© The Authors 2026

1. Introduction

The direct detection of gravitational waves (GWs) by the LIGO and Virgo collaborations has revolutionized our ability to study compact object populations in the Universe, most notably binary black holes (BBHs, Abbott et al. 2023a,b). As the catalog of GW events continues to grow, a pressing challenge is to identify and characterize the astrophysical environments in which such binaries form and evolve.

One particularly promising environment for BBH formation is the accretion disk of an active galactic nucleus (AGN), where a gas-rich and high-density environment provides ideal conditions for dynamical interactions among black holes (BHs, e.g., McKernan et al. 2012, 2025; Bartos et al. 2017; Stone et al. 2017; Yang et al. 2019a; Secunda et al. 2020; Tagawa et al. 2020a,b, 2023; Samsing et al. 2022; Vaccaro et al. 2024; Cook et al. 2025; Delfavero et al. 2025; Dittmann et al. 2025; Xue et al. 2025). Indeed, stellar-mass BHs can become embedded in the gaseous disk surrounding a supermassive BH (SMBH, Wang et al. 2024; Rowan et al. 2025a; Whitehead et al. 2025) and thus interact gravitationally with the disk, experiencing migration torques analogous to those governing planetary migration in protoplanetary disks (Paardekooper et al. 2010; McKernan et al. 2012; Bellovary et al. 2016; Masset 2017; Grishin et al. 2024).

As BHs migrate radially within the disk, the probability of close encounters and subsequent pair-up into binaries is enhanced (DeLaurentiis et al. 2023; Li et al. 2023; Rozner et al. 2023; Rowan et al. 2023, 2024; Whitehead et al. 2024b,a). Furthermore, the dissipation of orbital energy through gas drag and GW emission can efficiently harden such binaries, leading to mergers within the AGN’s lifetime (Peters 1964; Tiede et al. 2020; Ishibashi & Gröbner 2020, 2024; Rozner & Perets 2024; Rowan et al. 2025b).

BBH pair-up in AGN disks can occur either within migration traps – special radial locations where the total torque acting on a BH vanishes, halting its migration and promoting accumulation of several BHs – or throughout other regions of the disk, collectively referred to as “the bulk” (Wang et al. 2021). When two BHs migrate at different rates through the AGN disk, a process we refer to as differential migration, their radial paths may intersect. This increases the likelihood of close encounters that, under the influence of gas drag and GW emission, result in BBH formation. McKernan et al. (2020) find that more than 50% of BBH mergers originate in the bulk, but hierarchical mergers are efficient primarily within migration traps. Furthermore, Tagawa et al. (2020b) demonstrate that BBHs forming in the bulk through gas dynamical capture tend to migrate toward the migration trap as they harden.

Existing semi-analytical approaches focusing on the hierarchical merger scenario, including the FASTCLUSTER population synthesis tool developed by the authors (Mapelli et al. 2021, 2022; Vaccaro et al. 2024; Torniamenti et al. 2024), assign pair-ups to precomputed trap locations. While computationally convenient, this approach limits the physical realism of model outputs, especially for predicting the spatial distributions and formation times of merging binaries.

For this study, instead, we explicitly integrated BH migration trajectories within AGN disks. We used the Sirko & Goodman (2003) model to describe the disk structure and employed the publicly available pAGN module (Gangardt et al. 2024) to compute self-consistent radial profiles of gas surface density, temperature, and pressure. We evaluated migration torques using both a classical Type I prescription (Paardekooper et al. 2010; Bellovary et al. 2016) and an updated formulation including thermal effects (Masset 2017; Grishin et al. 2024). We accounted for transitions to Type II migration under gap-opening conditions (Lin & Papaloizou 1993; Bryden et al. 1999; McKernan et al. 2014; Vaccaro et al. 2024; Gilbaum et al. 2025).

We simulated millions of independent BH pairs, drawn from a mass distribution generated with population-synthesis calculations (Iorio et al. 2023), and tracked their one-dimensional migration under gravitational torques. We accounted for differential migration rates between BHs of varying masses, and identified binary formation using a combination of Hill-stability and orbit-crossing criteria. We further incorporated hierarchical evolution by seeding higher-generation (Ng) BHs based on the merger outcomes of first-generation (1g) binaries, including recoil kicks (Maggiore 2018).

Our main objective is to quantify the spatial distribution of BBH pair-ups across a broad parameter space. This includes variations in SMBH mass, disk viscosity, torque prescription, and BH mass. Our results provide realistic prescriptions for pair-up radii and timescales, which can be incorporated into semi-analytic population synthesis codes to improve merger rate predictions and GW data interpretation.

2. Methods

2.1. Disk structure

We modeled the physical properties of the accretion disk using the widely accredited steady-state analytic Sirko & Goodman (2003) disk model (hereafter SG). We solved the corresponding one-dimensional equations self-consistently via the pAGN module for the Python programming language (Gangardt et al. 2024) for given input parameters such as the mass of the central SMBH, M_•, and the viscosity coefficient, α. We fixed the mass accretion rate to a tenth of the Eddington rate, $\dot{M_{•}}=0.1{\dot{M_{•}}}_{,\mathrm{Edd}}$. We obtained radial profiles 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), the temperature, T, and the sound speed, c_s.

Self-consistent solutions are crucial, as approximate treatments can lead to large deviations in these profiles, which in turn dramatically affect predictions for BBH formation in AGN disks, as discussed in Appendix A.

2.2. Migration torques

We explicitly computed the torques felt by an embedded object of mass m at a radial distance R from a SMBH of mass M_•. We included the effects of Type I migration caused by the gravitational perturbation of embedded BHs in the disk (Paardekooper et al. 2010; Lyra et al. 2010; Jiménez & Masset 2017), accounting also for thermal migration, caused by the thermal response of the disk to the cold, overdense accretion disks surrounding embedded BHs (Masset 2017). The overall torque is

$$\begin{array}{c}\hfill \mathrm{\Gamma }={\mathrm{\Gamma }}_{\mathrm{type}\mathrm{I}}+{\mathrm{\Gamma }}_{\mathrm{th}}.\end{array}$$(1)

We investigated the difference between a classical Type I migration prescription, inspired by the work by Bellovary et al. (2016, hereafter B16), and one that also accounts for vertical and thermal effects, similar to Grishin et al. (2024, hereafter G24). The classical Type I migration torques were computed according to B16 as

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{type}\mathrm{I}}=\frac{{\mathrm{\Gamma }}_{\mathrm{ad}}{\mathrm{\Theta }}^2+{\mathrm{\Gamma }}_{\mathrm{iso}}}{{(\mathrm{\Theta }+1)}^2},\end{array}$$(2)

where the parameter Θ was computed according to Lyra et al. (2010) as¹

$$\begin{array}{c}\hfill \mathrm{\Theta }=\frac{c_{\mathrm{V}}{\mathrm{\Sigma }}_{\mathrm{g}}\mathrm{\Omega }{\tau }_{\mathrm{eff}}}{12\pi \sigma T^3},\end{array}$$(3)

and the isothermal Γ_iso and adiabatic Γ_ad torques are functions of the adiabatic index γ (assumed to be equal to 5/3 throughout this work) and the gradients of local surface density Σ_g, and temperature T

$$\begin{array}{c}\hfill a=-\frac{\partial log{\mathrm{\Sigma }}_{\mathrm{g}}}{\partial logR},b=-\frac{\partial logT}{\partial logR},c=b-(\gamma -1)a,\end{array}$$(4)

as

$$\begin{array}{cc}\hfill {\mathrm{\Gamma }}_{\mathrm{iso}}/{\mathrm{\Gamma }}_0& =-0.85-a-0.9b\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill \gamma {\mathrm{\Gamma }}_{\mathrm{ad}}/{\mathrm{\Gamma }}_0& =-0.85-a-1.7b+7.9c/\gamma .\hfill \end{array}$$(6)

Jiménez & Masset (2017) incorporated three-dimensional effects in their Type I migration torque,

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{type}\mathrm{I}}={\mathrm{\Gamma }}_{\mathrm{L}}+\left(\frac{0.46-0.96a+1.8b}{\gamma }\right){\mathrm{\Gamma }}_0,\end{array}$$(7)

and also accounted for the Lindblad torque Γ_L as

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{L}}/{\mathrm{\Gamma }}_0=(-2.34+0.1a-1.5b)f(\chi /{\chi }_c),\end{array}$$(8)

where

$$\begin{array}{c}\hfill f(x)=\frac{{(x/2)}^{1/2}+1/\gamma }{{(x/2)}^{1/2}+1},\chi =\frac{16\gamma (\gamma -1){\sigma }_{\mathrm{SB}}{T}^4}{3\kappa {\rho }_0^2{H}^2{\mathrm{\Omega }}^2},{\chi }_c=H^2\mathrm{\Omega }.\end{array}$$(9)

We took the electron-scattering opacity as κ = 0.34 cm²/g as in Gangardt et al. (2024).

The thermal torques were computed as (Masset 2017; Velasco Romero & Masset 2020)

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{th}}=1.61\frac{\gamma -1}{\gamma }\frac{x_c}{\lambda }(\frac{L}{L_c}\frac{4m_c}{m+4m_c}-\frac{2m_c}{m+2m_c})\frac{{\mathrm{\Gamma }}_0}{h}(1-e^{-\lambda {\tau }_V/h}),\end{array}$$(10)

where m_c = c_sχ/G is a critical mass, τ_V is the vertical optical depth, x_c is the corotational radius of the embedded body and λ is the typical size of overdense regions surrounding embedded BHs, evaluated as in G24:

$$\begin{array}{c}\hfill x_c=-\frac{H^2}{3\gamma R}\frac{\partial logP}{\partial logR},\lambda =\sqrt{\frac{2\chi }{3\gamma \mathrm{\Omega }}}·\end{array}$$(11)

We assumed that the embedded BH emits at Eddington luminosity, L = L_Edd = 4πGmc/κ, and took L_c = 4πGmρχ/γ as the critical luminosity (Masset 2017).

All migration torques are normalized by

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_0={\left(\frac{m}{M_{•}h}\right)}^2{\mathrm{\Sigma }}_{\mathrm{g}}{R}^4{\mathrm{\Omega }}^2.\end{array}$$(12)

We identified migration traps as locations in the disk where the total Type I migration torque Γ crosses zero, changing sign from positive to negative. These correspond to regions where BHs tend to accumulate. In contrast, locations where Γ changes sign from negative to positive are referred to as anti-traps, as they tend to be depleted of BHs.

We also accounted for the possibility of switching to Type II migration. Namely, we enforced a joint condition on the gas viscosity, migrator mass and local disk thinness, (Crida et al. 2006)

$$\begin{array}{c}\hfill K=\frac{m^2}{\alpha h^5{M}_{•}^2}>11.\end{array}$$(13)

If this condition is met, the migrator opens a gap in the disk and can no longer be subject to the migration torque Γ. Instead, conservatively assuming that no gas can cross the gap, its migration follows the viscous evolution of the disk’s gas, moving inwards on a timescale (McKernan et al. 2012)

$$\begin{array}{c}\hfill t_{\mathrm{type}\mathrm{II}}=t_{\mathrm{visc}}={(\alpha h^2\mathrm{\Omega })}^{-1}.\end{array}$$(14)

Our prescription for the disk profile does not allow us to self-consistently open a gap (i.e., set Σ_gas(R)≃0 locally), hence in case of Type II migration we artificially constructed the corresponding torque as the rate of change in angular momentum, L = mR²Ω,

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{type}\mathrm{II}}={\frac{\mathrm{d}L}{\mathrm{d}t}|}_{\mathrm{type}\mathrm{II}}\simeq \frac{L}{t_{\mathrm{type}\mathrm{II}}}=\alpha m{H}^2{\mathrm{\Omega }}^2.\end{array}$$(15)

We discuss the validity and impact of this assumption in Appendix B.

2.3. Integration scheme and pair-up criterion

We began each simulation by selecting two stellar-mass BHs that were assumed to already be embedded in the AGN disk. For stellar-origin BHs, henceforth referred to as first generation (1g) BHs, masses m are sampled independently from a distribution obtained via population-synthesis modeling SEVN for solar metallicity (Iorio et al. 2023). The initial radial location R_init, is drawn from a uniform distribution in logarithmic space. Specifically, the logarithm of the initial radius is uniformly sampled between the logarithm of the disk’s inner radius, assumed to be at the innermost stable circular orbit, R_min = 6 R_g in this context, and an outer radius R_max = 0.1 pc (M_•/10⁶ M_⊙)^1/2, beyond which the disk’s self-gravity becomes important (Goodman 2003; Yang et al. 2019b). This choice reflects the expectation that BHs can form or become embedded at various locations throughout the AGN disk, and provides a neutral baseline that does not bias our results toward any particular radial scale. However, the stellar-origin BHs we expect to find in an AGN disk come from the surrounding NSC, and the probability that they are captured by the disk depends on the properties of both the AGN disk and the NSC. As a result, embedded BHs are expected to accumulate preferentially in particular radial ranges, enhancing the probability of binary formation at those locations relative to our log-uniform assumption.

The subsequent orbital evolution is treated in one dimension, considering only the radial location as a function of time. The equation of motion is derived from angular momentum conservation in a Keplerian potential and takes the form:

$$\begin{array}{c}\hfill \dot{R}=\frac{2\mathrm{\Gamma }(R)}{m}\sqrt{\frac{R}{G{M}_{•}}},\end{array}$$(16)

where Γ(R) is obtained as in Eq. (1) or (15) depending on whether the conditions for Type II migration (Eq. 13) are met. The torque profile is precomputed from the AGN disk model and interpolated at each timestep to reflect the actual mass and location of the migrating body.

We solved this system of coupled differential equations using the second-order Runge–Kutta integration method with embedded third-order error estimation and an adaptive timestep. The integration spans the full duration of the AGN disk’s activity, which we fix at 2 Myr, and is initialized with a timestep of 10³ yr. To capture statistical trends in binary formation, this procedure was repeated for a large number of random initial conditions. We simulated 10⁶ independent realizations, each with different initial masses and orbital radii, resulting in a comprehensive catalog of pair-up events across a wide range of physical parameters.

A central goal of the simulation is to determine whether and when a given pair of migrating single BHs becomes a bound binary. To that end, we randomly select two BHs in our catalog and we label them as primary (mass m₁, radial location R₁) and secondary (m₂, R₂), with the convention that m₁ ≥ m₂. We applied two complementary criteria for binary formation in post-processing, after completing the orbital integration of each single BH.

The first and most stringent condition is based on Hill stability (Hill 1878). We computed the instantaneous separation between the two BHs, d = |R₁ − R₂|, at each integration step and compare it to their Hill radii. A binary is considered to have formed if the separation falls below the larger of the two Hill radii, hence the two BHs are close enough for their mutual gravitational attraction to overcome tidal forces from the central SMBH:

$$\begin{array}{c}\hfill d\le \underset{i=1,2}{max}{R}_{\mathrm{H}}(m_i,R_i),R_{\mathrm{H}}(m,R)=R{\left(\frac{m}{3(M_{•}+m)}\right)}^{1/3}.\end{array}$$(17)

To prevent the integrator from skipping over close encounters due to large timesteps or rapid migration, we also monitored the sign of the radial separation E_cross = R₁ − R₂. A sign change in E_cross indicates that the BHs have crossed each other’s orbits and may have come into sufficiently close proximity to allow for capture. If either the explicit Hill condition or this indirect orbit-crossing condition is met, a BBH is considered to have formed.

In either case, we recorded the BBH formation time and its pair-up radius, defined as the mass-weighted average position of the two BHs at the moment of binary formation:

$$\begin{array}{c}\hfill R_{\mathrm{pair}}=\frac{m_1{R}_1+m_2{R}_2}{m_1+m_2}·\end{array}$$(18)

This dual-criterion strategy balances accuracy and robustness, ensuring that bound systems are identified even when the BHs migrate rapidly or interact near the resolution limit of the integrator. However, it neglects the exact binary binding mechanism under the influence of gas torques and disk turbulence, investigated in hydrodynamical simulations (e.g. Whitehead et al. 2024b). The resulting catalog of binary formation events provides the foundation for our statistical analysis of BBH assembly in AGN disks.

2.4. Nth-generation BHs

For Nth-generation (Ng) BHs, which are the remnants of earlier mergers in the AGN disk, we estimated their orbital starting conditions using the outputs of 1g simulations. Hence, in the Ng simulations, the initial position R_init of the primary BH is no longer drawn randomly across the disk. Instead, it is sampled directly from the set of R_pair values obtained in the 1g simulation, accounting for the natal velocity kick v_kick that is imparted during the merger event. The recoil velocity v_kick is drawn from a precomputed FASTCLUSTER catalog (Vaccaro et al. 2024) that already incorporates the dependence on progenitor mass ratio and spin configuration as in Maggiore (2018). The kick direction is randomized assuming isotropic ejection, and contributions from enclosed gas and NSC members are neglected given their small mass relative to the SMBH. If the kick magnitude exceeds the local escape velocity, the remnant is assumed to leave the system. Otherwise, the object is recaptured and is then able to migrate in the disk following the same torque prescriptions as 1g BHs. The shifted radius $R_{\mathrm{init}}^{\mathrm{new}}$ is computed by means of simple orbital transfer calculations (Hohmann 1960) as

$$\begin{array}{c}\hfill R_{\mathrm{init}}^{\mathrm{new}}=\frac{G{M}_{•}}{v_{\mathrm{tot}}^2}\left[\frac{R_{\mathrm{init}}{v}_{\mathrm{tot}}^2/G{M}_{•}}{2-R_{\mathrm{init}}{v}_{\mathrm{tot}}^2/G{M}_{•}}\right],\end{array}$$(19)

where v_tot is the vector sum of the circular orbital velocity $v_{\mathrm{Kepl}}=\sqrt{G{M}_{•}/R}$ and the random kick v_kick.

In doing so, we implicitly assumed that 1g binaries will eventually merge at a radial location not significantly different from where they originally paired up.

We only allowed Ng BHs to pair up with first-generation companions (forming Ng − 1g BBHs). This choice is motivated by the low probability of drawing a secondary of generation M, given the steeply decreasing p(M)∝2^{−(M − 1)} (Zevin & Holz 2022). Accordingly, the secondary BH in each pair has mass and initial radial position sampled in the same way as in the 1g simulations.

This methodology ensures that Ng BHs are seeded in dynamically consistent locations, and allows us to model the spatial memory and redistribution of hierarchical mergers across the AGN disk. Other than the initial radius, the rest of the computation proceeds in an identical way to what previously described for 1g BHs.

2.5. Description of runs

We performed a suite of Monte Carlo simulations, running N = 10⁵ realizations for each combination of model parameters. We explore how BBH formation in AGN disks depends on the mass of the central SMBH, the disk’s viscosity, and the underlying migration prescription.

We vary the SMBH mass over the range log M_•/M_⊙ = 5.0 − 9.0, encompassing typical values expected across cosmic time (e.g. Greene & Ho 2007; Trinca et al. 2022).

The disk viscosity for a SG disk, and more generally for so-called α-disk models, is parametrized by coefficient α ∈ [0,  1]. We adopt α = 0.1 and α = 0.4 to bracket the range suggested by observations of fully ionized accretion disks, which typically yield α ∼ 0.1 − 0.4 (King et al. 2005; Martin et al. 2019). We also ran some models with α = 0.01 in order to align with conventional assumptions in theoretical models that often presume lower turbulence efficiency (e.g. standard disk models such as SG) and magneto-hydrodynamical simulations (e.g. Hawley et al. 2011). We are limited by the assumption of a spatially constant α across the disk radius, consistent with the approximation in the SG model, despite more recent simulations indicating potential radial variation in viscosity (e.g. Penna et al. 2012).

To investigate the impact of different assumptions on the migration mechanism, we distinguished between a B16 prescription, where we set Γ_th = 0 in Eq. (1) and computed Γ_type I as in Eq. (2) (as in Bellovary et al. 2016), and a G24 prescription where we set Γ_th as in Eq. (10) and Γ_type I as in Eq. (7) (as in Grishin et al. 2024).

3. Results

3.1. Migration torques and pair-up radii

We begin by illustrating the torque structure that governs BH migration in SG disks. Figure 1 shows representative profiles of the total torque acting on a BH as a function of radial distance.

Fig. 1. Examples of torque profiles as a function of radial distance (in Schwarzschild radii) for BHs embedded in a SG disk. Each plot corresponds to a different combination of torque prescription (top: G24; bottom: B16), disk viscosity parameter α (left: α = 0.01; right: α = 0.1) and SMBH mass (log M•/M⊙ = 5.0, 6.3, 7.7, 9.0, represented in different colors). Solid and dashed lines indicate positive and negative net torques, respectively, while shaded regions denote uncertainty in the torque determined by computing it over a range of migrator masses between 5 M⊙ and 50 M⊙. Vertical dotted lines mark the locations of migration traps.

The torque intensity and sign vary depending on the local properties of the disk such as temperature and density gradients as in Eqs. (2)–(9), which are themselves functions of α, R, and M_•. The mass of the migrator m only affects the intensity of the torque via the normalization factor in Eq. (12). Hence, migration trap locations are disk-specific and independent on the characteristics of the migrating BHs. Multiple torque sign reversals from Γ > 0 to Γ < 0 are generally present, corresponding to migration traps, whose number and position are of course also sensitive to M_• and α.

To connect this torque structure with the outcome of binary formation, we show in Figure 2 the full distribution of BBH pair-up locations across our catalog of simulations. Each panel displays a two-dimensional histogram in the radius–SMBH mass plane, indicating the logarithmic number of pair-ups in each bin. We also indicate the locations of migration traps, as extracted from the torque profiles. Most BHs that successfully form binaries do so while migrating under Type I migration, with fewer than 8% forming under Type II conditions.

Fig. 2. Two-dimensional histograms showing the distribution of BBH pair-up locations in the radius–mass plane (R vs. M•), for different combinations of torque prescription (top: G24; bottom: B16) and disk viscosity α (left: 0.01; right: 0.1) in a SG disk. The color map indicates the number of pair-ups in each bin Npair. Red crosses mark the locations of migration traps.

We find that the radial distribution of pair-ups is broadly consistent with the location of migration traps, supporting the idea that traps are favorable sites for binary formation. However, the distribution shows a significant spread around traps as well as some clear offsets, with efficient pair-ups also occurring at locations not identified as traps, while some traps contain no BBH pair-ups at all. These departures can be attributed to differential migration: mismatches in migration speed may cause pair of BHs to pair up upstream of the trap radius.

In some cases, such as at R ≃ 10³ R_S for M_• < 10⁶ M_⊙ in the B16 prescription, there is a clear overdensity of pair-ups at locations which are not marked as traps, caused by sharp changes in the slope of Γ. We refer to these locations as traffic jams.

In systems with smaller SMBHs (log M_•/M_⊙ ≲ 6), pair-ups are possible at almost all radial locations in the disk, whereas at higher SMBH masses (log M_•/M_⊙ ≳ 7), the distribution of pair-up radii becomes significantly narrower and more concentrated toward the inner disk regions (10² ≲ R/R_S ≲ 10⁴). This narrowing reflects the increased gravitational influence of more massive SMBHs which dominates over local disk perturbations, leading to a more stable torque structure. In particular, the radial torque Γ at large radii is less intense and uniformly negative at high M_•, reducing differential migration between BHs and suppressing their pairing in the outer disk region.

We observe relevant differences in the trap location based on the torque prescription, as shown in Figure 2. Using the full G24 migration torque gives rise to outer migration traps which are absent in the classical B16 case for SMBH masses log M_•/M_⊙ ∈ [5.6, ∼8.0] ([6.7, 9.0]) for α = 0.01 (0.1).

3.2. The relevance of migration traps

To quantify the role of migration traps in BBH formation, we track the radial locations at which BHs pair up and compare them to the positions of predefined migration traps in the disk. For each simulation, we identify a set of trap locations {R_trap} based on the local torque profile, as described in Section 2.2. A pair-up event is considered to occur “within a trap” if the radial distance between the pair-up location R_pair and any trap satisfies the following logarithmic tolerance condition:

$$\begin{array}{c}\hfill |log{R}_{\mathrm{pair}}-log{R}_{\mathrm{trap}}|<\delta ,\end{array}$$(20)

where δ is a fixed tolerance. In this work, we adopt δ = 0.05.

Based on this definition, we compute the fraction f_trap as

$$\begin{array}{c}\hfill f_{\mathrm{trap}}=\frac{N_{\mathrm{in}\mathrm{trap}}}{N_{\mathrm{pair}}},\end{array}$$(21)

where N_in trap is the number of pair-up events satisfying the condition in Eq. (20), and N_pair is the total number of pair-up events recorded in the simulation.

Figure 3 shows the resulting fraction f_trap, plotted as a function of the central SMBH mass. Overall, the majority of pair-ups occur within migration traps across most of the parameter space explored. For SMBH masses below a certain threshold $M_{•}^{\mathrm{thr}}$, f_trap has values of 0.7 to 0.9 for 1g BHs and it approaches unity for Ng BHs. As the SMBH mass increases, f_trap presents a pronounced transition from values close to unity to nearly zero, indicating that pair-up events are no longer connected with trap locations at high SMBH masses. This transition typically takes the form of an abrupt drop from f_trap ≃ 1 to f_trap ≃ 0, with specific value of $M_{•}^{\mathrm{thr}}$ depending on the disk and torque model adopted.

Fig. 3. Fraction of binary pair-ups occurring within migration traps as a function of the SMBH mass. The vertical axis shows the occurrence ftrap of pair-up events that take place in proximity to predefined migration trap locations, while the horizontal axis shows the logarithm of the central SMBH mass. Different marker fillings indicate the BH generation (filled: 1g; hollow: Ng), while color encodes the viscosity parameter α (orange: 0.01; teal: 0.1; navy blue: 0.4). The two torque formalisms are distinguished by marker shape and line style (squares and solid lines: G24; circles and dashed lines: B16). Trap proximity is defined using a logarithmic tolerance on radial separation (see Section 3.2 for details).

In the high-viscosity B16 models (α = 0.1, 0.4), the drop occurs early, at $M_{•}^{\mathrm{thr}}\simeq {10}^{7.4}{\mathrm{M}}_{\odot }$, whereas in the low-viscosity case (α = 0.01), $M_{•}^{\mathrm{thr}}\simeq {10}^{8.5}{\mathrm{M}}_{\odot }$. For the G24 model, the transition to low f_trap consistently occurs at large SMBH masses; specifically, at $M_{•}^{\mathrm{thr}}\simeq {10}^{8.8}{\mathrm{M}}_{\odot }$, > 10⁹ M_⊙, ≃10^8.5 M_⊙ for the cases with α = 0.01, 0.1, 0.4. Furthermore, f_trap curves present some outliers due to peculiar disk properties for specific combinations of M_• and α (see Figure 2 and additional |Γ| plots on Zenodo). For instance, in disks with log M_•/M_⊙ ∈ [5.0, 6.0] and α = 0.01, f_trap increases gradually from 0.2 to ∼0.8 for the B16 torque due to traffic jams: a sharp transition in the slope of Γ gives rise to many pair-up events outside of traps due to differential migration. Specifically, as visible in Figure 1, this happens at R ≲ 10² R_S for log M_•/M_⊙ = 5.0 − 5.1 and at R ∼ 10³ R_S for log M_•/M_⊙ < 6.0.

We show normalized histograms of binary pair-up radii R_pair, split by BH generation (1g vs Ng), in Figure 4 for selected SMBH masses and disk viscosity. Pair-up locations for Ng BHs are considerably more concentrated around peaks than for 1g BHs. This enhanced localization arises from the initial conditions: in our model, Ng remnants resume migration at a radius set by Eq. (19). Because recoil speeds v_kick are typically small relative to the local Keplerian velocity v_Kepl (Maggiore 2018; Vaccaro et al. 2024), the kick induces only a minor radial displacement, so Ng remnants remain close to the pair-up radius of their progenitor binary. Since 1g binaries preferentially form near peaks in the R_pair distribution, the ensuing Ng BHs tend to start their migration near those same radii. As a result, R_pair for Ng objects clusters tightly around peaks, which naturally explains the systematically higher values of f_trap seen for Ng in Figure 3 since most peaks are coincident with traps. In contrast, 1g BHs are initialized with randomly distributed radial positions within the disk, and their subsequent migration is governed entirely by the local torque profile. While traps still influence the dynamics, the broader range of initial conditions results in a more diffuse distribution of pair-up radii, and therefore a lower f_trap.

Fig. 4. Histograms of binary pair-up radii Rpair (in units of Schwarzschild radii) for first-generation (1g) and Nth-generation (Ng) BHs, across a range of SMBH masses (log M•/M⊙ = 6, 7, 8 in the three columns, respectively) and viscosity parameters (upper row: α = 0.01; lower row: α = 0.1). Panels show the pair-up radii probability density functions (PDFs), normalized separately in each panel. The 1g distributions are shown as step histograms, while Ng distributions are overlaid either as filled or hatched histograms for visual distinction. Vertical dotted lines indicate the positions of migration traps identified from the torque profile. Color represents the migration model used (green: G24; purple: B16).

For low viscosity (α = 0.01) under the B16 prescription, we also find a pronounced pile-up of pair-ups at R ≃ 10³ R_S for M_• ≲ 10⁶ M_⊙ that does not coincide with migration trap locations. This is a direct consequence of differential migration: as shown in Figure 1, the |Γ(R)| profiles exhibit an abrupt change in slope at these radii, producing a traffic-jam accumulation of pair-ups.

These findings emphasize the importance of migration traps as anchoring points for BBH formation in AGN disks. Their effectiveness is sensitive to both disk properties (through the viscosity parameter α) and central SMBH mass, and is particularly crucial for Ng pair-ups and their resulting hierarchical mergers.

3.3. Pair-up efficiency

We define the pair-up efficiency η as the fraction of BHs in a given simulation that successfully form a binary by the end of the AGN disk lifetime. That is,

$$\begin{array}{c}\hfill \eta =\frac{N_{\mathrm{pair}}}{N_{\mathrm{tot}}},\end{array}$$(22)

where N_pair is the number of BHs that participated in a pair-up event, and N_tot is the total number of simulated BHs. Figure 5 shows the behavior of η as a function of the central SMBH mass, for both 1g and Ng BHs. As expected, the efficiency tends to decrease with increasing SMBH mass, due to the suppression of differential migration in the outer regions of the disk discussed in Section 3.1. Nonetheless, Ng BHs consistently exhibit higher pair-up efficiency than their 1g analogs. This is a consequence of their more favorable initial conditions as previously observed for f_trap: they typically start closer to trap or traffic-jam locations, which increases their likelihood of encountering a suitable companion before the end of the disk lifetime.

Fig. 5. Pair-up efficiency as a function of the central SMBH mass. The vertical axis shows the fraction of BHs that successfully form a binary during the simulation, while the horizontal axis shows the logarithm of the SMBH mass. As in Figure 3, different marker fillings indicate the BH generation (filled: 1g; hollow: Ng), color encodes the viscosity parameter α (orange: 0.01; teal: 0.1; navy blue: 0.4), and the two torque prescriptions are distinguished by marker shape and line style (squares and solid lines: G24; circles and dashed lines: B16).

We also observe a systematic, although modest, dependence of η on the value of α and on the assumed migration model. For low-viscosity disks (α = 0.01), pair-up efficiency remains higher across a broader range of SMBH masses, whereas higher values of α (0.1, 0.4) slightly reduce the likelihood of BBH formation. Interestingly, these differences are more pronounced for the Ng BHs than for 1g, where the impact of α is comparatively marginal. This may be due to the inherently stochastic nature of pair-up dynamics in the outer disk regions at high viscosity (α ≥ 0.1), relatively low SMBH masses (log M_•/M_⊙ ≲ 6.5), where multiple nearby trap locations (noticeable in Figures 1 and 2) can lead to complex migration patterns.

We also note marked jumps in η which correlate with the temporary disappearance of outer migration traps in the corresponding M_• range, which allows outer BHs to move toward intermediate-radius (∼10³ R_S) locations, characterized by sudden changes in the slope of Γ which facilitate BBH captures by differential migration. This happens most notably at log M_•/M_⊙ ∼ 5.6 − 5.8 for α = 0.01, and at log M_•/M_⊙ ≳ 6.7 for α = 0.1.

3.4. Impact of initial conditions

Figure 6 is an analysis of how the pair-up efficiency η and the trap fraction f_trap depend on the initial conditions of the migrating BH, specifically its initial radial position R_init and mass m.

Fig. 6. Pair-up efficiency (top panels) and fraction of pair-ups occurring within migration traps (bottom panels) as functions of initial radius Rinit (left) and migrator mass m (right). Color encodes the viscosity parameter α (orange: 0.01; teal: 0.1; navy blue: 0.4), and the two torque prescriptions are distinguished by marker shape and line style (squares and solid lines: G24; circles and dashed lines: B16).

The pair-up efficiency η shows a strong radial dependence, peaking at intermediate radii (R_init ∼ 10² − 10⁵ R_S). This trend reflects the underlying structure of the torque profile: migration is inefficient near the inner and outer edges of the disk, while the central regions allow for efficient dynamical pairing. In the innermost disk, pair-ups are suppressed by limited radial volume and weaker torques, whereas in the outer disk, BHs typically migrate too slowly to encounter a companion within the AGN lifetime. A secondary dependence on mass is also evident: heavier BHs (m ≳ 20 M_⊙) tend to achieve higher pair-up efficiencies due to their stronger interaction with the gas and faster migration rates.

The trap fraction f_trap follows a similar, but more pronounced, trend to the pair-up efficiency η in terms of BH mass m. In particular, stellar-mass BHs with m ≲ 20 M_⊙ are generally unable to reach migration traps, with average f_trap remaining below 0.3, whereas BHs with m ≳ 20 M_⊙ consistently exhibit trap fractions exceeding 0.6 (0.9) in the B16 (G24), indicating that higher-mass objects are not only more efficient migrators but also more likely to pair up in the vicinity of traps.

However, the dependence of f_trap on the initial radius R_init appears more complex than that of the overall efficiency η. While the peak values of both quantities occur broadly in the same intermediate radial range (∼10² − 10⁵ R_S), the trap fraction displays sharper and more localized peaks. This suggests that proximity to a trap at birth plays a significant role: BHs that happen to become embedded sufficiently close to a stable trap are more likely to remain in its vicinity and undergo in-trap pair-up. Consequently, f_trap is subject to a degree of stochasticity tied to the initial orbital configuration, which is less evident in the smoother behavior of η.

3.5. Migration timescales

In semi-analytical models of BH dynamics in AGN disks, the migration timescale is often estimated using the rate of angular momentum change via gaseous torques. A common approximation (e.g. McKernan et al. 2012; Vaccaro et al. 2024) is to compute the migration time as the ratio between the angular momentum L of the BH and the torque normalization Γ₀:

$$\begin{array}{c}\hfill t_{\mathrm{migr},0}=\frac{L}{{\mathrm{\Gamma }}_0}=\frac{M_{•}^2{h}^2}{m{\mathrm{\Sigma }}_{\mathrm{g}}{R}^2\mathrm{\Omega }}\end{array}$$(23)

where L = mR²Ω and Γ₀ is as in Eq. (12). This simplified expression is widely used in the literature as it captures the characteristic timescale over which a BH would drift through the disk under idealized conditions. In Figure 7, the resulting timescales computed for a 10 M_⊙ BH are shown as dashed lines for different SMBH masses, viscosity parameters, and torque prescriptions.

Fig. 7. Two-dimensional distributions of the migration time tmigr as a function of the initial radius Rinit in SG AGN disks, for different combinations of the viscosity parameter α (left: 0.01; right: 0.1) and torque prescriptions (top: G24; bottom: B16). Each panel shows contours enclosing 40% and 68% of the migration time distribution, based on our simulations. Results are shown for four representative SMBH masses: log(M•/M⊙) = 5.0, 6.3, 7.7, 9.0, indicated by color. Solid (dashed) lines represent the migration timescale computed using the torque Γ in Eq. (1) (Γ0 in Eq. 12). Vertical dotted lines mark the locations of migration traps.

However, when using the full torque Γ (Eq. 1) the corresponding migration timescale, t_migr = L/Γ, deviates from the Γ₀-based estimate, especially in regions near torque sign-reversal points, i.e., in the vicinity of migration traps or anti-traps, where Γ = 0. As a result, the solid lines in Figure 7 exhibit sharp features and even divergences, reflecting the breakdown of the timescale approximation in those zones.

To assess how well either of these approximations tracks the actual dynamical evolution of BHs in the disk, we compare them to the pair-up times, t_pair, obtained directly from our simulations. These pair-up times, shown via the 40% and 68% contour levels in Figure 7, represent the elapsed time between the moment at which the single primary BH becomes embedded in the disk and the subsequent formation of a BBH. Unlike t_migr, they are measured durations rather than derived timescales, and naturally include the effects of stochasticity, mass differences, and local disk structure. Several important features emerge from this comparison. There is substantial scatter in the pair-up times at fixed radius, due in part to the variation in the mass of the migrating BH. There is a sharp upper limit at t_migr = 2 Myr, corresponding to the assumed AGN disk lifetime adopted in our model. Moreover, some radial regions are systematically disfavored for binary formation, as also noticeable in Figure 6, due to either negative torques in the innermost regions of the disk or to reduced differential migration in the outer disk regions for large M_•, as discussed in Section 3.1. Finally, pair-up times in our simulations are typically shorter than the estimated migration timescale in regions where the Γ profile has a large slope, and longer where the profile is flatter. This is a consequence of differential migration: pair-ups are set by relative drift rather than absolute torque. Where the migration speed profile is steep, nearby objects converge quickly and pair up fast. Vice versa, where the profile is flat, velocity differences are only mass-dependent and consequently pair-up times lengthen.

4. Caveats

Our current approach still presents several limitations.

First, we do not model gas accretion onto individual BHs. Accretion is expected to lead to spin alignment and mass growth, which in turn influences migration rates and pair-up probabilities. We also neglect any vertical stratification in the AGN disk and restrict all BHs to move strictly in the midplane, thereby omitting any effect related to vertical motion. In addition, the SG model is not fully realistic, as it assumes an unspecified heating source and does not account for gas consumption even in gravitationally unstable regions.

Our treatment of Type II migration is deliberately conservative, as gap-opening objects are assumed to effectively prevent any inflow across the gap. More realistically, the transition to Type II is not a hard switch: the gap would not open fully and gas inflow would not be completely suppressed (Malik et al. 2015). As a result, the Type II regime would more closely resemble a modified Type I migration with reduced local density (as shown in Kanagawa et al. 2018; Fung et al. 2014). We explore this alternative model for Type II migration in Appendix B.

Moreover, our model does not include resonant interactions, such as mean-motion resonances, which could alter or halt migration (Epstein-Martin et al. 2025; Moncrieff et al. 2026). We also do not consider dynamical formation via three-body scatterings, which is another viable mechanism of BBH formation in AGN environments (Tagawa et al. 2020b). Finally, while our pair-up criterion identifies candidate binaries based on proximity, it does not explicitly model binary binding or survivability under the influence of gas torques and disk turbulence.

We will cure these approximations in a forthcoming study.

5. Summary

Semi-analytical models such as FASTCLUSTER (Mapelli et al. 2021, 2022; Vaccaro et al. 2024; Torniamenti et al. 2024) and MCFACTS (McKernan et al. 2025; Cook et al. 2025; Delfavero et al. 2025) have proven to be useful tools for modeling hierarchical mergers of BHs in AGN disks. Their computational efficiency allows for rapid exploration of the parameter space and facilitates comparisons with observational data. However, these models often rely on simplified assumptions regarding the locations at which BBHs form. A common assumption, made for example by Vaccaro et al. (2024), is that all pair-ups occur at so-called migration traps. While this assumption captures part of the physical picture, it overlooks the diversity of pair-up locations that can arise from realistic differential migration.

In this work, we have achieved a more realistic distribution of pair-up radii via explicit integration of migration torques computed for self-consistent SG AGN disk models modeled with the software pAGN (Gangardt et al. 2024). We identify migration traps as any location in the disk where the torque changes sign from positive to negative. Our key findings are as follows:

• The radial distribution of BBH pair-ups mainly follows the locations of migration traps, but with some spread and offsets due to differential migration. The overall pair-up radii profile is mostly influenced by the mass of the SMBH M_•, specifically: (i) For log M_•/M_⊙ ≲ 6.5, BBH pair-ups occur broadly across the entire disk; (ii) for log M_•/M_⊙ ≳ 6.5, pair-ups concentrate in the inner disk (10² ≲ R/R_S ≲ 10⁴).

• The efficiency of BBH formation η is mildly influenced by the disk parameters, decreasing for increasing M_• and α, and it it mostly influenced by the initial conditions of the migrating BHs. Indeed, BBHs are formed more efficiently for intermediate initial radii (R_init ∼ 10² − 10⁵ R_S) and large migrator mass (m ≳ 20 M_⊙).

• A large fraction of BBH pair-ups (f_trap ≳ 80% for first-generation BHs, f_trap ∼ 100% for N-th generation BHs) occur in migration traps for low-to-intermediate SMBH masses $M_{•}<M_{•}^{\mathrm{thr}}$. The threshold value is model dependent: $M_{•}^{\mathrm{thr}}\simeq {10}^{7.4}{\mathrm{M}}_{\odot }$ for the high-viscosity (α ≥ 0.1) in the B16 case, or $M_{•}^{\mathrm{thr}}\ge {10}^{8.5}{\mathrm{M}}_{\odot }$ for all other models. At larger SMBH masses, f_trap drops to lower values.

• For low SMBH masses (log M_•/M_⊙ < 6.0), low viscosity (α = 0.01) and B16 torque prescriptions, pair-up events pronouncedly pile up at R ≃ 10³ R_S as consequence of differential migration: sharp changes in the slope of |Γ(R)| create a traffic-jam accumulation even in the absence of a trap.

• Estimating migration timescales as the inverse of the rate of angular momentum change ($t_{\mathrm{migr}}\simeq L/\dot{L}=L/\mathrm{\Gamma }$) tends to overestimate pair-up times in outer disk regions. True pair-up times (from simulations) are shorter due to differential migration and movement into regions of stronger torque.

Our results demonstrate that, although BBH formation predominantly happens in migration traps for a large portion of the explored parameter space, formation in the bulk is non-negligible.

By integrating the migration trajectories of individual BHs and applying physically motivated pair-up criteria, we recover a distribution of pair-up radii that reflects the full complexity of orbital evolution in AGN disks. This highlights the limitations of using predetermined locations to represent the entire BBH formation process.

Incorporating realistic pair-up radius distributions into semi-analytic population models will improve the physical accuracy of merger rate predictions. It allows for a more sensible representation of where and when binaries are likely to form, and it captures the contribution of dynamically overtaking encounters which is entirely missed when assuming a static trap location.

Data availability

Supplementary material, data, and code used in this study are available via Zenodo at this link.

Acknowledgments

We thank the anonymous referee for their constructive and insightful comments that improved this manuscript. We thank Kees Dullemond, Dominika Wylezalek, Ralf Klessen, Evgeni Grishin, Nick Stone, Alessandro A. Trani, K.E. Saavik Ford, Barry McKernan, Shawn Ray, Thomas Rometsch, Zoltan Haiman, Shmuel Gilbaum, and Lumen Boco for useful discussions. MPV and MM acknowledge financial support from the European Research Council for the ERC Consolidator grant DEMOBLACK, under contract no. 770017. They also acknowledge financial support from the German Excellence Strategy via the Heidelberg Cluster of Excellence (EXC 2181 – 390900948) STRUCTURES. The authors acknowledge support by 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. We thank Daria Gangardt for developing the pAGN software, which is publicly available via github following this link. This research made use of NUMPY (Harris et al. 2020), SCIPY (Virtanen et al. 2020), PANDAS (The pandas development Team 2024). Plots were produced using MATPLOTLIB (Hunter 2007).

References

Appendix A: Improved disk modeling and comparison to previous work

A proper and realistic description of the properties of AGN disks is crucial for modeling compact object dynamics in such environments. In this work, we model AGN disks using the analytical tool pAGN, developed by Gangardt et al. (2024), to numerically solve Sirko & Goodman (2003, SG) equations of disk structure. A similar approach has recently been adopted in most semi-analytical models for the AGN channel (e.g. Vaccaro et al. 2024; Cook et al. 2025; Delfavero et al. 2025). This provides smooth profiles for key quantities for any combination of system parameters in a manner fully consistent with the physical description of SG. The resulting profiles for selected SMBH masses M_• are shown in Figure A.1 and Figure A.2 (hereafter A.1 and A.1) for high (α = 0.1) and low (α = 0.01) viscosity respectively.

Fig. A.1. Radial profiles of surface density Σg (top panel) and aspect ratio h = H/R (bottom panel) of α = 0.1 disks, under different numerical treatments. The thin cyan line shows the original SG profile as published in Sirko & Goodman 2003, valid only for log M•/M⊙ = 8.0 and α = 0.01. Thicker colorful lines are extensions of the model to different SMBH masses (orange: log M•/M⊙ = 6.0; purple: 7.0; black: 8.0). Different line-styles distinguish the numerical method: solid lines are self-consistent solutions of the SG equations computed with pAGN, while dashed and dotted lines represent power-law scaling approximations from Vaccaro et al. (2024) and Grishin et al. (2024) respectively. The dot-dashed gray line is a Shakura & Sunyaev (1973) disk for comparison. Solid lines should be regarded as the most physically consistent profiles within this framework.

Fig. A.2. Same as Figure A.1, but for α = 0.01.

Fig. A.3. Two-dimensional histograms showing the distribution of BBH pair-up locations in the K18 model in the radius–mass plane (R vs M•), for different combinations of torque prescription (top: G24; bottom: B16) and disk viscosity α (left: 0.01; right: 0.1), for a SG disk model. The color map indicates the number of pair-ups in each bin, Npair. Pink dots mark the locations of migration traps in the K18 model, while red crosses mark the Type I trap positions (as in Figure 2). Pink contours show the values of the normalized contrast parameter Knorm.

In this appendix, we compare our current prescription with models of α-disks previously used in the literature. An α-disk is a model in which gas turbulence is assumed to be the cause of disk viscosity, characterized by the viscosity coefficient α ∈ [0,  1]. The prototypical α-disk is the one by Shakura & Sunyaev (1973), where crucial quantities such as the gas surface density Σ_g and aspect ratio h are modeled as power laws. We show the resulting profile in A.1 and A.1, using the Kocsis et al. (2011, eqs. 32–33) treatment and assuming $\dot{M_{•}}/{\dot{M_{•}}}_{,\mathrm{Edd}}=0.1$. The obvious problem with this disk model is that becomes unphysical in the outer regions, and it is hence only valid for R ≲ 10³R_S. Moreover, such disks are viscously, thermally, and convectively unstable to perturbations, meaning that they can easily fragment into smaller sub-clouds (Kocsis et al. 2011).

Sirko & Goodman (2003) extend the Shakura & Sunyaev (1973) model by assuming that star formation in the outskirts of the disk acts as stabilizing mechanism against fragmentation. Moreover, the SG model is also broadly consistent with observed quasar spectral energy distributions, as it more accurately reproduces the flat UV–optical slopes and IR excesses characteristic of luminous quasars. he full model consists of 16 coupled equations describing the radial structure of the disk, which were originally solved self-consistently for a fiducial set of parameters: log M_•/M_⊙ = 8.0, α = 0.01, and $\dot{M_{•}}/{\dot{M_{•}}}_{\mathrm{Edd}}=0.1$. Before the publication of the pAGN tool, various approximations were required to generalize this solution to other physical regimes. In the following, we illustrate two examples.

In Vaccaro et al. (2024, eqs. 1–2), we fit the original profile with broken power laws, and combine it with a first-order rescaling in M_• and α. As shown in A.1 and A.1, this provides quite a crude modeling of the relevant quantities, as it neglects the complex interdependency between parameters in the full system of equations. The approximation breaks down most noticeably at high α and high M_•, where it severely overestimates the surface density Σ_g in the inner regions of the disk.

Another approach is the one by Grishin et al. (2024, eqs. 3–11) who reengineer an analytical model which is similar in spirit to SG, but neglecting some of the physics and using simplified analytical expressions. The resulting profiles, also shown in A.1 and A.1 for $\dot{M_{•}}/{\dot{M_{•}}}_{\mathrm{Edd}}=0.1$, are qualitatively similar to the SG ones computed with pAGN (meaning that they feature an inner region where the density increases and the thickness decreases for growing radii, and an outer region where the opposite happens), but the transition radii are systematically one order of magnitude smaller than in the more accurate SG paradigm. Moreover, they also overestimate the density and underestimate the aspect ratio in the inner regions (R ≲ 10³R_S).

This highlights the importance of solving the full set of equations with pAGN, as the dependence of the disk structure on α and M_• is highly nonlinear and cannot be captured by simple rescaling or analytical fits, when constructing models of BH dynamics in AGN disks.

Appendix B: Type II migration

The classical picture of Type II migration assumes that, once a gap is opened, the embedded object becomes locked to the viscous inflow of the disk with migration timescale t_type II = t_visc, as assumed in subsection 2.2. This disk-dominated view is conservative: it assumes the gap to be completely cleared and gas inflow across the gap to be suppressed. Recent hydrodynamical simulations have shown that this picture is oversimplified: in reality, gap opening is often partial and gas can leak through the gap edges, leading to migration rates that are faster than the purely viscous timescale (e.g. Malik et al. 2015).

Kanagawa et al. (2018) (hereafter K18) proposed a modified prescription that interpolates between Type I and Type II regimes, so that migration is well described as a modified Type I regime with reduced local density

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{g}}\prime (R)=\frac{{\mathrm{\Sigma }}_{\mathrm{g}}(R)}{1+0.04K},\end{array}$$(B.1)

where the contrast parameter K is defined as in eq. 13. For low-mass migrators (K ≲ 1, Σ_g′≳0.96 Σ_g), the gap is shallow and standard Type I migration (eq. 1) remains valid. As the migrator mass increases, the contrast parameter K grows quadratically. The full transition to a Type II regime happens for K ≳ 20 (K18), as the disk surface density at the orbit is significantly reduced to Σ_g′≲0.55 Σ_g.

We calculate corotational gas torques in eqs. 3, 4, 5, 6, 7, and 12 using the reduced local density Σ_g′. In this regime, both the normalization of the torque Γ₀ (eq. 12), and the shape of the torque profiles (via the parameter a in eq. 4) are significantly modified. We include an exponential decay term which scales with the contrast factor K, as in K18. We leave the Lindblad torque (Γ_L, eq. 8) unperturbed.

The effect of gap opening of Lindblad torques is quite complex, as it depends on the exact profile of the gap (e.g. Kley & Nelson 2012; Petrovich & Rafikov 2012). A more comprehensive simulation would include modeling the gap’s density profile and computing Lindblad torques self-consistently. However, in the high-viscosity regime in which we are operating, we do not expect any nonlinear behavior that would strongly affect our results (Duffell 2015).

Explicitly, our torque equations read

$$\begin{array}{c}\hfill \mathrm{\Gamma }\prime =[\left(\frac{0.46-0.96a\prime +1.8b}{\gamma }\right)\mathrm{\Gamma }{\prime }_0+\mathrm{\Gamma }{\prime }_{\mathrm{th}}]e^{-K/20}+{\mathrm{\Gamma }}_{\mathrm{L}}\end{array}$$(B.2)

in the G24 case and

$$\begin{array}{c}\hfill \mathrm{\Gamma }\prime =\frac{\mathrm{\Gamma }{\prime }_{\mathrm{ad}}{(\mathrm{\Theta }\prime )}^2+\mathrm{\Gamma }{\prime }_{\mathrm{iso}}}{{(\mathrm{\Theta }\prime +1)}^2}{e}^{-K/20}\end{array}$$(B.3)

in the B16 case, where Γ′_i = Γ_i(Σ_g′) and a′= − ∂log Σ′_g/∂log R.

An important implication is that migration traps are no longer determined solely by the background disk structure, as the torques Γ′=Γ(Σ_g′) are now K-dependent and, consequently, m-dependent. This mass dependence adds an additional layer of complexity to the dynamics of BHs in AGN disks. In the following, the trap locations under the K18 prescription are defined as the radii where the total migration torque Γ′ acting on a migrator of average mass (m ≃ 13.5 M_⊙) changes sign from positive to negative.

While the gap contrast naturally depends on the properties of the migrator, a disk-specific contrast parameter K_norm can be defined by normalizing over the migrator mass:

$$\begin{array}{c}\hfill K_{\mathrm{norm}}=\frac{K}{m^2}={\alpha }^{-1}{h}^{-5}{\left(\frac{M_{•}}{{\mathrm{M}}_{\odot }}\right)}^{-2}\end{array}$$(B.4)

We recompute the binary formation catalog using the K18 prescription for gap opening. The resulting pair-up distributions are shown in Figure A.3. Compared to the results presented in subsection 3, the locations of migration traps differ substantially for low-mass SMBHs (M_• ≲ 10⁷ M_⊙), because the modified torques tend to shift the zero-torque radii away from regions of high contrast K_norm.

In the G24 scenario, the peaks in the pair-up radius distributions are generally located in proximity to the new migration traps, with the exception of systems with M_• ≲ 10⁶ M_⊙ which present traffic-jam accumulations of pair-up events in regions of low contrast (K_norm ≲ 1), even in the absence of traps.

In the B16 framework, the Lindblad torque is neglected (see eqs. 2 and B.3). Hence, the migration torque is highly suppressed in areas of high contrast (K_norm ≳ 1) and almost no pair-up events happen at radii R ∈ [10², 10⁵]R_g for M_• ≤ 10^5.8 M_⊙. Moreover, disks with SMBH masses up to 10^6.5 M_⊙ present prominent traffic-jams in the outer area of the disks.

Figure B.1 shows the resulting pair-up efficiency, η, and trap fraction, f_trap, obtained with the K18 model.

Fig. B.1. Fraction of pair-ups occurring within migration traps (left panels) and pair-up efficiency (right panels) as functions of SMBH mass M• (top row), migrator mass m (middle row), and initial radius Rinit (bottom row) in the K18 prescription. Color encodes the viscosity parameter α (orange: 0.01; teal: 0.1; navy blue: 0.4), and the two torque prescriptions are distinguished by marker shape and line style (squares and solid lines: G24; circles and dashed lines: B16).

The pair-up efficiency, η, is generally lower than in the fiducial runs, peaking at roughly 0.3 for M_• ≲ 10⁶ M_⊙ in the low-viscosity G24 case, while retaining similar trends in SMBH mass M_• and in initial radius R_init to those observed in subsection 3. On the other hand, the dependency on migrator mass m is flatter than in the fiducial case especially for disks with low viscosity (α = 0.01), implying that partial gap opening (K18) moderately suppresses the pair-up efficiency of larger mass migrators (with m ≳ 20 M_⊙).

The fraction of pair-ups occurring within migration traps, f_trap, retains a similar behavior than in Figure 3 for the high-viscosity G24 case, remaining large (> 0.6) for most of the parameter space. In the low-viscosity G24 case, f_trap is large for M_• ∈ [10⁷, 10^8.7] M_⊙ and ≲0.3 otherwise. Contrarily, in the B16 scenario, traffic-jam pair-ups are more relevant, especially at low viscosity for low SMBH mass. Therefore, f_trap < 0.2 with the exception of a certain α-dependent window: log M_•/M_⊙ ∈ [6.5, 8.6] for α = 0.01, log M_•/M_⊙ ∈ [6.0, 7.6] for α = 0.1, or log M_•/M_⊙ ∈ [5.8, 7.3] for α = 0.4. It retains a similar trend in both migrator mass m and initial radius R_init than discussed in subsection 3.4.

We compute pair-up times, t_pair, and evaluate how the K18 framework affects them. These pair-up times, shown via the 68% and 90% contour levels in Figure B.2, represent the elapsed time between the moment at which the single primary BH becomes embedded in the disk and the subsequent formation of a BBH. At low SMBH mass, M_• ≲ 3 × 10⁶ M_⊙, the normalized contrast parameter can reach larger values, K_norm ≳ 1, thus significantly reducing the local gas density Σ′ and slowing down the pair-up process, therefore increasing t_pair by almost two orders of magnitude for the lowest SMBH mass explored, M_• = 10⁵ M_⊙. In contrast, for M_• ≳ 3 × 10⁶ M_⊙, pair-up times are unaffected by the adopted Type II prescription.

Fig. B.2. Pair-up times tpair as a function of SMBH mass in SG AGN disks, for different combinations of the viscosity parameter α (left: 0.01; right: 0.1) and torque prescriptions (top: G24; bottom: B16). Each panel shows a median line and contours enclosing 68% and 90% of the pair-up time distribution, based on our simulations. Results are shown for two different Type II migration prescriptions, indicated by color and line style. In blue, dashed lines: full gap opening, as in the main text (eq. 15). In red, continuous lines: partial gap opening, as in this appendix (K18).

Overall, the adoption of the K18 prescription has different consequences depending on the torque model adopted.

In the G24 case, the position of traps is shifted toward regions of low contrast K_norm. For M_• ≲ 10⁶(10⁷) M_⊙ in the high-viscosity (low-viscosity) case, pair-up events happen in traffic jams in valleys of low contrast K_norm. For larger SMBH masses, both the pair-up efficiency and the relevance of migration traps are mostly unaffected.

In the B16 case, the pair-up process in disks with low SMBH mass (log M_•/M_⊙ ≲ 6) is strongly suppressed, and pair-up events often happen in the outskirts of the disk.

Consequently, the self-consistent treatment of gap opening lowers the pair-up efficiency of high-mass migrators and weakens the correlation between migration traps and binary formation sites.

All Figures

Fig. 1. Examples of torque profiles as a function of radial distance (in Schwarzschild radii) for BHs embedded in a SG disk. Each plot corresponds to a different combination of torque prescription (top: G24; bottom: B16), disk viscosity parameter α (left: α = 0.01; right: α = 0.1) and SMBH mass (log M•/M⊙ = 5.0, 6.3, 7.7, 9.0, represented in different colors). Solid and dashed lines indicate positive and negative net torques, respectively, while shaded regions denote uncertainty in the torque determined by computing it over a range of migrator masses between 5 M⊙ and 50 M⊙. Vertical dotted lines mark the locations of migration traps.

In the text

	Fig. 2. Two-dimensional histograms showing the distribution of BBH pair-up locations in the radius–mass plane (R vs. M•), for different combinations of torque prescription (top: G24; bottom: B16) and disk viscosity α (left: 0.01; right: 0.1) in a SG disk. The color map indicates the number of pair-ups in each bin Npair. Red crosses mark the locations of migration traps.
In the text

Fig. 3. Fraction of binary pair-ups occurring within migration traps as a function of the SMBH mass. The vertical axis shows the occurrence ftrap of pair-up events that take place in proximity to predefined migration trap locations, while the horizontal axis shows the logarithm of the central SMBH mass. Different marker fillings indicate the BH generation (filled: 1g; hollow: Ng), while color encodes the viscosity parameter α (orange: 0.01; teal: 0.1; navy blue: 0.4). The two torque formalisms are distinguished by marker shape and line style (squares and solid lines: G24; circles and dashed lines: B16). Trap proximity is defined using a logarithmic tolerance on radial separation (see Section 3.2 for details).

In the text

Fig. 4. Histograms of binary pair-up radii Rpair (in units of Schwarzschild radii) for first-generation (1g) and Nth-generation (Ng) BHs, across a range of SMBH masses (log M•/M⊙ = 6, 7, 8 in the three columns, respectively) and viscosity parameters (upper row: α = 0.01; lower row: α = 0.1). Panels show the pair-up radii probability density functions (PDFs), normalized separately in each panel. The 1g distributions are shown as step histograms, while Ng distributions are overlaid either as filled or hatched histograms for visual distinction. Vertical dotted lines indicate the positions of migration traps identified from the torque profile. Color represents the migration model used (green: G24; purple: B16).

In the text

Fig. 5. Pair-up efficiency as a function of the central SMBH mass. The vertical axis shows the fraction of BHs that successfully form a binary during the simulation, while the horizontal axis shows the logarithm of the SMBH mass. As in Figure 3, different marker fillings indicate the BH generation (filled: 1g; hollow: Ng), color encodes the viscosity parameter α (orange: 0.01; teal: 0.1; navy blue: 0.4), and the two torque prescriptions are distinguished by marker shape and line style (squares and solid lines: G24; circles and dashed lines: B16).

In the text

Fig. 6. Pair-up efficiency (top panels) and fraction of pair-ups occurring within migration traps (bottom panels) as functions of initial radius Rinit (left) and migrator mass m (right). Color encodes the viscosity parameter α (orange: 0.01; teal: 0.1; navy blue: 0.4), and the two torque prescriptions are distinguished by marker shape and line style (squares and solid lines: G24; circles and dashed lines: B16).

In the text

Fig. 7. Two-dimensional distributions of the migration time tmigr as a function of the initial radius Rinit in SG AGN disks, for different combinations of the viscosity parameter α (left: 0.01; right: 0.1) and torque prescriptions (top: G24; bottom: B16). Each panel shows contours enclosing 40% and 68% of the migration time distribution, based on our simulations. Results are shown for four representative SMBH masses: log(M•/M⊙) = 5.0, 6.3, 7.7, 9.0, indicated by color. Solid (dashed) lines represent the migration timescale computed using the torque Γ in Eq. (1) (Γ0 in Eq. 12). Vertical dotted lines mark the locations of migration traps.

In the text

Fig. A.1. Radial profiles of surface density Σg (top panel) and aspect ratio h = H/R (bottom panel) of α = 0.1 disks, under different numerical treatments. The thin cyan line shows the original SG profile as published in Sirko & Goodman 2003, valid only for log M•/M⊙ = 8.0 and α = 0.01. Thicker colorful lines are extensions of the model to different SMBH masses (orange: log M•/M⊙ = 6.0; purple: 7.0; black: 8.0). Different line-styles distinguish the numerical method: solid lines are self-consistent solutions of the SG equations computed with pAGN, while dashed and dotted lines represent power-law scaling approximations from Vaccaro et al. (2024) and Grishin et al. (2024) respectively. The dot-dashed gray line is a Shakura & Sunyaev (1973) disk for comparison. Solid lines should be regarded as the most physically consistent profiles within this framework.

In the text

	Fig. A.2. Same as Figure A.1, but for α = 0.01.
In the text

Fig. A.3. Two-dimensional histograms showing the distribution of BBH pair-up locations in the K18 model in the radius–mass plane (R vs M•), for different combinations of torque prescription (top: G24; bottom: B16) and disk viscosity α (left: 0.01; right: 0.1), for a SG disk model. The color map indicates the number of pair-ups in each bin, Npair. Pink dots mark the locations of migration traps in the K18 model, while red crosses mark the Type I trap positions (as in Figure 2). Pink contours show the values of the normalized contrast parameter Knorm.

In the text

Fig. B.1. Fraction of pair-ups occurring within migration traps (left panels) and pair-up efficiency (right panels) as functions of SMBH mass M• (top row), migrator mass m (middle row), and initial radius Rinit (bottom row) in the K18 prescription. Color encodes the viscosity parameter α (orange: 0.01; teal: 0.1; navy blue: 0.4), and the two torque prescriptions are distinguished by marker shape and line style (squares and solid lines: G24; circles and dashed lines: B16).

In the text

Fig. B.2. Pair-up times tpair as a function of SMBH mass in SG AGN disks, for different combinations of the viscosity parameter α (left: 0.01; right: 0.1) and torque prescriptions (top: G24; bottom: B16). Each panel shows a median line and contours enclosing 68% and 90% of the pair-up time distribution, based on our simulations. Results are shown for two different Type II migration prescriptions, indicated by color and line style. In blue, dashed lines: full gap opening, as in the main text (eq. 15). In red, continuous lines: partial gap opening, as in this appendix (K18).

In the text