The LVK collaboration has detected merging BBHs with masses extending into the pair-instability mass gap (
$\gtrsim 60\,M_\odot$
), suggesting a dynamical origin (Abbott et al. Reference Abbott, Abbott and Acernese2023). One promising environment for such mergers is the accretion disk surrounding a supermassive black hole (SMBH) in an AGN. Gas-assisted dynamics favor a high encounter rate between BHs, enhancing the likelihood of binary formation. Moreover, remnants of previous mergers are retained thanks to the deep gravitational potential well of the SMBH, leading to multiple generations of mergers (e.g. Tagawa et al. Reference Tagawa, Haiman and Kocsis2020).
We present a fast, flexible framework for simulating BBH mergers in AGN disks that reproduces key features of hierarchical mergers and can probe the origin of extreme-mass, high-spin BBHs. Moreover, the redshift evolution of AGN-driven merger rates can shed light on SMBH growth across cosmic time.
1. Semi-Analytical Model for BBH Mergers in AGN Disks
We model BBH formation and mergers in AGN disks by integrating the relevant physical processes into a consistent population-synthesis framework, based on the fastcluster code (Mapelli et al. Reference Mapelli, Dall’Amico and Bouffanais2021) extended with AGN-specific prescriptions (Vaccaro et al. Reference Vaccaro, Mapelli and Périgois2024). The AGN disk is described using the steady-state model of Sirko & Goodman (Reference Sirko and Goodman2003), which provides radial profiles of gas surface density
$\Sigma_\mathrm{g}$
and aspect ratio h. These profiles are computed numerically with the pAGN code (Gangardt et al. Reference Gangardt, Trani and Bonnerot2024) for given values of the SMBH mass
$M_\bullet$
and the disk viscosity parameter
$\alpha$
. We consider stellar-origin BHs formed in a coexisting nuclear star cluster, initialized at solar metallicity using the SEvN code (Iorio et al. Reference Iorio, Mapelli and Costa2023). BHs are distributed uniformly in radius with isotropic initial inclinations. We follow the evolution of BHs individually, computing key timescales analytically. Such BHs are captured by the disk via gas drag (timescale computed as in Wang et al. Reference Wang, Zhu and Lin2024). Once embedded, they undergo migration, which transports them towards locations of higher BH numerical density and favors binary formation (Bellovary et al. Reference Bellovary, Mac Low and McKernan2016). We model this analytically following McKernan et al. (Reference McKernan, Ford and Lyra2012). If the resulting timescale is shorter than the AGN lifetime
$\tau$
, we consider BBH formation and extract the mass
$m_2$
of the secondary BHs by physically-motivated distributions (Vaccaro et al. Reference Vaccaro, Mapelli and Périgois2024).
We show how the process of binary formation is sensitive to the component masses and the local disk properties. Migration traps, defined as locations where torques cancel and migration transitions from outwards to inwards, are favorable locations for BBH pair-up, but other locations in the disk are also viable sites for BBH assembly. The proportion between trap and non-trap mergers depends on the SMBH mass
$M_\bullet$
and the gas viscosity
$\alpha$
(Vaccaro et al. Reference Vaccaro, Seif and Mapelli2025).
After formation, BBHs evolve under three main processes: (i) Gas hardening, which exchanges angular momentum with the surrounding gas and is modeled following Ishibashi & Gröbner (Reference Ishibashi and Gröbner2024). In thin-disk regions (
$h \lesssim 0.04$
–
$0.16$
depending on the binary accretion rate) this shrinks the binary and excites its eccentricity; in thicker regions, it tends to widen it. (ii) Three-body encounters with single BHs, modeled by interfacing our code with tsunami (Trani et al. Reference Trani, Quaini and Colpi2024), can lead to exchanges, orbital tilts, and changes in semi-major axis and eccentricity (Samsing et al. Reference Samsing, Bartos and D’Orazio2022; Trani & Spera Reference Trani and Spera2023). These interactions generally shrink BBHs and increase eccentricity. (iii) Gravitational-wave emission, computed via the Peters (Reference Peters1964) formalism, drives the final inspiral and merger.
The remnant mass
$m_\mathrm{rem}$
and spin
$\chi_\mathrm{rem}$
are computed from numerical-relativity fits (Jiménez-Forteza et al. Reference Jiménez-Forteza, Keitel and Husa2017). Remnants may undergo further capture, migration, and binary formation, enabling hierarchical mergers, if three conditions are met: the recoil velocity
$v_\mathrm{kick}$
is below the local escape velocity
$v_\mathrm{esc}$
(met in
$99.4\%$
of cases); the cumulative accreted mass
$M_\mathrm{acc}=m_1+\sum m_2^{(N)}$
does not exceed the total captured BH mass in the disk (
$72.9\%$
); and the remnant’s formation timescale is shorter than
$\tau$
(
$41.3\%$
). When satisfied, this produces a chain of progressively more massive BBH mergers (Vaccaro et al. Reference Vaccaro, Mapelli and Périgois2024).
2. Results and Predictions
2.1. Population of hierarchical mergers in AGNs
As shown in Figure 1, our simulations show that AGN disks can efficiently produce a distinct population of hierarchical BBH mergers. This population extends the primary mass distribution well beyond
$100M_\odot$
, with a significant fraction reaching into and above the pair-instability mass gap (Vaccaro et al. Reference Vaccaro, Mapelli and Périgois2024). Mass ratios are often low (
$q \lesssim 0.1$
), reflecting the frequent pairing of massive merger remnants with lighter companions, and effective spins span extreme values (
$\chi_{\rm eff} \simeq \pm 0.8$
) due to repeated mergers and dynamical re-orientations. Spin and orbital properties are further diversified by three-body interactions, which can tilt spins, increase eccentricities, and in some cases exchange binary members.
Distributions of BBH merger properties from our AGN channel models. (a) Primary mass
$m_1$
, (b) mass ratio
$q=m_2/m_1$
, (c) effective spin
$\chi_\mathrm{eff}$
, and (d) precession spin
$\chi_\mathrm{p}$
. Solid blue and dashed orange lines show hierarchical mergers including all orbital orientations (“agnostic”) and restricted to prograde orbits at BBH formation (“prograde”), respectively. Light blue histograms show 1g mergers from our simulations. Grey curves and shaded regions indicate the median and 90% credible intervals from GWTC-3 (Abbott et al. Reference Abbott, Abbott and Acernese2023).

Figure 1 Long description
Panel A: A histogram shows the distribution of primary mass (m1) in solar masses (M☉). The x-axis ranges from 10^1 to 10^3 M☉, and the y-axis represents the probability density function (pdf). Solid blue and dashed orange lines represent hierarchical mergers with different orbital orientations, while light blue histograms show 1g mergers. Panel B: A histogram displays the distribution of mass ratio (q = m2/m1). The x-axis ranges from 10^-2 to 10^0, and the y-axis represents the pdf. Solid blue and dashed orange lines indicate hierarchical mergers, and light blue histograms show 1g mergers. Panel C: A histogram illustrates the distribution of effective spin (χeff). The x-axis ranges from -1.0 to 1.0, and the y-axis represents the pdf. Solid blue and dashed orange lines represent hierarchical mergers, while light blue histograms show 1g mergers. Panel D: A histogram shows the distribution of precession spin (χp). The x-axis ranges from 0.0 to 1.0, and the y-axis represents the pdf. Solid blue and dashed orange lines indicate hierarchical mergers, and light blue histograms show 1g mergers. Grey curves and shaded regions indicate the median and 90% credible intervals from GWTC-3.
2.2. Constraining the origin of observed BBH mergers
We assess the possible contribution of the AGN channel to the current GW detections by performing a hierarchical Bayesian mixture analysis including multiple formation channels. We use fastcluster predictions and evaluate transient event detectability for the LVK network as in Périgois (Reference Périgois2025). The mixing fraction
$f_\mathrm{AGN}$
is then inferred by comparing the model distributions to the posterior samples of BBH events in Abbott et al. (Reference Abbott, Abbott and Acernese2023). The inferred
$f_\mathrm{AGN}$
in our analysis has a median value of
$\sim 8\%$
with a 90% credible interval spanning from
$\lt1\%$
up to
$\sim 20\%$
(Vaccaro et al. Reference Vaccaro, Mapelli and Périgois2024). It is lower than for other dynamical channels (e.g. star clusters), but remains consistent with a non-negligible contribution of AGNs to the observed BBH population. A comparable upper limit (
$f_\mathrm{AGN} \lesssim 17\%$
for bright AGNs) has been independently obtained via a spatial correlation analysis between GW events and AGN catalogs by Veronesi et al. (Reference Veronesi, Rossi and van Velzen2023). The large uncertainties in our inference reflect both the limited number of detections and degeneracies between dynamical channels, underscoring that more stringent constraints will require the larger and more precise samples expected from third-generation detectors.
2.3. Tracing SMBH growth across cosmic time
The efficiency of hierarchical mergers in AGN disks depends sensitively on the properties of the host galaxy, particularly on the mass of its SMBH and the length of its cycles of active accretion. As a consequence, the cosmic evolution of the SMBH mass function directly impacts the redshift dependence of the AGN-channel BBH merger rate density. We explore two scenarios for SMBH growth, based on Trinca et al. (Reference Trinca, Schneider and Valiante2022): Eddington-limited accretion and merger-driven, episodically super-Eddington episodes. By combining the redshift-dependent SMBH mass function from these models with the cosmic evolution of the AGN number density
$n_\mathrm{AGN}(z)$
derived from cosmological simulations (Nelson et al. Reference Nelson, Springel and Pillepich2019), we predict the BBH merger rate density
$\mathcal{R}(z)$
for the AGN channel. In the merger-driven scenario, the early build-up of high-mass SMBHs leads to an excess of AGN-assisted mergers at
$z \gtrsim 3$
compared to the Eddington-limited case, whereas at lower redshift the two models behave in a similar way. Third-generation detectors such as the Einstein Telescope and Cosmic Explorer will extend the reach of BBH observations to high redshift, enabling direct measurements of the AGN-channel merger rate density in the epoch of SMBH assembly. A detection (or absence) of a high-z excess in
$\mathcal{R}(z)$
will thus provide a novel observational probe of SMBH growth, potentially distinguishing between slow Eddington-limited accretion and rapid, merger-driven growth across cosmic time.
Acknowledgments
The author thanks her collaborators M. Mapelli, A.A. Trani, Y. Seif, D. Wylezalek, A. Trinca, S. Torniamenti, and M. Dall’Amico. She gratefully acknowledges financial support from the MODEST-25 organizers. She also acknowledges support from the European Research Council for the ERC Consolidator grant DEMOBLACK, under contract no. 770017; from the German Excellence Strategy via the Heidelberg Cluster of Excellence (EXC 2181 - 390900948) STRUCTURES; from the state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through grants INST 35/1597-1 FUGG and INST 35/1503-1 FUGG. The public version of fastcluster is available via gitlab following this link.
m1
q=m2/m1
χeff
χp