EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 703 (November 2025) A&A, 703 (2025) A6 Full HTML
Open Access
Issue
A&A
Volume 703, November 2025
Article Number A6
Number of page(s) 10
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/202555879
Published online 31 October 2025

© The Authors 2025

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1. Introduction

Gravitational waves (GWs) are rapidly becoming a powerful tool to probe the complex dynamics of dense stellar environments and the physics of compact objects lurking at their centers (Mapelli 2021; Spera et al. 2022; The LIGO Scientific Collaboration et al. 2023; Abbott et al. 2023). Upcoming ground-based detectors such as the Einstein Telescope (ET, Punturo et al. 2010) and the Cosmic Explorer (CE, Reitze et al. 2019) will significantly extend sensitivity to lower frequencies and larger volumes, enabling the detection of black hole (BH) mergers out to high redshift and providing improved access to the early inspiral phase of stellar-mass binaries (Maggiore et al. 2020; Branchesi et al. 2023). Complementing these efforts, future space-borne GW detectors such as LISA (Amaro-Seoane et al. 2017) or TAIJI (Hu & Wu 2017) will open the sub-Hz band, allowing the detection of long-lived inspirals involving both supermassive BH (SMBH) binaries and stellar-mass BHs in extreme mass-ratio configurations (Hopman & Alexander 2005; Amaro-Seoane et al. 2007; Mandel et al. 2008; Babak et al. 2017).

As GW detectors reach higher sensitivity, attention is turning to the role of environmental effects in shaping binary evolution (Zwick et al. 2024, 2025; Takátsy et al. 2025). Among these, gas-rich environments such as active galactic nuclei (AGNs) offer particularly promising prospects. Interactions with AGN accretion disks can significantly alter the dynamics of embedded BH binaries through hydrodynamic drag, accretion torques, and disk-driven migration, potentially accelerating inspiral or exciting orbital eccentricity (Stone et al. 2016; Ishibashi & Gröbner 2020; Li et al. 2021, 2023; Li & Lai 2022, 2023, 2024; Dempsey et al. 2022; Samsing et al. 2022; Rowan et al. 2023, 2024a,b, 2025a; Whitehead et al. 2024a,b, 2025a,b; Dittmann et al. 2024, 2025; Dodici & Tremaine 2024; Mishra & Calcino 2024; Trani et al. 2024). These effects may leave detectable imprints on the GW signal, opening a window onto the astrophysical environments in which BH mergers occur. Other mechanisms, such as dynamical friction from dark matter spikes, have also been thought to influence binary evolution in galactic nuclei (e.g., Hannuksela et al. 2020; Montalvo et al. 2024; Mukherjee et al. 2024; Dosopoulou 2024; Fischer & Sagunski 2024; Kavanagh et al. 2025). However, AGN-assisted mergers are particularly compelling because they may also produce electromagnetic counterparts through shocks, variable accretion, or relativistic jets, thus enabling multi-messenger observations (Graham et al. 2020; Ren et al. 2022; Gayathri et al. 2023; Tagawa et al. 2023, 2024; Chen & Dai 2024).

Magnetohydrodynamical simulations have revealed that disk dynamics, ranging from protoplanetary to AGN scales, are largely governed by turbulence driven by magnetorotational instability (Balbus & Hawley 1998; Janiuk et al. 2004; Armitage 2011; Zubovas et al. 2024). This turbulence can induce substantial density and velocity fluctuations, which can have a strong impact on the orbital evolution of embedded objects (Papaloizou et al. 2004; Nelson & Papaloizou 2004; Nelson 2005; Johnson et al. 2006; Oishi et al. 2007; Yang et al. 2009, 2012; Nelson & Gressel 2010).

From the perspective of kinetic (i.e., particle based) simulations, only limited attention has been paid to the impact of the density fluctuations on the orbit of objects embedded in turbulent gaseous disks (Rein & Papaloizou 2009; Baruteau & Lin 2010; Secunda et al. 2019), and even less effort has been devoted to modeling the direct impact of turbulent velocity fields. A key obstacle lies in the broad range of scales involved: hydrodynamic turbulence evolves rapidly on small spatial scales, whereas orbital motion proceeds more slowly and over larger distances. This disparity makes it computationally challenging to simultaneously resolve both the trajectory of a compact object and the collective dynamics of the surrounding disk within a single hybrid simulation.

In this work, we investigated the role of turbulence in shaping the orbital evolution of a stellar-mass BH embedded in the accretion disk of an SMBH, focusing on orbits near the migration trap radius. To this end, we developed a semi-analytical framework based on a stochastic differential equation that captures turbulent velocity fluctuations through a time-varying drag coefficient.

The paper is organized as follows. In Sect. 2 we introduce the AGN disk model, the governing equations, and the numerical techniques employed. In Sect. 3 we present the results of our semi-analytical calculations and construct a reduced model to interpret the key features. Finally, in Sect. 4 we summarize our findings and discuss their astrophysical implications.

2. Models

2.1. The turbulent AGN disk model

In our kinetic model, we account for the effect of a turbulent AGN disk via an effective friction and diffusion process. The equation of motion for the secondary BH of stellar mass (mBH) orbiting in the gravitational field of the SMBH MBH reads

d 2 r d t 2 = − G M BH r 3 r − η ( v ∼ , R , z ) v ∼ , $$ \begin{aligned} \frac{\mathrm{d^2}\mathbf r }{\mathrm{d}t^2}=-\frac{GM_{\rm BH}}{r^3}\mathbf r -\eta (\tilde{v},R,z)\tilde{\mathbf{v }}, \end{aligned} $$(1)

where η is the local (i.e., dependent on the cylindrical radial position and vertical distance z) friction coefficient, accounting for the drag force exerted by the disk. The relative velocity, v ∼ $ \tilde{\mathbf{v}} $, is defined as

v ∼ = v − v circ − v turb . $$ \begin{aligned} \tilde{\mathbf{v }}=\mathbf v -\mathbf v _{\rm circ}-\mathbf v _{\rm turb}. \end{aligned} $$(2)

In the equation above, v = r ˙ $ \mathbf{v}=\dot{\mathbf{r}} $; vcirc is the gas circular velocity at cylindrical radius R, of magnitude v circ = G M BH / R $ v_{\mathrm{circ}}=\sqrt{GM_{\mathrm{BH}}/R} $. vturb is the stochastic fluctuation of vcirc induced by the turbulent gas motion in the disk. The drag coefficient is defined by

η = 4 π G 2 ρ m BH I ( M ) v ∼ 3 , $$ \begin{aligned} \eta = 4\pi G^2\rho m_{\rm BH} \frac{I(\mathcal{M} )}{\tilde{v}^3}, \end{aligned} $$(3)

where the quantity I is given as a function of the local Mach number M = v ∼ / c s $ \mathcal{M}=\tilde{v}/{{c}_\text{s}} $ by the Ostriker (1999) prescription

I = { 1 2 log 1 + M 1 − M − M M ≤ 0.95 , [ 4 p t ] 1 2 log ( 1 − 1 M 2 ) + 3.1 M ≥ 1.05 , [ 2 p t ] 0.88 + 10.3 ( M − 0.95 ) elsewhere . $$ \begin{aligned} I= {\left\{ \begin{array}{ll} \displaystyle \frac{1}{2}\log \frac{1+\mathcal{M} }{1-\mathcal{M} }-\mathcal{M} \quad \mathcal{M} \le 0.95,\\ [4pt] \displaystyle \frac{1}{2}\log \Big ( 1-\frac{1}{\mathcal{M} ^2}\Big )+3.1\quad \mathcal{M} \ge 1.05,\\ [2pt] \displaystyle 0.88+10.3\Big (\mathcal{M} -0.95\Big )\quad \mathrm{elsewhere}. \end{array}\right.} \end{aligned} $$(4)

In the definition above, the third expression is a linear interpolation between the subsonic and supersonic regimes1 to overcome the mathematical discontinuity at ℳ = 1 (see Fig. 1 and DeLaurentiis et al. 2023).

thumbnail Fig. 1.

Dependency of Ostriker (1999) the dynamical friction prescription as a function of the mach number ℳ (Eq. (4)). The dotted line marks the linear interpolation used to circumvent the mathematical discontinuity at ℳ = 1. The dashed line indicates the linear subsonic limit.

We note that the model given by Eq. (1) cannot be reduced in terms of a simple second order Langevin equation (see, e.g., the discussion in Sartorello et al. 2025) of the form

r ¨ = − ∇ Φ − η r ˙ + δ f , $$ \begin{aligned} \ddot{\mathbf{r }}=-\nabla \Phi -\eta \dot{\mathbf{r }}+\delta \mathbf f , \end{aligned} $$(5)

where δf is a fluctuating force per unit mass, as the velocity fluctuations induced by turbulence also enter explicitly the definition of the friction coefficient η.

In this work we neglected the gravitational field of the AGN disk and assumed an axisymmetric density ρ(R, z) of the form

ρ ( R , z ) = ρ R exp ( − z 2 2 H R 2 ) , $$ \begin{aligned} \rho (R,z) = \rho _R\exp \left(-\frac{z^2}{2H_R^2}\right), \end{aligned} $$(6)

where ρR and HR are the midplane density and scale height at cylindrical radius R, respectively. In the numerical simulations discussed hereafter, we interpolated its parameters from a 1D AGN model generated with pagn (Gangardt et al. 2024). We adopted the Sirko & Goodman (2003) model with MBH = 107 M⊙ and the alpha viscosity α = 0.01, the Eddington luminosity ratio lE = 0.5, and the radiative efficiency ε = 0.1. The density profile of the disk was truncated at R = 2 × 107 Rg. Throughout this work we assumed mBH = 10 M⊙. Figure 2 shows the radial density in the midplane of the AGN disk (top panel, corresponding to ρR in Eq. (6)), its sound speed and scale height profiles (mid-panels), as well as the magnitude of the so-called migration torque Γ (bottom panel), defined as the effective force exerted on mBH by the leading and trailing spiral density waves induced by the perturbations of the disk density. The migration torque is evaluated as in Gangardt et al. (2024, see also Masset 2017; Guilera et al. 2021; Grishin et al. 2024). The model disk has two migration traps, Rtrap, 1 ≃ 121 au ≃ 1230 Rg and Rtrap, 2 ≃ 939 au ≃ 9512 Rg, indicated in the figure by the vertical dashed lines. They are defined as the radii at which the torque Γ is such that, for R < Rtrap, i, the particle mBH is pushed outward, and inward otherwise. For reasons of simplicity, as we only considered initial conditions starting at migration traps, we neglected the effect of the migration torque in our model.

thumbnail Fig. 2.

Radial profiles of the AGN disk employed in this work. From top to bottom: midplane gas density ρ, sound speed cs, disk aspect ratio h/R, and magnitude of the migration torque Γ. The migration torque assumes a secondary BH mass of mBH = 20 M⊙. The outer migration trap lies within the star formation region, which begins at R ≃ 767 au.

2.2. Governing equations

Starting from Eq. (1) in Cartesian coordinates, we derived the Gaussian perturbation equations in terms of Keplerian orbital elements using standard perturbative techniques (see Beutler 2005). All numerical results in this work are based on the full form of these equations, which we independently validated by comparing with direct integrations of Eq. (3) in Cartesian coordinates (see Sect. 2.3). To develop analytical insight, however, we also considered simplified versions of the perturbation equations that capture the essential physics while remaining more tractable. Refer to Table 1 for the full list of symbols used in this manuscript.

Table 1.

Symbols used in the manuscript.

Rather than using the full expression for η ( v ∼ ) $ \eta(\tilde{v}) $, which depends on a complex nonlinear dependence on the relative velocity, we took advantage of the fact that the orbits under consideration are fully embedded within the AGN disk. These orbits typically have low eccentricities (e) and inclinations (ι). To first order in e and ι, the relative velocity (excluding contributions from the turbulence) are approximated as

v ∼ 2 ≃ v circ 2 [ e 2 ( 1 − 3 4 cos 2 ν ) + ι 2 cos 2 u ] , $$ \begin{aligned} \tilde{v}^2 \simeq v_\mathrm{circ} ^2 \left[ e^2 \left(1 - \frac{3}{4} \cos ^2{\upnu }\right) + \iota ^2 \cos ^2{u} \right], \end{aligned} $$(7)

where ν is the true anomaly and u = ω + ν is the argument of latitude. For a stellar-mass BH embedded in the disk with e, ι ≪ 1, the relative velocity satisfies v ∼ ≪ v circ = c s ( H / R ) − 1 $ \tilde{v} \ll v_{\mathrm{circ}} = c_{\mathrm{s}} (H / R)^{-1} $, where the expression for vcirc follows from vertical hydrostatic equilibrium.

Given that the disk aspect ratio satisfies H/R ≳ 0.01 everywhere (see Fig. 2), we obtain v ∼ ≪ c s $ \tilde{v} \ll {{c}_\text{s}} $, and the gas friction is always subsonic. In this regime, for Mach numbers ℳ ≪ 1 (a condition that holds up to ℳ ≃ 0.5; see Fig. 1), the drag coefficient in Eq. (3) becomes independent of v ∼ $ \tilde{v} $ and simplifies to

η ∗ = 4 π G 2 ρ m BH 3 c s 3 · $$ \begin{aligned} \eta _* = \frac{4\pi G^2 \rho \,{m}_\mathrm{BH} }{3 {c}_\mathrm{s} ^3}\cdot \end{aligned} $$(8)

Here, ρ retains its dependence on the radial and vertical coordinates R and z, as given by Eq. (6). For convenience, we also define the associated damping timescale τ*:

τ ∗ = 1 η ∗ = 3 c s 3 4 π G 2 ρ m BH · $$ \begin{aligned} \tau _* = \frac{1}{\eta _*} = \frac{3 {c}_\mathrm{s} ^3}{4\pi G^2 \rho \,{m}_\mathrm{BH} }\cdot \end{aligned} $$(9)

Consequently, in the subsonic regime, the drag force becomes linear in the velocity v ∼ $ \tilde{v} $, allowing us to decompose the acceleration appearing in Eq. (1) into the deterministic component

d 2 r d t 2 | det = − G M BH r 3 r − η ∗ ( v − v circ ) , $$ \begin{aligned} \left.\frac{\mathrm{d^2}\mathbf r }{\mathrm{d}t^2}\right|_{\rm det} = -\frac{GM_{\rm BH}}{r^3}\mathbf r -\eta _* (\mathbf v -\mathbf v _\mathrm{circ} ), \end{aligned} $$(10)

and the stochastic fluctuating term

d 2 r d t 2 | sto = η ∗ v turb . $$ \begin{aligned} \left.\frac{\mathrm{d^2}\mathbf r }{\mathrm{d}t^2}\right|_{\rm sto} = \eta _* \mathbf v _{\rm turb}. \end{aligned} $$(11)

In practice, one can collect η* and vturb in Eq. (11) into an effective fluctuating force δf*, thereby basically recovering Eq. (5).

Using standard first-order perturbation theory, the Gaussian equations governing the evolution of the semimajor axis (a), the eccentricity (e), the inclination (ι), the argument of pericenter (ω), the longitude of pericenter (ϖ = ω + Ω), and the longitude of ascending node (Ω) derived from the deterministic component of the equations Eq. (10) read

d a d t = − 2 a η ∗ 1 − e 2 ( 1 + e 2 + 2 e cos ν − cos ι ( 1 + e cos ν ) 3 / 2 ( cos 2 ι + cos 2 u sin 2 ι ) 3 / 4 ) , $$ \begin{aligned} \frac{\mathrm{d}a}{\mathrm{d}t}&= -\frac{2 a \eta _*}{1-e^2} \left(1 + e^2 + 2e\cos \upnu - \frac{\cos \iota \, (1 + e \cos \upnu )^{3/2}}{(\cos ^2\iota + \cos ^2{u} \sin ^2\iota )^{3/4}}\right),\end{aligned} $$(12)

d e d t = − η ∗ ( cos ν + e + cos ν 1 + e cos ν ) · ( 1 + e cos ν − cos ι ( 1 + e cos ν ) 1 / 2 ( cos 2 ι + cos 2 u sin 2 ι ) 3 / 4 ) − e η ∗ sin 2 ν , $$ \begin{aligned} \frac{\mathrm{d}e}{\mathrm{d}t}&= -\eta _* \left(\cos \upnu + \frac{e+\cos \upnu }{1+e \cos \upnu } \right)\cdot \nonumber \\&\quad \left(1+e\cos \upnu -\frac{\cos \iota \,(1 + e \cos \upnu )^{1/2} }{(\cos ^2\iota + \cos ^2{u} \sin ^2\iota )^{3/4}}\right) - e \eta _*\sin ^2\upnu , \end{aligned} $$(13)

d ι d t = − η ∗ cos 2 u sin ι ( 1 + e cos ν ) 1 / 2 ( cos 2 u + cos 2 ι sin 2 u ) 3 / 4 , $$ \begin{aligned} \frac{\mathrm{d}\iota }{\mathrm{d}t}&= -\eta _*\frac{\cos ^2{u}\,\sin \iota }{(1+e \cos \upnu )^{1/2} (\cos ^2{u} + \cos ^2{\iota } \sin ^2{u})^{3/4}}, \end{aligned} $$(14)

d ϖ d t = − η ∗ sin u e ( 2 − cos ι ( 2 + e cos ν ) ( 1 + e cos ν ) 1 / 2 ( cos 2 ι + cos 2 u sin 2 ι ) 3 / 4 ) , $$ \begin{aligned} \frac{\mathrm{d}\varpi }{\mathrm{d}t}&= -\eta _* \frac{\sin {u}}{e}\left(2 - \frac{\cos \iota \,(2 + e\cos \upnu )}{(1 + e \cos \upnu )^{1/2} (\cos ^2\iota + \cos ^2{u} \sin ^2\iota )^{3/4}}\right), \end{aligned} $$(15)

d Ω d t = − η ∗ sin u cos u ( 1 + e cos ν ) 1 / 2 ( cos 2 u + cos 2 ι sin 2 u ) 3 / 4 · $$ \begin{aligned} \frac{\mathrm{d}\Omega }{\mathrm{d}t}&= -\eta _* \frac{\sin {u}\cos {u}}{(1+e\cos \upnu )^{1/2}(\cos ^2{u} + \cos ^2{\iota } \sin ^2{u})^{3/4}}\cdot \end{aligned} $$(16)

This system is closed by the usual two-body expressions for dν/dt and dω/dt, while ensuring that the reference frame contributions are taken into account (e.g., Burns 1976):

d ν d t = μ a 3 ( 1 + e cos ν ) 2 ( 1 − e 2 ) 3 / 2 − d ϖ d t , $$ \begin{aligned} \frac{\mathrm{d}\upnu }{\mathrm{d}t}&= \sqrt{\frac{\mu }{a^3}} \frac{(1 + e \cos {\upnu })^2}{(1-e^2)^{3/2}} - \frac{\mathrm{d}\varpi }{\mathrm{d}t},\end{aligned} $$(17)

d ω d t = d ϖ d t − cos ι d Ω d t · $$ \begin{aligned} \frac{\mathrm{d}\omega }{\mathrm{d}t}&= \frac{\mathrm{d}\varpi }{\mathrm{d}t} -\cos {\iota }\frac{\mathrm{d}\Omega }{\mathrm{d}t}\cdot \end{aligned} $$(18)

Here, μ = GMBH denotes the standard gravitational parameter.

To gain analytical insight, we expanded the perturbation equations to first order in e and ι. The resulting linearized system for the orbital elements reads:

d a d t = − η ∗ a e cos u , $$ \begin{aligned} \frac{\mathrm{d}a}{\mathrm{d}t}&= - \eta _* a e \cos {u}, \end{aligned} $$(19)

d e d t = − η ∗ e , $$ \begin{aligned} \frac{\mathrm{d}e}{\mathrm{d}t}&= - \eta _* e, \end{aligned} $$(20)

d ι d t = − η ∗ ι cos 2 u , $$ \begin{aligned} \frac{\mathrm{d}\iota }{\mathrm{d}t}&= -\eta _* \iota \cos ^2{u}, \end{aligned} $$(21)

d ϖ d t = 1 4 η ∗ e cos 2 ν sin ν , $$ \begin{aligned} \frac{\mathrm{d}\varpi }{\mathrm{d}t}&= \frac{1}{4} \eta _* e \cos ^2{\upnu } \sin {\upnu }, \end{aligned} $$(22)

d Ω d t = 1 2 η ∗ ( cos ν − 2 ) sin u cos u . $$ \begin{aligned} \frac{\mathrm{d}\Omega }{\mathrm{d}t}&= \frac{1}{2} \eta _* (\cos {\upnu } - 2) \sin {u}\cos {u}. \end{aligned} $$(23)

To model the stochastic component of the orbital evolution, we assumed that the turbulent velocity perturbation, vturb, follows a stationary stochastic Gaussian process with an autocorrelation function given by

⟨ v turb ( t 1 ) v turb ( t 2 ) ⟩ = G ( t 1 − t 2 ) = σ turb 2 exp ( − | t 1 − t 2 | τ c ) , $$ \begin{aligned} \langle v_{\rm turb}(t_1)\,v_{\rm turb}(t_2) \rangle = G(t_1-t_2) = \sigma _{\rm turb}^2\, \exp {\left(-\frac{|t_1 - t_2|}{\tau _\mathrm{c} }\right)}, \end{aligned} $$(24)

where σturb2 ≡ ⟨vturb2⟩ is the variance of the turbulent velocity field, and τc = 1/Ωcirc is the autocorrelation time, with Ωcirc denoting the local gas circular frequency at the position of the stellar-mass BH mBH. This assumption is consistent with previous models of turbulence in protoplanetary disks (e.g., Rein & Papaloizou 2009; Picogna et al. 2018). The amplitude of vturb is directly related to the local kinematic viscosity of the disk, which within the framework of α-disks, is given by

ν kin = α c s H . $$ \begin{aligned} \nu _{\rm kin} = \alpha {c}_\mathrm{s} H. \end{aligned} $$(25)

It follows that the 1D variance of the turbulent velocity field σturb2 field reads

σ turb 2 = ν kin τ c = α c s H Ω circ = α c s 2 , $$ \begin{aligned} \sigma _\mathrm{turb} ^2 = \frac{\nu _{\rm kin}}{\tau _\mathrm{c} } = \alpha {c}_\mathrm{s} H \,\Omega _{\rm circ} = \alpha {c}_\mathrm{s} ^2, \end{aligned} $$(26)

where the last identity uses H = cs/Ωcirc, which follows from vertical hydrostatic equilibrium.

The turbulent velocity vector in Cartesian coordinates is expressed as

v turb = ( v turb , r cos ϕ − v turb , ϕ sin ϕ v turb , r sin ϕ + v turb , ϕ cos ϕ v turb , z ) , $$ \begin{aligned} \mathbf v _{\rm turb}&= \begin{pmatrix} v_{\mathrm{turb}, r} \cos {\phi } - v_{\mathrm{turb}, \phi } \sin {\phi } \\ v_{\mathrm{turb}, r} \sin {\phi } + v_{\mathrm{turb}, \phi } \cos {\phi } \\ v_{\mathrm{turb}, z} \end{pmatrix}, \end{aligned} $$(27)

where ϕ is the cylindrical azimuth of the BH, i.e., the angle in the disk plane from the +x axis to the projection of the BH position onto the x–y plane, such that vturb, r and vturb, ϕ are the radial and azimuthal components of the turbulent velocity at the BH location, and vturb, z is the vertical component. We assumed that the radial, azimuthal, and vertical components of vturb were independent, stationary Gaussian processes, each with variance σturb2.

Deriving the Gaussian perturbation equations corresponding to the stochastic acceleration term in Eq. (11) is cumbersome but straightforward. The full set of equations reads:

d a d t = η ∗ a 3 μ ( 1 − e 2 ) { [ e sin ν ( cos 2 ι + cos 2 u sin 2 ι ) 1 / 2 ] v turb , r + cos ι ( 1 + e cos ν ) ( cos 2 u + cos 2 ι sin 2 u ) 1 / 2 v turb , ϕ + sin ι ( e cos ω + cos u ) v turb , z } , $$ \begin{aligned} \frac{\mathrm{d}a}{\mathrm{d}t}&= \eta _* \sqrt{\frac{ a^3}{\mu (1-e^2)}} \Bigg \{ \left[ e \sin \upnu \,(\cos ^2\iota + \cos ^2{u} \sin ^2\iota )^{1/2} \right] v_{\mathrm{turb}, r} \nonumber \\&\quad + \frac{\cos \iota \, (1 + e\cos \upnu )}{(\cos ^2{u} + \cos ^2\iota \sin ^2{u})^{1/2}} v_{\mathrm{turb}, \phi } \nonumber \\&\quad + \sin \iota \,(e\cos \omega + \cos {u}) \,v_{\mathrm{turb}, z} \Bigg \}, \end{aligned} $$(28)

d e d t = η ∗ a ( 1 − e 2 ) μ { [ sin ν ( cos 2 ι + cos 2 u sin 2 ι ) 1 / 2 + − ( cos ν + e + cos ν 1 + e cos ν ) sin 2 ι sin u cos u ( cos 2 u + cos 2 ι sin 2 u ) 1 / 2 ] v turb , r + [ cos ι ( cos 2 u + cos 2 ι sin 2 u ) 1 / 2 ( cos ν + e + cos ν 1 + e cos ν ) ] v turb , ϕ + [ sin ν sin u + ( cos ν + e + cos ν 1 + e cos ν ) cos u ] sin ι v turb , z } , $$ \begin{aligned} \frac{\mathrm{d}e}{\mathrm{d}t}&= \eta _* \sqrt{\frac{a (1-e^2)}{\mu }} \Bigg \{ \Bigg [ \sin \upnu \,(\cos ^2\iota + \cos ^2{u} \sin ^2\iota )^{1/2} + \nonumber \\ &\quad - \left(\cos \upnu + \frac{e+\cos \upnu }{1+e\cos \upnu }\right) \frac{\sin ^2\iota \sin {u}\cos {u}}{(\cos ^2{u} + \cos ^2\iota \sin ^2{u})^{1/2}} \Bigg ] v_{\mathrm{turb}, r} \nonumber \\&\quad + \left[\frac{\cos \iota }{(\cos ^2{u} + \cos ^2\iota \sin ^2{u})^{1/2}} \, \left(\cos \upnu + \frac{e+\cos \upnu }{1+e\cos \upnu }\right)\right] v_{\mathrm{turb}, \phi } \nonumber \\&\quad + \Bigg [\sin \upnu \sin {u} + \left(\cos \upnu +\frac{e+\cos \upnu }{1+e\cos \upnu }\right)\cos {u}\Bigg ] \sin \iota \, v_{\mathrm{turb}, z} \Bigg \}, \end{aligned} $$(29)

d ι d t = − η ∗ a ( 1 − e 2 ) μ cos u ( 1 + e cos ν ) ( cos 2 u + cos 2 ι sin 2 u ) 1 / 2 · { sin ι cos ι sin u v turb , r + sin ι cos u v turb , ϕ + − ( cos 2 u + cos 2 ι sin 2 u ) 1 / 2 cos ι v turb , z } , $$ \begin{aligned} \frac{\mathrm{d}\iota }{\mathrm{d}t}&= -\eta _* \sqrt{\frac{a (1-e^2)}{\mu }} \frac{\cos {u}}{(1+e\cos \upnu )(\cos ^2{u} + \cos ^2\iota \sin ^2{u})^{1/2}} \cdot \nonumber \\&\quad \Bigg \{\sin \iota \cos \iota \sin {u}\,v_{\mathrm{turb}, r} + \sin \iota \cos {u} \,v_{\mathrm{turb}, \phi } + \nonumber \\&\quad - (\cos ^2{u} + \cos ^2\iota \sin ^2{u})^{1/2}\cos \iota \, v_{\mathrm{turb}, z} \Bigg \}, \end{aligned} $$(30)

d ϖ d t = η ∗ e a ( 1 − e 2 ) μ { − [ cos ν ( cos 2 ι + cos 2 u sin 2 ι ) 1 / 2 $$ \begin{aligned} \frac{\mathrm{d}\varpi }{\mathrm{d}t}&= \frac{\eta _*}{e} \sqrt{\frac{a (1-e^2)}{\mu }} \Bigg \{-\Bigg [\cos \upnu (\cos ^2\iota + \cos ^2{u} \sin ^2\iota )^{1/2} \end{aligned} $$(31)

+ ( 2 + e cos ν 1 + e cos ν ) sin 2 ι sin ν sin u cos u ( cos 2 u + cos 2 ι sin 2 u ) 1 / 2 ] v turb , r + [ ( 2 + e cos ν 1 + e cos ν ) cos ι sin ν ( cos 2 u + cos 2 ι sin 2 u ) 1 / 2 ] v turb , ϕ + [ ( 2 + e cos ν 1 + e cos ν ) cos u sin ν − cos ν sin u ] sin ι v turb , z } , $$ \begin{aligned}&\quad + \left(\frac{2+e\cos \upnu }{1+e\cos \upnu }\right)\frac{\sin ^2\iota \sin \upnu \sin {u}\cos {u}}{(\cos ^2{u} + \cos ^2\iota \sin ^2{u})^{1/2}}\Bigg ]v_{\mathrm{turb}, r} \nonumber \\&\quad + \left[\left(\frac{2+e\cos \upnu }{1+e\cos \upnu }\right)\frac{\cos \iota \sin \upnu }{(\cos ^2{u} + \cos ^2\iota \sin ^2{u})^{1/2}}\right] v_{\mathrm{turb}, \phi } \nonumber \\&\quad +\left[\left(\frac{2+e\cos \upnu }{1+e\cos \upnu }\right)\cos {u}\sin \upnu - \cos \upnu \sin {u}\right] \sin \iota \,v_{\mathrm{turb}, z} \Bigg \}, \end{aligned} $$(32)

d Ω d t = − η ∗ a ( 1 − e 2 ) μ sin u 1 + e cos ν [ cos ι sin u v turb , r ( cos 2 u + cos 2 ι sin 2 u ) 1 / 2 + cos u ( cos 2 u + cos 2 ι sin 2 u ) 1 / 2 v turb , ϕ − cos ι sin ι v turb , z ] . $$ \begin{aligned} \frac{\mathrm{d}\Omega }{\mathrm{d}t}&= -\eta _* \sqrt{\frac{a (1-e^2)}{\mu }} \frac{\sin {u}}{1+e\cos \upnu } \Bigg [ \frac{\cos \iota \sin {u} \,v_{\mathrm{turb}, r} }{(\cos ^2{u} + \cos ^2\iota \sin ^2{u})^{1/2}} \nonumber \\&\quad + \frac{\cos {u}}{(\cos ^2{u} + \cos ^2\iota \sin ^2{u})^{1/2}}v_{\mathrm{turb}, \phi } - \frac{\cos \iota }{\sin \iota }v_{\mathrm{turb}, z} \Bigg ]. \end{aligned} $$(33)

We could, in principle, linearize Eqs. (28–33) in e and ι, following the same procedure used for the deterministic system in Eqs. (12–16). However, the turbulent velocity dispersion σturb is itself of order 𝒪(e, ι), as shown in Eqs. (7) and (26). Therefore, to ensure consistency with the perturbative order of the deterministic equations, we expand Eqs. (28–33) to zeroth order in eccentricity and inclination:

d a d t = 2 a η ∗ a μ v turb , ϕ , $$ \begin{aligned} \frac{\mathrm{d}a}{\mathrm{d}t}&= 2a \, \eta _* \sqrt{\frac{a}{\mu }} v_{\mathrm{turb}, \phi },\end{aligned} $$(34)

d e d t = η ∗ a μ ( v turb , r sin ν + 2 v turb , ϕ cos ν ) , $$ \begin{aligned} \frac{\mathrm{d}e}{\mathrm{d}t}&= \eta _* \sqrt{\frac{a}{\mu }} \left(v_{\mathrm{turb}, r} \sin \upnu + 2 v_{\mathrm{turb}, \phi } \cos \upnu \right), \end{aligned} $$(35)

d ι d t = η ∗ a μ cos u v turb , z , $$ \begin{aligned} \frac{\mathrm{d}\iota }{\mathrm{d}t}&= \eta _* \sqrt{\frac{a}{\mu }} \cos {u} \, v_{\mathrm{turb}, z}, \end{aligned} $$(36)

d ϖ d t = η ∗ e a μ ( − cos ν v turb , r + 2 sin ν v turb , ϕ ) , $$ \begin{aligned} \frac{\mathrm{d}\varpi }{\mathrm{d}t}&= \frac{\eta _*}{e} \sqrt{\frac{a}{\mu }} \left(- \cos \upnu \,v_{\mathrm{turb}, r} + 2 \sin \upnu \, v_{\mathrm{turb}, \phi } \right), \end{aligned} $$(37)

d Ω d t = η ∗ i a μ sin u v turb , z . $$ \begin{aligned} \frac{\mathrm{d}\Omega }{\mathrm{d}t}&= \frac{\eta _*}{i} \sqrt{\frac{a}{\mu }} \sin {u} \,v_{\mathrm{turb}, z}. \end{aligned} $$(38)

A numerical integration of the reduced stochastic system given by Eqs. (34–38) alongside the linearized deterministic system in Eqs. (19–23) yields results that are virtually indistinguishable from those obtained using the full equations, Eqs. (12–16) and (28–33). Furthermore, the latter reduced equations remain valid across the entire Mach number range, provided that the drag coefficient η* in Eq. (8) is replaced with its general expression η from Eq. (3).

2.3. Numerical methods

In the numerical simulations discussed in the next section, we integrated the system of first-order differential equations (Eqs. 12–16) governing the evolution of the Keplerian elements using an eighth-order Dormand–Prince scheme (Hairer & Wanner 1993). The drag coefficient η* was replaced with its general form η ( v ∼ , R , z ) $ \eta(\tilde{v}, R, z) $ to account for full velocity and spatial dependence. We employed an adaptive timestep, δt, with a relative tolerance of 10−8. Turbulent velocity fluctuations were modeled as a 3D Gaussian random process with isotropic velocity dispersion σ turb = α c s $ \sigma_{\mathrm{turb}} = \sqrt{\alpha}\, {{c}_\text{s}} $. The temporal correlation of the fluctuations is characterized by the correlation time τc, as introduced above (see Terzic & Kandrup 2003; Sideris & Kandrup 2004 for further discussion). To implement this, we defined a decay factor fδt = exp(−δt/τc), such that at each timestep k + 1, each component of the turbulent velocity vector vturb was independently updated from its previous value at step k according to

v turb k + 1 = v turb k f δ t + 1 − f δ t 2 G ( σ turb ) , $$ \begin{aligned} v_{\rm turb}^{k+1} = v_{\rm turb}^k f_{\delta t}+\sqrt{1-f_{\delta t}^2}\, \mathcal{G} (\sigma _{\rm turb}), \end{aligned} $$(39)

where 𝒢(σturb) denotes a random variate drawn from a Gaussian distribution with zero mean and standard deviation σturb, generated at each timestep using the Box–Muller algorithm (see, e.g., Press et al. 2002).

For comparison, we also integrated the full equation of motion (Eq. (1)) in Cartesian coordinates for the same AGN model and stochastic prescription. In this case, we adopted a second-order modified midpoint leapfrog integrator with a fixed timestep δt = 10−4P, where P is the orbital period of mBH around MBH (see, e.g., Mikkola & Merritt 2006). After recovering the Keplerian elements from the Cartesian trajectories using standard formulae (see, e.g., Roy 2005), we compared the statistical outcomes of both approaches. For an ensemble of 2 × 103 orbits evolved over 20 unperturbed periods from identical initial conditions, we find that the spread in eccentricity and semimajor axis differs by approximately 5% and 25%, respectively. This discrepancy is attributed to the perturbative nature of the Keplerian-element equations. In contrast, the average values of a, e, and ι remain in good agreement between the two methods, with deviations consistently below 1%.

3. Simulations and results

3.1. Evolution of the orbital elements

Figure 3 shows the evolution of the orbital elements of a 20 M⊙ BH embedded in the AGN disk, comparing cases with and without turbulent forcing. The BH is initially placed at the migration trap radius, R ∼ 104 Rg ∼ 103 au, where the net migration torque vanishes. It begins with an eccentricity of e = 0.1 and an inclination of ι = 0.5°, ensuring that the orbit remains within a few disk scale heights. At this radius the disk has an aspect ratio of H/R = 8.5 × 10−3, corresponding to an inclination of 0.487°. The BH has an initial period of P = 9.1 yr.

thumbnail Fig. 3.

Evolution of the semimajor axis (top), the eccentricity (middle), and the inclination (bottom) for an inclined (ι0 = 0.5°) and eccentric (e0 = 0.1) BH (mBH = 20 M⊙) undergoing dynamical friction in the AGN disk. The thick blue line indicates the evolution without the turbulent velocity field, while each thin red line includes a different realization of the turbulent velocity field. The BH is placed at the migration trap within an AGN disk around a 107 M⊙ SMBH. The dotted black line is the mean evolution of the realizations including turbulence. In the top panel, the dot-dashed lines indicate the mean square change in semimajor axis ⟨Δa2⟩ (Eq. (44)). In the middle and bottom panels, the dot-dashed green line is the ratio between the turbulent velocity dispersion σturb and the circular velocity vcirc. The vertical dashed lines indicate the first three midplane crossings. The turbulence prevents the full circularization and alignment of the embedded BH.

In the absence of turbulence (blue curve), the BH rapidly circularizes and aligns with the disk midplane. During the first two orbital periods, dynamical friction is enhanced near the midplane crossings (vertical lines) due to the steep vertical density gradient. At each crossing, the semimajor axis increases when the BH passes through the disk from above and decreases when crossing from below. This asymmetry arises because the first crossing occurs near apocenter (corresponding to the negative cosu in Eq. (19)) where the BH’s orbital velocity is lower than the local circular velocity, allowing it to gain energy from the gas flow. After approximately three crossings (i.e., ∼1.5 orbital periods), the inclination drops below ∼0.18° (zmax/H ∼ 0.4), such that the vertical variations in gas density become negligible. From this point on, the semimajor axis stabilizes at the final value of afin ≃ 975 au ≃ 9.9 × 104 Rg, while both eccentricity and inclination decay exponentially, consistent with the analytic expectations from Eqs. (20–21) and the results of hydrodynamical simulations (Rowan et al. 2025b; Whitehead et al. 2025b).

Notably, the inclination decays at an exponential rate of ∼η*/2.6, which is slightly slower than the rate ∼η*/2 obtained by averaging the factor cos2u in Eq. (21) over the argument of latitude u. This discrepancy likely arises because the inclination decays too rapidly for u to circulate during the alignment phase. Given the eccentricity damping timescale τ* = 1/η* ≃ 0.75 yr ≃ 0.08P, the orbit settles into the disk while u librates around u* = 5.38 such that 1/cos2u* ≈ 2.6, making the decay rate sensitive to the librating argument of latitude rather than its time-averaged behavior.

In the joint limit of low eccentricity and inclination, the fixed point u* can be estimated by setting du/dt = 0, yielding

n [ 1 + 2 e cos ( u − ω 0 ) ] = 1 4 η ∗ [ cos ( u − ω 0 ) − 2 ] sin ( 2 u ) , $$ \begin{aligned} n \left[1+2e\cos {(u-\omega _0)}\right] = \frac{1}{4} \eta _* \left[\cos (u-\omega _0) - 2\right] \sin {(2u)}, \end{aligned} $$(40)

where n = μ / a 3 $ n=\sqrt{{\mu}/{a^3}} $ is the mean motion and ω0 is the argument of pericenter at the onset of the exponential decay. This transcendental equation admits real roots only in the fast damping regime η* ≫ n, i.e., when the damping timescale is shorter than the orbital period. As n/η* → 0, the solutions collapse onto the four zeros of sin(2u) = 0. However, only the zeros at 0 and π are physical because the inclination damping in Eq. (21) would otherwise vanish. Setting e = 0, a first-order expansion in n/η* about each root then gives

u ∗ ( k ) = k π + n η ∗ 2 cos ( k π − ω 0 ) − 2 , k ∈ { 0 , 1 } . $$ \begin{aligned} u_*^{(k)} = k \pi + \frac{n}{\eta _*} \frac{2}{\cos (k\pi - \omega _0) -2},\quad k\in \{0,1\}. \end{aligned} $$(41)

We obtain u*(3) = 5.48, in good agreement with the numerically determined libration of the argument of latitude. This level of agreement is notable given the relatively large value of the ratio n/η* ≃ 0.48, which pushes the limits of the fast damping approximation used in the derivation.

In the turbulent case (red curves), the evolution of the BH becomes stochastic once it settles into the disk at ∼1.5 orbital periods. The semimajor axis begins to diffuse around its final unperturbed value afin, as predicted by Eq. (34). Both eccentricity and inclination deviate from the deterministic case, exhibiting sustained fluctuations above zero. Over time, they reach a steady state in which stochastic forcing balances dynamical damping. When ensemble-averaged (dotted black line), e and ι converge to equilibrium values of (e, ι)eq ≃ (1.0 × 10−3, 1.3 × 10−3). This behavior supports the statistical description we introduced in Sect. 3.3, where the interplay between linear damping and stochastic driving yields steady-state Rayleigh distributions for both e and ι.

The dash-dotted green lines in the lower panels represent the order-of-magnitude theoretical equilibrium root mean square values of e and ι based on the ratio between the turbulent velocity dispersion and the local gas circular velocity, σturb/vcirc ≃ 8.6 × 10−4 (see Eq. (26)). The good agreement with the mean values supports the idea that the final state reflects a statistical balance between stochastic driving and dissipative damping.

3.2. Steady-state distributions of e and ι

The turbulent velocity field hampers the full orbital circularization and alignment of the embedded objects. This occurs because dynamical friction tends to bring the BH at rest with respect to the local velocity field, which is dominated by large-scale turbulent flows. Consequently, the BH’s orbit does not become fully circular or aligned with the midplane but retains a residual eccentricity and inclination. On average, the eccentricity and inclination settle onto a plateau, whose characteristic magnitude is determined by the ratio of the turbulent velocity dispersion σt and the local circular velocity vcirc. This is shown in Fig. 4, which depicts the distribution of the eccentricity and inclination of the embedded BHs after 20 orbital periods, roughly corresponding to 244 τ*.

thumbnail Fig. 4.

Distributions (normalized to one) of the semimajor axis (top), the eccentricity (middle) and the inclination (bottom) for 2.5 × 104 realizations of an embedded BH, evaluated after 20 initial orbital periods. The initial conditions are identical to those in Fig. 3. In the top panel, the blue line marks the final semimajor axis from the nonturbulent simulation, while the green curve shows a Gaussian distribution with standard deviation given by Eq. (44). In the middle and bottom panels, the dot-dashed green line is the ratio between the turbulent velocity dispersion σturb and the circular velocity vcirc, while the green curves indicate Rayleigh distributions with mean value σturb/vcirc. The insets show the same distributions in logarithm scale, highlighting the presence of fatter tails with respect to the Rayleigh prediction.

The semimajor axis undergoes a random walk due to the stochastic forcing, whose mean square amplitude can be estimated from statistical considerations. Given the stochastic derivative in Eq. (34), the mean squared change ⟨Δa2⟩ is

〈 d a d t , d a d t 〉 = ( 2 a η ∗ v circ ) 2 ∫ 0 t ∫ 0 t v turb ( t 1 ) v turb ( t 2 ) d t 1 d t 2 , $$ \begin{aligned} \left\langle \frac{\mathrm{d}a}{\mathrm{d}t}, \frac{\mathrm{d}a}{\mathrm{d}t} \right\rangle&= \left( \frac{2 a\eta _*}{v_\mathrm{circ} } \right)^2 \int _0^t\int _0^t v_\mathrm{turb} (t_1) v_\mathrm{turb} (t_2) \mathrm{d}t_1 \mathrm{d}t_2,\end{aligned} $$(42)

= ( 2 a η ∗ σ turb v circ ) 2 ∫ 0 t ∫ 0 t exp ( − | t 1 − t 2 | τ c ) d t 1 d t 2 , $$ \begin{aligned}&= \left( \frac{2 a\eta _* \sigma _{\rm turb}}{v_\mathrm{circ} } \right)^2 \int _0^t\int _0^t \exp {\left(-\frac{|t_1 - t_2|}{\tau _{\rm c}}\right)}\mathrm{d}t_1 \mathrm{d}t_2, \end{aligned} $$(43)

which for t ≫ τc reduces to

⟨ Δ a 2 ⟩ ( t ) = 8 η ∗ 2 a 2 σ turb 2 v circ 2 t , $$ \begin{aligned} \langle \Delta a^2 \rangle (t) = 8 \eta _*^2 a^2 \frac{\sigma _\mathrm{turb} ^2}{v_{\rm circ}^2}\, t, \end{aligned} $$(44)

where we used vcirc2 = μ/a and τc2 = a3/μ. We plot Eq. (44) as a function of time in Fig. 3 and show the corresponding normal distribution in Fig. 4. Our analytic estimate matches the numerical results very well, despite the many assumptions, with small deviations likely due to neglecting the deterministic derivative in Eq. (12). In fact, even though any pure deterministic term would average to zero in Eq. (42), Eq. (12) depends on the eccentricity, whose mean square value is nonzero, albeit small.

Contrary to what happens for the semimajor axis, the distributions of e and ι settle to a steady state, reasonably well described in terms of Rayleigh distributions. In the middle and bottom panel of Fig. 4 we compare the distributions obtained from the numerical integrations (red histograms) with the matching Rayleigh distribution (solid green lines) with a choice of scale parameter so that the mean value is equal to σturb/vcirc. For both quantities, the empirical distribution and its semi-analytical estimate differ by less than the 15% over the significant range of e and ι. We note that, for the distribution of eccentricity, an even better match is obtained by setting the mean value to 1.2σturb/vcirc.

3.3. Stochastic equilibrium of e and ι

The equations for de/dt and dι/dt exhibit similar dynamics. The deterministic terms in Eqs. (20–21) include linear damping components proportional to −η*e and −η*ι, while the stochastic terms in Eqs. (35–36) account for turbulent forcing, proportional to vturb. The deterministic terms drive e and ι toward zero, while the stochastic terms introduce diffusion. We stress the fact that the above equations are valid in the subsonic regime where the friction coefficient reduces to η* (see Eq. (8)). Over long timescales, the probability distributions for the eccentricity, 𝒫(e), and inclination, 𝒫(ι), reach steady-state distributions dictated by the balance between frictional damping and stochastic forcing. These steady-state distributions can be derived analytically by exploiting the fact that Eqs. (35–38) and (20–23) are closely related to the general stochastic differential equation for a Ornstein–Uhlenbeck process Xt (hereafter OU; see, e.g., Risken 1989; Gardiner 1994) that reads

d X t = − θ X t d t + σ d W t . $$ \begin{aligned} \mathrm{d}X_t = - \theta X_t \mathrm{d}t + \sigma \mathrm{d}W_t. \end{aligned} $$(45)

In the equation above, θ is the damping rate, σ is the noise intensity, and Wt represents the underlying Wiener process. In our case η* and σturb correspond to θ and σ, respectively. At variance with the standard OU process, here we consider a Markovian stochastic forcing with an exponentially decaying autocorrelation function rather than a delta function.

We can explicitly recast the eccentricity evolution equations into the form of an OU process. While the equation for de/dt might appear OU-like, it is coupled to the evolution of the longitude of pericenter dϖ/dt. To isolate the true linear stochastic equation, we switch variables to the components of the eccentricity vector:

g x = e cos ϖ , $$ \begin{aligned} g_x&= e \cos {\varpi },\end{aligned} $$(46)

g y = e sin ϖ . $$ \begin{aligned} g_y&= e \sin {\varpi }. \end{aligned} $$(47)

After algebraic manipulation, we obtain

d g x d t = − η ∗ g x + η ∗ a μ ( 2 v turb , ϕ cos ℓ + v turb , r sin ℓ ) , $$ \begin{aligned} \frac{\mathrm{d}g_x}{\mathrm{d}t}&= - \eta _* g_x+ \eta _* \sqrt{\frac{a}{\mu }} \left( 2v_{\mathrm{turb}, \phi } \cos \ell + v_{\mathrm{turb}, r} \sin \ell \right), \end{aligned} $$(48)

d g y d t = − η ∗ g y + η ∗ a μ ( 2 v turb , ϕ sin ℓ − v turb , r cos ℓ ) , $$ \begin{aligned} \frac{\mathrm{d}g_y}{\mathrm{d}t}&= - \eta _* g_y + \eta _* \sqrt{\frac{a}{\mu }} \left( 2v_{\mathrm{turb}, \phi } \sin \ell - v_{\mathrm{turb}, r} \cos \ell \right), \end{aligned} $$(49)

which describe a 2D OU-like process for the eccentricity vector. Here ℓ = ϖ + ν = Ω + ω + ν represents the true longitude. In the two equations above, the first term represents the deterministic linear damping θ, and the second term is the stochastic forcing σ due to the radial and azimuthal components of the velocity.

A 1D OU process has a steady-state distribution that is a Gaussian with zero mean. The variance of this distribution can be derived using Itô’s calculus (e.g., see Risken 1989; Gardiner 1994). We first expressed the solution of Eq. (48) as Itô’s stochastic integral:

g x ( t ) = D ∫ − ∞ t exp ( − η ∗ ( t − s ) ) Σ x ( s ) d s , $$ \begin{aligned} g_x(t) = D \int ^{t}_{-\infty } \exp {\Big (-\eta _* (t-s)\Big )} \,\Sigma _x(s) \,\mathrm{d}s, \end{aligned} $$(50)

where D = η*/vcirc and Σx(t) = 2vturb, ϕcosu(t) + vturb, r sin u(t). The variance of gx is then

⟨ g x 2 ⟩ = D 2 ∫ 0 ∞ ∫ 0 ∞ exp ( − η ∗ ( t 1 + t 2 ) ) · 〈 Σ x ( t − t 1 ) , Σ x ( t − t 2 ) 〉 d t 1 d t 2 = 2 η ∗ σ turb 2 v circ 2 ∫ 0 ∞ exp ( − η ∗ t ) S ( t ) d t = 2 η ∗ σ turb 2 v circ 2 A ( η ∗ + Ω circ ) − B n ( η ∗ + Ω circ ) 2 + n 2 , $$ \begin{aligned} \langle g_x^2 \rangle&= D^2 \int _0^\infty \int _0^\infty \exp \Big ( -\eta _*(t_1 + t_2) \Big ) \cdot \nonumber \\&\qquad \qquad \qquad \left\langle \Sigma _x(t - t_1), \Sigma _x(t - t_2) \right\rangle \, \mathrm{d}t_1 \, \mathrm{d}t_2 \nonumber \\&= \frac{2 \eta _* \sigma ^2_{\rm turb}}{v^2_{\rm circ}} \int _0^\infty \exp (-\eta _* t) S(t) \, \mathrm{d}t \nonumber \\&= \frac{2 \eta _* \sigma ^2_{\rm turb}}{v^2_{\rm circ}} \frac{A \,(\eta _* + \Omega _{\rm circ}) - B \,n}{(\eta _* + \Omega _{\rm circ})^2 + n^2}, \end{aligned} $$(51)

where

S ( t ) = exp ( − t / τ c ) [ A cos ( n t ) − B sin ( n t ) ] , $$ \begin{aligned} S(t)&= \exp {(-t/\tau _c)} \left[A \cos (nt) - B \sin (nt) \right],\end{aligned} $$(52)

A = 4 cos 2 ℓ + sin 2 ℓ , $$ \begin{aligned} A&= 4\cos ^2\ell + \sin ^2\ell , \end{aligned} $$(53)

B = 3 cos ℓ sin ℓ . $$ \begin{aligned} B&= 3\cos \ell \sin \ell . \end{aligned} $$(54)

We distinguish τc = 1/Ωcirc, the correlation time of turbulent forcing (set by the gas), from the BH’s mean motion n. Although they are equal in our model, this separation allows us to keep the generality for more complex models.

The A and B terms depend on the true longitude and need to be estimated depending on the ratio between the damping rate η* and the orbital and gas correlation frequencies n and Ωcirc. Figure 5 shows how the damping timescale τ* compares to the orbital period in the AGN disk.

thumbnail Fig. 5.

Ratio between the damping timescale τ* (Eq. (9)) and the orbital period P, as a function of the distance from the central SMBH. The estimate assumes a BH with mass mBH = 20 M⊙ based on the AGN disk profile for SMBHs with masses MSMBH = 106, 3 × 107, 7 × 106, and 2 × 107 M⊙, with α = 0.1. The vertical lines indicate the location of the migration traps.

In the region of the migration trap at R ∼ 104 Rg, we obtain n/η* ≃ 0.48. In this case, the damping acts moderately faster than the orbital timescale, and we consider the true longitude ℓ to be constant over one damping timescale. However, since ℓ continues to circulate and is not dynamically driven to a fixed point (unlike the argument of latitude u, which librates around u*, see Eq. (41)), we do not expect any preferred value of ℓ during the damping. We can thus ensemble average A over ℓ, yielding ⟨A⟩ = 5/2. Using σe2 = 2⟨gx2⟩. This holds due to the statistical symmetry between gx and gy. The steady-state variance of the eccentricity in the fast damping regime reads:

σ e 2 ≃ 10 ( σ turb v circ ) 2 (fast damping) . $$ \begin{aligned} \sigma ^2_e \simeq 10 \left(\frac{\sigma _{\rm turb}}{v_{\rm circ}}\right)^2 \quad \text{(fast} \text{ damping)}. \end{aligned} $$(55)

In other regions of the AGN disk, or for smaller BHs, the system can enter the slow damping regime where η* ≪ n. In this case, the true longitude circulates rapidly compared to the damping timescale, and we can average over its values. As a result, the stochastic forcing averages over all phases, leading to the orbit-averaged values of ⟨A⟩ = 5/2 and ⟨B⟩ = 0. The resulting eccentricity variance in the slow damping regime is

σ e 2 ≃ 5 η ∗ n ( σ turb v circ ) 2 (slow damping) , $$ \begin{aligned} \sigma ^2_e \simeq 5\,\frac{\eta _*}{n} \left(\frac{\sigma _{\rm turb}}{v_{\rm circ}}\right)^2 \quad \text{(slow} \text{ damping)}, \end{aligned} $$(56)

where we set Ωcirc = n. This expression shows that the stochastic excitation is suppressed by the small parameter η*/n ≪ 1, reflecting the fact that the orbital motion rapidly de-correlates each turbulent kick before it can accumulate significant eccentricity.

The variance of the inclination, σι2, can be derived by applying the same OU process analysis to the inclination vector components obtained via the change of variables:

h x = i cos Ω , $$ \begin{aligned} h_x&= i \cos {\Omega },\end{aligned} $$(57)

h y = i sin Ω . $$ \begin{aligned} h_y&= i \sin {\Omega }. \end{aligned} $$(58)

In this formulation, hx and hy follow stochastic differential equations analogous to those for eccentricity. The resulting steady-state variance of ι takes the form:

σ ι 2 = { 2 3 ( σ turb v circ ) 2 , if η ∗ ≫ n (fast damping) , [ 2.5 e x ] 2 η ∗ n ( σ turb v circ ) 2 , if η ∗ ≪ n (slow damping) . $$ \begin{aligned} \sigma ^2_\iota = {\left\{ \begin{array}{ll} \dfrac{2}{3} \left(\dfrac{\sigma _{\rm turb}}{v_{\rm circ}}\right)^2,&\text{ if } \eta _* \gg n \quad \text{(fast} \text{ damping)},\\ [2.5ex] 2 \dfrac{\eta _*}{n} \left(\dfrac{\sigma _{\rm turb}}{v_{\rm circ}}\right)^2,&\text{ if } \eta _* \ll n \quad \text{(slow} \text{ damping)}. \end{array}\right.} \end{aligned} $$(59)

These analytical results are consistent with the numerical simulations shown in Fig. 4, and they account for the observed difference between σι and σe. In particular, they explain why the steady-state variance of eccentricity exceeds that of inclination. For the simulations considered here, the ratio n/η* ≃ 0.48 places the system between the fast and slow damping regimes, leading to intermediate values of σι and σe.

4. Summary and conclusions

Using semi-analytical stochastic methods and numerical simulations, we investigated the drag exerted by a turbulent AGN disk surrounding an SMBH on a stellar-mass BH initially located at the migration trap. We introduced a numerical technique, inspired by analogous studies in protoplanetary disks, that incorporates turbulence, modeled as a stochastic velocity field, without solving the full hydrodynamical equations. Instead, the method relies on effective stochastic differential equations.

We implemented a full solver in Cartesian coordinates and a reduced solver based on Keplerian orbital elements, and find excellent agreement between the two approaches across a range of initial conditions. In the subsonic regime (ℳ ≪ 1), we derived approximate Langevin equations governing the orbital element evolution, enabling us to compute the stationary distributions of eccentricity e and inclination ι. We compared the resulting distributions from N = 2.5 × 104 independent realizations to analytical Rayleigh probability density functions, finding good agreement.

We analytically estimated the variance of the steady-state Rayleigh distributions using the linearized stochastic equations and find good agreement with the results of the full numerical simulations. By expressing the relevant quantities in terms of local disk properties, we obtain the following expressions:

σ e 2 = 15 σ ι 2 = 10 α ( H R ) 2 , (fast damping, η ∗ ≫ n ) , σ e 2 = 5 2 σ ι 2 = 5 α η ∗ n ( H R ) 2 , (slow damping, η ∗ ≪ n ) . $$ \begin{aligned}&\sigma _e^2 = 15\, \sigma _\iota ^2 = 10\, \alpha \left(\frac{H}{R}\right)^2, &\quad (\text{fast} \text{ damping,} \eta _* \gg n),\\ &\sigma _e^2 = \frac{5}{2}\, \sigma _\iota ^2 = 5\, \alpha \frac{\eta _*}{n} \left(\frac{H}{R}\right)^2, &\quad (\text{slow} \text{ damping,} \eta _* \ll n). \end{aligned} $$(60)

Here, α is the alpha-viscosity parameter characterizing the strength of turbulence, H is the disk scale height, n = Ω circ = G M SMBH / R 3 $ n = \Omega_{\mathrm{circ}} = \sqrt{G M_{\mathrm{SMBH}}/R^3} $ is the circular orbital frequency, and η* is the subsonic damping rate determined by the local gas properties (see Eq. (8)).

Equation (60) can be readily used to sample the eccentricity and inclination of BHs embedded in AGN disks, both in Monte-Carlo simulations (McKernan et al. 2020, 2025; Tagawa et al. 2021; Rowan et al. 2024a; Cook et al. 2024; Delfavero et al. 2024) and as initial conditions for BH scattering hydrodynamical simulations (Whitehead et al. 2025a; Rowan et al. 2025a), which to date have assumed perfectly circular and aligned orbits. Notably, Trani et al. (2024) show that the initial eccentricity and inclination of BHs in a disk configuration strongly affect the merger rate from three-body encounters. Specifically, they find that the merger rate in a dispersion-dominated disk is suppressed by a factor of ∼42 relative to a shear-dominated disk (see their Figure 5). A similar suppression is expected for the formation of binaries via Jacobi captures mediated by dynamical friction, since even modest values of e and ι increase the relative velocity between interacting BHs, significantly raising the amount of energy that must be dissipated for capture to occur (see Dodici & Tremaine 2024).

In our current formulation and numerical implementation, we neglected migration torques arising from the perturbation of the disk’s density by the orbiting BH. As a result, our simulations and theoretical predictions likely overestimate the diffusion in semimajor axis. In reality, such torques would act as a restoring force, counteracting orbital changes from the gas drag as well as the stochastic fluctuations in a. Consequently, the BH would remain near the migration trap, and the semimajor axis distribution would settle into a steady-state centered on the trap location. Nonetheless, our simulations show that the radial wandering remains sufficiently small to justify the assumption of constant a used in the derivation of Eq. (60).

We note that simplified hydrodynamical simulations by Wu et al. (2024) suggest that turbulence may also reduce the magnitude of the migration torque itself, thereby diminishing its ability to stabilize the semimajor axis (see also Baruteau & Lin 2010; Pierens et al. 2012; Stoll et al. 2017). We stress that the framework presented here relies on the linear theory of dynamical friction, which neglects nonlinear interactions between the BH’s wake and the turbulent gas. Furthermore, our treatment neglected the back-reaction of the embedded BH on the surrounding gas. In reality, accretion onto the BH can heat the local disk, drive local turbulence, or create a low-density cavity, potentially reducing the drag and altering the migration trap structure (e.g., Gilbaum & Stone 2022; Epstein-Martin et al. 2025; Chen & Dai 2025). Quantifying these effects lie beyond the scope of this work and should be investigated using hydrodynamical simulations that self-consistently model both turbulence and gas–BH coupling.

1

Note that, in the limit of ℳ ≫ 1, one recovers the classical expression for the Chandrasekhar (1943) dynamical friction in a collisionless gravitational system of mean mass density ρ.

Acknowledgments

A.A.T. would like to thank Martin Pessah and Troels Haugbølle for helpful discussions, which provided valuable directions for this research. A.A.T. further acknowledges the hospitality of Firenze University’s physics & astronomy department at Arcetri, where this work was started, and acknowledges support from the Horizon Europe research and innovation programs under the Marie Skłodowska-Curie grant agreement no. 101103134. P.F.D.C. wishes to acknowledge funding by “Fondazione Cassa di Risparmio di Firenze” under the project HIPERCRHEL for the use of high performance computing resources at the University of Firenze. We thank the referee for their insightful comments, which helped improve the clarity of this work.

References

  1. Abbott, R., Abbott, T. D., Acernese, F., et al. 2023, Phys. Rev. X, 13, 041039 [Google Scholar]
  2. Amaro-Seoane, P., Gair, J. R., Freitag, M., et al. 2007, Classical Quantum Gravity, 24, R113 [NASA ADS] [CrossRef] [Google Scholar]
  3. Amaro-Seoane, P., Audley, H., Babak, S., et al. 2017, ArXiv e-prints [arXiv:1702.00786] [Google Scholar]
  4. Armitage, P. J. 2011, ARA&A, 49, 195 [Google Scholar]
  5. Babak, S., Gair, J., Sesana, A., et al. 2017, Phys. Rev. D, 95, 103012 [NASA ADS] [CrossRef] [Google Scholar]
  6. Balbus, S. A., & Hawley, J. F. 1998, Rev. Mod. Phys., 70, 1 [Google Scholar]
  7. Baruteau, C., & Lin, D. N. C. 2010, ApJ, 709, 759 [NASA ADS] [CrossRef] [Google Scholar]
  8. Beutler, G. 2005, Methods of Celestial Mechanics. Vol. I: Physical, Mathematical, and Numerical Principles [Google Scholar]
  9. Branchesi, M., Maggiore, M., Alonso, D., et al. 2023, JCAP, 2023, 068 [CrossRef] [Google Scholar]
  10. Burns, J. A. 1976, Am. J. Phys., 44, 944 [NASA ADS] [CrossRef] [Google Scholar]
  11. Chandrasekhar, S. 1943, ApJ, 97, 255 [Google Scholar]
  12. Chen, K., & Dai, Z.-G. 2024, ApJ, 961, 206 [Google Scholar]
  13. Chen, K., & Dai, Z.-G. 2025, ApJ, 987, 214 [Google Scholar]
  14. Cook, H. E., McKernan, B., Ford, K. E. S., et al. 2024, ApJ, submitted [arXiv:2411.10590] [Google Scholar]
  15. DeLaurentiis, S., Epstein-Martin, M., & Haiman, Z. 2023, MNRAS, 523, 1126 [NASA ADS] [CrossRef] [Google Scholar]
  16. Delfavero, V., Ford, K. E. S., McKernan, B., et al. 2024, ApJ, 989, 67 [Google Scholar]
  17. Dempsey, A. M., Li, H., Mishra, B., & Li, S. 2022, ApJ, 940, 155 [Google Scholar]
  18. Dittmann, A. J., Dempsey, A. M., & Li, H. 2024, ApJ, 964, 61 [Google Scholar]
  19. Dittmann, A. J., Dempsey, A. M., & Li, H. 2025, ApJ, 990, 137 [Google Scholar]
  20. Dodici, M., & Tremaine, S. 2024, ApJ, 972, 193 [Google Scholar]
  21. Dosopoulou, F. 2024, Phys. Rev. D, 110, 083027 [Google Scholar]
  22. Epstein-Martin, M., Tagawa, H., Haiman, Z., & Perna, R. 2025, MNRAS, 537, 3396 [Google Scholar]
  23. Fischer, M. S., & Sagunski, L. 2024, A&A, 690, A299 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Gangardt, D., Trani, A. A., Bonnerot, C., & Gerosa, D. 2024, MNRAS, 530, 3689 [Google Scholar]
  25. Gardiner, C. W. 1994, Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences (Berlin: Springer) [Google Scholar]
  26. Gayathri, V., Wysocki, D., Yang, Y., et al. 2023, ApJ, 945, L29 [Google Scholar]
  27. Gilbaum, S., & Stone, N. C. 2022, ApJ, 928, 191 [NASA ADS] [CrossRef] [Google Scholar]
  28. Graham, M. J., Ford, K. E. S., McKernan, B., et al. 2020, Phys. Rev. Lett., 124, 251102 [NASA ADS] [CrossRef] [Google Scholar]
  29. Grishin, E., Gilbaum, S., & Stone, N. C. 2024, MNRAS, 530, 2114 [NASA ADS] [CrossRef] [Google Scholar]
  30. Guilera, O. M., Miller Bertolami, M. M., Masset, F., et al. 2021, MNRAS, 507, 3638 [NASA ADS] [CrossRef] [Google Scholar]
  31. Hairer, E., & Wanner, G. 1993, Runge-Kutta and Extrapolation Methods (Berlin, Heidelberg: Springer), 129 [Google Scholar]
  32. Hannuksela, O. A., Ng, K. C. Y., & Li, T. G. F. 2020, Phys. Rev. D, 102, 103022 [Google Scholar]
  33. Hopman, C., & Alexander, T. 2005, ApJ, 629, 362 [NASA ADS] [CrossRef] [Google Scholar]
  34. Hu, W.-R., & Wu, Y.-L. 2017, Natl. Sci. Rev., 4, 685 [CrossRef] [Google Scholar]
  35. Ishibashi, W., & Gröbner, M. 2020, A&A, 639, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. Janiuk, A., Czerny, B., Siemiginowska, A., & Szczerba, R. 2004, ApJ, 602, 595 [NASA ADS] [CrossRef] [Google Scholar]
  37. Johnson, E. T., Goodman, J., & Menou, K. 2006, ApJ, 647, 1413, publisher: IOP Publishing [Google Scholar]
  38. Kavanagh, B. J., Karydas, T. K., Bertone, G., Di Cintio, P., & Pasquato, M. 2025, Phys. Rev. D, 111, 063071 [Google Scholar]
  39. Li, R., & Lai, D. 2022, MNRAS, 517, 1602 [Google Scholar]
  40. Li, R., & Lai, D. 2023, MNRAS, 522, 1881 [Google Scholar]
  41. Li, R., & Lai, D. 2024, MNRAS, 529, 348 [Google Scholar]
  42. Li, Y.-P., Dempsey, A. M., Li, S., Li, H., & Li, J. 2021, ApJ, 911, 124 [NASA ADS] [CrossRef] [Google Scholar]
  43. Li, J., Dempsey, A. M., Li, H., Lai, D., & Li, S. 2023, ApJ, 944, L42 [NASA ADS] [CrossRef] [Google Scholar]
  44. Maggiore, M., Van Den Broeck, C., Bartolo, N., et al. 2020, JCAP, 03, 050 [CrossRef] [Google Scholar]
  45. Mandel, I., Brown, D. A., Gair, J. R., & Miller, M. C. 2008, ApJ, 681, 1431 [Google Scholar]
  46. Mapelli, M. 2021, Handbook of Gravitational Wave Astronomy, 4 [Google Scholar]
  47. Masset, F. S. 2017, MNRAS, 472, 4204 [NASA ADS] [CrossRef] [Google Scholar]
  48. McKernan, B., Ford, K. E. S., O’Shaugnessy, R., & Wysocki, D. 2020, MNRAS, 494, 1203 [NASA ADS] [CrossRef] [Google Scholar]
  49. McKernan, B., Ford, K. E. S., Cook, H. E., et al. 2025, ApJ, 990, 217 [Google Scholar]
  50. Mikkola, S., & Merritt, D. 2006, MNRAS, 372, 219 [NASA ADS] [CrossRef] [Google Scholar]
  51. Mishra, B., & Calcino, J. 2024, ArXiv e-prints [arXiv:2409.05614] [Google Scholar]
  52. Montalvo, D., Smith-Orlik, A., Rastgoo, S., et al. 2024, Universe, 10, 427 [Google Scholar]
  53. Mukherjee, D., Holgado, A. M., Ogiya, G., & Trac, H. 2024, MNRAS, 533, 2335 [NASA ADS] [CrossRef] [Google Scholar]
  54. Nelson, R. P. 2005, A&A, 443, 1067 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Nelson, R. P., & Gressel, O. 2010, MNRAS, 409, 639 [Google Scholar]
  56. Nelson, R. P., & Papaloizou, J. C. B. 2004, MNRAS, 350, 849 [NASA ADS] [CrossRef] [Google Scholar]
  57. Oishi, J. S., Mac Low, M.-M., & Menou, K. 2007, ApJ, 670, 805 [Google Scholar]
  58. Ostriker, E. C. 1999, ApJ, 513, 252 [Google Scholar]
  59. Papaloizou, J. C. B., Nelson, R. P., & Snellgrove, M. D. 2004, MNRAS, 350, 829, publisher: OUP ADS Bibcode: 2004MNRAS.350.829P [NASA ADS] [CrossRef] [Google Scholar]
  60. Picogna, G., Stoll, M. H. R., & Kley, W. 2018, A&A, 616, A116 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  61. Pierens, A., Baruteau, C., & Hersant, F. 2012, MNRAS, 427, 1562 [Google Scholar]
  62. Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 2002, Numerical Recipes in C++ : The Art of Scientific Computing [Google Scholar]
  63. Punturo, M., Abernathy, M., Acernese, F., et al. 2010, Class. Quant. Grav., 27, 194002 [Google Scholar]
  64. Rein, H., & Papaloizou, J. C. B. 2009, A&A, 497, 595 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  65. Reitze, D., Adhikari, R. X., Ballmer, S., et al. 2019, Bull. Am. Astron. Soc., 51, 035 [Google Scholar]
  66. Ren, J., Chen, K., Wang, Y., & Dai, Z.-G. 2022, ApJ, 940, L44 [Google Scholar]
  67. Risken, H. 1989, The Fokker-Planck Equation. Methods of Solution and Applications [Google Scholar]
  68. Rowan, C., Boekholt, T., Kocsis, B., & Haiman, Z. 2023, MNRAS, 524, 2770 [NASA ADS] [CrossRef] [Google Scholar]
  69. Rowan, C., Whitehead, H., & Kocsis, B. 2024a, ArXiv e-prints [arXiv:2412.12086] [Google Scholar]
  70. Rowan, C., Whitehead, H., Boekholt, T., Kocsis, B., & Haiman, Z. 2024b, MNRAS, 527, 10448 [Google Scholar]
  71. Rowan, C., Whitehead, H., Fabj, G., et al. 2025a, MNRAS, 539, 1501 [Google Scholar]
  72. Rowan, C., Whitehead, H., Fabj, G., et al. 2025b, MNRAS, 543, 132 [Google Scholar]
  73. Roy, A. E. 2005, Orbital Motion (Bristol, UK: Institute of Physics Publishing) [Google Scholar]
  74. Samsing, J., Bartos, I., D’Orazio, D. J., et al. 2022, Nature, 603, 237 [NASA ADS] [CrossRef] [Google Scholar]
  75. Sartorello, S., Di Cintio, P., Trani, A. A., & Pasquato, M. 2025, A&A, 698, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  76. Secunda, A., Bellovary, J., Mac Low, M.-M., et al. 2019, ApJ, 878, 85 [NASA ADS] [CrossRef] [Google Scholar]
  77. Sideris, I. V., & Kandrup, H. E. 2004, ApJ, 602, 678 [Google Scholar]
  78. Sirko, E., & Goodman, J. 2003, MNRAS, 341, 501 [NASA ADS] [CrossRef] [Google Scholar]
  79. Spera, M., Trani, A. A., & Mencagli, M. 2022, Galaxies, 10, 76 [NASA ADS] [CrossRef] [Google Scholar]
  80. Stoll, M. H. R., Picogna, G., & Kley, W. 2017, A&A, 604, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  81. Stone, N. C., Metzger, B. D., & Haiman, Z. 2016, MNRAS, 464, 946 [Google Scholar]
  82. Tagawa, H., Kocsis, B., Haiman, Z., et al. 2021, ApJ, 907, L20 [NASA ADS] [CrossRef] [Google Scholar]
  83. Tagawa, H., Kimura, S. S., Haiman, Z., Perna, R., & Bartos, I. 2023, ApJ, 950, 13 [NASA ADS] [CrossRef] [Google Scholar]
  84. Tagawa, H., Kimura, S. S., Haiman, Z., Perna, R., & Bartos, I. 2024, ApJ, 966, 21 [Google Scholar]
  85. Takátsy, J., Zwick, L., Hendriks, K., et al. 2025, ArXiv e-prints [arXiv:2505.09513] [Google Scholar]
  86. Terzic, B., & Kandrup, H. E. 2003, Phys. Rev. E, submitted [arXiv:astro-ph/0312434] [Google Scholar]
  87. The LIGO Scientific Collaboration, the Virgo Collaboration, the KAGRA Collaboration, et al. 2023, Phys. Rev. X, 13, 011048 [NASA ADS] [Google Scholar]
  88. Trani, A. A., Quaini, S., & Colpi, M. 2024, A&A, 683, A135 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  89. Whitehead, H., Rowan, C., Boekholt, T., & Kocsis, B. 2024a, MNRAS, 531, 4656 [Google Scholar]
  90. Whitehead, H., Rowan, C., Boekholt, T., & Kocsis, B. 2024b, MNRAS, 533, 1766 [Google Scholar]
  91. Whitehead, H., Rowan, C., & Kocsis, B. 2025a, MNRAS, 542, 1033 [Google Scholar]
  92. Whitehead, H., Rowan, C., & Kocsis, B. 2025b, MNRAS, submitted [arXiv:2505.23899] [Google Scholar]
  93. Wu, Y., Chen, Y.-X., & Lin, D. N. C. 2024, MNRAS, 528, L127 [Google Scholar]
  94. Yang, C.-C., Mac Low, M.-M., & Menou, K. 2009, ApJ, 707, 1233, publisher: The American Astronomical Society [Google Scholar]
  95. Yang, C.-C., Low, M.-M. M., & Menou, K. 2012, ApJ, 748, 79 , arXiv:1103.3268 [astro-ph] [Google Scholar]
  96. Zubovas, K., Tartènas, M., & Bourne, M. A. 2024, A&A, 691, A151 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  97. Zwick, L., Tiede, C., Trani, A. A., et al. 2024, Phys. Rev. D, 110, 103005 [Google Scholar]
  98. Zwick, L., Takátsy, J., Saini, P., et al. 2025, ApJ, 991, 131 [Google Scholar]

All Tables

Table 1.

Symbols used in the manuscript.

In the text

All Figures

thumbnail Fig. 1.

Dependency of Ostriker (1999) the dynamical friction prescription as a function of the mach number ℳ (Eq. (4)). The dotted line marks the linear interpolation used to circumvent the mathematical discontinuity at ℳ = 1. The dashed line indicates the linear subsonic limit.
In the text
thumbnail Fig. 2.

Radial profiles of the AGN disk employed in this work. From top to bottom: midplane gas density ρ, sound speed cs, disk aspect ratio h/R, and magnitude of the migration torque Γ. The migration torque assumes a secondary BH mass of mBH = 20 M⊙. The outer migration trap lies within the star formation region, which begins at R ≃ 767 au.
In the text
thumbnail Fig. 3.

Evolution of the semimajor axis (top), the eccentricity (middle), and the inclination (bottom) for an inclined (ι0 = 0.5°) and eccentric (e0 = 0.1) BH (mBH = 20 M⊙) undergoing dynamical friction in the AGN disk. The thick blue line indicates the evolution without the turbulent velocity field, while each thin red line includes a different realization of the turbulent velocity field. The BH is placed at the migration trap within an AGN disk around a 107 M⊙ SMBH. The dotted black line is the mean evolution of the realizations including turbulence. In the top panel, the dot-dashed lines indicate the mean square change in semimajor axis ⟨Δa2⟩ (Eq. (44)). In the middle and bottom panels, the dot-dashed green line is the ratio between the turbulent velocity dispersion σturb and the circular velocity vcirc. The vertical dashed lines indicate the first three midplane crossings. The turbulence prevents the full circularization and alignment of the embedded BH.
In the text
thumbnail Fig. 4.

Distributions (normalized to one) of the semimajor axis (top), the eccentricity (middle) and the inclination (bottom) for 2.5 × 104 realizations of an embedded BH, evaluated after 20 initial orbital periods. The initial conditions are identical to those in Fig. 3. In the top panel, the blue line marks the final semimajor axis from the nonturbulent simulation, while the green curve shows a Gaussian distribution with standard deviation given by Eq. (44). In the middle and bottom panels, the dot-dashed green line is the ratio between the turbulent velocity dispersion σturb and the circular velocity vcirc, while the green curves indicate Rayleigh distributions with mean value σturb/vcirc. The insets show the same distributions in logarithm scale, highlighting the presence of fatter tails with respect to the Rayleigh prediction.
In the text
thumbnail Fig. 5.

Ratio between the damping timescale τ* (Eq. (9)) and the orbital period P, as a function of the distance from the central SMBH. The estimate assumes a BH with mass mBH = 20 M⊙ based on the AGN disk profile for SMBHs with masses MSMBH = 106, 3 × 107, 7 × 106, and 2 × 107 M⊙, with α = 0.1. The vertical lines indicate the location of the migration traps.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.