Introduction

Astrophysical black holes (BHs) are classified according to their mass as stellar-mass BHs (10 2 M e ), intermediate-mass black holes (IMBHs; ∼10 2 -10 5 M e ), or supermassive black holes (SMBHs; 10 5 M e ). IMBHs could play a fundamental role in a wide variety of contexts: They could be the seeds that grow into supermassive BHs at the centers of galaxies, providing feedback on galaxy evolution, and a possible source of reionization (e.g., Madau & Rees 2001;Silk 2017;Tagawa et al. 2020;Natarajan 2021); they can produce observable tidal disruption events (e.g., Rosswog et al. 2009;Chen et al. 2011;MacLeod et al. 2016;Fragione & Leigh 2018a); they can be in binaries that become gravitational wave (GW) sources with unique properties (e.g., Miller 2002;Mandel et al. 2008;Gair et al. 2011;Fragione & Leigh 2018b;Arca Sedda et al. 2021); and they could be accreting from a stellar binary companion, producing ultraluminous X-ray sources (e.g., Kaaret et al. 2017). Therefore, characterizing the population of IMBHs represents a crucial step toward our understanding of the universe (see Greene et al. 2020, for a review)

There are four main observational methods to detect IMBHs. The first two consist of tracking stellar and gas dynamics and modeling gas accretion around a massive BH, respectively. Both methods have provided several massive ( ∼10 4 -10 5 M e ) nearby candidates (e.g., Baldassare et al. 2015;Chilingarian et al. 2018;Mićić et al. 2022;Pechetti et al. 2022), but none of them have confirmed the existence of an IMBH beyond reasonable doubt. The third method is based on looking for tidal disruption events consistent with an IMBH, similarly to what is typically done for supermassive BHs in galactic nuclei (e.g., Lin et al. 2018;Peng et al. 2019;Shen 2019). Detecting a signal in the outskirts of a galaxy and from the disruption of a white dwarf would be promising evidence for the existence of an IMBH (e.g., Rosswog et al. 2008Rosswog et al. , 2009;;MacLeod et al. 2016), and new facilities, such as JWST and LSST/VRO, may be able to detect tens of these events out to redshift z ∼ 1. The fourth method consists of looking for GW events, where one of the components in the merging binary is in the IMBH mass range. While LIGO/Virgo/KAGRA can detect the mergers of IMBHs with masses ∼100 M e out to z ∼ 1 (e.g., Abbott et al. 2019), as already done in the notable case of GW190521 (Abbott et al. 2020); the upcoming LISA, Einstein Telescope (ET), TianQin, and DECIGO have the potential to detect IMBHs throughout the universe (e.g., Luo et al. 2016;Amaro-Seoane et al. 2017;Jani et al. 2020). Current LIGO/Virgo/ KAGRA upper limits for mergers of IMBHs with masses up to about 500 M e are ∼0.1-10 Gpc -3 yr -1 (Abbott et al. 2017(Abbott et al. , 2019(Abbott et al. , 2022)).

Three main pathways to form IMBHs have been discussed in the literature. The first two channels involve the direct collapse of a gas cloud of pristine gas (Loeb & Rasio 1994;Bromm & Loeb 2003;Begelman et al. 2006) or of massive Population III stars (Fryer et al. 2001;Madau & Rees 2001;Bromm & Larson 2004) at high redshift, which might form an IMBH of ∼10 4 -10 5 M e and ∼100 M e , respectively. IMBHs with masses in between these two extremes can be produced at any redshift through repeated mergers either of massive mainsequence stars, later collapsing to form a BH (Portegies Zwart & McMillan 2002;Gürkan et al. 2004;Freitag et al. 2006;Pan et al. 2012;Giersz et al. 2015;Tagawa et al. 2020;Das et al. 2021;Di Carlo et al. 2021), or of stellar-mass BHs (Miller & Hamilton 2002;O'Leary et al. 2006;Antonini & Rasio 2016;Antonini et al. 2019;González et al. 2021;Mapelli et al. 2021;Weatherford et al. 2021;Fragione et al. 2022). Other possibilities include the fragmentation of active galactic nucleus (AGN) disks (McKernan et al. 2012(McKernan et al. , 2014)), super-Eddington accretion onto stellar BHs embedded in AGN disks (e.g., Kocsis et al. 2011), and repeated mergers of BHs with stars in galactic nuclei (e.g., Stone et al. 2017;Rose et al. 2022).

GWs provide the most secure way to weigh IMBHs since their mass is strictly encoded in the chirp mass and the characteristic strain of the signal. Only a few studies systematically explore the growth and formation of IMBHs in the context of repeated mergers (e.g., Antonini et al. 2019;Mapelli et al. 2021;Fragione et al. 2022), and even fewer predict their merger rates to compare to present and future constraints from GW detectors (e.g., Fragione et al. 2018aFragione et al. , 2018b;;Rasskazov et al. 2020). However, IMBH merger rates are typically predicted in the range ∼0.01-10 Gpc -3 yr -1 , which would imply several detectable mergers per year with current and upcoming detectors. In this paper, we use the results from a large set of semianalytic calculations, performed with the code developed in Fragione et al. (2022), to study the formation of IMBHs via repeated mergers, and we compute for the first time the cosmic rate of IMBH mergers in nuclear star clusters (NSCs). We present merger rates as a function of the initial seed mass and the characteristic BH spin at birth, and we make predictions about the multiband detectability of mergers of IMBH binaries.

This paper is organized as follows. In Section 2, we discuss our numerical framework to study the formation and mergers of IMBHs. In Section 3, we present and discuss our results for merger rates and multiband detectability. Finally, in Section 4, we discuss the implications of our results and draw our conclusions.

Method

We consider NSCs since they represent the ideal environment to form IMBHs via repeated mergers owing to their large escape speeds (e.g., Antonini et al. 2019;Mapelli et al. 2021;Fragione et al. 2022). In what follows, we summarize the semianalytical scheme developed in Fragione et al. (2022) that we adopt here.

We assume that the NSC density is described by a threeparameter potential-density pair (for details see Stone & Ostriker 2015). To generate a population of NSCs, we start from sampling galaxy masses (M *,gal ) from a Schechter function

where M c = 10 11.14 M e and α c = 1.43 (Furlong et al. 2015). To scale the galaxy mass to the NSC mass (M NSC ), we use scaling relations for late-type galaxies from Georgiev et al. (2016),

and, to sample their half-mass-radius (r h ), we use

where c 1 = 2.78 × 10foot_1 M e , c 2 = 3.94 × 10 9 M e , 1.001 0.067 andc 3 =3.31 pc, c 4 = 3.60 × 10 6 M e , 0 . Note that we consider the scatter in the fit parameters in sampling from Equations (2) to (3).foot_0 From the NSC mass and size, we compute the escape velocity from the center (e.g., Fragione & Silk 2020;Fragione et al. 2022)

Hereafter, we summarize our numerical procedure to follow the growth and mergers of IMBHs, starting from a BH seed of mass M seed , which undergoes mergers with other stellar-mass BHs (for details see Fragione et al. 2022).

1. Initial seed mass. We either fix the mass of the growing BH seed to some initial value (e.g., produced through repeated mergers of massive stars; González et al. 2021;Weatherford et al. 2021) or we derive it directly from stellar evolution. In the latter case, we sample the stellar mass of the seed progenitor from a Kroupa (2001) initial mass function 6 and evolve it using SSE (Hurley et al. 2000;Banerjee et al. 2020). If the seed is not ejected by natal kicks, v natal , as a result of asymmetric supernova explosion, that is v natal < v esc , we compute the timescale for the seed to sink to the cluster center via dynamical friction (Chandrasekhar 1943). 2. Formation of binaries and mergers. We estimate the timescale for the seed to find a BH companion by accounting for three-body binary formation and binarysingle encounters, whichever is the faster depending on the NSC mass and density (e.g. For details about the various parameters and important quantities adopted in our model, see Table 1 in Fragione et al. (2022).

We run three different models corresponding to different assumptions for the initial mass of the BH seed. In the first model, we randomly sample the initial seed mass from a realistic mass spectrum for stellar BH remnants. In this case, the average initial seed mass would be about 10 M e , with little dependence on the cluster metallicity. In the second and third models, we assume an initial seed mass of 50 and 100 M e , respectively. The latter two models could represent the case where a massive seed is produced in the beginning of the cluster lifetime as a result of the collapse of a very massive star, formed through repeated stellar mergers (e.g., Portegies Zwart & McMillan 2002;González et al. 2021). The birth spin of BHs is quite uncertain. To encompass all the possibilities, we adopt four different fixed values for the initial spin parameter of BHs, namely χ = 0, 0.2, 0.5, 0.8. For each model, we consider eight different metallicities (Z = 0.0001, 0.0002, 0.0005, 0.001, 0.002, 0.005, 0.01, and 0.02). For each combination, we run 50 k simulations, for a total of 4.8 × 10foot_3 .

Results

Merger Rates

We compute the differential merger rate of IMBH binaries as follows. We start with computing the merger rate as

where z 6 max = , t lb is the look-back time at redshift z 7 , Ψ NSC (z) is the cosmic NSC formation history, Φ(z, Z) is the merger efficiency at a given metallicity Z, and Π(z, Z) is the metallicity distribution at a given redshift, which we assume is described by a log-normal distribution with mean given by Madau & Fragos (2017) Z z log Z 0.153 0.074 , 6

1.34

and a standard deviation of 0.5 dex (Dvorkin et al. 2015). The formation history of NSCs is highly uncertain (e.g., Neumayer et al. 2020). Here, we model it as

In the previous equation, we assume et al. 2021), chosen so that the NSC density in the local universe is consistent with the observed one. We also fix z NSC = 3.2 and σ NSC = 1.5, respectively, following the cosmic formation history of globular clusters (e.g., Gratton et al. 1997Gratton et al. , 2003;;VandenBerg et al. 2013;El-Badry et al. 2019), under the assumption that NSCs form mostly from the inspiral of globular clusters into the galaxy through dynamical friction (e.g., Capuzzo-Dolcetta & Miocchi 2008;Antonini 2014). Finally, we model the merger efficiency at a given metallicity Z as (e.g., Mapelli et al. 2021)

where ò(z, Z) is the number of mergers of IMBH binaries per simulated system at a given redshift, N BH (Z) the total number of BHs at a given metallicity, assuming a Kroupa (2001) initial mass function, and M tot (Z) is the total simulated NSC mass at metallicity Z. We compute Equation (5) for different IMBH mass bins, that is R i (z).

Figure 1 shows the volumetric merger rate for IMBH binaries for different initial choices of the seed mass, assuming that the initial BH spins are zero. In the first model, we find that IMBH merger rates are peaked at z  1, and are about 1.5 Gpc -3 yr -1 , 1 Gpc -3 yr -1 , 2 Gpc -3 yr -1 , 0.4 Gpc -3 yr -1 , 0.5 Gpc -3 yr -1 for IMBH in the mass bins 100 M e -200 M e , 200 M e -500 M e , 500 M e -800 M e , 800 M e -1000 M e , 1000 M e -2000 M e , 2000 M e -5000 M e , respectively, in the local universe. The assumption on the initial seed mass affects the overall merger rate and the relative rates of different mass bins. For M seed = 50 and 100 M e , the overall merger rate increases by a factor of about 5 and 25, respectively, and the merger rate of IMBHs with masses 500 M e becomes comparable to the merger rate of IMBHs with masses 500 M e . This trend can be explained considering that the recoil kick imparted to the growing seed becomes smaller for larger initial seed masses, therefore rendering less likely the ejection of the growing IMBH. Our estimates are consistent (within current uncertainties) with LIGO/Virgo/KAGRA upper limits for mergers of GW190521-like systems (The LIGO Scientific Collaboration & the Virgo Collaboration 2020) and IMBH binaries with masses up to about 500 M e (Abbott et al. 2017(Abbott et al. , 2019(Abbott et al. , 2022)).

We show the effect of the initial BH spin on the expected merger rate in Figure 2. Since recoil kicks critically depend on BH spins (as long as the mass ratio of a merger is 0.1), the rates are significantly affected by spins. If BHs are born with nonzero spins, the growing seed can be more efficiently ejected from the host star cluster, halting further growth. Since recoil kicks also depend on the mass ratio of the merging binary, merger rates are less affected by a nonzero initial BH spin for larger initial seed masses. For example, for masses in the range 100-200 M e , the merger rates become almost two orders of magnitude smaller when the initial BH spin parameter goes from 0.2 to 0.8 for an initial random seed, while they go down by a much smaller factor of about four when the initial seed mass is 100 M e . Finally, the effect of the initial BH spin is more important on the merger rates of larger IMBHs. Indeed, larger initial spins imply larger recoil kicks, which can eject the growing seed from the parent cluster, reducing the merger rates of large IMBHs.

Dependence on NSC Formation History

The formation history of NSCs is highly uncertain. As shown in Equation (7), we model it as a Gaussian distribution with mean z NSC = 3.2 and variance σ NSC = 1.5. To analyze how our results depend on these assumptions, we run two additional models where we set z NSC = 1.6 ( M 4.9 10 Mpc yr

=

´---



), respectively. These two additional models represent a later and earlier NSC formation peak, respectively, with respect to our main model.

We show the results of our additional runs in Figure 3, in the case of a random initial seed mass in the mass spectrum of stellar-mass BHs. We find that our various formation histories do not have a significant impact on the magnitude of the IMBH merger rate, while affecting its shape for masses 500 M e .

Number of Mergers

The number of mergers we observed per year is given by the cumulative merger rate

where (dV c (z)/dz) is the amount of comoving volume in a slice of the universe at redshift z and (1 + z) -1 is the difference in comoving time between the merger redshift and the observer at z = 0.

In Figure 4, we show the cumulative merger rate as a function of the IMBH mass for different choices of the initial seed mass and initial BH spins, at redshift z = 1 (squares) and z = 3 (stars). As expected from the analysis of Figures 12, most of the binary mergers take place at low redshift, z  1, with the exception of the smallest IMBHs with masses 300 M e . In this case, the number of events is higher by a factor of a few at z = 3 with respect to z = 1, with the largest differences attained for larger initial seed masses and smaller initial BH spins. Mergers of IMBHs with masses 300 M e occur at z < 1, as a result of the progressive building up of massive IMBHs.

Multiband Detections

We compute the signal-to-noise ratio (S/N) of our merging systems in a detector frequency band, to understand the prospects for multiband detections by ground-based and spacebased observatories. Given the masses of the merging BHs, m 1 and m 2 , we compute the average S/N as

where f min and f max are the minimum and maximum frequency of the binary in the detector band, respectively, S n ( f ) is the effective noise power spectral density, and h f | ˜( )| is the frequency-domain waveform amplitude for a face-on binary, approximated with a PhenomA waveform (e.g., Equation (20) in Robson et al. 2019)

where the values of {f k , a k , b k , c k } are taken from Table 2 in Robson et al. (2019), f is the detector-frame frequency, related to the binary orbital frequency by f = f orb (1 + z) -1 , M c,z is the redshifted chirp mass, related to the rest-frame chirp mass by M c = M c,z /(1 + z), and

is the luminosity distance, with c and H 0 being the velocity of light and Hubble constant, respectively, and Ω M = 0.286 and Ω Λ = 0.714 (Planck Collaboration et al. 2016), respectively. We compute the power spectral density of LISA as in Robson et al. (2019), of ET as in Hild et al. (2011), of DECIGO as in Yagi & Seto (2011), and of LIGO at design sensitivity as in Ajith (2011). In order to assess the possibility of multiband detections, we choose f f f max , min mb min,j = ( ) , where f min,j is the minimum frequency in a detector band (10 -5 Hz,

is the minimum initial frequency for a binary to merge within the LISA mission lifetime, which we set to T LISA = 4 yr. We set as a minimum detection threshold 〈S/N〉 thr = 8. We show in Figure 5 the expected average S/N of IMBH binaries as a function of redshift for different detectors and various masses of the IMBH and its companion. The average S/N in both LIGO and ET is generally higher for larger companion masses and smaller IMBH masses, as expected since binaries with more massive IMBHs will spend fewer cycles in band. With its smaller characteristic noise, ET will be able to detect these binaries further than LIGO. For instance, a merging IMBH binary is expected to be detected with an S/N  10 by ET even at z = 2, while the corresponding S/N in LIGO is 10. On the contrary, DECIGO and LISA will be able to detect merging binaries to larger distances for a larger IMBH and secondary masses. While DECIGO is expected to detect such systems with S/N  10 3 even at redshift z = 2, LISA will essentially be able to detect merging binaries with IMBH mass 800 M e and companion masses 30 M e at z  0.5. 8We show in Figure 6 the maximum number of detected merger events per year as a function of the IMBH mass for different GW observatories, for various choices of the initial seed mass and BH spin. As discussed, most of the binary mergers take place at low redshift, z  1, with the exception of the smallest IMBHs with masses 300 M e , for which a nice fraction of events comes from larger redshifts. Our model predicts that a few merger events per year should be detectable with LISA, DECIGO, ET, and LIGO for IMBHs with masses 1000 M e , and a few tens of merger events per year with DECIGO, ET, and LIGO only.

Discussion and Conclusions

GW detectors are expected to revolutionize our understanding of the elusive IMBHs, after establishing firmly their existence.

We have analyzed the statistical results of a wide range of simulations, which encompass different metallicities, initial seed masses, initial BH spins, run with the semianalytical code developed in Fragione et al. (2022). We have computed the merger rate of IMBH binaries and have found that rates are in the range 0.01 Gpc -3 yr -1 -10 Gpc -3 yr -1 for different IMBH masses. We have also computed the number of multiband detections in ground-based and space-based observatories. Our model predicts that a few merger events per year should be detectable with LISA, DECIGO, ET, and LIGO for IMBHs with masses  1000 M e , and a few tens of merger events per year with DECIGO, ET, and LIGO only. Our estimates are consistent (within current uncertainties) with LIGO/Virgo/ KAGRA upper limits for mergers of GW190521-like systems (The LIGO Scientific Collaboration & the Virgo Collaboration 2020) and IMBH binaries with masses up to about 500 M e (Abbott et al. 2017(Abbott et al. , 2019(Abbott et al. , 2022)). As the sensitivity of current detectors is improved and new observatories start operating, our models predict tens or hundreds of merging binary signals in the next few years. In the case of fewer detections, the inferred observational merger rates can be used to constrain the model parameters of NSCs, their formation histories, and IMBH spins.

As already widely discussed in Fragione et al. (2022), there are some limitations to our approach. For example, we have assumed that the masses and sizes of NSCs do not evolve significantly during their lifetime and that NSCs are isolated with respect to the rest of their host galaxy history. In reality, NSCs are not isolated from their environments and evolve on a cosmological timescale, as there can be a continuous supply of stars and gas from the rest of the galaxy, with ongoing star formation and accretion of smaller star clusters. The NSC