Impact of gas hardening on the population properties of hierarchical black hole mergers in active galactic nucleus disks (cid:63)

Hierarchical black hole (BH) mergers in active galactic nuclei (AGNs) are unique among formation channels of binary black holes (BBHs) because they are likely associated with electromagnetic counterparts and can e ﬃ ciently lead to the mass growth of BHs. Here, we explore the impact of gas accretion and migration traps on the evolution of BBHs in AGNs. We have developed a new fast semi-analytic model, that allows us to explore the parameter space while capturing the main physical processes involved. We ﬁnd that an e ﬀ ective exchange of energy and angular momentum between the BBH and the surrounding gas (i.e., gas hardening) during inspiral greatly enhances the e ﬃ ciency of hierarchical mergers, leading to the formation of intermediate-mass BHs (up to 10 4 M (cid:12) ) and triggering spin alignment. Moreover, our models with e ﬃ cient gas hardening show both an anticorrelation between the BBH mass ratio and the e ﬀ ective spin and a correlation between the primary BH mass and the e ﬀ ective spin. In contrast, if gas hardening is ine ﬃ cient, the hierarchical merger chain is already truncated after the ﬁrst two or three generations. We compare the BBH population in AGNs with other dynamical channels as well as isolated binary evolution.

Stellar dynamics provides some of the most straightforward channels to explain the birth of such oversized BHs, via star-star collisions (Di Carlo et al. 2019, 2020;Kremer et al. 2020;Renzo et al. 2020;Torniamenti et al. 2022;Costa et al. 2022;Ballone et al. 2023) or repeated mergers of stellar- The data underlying this article will be shared upon reasonable request to the corresponding authors.
In AGNs, a central supermassive black hole (SMBH) is surrounded by a dense gaseous accretion disk.AGNs can appear as extremely luminous objects called quasars, but also as Seyfert galaxies, radio galaxies, or blazars, depending on their luminosity and our viewing angle (see Netzer 2015 for a review of the unification scheme and its controversy).Stars and stellar-sized BHs orbiting the SMBH are subject to gas torques, which can bend their orbits, aligning them with the disk.This is expected to lead to a large overdensity of BHs with similar orbits in small areas of the disk called migration traps (McKernan et al. 2012;Bellovary et al. 2016), where BHs can pair up efficiently via gas capture (DeLaurentiis et al. 2023;Li et al. 2023;Rowan et al. 2023Rowan et al. , 2024;;Whitehead et al. 2023).Therefore, AGNs are potential factories of BBH mergers with masses in the pair-instability mass gap and above (Yang et al. 2019a;Tagawa et al. 2020a,b).
This channel has attracted great interest because GW signals from mergers could be accompanied by some electromagnetic emission (Bartos et al. 2017;Tagawa et al. 2023), possibly anticipated by a neutrino detection (Zhu 2024;Zhou & Wang 2023).To explore this possibility, electromagnetic follow-up observations have been carried out for several GW events, but so far all associations with BBH mergers remain controversial (Greiner et al. 2016;Coughlin et al. 2020;Calderón Bustillo et al. 2021).
Several studies have investigated the possibility that the high-mass BBH merger event GW190521 is associated with an AGN disk (Tagawa et al. 2021;Samsing et al. 2022) because of its large mass, high spin (Abbott et al. 2020), and claimed electromagnetic counterpart (Graham et al. 2020;Calderón Bustillo et al. 2021;Morton et al. 2023).Yang et al. (2019a) explored the AGN scenario for another high-mass BBH, GW170729, for which there is evidence of nonzero effective and precessing spins.Some authors, such as Gayathri et al. (2021) and Ford & McKernan (2022), have predicted that a sizeable fraction of the observed BBH mergers (∼20% up to 80%) may originate in AGNs.However, according to an analysis based on the sky localization of GW signals, the fraction of detected BBH mergers that originated in bright AGNs (L ≥ 10 46 erg s −1 ) cannot be higher than 17% (Veronesi et al. 2023).
Binary BH pair-up in AGNs can happen either in the migration traps or at other locations in the disk (Wang et al. 2021) collectively referred to as 'the bulk'.McKernan et al. (2020) find that, although more than 50% of mergers happen in the bulk, hierarchical mergers are only efficient in migration traps.Moreover, Tagawa et al. (2020b) find that BBHs assembled in the bulk via gas capture migrate toward the migration trap while hardening.
The complex physics of BBHs in AGN disks has been extensively explored in the literature, with studies accounting for the effects of binary-single interactions (Stone et al. 2017;Leigh et al. 2017) and gas torques (using models borrowed from protoplanetary physics; e.g., McKernan et al. 2012; Bartos et al. 2017;Yang et al. 2019a,b;Secunda et al. 2020;or hydrodynamical simulations;e.g., Li et al. 2023;Kaaz et al. 2023;Rowan et al. 2023Rowan et al. , 2024;;Whitehead et al. 2023).For example, Ishibashi & Gröbner (2020) explore the effects of energy and angular momentum exchange between the BBH and the surrounding gas.In particular, they assume that the BBH evolves surrounded by a circumbinary disk, which induces orbital decay and increases its orbital eccentricity (i.e., gas hardening).
We studied hierarchical BBH mergers in the migration traps of AGN disks and compared our results to mergers in YSCs, GCs, and NSCs (Mapelli et al. 2021(Mapelli et al. , 2022)).Specifically, we tested the impact of gas hardening (GH) on the BBH merger population in AGNs (Ishibashi & Gröbner 2020).We find that the hierarchical merger process is pronouncedly more efficient when accounting for gas hardening.
We have developed a new semi-analytic code for the simulation of hierarchical mergers in AGNs.It is effective at exploring the parameter space, and it models the relevant physical processes much faster than an N-body or hydrodynamical code.The new code is publicly available within the Fastcluster software environment 1 (Mapelli et al. 2021(Mapelli et al. , 2022)).

AGN disk model
We assumed the AGN disk to be as described by Sirko & Goodman (2003, hereafter, SG).These authors introduce a hydro-dynamical model for a geometrically thin and optically thick disk with steady-state accretion of 0.1 ṀEdd onto the central SMBH.Gas turbulence is assumed to be the cause of disk viscosity, characterized by the viscosity coefficient α ∈ [0, 1].This model neglects any effects due to magnetic fields and general relativity.The physical parameters of the disk are functions of the mass of the SMBH M SMBH , the distance R from the SMBH, and the viscosity parameter α.We slightly simplified the radial dependence of the SG model for numerical ease and re-scaled their functions to allow for different M SMBH and α parameters, as shown in Fig. 1.The resulting expressions for the gas surface density, Σ gas , the disk aspect ratio, h = H/R (defined as the ratio between the height of the disk, H (R), and its radius, R), and the sound speed, c s , are the following: In the above equations, we define the gravitational radius as R g = GM SMBH /c 2 , where G is the gravity constant and c the speed of light.
The viscosity coefficient α is a free parameter.We assumed it to be constant over the whole extension of the disk and to have a constant value independent of the other physical properties of the disk.The assumed value is α = 0.1 which is consistent with observations (King et al. 2007).
Migration traps are locations in the disk where migration stalls and BHs pile up.Bellovary et al. (2016) find that migration traps are located where the slope of the gas surface density profile changes sign from positive to negative (i.e. at local maxima).In a SG disk there are two local maxima in the gas surface density and therefore two migration traps: an inner trap at ≈100 R g and an outer one at ≈1300 R g .
We see in Fig. 1 that the inner migration trap coincides with a local maximum in the density profile of the original SG model, while the outer migration trap corresponds to a global maximum.For simplicity, we ignored the local overdensity in the disk at 10 2 R g and we only assumed the existence of the outer migration trap.We thus defined the trap radius as R trap = 10 3.1 R g , although its location and existence may be affected by other effects (e.g.Grishin et al. 2024;Pan & Yang 2021 more details in Sect. 4.4).
The disk is assumed to have radial extension between the innermost stable circular orbit (ISCO) radius for a nonrotating BH, which we call R min = 6 R g in this context, and an outer radius R max = 0.1 pc(M SMBH /10 6 M ) 1/2 , beyond which the disk's self-gravity becomes important (Goodman 2003;Yang et al. 2019b).For R > R max , the disk is expected to fragment and experience star formation.Energetic feedback from the newly formed stars may keep the disk vertically supported (Sirko & Goodman 2003), but the outcome of viscous interactions with the BHs in this region of the disk is highly uncertain, so we conservatively neglected it.Moreover, we neglected the contribution of stars formed in AGN disks, which are expected to be massive and whose compact remnants might also become GW sources (Cantiello et al. 2021;Wang et al. 2023).
We assumed a non-spinning SMBH.We randomly generated the SMBH mass from the observational distribution derived by Greene & Ho (2007) in the Local Universe, which is a Gaussian distribution with mean log M SMBH /M = 6.576 ± 0.591.
A mass-accretion episode onto a SMBH lasts for a finite amount of time, which we refer to as the AGN disk lifetime.The lifetime of AGN disks is subject to large uncertainty: different estimates span several orders of magnitude in the range of 10 −2 −10 3 Myr (Khrykin et al. 2021).Also, it is not clear whether accretion onto SMBHs happens continuously over a given time span, or episodically through many cycles of efficient accretion.Here, we used the estimate by Khrykin et al. (2021), based on observations of quasars' proximity effect.They find that the quasar lifetime τ is distributed according to a Gaussian distribution with mean log τ/Myr = 0.22 ± 0.80.We randomly extracted the lifetime for AGN disks in our model from this distribution.
There are other models for stable AGN accretion disks, with notably different features compared to the one we adopt here.For example, the model by Thompson et al. (2005) differs from the SG model in both density and aspect ratio (McKernan et al. 2022).We will explore the impact of different disk models in a follow-up study.

Nuclear star cluster (NSC)
Supermassive black holes and NSCs commonly inhabit galactic spheroids with stellar masses ranging from 10 8 M to 10 11 M (Graham & Spitler 2009).For M SMBH ≤ 5 × 10 7 M , the mass of the SMBH and that of the NSC scale as log For M SMBH > 5 × 10 7 M , the SMBH generally does not coexist with a NSC and we set M NSC = 0.The effective radius of a NSC mildly correlates to its mass (Neumayer et al. 2020).The best-fit relation between the mass of the NSC and its effective radius is log (5) To account for the spread in the data, we sampled log r h /pc from a Gaussian distribution with the mean centered on the corresponding value determined by Eq. ( 5) and width σ = 0.1 (0.2) for M NSC < 10 6 M (≥10 6 M ) so that most of the data fell under the ±2σ dispersion.
We approximated the spatial distribution of stars in the NSC with a Plummer (1911) model, with mass where R is the distance from the SMBH, and a PL = r h /(1.3 pc) is the scale parameter for the Plummer model.The mass fraction of stellar-origin BHs in the AGN disk f BH is sampled from a Gaussian with mean 0.04 and standard deviation 0.01 to account for mass segregation in the NSC (Bartos et al. 2017).The number of BHs in the AGN disk and their cumulative mass are thus given by the following equations: A51, page 3 of 20 where M NSC (R max ) is defined in Eq. ( 6), the mean stellar-origin BH mass m BH is computed from data obtained from the population synthesis simulation code sevn2 (Iorio et al. 2023) at solar metallicity (Appendix A), the distance R max is the outer radius of the disk (Sect.2.1), the factor 1/2 accounts for prograde orbiters3 only, and the mean stellar mass m * 1 M is computed using a Kroupa (2001) initial mass function.
The parameter f trap is the ratio between the number of BHs that are able to reach the migration trap on a timescale shorter than the disk lifetime τ and the total number of BHs that interact with the disk: this study focuses on BBH pair-ups in the migration trap, and therefore we are not interested in any BHs situated outside of that location.An operational definition of f trap is given in Sect.2.3.3.
We assumed the velocity dispersion of stars to scale with the SMBH mass as (Merritt & Ferrarese 2001) However, it has been shown (e.g.Scott & Graham 2013;Sahu et al. 2019;Graham 2022) that this result depends on the morphology of galaxies included in the sample.For example, Sahu et al. (2019) find an exponent of ∼1/6 and show that it is caused by the combined contribution of Sérsic (i.e., following the Sérsic 1963 brightness profile) and core-Sérsic (i.e., centrally depleted) galaxies following two different M SMBH − σ relations with exponents (5.75 ± 0.34) −1 and (8.64 ± 1.10) −1 , respectively.

First-generation (1g) BHs
We randomly drew first-generation (1g, i.e., stellar-origin) BH masses m BH from a catalog obtained with the population synthesis code sevn (Spera & Mapelli 2017;Spera et al. 2019;Mapelli et al. 2020;Iorio et al. 2023).Sevn relies on up-to-date stellar tracks (Bressan et al. 2012;Costa et al. 2019;Nguyen et al. 2022) and models the formation of compact objects by taking into account electron-capture (Giacobbo & Mapelli 2019), core-collapse (Iorio et al. 2023) and pair-instability supernovae (Mapelli et al. 2020).In particular, here we used the rapid core-collapse supernova model from Fryer et al. (2012), which enforces the existence of a mass gap between the maximum neutron star mass and the minimum BH mass (Özel et al. 2010).We used the fiducial model from Iorio et al. (2023) and considered single stellar evolution only, as described in Appendix A. We assumed metallicity Z = 0.02 (i.e., approximately solar), which matches the typical metallicity at the center of massive galaxies in the Local Universe (Gallazzi et al. 2008).We randomly drew an initial radial position for each BH.The radial extension of the accretion disk is small compared to the typical dimension of a NSC, so we neglected mass segregation on this scale and considered the numerical density of objects at radii R < R max to be uniform in radius: R ∼ U (R min , R max ).
We drew the dimensionless spin magnitude, χ, from a Maxwellian distribution with one-dimensional root-meansquare σ χ = 0.1, truncated at χ = 1.We chose σ χ = 0.1 because it is reminiscent of the spins inferred from the third GW transient catalog (GWTC-3, Abbott et al. 2023b).This assumption does not take into account that the BBH population in GWTC-3 likely comes from multiple formation channels, including the AGN disk scenario.Current data are not sufficiently informative to differentiate between formation channels.We set the primary spin tilt as in Appendix B.

Gas capture
After setting up the properties of 1g BHs, we followed their evolution in the disk.When NSC objects orbit around the central SMBH, their orbits can cross the disk and gather some of the disk gas, subjecting them to strong gas drag.This is expected to dampen both the inclination i and the eccentricity e of their orbit (Cresswell et al. 2007).Therefore, after a sufficient number of laps, these objects will have circular orbits embedded in the disk.This process is called gas capture or orbital damping.
The gas accretion and subsequent gas drag are significant only for prograde orbiters.Hence, we neglected any variation in the orbits of retrograde orbiters.
We defined the inclination damping timescale, t damp , for a BH of mass, m BH , on an initial orbit of semimajor axis4 , A, around a SMBH of mass, M SMBH , as (Wang et al. 2024) where Σ gas (A) is the surface density of the gas.In this model, we neglect the mass increase due to gas accretion.Hence, during the orbital evolution, m BH is a constant quantity.

Migration
Once a BH is embedded in the disk, it exchanges angular momentum with the surrounding gas and is subject to gas torques.Torques can be both positive or negative, leading to outward or inward migration, respectively.Similar to what happens to planet seeds in protoplanetary disks, migration can happen in two different ways called Type I and Type II.Small to medium-mass objects are subject to Type I migration, meaning that they change their radial position in the disk without significantly perturbing the density distribution of the disk itself.For a BH of mass m BH on a circular orbit with radius R, this happens on the timescale (Lyra et al. 2010;McKernan et al. 2012;Baruteau et al. 2014) where h is the aspect ratio of the disk and Ω is the Keplerian angular velocity around the SMBH.Differently from Eq. ( 9), here we are considering a radius, R, rather than a semimajor axis, A, because gas capture happens necessarily before migration5 , so the orbits have already been circularized when migration sets in.
In our disk model (Fig. 1), torques are positive in the inner region of the disk, where the slope of the surface density is positive, and they are negative in the outer region, where the slope of the surface density is negative (Bellovary et al. 2016).Therefore, Type I migration is directed outward in the inner disk and inward in the outer disk.At the location where the torques change sign, called a migration trap, migration will stall leading to a large accumulation of objects.Hence, after a timescale t migr, I (Eq.( 10)), the migrating object will be in the migration trap.
Larger objects, on the other hand, can open gaps in the disk.This happens because the motion of a massive object exerts an intense tidal perturbation on the disk, which effectively pushes material away from the orbit's trail (Bryden et al. 1999).This is called Type II migration.An object of mass m BH can open a gap in the disk if (McKernan et al. 2014) where q = m BH /M SMBH is the mass ratio with respect to the central SMBH, α is the viscosity parameter, and h = h (R) is the aspect ratio of the disk at radius R.
If an object opens a gap in the disk, assuming that no gas can cross the gap, its migration follows the viscous evolution of the disk's gas.Hence, the timescale for Type II migration is the timescale for the viscous evolution of the disk (McKernan et al. 2012): However, pressure forces in the disk push to close the gap.So, even if an object is massive enough to open a gap, such a gap can stay open against pressure forces only if (Bryden et al. 1999;McKernan et al. 2014) Type II migrators typically have high mass: taking as fiducial values α = 0.1 and h = 0.02, Eqs. ( 11) and ( 13) entail m BH 0.1 M SMBH .They are bound to their radial location in the disk and can only move with the disk on its viscous timescale.Hence, they will never reach the migration trap.Moreover, the gaps they create will prevent some Type I migrators from reaching the trap: they will intercept inward-moving migrators if they are located at a radius greater than the trap's, or they will intercept outward moving migrators if the opposite is true.These intercepted BHs can potentially pair up and merge with the Type II migrator, although their merger would not be assisted by GH.However, we expect such a massive BH in an AGN disk only in two cases: a high-generation hierarchical BH or the central BH of a dwarf galaxy dragged into the AGN disk after a galaxy-galaxy merger (e.g., Di Matteo et al. 2008).Describing galaxy mergers is out of the scope of this paper, while a highgeneration hierarchical BH can form only in the late stages of the disk's lifetime.We can assume that by that point most of the BHs have already migrated in the migration trap, and we can neglect any further pair-up event with a Type II migrator.
Therefore, in our model we consider the onset of Type II migration to be one of the processes that can halt hierarchical mergers.Moreover, the presence of gaps has noteworthy consequences on disk structure.Nevertheless, we bluntly neglect any evolution of the disk density profile.

Pair-up
Depending on the physical properties of the AGN (such as viscosity, gas density, aspect ratio, and SMBH mass), gas capture and migration can happen on short timescales.When these processes are efficient, they can lead to a large abundance of BHs in the narrow region of the migration trap.Also, all BHs in the migration traps are on similar orbits (prograde and quasi-Keplerian), so their relative velocities of encounter are small.Under such conditions, it is easy for two BHs to become gravitationally bound in a binary.Therefore, efficient damping and migration lead to efficient binary pair-up.
In this work, we assume that the pair-up of a primary and secondary BH is immediate as soon as the primary reaches the migration trap.This assumption is fully justified by the efficiency of dynamical friction in the migration trap (see Qian et al. 2024 and Sect. 4).The pairing timescale of a BBH is therefore where t in is the formation time of the primary BH since the time of formation of the disk.We considered each BBH to be in a circular orbit in the migration trap at a radius R trap (Sect.2.1).We computed the fraction, f trap , of BHs that reach the migration trap by counting the number of 1g BHs for which t pair < τ, where τ is the disk lifetime, and dividing it by the total number, N, of simulated firstgeneration BHs.We used this parameter for the computation of the maximum mass that can be accreted by a single BH, as in Eq. ( 7).

BBH properties
We determined the secondary BH mass as in Appendix C, and we set the secondary BH spin in the same way as for the primary BH.We set the secondary BH spin tilt as in Appendix B. We assigned the initial semimajor axis, a, of the binary sampling from a distribution p (a) ∝ a 9/2 (Binney & Tremaine 2008;Tagawa et al. 2020b) for a ∈ [a min , a max ], where a min = 1 R and a max is equal to the Hill radius, computed as where m 1 and m 2 are the primary and secondary BH mass, respectively.We set the initial eccentricity e following a thermal distribution p(e) ∝ 2e for e between 0 and 1 (Jeans 1919).
According to the Heggie (1975) law, a binary can survive in a star cluster only if it is hard, meaning that its binding energy is greater than the average kinetic energy of a field star where a is the semimajor axis of the binary, m * is the average mass of a star in the NSC, and σ is the three-dimensional velocity dispersion.Hence, we dynamically evolveed hard binaries only.

Orbital evolution and merger
Once the binary is in the migration trap, we set R = R trap and update the quantities in Eqs. ( 1)-( 3) accordingly.The semimajor axis, a, and eccentricity, e, of a binary with component masses m 1 and m 2 , embedded in the disk, evolve due to GH, irrespective of whether the binary is prograde or retrograde, as A51, page 5 of 20 (Ishibashi & Gröbner 2020) where µ = m 1 m 2 /(m 1 + m 2 ) is the reduced mass and Ω b = G(m 1 + m 2 )/a 3 is the Keplerian orbital frequency.
The binary also hardens due to the effect of GW emission, which will govern the evolution at small semimajor axes.The evolution of the semimajor axis, ȧGW , and eccentricity, ėGW , due to GW hardening proceeds as in Peters (1964).
The interaction with other objects also contributes to the hardening of hard binaries according to Heggie (1975).We neglected the hardening effect due to three-body interactions because they typically occur on a timescale longer than that of GH (Leigh et al. 2017, see discussion in Sect.4.3).The overall evolution of the binary is thus described by The equations for GH (Eqs.( 17) and ( 18)) are valid under the assumption that gravitational torques from the binary act axisymmetrically upon the disk so that they clear a cavity in the surrounding gas distribution, which remains circular throughout its inspiral.This is an idealized model since cavities in AGN disks can become eccentric (MacFadyen & Milosavljević 2008;Cimerman & Rafikov 2024) and can lead the orbital separation to grow in time (Miranda et al. 2016).Hence, we considered an "optimistic" model in which the evolution is given by Eq. ( 19) (hereafter, the GH model), and a "pessimistic" model in which we neglect GH and only consider the effect of GW emission (hereafter, the no-GH model).In both the optimistic and pessimistic cases, we integrated the hardening equations using the Euler method and an adaptive timestep (Mapelli et al. 2021).
We refer to the delay time between the pair-up and the merger as t del .If t pair + t del is longer than the lifetime of the disk, the disk has evaporated before the binary could merge.In this case, the BBH keeps hardening due to GW emission only.
We assumed that the BBH merges when its members cross the ISCO radius of a non-spinning BH with mass equal to the total mass of the binary system, r ISCO = 6 G (m 1 + m 2 )/c 2 , with a tolerance of 0.1 r ISCO .This happens on a merger timescale We modeled the mass and spin of the merger remnant using fitting formulas from numerical relativity, as described by Jiménez-Forteza et al. (2017).At birth, merger remnants receive a relativistic kick v kick because of the transfer of linear momentum caused by asymmetries in GW emission.We used the model of Maggiore (2018, Eq. (14.202)) for the magnitude and direction of the kick velocity, u kick .
The relativistic kick u kick pushes the merger remnant out of the migration trap.We computed the new velocity as u fin = u Kepl (R trap ) + u kick , where the Keplerian velocity in the migration trap is computed accounting for the mass of the SMBH M SMBH and the inner part of the NSC M NSC , while neglecting the mass of the gas disk: The final semimajor axis A fin of the remnant orbit, computed by means of simple orbital transfer calculations (Hohmann 1960), is The quantities in Eqs. ( 1)-(3) were updated accordingly.As a safety check, we ensured that the new orbital semimajor axis, A fin , is smaller than the maximum radius of the disk, R max , meaning that the remnant can experience damping and become embedded in the disk.Otherwise, we did not consider the remnant for future generations.

Nth-generation (Ng) BHs
A seed BH can only go through a finite number of hierarchical mergers before it comes across one of the following scenarios: 1.The disk has evaporated.Therefore, damping, migration, and any other effects due to gas torques stop.We evaluated this by checking if, at generation N, the merger timescale, t (N) merg , is shorter than the disk lifetime, τ.In our formalism, we defined the merger timescale of a Ng BH as t (N) merg = t (N)  in + t damp + t migr, I + t del (N) , where t (N) in is the evolutionary time of the previous generations.2. The relativistic kick received at merger is so strong that the remnant is ejected from the AGN.We computed the escape velocity considering only the gravitational potential of the SMBH and the inner NSC, while neglecting the mass of the gaseous disk, as Then we computed the final velocity after kick u fin = u in + u kick and ensured that it is smaller than the escape velocity, v esc .We also ensured that the merger remnant is on a diskcrossing orbit by comparing its semimajor axis, A fin , with the outermost disk radius, R max .3. The number of BHs in the NSC is finite.Therefore, the mass that can be accreted is limited and the BH may not find any more companions to pair up with.We kept track of the mass accreted by a single BH, which is the sum of its initial mass, m 1 , and of the masses of all the secondaries, m (N) 2 , it pairs up and merges with where the index N represents the generation number.
We ensure that the BH does not accrete more mass than what is available in the inner NSC in the form of other BHs, namely M acc ≤ M max BH , where M max BH is obtained as in Eq. ( 7). 4. The BH is so massive that it opens a gap in the disk and can only move from its radial location due to Type II migration.We checked whether both conditions in Eqs. ( 11) and ( 13) are respected.For typical values of viscosity and aspect ratio, these conditions entail m BH 10 −1 M SMBH .If any of conditions (1)-( 4) was met, we did not consider the merger remnant for future generations.We followed the evolution of Nth generation BHs with a procedure similar to the one outlined for first-generation BHs.In hierarchical mergers, A51, page 6 of 20 the remnant of an (N − 1)th-generation merger acts as the primary BH for the Nth-generation.So, the primary mass and spin are simply set as the remnant mass and spin of the previous generation, computed according to Jiménez-Forteza et al. (2017).Similarly, the initial position of the primary component of the Nth-generation BBH is set as the position of the merger remnant of the previous generation, set in Eq. ( 22).This value is used to compute the pairing time t pair = t damp + t migr, I + t in , where for Nth generations t (N)  in = t (N−1) merg .The migration and pair-up physics are the same as described for the first generation.
We modeled the secondary component of each Nthgeneration BBH as described in Appendix C. We define an Nthgeneration (Ng) BH to be the result of the repeated merger of N stellar-origin BHs.For instance, an Ng BH can be the result of either an Mg − 1g merger (where M + 1 = N) or an Mg − Lg merger (where M + L = N, L > 1).We kept iterating the same procedure until the considered BH met at least one of requirements (1)-( 4).

Cosmic evolution
We modeled the cosmic evolution of AGNs using data from the IllustrisTNG100, a magneto-hydrodynamical cosmological simulation adopting a cubic box of size 110.7 Mpc 3 with a resolution of roughly 6 Kpc 3 (see Pillepich et al. 2018 andSpringel et al. 2018 for further details).
We calibrated the AGN density distribution n AGN (z) in the IllustrisTNG data by choosing a threshold in the mass accretion rate such that the value at redshift zero, n AGN (0), is compatible with the observational value of ∼6 × 10 −6 Mpc −3 obtained from a sample of X-ray-selected AGNs (Buchner et al. 2015).As displayed in Fig. 2, we find that setting this threshold to a fifth of the Eddington mass accretion rate, namely ṀSMBH ≥ 0.2 ṀEdd , provides an appropriate normalization.We employed n AGN (z) to compute the BBH merger rate properties in various redshift bins as explained in Appendix D.
The SMBH mass distribution is expected to evolve as a function of redshift (Weinberger et al. 2018).Nevertheless, since the mass evolution cannot be constrained by current data, we used the observational distribution from Greene & Ho (2007) at every redshift for simplicity.

Description of runs
We considered two different models for the AGN channel: with and without GH (hereafter, GH and no-GH).We ran N = 10 5 realizations of our two models, each with different SMBH mass and AGN lifetime randomly extracted as described at the end of Sect.2.1.In each run, we simulated all the BHs that reach the migration trap within the AGN disk lifetime.Furthermore, we simulated four other channels (isolated, YSC, GC, and NSC) with the same code (Fastcluster; Mapelli et al. 2021Mapelli et al. , 2022) ) and using the same underlying initial conditions, such as the stellar evolution model determining the 1g BHs mass distribution (Appendix A).This allowed us to filter out any bias that might arise by using different numerical codes for different environments: the differences we see are not due to the numerical approach adopted but to intrinsic differences among channels.We performed a Bayesian analysis to determine the AGN mixing fraction as described in Appendix E.

Impact of gas hardening
We find that the efficiency of the hierarchical merger process is significantly higher in the presence of GH: in our GH model, seed BHs can go through up to roughly 500 merger episodes, whereas in the no-GH case the hierarchical chain usually stops after a few generations (Fig. 3).In the GH (no-GH) model, the hierarchical merger chain stops because of the condition on timescales, maximum mass, Type II migration, or ejection in the A51, page 7 of 20 58.7% (81.6%), 27.1% (18.0%), 13.6% (0%), and 0.6% (0.4%) of simulations, respectively.Figure 4 shows the main properties of dynamically assembled BBHs in our AGN simulations.We display only BBHs that merge within a Hubble time.In both models, hierarchical mergers in the migration trap are successful in producing BBHs with primary mass in the pair-instability mass gap and above, but there is a major difference between the GH and no-GH scenarios: the primary BH mass extends only up to 300 M in the no-GH scenario, while it reaches ≈5 × 10 3 M in the GH scenario, because GH dramatically increases the efficiency of hierarchical mergers.Hence, the GH mechanism produces a recognizable feature in the high-mass end (m 1 100 M ) of the mass spectrum.This is caused by the stark difference in delay timescales (Fig. 4f) due to the different GH prescriptions: in the GH scenario, the relative distribution has a pronounced peak between 1 yr and 1 Myr, which in the no-GH case is absent.Information on other relevant timescales is illustrated in Sect.3.4.Secondary BHs can also have masses in the pair-instability mass gap but only up to a few times 10 2 M in the GH scenario, as the AGN channel tends to favor mergers with low mass ratio (Fig. 4b).
The GH model has a sharp peak in the distribution of the primary spin magnitude at χ 1 = 1 as well as a peak in the distribution of the secondary spin magnitude at χ 2 0.9 associated with high-generation mergers (Fig. 4c,d).Both models GH and no-GH also have a peak at χ 1 ≈ 0.75 (also χ 2 ≈ 0.75 for the GH case), corresponding to the second generation.
The distribution of spin-tilt angles is influenced by GH, as shown in Fig. 5 and explained in Appendix B. In the no-GH scenario, the spin tilts θ 1,2 are isotropically distributed with respect to the orbital angular momentum L. Instead, in the presence of GH, χ 1 , χ 2 , and L are all preferentially aligned with each other.The distribution of the secondary spin tilt θ 2 is slightly wider than the primary's spin tilt, because of a cumulative effect: the misalignment between χ 2 and L is set keeping into account the misalignment between χ 1 and L (see Appendix B).
Both scenarios preferentially produce BBHs with low eccentricity when6 the GW frequency is 10Hz.Nevertheless, they can produce BBHs with eccentricity e ≥ 10 −2 (Fig. 4e), which is potentially detectable by the LIGO-Virgo-KAGRA (LVK) interferometers (Romero-Shaw et al. 2021).In our GH model, there are two processes at play: GH pumps the eccentricity (Eq.( 18)), whereas GW emission damps it (Peters 1964).Hence, the GH scenario is more likely to produce BBHs with eccentricity in the range e ∈ [10 −5 , 10 −2 ] for f GW = 10 Hz than the no-GH one.In a more accurate model, we should also include the effect of three-body scatterings (Samsing et al. 2022) and tidal forces exerted by the SMBH (Rom et al. 2024), both of which pump BBH eccentricity.
Figure 6 shows the relation between the BBH mass m 1 + m 2 and mass ratio m 2 /m 1 in our models, compared to GW data (Abbott et al. 2023a,b).The AGN channel tends to favor mergers with low mass ratio (Appendix C).This effect is particularly important for high primary BH masses m 1 10 2 M when accounting for GH.
The effective spin distribution in the GH model (Fig. 7) peaks at χ eff ∼ 1, corresponding to maximum alignment, and shows lower peaks at χ eff 0.2 and χ eff 0.5.From Fig. 8, we see that 1g mergers populate the peak at χ eff 0.2, 2g mergers create a peak at χ eff 0.5, while third-and higher-generation BBHs contribute the main peak at χ eff 1.The precession spin distribution, instead, has no sharp features in this scenario.
In the no-GH scenario, the orbital angular momentum is isotropically oriented with respect to the AGN disk (Appendix B).This produces a bell-shaped effective spin distribution centered on χ eff = 0 with half-width at half-maximum HWHM 0.3 and a precession spin distribution that peaks at χ p 0.2 and χ p 0.75 (Fig. 7).As the BBH generation increases, so does the magnitude of its BH spins, making the corresponding effective spin distribution progressively wider as shown in Fig. 8.The peaks at χ p 0.2 and 0.7 are populated mostly by 1g and 2g mergers, respectively, whereas the tail at χ p > 0.8 results from higher-generation mergers.

Anticorrelation between q and χ eff
Figure 9 shows the relationship of the effective spin with primary mass and mass ratio of BBH mergers in our simulations.In the GH scenario, there is a clear correlation between effective spin and BBH mass, as well as a clear anticorrelation between effective spin χ eff and mass ratio q.This is a possible explanation for the anticorrelation between χ eff and q found in the LVK data (Callister et al. 2021;Santini et al. 2023).If this anticorrelation stems from hierarchical mergers in the GH regime, our simulations suggest that it should extend to lower mass ratios and higher effective spins than currently observed by LIGO and Virgo.
The anticorrelation is particularly noticeable for χ eff 0.4 because all higher-generation mergers display the following features: high mass, small mass ratio, and high effective spin.In the no-GH model, both the spin alignment and the hierarchical merger process are suppressed; so there is no clear correlation between the aforementioned quantities.

Comparison with other channels
Figure 10 compares the main properties of BBH mergers in five different formation channels, namely AGN disks, NSCs, GCs, YSCs, and isolated binary evolution (iso).The simulations for NSCs, GCs, YSCs, and iso adopt the same setup as model B of Mapelli et al. (2022), and use as initial conditions catalogs of BH masses derived with Sevn (Iorio et al. 2023), for consistency with the catalogs of the AGN channel.We provide more details about NSC, GC, YSC, and iso simulations in Appendix A.
All dynamical formation channels can produce BBH mergers with masses and spin magnitudes higher than the isolated channel.The no-GH AGN disk scenario resembles most closely the results of the other dynamical channels, whereas the GH AGN disk channel is completely different from the others.
The population of dynamical channels appears as follows: the mass distributions show one or two peaks at a few tens of solar masses and a tail that extends up to a few hundred solar masses; the effective spin distribution is symmetric, peaks at χ eff = 0, and has a half width at half maximum of ∼0.5 (∼0.3 for the YSC channel); and the precession spin distribution has two peaks, at χ p 0.2 and χ p 0.7.This contrasts with the isolated channel population, which features a primary (secondary) BH mass distribution that does not extend any higher than 50 M (40 M ).
The primary BH mass distribution for the GH AGN channel, instead, extends up to ∼5 × 10 3 M .The effective spin distribution has two sharp peaks at χ eff 0.2 and χ eff 1, as well as a minor peak at χ eff 0.5, whereas the precession spin distribution has no evident peaks.
We computed the mixing fractions, f mix , for our five channels, as described in Appendix E. To derive f mix , we marginalized over the merger rate density, to avoid this extremely uncertain quantity affecting our results.As shown in Fig. 11, the NSC and YSC channels are associated with larger median mixing fractions than the other channels.This happens because, in our models, the NSC channel can account for the peak at low BH mass (m 1 ≈ 10 M ) found in the LVK data, while the YSC channel contributes mostly to the high-mass peak at ∼35 M (Abbott et al. 2023b).In fact, in our model, the YSC scenario produces larger median BBH masses than the NSC scenario, because the escape velocity from YSCs is much lower, preventing the low-mass BHs from merging in YSCs (they are ejected by supernova kicks, see the discussion in Mapelli et al. 2021).The isolated channel also produces BBH mergers with a peak of the primary BH mass at ∼10 M , but has too little support for zero and negative values of χ eff with respect to the observed ones.The mass and spin properties of the GC scenario are intermediate between YSCs and NSCs.
There is a small difference between the mixing-fraction distribution of the GH and no-GH scenarios because this analysis only considers detectable events: the most massive BHs of the GH population (which are the main difference with respect to the no-GH model) have no impact on f mix .Overall, our results A51, page 9 of 20 Fig. 6.Probability density function (pdf) of the binary mass and mass ratio for all BBH mergers in our simulations (blue color palette).We do not apply any selection effects to our simulations.Left-hand (right-hand) panels: GH (non-GH) scenario.The magenta dots show the median values of the parameters for LVK BBH merger event candidates with p astro > 0.9 from GWTC-3 (Abbott et al. 2023a,b).We do not include error bars for readability purposes.We do not apply any selection effects to our simulations.Solid black lines show the posterior distribution inferred from LVK data, after applying a parametric-model description of the intrinsic BBH population with hierarchical Bayesian inference (Abbott et al. 2023b).
confirm that the current LVK sample of BBH mergers is too small and the uncertainties on theoretical models are too large to draw an informative mixing-fraction analysis with five channels (see Sect. 3.4 of Mapelli et al. 2022).

Timescales
Figure 12 shows the distribution of all relevant timescales with and without GH.Both the damping time t damp and the migration time t migr, I span a large range, which is representative of the large range in SMBH mass and initial position R.The damping timescale for the first hierarchical merger generation has a flat distribution in the range 10 −4 −10 2 Myr, while the migration timescale ranges between 10 −2 and 10 3 Myr.In the GH scenario, t damp and t migr, I display additional peaks for subsequent generations at 10 −5 Myr and 10 −2 Myr, respectively.The delay timescale t del spans a large range of values, because of the large range encompassed in BH masses and initial BBH semimajor axes.In the GH scenario, it has a pronounced peak at ∼1 yr, whereas in the no-GH case it is typically as large as the Hubble timescale (∼14 Gyr).Indeed, the longest timescale in the no-GH scenario is the delay timescale.In the GH scenario, instead, the evolution is overall governed by the migration timescale.We also display the timescale for gas dynamical friction τ 0 , responsible for BBH pair-up, which we discuss in Sect.4.2.It is typically on the order of 10 −13 Myr and can be as short as 10 −20 Myr for the disk parameters in our model; hence, it is negligible compared to other timescales at play.

Outlook for ET and LISA
We studied the evolution of hierarchical mergers in the migration traps of AGN disks and explored the role of GH.Our models, especially the GH scenario, predict the formation of merger events with higher mass than is detectable by current ground-based detectors.Indeed, the frequency emitted by a BBH A51, page 10 of 20 Fig. 8. Effective and precession spin distributions for the GH (left-hand panels) and no-GH (right-hand panels) models, where the shade represents different hierarchical merger generations.The lightest hue refers to the first generation (stellar-origin BHs), while the darkest hue refers to generations higher than the third.
depends on its mass and steadily increases during its inspiral.The frequency of the ISCO is equal to (Maggiore 2008) Figure 13 shows the probability distribution function of the maximum emitted frequency by AGN BBHs.The LVK interferometers are only sensitive in the frequency range from a few tens of hertz up to ∼1 KHz (Abbott et al. 2023a).Hence, only a fraction of our predicted merger events in the GH scenario are observable with existing detectors.
The next generation ground-based interferometers, the Einstein Telescope (ET) and Cosmic Explorer (Punturo et al. 2010;Maggiore et al. 2020;Evans et al. 2021;Branchesi et al. 2023), and space-borne interferometers Laser Interferometer Space Antenna (LISA), Deci-hertz Interferometer GW Observatory (DECIGO) and TianQin (Amaro Seoane et al. 2013, 2017;Luo et al. 2016;Kawamura et al. 2019) will be able to detect GW signals with lower frequency than currently possible.For example, the frequency range of the ET and Cosmic Explorer (≈1−10 4 Hz) has substantial overlap with the GW frequency of BBH mergers from the GH AGN channel, with the exception of the very high-mass tail.LISA will instead be sensitive to frequencies lower than ∼1 Hz (Robson et al. 2019), so it would be able to detect the highest-mass end of the synthetic mergers predicted by our GH AGN model.We will explore detectability by ET and LISA in detail in a follow-up work.

BBH pair-up
In our model, we assume that the pairing of a primary and a secondary BH is immediate as soon as the primary reaches the migration trap.A realistic model for BBH pair-up in AGN disks should keep into account GW two-body capture, threebody encounters, and gas dissipation.Whitehead et al. (2023) run hydro-dynamical simulations of close encounters between pairs of BHs embedded in the gaseous disk and find that dissipation by gas gravitation is not always efficient for the formation of bound BBHs.Specifically, they find that it is usually an effective mechanism for BBH formation for gas densities ρ g in the range −4.5 log ρ g R 3 H / (m 1 + m 2 ) −2.5, where R H is the Hill radius of the BBH.In our simulations, we typically have log ρ g R 3 H / (m 1 + m 2 ) −6 for our BBHs because the high gas density in the migration trap is compensated for by small Hill radii due to the proximity to the SMBH.Hence, gas dissipation is expected to be inefficient.
On the other hand, in their recent N-body simulations, Qian et al. (2024) find that BBH formation from two BHs on similar orbits embedded in the disk is efficient for Ω τ 0 10, where Ω is the Keplerian angular velocity of the BHs and is the dynamical friction timescale (Ostriker 1999), which, as shown in Qian et al. (2024), is linearly related to the timescale of BBH pair-up.In our model typically Ω τ 0 ∼ 10 −4 −10 −3 due to the high gas density ρ g in the migration trap, which justifies our assumption.Moreover, we expect a high number density n BH of BHs in the migration trap due to efficient migration.This may further aid BBH pair-up because the timescales for two-body and threebody captures scale respectively as n −1 BH (Quinlan & Shapiro 1990) and n −2 BH (Fragione & Silk 2020).In future work, we will further explore the process that leads to BBH pair-up, accounting for these additional effects.
In this work we neglect BBH formation in the bulk (i.e., outside of the migration trap), which is expected to be efficient due A51, page 11 of 20 to the high number of Type I migrators embedded in the disk (e.g.Tagawa et al. 2020a,b).Interestingly, gas-assisted bulkassembled BBHs migrate toward the migration trap while hardening (Tagawa et al. 2020b, Fig. 8).Hence, some binaries may merge in the migration trap even if they were assembled in the bulk.

Three-body encounters
We neglect three-body interactions and their effects on BH scattering and BBH hardening.We make this approximation because of the dearth of semi-analytical models for these interactions: in the literature, there are some works (e.g., Coleman Miller & Lauburg 2009;Fragione & Silk 2020) that assume an isotropic distribution of velocities and are appropriate for spherical star clusters.In a Keplerian disk geometry, the distribution of velocities is much different and these models are not appropriate (McKernan et al. 2022).Nevertheless, a hard binary hardens via binary-single encounters (Heggie 1975), and this has a twofold effect: the inspiral of BBHs is accelerated via three-body effects, and the third intruding body receives a recoil kick that may prevent it from reaching the migration trap.Also, binary-single interactions are expected to increase the eccentricity of all BBHs to e (10 Hz) ≥ 10 −4 and potentially flip the A51, page 12 of 20 orientation of the BBH orbital angular momentum, leaving key signatures in the effective spin distribution (Samsing et al. 2022).Additionally, Rowan et al. (2023) show that the initial eccentricity of BBHs formed in AGN disks is high, which is not well modeled by our assumption of a thermal distribution (Jeans 1919).Tagawa et al. (2020b) simulated the evolution of the compact object population in AGN disks using a one-dimensional N-body simulation combined with a semi-analytical model for the formation, disruption, and evolution of binaries.They include the effects of three-body interactions and employ a Thompson et al. (2005) disk model rather than the SG model that we use.Nevertheless, comparing our Fig.6 with Tagawa et al. (2020b, Fig. 12a), we can point out that the output from our fiducial model is compatible with the output from theirs.This suggests that three-body effects have a relatively minor impact on the population of merging BBHs.Leigh et al. (2017) find that the approximate encounter timescale between a BH in the migration trap and an intruder is a fraction of the Type I migration timescale (Eq.( 10)) of an object of mass m * starting its migration from the outer skirts of the disk, namely where N * is the number of Type I migrators in the disk.
Assuming m * = 1 M and N * = 100, we find that the encounter timescale is typically larger than the delay timescale A51, page 13 of 20 10 13 10 10 10 7 10 4 10 1 10 2 10 5  in the GH scenario (Fig. 14).This justifies our disregarding of three-body effects in this scenario.Nevertheless, they should be included in the no-GH model.
This implies that Heggie's law (employed in Eq. ( 16)) may not be valid in the GH case, as even soft binaries may be hardened considerably on a timescale shorter than the one on threebody interactions.We find that forming soft BBHs in our model is an extremely rare event (fewer than 1 in 10 000 BBHs), so including them in our computation would not have affected the results significantly.
Furthermore, efficient damping, migration, and BBH pairup in AGN disks are expected to give rise to a high concentration of BBHs in migration traps.Hence, we should take into account the effect of binary-binary interactions.Such encounters are expected to efficiently ionize one of the binaries involved and lead to the formation of a stable triple (Marín Pina & Gieles 2024).

Position of migration traps
The position of migration traps is highly uncertain: it strongly depends on the disk model and the migration prescription used.For instance, Pan & Yang (2021) compute Type I migration torques in the locally isothermal approximation as well as torques due to winds and to the gravitational interaction with the SMBH for three different disk models.In their approximation, no disk model develops migration traps.Moreover, they find that in the inner part of the disk (R 10 2 R g ) Type I migration torques are always overpowered by the other effects, which prompt embedded BHs to eventually merge with the SMBH.This justifies our assumption that we can neglect the inner migration trap for the SG model (Sect.2.1), but it implies that we somewhat overestimate the number of Type I migrators.
In their recent work, Grishin et al. (2024) account for both Type I migration, caused by the gravitational perturbation of embedded BHs in the disk, and thermal migration, caused by the thermal response of the disk to the small and overdense accretion disks surrounding embedded BHs.They find that the resulting migration trap position is at much larger distances than previously identified by considering Type I migration only.
A51, page 14 of 20 In our fiducial model, we followed the prescription by Bellovary et al. (2016) and identified the position of the migration trap as R trap = 1324 R g , whereas Grishin et al. (2024) find that the updated position of the migration trap is typically larger R trap = 10 3 −10 6 R g , and has a steep dependence on the SMBH mass as log R trap /R g − log (M SMBH /M ) + 11, for M SMBH 10 8 M .Hence, for a typical SMBH of mass 10 6 M , the migration trap is at ∼10 5 R g .
Migration is also affected by the additional thermal torque, which reduces the migration timescale by a factor of where t migr, I is defined as in Eq. ( 10).
Figure 15 shows a comparison between the results of our fiducial model and the results we obtain assuming the position of migration traps and the migration timescales predicted by Grishin et al. (2024).We find that the hierarchical merger process is suppressed when including a treatment for thermal torques, as only a handful of seed BHs reach the second generation in our simulations.This happens because the gaseous disk is significantly thicker and more diluted in the outer area of the disk where the new migration trap is located (see Sect. 2.1).Therefore, as the migration timescale is proportional to h 3 and Σ −1 gas , the migration process is much slower in this scenario.Moreover, even when BHs successfully reach the migration trap, the delay timescale ends up being too long because of the low gas surface density.

General caveats
In our model, we disregard dynamical interactions with the SMBH and gravitational perturbations caused by intermediate-mass BHs.For example, Deme et al. (2020) show that the presence of intermediate-mass BHs in the disk may enhance the ionization of BBHs and consequently decrease the merger rate.
Moreover, in Sect.2.3.1 we assumed that eccentricity and inclination of a BH orbit are damped on similar timescales due to gas torques.Wang et al. (2024, Fig. 4) show that this is not always accurate: if the initial orbit is highly eccentric (e 0.9), the eccentricity-damping timescale can be up to five times longer than the inclination-damping one.A BH on a gas-embedded eccentric orbit (i ∼ 0, e h) would be subject to spin-down (McKernan & Ford 2023), which would leave an imprint on the expected spin distribution.
Finally, we ignored the evolution of AGN disks in time.This can happen slowly over the AGN lifetime as BHs are embedded in the disk (Tagawa et al. 2022) and gas is accreted by the SMBH.The efficiency of all dynamical processes strongly depends on disk density and aspect ratio.If these quantities evolve over time, the resulting BBH population will be affected as well.

GW190521 and other transient events
It has previously been suggested that the transient events GW170729 and GW190521 could have originated in AGN disks (Yang et al. 2019a;Tagawa et al. 2021;Samsing et al. 2022;Graham et al. 2020;Morton et al. 2023).In particular, Yang et al. (2019a) found that it is five times more likely for GW170729 to arise from hierarchical mergers in AGNs than assuming that all events in the GWTC-1 catalog (Abbott et al. 2019) arise from the same channel, whereas Morton et al. (2023) show that the association between GW190521 and the AGN flare ZTF19abanrhr (Graham et al. 2020) is highly preferred over the lack of association, suggesting that the GW transient event was generated in an AGN.
A51, page 15 of 20 We compare the posterior contours for the primary mass and effective spins of such events with the output of our AGN simulations, as shown in Fig. 9.We also include the posterior contour of the transient event GW190403_051519, although it has a low S/N of 7.6 +0.6  −1.1 and a low probability of having an astrophysical origin, p astro = 0.61 (Abbott et al. 2021).
We find that the contours of GW170729 have some overlap with both the GH and the no-GH AGN models, and thus might be compatible with an origin in an AGN environment.In contrast, GW190521 has no significant overlap with our models in the χ eff − m 1 space, and negligible overlap with our no-GH model in the χ eff − χ p space.Indeed, the posterior probability distribution of GW190521 shows moderate support for high precession spin and low effective spin, suggesting that its BH spin vectors are large but misaligned.In contrast, in our GH model, we predict that they should be aligned for a BBH of primary mass ∼100 M .This does not rule out the hypothesis of an AGN origin for GW190521, as we might speculate that the BBH may have been perturbed by three-body encounters, which we do not account for in our GH model (Appendix B).This would alter the orientation of the orbital angular momentum and decrease the effective spin, as well as increase the BBH eccentricity (Samsing et al. 2022).
Finally, the posterior contours of GW190403_051519 significantly overlap with our GH AGN model.However, this event candidate has a high false alarm rate and may not have astrophysical origin (Abbott et al. 2021).

Summary
We explored the formation of BBH mergers in AGNs employing a new semi-analytical model.our key findings are as follows: -The presence of GH significantly increases the efficiency of hierarchical mergers in AGN disks, allowing seed BHs to go through up to 500 merger episodes.This leads to the formation of BBHs with high masses (up to a few thousand solar masses) and low mass ratios (q 10 −2 ).-In contrast, if GH is not efficient, the hierarchical merger chain is truncated after a few generations, leaving a BBH population with no components more massive than ∼10 2 M .-The distribution of spin tilt angles is influenced by GH.
There is a preferential alignment of spins in the GH scenario and isotropic distribution in the no-GH scenario.Hence, GH leads to distinct features in the effective spin distributions: a main peak on χ eff 1 corresponding to maximum alignment and smaller peaks at χ eff 0.2 and χ eff 0.5.-In the GH scenario, we find an anticorrelation between q and χ eff that might even extend to lower values of q and higher values of χ eff than currently observed by LVK.-Comparison with other formation channels shows that AGNdriven mergers in the GH scenario result in higher primary BH masses and higher effective spins.
In summary, we find that efficient GH in AGN disks enhances the formation of BBH mergers with high primary masses and low mass ratios, and leads to a strong preference for χ eff ≈ 1.Given the high masses of BHs, next-generation ground-based detectors such as the ET are an ideal test bed for such unique features.
for M ≤ N, so that a 1g primary BH (N = 1) will necessarily pair up with a 1g secondary (M = 1).For M = 1, the secondary BH mass m 2 is randomly drawn with probability distribution = m 1 .For M > 1, we determined the secondary mass as follows.First, we generated a 1g BH determining its mass and spin as in Eq.C.2.Then, we let it go through a certain number of Ng − 1g mergers 7 until it created an Mg remnant.
At each step, the primary will be the remnant of a previous merger event and its mass, m 1 , and spin, χ With this choice, if the primary has mass compatible with a 1g BH (m 1 ≤ m max 1 ), we sample the secondary in Eq. (C.2) (model from O' Leary et al. 2016).
Otherwise, if the primary is a higher-generation BH (m 1 > m max 1 ), we keep the same analytical description as in eq.C.2 but we force the secondary to have mass compatible with a 1g BH (m 2 ≤ m max 1 ).This calculation is quite fast because we only compute the remnant mass and spin at each step, neglecting the merger times.

Fig. 1 .
Fig.1.Surface density, Σ gas , aspect ratio, h and sound speed, c s , profiles as a function of the scale distance, R/R g .Dashed maroon lines represent the SG model for a disk with viscosity α = 0.01 around a 10 8 M SMBH.Solid blue lines show our broken power-law best fit to the SG model for α = 0.01 and M SMBH = 10 8 M .The dash-dotted green and dotted orange lines show our fits for a 10 7 M and a 10 6 M SMBH, respectively (both with α = 0.01).

Fig. 3 .
Fig.3.Number of mergers and characteristic masses at each merger generation.The blue histogram and the left-hand y-axis show the number of BBH merger events for each hierarchical merger generation, N g .The orange (green) shaded area and the upper right-hand y-axis show the 25%-75% percentile of the chirp mass (primary mass) for merging BBHs of each generation.The orange (green) dots and lower right-hand y-axis show the average chirp mass (primary mass) for merging BBHs of each generation.

Fig. 4 .
Fig. 4. Main properties of dynamical BBH mergers in our GH and no-GH AGN models.The unfilled solid blue (dashed orange) histogram shows BBHs in the GH (no-GH) AGN scenario.The filled blue (orange) histogram shows the 1g BBHs in the GH (no-GH) AGN scenario.If 1g BHs in the no-GH scenario are identical to the GH scenario, they are not shown.Panels (a) and (b): probability density function (pdf) of the primary and secondary BH mass (m 1 and m 2 ).Panels (c) and (d): primary and secondary spin magnitudes (χ 1 and χ 2 ).Panel (e): orbital BBH eccentricity e(10 Hz) when the GW frequency is f GW = 10 Hz.In gray we show the thermal distribution p(e) = 2e for comparison.Panel (f): timescales of delay time, t del , between the BBH pair-up and the merger.

Fig. 5 .
Fig. 5. Probability density function (pdf) of spin tilt angles θ 1 and θ 2 .The blue (green) histogram shows the primary (secondary) spin tilt in the GH scenario.The unfilled orange histogram shows the primary and secondary tilts in the no-GH scenario.The first generation is not shown separately in this plot because the distributions are identical at each generation.

Fig. 7 .
Fig. 7. Effective (left) and precession spin (right) probability density functions (pdf) GH (solid blue lines) and no-GH (dashed orange lines) model.We do not apply any selection effects to our simulations.Solid black lines show the posterior distribution inferred from LVK data, after applying a parametric-model description of the intrinsic BBH population with hierarchical Bayesian inference(Abbott et al. 2023b).

Fig. 9 .
Fig. 9. Effective spin correlations with relevant BBH parameters.The top row shows the probability density distribution for the primary BH mass, m 1 , and the effective spin, χ eff , compared with posterior contour plots at credibility levels 50% and 90% of LVK BBH merger events GW190521 (dash-dotted red line, Abbott et al. 2020) and GW170729 (solid orange line, Abbott et al. 2023b), and LVK transient event GW190403_051519 (dashed yellow line, Abbott et al. 2021).The central (bottom) row has the same legend as the first row, but shows the mass ratio, q = m 2 /m 1 , (precession spin, χ p ) and the effective spin, χ eff .

Fig. 10 .Fig. 11 .
Fig.10.Probability density distribution (pdf) of the primary and secondary BH masses, m 1,2 , the effective spin, χ eff , and the precession spin, χ p , in the five channels modeled with Fastcluster: active galactic nuclei (AGNs) with GH (solid dark red line) or with no-GH (dashed bright red line), NSCs (dark gray line), GCs (light gray line), YSCs (dark gold line), and isolated binary evolution (light gold line).

Fig. 12 .Fig. 13 .
Fig. 12. Relevant timescales in our model.The black line shows the overall merger timescale, t merg .The dash-dotted green line shows the gas capture timescale, t damp , (the shaded green histogram shows t damp for 1g BBHs only).The dotted navy line shoiws the Type I migration timescale, t migr, I , (the shaded navy histogram shows t migr, I for 1g BBHs only).The dashed pink line shows the delay timescale, t del , between the pair-up and the merger (the shaded pink histogram shows t del for 1g BBHs only).The light blue line shows the gas dynamical friction timescale, τ 0 , related to BBH pairing (Qian et al. 2024).

Fig. 14 .
Fig.14.Delay timescale, t del , in the GH model compared against the three-body encounter timescale, t enc , estimated as inLeigh et al. (2017).The red line is t del = enc .

Fig. 15 .
Fig. 15.Main properties of dynamical BBH mergers with two different models for the migration timescale and the location of migration traps.The solid blue line shows the fiducial GH model.The dashed red line shows the results of runs where we used the model by Grishin et al. (2024) for both the migration timescale and the location of the migration trap.Filled light blue histograms refer to 1g mergers in the fiducial model.Panels (a) and (b): primary and secondary BH masses.Panels (c) and (d): primary and secondary BH spin magnitudes.
(O'Leary et al. 2016)p (m 2 | m 1 ) ∝ (m 1 + m 2 1 , are computed according to Jiménez-Forteza et al. (2017).The secondary BH mass is sampled from p(m 2 | m 1 ) ∝ (m 1 + m max 2 on the value of m 1 .In particular, if m 1 is smaller than the maximum mass m max 1 of a 1g BH in the input sample coming from the population synthesis code sevn, the secondary mass cannot exceed the mass of the primary: m max 2 = m 1 .Otherwise, we set m max 2 = m max 1 .This is a modification of Eq. (C.2).