Supermassive Black Hole Growth in Hierarchically Merging Nuclear Star Clusters Konstantinos Kritos1aa, Ricarda S. Beckmann2,3aa, Joseph Silk1,4,5aa, Emanuele Berti1aa, Sophia Yi1aa, Marta Volonteri6aa, Yohan Dubois6aa, and Julien Devriendt7aa 1William H. Miller III Department of Physics and Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA; kkritos1@jhu.edu 2 Institute of Astronomy and Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK 3 Institute for Astronomy, University of Edinburgh, Royal Observatory, Edinburgh EH9 3HJ, UK 4 Institut d’Astrophysique de Paris, UMR 7095 CNRS and UPMC, Sorbonne Université, F-75014 Paris, France 5 Department of Physics, Beecroft Institute for Particle Astrophysics and Cosmology, University of Oxford, Oxford OX1 3RH, UK 6 Institut d’Astrophysique de Paris, CNRS, Sorbonne Université, UMR7095, 98bis bd Arago, 75014 Paris, France 7 Department of Physics, University of Oxford, Keble Road, Oxford, OX1 3RH, UK Received 2025 January 7; revised 2025 May 21; accepted 2025 June 6; published 2025 September 16 Abstract Supermassive black holes are prevalent at the centers of massive galaxies, and their masses scale with galaxy properties, increasing evidence suggesting that these trends continue to low stellar masses. Seeds are needed for supermassive black holes, especially at the highest redshifts explored by the James Webb Space Telescope. We study the hierarchical merging of galaxies via cosmological merger trees and argue that the seeds of supermassive black holes formed in nuclear star clusters via stellar black hole mergers at early epochs. Observable tracers include intermediate-mass black holes, nuclear star clusters, and early gas accretion in host dwarf galaxies, along with a potentially detectable stochastic gravitational-wave background, ejection of intermediate and supermassive black holes, and consequences of a significant population of early tidal disruption events and extreme mass ratio inspirals. Unified Astronomy Thesaurus concepts: Supermassive black holes (1663) 1. Introduction Supermassive black hole (SMBH) formation is an out- standing problem in contemporary cosmology and astrophy- sics. SMBHs often coexist with nuclear star clusters (NSCs), the densest known stellar systems in the Universe, found in galactic nuclei (N. Neumayer et al. 2020). Moreover, scaling relations between the SMBH’s and NSC’s mass have been suggested (I. Y. Georgiev et al. 2016), which may imply a coevolution history, and the NSC might have contributed to the formation of black hole (BH) seeds that further grow into SMBHs. Thanks to its infrared capabilities, the James Webb Space Telescope (JWST) has revolutionized our understanding of the high-redshift Universe, and a large population of distant SMBHs has been revealed (J. E. Greene et al. 2024; I. Juodžbalis et al. 2024; O. E. Kovács et al. 2024; R. Maiolino et al. 2024a, 2024b). It is generally believed that SMBHs have grown out of smaller BH seeds that may be light or massive. The former designates seeds of stellar origin with masses up to 100M⊙, and the latter refers to seeds of intermediate BH mass (IMBH) in the ∼103–105M⊙ range regime (J. E. Greene et al. 2020). The existence of stellar-mass BHs is observationally well established, and that of massive seeds has theoretical backing from Population III studies motivated by tentative detections of spectroscopic signatures from a metal-free stellar population at high redshift, including strong HHEII emission, extremely blue UV continua, and enhanced [N/O] (B. Liu et al. 2024). The observation of SMBHs at high redshift (z > 6) implies that seed formation and growth must have occurred rapidly within a few hundred Myr (at least for a subset of the observed massive BHs). Moreover, a number of strongly lensed magnified massive (masses ∼106M⊙) compact (radii <1 pc) young star clusters have been detected by JWST at z ∼ 10 (A. Adamo et al. 2024), and these could be the sites of massive BH seed formation. According to simulations, the assembly of such massive star clusters is much more efficient at z > 7 (L. Mayer et al. 2025). Theoretical works that study the growth of SMBHs in the cosmological context rely on either expensive cosmological simulations or semianalytic methods, which are faster to implement. A number of previous works have studied the formation and evolution of massive BHs in the cosmological context (T. Di Matteo et al. 2008; S. Bonoli et al. 2009; P. Natarajan 2011; Y. Dubois et al. 2014a, 2014b; D. Sijacki et al. 2015; M. Tremmel et al. 2015; M. Volonteri et al. 2020; D. Izquierdo-Villalba et al. 2021; C. A. Dong-Páez et al. 2023a, 2023b, 2025; A. K. Bhowmick et al. 2024a, 2025; J. Ellis et al. 2024; K. Li et al. 2024). Cosmological simulations follow the three-dimensional coupled evolution of dark matter, gas, stars, and massive BHs. Due to the many orders of magnitude of length scale required to be simulated, simulations are typically too expensive to capture scales below 100 pc over long periods of time. Therefore, simulating the formation process of massive BH seeds is nontrivial because they occur on spatial and timescales much smaller than can be resolved. Consequently, these cosmological simulations must rely on extrapolated models to place seed BHs in galaxies. At the same time, gravitational runaway processes occur in the cores of dense star clusters at the subparsec scale and operate on timescales of a few Myr. Therefore, in this work, we simulate the growth of these light BH seeds assuming that they are initially formed by the usual astrophysical process of massive star evolution. In previous work, K. Kritos et al. (2024a) developed NUCE, a code that self-consistently evolves NSCs, and studied the The Astrophysical Journal, 991:58 (18pp), 2025 September 20 https://doi.org/10.3847/1538-4357/adeb44 © 2025. The Author(s). Published by the American Astronomical Society. aaaaaaa Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. 1 https://orcid.org/0000-0002-0212-3472 https://orcid.org/0000-0002-2850-0192 https://orcid.org/0000-0002-1566-8148 https://orcid.org/0000-0003-0751-5130 https://orcid.org/0000-0002-9104-1734 https://orcid.org/0000-0002-3216-1322 https://orcid.org/0000-0003-0225-6387 https://orcid.org/0000-0002-8140-0422 mailto:kkritos1@jhu.edu http://astrothesaurus.org/uat/1663 https://doi.org/10.3847/1538-4357/adeb44 https://crossmark.crossref.org/dialog/?doi=10.3847/1538-4357/adeb44&domain=pdf&date_stamp=2025-09-16 https://creativecommons.org/licenses/by/4.0/ formation of SMBHs in NSCs, simulating the seeding process. The only true seeds that must be put by hand in NSCs are stellar-mass BHs with mass ∼10M⊙ that form after the evolution of massive stars. A number of works have emphasized the importance of a dense stellar cluster environment in rapidly forming massive BH seeds through mergers and accretion (G. D. Quinlan & S. L. Shapiro 1987; S. F. Portegies Zwart & S. L. W. McMillan 2002; K. Omukai et al. 2008; B. Devecchi et al. 2010, 2012; T. Alexander & P. Natarajan 2014; M. Mapelli 2016; B. Reinoso et al. 2018; F. Antonini et al. 2019; P. Kroupa et al. 2020; P. Natarajan 2021; K. Kritos et al. 2023; M. C. Vergara et al. 2023; D. Mehta et al. 2024; L. R. Prole et al. 2024; A. Rantala et al. 2024; A. Dekel et al. 2025; C. Partmann et al. 2025). Alternative scenarios also exist, such as primordial BHs (D. Hooper et al. 2024; F. Ziparo et al. 2025) and the growth of Population III stars (K. Freese et al. 2010; V. Cammelli et al. 2025). In this paper, we present a semianalytic model for the evolution of NSCs and their central massive BHs based on galaxy merger trees extracted from the NEWHORIZON simula- tion (Y. Dubois et al. 2021); see, e.g., M. Polkas et al. (2024) for another semianalytic model of NSCs. In Section 2, we present our methodology for coupling the NUCE code with those merger trees; in Section 3, we show our results on the generation of massive BH and coexisting NSC populations; in Section 4, we suggest that associated stellar tidal disruption events (TDEs) in NSCs will be an inevitable corollary of our model and discuss more generally the caveats of our study; and in Section 5, we present our conclusions.8 2. Methods In this section, we present our assumptions and methodol- ogy. In Section 2.1, we introduce the NEWHORIZON simulation from which we extract merger trees. In Section 2.2, we discuss NSC initial conditions. In Section 2.3, we describe the formation of massive BH seeds in NSCs. Finally, in Section 2.4, we present our algorithm for treating NSC–NSC coalescence and massive BH binary formation and merger during major galaxy–galaxy mergers. 2.1. Merger Trees In this paper, we build our model of NSCs and IMBHs on galaxy evolution histories extracted from the NEWHORIZON simulation. NEWHORIZON has been presented in detail in Y. Dubois et al. (2021), so we only briefly recap its relevant features here. NEWHORIZON is a high-resolution resimulation of a sphere of radius ∼16Mpc of the HORIZON-AGN simulation. NEW- HORIZON was performed with RAMSES (R. Teyssier 2002) and uses a ΛCDM cosmology consistent with WMAP-7 data (E. Komatsu et al. 2011). The simulation reached a maximum resolution of 34 comoving pc using a quasi-Lagrangian and super-Lagrangian refinement scheme. Gas follows an ideal monoatomic equation of state with an adiabatic index of γad = 5/3. Cooling is modeled down to 104 K using cooling curves from R. S. Sutherland & M. A. Dopita (1993) assuming equilibrium chemistry. After redshift zreion = 10, heating from a uniform UV background is included following F. Haardt & P. Madau (1996). Star formation takes place following a Schmidt relation (T. Kimm et al. 2017; M. Trebitsch et al. 2017, 2020) in cells with a gas number density above n0 = 10 H cm−3 using a star formation efficiency that depends on the properties of the interstellar medium such as the Mach number and virial parameters. Stars have a stellar mass resolution of 1.3 × 104M⊙ and are assumed to have a Chabrier (G. Chabrier 2005) initial mass function with cutoffs at 0.1M⊙ and 150M⊙. Stellar feedback separately tracks the momentum and energy-conserving phase of the explosion (T. Kimm et al. 2015). NEWHORIZON contained on-the-fly evolution for BH formation and evolution, but as its properties are not included in this analysis, we omit discussing this further. Interested readers are referred to R. S. Beckmann et al. (2023). Halos and galaxies were identified in the simulation using the HOP structure finding algorithm (D. Aubert et al. 2004), using an overdensity of δhalos = 80 (for halos) and δgalaxies = 160 for galaxies. Uncontaminated halos only contain high-resolution dark matter particles within their virial radius. Our sample of target galaxies consists of all central galaxies with a stellar mass above 6.5 × 105M⊙ (equivalent to at least 50 star particles) that are located within the central 0.1 virial radii of a dark matter halo at z = 0.25. In Figure 1, we show the distribution and range of stellar galaxy masses of all 451 galaxies within the simulation box extracted from the NEWHORIZON simulation at z = 0.25. The sample contains a large population of dwarf galaxies with M� < 1010M⊙ and 17 massive galaxies. From this sample of galaxies, we construct merger trees across cosmic time (see Figure 2 for an example) that show the stellar mass evolution history of the galaxy’s progenitors across cosmic time. Merger trees are constructed by analyzing a set of 997 simulation snapshots in the redshift range z = 10–0.25. We only retain galaxies that attain a mass of at least 107M⊙ in stellar mass before the merger and discard all others. Throughout the 105 106 107 108 109 1010 1011 1012 M?/M� 0 5 10 15 20 25 30 C ou n t z = 0.25 Figure 1. Distribution of all galaxy stellar masses extracted from the NEWHORIZON data set at redshift z = 0.25. The sample contains 451 galaxies in a box of comoving volume (16 Mpc)3. 8 In this manuscript, we use massive BH to refer to any BH with a mass >100 M⊙. However, since we focus on a mass range with more than 10 orders of magnitude, it is common practice to further categorize massive BHs into IMBHs in the range 100–106 M⊙ and SMBHs with masses >106 M⊙. 2 The Astrophysical Journal, 991:58 (18pp), 2025 September 20 Kritos et al. simulations carried out in this study, the evolution of galaxy stellar mass is tracked from the snapshots of the NEWHORIZON simulation. The NEWHORIZON simulation originally included massive BH-galaxy evolution; however, in this study, we evolve BH growth within NSCs along the static branches and nodes of the extracted merger trees. 2.2. NSC Initial Conditions In our NSC paradigm, SMBH seeds form from repeated merger processes and accretion in the dense centers of protogalaxies. We assume that each protogalaxy forms an NSC in a star formation burst when its stellar mass reaches a critical value,Mcr. The initial stellar massM�,0 contained in the NSC is then taken to be a random fraction of Mcr, denoted by [ ]f f0, ,max , and we write M�,0 = f�Mcr. Assuming a Kroupa initial mass function in the range 0.08–150M⊙ with the mean parameters from P. Kroupa (2002; which we assume to be universal), the initial average stellar mass is =m M0.6,0 . The initial number of stars is then computed as /=N M m,0 ,0 ,0, which is equivalent to ( )/=N f M m . 1,0 cr ,0 An amount of gas is also initially added into the NSC such that a random fraction [ ]f f0,g g,max of the cluster’s mass is in the form of gas that did not fragment into stars. Thus, if Mcl,0 = M�,0 + Mg,0 is the total initial mass of the NSC, then Mg,0 = fgMcl,0 is the initial amount of gas in the system not turned into stars assumed to be in the form of ionized hydrogen with a sound speed of cs = 10 km s−1. We compute the initial gas mass as ( ) ( )/=M M f f1 , 2g,0 ,0 g g and the star formation efficiency is given by η� = 1 − fg. The gas is completely removed from the NSC with an exponential time dependence on a timescale τge, which is a free parameter of our model (H. Baumgardt & P. Kroupa 2007). The prior on the initial half-mass radius rh,0 is drawn from a log-uniform distribution from 0.1 to 10 pc. This choice is supported by young massive star cluster properties in the local Universe (S. F. Portegies Zwart et al. 2010) and by recent observations of high-redshift young star clusters with effective radii on parsec and subparsec scales (E. Vanzella et al. 2023; A. Adamo et al. 2024). Moreover, according to B. Devecchi et al. (2010), NSCs may form in protogalaxies with rh,0 < 0.5 pc. Since most NSCs observed in the local Universe have effective radii of a few parsecs (N. Neumayer et al. 2020), and, by Hénon’s principle, energy generation in their core causes them to expand during balanced evolution, these systems must have been more compact in the past when they formed. The diversity of radii in present-day NSCs indicates a variety of initial radii and possibly different evolutionary conditions. We evolve NSCs and simulate the formation of SMBH seeds using the numerical code NUCE (for “nuclear cluster evolution”), a semianalytic model that relies on Hénon’s principle, described in K. Kritos et al. (2024a). 2.3. Massive BH Seeds In this subsection, we describe the BH growth channels considered in our model (Section 2.3.1). Then we show the evolution of the mass and spin of the seeds (Section 2.3.2), and finally, we present the dependence of seed formation on the NSC initial conditions (Section 2.3.3). 2.3.1. BH Growth Channels Our NSC model NUCE implements the following seeding channels. 1. Stellar collisions. Successive physical collisions among stars in the cluster lead to the formation of a very massive star (>150 M⊙), which then directly collapses into a massive BH seed. This runaway process initiates provided that 0.2τrh,0 < τse = 3Myr, where τrh,0 is the initial half-mass relaxation timescale and τse is the massive star evolution timescale. The mass of the massive object is a fraction ∼10−3 of the initial stellar mass (S. F. Portegies Zwart & S. L. W. McMillan 2002; F. A. Rasio et al. 2004). Since runaway stellar collisions rapidly occur in the first few Myr if the previous condition is met (M. S. Fujii et al. 2024), we do not simulate them; rather, we assume that a massive BH seed has formed in the cluster with mass =M N m10BH 3 ,0 ,0. This BH may then grow further through the following channels. If the runaway condition is not met, then the seed BH mass is set to be 10M⊙. 2. Gas accretion. We assume that gas is accreted onto BHs at the Bondi rate (H. Bondi 1952) capped to the Eddington limit and that accretion operates only for as long as there is remaining gas in the cluster. For the gas density, we take ( )/= M r2g g h 3 , and we assume a sound speed of cs = 10 km s−1 (K. Kritos et al. 2024a). We recall that the gas is removed exponentially with decay timescale τge (H. Baumgardt & P. Kroupa 2007). 3. BH mergers. Binary BHs form after the core collapse of the BH subsystem, and if a binary does not receive an interaction recoil that exceeds the escape velocity vesc, it then merges within the host cluster as a consequence of hardening. The binary is ejected if the gravitational-wave 0 2 4 6 8 10 12 t/Gyr 106 107 108 109 1010 M ? /M � 8 4 2 1 0.5 0.25 z Figure 2. Example galaxy merger tree showing the galaxy stellar mass evolution of the main progenitor (black), as well as all other progenitors (colored lines) that reach a mass of M� > 107 M⊙ before the galaxy merger. Colored markers indicate galaxy mergers. 3 The Astrophysical Journal, 991:58 (18pp), 2025 September 20 Kritos et al. (GW) kick vGW imparted to the postmerger remnant is such that vGW > vesc; otherwise, it is retained in the cluster, and it may merge with another BH. These repeated BH mergers lead to BH growth. 4. Star consumption. The contribution of repeated stellar consumption events becomes nonvanishing in our model only after the BH subsystem has evaporated, within approximately 10 initial half-mass relaxation times (P. G. Breen & D. C. Heggie 2013). We compute the full-loss-cone stellar consumption rate implementing Equation (6.15) from D. Merritt (2013) assuming a Bahcall–Wolf cusp. The loss-cone radius is set to double the tidal radius to account for tidal captures (F. P. Rizzuto et al. 2023). We assume that a mass of m0.05 ,0 is accreted. 2.3.2. Time Evolution of MBH and χBH In this subsection, we test the massive BH growth model in NSCs without using the merger trees described in Section 2.1. In Figure 3, we show the time evolution of the mass (MBH) and spin (χBH) of massive BH seeds when we consider different initial conditions for (N�,0, rh,0, η�), i.e., the initial number of stars, half- mass radius, and star formation efficiency, respectively. We compute the time evolution of the BH spin from the relation ( )= d dt d dM dM dt , 3BH BH BH BH where dχBH/dMBH is computed by following the classic work of J. M. Bardeen (1970); see Equation (1) in K. Kritos et al. (2024b). The factor dMBH/dt contains contributions from gas accretion and consumption of stars. We assume that accretion (of gas or stars) is coherent along the plane orthogonal to the direction of the BH spin, and we randomize the direction of the accretion (prograde or retrograde) with a probability of 50%. Below, we describe in detail the evolution of MBH and χBH in the specific cases shown in Figure 3. The case (106, 0.2 pc, 1.0). The cluster does not satisfy the conditions for runaway stellar collisions or repeated BH mergers to assemble an IMBH. Second-generation BHs (≃20M⊙) form and merge with other first-generation BHs (10M⊙). Still, the third-generation remnants of these mergers are ejected promptly due to the cluster’s relatively small escape speed, which is initially ≃185 km s−1. In particular, since second-generation merger products of comparable-mass BHs have a spin of ∼0.7 (V. Baibhav et al. 2021), and first- generation BHs are assumed to be nonspinning, the probability of retaining the third-generation product is about 4%. This repeated formation of second-generation BHs and subsequent ejections results in the oscillatory behavior seen in the mass and spin evolution plots. Once 90% of the BHs are ejected at ≃200Myr, repeated consumption of stars initiates and slowly leads to the formation of a ≈200M⊙ massive BH (within 1 Gyr). At the same time, the spin asymptotes to zero. The case (106, 0.1 pc, 1.0). The cluster has the required initial conditions to undergo runaway stellar collisions, and a ≈600M⊙ massive BH forms by t = 3Myr. We assume that the very massive stellar progenitor directly collapses into a BH without mass loss, and we (arbitrarily) set the spin of the direct-collapse BH to zero. By t = 1 Gyr, the IMBH further grows to ≈2000M⊙ by merging with other BHs and consuming stars, while its spin hovers around χBH ∼ 0.2 within the simulated time window. The case (107, 0.2 pc, 1.0). The cluster is massive and compact enough for the initial escape velocity to exceed 300 km s−1, which is the threshold for hierarchical mergers to happen. A 103M⊙ massive BH is formed within a few tens of Myr and grows to ≈2 × 104M⊙ by t = 1 Gyr. The spin of this BH increases to ∼0.7 after the first merger and drops below χBH ∼ 0.1 by the end of the simulation as a result of the isotropic spin directions in merger events. The case (106, 0.5 pc, 0.5). As in the first case above, the probability of forming a massive BH from repeated BH mergers is small due to the insufficient escape velocity of the system. The system contains a lot of gas, but the gas accretion rate is too low to grow the BH, because the mass-doubling timescale (based on the Bondi rate) is smaller than the gas expulsion timescale. Only after the initiation of stellar consumption at ≈530Myr (following the evaporation of the BH subsystem) does the gas accretion lead to BH growth and spins it up to the maximal value. At ≈700Myr, the remnant gas has been removed completely, and further growth proceeds through stellar consumption. During this final phase, the BH spin drops, asymptoting to zero on the mass-doubling time as a consequence of repeated consumption of stars. The cases (107, 0.2 pc, 0.5) and (107, 0.2 pc, 0.1). We consider the same structural initial conditions as in the case (107, 0.2 pc, 1.0) above, but we consider different values for the star formation efficiency. In both cases, the seed BH grows by mergers and gas accretion that act simultaneously but on different timescales. Depending on the amount of gas present in the system, the final BH reaches ≃7 × 105M⊙ (η� = 0.5) or ≃4 × 106M⊙ (η� = 0.1) within 300Myr. The BH starts to slowly grow by gas accretion in the first ≃50Myr and spins up to the maximal value, but the accretion rate is enhanced once repeated BH mergers lead to BH masses of ∼103–104M⊙. During the hierarchical BH merger phase, the spin drops to 101 102 103 104 105 106 107 M B H /M � (106, 0.2 pc, 1.0) (106, 0.1 pc, 1.0) (107, 0.2 pc, 1.0) 100 101 102 103 t/Myr 0.00 0.25 0.50 0.75 1.00 χ B H (106, 0.5 pc, 0.5) (107, 0.2 pc, 0.5) (107, 0.2 pc, 0.1) Figure 3. The time evolution of BH mass (MBH; upper panel) and dimensionless spin (χBH; lower panel). Different colors correspond to various initial conditions for (N�,0, rh,0, η�), as shown in the legends. 4 The Astrophysical Journal, 991:58 (18pp), 2025 September 20 Kritos et al. low values of ∼0.1 (see the dips in the evolution of the spin in the orange and green lines in the lower panel of Figure 3). Once gas accretion takes over the evolution of the BH, it steadily and rapidly grows into the SMBH regime; the spin magnitude grows back to the maximal value within a mass- doubling time and remains at the maximal value, since we assume coherent gas accretion. We remind the reader that we have chosen to isotropically accrete stars and stellar-mass BHs onto the growing BH, while gas is accreted coherently. 2.3.3. Mapping BH Masses in the NSC Parameter Space In Figure 4, we show that the mass of the heavier BH formed within 1 Gyr for 104 randomly chosen initial condi- tions in the ranges 104� N�,0� 108 and 0.1� rh,0/pc� 10. The boundaries of massive BH formation are easier to understand in the dry (η� = 100%) case without gas accretion. Clusters with initial conditions below the lime dashed line of Figure 4 satisfy the condition for runaway stellar collisions and form a massive BH within 3Myr. Clusters with initial conditions below the orange dashed line eject their stellar- mass BHs and undergo runaway stellar consumption. The normalization constant of the orange dashed boundary depends on the time since the beginning of the simulation, and in all panels of Figure 4, we show the BH mass at t = 1 Gyr. Since the BH subsystem evaporates within around 10 initial half- mass relaxation times (P. G. Breen & D. C. Heggie 2013), the condition 10τrh,0 < t is required for runaway stellar consumption to have an effect. Otherwise, the BH subsystem Figure 4. Masses of the BH seeds MBH formed through the combined processes of successive mergers and gas accretion within 1 Gyr in NSCs as a function of the initial number of stars N�,0 and initial half-mass radius rh,0. Each point corresponds to an NSC simulation with NUCE. A few pairs of the star formation efficiency (η�) and gas expulsion timescale (τge) have been assumed and shown in the title of each panel. Note that in the lower left panel, no value of τge is given, because there is no gas to be ejected. The colored stars in the panels represent the set of initial conditions for the simulations of Figure 3. 5 The Astrophysical Journal, 991:58 (18pp), 2025 September 20 Kritos et al. dominates the core, and no massive BH forms unless repeated BH mergers occur. Thus, the parameter space of massive BHs increases with time. If the cluster’s initial conditions lie below the cyan dashed line—corresponding to an initial escape velocity of 300 km s−1—then repeated BH mergers can form a massive BH. This result agrees with the findings of F. Antonini et al. (2019) for the threshold for massive BH formation via repeated mergers. Now we discuss the rest of the panels of Figure 4. The presence of gas in the system boosts the growth of the massive BH, and the available parameter space expands as long as the gas remains in the system long enough, as we argue below. In the region below the yellow lines in the upper panels of Figure 4, BH interactions become efficient, and the binary formation rate is much smaller than 100Myr. The growing BH merges with other stellar-mass BHs in the core and eventually gets ejected due to GW recoil, so it does not have time to grow. Only after most BHs have been ejected is the successive consumption of stars turned on, which contributes to the formation of massive BHs. In the high-rh,0 region, the gas density is not high enough; hence, accretion at the Bondi rate becomes insufficient to grow the BH in the available 100Myr before it is removed, because we assume the gas density to be ( )/M r2g g h 3 (K. Kritos et al. 2024a). In the case τge = 10Myr (lower right panel of Figure 4), the gas is removed faster, and the cluster expands by a larger factor than if it was ejected at a slower rate (K. Kritos et al. 2024a). Consequently, when the BH subsystem forms, the escape velocity has dropped by a factor of ∼10, and the probability of ejection due to a GW kick is much higher (although the interaction rate also drops). Forming a massive BH through BH mergers requires a much higher escape velocity in this scenario, and the initial escape velocity of the cluster has to be larger (≳400 km s−1) because of the rapid expansion, which lowers vesc. Thus, most of the parameter space in this case does not lead to the production of massive BHs. 2.4. Major Merger Episodes We keep track of the merger times of galaxies in the merger tree. We classify the merger as major or minor based on a decision boundary depending on the mass ratio Q of their stellar masses (defined here as the secondary over the primary galaxy stellar mass). If the mass ratio is larger than some threshold value Qth, we assume that the merger is major, and that if both galaxies contain an NSC, their corresponding NSCs will coalesce. If both NSCs contain massive BHs, these massive BHs will pair up and eventually merge. On the other hand, if the asymmetry in the galaxy masses is smaller, we assume that the lighter galaxy dissolves in the primary halo, and the NSC becomes an ultracompact dwarf galaxy with a central massive BH. Following the merger of two galaxies, in the case of a major merger, the two NSCs are merged into a single NSC, assuming conservations of mass and energy. Hence, if M1, r1 and M2, r2 are the total masses and half-mass radii of the two progenitor NSCs, then the mass M and half-mass radius r of the merged NSC will be ( )= +M M M , 4a1 2 ( )= + r M , 4b M r M r 2 1 2 1 2 2 2 respectively, as in A. Lupi et al. (2014). Moreover, we add an amount of gas given by ( )/f M f1g g , where fg is randomly sampled from [ ]f0, g,max , and the choice of fg,max will be discussed later. In this model, we ignore the effect of any time delays between galaxy mergers, NSC mergers, and massive BH binary mergers. Following the coalescence of two massive BHs, we compute the remnant mass, remnant spin, and GW kick vGW imparted to the merger remnant using the fitting formulae of Section V of D. Gerosa & M. Kesden (2016), which are based on interpolation between numerical relativity fits and the test particle limit. The spin orientations of the progenitor BHs that coalesce following a major merger are randomized. The assumption of spin orientation is important because in certain configurations, when the relative angle of both spins with the orbital angular momentum is ≈45° and their projections on the orbital plane are antialigned, the merger remnant may receive a superkick of thousands of km s−1, which may well exceed the escape velocity of most galaxies (M. Campanelli et al. 2007; J. A. Gonzalez et al. 2007). The postmerger remnant BH is retained (ejected) if the magnitude of the GW recoil velocity, vGW, is smaller (larger) than the escape velocity of the merged NSC, vesc. We do not distinguish between massive BHs that receive a superkick (in thousands of km s−1) and those that are marginally ejected with vGW barely larger than vesc: if vGW > vesc, we assume that the ejected massive BH never returns to the center but rather wanders freely. To summarize, massive BHs in our simulations are characterized as being of the following kinds. 1. Nuclear. The massive BH lies in the galaxy’s center and is surrounded by an NSC. 2. Off-nuclear/wandering. These BHs, in turn, can be either of the following. (a) Satellite. The massive BH lies in the center of an ultracompact dwarf hosted in the halo of the primary galaxy following a minor merger episode. (b) Ejected. The massive BH is ejected from the nuclear region following a sufficiently large GW kick. 3. Results In this section, we present the model’s predictions for the properties of massive BHs and NSCs in the local Universe (Sections 3.1 and 3.2, respectively). In Section 3.3, we discuss transient events, such as massive and stellar-mass BH binary mergers and TDEs. We use the merger trees discussed in Section 2.1. Each merger tree terminates at z = 0.25, and our results are shown at that epoch. In our simulations, we choose Mcr = 107M⊙. This is motivated partly by local observations of NSC masses (see Figure 6 in N. Neumayer et al. 2020) and partly by the constraint that the initial stellar NSC mass must be a fraction of the total galaxy stellar mass. In our model, heavier NSCs at low redshift result from further evolution and mass growth through hierarchical mergers. The choice of threshold for major galaxy mergers was motivated by A. Lupi et al. (2014). We further choose a gas expulsion timescale of 100Myr, which, as K. Kritos et al. (2024a) showed, is long enough for Eddington-limited gas accretion to boost BH growth. The lifetime of accretion disks may be several tens of Myr (e.g., see M. Chiaberge et al. 2025). Moreover, we set = =f f 1g,max ,max 6 The Astrophysical Journal, 991:58 (18pp), 2025 September 20 Kritos et al. to allow for the possibility that a major fraction of the stellar mass is contained within the central collisional component of the galaxy and is gas-dominated, motivated by the recent little red dots observed by the JWST (J. Matthee et al. 2024) and the inside-out formation paradigm of galaxies (J. Silk et al. 2024). In Table 1, we summarize the set of fiducial hyperparameters used in the simulations. While throughout this section, we present a comparison of our results with observational data, we nevertheless have not attempted to fit the hyperparameters of our model to those data. We run the model for each merger tree once, using a different seed number. See the Appendix for the variations of our results on different realizations. 3.1. Local Massive BH Properties In this subsection, we discuss the properties of all massive BHs in our (16Mpc)3 comoving simulation box at redshift z = 0.25. We present, in order, the MBH–M� relation (Section 3.1.1), the massive BH spin (Section 3.1.2), the occupation fraction of massive BHs as a function of the galaxy’s stellar mass (Section 3.1.3), and finally, the number densities of IMBHs and SMBHs and their mass function (Section 3.1.4). 3.1.1. The MBH–M� Relation In Figure 5, we show all nuclear massive BHs in our simulation box in the MBH–M� plane, where M� is the stellar mass of the host galaxy (this is just the result of one realization of NSCs, because f� and fg are drawn randomly). We also include observational data (red points for detections and red triangles for limits) from J. E. Greene et al. (2020) as well as data from the catalog of dwarf active galactic nuclei of M. Mezcua et al. (2023; magenta squares with error bars). Notice that the sample is intrinsically biased: the heavier BHs have a higher probability of detection in dwarfs, and a direct comparison between our results and the magenta points requires a proper consideration of selection effects. Our results do not contain galaxies with M� < 107M⊙ because we set the critical stellar mass for NSC formation Mcr = 107M⊙ in the merger tree. We remind the reader that galaxy stellar masses M� are taken from the NEWHORIZON simulation. If we fit a linear model to our simulation data for ( )/ /M M M Mlog , log 10BH 10 using the linregress func- tion from SciPy, and we include only massive BHs with MBH > 100M⊙ (i.e., we restrict the fit to galaxies that contain a nuclear massive BH), we obtain (note that the logarithm is base 10) ( ) ( ) ( ) = ± + ± M M M M log 7.8 0.1 1.10 0.05 log 10 , 5 BH 10 with a Pearson correlation coefficient of ≃0.74, which implies a strong correlation. The cyan line in Figure 5 corresponds to the mean slope and intercept from this simple model. For comparison, orange lines are similar fits of the observational data from J. E. Greene et al. (2020) that either include upper limits (solid orange) or consider only estimated masses (dashed orange—i.e., in this case, we exclude systems for which we only have upper limits). On average, our model predicts larger BH masses relative to the observational data, although there is significant scatter in the data. It also predicts a shallower correlation, which roughly follows the MBH ∼ 10−2 M� relation. Thus, our population of massive BHs is heavier in the dwarf galaxy regime M� < 109M⊙ (bottom-heavy). Moreover, the majority of our galaxy sample are dwarfs with M� < 109M⊙. Our MBH–M� relation has more scatter in the IMBH mass range (102M⊙ < MBH < 106M⊙), and the MBH/M� ratio is in the range 10−3–10−1. This is consistent with overmassive BHs observed in nearby active dwarfs (M. Mezcua et al. 2023). The relation becomes tighter in the SMBH regime (MBH > 106M⊙) and for massive galaxies (M� > 109M⊙). For context, the Milky Way has M� ∼ 1011M⊙ with a central BH of ∼106M⊙, and it is therefore an atypical massive galaxy with a relatively small SMBH compared to the average. 3.1.2. Massive BH Spin In Figure 6, we show our simulated massive BH population in the mass–spin plane. We decompose our data set into nuclear BHs (gray filled circles), satellite BHs (gray plus signs), and ejected BHs (gray crosses). For comparison, we also overplot spin estimates (blue diamonds) and lower limits (red diamonds) made with the X-ray reflection technique, as listed in Table 1 of C. S. Reynolds (2021). Notice that these published measurements are only in the SMBH regime Table 1 Set of Fiducial Hyperparameter Values Mcr Qth τge fg,max f ,max (M⊙) (Myr) 107 1:10 100 1.0 1.0 106 107 108 109 1010 1011 1012 M?/M� 102 103 104 105 106 107 108 109 1010 1011 M B H /M � 100 10−2 10−4 galaxy, limits galaxy, detections dwarf AGN Figure 5. The mass of nuclear massive BHs and the corresponding stellar mass of their host galaxies at redshift z = 0.25 from our model (gray points). The red points with error bars show BH mass estimates in nearby galaxies, while the red triangles are upper limits on BH masses from the catalog of J. E. Greene et al. (2020). The magenta squares and associated error bars correspond to the catalog of dwarf active galactic nuclei from M. Mezcua et al. (2023). The cyan solid and orange solid (dashed) lines are power-law fits of the gray and red points including (excluding) upper limits. For reference, the thin lime lines correspond to different values of MBH/M�. 7 The Astrophysical Journal, 991:58 (18pp), 2025 September 20 Kritos et al. (MBH > 106M⊙), and they represent a biased population due to selection effects: highly spinning SMBHs are easier to measure with the X-ray reflection technique because they are more luminous. At the time of writing, no confident spin measurements of IMBH candidates have been published (J. E. Greene et al. 2020). We find that nuclear SMBHs (MBH > 106M⊙) are highly spinning (χBH > 0.8). Most of our nuclear SMBH population has χBH ∼ 1. This is because SMBHs grow primarily by accreting significant amounts of gas during major galaxy– galaxy mergers. At each gas accretion episode, we model the accretion phase as coherent, and we assume alignment between the BH spin and the angular momentum of the disk. If a BH coherently accretes an amount of gas at least of the order of its mass, then its spin asymptotes to its maximal value in a short finite time (J. M. Bardeen 1970). For context, the mass increases by one e-fold in about 50 Myr, assuming the Eddington accretion rate and a radiative efficiency of 10%. The spins of IMBHs are generally not as high. The distribution is bimodal, with one cluster of values around χBH ∼ 0 and another one at χBH ∼ 1. The spin of IMBHs depends on their formation channel and host galaxy merger history. If most of the IMBH’s mass is assembled through repeated BH mergers or consumption of stars, the spin asymptotes to zero (S. A. Hughes & R. D. Blandford 2003). This results from the asymmetry of the innermost stable circular radius for prograde versus retrograde orbits and the assumption of isotropy. On the other hand, coherent gas accretion events spin up the BH to unity as long as the BHs accrete an amount of gas ∼MBH (J. M. Bardeen 1970; E. Berti & M. Volonteri 2008; M. Dotti et al. 2013; Y. Dubois et al. 2014b; K. Kritos et al. 2024b; R. S. Beckmann et al. 2025). 3.1.3. Occupation Fraction In Figure 7, we show the occupation fraction of massive nuclear BHs with mass >103M⊙ (gray dashed), >106M⊙ (gray dashed–dotted), and >108M⊙ (gray dotted), as well as the occupation fraction fNSC of NSCs with half-mass radii <50 pc (thick gray solid line). We assume any NSC with a half-mass radius that exceeds 50 pc to have dissolved in the host galaxy and contributed to the bulge component. For comparison, we show limits on the BH occupation fraction in low-mass galaxies with masses 109–1010M⊙ from dynamical (D. D. Nguyen et al. 2019) and X-ray (B. P. Miller et al. 2015) surveys (see also Figure 5 in J. E. Greene et al. 2020). For reference, the blue line shows the nucleation fraction of galaxies in the Virgo cluster from the catalog of R. Sánchez- -Janssen et al. (2019), and the red line is the linear model for the occupation fraction of BHs with MBH > 105M⊙ from J. E. Greene et al. (2020). Our model predicts that ∼80% of the dwarf galaxies with stellar mass ≲109M⊙ contain IMBHs with mass >103M⊙. A fraction, <20%, of these dwarf galaxies contain SMBHs (MBH > 106M⊙). Massive BHs with MBH > 108M⊙ are hosted by galaxies with stellar mass >109M⊙. The occupation fraction converges to unity in massive galaxies with stellar mass ≳1010M⊙. For comparison, the value of the occupation fraction for BHs with MBH > 105M⊙ predicted by other studies is in the range 0.1–0.7 (0.1–1.0) for 108M⊙ (109M⊙) stellar-mass galaxies (see Table 1 of J. E. Greene et al. 2020). Our predictions tend to be on the higher end of these ranges. The NSC occupation fraction approaches unity for dwarfs with stellar mass below 108M⊙ and drops to zero in massive galaxies above 1010M⊙. In contrast, as shown by the right 102 103 104 105 106 107 108 109 1010 1011 MBH/M� 0.0 0.2 0.4 0.6 0.8 1.0 χ B H measurements limits Figure 6. The population of simulated massive BHs at redshift z = 0.25 in the mass–spin plane. The gray filled circles correspond to nuclear BHs, the plus signs to satellite BHs, and the crosses to ejected massive BHs. We also show measurements with a meaningful upper bound (blue diamonds) and lower limits on spin measurements (red diamonds) made with the X-ray reflection spectroscopy technique. All data are from C. S. Reynolds (2021). 107 108 109 1010 1011 M?/M� 0.0 0.2 0.4 0.6 0.8 1.0 O cc u p at io n fr ac ti on L in ea r V irgo N S C s X-ray Dynamical fNSC MBH/M� > 103M� 106M� 108M� Figure 7. Occupation fraction of massive BHs with mass >103 M⊙ (gray dashed), >106 M⊙ (gray dashed–dotted), and >108 M⊙ (gray dotted), as well as the NSC occupation fraction ( fNSC; thick gray solid). The linear model of J. E. Greene et al. (2020) for the MBH > 105 M⊙ occupation fraction and the NSC occupation fraction in the Virgo cluster from R. Sánchez-Janssen et al. (2019) are shown with red and blue lines, respectively. The cyan and yellow shaded regions correspond to massive BH occupation fraction limits from the dynamical and X-ray surveys in D. D. Nguyen et al. (2019) and B. P. Miller et al. (2015), respectively. 8 The Astrophysical Journal, 991:58 (18pp), 2025 September 20 Kritos et al. panel of Figure 3 based on observational data in N. Neumayer et al. (2020), the occupation declines for dwarfs. However, it should be noted that the results of the NSC occupation fraction presented in N. Neumayer et al. (2020) are lower limits, due to the difficulty of identifying NSCs. Our predictions for high occupation fractions are impacted by our choice to form an NSC in every galaxy with stellar mass larger than 107M⊙. We have found that the final stellar mass of a galaxy (and thus the mass of the central BH) depends on the number of merger nodes in the tree, and that these two quantities are correlated. Thus, massive galaxies have typically undergone tens of major merger episodes during their evolution, while some dwarf galaxies (M� < 1010M⊙) may have had no merger history. At each merger node, NSCs coalesce and the radius of the remnant NSC increases according to energy conservation (see Equation (4b)). This effect contributes to the expansion of the NSC, in addition to expansion due to internal heating determined by Hénon’s principle. Eventually, the NSCs of all massive galaxies in our sample have expanded beyond 50 pc; hence, the NSC occupation fraction drops for M� > 1010M⊙. Despite the high NSC occupation fraction in dwarf galaxies with no merger history, a few NSCs in our sample have still expanded beyond 50 pc (see the open circles in Figure 8). We identified these systems as containing maximally spinning ∼104M⊙ IMBHs that have grown by accretion in NSCs that had a seed mass of ∼107M⊙ and an initial half-mass radius of several parsecs. These systems expanded by a factor of ∼100 over the span of 10 Gyr. Observationally, these systems would show up as low-luminosity active galactic nuclei hosted in dwarfs with a diffuse bulge lacking an NSC component. 3.1.4. Number Density In Figure 9, we show the local massive BH mass function at z = 0.25. We decompose the total number density into satellite, nuclear, and ejected BH components. For compar- ison, we show the lower-limit estimates of the mass function in Figure 6 from J. E. Greene et al. (2020), which assumes two models for the occupation fraction: the “Linear” and “NSC” model. In total, the simulation box contains 522 (201) IMBHs (SMBHs) in the (16Mpc)3 comoving volume at z = 0.25, corresponding to a local number density of 0.13Mpc−3 (4.9 × 10−2 Mpc−3). Moreover, the number density of ejected IMBHs (SMBHs) is 1.1 × 10−2 Mpc−3 (1.4 × 10−2 Mpc−3). We find roughly equal contributions (≈6.4 × 10−2 Mpc−3) of nuclear and satellite IMBHs. For comparison, according to other models discussed in Section 2.4 of J. E. Greene et al. (2020), the number density of IMBHs in nuclei is in the range 0.02–0.25Mpc−3. This is consistent with our result for the density of nuclear IMBHs. On the other hand, according to other studies (see J. E. Greene et al. 2020), the number density of wandering IMBHs—what we term “off-nuclear,” which combines the satellite and ejected components—is predicted to be >0.3 Mpc−3. Consequently, our model predicts a factor of at least 4 lower number densities of off-nuclear IMBHs than other studies. The BH mass function turns over at MBH ∼ 104M⊙ and indicates that our NSC scenario produces a larger abundance of heavy seeds with masses ∼105M⊙. Most of these IMBHs are assembled in dwarf galaxies that do not have a rich merger history. The IMBHs stay in the center of their hosts until the present era. There is an underdensity of light seeds in the range 102–104M⊙, while the majority of BHs in the Universe are stellar-mass BHs, remnants of massive stars that did not grow. We summarize the number densities discussed in this section in Table 2. Notice that our IMBH densities have been derived assuming that they assemble solely in NSCs. In contrast, several IMBHs might have been formed in off-nuclear globular star clusters or other environments. Therefore, our IMBH number density estimates represent lower limits. Finally, given the small volume of NEWHORIZON, there are 107 108 109 1010 1011 M?/M� 100 101 102 103 r h /p c 100 101 102 Nnodes Figure 8. Our sample of half-mass radii of NSCs vs. the stellar mass of their host galaxy at redshift z = 0.25. The color denotes the number of nodes in each host’s merger tree, and the open circles correspond to systems with no merger history (Nnodes = 0). 102 103 104 105 106 107 108 109 1010 1011 MBH/M� 10−5 10−4 10−3 10−2 10−1 n B H /M p c− 3 Total Linear NSC Satellite Nuclear Ejected Figure 9. Local massive BH mass function, decomposed into satellite (dashed–dotted), nuclear (dashed), and ejected (dotted) BH components. The thick solid line corresponds to the total number density. The lower-limit estimates “Linear” and “NSC” from J. E. Greene et al. (2020) for the nuclear contribution of the mass function are shown in red and blue, respectively. 9 The Astrophysical Journal, 991:58 (18pp), 2025 September 20 Kritos et al. no massive galaxies; therefore, the number density is a lower limit at the high SMBH mass end. 3.2. Local NSC Properties In this subsection, we discuss the present-day properties of NSCs within our simulation box at redshift z = 0.25. We consider only NSCs with a final half-mass radius smaller than 50 pc, as we assume that sparser NSCs have dissolved and contribute to the bulge component of their host. In total, the (16Mpc)3 comoving box contains 1133 NSCs with rh < 50 pc, corresponding to a local NSC number density of roughly 0.28Mpc−3. This is comparable to (but a factor of 2 larger than) the predicted total local number density of massive BHs (MBH > 100M⊙). In the following, we present the MBH–MNSC relation (Section 3.2.1) as well as the rh–MNSC and MNSC–M� relations (Section 3.2.2) predicted by the model. Once again, the evolution of M� is derived from the NEWHORIZON simulation. 3.2.1. The MBH–MNSC Relation In Figure 10, we show the local MBH–MNSC correlation. For comparison, we include BH and NSC data from J. E. Greene et al. (2020) corresponding to either claimed detections or upper limits in the centers of galaxies and ultracompact dwarfs. Despite the large scatter in the data and in the simulation results, we observe a correlation between the BH mass and the NSC mass, because heavier clusters harbor heavier BHs. In particular, BHs with mass >106M⊙ and >108M⊙ are only hosted in NSCs with mass≳106M⊙ and≳107M⊙, respectively. The scatter increases in NSCs with mass ≲107M⊙, and the results span a wide range of values of MBH/MNSC ∈ [10−4, 10]. The observational data also suggest a positive correlation between BH and NSC mass, albeit with a shallower slope. However, due to selection biases, the comparison is only qualitative. The heaviest NSCs in our catalog (MNSC ≳ 108M⊙) contain overmassive SMBHs, in the sense that the BH mass is larger than the NSC mass, in some cases up to a factor of ∼10. Most of these NSCs have half-mass radii larger than 50 pc, which we define as the threshold for NSC detectability. A special case is NGC 1023, which may contain an SMBH whose mass exceeds the host NSC mass by a factor of ∼10 or more. However the error bar on MBH is large, and the MBH/MNSC ratio could be as high as ∼1. Despite the low statistics and observational challenges at the high NSC mass above 108M⊙, we find that our model does not seem to populate lower values of MBH/MNSC. This happens because the gas mass replenished at each merger episode is proportional to the NSC mass. The fact that this assumption is unable to predict lower values of MBH/MNSC suggests that replenishing gas events for high- mass NSCs should become less efficient. 3.2.2. The rh–MNSC and MNSC–M� Relations In Figure 11, we show the structural properties of our NSC population in the local Universe. We mark with open gray circles all NSCs with a half-mass radius that has expanded past 50 pc. We further compare our results against the NSC catalog of I. Y. Georgiev et al. (2016) and the population of Milky Way globular clusters (in the left panel). Our NSC population is heavier than the sample from I. Y. Georgiev et al. (2016). Moreover, NSCs are heavier than globular clusters by at least an order of magnitude (on average) but have similar half-mass radii, as confirmed by observations (N. Neumayer et al. 2020). We notice that our NSC population is more massive in the dwarf galaxy regime, with NSC masses roughly uniformly populated in the range 106–107M⊙ (see right panel of Figure 11). This is related to our choice of initial cluster mass being a random fraction of Mcr = 107M⊙. Points with MNSC > 107M⊙ have grown their mass through subsequent mergers. One of the assumptions that lead to overmassive NSCs is the NSC mass conservation during NSC–NSC mergers, which suggests that NSCs should lose a nonnegligi- ble fraction of their mass during major merger episodes. 3.3. Transient Events In this subsection, we compute the merger rates of massive BH mergers (Section 3.3.1), stellar-mass BH–BH captures in NSCs (Section 3.3.2), extreme mass ratio inspirals (Section 3.3.3), and TDEs (Section 3.3.4). 3.3.1. Massive BH Binary Merger Rate Massive BH binaries form when galaxies merge with a minor-to-major galaxy mass ratio that does not exceed a certain threshold, for which we take the fiducial value Table 2 Predicted Number Densities of Massive BHs (IMBHs and SMBHs) at z = 0.25 Massive BH nIMBH nSMBH Type (Mpc−3) (Mpc−3) Nuclear 0.064 0.032 Satellite 0.064 0.017 Ejected 0.011 0.014 Total 0.130 0.049 105 106 107 108 109 MNSC/M� 102 103 104 105 106 107 108 109 1010 1011 M B H /M � 101 100 10−1 10−2 NGC 1023 galaxy, limits UCD, limits galaxy, detections UCD, detections Figure 10. The BH mass vs. the total mass of the host NSC for the population of massive BHs at redshift z = 0.25. The open circles correspond to NSCs whose half-mass radius has expanded past 50 pc. We also plot observational data from J. E. Greene et al. (2020), including claimed detections (points with error bars) and upper mass limits (triangles) in galactic centers (red) and ultracompact dwarfs (UCDs; blue). For reference, the thin lime lines show different values of MBH/MNSC. 10 The Astrophysical Journal, 991:58 (18pp), 2025 September 20 Kritos et al. Qth = 0.1. Thus, the massive BH binary merger rate in our model follows the major galaxy–galaxy merger rate. In Figure 12, we show the source-frame massive BH binary merger rate as a function of redshift, ( )R zm , as well as its decomposition into the three subchannels: IMBH–IMBH, IMBH–SMBH, and SMBH–SMBH. The inset shows the observer-frame cumulative merger rate up to redshift z, computed using Equation (17) from C. L. Rodriguez et al. (2016), ( ) ( ) ( ) ( )< = + RR z z dV z dz dz z1 , 6 z c 0 m com where dVcom(z) is the all-sky differential comoving volume element at redshift z. We find that the total merger rate attains its maximum at z ≃ 3, with a peak value of ≃1.4 × 10−2 yr−1 Gpc−3. The corresponding values for each subchannel are ≃0.8 × 10−2 yr−1 Gpc−3 at z ≃ 4.0 for IMBH– IMBH, ≃0.4 × 10−2 yr−1 Gpc−3 at z ≃ 3.5 for IMBH–SMBH, and ≃0.5 × 10−2 yr−1 Gpc−3 at z ≃ 2.5 for SMBH–SMBH mergers. Hence, the total source-frame mass decreases with redshift, as expected from a hierarchically merging massive BH population as larger structures assemble with time. In the local Universe, the SMBH–SMBH merger rate is predicted to be of the order of a few× 10−2 yr−1 Gpc−3. The observer- frame merger rate out to z = 5 is ≃2.7 yr−1 for IMBH–IMBH, ≃1.5 yr−1 for IMBH–SMBH, and ≃1.1 yr−1 for SMBH– SMBH mergers, adding up to a total of ≃5.3 yr−1. The merger rate decreases for redshifts above z = 4 because there are very few sufficiently massive galaxies with large BHs present in the simulation box at those redshifts. In our model, we set the critical mass for NSC formation (and, consequently, for massive BH formation) at 107M⊙. While many major galaxy– galaxy mergers occur beyond z = 4, the merging galaxies have not yet developed massive BHs. If we consider that massive BHs could also form through alternative scenarios unrelated to their seeding in an NSC environment (e.g., as direct-collapse BHs), then our merger rate at these high redshifts may be underestimated. In Figure 13, we show the properties of our massive BH binary merger catalog. These GW sources are observational targets for spaceborne GW observatories such as LISA (P. Amaro-Seoane et al. 2017). In total, the catalog contains 192 events. For each event with a primary-to-secondary mass 104 105 106 107 108 109 MNSC/M� 10−1 100 101 102 103 r h /p c Globular clusters Georgiev+2016 107 108 109 1010 1011 1012 M?/M� 104 105 106 107 108 109 M N S C /M � Georgiev+2016 Figure 11. Left: present-day structural properties of NSCs (gray points) in the half-mass–radius (rh) vs. NSC mass (MNSC) plane. Also shown are Milky Way globular clusters (blue plus signs) and the NSC catalog of I. Y. Georgiev et al. (2016; red points with error bars). Right: same as the left panel but in the galaxy stellar mass (M�) vs. NSC mass (MNSC) plane. 0 1 2 3 4 5 6 z 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 R m × 10 2 (y r− 1 G p c− 3 ) Total IMBH-IMBH IMBH-SMBH SMBH-SMBH 0 2 4 6 0 2 4 6 Rc( 10 out to redshift z = 6. The loudest events have total source-frame masses in the range 106–107M⊙, roughly equal-mass components (q ≃ 1), and z ≲ 2. 3.3.2. Stellar-mass BH–BH Capture Rate Interactions among stellar-mass BHs around massive BHs in the cores of NSCs lead to the formation of BH binaries that merge. Here, we compute the rate of BH–BH captures through the emission of GWs in NSCs. This rate corresponds to the merger rate, because captured pairs merge promptly, on a timescale smaller than the disruption timescale, as long as the density is ( )/v10 pc 100 km s12 3 rel 1 6 (R. M. O’Leary et al. 2009). The two-body GW capture cross section, Σcap, is given by Equation (4) from H. Mouri & Y. Taniguchi (2002), where we set m1 = m2 = 10M⊙ and =v v2rel BH (where vBH is the three-dimensional velocity dispersion in the BH subsystem). Assuming for simplicity uniform conditions within the half- mass radius of the BH subsystem, rBH, which is related to vBH through the virial theorem, we compute the capture rate per NSC as ( ) / × N r v N v 4 2.6 10 yr 1000 100 km s . 7 cap BH 2 BH 3 cap rel 5 1 BH BH 1 31 7 Here NBH is the number of BHs in the BH subsystem. As the NSC expands, the velocity dispersion drops, and the capture rate declines with time. To compute the total source-frame capture rate density at redshift z, we add up the contributions from each NSC within the simulated comoving volume and divide by that volume, ( ) ( ) ( ) ( )=R z z 1 16 Mpc , 8 i icap 3 cap, where the index i runs over all NSCs for which the time- integrated capture rate does not exceed NBH/2. In Figure 14, we plot the source-frame rate density of BH– BH captures as a function of redshift (solid blue), as well as the observed intrinsic cumulative merger rate (left inset) and number density of contributing NSCs, nNSC(z) (right inset). Due to finite number statistics, the estimated rate density shows fluctuations; nevertheless, on average, the source-frame volumetric capture rate is ≈10 yr−1 Gpc−3 out to redshift z ≃ 5. The rate decreases at z > 5 due to a declining number density of NSCs. Predictions of the local BH–BH merger rate from NSCs in other studies range from ∼10−2 up to a few tens of events yr−1 Gpc−3 (see Figure 3 in I. Mandel & F. S. Broekgaarden 2022 and the blue vertical band in Figure 14), while the inferred merger rate from GW observations during O3 is in the range 18–44 yr−1 Gpc−3 (R. Abbott et al. 2023; see the black error bar “LVK” in Figure 14). We compute the cumulative merger rate Rc from Equation (6) assuming the average source-frame merger rate density. Within z = 5, we find the cumulative observed intrinsic rate of BH–BH captures to be ≃6000 yr−1. For comparison, it is expected that several hundred thousand 102 103 104 105 106 107 108 109 1010 1011 MMBH,tot/M� 0 1 2 3 4 5 6 z q = 1 q = 10 q = 100 q > 1000 100 101 102 103 104 105 LISA SNR Figure 13.Massive BH binary mergers in our generated catalog in the plane of merger redshift z and total source-frame binary mass MMBH,tot. The color shows the SNR of each merger event in the LISA detector. The point size represents q , where q = m1/m2 > 1 is the binary mass ratio. Benchmark point sizes for q = [1, 10, 100] are given in the legend, and crosses correspond to merger events for which q > 1000, such that our waveform models cannot be trusted for the SNR calculation. 0 1 2 3 4 5 6 7 z 10−5 10−4 10−3 10−2 10−1 100 101 102 103 104 105 106 R (y r− 1 G p c− 3 ) LVK Rcap RTDE REMRI 0.001 0.1 10 100 102 104 106 108 Rc(104M⊙ (see the green line in Figure 14). We find an intrinsic merger rate density comparable to the capture rate between stellar-mass BHs and a cumulative rate of ∼4000 inspirals per year. For comparison, we show the range of predictions in Table I from S. Babak et al. (2017; green band in the left inset of Figure 14), where extreme mass ratio inspiral rates lie in the range 13–20,000 yr−1, depending on model assumptions. 3.3.4. TDE Rate Stars approaching BHs within the tidal radius rt may be tidally disrupted. The condition for generating an observable electromagnetic counterpart is that rs < rt < r�, where rs is the BH Schwarzschild radius and r� is the star’s physical radius, which constrains the BH mass to be in the range 105M⊙ ≲ MBH ≲ 108M⊙ (S. Gezari 2021). If these conditions are met, the TDE cross section in the gravitational focusing regime is given by ( )r r c v , 9s tTDE 2 where v� is the three-dimensional stellar velocity dispersion. Using the half-mass number density of stars, the TDE rate per NSC in the full loss cone is computed as (H. Pfister et al. 2020) ( ) / × × N r v N M M v 4 3 10 yr 10 10 100 km s . 10 TDE h 3 TDE 4 1 7 2 BH 6 4 3 1 5 Again, the half-mass radius rh is related to v� through the virial theorem. The source-frame volumetric TDE rate density ( )R zTDE and the observer-frame cumulative rate are computed as in Equations (6) and (8). The TDE rate density is in the range 104–105 yr−1 Gpc−3 out to z = 5, about 3–4 orders of magnitude above the BH–BH capture rate density. Based on the left inset of Figure 14, we expect thousands of TDEs and a single BH–BH capture in an NSC at a luminosity distance of 500Mpc (z ≃ 0.1). 4. Caveats In this section, we list the main limitations of our study (Section 4.1), and we qualitatively discuss the dependence of our results on the choice of our model hyperparameters (Section 4.2). 4.1. Limitations Below, we list how the assumptions made in our model may affect our results. The NSC initial conditions presented in Section 2.2 were arbitrarily chosen. While a few strongly magnified star clusters have been identified with the JWST (L. Mowla et al. 2022; E. Vanzella et al. 2023; B. Welch et al. 2023; A. Adamo et al. 2024), there is no consensus on the initial conditions for high- redshift NSCs. Compact NSCs on subparsec scales could theoretically form in the high-redshift Universe (z > 10; B. Devecchi et al. 2010; T. Kimm et al. 2016; L. Mayer et al. 2025). This assumption is essential for assembling SMBH seeds in the NSC paradigm (see Figure 3). We consider a single burst of NSC formation when the stellar mass of each galaxy reaches the critical value of Mcr = 107M⊙. This assumption may lead to an overabundance of NSCs in the local Universe if not every galaxy forms an NSC in its core. While we consider the in situ formation of NSCs in protogalaxies, we note that globular clusters that sink to the center due to dynamical friction and form near the nuclear region of the host galaxy may contribute to the mass assembly of the NSC (see, e.g., O. Y. Gnedin et al. 2014). In our model, NSCs lose mass as they evolve and gain mass through NSC mergers and through gas inflows during major mergers. We simplify our analysis by neglecting the contrib- ution from this ex situ channel of inspiralling globular clusters. Furthermore, NSCs have been observed with multiple stellar populations (N. Neumayer et al. 2020), indicating complex star formation histories rather than a monolithic formation at some high redshift. The larger number of NSC mergers also contributes to the presence of overmassive BHs in our model, because BHs grow by mergers and gas accretion. The off-nuclear number density we estimated in Section 3.1.4 is likely a lower bound, because some off-NSCs could form IMBHs in their centers. If some of these IMBHs are ejected from their cores (as suggested by some simulations, especially for IMBHs with low mass, MBH ≲ 1000M⊙; K. Holley-Bockelmann et al. 2008; E. G. Prieto et al. 2022), our estimated number density of ejected massive BHs is also a lower bound. We have not included the time delay in massive BH binary coalescence following major galaxy mergers (but see K. Li et al. 2024 for subgrid modeling). This time may be decomposed into three contributions: the dynamical friction timescale (τdf; from galaxy merger to formation of the tight massive BH binary), the binary hardening timescale (τhar; from binary formation until its evolution becomes dominated by GW emission), and the GW merger timescale (τgw), so that τdelay = τdf + τhar + τgw (see, e.g., Equation (6) from V. Langen et al. 2025). We assume instantaneous binary formation in the postmerger galaxy by ignoring τdf. Further- more, we do not consider cases where a third massive BH is 13 The Astrophysical Journal, 991:58 (18pp), 2025 September 20 Kritos et al. brought into the vicinity of the evolving massive BH pair following a third subsequent galaxy–galaxy merger episode. Assuming that the final parsec problem is not, in fact, a problem (F. M. Khan et al. 2013), we can estimate the time delay from binary formation until binary merger as (see Equations (30) and (31) from G. D. Quinlan 1996) ( ) / / / / + × H M c M M 44 Myr 100 km s 16 10 p 10 . 11 har gw 1 4 5 4 5 5 3 4 5 BH 6 3 5 This timescale can thus be much smaller than the timescale between successive galaxy mergers, and it can be neglected. Our simulations are built on merger tree outputs from the NEWHORIZON results. Hence, no feedback effect from the massive BH or the NSC on the merger tree has been implemented. In our model, the parameters f� and fg are both independent and randomly sampled from the ranges [ ]f0, ,max and [ ]f0, g,max , respectively (see Section 2.2). According to our criterion, the parameter Mcr represents the critical stellar mass that the host galaxy must reach for NSC formation to begin in the center. In some cases, it is assumed that millions of solar masses of residual gas can accumulate in the center without fragmenting into stars at the time NSC formation occurs, particularly when fg is close to unity, indicating low star formation efficiency. It is important to note that the condition Mcl,0 ≫ Mcr is not contradictory. This is becauseMcl,0 includes both stellar and gas contributions, while Mcr is a threshold specifically for the critical stellar mass of the galaxy. This situation indeed facilitates the formation of early massive BHs. We found that accumulating gas mass in the center is essential for the rapid growth of BH seeds throughout cosmic time. However, the ratio MBH/M� never exceeds unity; thus, we do not observe unrealistically overmassive BHs. Moreover, M. Mezcua et al. (2023) identified an overmassive population of BHs in nearby dwarf galaxies, and our findings are not in substantial conflict with these observations, as illustrated in Figure 5. Our population of SMBHs is highly spinning. In principle, the presence of strong magnetic fields in the vicinity of SMBHs may efficiently extract angular momentum from the BH and slow it down through the Blandford–Znajek process (R. D. Blandford & R. L. Znajek 1977). We have not simulated this process. The highly spinning SMBHs in our simulations are consistent with our assumed coherent gas accretion episodes. On the other hand, if the accretion process is chaotic, then the spin asymptotes to zero (A. R. King et al. 2008). We have limited the accretion rate to the Eddington limit. A subpopulation of super-Eddington accreting SMBHs has been identified at large distances (J. E. Greene et al. 2024; R. Mai- olino et al. 2024a; H. Suh et al. 2025). Moreover, the observation of a dormant SMBH at redshift z = 6.68 suggests that massive BHs are likely assembled through a series of super-Eddington burst episodes (I. Juodžbalis et al. 2024). Our analysis is limited to redshifts below z = 6 and to a typical comoving volume of the Universe, in the sense that we ignore rare peaks in the density field, which may give rise to ultraluminous quasars (E. O. Nadler et al. 2023). Last but not least, we have not included selection effects when comparing with observational data. For this reason, all of our comparisons with observations should be regarded as qualitative. It is for this reason we have not attempted to fit the hyperparameters of our model, as that would require careful inclusion of observational biases. 4.2. Dependence on Hyperparameters We do not attempt to perform a comprehensive study of the dependence of our predictions on the model hyperparameters, as this would go beyond the scope of the present work. Therefore, this section focuses on a qualitative comparison between the fiducial results and the effect of varying the main model hyperparameters, one at a time. 4.2.1. Lower Mass Threshold for NSC Formation We have repeated simulations with Mcr = 106M⊙, i.e., an order of magnitude smaller than the fiducial value used in our model, fixing all other hyperparameters to their fiducial values. In this case, the NSCs are typically much lighter (by about an order of magnitude), and not heavier than 108M⊙. As a consequence, massive BH masses are also smaller because the masses of BH seeds strongly depend on the initial conditions of NSCs (see Figure 4). Moreover, at each galaxy merger episode, an amount of gas proportional to the NSC mass is added to the system (see Section 2.4). Since, in this case, NSCs are lighter than in the results presented in Section 3, BHs accrete less gas as they evolve in the merger tree. On the other hand, using the fiducial hyperparameters, we have shown that SMBHs can become overmassive relative to the observed local MBH–M� relation. 4.2.2. Lower Mass Ratio Threshold for Major Galaxy Mergers We have checked the effect of setting the mass ratio threshold to be Qth = 0.01, instead of the fiducial value of Qth = 0.1 used in Section 3. In this case, the merger rate increases, resulting in more ejections and a smaller density of massive satellite BHs. Nevertheless, this parameter is not critical, in the sense that it does not strongly affect the other observed correlations. During major gas inflows into the NSC, BHs accrete nearly at the Eddington rate. We have computed the ratio between the timescale of the accretion phase (set to 300Myr) and the time between successive major merger episodes, also known as the duty cycle. We find that the duty cycle always nearly peaks around 5% for both choices of Qth we have simulated. Finally, the total number of major mergers occurring in the simulated comoving box is 160 and 262 when we set Qth = 0.1 and 0.01. Such an increase by a factor of ∼2 does not lead to a significant gas supply increase; thus, the emergent BH population is not significantly altered. 4.2.3. Smaller Gas Expulsion Timescale We have carried out a simulation in which we have set τge = 10Myr while fixing all other hyperparameters to their fiducial values. Lowering the gas expulsion timescale by an order of magnitude limits the amount of gas accretion onto massive BHs—their dominant growth channel. Consequently, we find a larger number density of IMBHs compared with the fiducial results, and we find only a few SMBHs in the 14 The Astrophysical Journal, 991:58 (18pp), 2025 September 20 Kritos et al. (16Mpc)3 comoving volume. In particular, the heaviest BH in this volume is found to be ∼108M⊙. In addition, the occupation fraction of SMBHs drops, the slope of the MBH–M� relation decreases, and we end up underestimating the masses of massive BHs. 4.2.4. =f 0.1g,max We have repeated our simulations with a lower value of =f 0.1g,max . This choice limits the amount of gas in NSCs, so we find a larger population of slowly spinning BHs. Massive BHs have lower masses (below 108M⊙) compared to our fiducial model. This is because BHs have less gas to accrete from and therefore they do not grow as much from their seed mass value; it is not due to a lower density of NSCs that can produce heavy seeds. Moreover, the MBH–M� relation in this case is flatter, with a slope of ≃0.35. 4.2.5. =f 0.1,max In the case where we set =f 0.1,max , we limit the initial mass of the formed NSC relative to the host galaxy to be at most 106M⊙. We find anMBH–M� correlation that agrees more closely with the observed correlation, with a slope of ≃1.13. In this case, massive BHs are lighter because the seed masses are smaller, not because there is not enough gas accretion. 5. Conclusions SMBHs are prevalent at the centers of massive galaxies, and their masses scale with galaxy properties (including the velocity dispersion and mass of old stars and the dark halo mass). Increasing evidence suggests that these trends continue to low stellar masses and toward high redshift. Seeds are needed for SMBH formation, and this requirement is especially compelling at the highest redshifts explored by JWST. Our study of the hierarchical merging of galaxies over cosmic time, a central tenet of current cosmology, suggests that the central SMBHs found in nearby dwarf spheroidal galaxies may provide evidence for surviving IMBHs that once were prevalent as SMBH seeds at early epochs. Our key hypothesis is that such seeds formed in NSCs at early epochs, and we explore this model using cosmological merger trees. Our model succeeds despite adding ejections. A diverse population of SMBHs emerges, with masses up to (rarely) ∼1011M⊙. Even if a galaxy loses an SMBH due to a GW kick, a subsequent major merger event with another galaxy usually brings another SMBH into the center, maintaining the high occupation fraction. Gas is essential for forming massive BH seeds and for their subsequent growth. Inefficient star-forming NSCs in which the gas is ejected within a few tens of Myr do not contribute much to forming massive BH seeds. On the other hand, massive BH seed formation is efficient if the gas is retained longer in the cluster’s core, in the sense that in a large region of the parameter space, seeds form with a mass of ∼105 M⊙, independently of the star formation efficiency. In total, with our fiducial choice of hyperparameters, we find a massive BH mass of ≈1.7 × 1011M⊙ (not including stellar-mass BHs) within the (16Mpc)3 cosmological box. This translates to ≈4.2 × 1016M⊙Gpc −3, which is a factor of ≈10−4 less than the critical density to close the Universe. For comparison, according to A. Soltan (1982), the total mass accumulated by quasars is at least 4.7 × 1013M⊙Gpc −3, where the accretion efficiency is set to 10%. Notice that a galaxy that has lost its nuclear SMBH may regain it if it merges with another galaxy that hosts a central SMBH, provided that it has a rich merger history. Therefore, galaxies that have undergone several major merger episodes will likely contain nuclear massive BHs and a few wandering off-nuclear massive BHs. The off-nuclear BHs could either be surrounded by a stellar system (classified as an ultracompact dwarf), indicating a minor merger episode, or completely isolated, if they received a large enough GW kick. Such isolated SMBHs could only be identified either by their microlensing effect on foreground stars or by their electro- magnetic signatures when accreting from the interstellar medium. Our model allows us to infer a population of massive NSCs at high redshift and predicts that NSCs should generally contain massive central BHs. However, especially the most massive SMBHs are not invariably associated with NSCs, as the latter can be dynamically disrupted. Our code has many free parameters that we have chosen to give a crude fit to local scaling laws involving NSCs. We caution the reader that this is in the same spirit as other discussions of semianalytical models of galaxy formation, e.g., E. O. Nadler et al. (2023). Independent studies using different prescriptions are needed to test the robustness of our results. Our code and data set are publicly available on Zenodo (K. Kritos 2025). We anticipate future applications that will include the use of faster codes capable of generating merger tree histories for rare halos hosting extremely luminous high-redshift galaxies (as seen by JWST), follow cosmological simulations of halo assembly history from z > 10 down to mass scales of ∼105M⊙, and incorporate dynamical modeling of NSCs. We have also predicted GW merger rates, including massive BH binaries, extreme mass ratio inspirals, and stellar-mass BH binaries. Future experiments such as LISA, the Einstein Telescope, and Cosmic Explorer could detect these sources, either as a stochastic GW background or through the observation of individual merger events. Acknowledgments We thank Francesco Iacovelli and Luca Reali for discus- sions. K.K., E.B., and S.Y. are supported by NSF grant Nos. AST-2307146, PHY2207502, PHY-090003, and PHY-20043; NASA grant Nos. 20-LPS20-0011 and 21-ATP21-0010; the John Templeton Foundation grant 62840; the Simons Founda- tion; and the Italian Ministry of Foreign Affairs and International Cooperation grant No. PGR01167. K.K. is supported by the Onassis Foundation—Scholarship ID: F ZT 041-1/2023-2024. S.Y. is supported by the NSF Graduate Research Fellowship Program under grant No. DGE2139757. R.B. acknowledges support from a UKRI Future Leaders Fellowship (grant code: MR/Y015517/1). This work was carried out at the Advanced Research Computing at Hopkins (ARCH) core facility (https://www.arch.jhu.edu/), which is supported by the NSF grant No. OAC-1920103. For the purpose of open access, the author has applied a creative commons attribution (CC BY) license to any author accepted manuscript version arising from this submission. 15 The Astrophysical Journal, 991:58 (18pp), 2025 September 20 Kritos et al. http://www.arch.jhu.edu/ Appendix In this appendix, we explore the variations of our results on different realizations, assuming the fiducial hyper- parameters of Table 1. We show the variation of the local BH mass and spin distribution in Figure 15, NSC masses and radii in Figure 16, and the primary masses and mass ratios of all massive binary BH mergers in Figure 17. Note that the variation is the largest in the properties of ejected BHs, due to the small number of events per simulation. 100 101 102 103 104 105 106 107 108 109 1010 1011 1012 MBH/M� 100 101 102 C ou n t Nuclear 100 101 102 103 104 105 106 107 108 109 1010 1011 1012 MBH/M� 100 101 C ou n t Satellite 100 101 102 103 104 105 106 107 108 109 1010 1011 1012 MBH/M� 100 101 C ou n t Ejected 0.0 0.2 0.4 0.6 0.8 1.0 χBH 100 101 102 C ou n t Nuclear 0.0 0.2 0.4 0.6 0.8 1.0 χBH 100 101 102 C ou n t Satellite 0.0 0.2 0.4 0.6 0.8 1.0 χBH 100 101 C ou n t Ejected Figure 15. Eleven different realizations of the BH mass (top panels) and spin (lower panels) distributions at z = 0.25. The thick red histograms correspond to the data set in Section 3. The three columns correspond to nuclear (left), satellite (middle), and ejected (right) BH components. 104 105 106 107 108 109 MNSC/M� 100 101 102 C ou n t 10−1 100 101 102 103 104 rh/pc 100 101 102 C ou n t Figure 16. The same 11 realizations as in Figure 15 for the NSC masses (left) and half-mass radii (right). 16 The Astrophysical Journal, 991:58 (18pp), 2025 September 20 Kritos et al. ORCID iDs Konstantinos Kritosaa https://orcid.org/0000-0002-0212-3472 Ricarda S. Beckmannaa https://orcid.org/0000-0002-2850-0192 Joseph Silkaa https://orcid.org/0000-0002-1566-8148 Emanuele Bertiaa https://orcid.org/0000-0003-0751-5130 Sophia Yiaa https://orcid.org/0000-0002-9104-1734 Marta Volonteriaa https://orcid.org/0000-0002-3216-1322 Yohan Duboisaa https://orcid.org/0000-0003-0225-6387 Julien Devriendtaa https://orcid.org/0000-0002-8140-0422 References Abbott, R., Abbott, T. D., Acernese, F., et al. 2023, PhRvX, 13, 011048 Adamo, A., Bradley, L. D., Vanzella, E., et al. 2024, Natur, 632, 513 Alexander, T., & Natarajan, P. 2014, Sci, 345, 1330 Amaro-Seoane, P., Audley, H., Babak, S., et al. 2017, arXiv, arXiv:1702. 00786 Antonini, F., Gieles, M., & Gualandris, A. 2019, MNRAS, 486, 5008 Aubert, D., Pichon, C., & Colombi, S. 2004, MNRAS, 352, 376 Babak, S., Gair, J., Sesana, A., et al. 2017, PhRvD, 95, 103012 Baibhav, V., Berti, E., Gerosa, D., et al. 2019, PhRvD, 100, 064060 Baibhav, V., Berti, E., Gerosa, D., Mould, M., & Wong, K. W. K. 2021, PhRvD, 104, 084002 Bardeen, J. M. 1970, Natur, 226, 64 Baumgardt, H., & Kroupa, P. 2007, MNRAS, 380, 1589 Beckmann, R. S., Dubois, Y., Volonteri, M., et al. 2023, MNRAS, 523, 5610 Beckmann, R. S., Dubois, Y., Volonteri, M., et al. 2025, MNRAS, 536, 1838 Berti, E., & Volonteri, M. 2008, ApJ, 684, 822 Bhowmick, A. K., Blecha, L., Torrey, P., et al. 2024a, MNRAS, 529, 3768 Bhowmick, A. K., Blecha, L., Torrey, P., et al. 2025, MNRAS, 538, 518 Blandford, R. D., & Znajek, R. L. 1977, MNRAS, 179, 433 Bondi, H. 1952, MNRAS, 112, 195 Bonoli, S., Marulli, F., Springel, V., et al. 2009, MNRAS, 396, 423 Breen, P. G., & Heggie, D. C. 2013, MNRAS, 432, 2779 Cammelli, V., Monaco, P., Tan, J. C., et al. 2025, MNRAS, 536, 851 Campanelli, M., Lousto, C. O., Zlochower, Y., & Merritt, D. 2007, ApJL, 659, L5 Chabrier, G. 2005, Astrophys. Space Sci. Libr., 327, 41 Chiaberge, M., Morishita, T., Boschini, M., et al. 2025, arXiv:2501.18730 Dekel, A., Stone, N. C., Dutta Chowdhury, D., et al. 2025, A&A, 695, A97 Devecchi, B., Volonteri, M., Colpi, M., & Haardt, F. 2010, MNRAS, 409, 1057 Devecchi, B., Volonteri, M., Rossi, E. M., Colpi, M., & Portegies Zwart, S. 2012, MNRAS, 421, 1465 Di Matteo, T., Colberg, J., Springel, V., Hernquist, L., & Sijacki, D. 2008, ApJ, 676, 33 Dong-Páez, C. A., Volonteri, M., Beckmann, R. S., et al. 2023a, A&A, 676, A2 Dong-Páez, C. A., Volonteri, M., Beckmann, R. S., et al. 2023b, A&A, 673, A120 Dong-Páez, C. A., Volonteri, M., Dubois, Y., Beckmann, R. S., & Trebitsch, M. 2025, A&A, 695, A231 Dotti, M., Colpi, M., Pallini, S., Perego, A., & Volonteri, M. 2013, ApJ, 762, 68 Dubois, Y., Beckmann, R., Bournaud, F., et al. 2021, A&A, 651, A109 Dubois, Y., Volonteri, M., & Silk, J. 2014a, MNRAS, 440, 1590 Dubois, Y., Volonteri, M., Silk, J., Devriendt, J., & Slyz, A. 2014b, MNRAS, 440, 2333 Ellis, J., Fairbairn, M., Urrutia, J., & Vaskonen, V. 2024, arXiv:2410.24224 Freese, K., Ilie, C., Spolyar, D., Valluri, M., & Bodenheimer, P. 2010, ApJ, 716, 1397 Fujii, M. S., Wang, L., Tanikawa, A., Hirai, Y., & Saitoh, T. R. 2024, Sci, 384, adi4211 Fumagalli, G., Romero-Shaw, I., Gerosa, D., et al. 2024, PhRvD, 110, 063012 García-Quirós, C., Colleoni, M., Husa, S., et al. 2020, PhRvD, 102, 064002 Georgiev, I. Y., Böker, T., Leigh, N., Lützgendorf, N., & Neumayer, N. 2016, MNRAS, 457, 2122 Gerosa, D., & Kesden, M. 2016, PhRvD, 93, 124066 Gezari, S. 2021, ARA&A, 59, 21 Gnedin, O. Y., Ostriker, J. P., & Tremaine, S. 2014, ApJ, 785, 71 Gonzalez, J. A., Hannam, M. D., Sperhake, U., Bruegmann, B., & Husa, S. 2007, PhRvL, 98, 231101 Greene, J. E., Labbe, I., Goulding, A. D., et al. 2024, ApJ, 964, 39 Greene, J. E., Strader, J., & Ho, L. C. 2020, ARA&A, 58, 257 Haardt, F., & Madau, P. 1996, ApJ, 461, 20 Holley-Bockelmann, K., Gultekin, K., Shoemaker, D., & Yunes, N. 2008, ApJ, 686, 829 Hooper, D., Ireland, A., Krnjaic, G., & Stebbins, A. 2024, JCAP, 04, 021 Hughes, S. A., & Blandford, R. D. 2003, ApJL, 585, L101 Izquierdo-Villalba, D., Sesana, A., Bonoli, S., & Colpi, M. 2021, MNRAS, 509, 3488 Juodžbalis, I., Maiolino, R., Baker, W. M., et al. 2024, Natur, 636, 594 Khan, F. M., Holley-Bockelmann, K., Berczik, P., & Just, A. 2013, ApJ, 773, 100 Kimm, T., Cen, R., Devriendt, J., Dubois, Y., & Slyz, A. 2015, MNRAS, 451, 2900 Kimm, T., Cen, R., Rosdahl, J., & Yi, S. K. 2016, ApJ, 823, 52 Kimm, T., Katz, H., Haehnelt, M., et al. 2017, MNRAS, 466, 4826 King, A. R., Pringle, J. E., & Hofmann, J. A. 2008, MNRAS, 385, 1621 Komatsu, E., Smith, K. M., Dunkley, J., et al. 2011, ApJS, 192, 18 Kovács, O. E., Bogdán, Á., Natarajan, P., et al. 2024, ApJL, 965, L21 Kritos, K. 2025, Supermassive Black Hole Growth in Hierarchically Merging Nuclear Star Clusters: Data and Code Release, Version v1, Zenodo, doi:10. 5281/zenodo.15776363 Kritos, K., Berti, E., & Silk, J. 2023, PhRvD, 108, 083012 Kritos, K., Berti, E., & Silk, J. 2024a, MNRAS, 531, 133 Kritos, K., Reali, L., Gerosa, D., & Berti, E. 2024b, PhRvD, 110, 123017 Kroupa, P. 2002, Sci, 295, 82 Kroupa, P., Subr, L., Jerabkova, T., & Wang, L. 2020, MNRAS, 498, 5652 Langen, V., Tamanini, N., Marsat, S., & Bortolas, E. 2025, MNRAS, 536, 3366 Li, K., Volonteri, M., Dubois, Y., Beckmann, R., & Trebitsch, M. 2024, arXiv:2410.07856 Liu, B., Gurian, J., Inayoshi, K., et al. 2024, MNRAS, 534, 290 Lupi, A., Colpi, M., Devecchi, B., Galanti, G., & Volonteri, M. 2014, MNRAS, 442, 3616 Maiolino, R., Scholtz, J., Curtis-Lake, E., et al. 2024a, A&A, 691, A145 Maiolino, R., Scholtz, J., Witstok, J., et al. 2024b, Natur, 627, 59 Mandel, I., & Broekgaarden, F. S. 2022, LRR, 25, 1 Mapelli, M. 2016, MNRAS, 459, 3432 Marsat, S., Baker, J. G., & Dal Canton, T. 2021, PhRvD, 103, 083011 Matthee, J., et al. 2024, ApJ, 963, 129 Mayer, L., van Donkelaar, F., Messa, M., Capelo, P. R., & Adamo, A. 2025, ApJL, 981, L28 Mehta, D., Regan, J. A., & Prole, L. 2024, OJAp, 7, 107 Merritt, D. 2013, Dynamics and Evolution of Galactic Nuclei (Princeton, NJ: Princeton Univ. Press) Mezcua, M., Siudek, M., Suh, H., et al. 2023, ApJL, 943, L5 Miller, B. P., Gallo, E., Greene, J. E., et al. 2015, ApJ, 799, 98 Mouri, H., & Taniguchi, Y. 2002, ApJL, 566, L17 Mowla, L., Iyer, K. G., Desprez, G., et al. 2022, ApJL, 937, L35 Nadler, E. O., Benson, A., Driskell, T., Du, X., & Gluscevic, V. 2023, MNRAS, 521, 3201 Natarajan, P. 2011, BASI, 39, 145 Natarajan, P. 2021, MNRAS, 501, 1413 Neumayer, N., Seth, A., & Boeker, T. 2020, A&ARv, 28, 4 Neumayer, N., Seth, A., & Böker, T. 2020, A&ARv, 28, 4 Nguyen, D. D., Seth, A. C., Neumayer, N., et al. 2019, ApJ, 872, 104 O’Leary, R. M., Kocsis, B., & Loeb, A. 2009, MNRAS, 395, 2127 Omukai, K., Schneider, R., & Haiman, Z. 2008, ApJ, 686, 801 Partmann, C., Naab, T., Lahén, N., et al. 2025, MNRAS, 537, 956 Pfister, H., Dai, J., Volonteri, M., et al. 2020, MNRAS, 500, 3944 Polkas, M., Bonoli, S., Bortolas, E., et al. 2024, A&A, 689, A204 102 103 104 105 106 107 108 109 1010 1011 1012 m1/M� 100 101 C ou n t 10−8 10−7 10−6 10−5 10−4 10−3 10−2 10−1 100 q 100 101 C ou n t Figure 17. The same 11 realizations as in Figure 16 for the primary masses (left) and mass ratios (right). 17 The Astrophysical Journal, 991:58 (18pp), 2025 September 20 Kritos et al. https://orcid.org/0000-0002-0212-3472 https://orcid.org/0000-0002-2850-0192 https://orcid.org/0000-0002-1566-8148 https://orcid.org/0000-0003-0751-5130 https://orcid.org/0000-0002-9104-1734 https://orcid.org/0000-0002-3216-1322 https://orcid.org/0000-0003-0225-6387 https://orcid.org/0000-0002-8140-0422 https://doi.org/10.1103/PhysRevX.13.011048 https://ui.adsabs.harvard.edu/abs/2023PhRvX..13a1048A/abstract https://doi.org/10.1038/s41586-024-07703-7 https://ui.adsabs.harvard.edu/abs/2024Natur.632..513A/abstract https://doi.org/10.1126/science.1251053 https://ui.adsabs.harvard.edu/abs/2014Sci...345.1330A/abstract http://arXiv.org/abs/1702.00786 http://arXiv.org/abs/1702.00786 https://doi.org/10.1093/mnras/stz1149 https://ui.adsabs.harvard.edu/abs/2019MNRAS.486.5008A/abstract https://doi.org/10.1111/j.1365-2966.2004.07883.x https://ui.adsabs.harvard.edu/abs/2004MNRAS.352..376A/abstract https://doi.org/10.1103/PhysRevD.95.103012 https://ui.adsabs.harvard.edu/abs/2017PhRvD..95j3012B/abstract https://doi.org/10.1103/PhysRevD.100.064060 https://ui.adsabs.harvard.edu/abs/2019PhRvD.100f4060B/abstract https://doi.org/10.1103/PhysRevD.104.084002 https://ui.adsabs.harvard.edu/abs/2021PhRvD.104h4002B/abstract https://doi.org/10.1038/226064a0 https://ui.adsabs.harvard.edu/abs/1970Natur.226...64B/abstract https://doi.org/10.1111/j.1365-2966.2007.12209.x https://ui.adsabs.harvard.edu/abs/2007MNRAS.380.1589B/abstract https://doi.org/10.1093/mnras/stad1544 https://ui.adsabs.harvard.edu/abs/2023MNRAS.523.5610B/abstract https://doi.org/10.1093/mnras/stae2595 https://ui.adsabs.harvard.edu/abs/2025MNRAS.536.1838B/abstract https://doi.org/10.1086/590379 https://ui.adsabs.harvard.edu/abs/2008ApJ...684..822B/abstract https://doi.org/10.1093/mnras/stae780 https://ui.adsabs.harvard.edu/abs/2024MNRAS.529.3768B/abstract https://doi.org/10.1093/mnras/staf269 https://ui.adsabs.harvard.edu/abs/2025MNRAS.538..518B/abstract https://doi.org/10.1093/mnras/179.3.433 https://ui.adsabs.harvard.edu/abs/1977MNRAS.179..433B/abstract https://doi.org/10.1093/mnras/112.2.195 https://ui.adsabs.harvard.edu/abs/1952MNRAS.112..195B/abstract https://doi.org/10.1111/j.1365-2966.2009.14701.x https://ui.adsabs.harvard.edu/abs/2009MNRAS.396..423B/abstract https://doi.org/10.1093/mnras/stt628 https://ui.adsabs.harvard.edu/abs/2013MNRAS.432.2779B/abstract https://doi.org/10.1093/mnras/stae2663 https://ui.adsabs.harvard.edu/abs/2025MNRAS.536..851C/abstract https://doi.org/10.1086/516712 https://ui.adsabs.harvard.edu/abs/2007ApJ...659L...5C/abstract https://ui.adsabs.harvard.edu/abs/2007ApJ...659L...5C/abstract https://doi.org/10.1007/978-1-4020-3407-7_5 https://ui.adsabs.harvard.edu/abs/2005ASSL..327...41C/abstract http://arXiv.org/abs/2501.18730 https://doi.org/10.1051/0004-6361/202452393 https://ui.adsabs.harvard.edu/abs/2025A&A...695A..97D/abstract https://doi.org/10.1111/j.1365-2966.2010.17363.x https://ui.adsabs.harvard.edu/abs/2010MNRAS.409.1057D/abstract https://ui.adsabs.harvard.edu/abs/2010MNRAS.409.1057D/abstract https://doi.org/10.1111/j.1365-2966.2012.20406.x https://ui.adsabs.harvard.edu/abs/2012MNRAS.421.1465D/abstract https://doi.org/10.1086/524921 https://ui.adsabs.harvard.edu/abs/2008ApJ...676...33D/abstract https://doi.org/10.1051/0004-6361/202346435 https://ui.adsabs.harvard.edu/abs/2023A&A...676A...2D/abstract https://ui.adsabs.harvard.edu/abs/2023A&A...676A...2D/abstract https://doi.org/10.1051/0004-6361/202346295 https://ui.adsabs.harvard.edu/abs/2023A&A...673A.120D/abstract https://ui.adsabs.harvard.edu/abs/2023A&A...673A.120D/abstract https://doi.org/10.1051/0004-6361/202453070 https://ui.adsabs.harvard.edu/abs/2025A&A...695A.231D/abstract https://doi.org/10.1088/0004-637X/762/2/68 https://ui.adsabs.harvard.edu/abs/2013ApJ...762...68D/abstract https://ui.adsabs.harvard.edu/abs/2013ApJ...762...68D/abstract https://doi.org/10.1051/0004-6361/202039429 https://ui.adsabs.harvard.edu/abs/2021A&A...651A.109D/abstract https://doi.org/10.1093/mnras/stu373 https://ui.adsabs.harvard.edu/abs/2014MNRAS.440.1590D/abstract https://doi.org/10.1093/mnras/stu425 https://ui.adsabs.harvard.edu/abs/2014MNRAS.440.2333D/abstract https://ui.adsabs.harvard.edu/abs/2014MNRAS.440.2333D/abstract http://arXiv.org/abs/2410.24224 https://doi.org/10.1088/0004-637X/716/2/1397 https://ui.adsabs.harvard.edu/abs/2010ApJ...716.1397F/abstract https://ui.adsabs.harvard.edu/abs/2010ApJ...716.1397F/abstract https://doi.org/10.1126/science.adi4211 https://ui.adsabs.harvard.edu/abs/2024Sci...384.1488F/abstract https://ui.adsabs.harvard.edu/abs/2024Sci...384.1488F/abstract https://doi.org/10.1103/PhysRevD.110.063012 https://ui.adsabs.harvard.edu/abs/2024PhRvD.110f3012F/abstract https://doi.org/10.1103/PhysRevD.102.064002 https://ui.adsabs.harvard.edu/abs/2020PhRvD.102f4002G/abstract https://doi.org/10.1093/mnras/stw093 https://ui.adsabs.harvard.edu/abs/2016MNRAS.457.2122G/abstract https://doi.org/10.1103/PhysRevD.93.124066 https://ui.adsabs.harvard.edu/abs/2016PhRvD..93l4066G/abstract https://doi.org/10.1146/annurev-astro-111720-030029 https://ui.adsabs.harvard.edu/abs/2021ARA&A..59...21G/abstract https://doi.org/10.1088/0004-637X/785/1/71 https://ui.adsabs.harvard.edu/abs/2014ApJ...785...71G/abstract https://doi.org/10.1103/PhysRevLett.98.231101 https://ui.adsabs.harvard.edu/abs/2007PhRvL..98w1101G/abstract https://doi.org/10.3847/1538-4357/ad1e5f https://ui.adsabs.harvard.edu/abs/2024ApJ...964...39G/abstract https://doi.org/10.1146/annurev-astro-032620-021835 https://ui.adsabs.harvard.edu/abs/2020ARA&A..58..257G/abstract https://doi.org/10.1086/177035 https://ui.adsabs.harvard.edu/abs/1996ApJ...461...20H/abstract https://doi.org/10.1086/591218 https://ui.adsabs.harvard.edu/abs/2008ApJ...686..829H/abstract https://doi.org/10.1088/1475-7516/2024/04/021 https://ui.adsabs.harvard.edu/abs/2024JCAP...04..021H/abstract https://doi.org/10.1086/375495 https://ui.adsabs.harvard.edu/abs/2003ApJ...585L.101H/abstract https://doi.org/10.1093/mnras/stab3239 https://ui.adsabs.harvard.edu/abs/2022MNRAS.509.3488I/abstract https://ui.adsabs.harvard.edu/abs/2022MNRAS.509.3488I/abstract https://doi.org/10.1038/s41586-024-08210-5 https://ui.adsabs.harvard.edu/abs/2024Natur.636..594J/abstract https://doi.org/10.1088/0004-637X/773/2/100 https://ui.adsabs.harvard.edu/abs/2013ApJ...773..100K/abstract https://ui.adsabs.harvard.edu/abs/2013ApJ...773..100K/abstract https://doi.org/10.1093/mnras/stv1211 https://ui.adsabs.harvard.edu/abs/2015MNRAS.451.2900K/abstract https://ui.adsabs.harvard.edu/abs/2015MNRAS.451.2900K/abstract https://doi.org/10.3847/0004-637X/823/1/52 https://ui.adsabs.harvard.edu/abs/2016ApJ...823...52K/abstract https://doi.org/10.1093/mnras/stx052 https://ui.adsabs.harvard.edu/abs/2017MNRAS.466.4826K/abstract https://doi.org/10.1111/j.1365-2966.2008.12943.x https://ui.adsabs.harvard.edu/abs/2008MNRAS.385.1621K/abstract https://doi.org/10.1088/0067-0049/192/2/18 https://ui.adsabs.harvard.edu/abs/2011ApJS..192...18K/abstract https://doi.org/10.3847/2041-8213/ad391f https://ui.adsabs.harvard.edu/abs/2024ApJ...965L..21K/abstract https://doi.org/10.5281/zenodo.15776363 https://doi.org/10.5281/zenodo.15776363 https://doi.org/10.1103/PhysRevD.108.083012 https://ui.adsabs.harvard.edu/abs/2023PhRvD.108h3012K/abstract https://doi.org/10.1093/mnras/stae1145 https://ui.adsabs.harvard.edu/abs/2024MNRAS.531..133K/abstract https://doi.org/10.1103/PhysRevD.110.123017 https://ui.adsabs.harvard.edu/abs/2024PhRvD.110l3017K/abstract https://doi.org/10.1126/science.1067524 https://ui.adsabs.harvard.edu/abs/2002Sci...295...82K/abstract https://doi.org/10.1093/mnras/staa2276 https://ui.adsabs.harvard.edu/abs/2020MNRAS.498.5652K/abstract https://doi.org/10.1093/mnras/stae2694 https://ui.adsabs.harvard.edu/abs/2025MNRAS.536.3366L/abstract https://ui.adsabs.harvard.edu/abs/2025MNRAS.536.3366L/abstract http://arXiv.org/abs/2410.07856 https://doi.org/10.1093/mnras/stae2066 https://ui.adsabs.harvard.edu/abs/2024MNRAS.534..290L/abstract https://doi.org/10.1093/mnras/stu1120 https://ui.adsabs.harvard.edu/abs/2014MNRAS.442.3616L/abstract https://doi.org/10.1051/0004-6361/202347640 https://ui.adsabs.harvard.edu/abs/2024A&A...691A.145M/abstract https://doi.org/10.1038/s41586-024-07052-5 https://ui.adsabs.harvard.edu/abs/2024Natur.627...59M/abstract https://doi.org/10.1007/s41114-021-00034-3 https://ui.adsabs.harvard.edu/abs/2022LRR....25....1M/abstract https://doi.org/10.1093/mnras/stw869 https://ui.adsabs.harvard.edu/abs/2016MNRAS.459.3432M/abstract https://doi.org/10.1103/PhysRevD.103.083011 https://ui.adsabs.harvard.edu/abs/2021PhRvD.103h3011M/abstract https://doi.org/10.3847/1538-4357/ad2345 https://ui.adsabs.harvard.edu/abs/2024ApJ...963..129M/abstract https://doi.org/10.3847/2041-8213/adadfe https://ui.adsabs.harvard.edu/abs/2025ApJ...981L..28M/abstract https://doi.org/10.33232/001c.126629 https://ui.adsabs.harvard.edu/abs/2024OJAp....7E.107M/abstract https://doi.org/10.3847/2041-8213/acae25 https://ui.adsabs.harvard.edu/abs/2023ApJ...943L...5M/abstract https://doi.org/10.1088/0004-637X/799/1/98 https://ui.adsabs.harvard.edu/abs/2015ApJ...799...98M/abstract https://doi.org/10.1086/339472 https://ui.adsabs.harvard.edu/abs/2002ApJ...566L..17M/abstract https://doi.org/10.3847/2041-8213/ac90ca https://ui.adsabs.harvard.edu/abs/2022ApJ...937L..35M/abstract https://doi.org/10.1093/mnras/stad666 https://ui.adsabs.harvard.edu/abs/2023MNRAS.521.3201N/abstract https://ui.adsabs.harvard.edu/abs/2011BASI...39..145N/abstract https://doi.org/10.1093/mnras/staa3724 https://ui.adsabs.harvard.edu/abs/2021MNRAS.501.1413N/abstract https://doi.org/10.1007/s00159-020-00125-0 https://ui.adsabs.harvard.edu/abs/2020A&ARv..28....4N/abstract https://doi.org/10.1007/s00159-020-00125-0 https://ui.adsabs.harvard.edu/abs/2020A&ARv..28....4N/abstract https://doi.org/10.3847/1538-4357/aafe7a https://ui.adsabs.harvard.edu/abs/2019ApJ...872..104N/abstract https://doi.org/10.1111/j.1365-2966.2009.14653.x https://ui.adsabs.harvard.edu/abs/2009MNRAS.395.2127O/abstract https://doi.org/10.1086/591636 https://ui.adsabs.harvard.edu/abs/2008ApJ...686..801O/abstract https://doi.org/10.1093/mnras/staf002 https://ui.adsabs.harvard.edu/abs/2025MNRAS.537..956P/abstract https://doi.org/10.1093/mnras/staa3471 https://ui.adsabs.harvard.edu/abs/2021MNRAS.500.3944P/abstract https://doi.org/10.1051/0004-6361/202449470 https://ui.adsabs.harvard.edu/abs/2024A&A...689A.204P/abstract Portegies Zwart, S. F., & McMillan, S. L. W. 2002, ApJ, 576, 899 Portegies Zwart, S. F., McMillan, S. L. W., & Gieles, M. 2010, ARA&A, 48, 431 Prieto, E. G., Kremer, K., Fragione, G., et al. 2022, ApJ, 940, 131 Prole, L. R., Regan, J. A., Whalen, D. J., Glover, S. C. O., & Klessen, R. S. 2024, A&A, 692, A213 Quinlan, G. D. 1996, NewA, 1, 35 Quinlan, G. D., & Shapiro, S. L. 1987, ApJ, 321, 199 Rantala, A., Naab, T., & Lahén, N. 2024, MNRAS, 531, 3770 Rasio, F. A., Freitag, M., & Gürkan, M. A. 2004, in Carnegie Observatories Centennial Symp., Coevolution of Black Holes and Galaxies, ed. L. C. Ho (Cambridge: Cambridge Univ. Press), 138 Reinoso, B., Schleicher, D. R. G., Fellhauer, M., Klessen, R. S., & Boekholt, T. C. N. 2018, A&A, 614, A14 Reynolds, C. S. 2021, ARA&A, 59, 117 Rizzuto, F. P., Naab, T., Rantala, A., et al. 2023, MNRAS, 521, 2930 Rodriguez, C. L., Chatterjee, S., & Rasio, F. A. 2016, PhRvD, 93, 084029 Sánchez-Janssen, R., Côté, P., Ferrarese, L., et al. 2019, ApJ, 878, 18 Sijacki, D., Vogelsberger, M., Genel, S., et al. 2015, MNRAS, 452, 575 Silk, J., Begelman, M. C., Norman, C., Nusser, A., & Wyse, R. F. G. 2024, ApJL, 961, L39 Soltan, A. 1982, MNRAS, 200, 115 Suh, H., Scharwächter, J., Farina, E. P., et al. 2025, NatAs, 9, 271 Sutherland, R. S., & Dopita, M. A. 1993, ApJS, 88, 253 Teyssier, R. 2002, A&A, 385, 337 Trebitsch, M., Blaizot, J., Rosdahl, J., Devriendt, J., & Slyz, A. 2017, MNRAS, 470, 224 Trebitsch, M., Volonteri, M., & Dubois, Y. 2020, MNRAS, 494, 3453 Tremmel, M., Governato, F., Volonteri, M., & Quinn, T. R. 2015, MNRAS, 451, 1868 Vanzella, E., Claeyssens, A., Welch, B., et al. 2023, ApJ, 945, 53 Vergara, M. C., Escala, A., Schleicher, D. R. G., & Reinoso, B. 2023, MNRAS, 522, 4224 Volonteri, M., Pfister, H., Beckmann, R. S., et al. 2020, MNRAS, 498, 2219 Welch, B., Coe, D., Zitrin, A., et al. 2023, ApJ, 943, 2 Zevin, M., Romero-Shaw, I. M., Kremer, K., Thrane, E., & Lasky, P. D. 2021, ApJL, 921, L43 Ziparo, F., Gallerani, S., & Ferrara, A. 2025, JCAP, 2025, 040 18 The Astrophysical Journal, 991:58 (18pp), 2025 September 20 Kritos et al. https://doi.org/10.1086/341798 https://ui.adsabs.harvard.edu/abs/2002ApJ...576..899P/abstract https://doi.org/10.1146/annurev-astro-081309-130834 https://ui.adsabs.harvard.edu/abs/2010ARA&A..48..431P/abstract https://ui.adsabs.harvard.edu/abs/2010ARA&A..48..431P/abstract https://doi.org/10.3847/1538-4357/ac9b0f https://ui.adsabs.harvard.edu/abs/2022ApJ...940..131G/abstract https://doi.org/10.1051/0004-6361/202452486 https://ui.adsabs.harvard.edu/abs/2024A&A...692A.213P/abstract https://doi.org/10.1016/S1384-1076(96)00003-6 https://ui.adsabs.harvard.edu/abs/1996NewA....1...35Q/abstract https://doi.org/10.1086/165624 https://ui.adsabs.harvard.edu/abs/1987ApJ...321..199Q/abstract https://doi.org/10.1093/mnras/stae1413 https://ui.adsabs.harvard.edu/abs/2024MNRAS.531.3770R/abstract https://ui.adsabs.harvard.edu/abs/2004cbhg.symp..138R/abstract https://doi.org/10.1051/0004-6361/201732224 https://ui.adsabs.harvard.edu/abs/2018A&A...614A..14R/abstract https://doi.org/10.1146/annurev-astro-112420-035022 https://ui.adsabs.harvard.edu/abs/2021ARA&A..59..117R/abstract https://doi.org/10.1093/mnras/stad734 https://ui.adsabs.harvard.edu/abs/2023MNRAS.521.2930R/abstract https://doi.org/10.1103/PhysRevD.93.084029 https://ui.adsabs.harvard.edu/abs/2016PhRvD..93h4029R/abstract https://doi.org/10.3847/1538-4357/aaf4fd https://ui.adsabs.harvard.edu/abs/2019ApJ...878...18S/abstract https://doi.org/10.1093/mnras/stv1340 https://ui.adsabs.harvard.edu/abs/2015MNRAS.452..575S/abstract https://doi.org/10.3847/2041-8213/ad1bf0 https://ui.adsabs.harvard.edu/abs/2024ApJ...961L..39S/abstract https://doi.org/10.1093/mnras/200.1.115 https://ui.adsabs.harvard.edu/abs/1982MNRAS.200..115S/abstract https://doi.org/10.1038/s41550-024-02402-9 https://ui.adsabs.harvard.edu/abs/2025NatAs...9..271S/abstract https://doi.org/10.1086/191823 https://ui.adsabs.harvard.edu/abs/1993ApJS...88..253S/abstract https://doi.org/10.1051/0004-6361:20011817 https://ui.adsabs.harvard.edu/abs/2002A&A...385..337T/abstract https://doi.org/10.1093/mnras/stx1060 https://ui.adsabs.harvard.edu/abs/2017MNRAS.470..224T/abstract https://doi.org/10.1093/mnras/staa1012 https://ui.adsabs.harvard.edu/abs/2020MNRAS.494.3453T/abstract https://doi.org/10.1093/mnras/stv1060 https://ui.adsabs.harvard.edu/abs/2015MNRAS.451.1868T/abstract https://ui.adsabs.harvard.edu/abs/2015MNRAS.451.1868T/abstract https://doi.org/10.3847/1538-4357/acb59a https://ui.adsabs.harvard.edu/abs/2023ApJ...945...53V/abstract https://doi.org/10.1093/mnras/stad1253 https://ui.adsabs.harvard.edu/abs/2023MNRAS.522.4224V/abstract https://doi.org/10.1093/mnras/staa2384 https://ui.adsabs.harvard.edu/abs/2020MNRAS.498.2219V/abstract https://doi.org/10.3847/1538-4357/aca8a8 https://ui.adsabs.harvard.edu/abs/2023ApJ...943....2W/abstract https://doi.org/10.3847/2041-8213/ac32dc https://ui.adsabs.harvard.edu/abs/2021ApJ...921L..43Z/abstract https://doi.org/10.1088/1475-7516/2025/04/040 https://ui.adsabs.harvard.edu/abs/2025JCAP...04..040Z/abstract 1. Introduction 2. Methods 2.1. Merger Trees 2.2. NSC Initial Conditions 2.3. Massive BH Seeds 2.3.1. BH Growth Channels 2.3.2. Time Evolution of MBH and χBH 2.3.3. Mapping BH Masses in the NSC Parameter Space 2.4. Major Merger Episodes 3. Results 3.1. Local Massive BH Properties 3.1.1. The MBH–M⋆ Relation 3.1.2. Massive BH Spin 3.1.3. Occupation Fraction 3.1.4. Number Density 3.2. Local NSC Properties 3.2.1. The MBH–MNSC Relation 3.2.2. The rh–MNSC and MNSC–M⋆ Relations 3.3. Transient Events 3.3.1. Massive BH Binary Merger Rate 3.3.2. Stellar-mass BH–BH Capture Rate 3.3.3. Extreme Mass Ratio Inspirals 3.3.4. TDE Rate 4. Caveats 4.1. Limitations 4.2. Dependence on Hyperparameters 4.2.1. Lower Mass Threshold for NSC Formation 4.2.2. Lower Mass Ratio Threshold for Major Galaxy Mergers 4.2.3. Smaller Gas Expulsion Timescale 4.2.4. fg,max=0.1 4.2.5. f⋆,max=0.1 5. Conclusions Appendix References