Probing Formation Channels of Extreme Mass-Ratio Inspirals

TL;DR

This paper explores how measuring eccentricity, inclination, and mass in extreme mass-ratio inspirals can reveal their formation channels and underlying astrophysical environments, aiding future gravitational wave observations.

Contribution

It establishes theoretical expectations for observable properties in different formation channels of extreme mass-ratio inspirals and discusses their use in probing astrophysical environments.

Findings

Distinct distributions of observables can differentiate formation channels.

Eccentricity and inclination measurements inform about galactic nucleus environments.

Component mass data helps characterize black hole populations.

Abstract

The population study of stellar-mass black hole (sBH) binaries with ground-based gravitational wave detection has achieved tremendous success in recent years. Future observation of extreme mass-ratio inspirals will similarly require proper population analysis that identify the formation channels, measuring the branch ratio(s) and characterizing major properties within each major channel. In this work, we propose that the measurement of eccentricity, inclination, and component mass provides critical information to distinguish different formation channels and probe detailed formation mechanisms. Focusing on the dry and wet extreme mass-ratio inspirals, we establish the theoretical expectation of these observables in each formation channel. We also discuss how their distributions can be used to probe lifetime and turbulence level of active galactic nuclei disks, accretion patterns of supermassive black holes and population properties of sBHs within nuclear star clusters.

Tables

Table 1. Table 1: Comparison between Wet EMRIs and Dry EMRIs at r p = 10 ​ M ∙ r_{\rm p}=10M_{\bullet}

	Wet EMRI	Dry EMRI
Eccentricity	Multi-body dynamics: $\lesssim 0.01$, Fig. 7 Turbulence: $\simeq {10}^{-4}$, Fig. 16	$\gtrsim 0.01$, Fig. 2
Inclination	Coherent accretion: aligned orbit Chaotic accretion: tilted distribution, Fig. 19	Tilted distribution, Fig. 21
sBH Mass	Capture: larger 35$M_{\odot }$ peak Accretion: component mass increases Mergers: new peaks around 20$M_{\odot }$ and 45$M_{\odot }$ Fig. 28	Larger $35M_{\odot }$ peak, Fig. 25

Table 2. Table 2: Symbols

Physical Quantity	Physical Meaning	Physical Quantity	Physical Meaning
$M_{\bullet }$	Mass of center MBH	$t_{\mathrm{GW}}$	Gravitational wave damping timescale for semi-major axis
$M$	Total mass of the MBH and sBH	$t_{J,\mathrm{sca}}$	Gravitational scattering timescale
$\eta$	Symmetrical mass ratio	$\tilde{\gamma }$	Power-law index of MBH mass function
$m$	Mass of sBH	${\dot{M}}_{\mathrm{Edd}}$	Eddington accretion rate for the MBH
$q$	Mass ratio between sBH and MBH	${\dot{M}}_{\bullet }$	Mass accretion rate for the MBH
$\alpha$	Viscosity parameter	$h$	Disk aspect ratio
$E_{\mathrm{cr}}$	Critical energy distinguishing EMRI and Plunge	$H$	Disk scale height
$E_{\max }$	Maximum energy boundary in FP equation	$\Sigma \left(r\right),{\Sigma }_0$	Disk surface density at $r$ and $r_0$
$J_c\left(E\right)$	Angular momentum on circular orbit at given energy	${\Phi }_{\mathrm{turb}}$	Potential for the turbulent disk
$J_{\mathrm{lc}}$	Angular momentum on loss cone boundary	$\gamma$	Dimensionless strength for turbulence
$J$	Angular momentum of sBH	$\omega$ or $\Omega$	Angular velocity
$E$	Orbital energy of sBH	$\chi$	Dimensionless spin parameter
$R$	Dimensionless angular momentum variable in FP equation	$n_k$	Wavenumber of the turbulence mode
$N(E,R),f(E,R),C(E,R)$	Number density in $(E,R)$ phase-space ($N$) and position-velocity space ($f$), weight function ($C$)	$r_k,{\varphi }_k$	Radial position and azimuthal angle of each turbulence mode
$F_E,F_R$	Phase space flux along $E/R$ direction in FP equation	${\gamma }_e,{\tau }_e$	Eccentricity damping rate, timescale for the turbulent disk
$D_{EE},D_{RR},D_{ER},D_E,D_R$	Diffusion and advection coefficients	$S_{\delta r},S_{f_r},S_{{\tau }_{\varphi }}$	Power spectrum of radial motion, radial force, and azimuthal torque in the turbulent disk
$a,a_{\mathrm{in}},a_{\mathrm{out}}$	Semi-major axis of the (inner/outer) object	$t_J,t_{J,\mathrm{in}},t_{J,\mathrm{out}}$	Disk-driven damping timescale for angular momentum of (inner/outer) object
$e,e_{\mathrm{in}},e_{\mathrm{out}}$	Eccentricity of the (inner/outer) object	$t_e,t_{e,\mathrm{in}},t_{e,\mathrm{out}}$	Disk-driven damping timescale for eccentricity of (inner/outer) object
${\Gamma }_{\mathrm{EMRI}}$	EMRI rate	$\iota$	Inclination angle
${\varphi }_{\mathrm{in}},{\varphi }_{\mathrm{out}}$	Resonance angle of inner and outer objects	$f_r,{\tau }_{\varphi }$	Specific gravitational radial force and torque for the turbulent disk
${\varpi }_{\mathrm{in}},{\varpi }_{\mathrm{out}}$	Longitude of pericenter of inner and outer objects	$f_r\left(\omega \right),{\tau }_{\varphi }\left(\omega \right)$	Specific gravitational radial force and torque in the frequency domain for the turbulent disk
$r_a,r_p$	Apocenter, Pericenter of the sBH	$\delta F_r,\delta \Gamma$	Total stochastic radial force and torque for the turbulent disk
$ϵ$	Radiative efficiency of the AGN disk	$F_r,\Gamma$	Total radial force and torque for the turbulent disk
$R_{\mathrm{warp}}$	Warped disk radius	$P_{\mathrm{LVK}}$	LVK primary mass distribution
$R_{\mathrm{dec}}$	Decoupling radius where LT torque dominates over disk torques	$S\left(x\right)$	Smooth function
${\tau }_{\mathrm{BP}}$	Bardeen-Petterson timescale	$g\left(q\right)$	Pairing function
$J_d$	Disk angular momentum	$F(m_1,m_2)$	Joint distribution of binary mass
$J_{\bullet }$	MBH spin angular momentum	$\tilde{\beta }$	Power-law index of pairing function
${\alpha }_p,{\mu }_m,{\sigma }_m,{\delta }_m,{\lambda }_{\mathrm{peak}}$	Hyper-parameters of LVK primary BH mass distribution	${\tilde{\alpha }}_p,{\tilde{\mu }}_m,{\tilde{\sigma }}_m,{\tilde{\delta }}_m,{\tilde{\lambda }}_{\mathrm{peak}}$	Hyper-parameters of single BH mass distribution
$P_{\mathrm{cap}}$	Wet EMRI mass distribution after being captured	$P_{\mathrm{acc}}$	Wet EMRI mass distribution after accreting
$f_{\mathrm{mer}}$	Merger fraction	$P_{f=1}$	Wet EMRI mass distribution after merging for $f_{\mathrm{mer}}=1$
$q_{\mathrm{crit}}$	Critical mass ratio between type-I and type-II migration	$P_{\mathrm{mer}}$	Wet EMRI mass distribution after merging
${\dot{J}}_{\mathrm{mig},I}$	Type-I migration torque	${\dot{J}}_{\mathrm{gw}}$	Gravitational radiation torque
${\dot{J}}_{\mathrm{wind}}$	Wind torque	${\dot{m}}_{\mathrm{wind}}$	Effective headwind mass flux
${\dot{J}}_{\mathrm{mig},\mathrm{II}}$	Type-II migration torque	${\dot{m}}_{\mathrm{in}}$	Mass inflow at outer boundary of the sBH
${\dot{m}}_{\mathrm{std}}$	Steady accretion rate of AGN disk	$T_{\mathrm{disk}}$	Disk lifetime for AGN
$\Gamma (m,t)$	Per-object EMRI rate	${\Gamma }_{\mathrm{cap}}(m,t)$	Per-object capture rate

Equations

$$N(E,R)dEdR={\int }_{r_p}^{r_a}{d}^{3}r{d}^{3}vf(E,R)$$

$$N(E,R)$$

$$:$$

$$C(E,R)\frac{\partial f}{\partial t}=-\frac{\partial F_E}{\partial E}-\frac{\partial F_R}{\partial R}.$$

$$-F_E$$

$$-F_R$$

$$t_{\mathrm{GW}}$$

$$t_{J,\mathrm{sca}}$$

$${\Gamma }_{\text{EMRI}}={\int }_{E_{\text{cr}}}^{E_{\text{max}}}dE\left(-F_R+F_E\frac{dR}{dE}\right),$$

$$P(\log e_0)={\Gamma }_{\mathrm{EMRI}}^{-1}\left(-F_R+F_E\frac{dR}{dE}\right)\frac{dE}{d\log e_0}.$$

$$\frac{dn}{d\log M_{\bullet }}=\tilde{C}{\left(\frac{M_{\bullet }}{10^6{M}_{\odot }}\right)}^{\tilde{\gamma }}{\text{Mpc}}^{-3},$$

$$⟨P(e_0){⟩}_{M_{\bullet }}=\frac{\int d{M}_{\bullet }{\Gamma }_{\mathrm{EMRI}}(M_{\bullet })P(e_0)\frac{dn}{d{M}_{\bullet }}}{\int d{M}_{\bullet }{\Gamma }_{\mathrm{EMRI}}(M_{\bullet })\frac{dn}{d{M}_{\bullet }}},$$

$$h\simeq 0.7{\dot{M}}_{\bullet ,0.1}\frac{M_{\bullet }}{r}$$

$$\Sigma (r)={\Sigma }_0{\left(\frac{r}{100M_{\bullet }}\right)}^n.$$

$$t_J$$

$$t_e$$

$${\vec{a}}_{\text{PN}}=$$

$$+\frac{8}{5}\eta \frac{M}{r^2}\frac{M}{r}\left[v_r\frac{\vec{r}}{r}{\text{A}}_{\text{2.5PN}}-\vec{v}{\text{B}}_{\text{2.5PN}}\right],$$

$${\vec{a}}_{\text{mig}}=-\frac{\vec{v}}{t_J}-2\frac{v_r}{t_e}\frac{\vec{r}}{r}.$$

$${\varphi }_{\text{in}}$$

$${\varphi }_{\text{out}}$$

$$(p+q){\omega }_{\varphi ,\mathrm{out}}-(p+q){\omega }_{\varphi ,\mathrm{in}}+q{\omega }_{r,\mathrm{in}}=0,$$

$$p{\omega }_{\varphi ,\mathrm{out}}-p{\omega }_{\varphi ,\mathrm{in}}+q{\omega }_{r,\mathrm{out}}=0,$$

$$e=\frac{r_a-r_p}{r_a+r_p}$$

$$k_1{\omega }_{\varphi ,\text{in}}+k_2{\omega }_{r,\text{in}}+k_3{\omega }_{\varphi ,\text{out}}+k_4{\omega }_{r,\text{out}}=0$$

$${\Phi }_{\mathrm{turb}}(r,\varphi ,t)=\gamma r^2{\Omega }^2\sum _{k=1}^{50}{\Lambda }_k(n_k,r,\phi ,t),$$

$${\Lambda }_k={\xi }_k{e}^{-(r-r_k{)}^2/{\sigma }_k^2}\cos (n_k\varphi -{\varphi }_k-{\Omega }_k{\tilde{t}}_k)\sin (\pi {\tilde{t}}_k/\Delta t_k),$$

$$\alpha \simeq 35(\gamma /h_0{)}^2.$$

$$\ddot{r}-r{\dot{\varphi }}^2+\frac{M_{\bullet }}{r^2}=f_r,$$

$$\ddot{\varphi }{r}^2+2\dot{\varphi }\dot{r}r=t_{\varphi },$$

$$r=r_0+\delta r,\varphi ={\Omega }_{0}t+\delta \varphi ,$$

$$\ddot{r}-3{\Omega }_0^2\delta r-2r_0{\Omega }_0\delta \dot{\varphi }=f_r,$$

$$\delta \ddot{\varphi }{r}_0^2+2{\Omega }_0{r}_0\delta \dot{r}=t_{\varphi }.$$

$$-{\omega }^2\delta r-3{\Omega }_0^2\delta r+2i{r}_0{\Omega }_0\omega \delta \varphi =f_r(\omega ),$$

Full text

I Introduction

The study of gravitational wave source populations generally provides a demographic description of the source distribution for each formation channel and offers critical insights into key aspects of individual formation mechanisms. One successful example is the measurement of binary black hole (BBH) merger events [1, 2, 3] , which utilizes information about distances, component masses, and spins to begin inferring the source redshift distributions, branching ratios of various formation mechanisms, delay times, etc. With the rapidly increasing number of detection events from LIGO-Virgo-KAGRA observations, many important questions related to the formation of BBHs are expected to see significant progress in the near future [4, 5, 6, 7, 8].

For space-borne gravitational wave detectors [9, 10], extreme mass-ratio inspirals (EMRIs) are among the primary extragalactic sources. Identifying and characterizing their formation channels are key scientific goals of these missions. Several major formation channels have been proposed. First, within the nuclear star cluster surrounding a massive black hole (MBH), stellar-mass black holes (sBHs) can be gravitationally scattered into extremely low-angular-momentum orbits. In such cases, gravitational wave emission during pericenter passages may efficiently damp these highly eccentric orbits before two-body scattering disrupts their trajectory toward the so-called “loss cone” [11]. These “dry EMRIs” may contribute $\mathcal{O}(1)-\mathcal{O}(10^{3})$ events to space-borne detectors, with main uncertainties arising from the mass function of MBHs, the composition and distribution of nuclear star clusters [9]. Second, in a small fraction ($\mathcal{O}(1\%)-\mathcal{O}(10\%)$), according to observations [12] of galaxies that host active galactic nuclei (AGNs), sBHs may be captured through interactions with the accretion disks and subsequently migrate toward the central MBHs. It turns out that disk-driven migration is a more efficient mechanism for generating EMRIs than the “loss-cone” mechanism, such that the overall number of “wet EMRIs” may be comparable to or even greater than that of “dry EMRIs” [13, 14, 15, 16], despite the small fraction of galaxies hosting AGNs. The main uncertainties in the rate predictions arise from the modeling of nuclear star clusters, the mass function of MBHs, and the structure of accretion disks. Third, additional formation channels for EMRIs have been proposed, including the tidal disruption of sBH binaries near MBHs (the “Hills” mechanism [17]) and the Kozai-Lidov mechanism in systems with supermassive black hole (SMBH) binaries [18]. Determining the rate of the Hills channel requires a quantitative population study of the number and distribution of sBH binaries within nuclear star clusters, as well as the probability of binary tidal disruption and subsequent EMRI formation. Previous simulations show that Hills EMRIs are actually quite similar to dry EMRIs in the sensitivity band of space-borne gravitational wave missions [19]. On the other hand, the EMRI formation rate in systems with SMBH binaries remains under debate [20, 21]. Therefore, in this work, we focus on the dry and wet EMRI channels, leaving the discussion of additional formation pathways to future studies when more quantitative predictions become available.

To diagnose EMRI formation channels, a primary objective is to identify the key observables and their distributions for each channel. For example, the spin magnitudes and orientations of sBHs provide crucial information for distinguishing field binaries from those formed in dynamical environments, particularly in the context of ground-based detectors [22]. For EMRIs, we will show that while dry EMRIs naturally retain eccentricities of $\geq 0.01$ in the Laser Interferometer Space Antenna (LISA) band, wet EMRIs are not perfectly circular either. Two possible mechanisms for exciting eccentricity in wet EMRIs are likely to operate: turbulence in the AGN disk and multi-body resonance effects within the disk. As a result, the eccentricities of wet EMRIs may be pumped to $\mathcal{O}(10^{-4})$-$\mathcal{O}(10^{-2})$ in the LISA band, depending on the level of turbulence and whether multiple sBHs migrate together. It should be noted that the eccentricity measurement precision for LISA, Taiji, and TianQin is expected to reach $\mathcal{O}(10^{-5})$ [9]. Such high precision may not only help distinguish different formation channels, but also enable the characterization of key channel properties, such as the relevant surface density and average turbulence strength in AGN disks.

Recent studies have explored the possibility that accretion disks may directly introduce observable phase shifts into EMRI waveforms, which potentially enables the direct detection of disk effects from wet EMRIs [23, 24]. This necessarily requires a rather high disk surface density ($\gtrsim 10^{7}\ \mathrm{g\ cm^{-2}}$) in order to produce sufficiently strong effects [23]. This requirement cannot be satisfied by the canonical $\alpha$-disk scenario, but may be fulfilled by one of its variants, the so-called “$\beta$-disk” model [25]. The plausibility of this scenario remains to be tested by further numerical simulations and future observations. However, population studies involving eccentricity and inclination angle do not require a direct disk imprint on the waveform. Instead, they infer disk effects during the earlier stages of wet EMRIs by measuring orbital parameters. This type of indirect measurement offers a unique opportunity to probe the formation and early dynamical evolution of EMRIs.

In addition to eccentricity, the inclination angle relative to the spin axis of the central MBH is also a key observable in population studies. Our results show that dry EMRIs tend to favor prograde orbits, with the event rate varying by more than an order of magnitude across different inclination angles. The detailed distribution of the inclination angle is derived and can be directly compared with future observational data. For wet EMRIs, the inclination angle distribution depends on whether the MBH’s spin is more chaotic or coherent [26], as well as on how the disk lifetime compares to the Bardeen-Petterson timescale [27].

We further investigate the mass distribution of sBHs associated with different formation channels, each characterized by distinct signatures. In particular, for BBHs detected by the LIGO-Virgo Collaboration, the component mass distribution shows prominent peaks near $10\ M_{\odot}$ and $35\ M_{\odot}$. We adopt a phenomenological pairing function that characterizes the efficiency of binary BH formation from isolated sBHs, and use it to infer the underlying mass function of isolated sBHs. For dry EMRIs, the mass function of their sBHs exhibits more prominent peaks around $35\,M_{\odot}$ compared to the mass function of isolated sBHs, due to the more efficient gravitational scattering of more massive sBHs and the effect of mass segregation. For wet EMRIs, the mass function is further shaped by additional processes, including the capture of sBHs by the AGN disk, mergers of sBHs in the AGN disk, and mass growth due to gas accretion. Multiple new features on the mass function may appear, as summarized in Table 1, along with the main predictions of the distributions of the observables considered in the wet and dry formation channels. Other observables, such as the spins of secondary black holes, have also been proposed to probe different formation channels [28]. However, the in situ distribution of isolated sBHs remains unknown, and the accretion processes of these objects within AGN disks still require systematic characterization through numerical simulations. Moreover, secondary spin measurement precision degrades dramatically for low-eccentricity systems (such as wet EMRIs) [28]. As a result, we do not perform quantitative analysis in spin distributions in this work.

This paper is organized as follows. In Section II, we discuss the eccentricity distribution of EMRIs, covering dry EMRIs in Section II.1 and wet EMRIs in Section II.2. For the latter, we examine two mechanisms that can excite eccentricity: multi-body resonances (Section II.2.1) and turbulence (Section II.2.2). Section III focuses on the distribution of inclination angles, with dry and wet EMRIs discussed in Sections III.2 and III.1, respectively. In Section IV, we analyze the mass distribution, with corresponding discussions for dry and wet EMRIs in Sections IV.1 and IV.2. Finally, Section V summarizes our conclusions and outlines the prospects for future studies of EMRI populations. If no special instructions are given, in this paper we use the $G=c=1$ units, and we adopt the notation $M_{\bullet}$ for the mass of the MBH. Table LABEL:tab:symbols (in Appendix. E) provides a summary of the physical quantities employed in this study.

II Eccentricity

In this section, we investigate the expected eccentricity distributions associated with the dry and wet EMRI formation channels. In particular, since EMRI waveforms typically consist of a large number of cycles ($\gtrsim 10^{4}$), even small deviations in the intrinsic parameters of the system can accumulate and result in substantial phase differences at the end of the waveform. This implies that the measurement precision of many system parameters for EMRIs is significantly higher than that for comparable-mass binary systems. For example, as noted in [9], the eccentricity of a typical EMRI event can be measured with precision of the order of $10^{-5}$. This exceptional accuracy not only enables a clear distinction between EMRIs originating from different formation channels but also opens the possibility of probing properties of the EMRI environment during the early stages of its evolution, as will be discussed later in this section.

II.1 Dry EMRI

The dry EMRI channel, originating from two-body gravitational scatterings in nuclear star clusters, was the first major EMRI formation mechanism proposed for detection by space-based gravitational wave observatories. Several studies have been conducted on the formation rates and eccentricity distributions of EMRI [29, 30, 31, 32, 9]. For completeness, in this section we briefly outline the main procedure for computing the distribution of dry EMRIs using the Fokker-Planck (FP) method, with particular emphasis on the eccentricity distribution. A more detailed discussion of the formalism can be found in [14].

Within the FP formalism, the phase-space distribution of stars and sBHs in a nuclear star cluster is modeled under the assumption of spherical symmetry. Each object’s orbit is characterized by its energy $E$ and angular momentum $J$. For convenience, the angular momentum is often re-parameterized using the dimensionless variable $R=J^{2}/J_{c}^{2}(E)$, where $J_{c}(E)$ is the angular momentum of a circular orbit with energy $E$. This choice simplifies the treatment of diffusion in angular momentum space, as $R\in[0,1]$.

There is a distinct region in the phase space $(E,R)$, bounded by the energy-dependent angular momentum $J_{\rm lc}(E)$, as shown by the orange area in Fig. 1. Objects such as stellar-mass black holes that initially occupy this region will plunge into the MBH within a single orbital period. This region is known as the “loss cone” [11], and its boundary, $J_{\rm lc}(E)$, referred to as the “loss-cone boundary” (i.e., the last stable orbit), is represented by the orange line in Fig. 1.

Since dry EMRIs typically have eccentricities close to unity, it is a good approximation to define the loss cone boundary using $J_{\rm lc}=4M_{\bullet}$ [11], assuming a non-spinning (Schwarzschild) MBH. This expression becomes exact in the limit of extreme eccentricity. Note that in this figure we rescale the specific energy using $\sigma$, defined as the stellar velocity dispersion following [14, 33].

Assume that $f$ denotes the number density in position-velocity space, and $N$ represents the distribution function in the $(E,R)$ phase space. The two are connected through the following relation [34, 35]:

[TABLE]

in which $r_{p}(E,R)$ and $r_{a}(E,R)$ are the orbital pericenter and apocenter, respectively [14]. This number density relation can be translated to

[TABLE]

where $P(E,R)$ is the orbital period and $C(E,R)$ is often referred to as the weight function. The FP equation that describes the time-dependent evolution of the distribution is

[TABLE]

Here, $F_{E,R}$ is the phase-space flux along the $E/R$ direction. They can be obtained from

[TABLE]

The diffusion coefficients $\{D_{EE},D_{ER},D_{RR}\}$ and the advection coefficients $\{D_{E},D_{R}\}$ can be found in [14, 36].

There are two key timescales that govern the evolution of a nuclear star cluster: the gravitational scattering timescale, $t_{J,\rm sca}$, which characterizes the rate at which an object’s angular momentum changes due to interactions with other objects in the cluster; and the gravitational wave radiation timescale, $t_{\rm GW}$, which describes the orbital evolution of a two-body system composed of the object and the central MBH, driven by gravitational wave emission. Based on [36, 29, 37, 38, 39], these two timescales are

[TABLE]

where $t_{\rm GW}$ above is modified Peter’s formula for a highly eccentric binary system with semi-major axis $a$, eccentricity $e$, and mass components $M_{\bullet}$ and $m$.

In the $(E,R)$ phase space, it is instructive to define a line where the gravitational wave inspiral timescale equals the angular momentum relaxation timescale, i.e., $t_{\rm GW}=t_{J,\rm sca}$. This condition is shown as the blue line in Fig. 1. It also intersects the loss-cone boundary, and the corresponding energy at the intersection is denoted by $E_{\rm cr}$. It is evident that in the region between the loss cone boundary and the blue line, the gravitational-wave timescale $t_{\rm GW}$ is shorter than the angular momentum diffusion timescale $t_{J,\rm sca}$, implying that gravitational-wave emission dominates the secular orbital evolution. This region is commonly referred to as the “EMRI regime”, where a highly eccentric EMRI initially forms and subsequently evolves toward the loss cone along the dashed trajectories shown in Fig. 1. To compute the EMRI formation rate, one may evaluate the integrated flux along either the blue line or at the loss-cone boundary. In a quasi-steady state, both approaches yield the same result due to flux conservation:

[TABLE]

where $E_{\rm max}$ corresponds to the point where the loss-cone boundary intersects the boundary $R=1$, i.e. the circular orbit with angular momentum $J_{c}(E)=J_{\rm lc}=4M_{\bullet}$.

For a typical dry EMRI, the eccentricity is initially close to unity. During repeated pericenter passages, the system emits substantial gravitational radiation, resulting in the loss of both energy and angular momentum. Consequently, the eccentricity decreases rapidly from cycle to cycle. Because the initial semi-major axis is large compared to the size of the central MBH (and thus the system is only weakly bound), the semi-major axis $a$ is more strongly affected by gravitational radiation than the pericenter distance $r_{p}$. As a result, in the early stages of evolution, $a$ can decrease by orders of magnitude while $r_{p}$ experiences only modest fractional changes. Once the eccentricity drops below approximately 0.9, the waveform enters the valid parameter range of state-of-the-art waveform models, such as fastemriwaveforms (FEW) [40].

The eccentricity distribution of dry EMRIs can be derived from the phase-space distribution function in $(E,R)$, obtained by solving the FP equation. At each point in the phase space characterized by $(E,R)$, the corresponding pericenter distance $r_{p}(E,R)$ and apocenter distance $r_{a}(E,R)$ can be calculated. The semimajor axis a and eccentricity e are then given by $a=(r_{a}+r_{p})/2,\quad e=(r_{a}-r_{p})/(r_{a}+r_{p})$. Within the EMRI regime, as the orbital parameters $(e,a)$ evolve under gravitational wave-driven dynamics [41], we define the reference eccentricity of each trajectory as $e_{0}:=e|_{r_{p}=10M_{\bullet}}$. The eccentricity distribution is then obtained by evaluating the EMRI rate -measured at the boundary where $t_{\rm GW}=t_{J,\rm sca}$ - for different trajectories such as

[TABLE]

For an MBH with mass $M_{\bullet}$, the resulting time-averaged distribution of $e_{0}$ is shown in the left panel of Fig. 2. We find that the eccentricities of dry EMRIs fall predominantly within the range $10^{-2}$ to $1$, representing more than $95\%$ of the total population. This is consistent with the results in [29].

To assess the overall statistical distribution of the dry EMRI eccentricity across the population of MBHs in the universe, we incorporate the MBH mass function. Specifically, we assume a parametric form for the comoving number density of MBHs as a function of the MBH mass given by

[TABLE]

where $\tilde{\gamma}$ is a power law index that characterizes the slope of the MBH mass function and $\tilde{C}$ is a constant that does not affect the eccentricity distribution. We consider the range of index from $\tilde{\gamma}=-0.3$ to $\tilde{\gamma}=0.3$ [9] to reflect plausible uncertainties in the MBH population models inferred from observations and simulations.

The mass-weighted average eccentricity distribution for a given mass function is

[TABLE]

in which $\Gamma_{\rm EMRI}(M_{\bullet})$ calculated from the FP equation, and we have only considered MBHs in the mass range $\in(10^{5},10^{7})M_{\odot}$, relevant for spaceborne GW detection. The overall distribution is shown in the right panel of Fig. 2. In general, the change in the power law index $\tilde{\gamma}$ only slightly shifts the distribution around $e_{0}\sim 10^{-2}$.

We also compare the eccentricity distribution computed here with a recent study by [42]. Although there are a few differences in the setup, the results are broadly consistent. More detailed comparison can be found in Appendix. D.

II.2 Wet EMRI

Wet EMRIs, on the other hand, may exhibit diminishing eccentricities according to standard disk migration theory. However, we will show that two novel mechanisms can potentially excite the eccentricities of wet EMRIs: one arising from multibody interactions and resonances within the accretion disk (see Sec. II.2.1), and the other from turbulent fluctuations in the disk (see Sec. II.2.2). Accurate eccentricity measurements offer a valuable opportunity to probe the evolutionary history and environmental conditions of wet EMRIs.

We consider the orbital evolution of EMRI using parametrized AGN disk models. For an AGN disk, the disk profiles within $100\ M_{\bullet}$ can be parameterized as [23]

[TABLE]

and

[TABLE]

For the $\alpha$ model with $\alpha=0.1$ and $\dot{M}_{0.1,\bullet}=\dot{M}_{\bullet}/0.1\dot{M}_{\rm Edd}$ where $\dot{M}_{\rm Edd}$ is the Eddington accretion rate for the MBH, we have $n=1.5,\Sigma_{0}=7\times 10^{4}\ {\rm g\ cm^{-2}}$. For the $\beta$ model with the same model parameters, we have $n\simeq-0.6,\Sigma_{0}=6\times 10^{5}\ {\rm g\ cm^{-2}}$. The disk profiles with different $\alpha$ and $\dot{M}_{\bullet}$ for $\alpha$ and $\beta$ models are shown in Fig. 3. The real disk profiles may deviate from both $\alpha$- and $\beta$-disk prescriptions, so it is instructive to consider the parametrized form of the disk, especially for the discussion in Sec. II.2.1 and Sec. II.2.2.

II.2.1 Multi-body resonance

An isolated sBH embedded in a thin accretion disk generically excites density waves that carry away energy and angular momentum, analogous to the behavior of planets in protoplanetary disks. The corresponding damping timescales for angular momentum ($t_{J}$) and eccentricity ($t_{e}$) are

[TABLE]

where $m$ is the sBH mass, $\Sigma$ is the disk surface density, $h$ is the disk aspect ratio, and $\Omega$ is Keplerian angular velocity. All these quantities related to the disk are measured at the semi-major axis $a$ of the sBH orbit. For a thin disk the damping timescale of eccentricity is much smaller than that of the energy and angular momentum (by a factor $h^{2}\ll 1$), so the eccentricity would quickly decrease below $10^{-5}$, as can be seen from hydrodynamical simulations [43, 44]. However, it is instructive to examine the underlying assumption of this argument, e.g., an isolated sBH, as multibody interactions may significantly modify the dynamics of sBHs.

In fact, since migration within the AGN is dominated by disk interactions at large radii and gravitational wave back reaction at small radii, the effective lifetime (measured by $a/\dot{a}$) is maximized at a few hundred gravitational radii, according to the analysis in [14] (see also Fig. 17 below). The refined population study (following a similar formalism) in [45] predicts that there are on average one or more sBH staying at a few hundred gravitational radii. This number is larger for shorter disk lifetimes, as the wet EMRI rate generally decreases over time due to decreasing supply of sBHs [14, 36]. Therefore, it is natural to expect that a fraction of wet EMRIs are accompanied by close companions. We will need to jointly evolve the multibody system to determine its long-term dynamics. In this work, for the sake of simplicity, we focus on three-body systems with two sBHs initially placed at two different radii within a thin accretion disk. The motion is assumed to be planar.

To numerically evolve the three-body system, we adopt the N-body code REBOUND [46] with the following modifications. First, since the later stages of the EMRI trajectory enter the strong-gravity regime, we include the 1PN, 2PN, and 2.5PN Post-Newtonian (PN) corrections to the equations of motion of the sBHs [47]:

[TABLE]

where $M=M_{\bullet}+m$ is the total mass of center MBH and sBH, $\eta=M_{\bullet}m/(M_{\bullet}+m)^{2}$ is symmetrical mass ratio, $v_{r}$ is the radial velocity of the sBH, $\vec{r}$ and $\vec{v}$ are the position and velocity vector of the sBH. A and B are the PN coefficients that are presented in Appendix. A. Notice that we have not included PN corrections to the gravitational force between the two sBHs, because this force at Newtonian order is already $m/M_{\bullet}$ times smaller than the force between an sBH and the MBH.

Second, we also include the disk force in the equation of motion. Assuming the characteristic timescales $t_{\rm J},t_{\rm e}$, the effective acceleration term can be written as

[TABLE]

We apply the approximation formula in Eq. (II.2.1) to implement such disk-force component.

In the following analysis, we adopt a central MBH mass of $M_{\bullet}=4\times 10^{6}M_{\odot}$ and use $\alpha$-disk parameters with $\alpha=0.1$ and $\dot{M}_{\bullet}=0.1\dot{M}_{\rm Edd}$ [14].

II.2.1.1 Dynamics.

Depending on initial conditions of the system, the multibody dynamics can be classified into resonance, non-resonant migration, and scattering regimes. In principle, dynamical formation of sBH binaries can also be classified as the scattering case, as it typically involves close gravitational encounters. However, as simulations [48, 49, 50] suggest that it often requires the modeling of circum-single disks, we shall not consider this scenario in this study.

Resonance locking generally happens if the outer object migrates faster (e.g., being the more massive object) than the inner object, so that the relative motion is convergent. Similarly to mean-motion resonances in planetary dynamics [51], we can define the resonant angles for inner and outer resonances as

[TABLE]

where the $\lambda$ is mean longitude and $\varpi$ is longitude of pericenter of true physical orbit. A system is trapped within a $(p+q):p$ inner/outer mean-motion resonance if $\phi_{\rm in}/\phi_{\rm out}$ is bound around a fixed value. In this case, notice that the time-averaged derivative of $\lambda$ corresponds to $\omega_{\phi}$ and the time-averaged derivative of $\lambda-\varpi$ corresponds to $\omega_{r}$, where $\omega_{\phi}$ and $\omega_{r}$ denote the angular frequency in the $\phi$ and $r$ direction, respectively. These two frequencies are the same in Newtonian theory, but start to differ from each other when PN corrections are included, because of the PN pericenter precession effect. The bounded value of $\phi_{\rm in}/\phi_{\rm out}$ also implies that

[TABLE]

for the inner resonance and

[TABLE]

for the outer resonance. In the following discussion, we will use $(p,q)_{\rm in}$ and $(p,q)_{\rm out}$ to denote these resonances.

In reality, because of the dissipation due to disk density waves and gravitational wave emission, the system gradually loses energy and angular momentum, causing sBHs to migrate towards the MBH. If the migration timescale is longer than the liberation timescale of the particular resonance, the resonance state tends to be maintained. An example evolution is shown in Fig. 4, for which a pair of sBHs are trapped in both $(8,1)_{\rm in}$ and $(8,1)_{\rm out}$ resonances (as indicated by the resonant angles $\phi_{\rm in,out}$) and jointly migrate inward. The angular period ratio remains approximately constant.

The evolution of eccentricities of the inner and outer sBHs also shows distinctive features. They are defined as

[TABLE]

computed with the pericenter $r_{p}$ and apocenter $r_{a}$ of an orbit. Unlike isolated sBH migrations, we observe significant eccentricities after the initial stage without resonance locking. As shown by the analytical theory in Appendix. B, the “equilibrium” eccentricity at any time depends on the particular resonance trapped in, and also on the decay time scales $t_{J},t_{e}$. The migration dynamics competes with the resonance dynamics between the two sBHs, as the migration timescale becomes smaller compared with the libration timescale, the eccentricity increases in time. This behavior of eccentricity evolution was previously seen in dynamical theories developed for planet migrations [51], and is still applicable to the sBH migration in AGN disks [52]. When the sBH migrates inward further, the PN pericenter precession starts to modify the resonance dynamics. As we include the PN precession term in the analytical theory, we observe that the inner “equilibrium eccentricity” indeed starts to decrease within a certain radius, which matches the numerical evolution of $e_{\rm in}$ in Fig. 4. The difference between analytical theory and numerical evolution between $e_{\rm in}$ and $e_{\rm in,equi}$ ($e_{\rm out}$ and $e_{\rm out,equi}$ show remarkable agreement) may arise from the simplified treatment in analytical theory, as only the 1PN precession term is included therein. However, this is the first time that a PN effect on the eccentricity of mean-motion resonances is identified, leading to eccentricity $\sim 10^{-3}$ as the inner sBH reaches $10M_{\bullet}$ for this particular system.

When the two sBHs are not locked in resonance, their period ratios no longer remain constant, and eccentricities decrease over time because of dissipation from the disk and gravitational wave emission. The mutual gravitational force between sBHs in general does not contribute secular effects at the leading order, except during the tidal resonances with [52]

[TABLE]

where $(k_{1},k_{2},k_{3},k_{4})$ are integers. Notice that in the most general cases all six frequencies of $\omega_{r},\omega_{\theta},\omega_{\phi}$ are allowed in the above relation, if the restriction on planar motion is removed. The conserved quantities of motion may undergo finite changes within any individual tidal resonance, so that one observes eccentricity excitations across tidal resonances.

It is instructive to show a few such examples. In Fig. 5 a pair of sBHs migrate out of resonance except briefly locked into the $5:2$ resonance around $r_{\rm in}\sim 100M_{\bullet}$. During the non-resonant phase, the eccentricities decrease over time, but significant eccentricity excitations are also observed across at least three transient tidal resonances. The inner sBH’s eccentricity is around $10^{-3}$ at $r=10M_{\bullet}$.

In Fig. 6 another pair of sBHs is shown, exhibiting a clear separation between the mean-motion resonance phase and non-resonant phase. During the $4:3$ resonance phase we observe behavior similar to that in Fig. 4. During the non-resonant phase we observe both the eccentricity damping when the system is away from major tidal resonances and eccentricity pumping across tidal resonances. The final eccentricity of the inner object at $r=10M_{\bullet}$ is around $0.01$.

In addition to the resonant and non-resonant regimes, a scattering regime arises when the two sBHs are in close proximity. As illustrated in Fig. 4, at very small semi-major axes the sustained resonance further reduces their separation, leading to strong mutual scattering. This interaction excites the eccentricity of the inner sBH to relatively high values, up to $\sim 0.1$, and breaks their sustained resonance.

II.2.1.2 Population.

To study the rich phenomena of multibody dynamics, we perform a statistical simulation of 100 examples. The central MBH has mass $M_{\bullet}=4\times 10^{6}M_{\odot}$, the surface density of disk $\Sigma(r)$ and aspect ratio $h(r)$ are chosen as the same setup as the $\alpha$-disk in [14] with $\dot{M}_{\bullet}=0.1\dot{M}_{\rm Edd}$, $\alpha=0.1$ as shown in Fig. 3. The component mass of sBHs are randomly chosen in the range $\{m_{\text{in}},m_{\text{out}}\}\in(10,30)M_{\odot}$, the initial semi-major axis $(a_{\text{in}},a_{\text{out}})$ in the range $(400,1000)M_{\bullet}$, the initial eccentricities $(e_{\text{in}},e_{\text{out}})$ in the range $(0.0,0.2)$, the initial true anomaly $(f_{\text{in}},f_{\text{out}})$ in the range $(0,2\pi)$. The initial longitudes of pericenter are set to zero, and the ascending nodes are also zero due to planar motion. The simulation is terminated when the inner sBH migrates to $r=8M_{\bullet}$. In order to quantity the eccentricity distribution, we measure the benchmark eccentricity for each system when the inner object migrates to $r=10M_{\bullet}$.

Under this set of initial conditions, there are cases where the two sBHs start in close proximity and undergo strong gravitational scattering, resulting in one of them being ejected and becoming unbound. Such scenarios are excluded from the population study, as our focus is on bound states of two sBHs.

The configurations of the systems in their final state are shown in Fig. 7. First, we observe a gap in the period ratios for the sBH pair. The right branch corresponds to systems that remain in mean motion resonances for most of the evolution time, with period ratios closely matching the resonance values. Since locking into and maintaining resonance generally requires faster intrinsic migration of the outer object, we find that in the right branch, the outer objects are typically more massive than the inner ones. The left branch, on the other hand, is contributed by pairs that have never been locked into a mean-motion resonance or pairs that are locked into a resonance briefly and escape the resonance midway, so that the period ratios continuously change over time as the inner object migrates faster than the outer object. For systems with large $m_{\rm in}/m_{\rm out}$, the influence of tidal resonances on the outer object is small (though not necessarily on the inner object), allowing the outer eccentricity to be damped to very small values by the end of the evolution. The gap between the left and right branches results from the relatively short time systems spend in these period ratios during the evolution.

The complete sample of eccentricities for the inner object at $r=10M_{\bullet}$ is summarized in the histogram shown in the right panel of Fig. 7. We find that resonance effects generally lead to $e_{\rm in}\sim(10^{-3},10^{-2})$ when the inner object enters the sensitivity band of space-borne detectors.

To examine the influence of the disk model on the final eccentricity, we also consider a $\beta$-disk with $\alpha=0.1$ and $\dot{M}_{\bullet}=0.1\dot{M}_{\rm Edd}$, using the same initial conditions for the 100 realizations of $\alpha$-disk model. The result, shown in Fig. 8, exhibits a distribution similar to that of the $\alpha$-disk in the right panel of Fig. 7, suggesting that the disk model has only a minor influence compared to the PN effect.

II.2.2 Turbulence

AGN disks are expected to exhibit ubiquitous strong turbulence driven by magneto-rotational instability [53]. The turbulent eddies associated with such disks exert stochastic forces on embedded sBHs, thereby exciting their orbital eccentricities. In order to quantify the eccentricity evolution of sBHs embedded in AGN disks, we first perform a series of hydrodynamical simulations tracking the evolution of a single sBH in a turbulent disk. Based on these hydrodynamical simulations, we can extract the scaling relations for the turbulence forces acting on the sBH. By incorporating these prescribed scaling relations into an N-body code, we can then investigate the long-term evolution of the sBH and infer the eccentricity when the EMRI enters the sensitivity band of detectors, for example $r=10\ M_{\bullet}$. The inferred eccentricity of the EMRI typically depends on the properties of the turbulent disks. Conversely, these eccentricities can be used to constrain key disk properties, such as the turbulence level, surface density, and scale height, etc.

II.2.2.1 Numerical Setup for Hydrodynamical Simulations.

To explore EMRI dynamics within an AGN disk, we use the hydrodynamic grid-based code FARGO3D [54] in our simulations.

For the disk model, we adopt a locally isothermal disk model with a constant aspect ratio $h\equiv H/r=h_{0}$. We also set the initial disk gas surface density profile to be $\Sigma=\Sigma_{0}(r/r_{0})^{-1/2}\left[M_{\bullet}/r_{0}^{2}\right]$. Here $\Sigma_{0}$ is a normalization constant relevant to disk mass, $r_{0}$ is the code unit length and corresponds to the initial inner planet’s orbital radius, $M_{\bullet}$ is the mass of the central MBH and the code unit for mass, and $\Omega=\sqrt{GM_{\bullet}/r^{3}}$ is the Keplerian angular frequency. Note that even though the slopes for the surface density and aspect ratio are not exactly the same as the $\alpha$- or $\beta$-disk models presented in Fig. 3, they should not significantly influence our hydrodynamical simulation results. Given the low disk masses considered in our simulations, we neglect the self-gravity of the disk. The indirect term associated with the motion of the central object is included, as the central MBH is fixed at the origin of the reference frame.

We do not include the classical $\alpha$ viscosity [55] in our disk model. Instead, we incorporate a phenomenological turbulence prescription into the accretion disk [56, 57, 58, 59, 60, 61] based on previous MHD simulations [62] to study the orbital evolution in turbulent disks. Specifically, we add a fluctuating potential $\Phi_{\rm turb}$ to the momentum equation consisting of 50 stochastic modes at each time step:

[TABLE]

where $\gamma$ is a dimensionless characteristic amplitude of turbulence. Each mode, denoted as $\Lambda_{k}$, is given by

[TABLE]

which is associated with a wavenumber $n_{k}$ drawn from a logarithmically uniform distribution between $n_{k}=1$ and the maximum value $n_{\rm max}$, corresponding to the azimuthal grid scale. The initial radial position $r_{k}$ and azimuthal angle $\phi_{k}$ of each mode are sampled from a uniform distribution. The radial extent of each mode is defined as $\sigma_{k}=\pi r_{k}/4n_{k}$. The modes are activated at time $t_{0,k}$ and persist for a duration $\Delta t_{k}=0.2\pi r_{k}/n_{k}c_{s}$, where $c_{s}$ denotes the local sound speed. $\Omega_{k}$ represents the Keplerian frequency at $r_{k}$, $\tilde{t}_{k}=t-t_{0,k}$, and $\xi_{k}$ is a dimensionless constant sampled from a Gaussian distribution with unit width. Following Baruteau and Lin [57], the parameter choices for this turbulence driver emulate a Kolmogorov cascade power spectrum, maintaining a Scale Law of $n_{k}^{-5/3}$.

Based on the velocity fluctuations induced by this turbulence driver, an effective time-averaged Reynolds stress parameter $\langle\alpha_{\rm R}\rangle$ (hereafter referred to as the dimensionless viscosity parameter $\alpha$) in the vicinity of the sBH can be calibrated as [57]

[TABLE]

This corresponds to an effective kinematic viscosity of $\nu=\alpha H^{2}\Omega$, which is consistent with previous calibrations [59, 60, 61].

Our simulations are performed in a 2D $\left(r,\phi\right)$ coordinate system, with a computational domain ranging radially from 0.4 $r_{0}$ to 3.2 $r_{0}$ and azimuthally from 0 to $2\pi$. The domain is resolved by 512 logarithmic-spaced grid cells in the radial direction and 1536 uniformly spaced cells in the azimuthal direction. We apply wave-damping radial boundary conditions [63] to both the inner edge ($[0.4-0.5]r_{0}$) and the outer ($[2.8-3.2]r_{0}$) edge to avoid wave reflection. Tests with different numerical resolutions have been performed to confirm the convergence of our simulation results. To avoid numerical artifacts caused by damping of turbulence at the boundaries, we limit the turbulent region to a radial range between 0.5 $r_{0}$ and 2.8 $r_{0}$, instead of extending turbulence across the entire disk.

In all simulations presented in this work, the sBH is initially located at an orbital radius of $r=1.0~r_{0}$ with a mass ratio of sBH to MBH of $q\equiv m/M_{\bullet}$. We fully account for the orbital evolution of the sBH because of the gravitational interaction with the disk. We will explore the orbital evolution of the sBH with different turbulence strengths $\gamma$, mass ratios $q$, disk surface densities $\Sigma_{0}$, and different disk aspect ratios $h_{0}$.

II.2.2.2 Orbital Perturbation with Turbulent Forces.

Before presenting the results of the hydrodynamical simulation, we first discuss how turbulent disk forces induce orbital perturbations. Considering the sBH perturbed by the (turbulent) fluctuating gravitational force from the disk, the equation of motion is given by [64]

[TABLE]

with $f_{r}$ being the specific force in the radial direction and $t_{\phi}$ the specific torque in the azimuthal direction. We consider the case where sBH is weakly perturbed from a circular orbit.

[TABLE]

with $\Omega_{0}^{2}=M_{\bullet}/r^{3}_{0}$. The equation of motion can be simplified as

[TABLE]

In the frequency domain, we have the replacement rules $\partial_{t}\rightarrow-i\omega,\partial^{2}_{t}\rightarrow-\omega^{2}$, so the equation of motion further becomes

[TABLE]

The solution is given by

[TABLE]

The induced eccentricity, which is described by the magnitude of radial fluctuations, can be estimated as

[TABLE]

This equation tends to yield divergent results if $S_{f_{r}}(\Omega_{0})$ or $S_{t_{\phi}}(\Omega_{0})$ is nonzero. To regularize the equilibrium eccentricity, we recognize that the co-orbital Lindblad resonances and density waves in the presence of turbulence tend to damp the eccentricity, say at a rate $\gamma_{e}=1/\tau_{e}\ll\omega$. Therefore, we effectively modify the first equation of Eq. (II.2.2.2) by

[TABLE]

so that

[TABLE]

The pole or resonance part of this integral is given by (assuming $\gamma_{e}\ll\Omega_{0}$)

[TABLE]

The equilibrium eccentricity can be estimated by combining the smooth component of Eq. (30) and the pole component (Eq. 31) of the power spectrum. In other words, the underlying eccentricity damping rate, $\gamma_{e}$, can be constrained based on the simulated eccentricity (or the orbital perturbation $\delta r/r$) of the embedded object and the turbulent force acting on the embedded object. In most cases, we find that the radial forces dominate over the azimuthal forces; thus it is the first term that mainly contributes to the eccentricity evolution of the EMRI.

II.2.2.3 Hydrodynamical Evolution of EMRI.

We first present the simulation results for an sBH embedded in a disk with a typical model parameters. The mass ratio of the sBH to the MBH is set to $q=m/M_{\bullet}=5\times 10^{-6}$, where the mass of the sBH is $m=20\ M_{\odot}$ and the mass of the MBH is $M_{\bullet}=4\times 10^{6}\ M_{\odot}$. The other disk parameters are $\Sigma_{0}=10^{-4}\ M_{\bullet}/r_{0}^{2}$ and $h_{0}=0.03$ at $r_{0}$. We consider a turbulent disk with a turbulence strength of $\gamma=1.6\times 10^{-3}$, which corresponds to the viscosity parameter $\alpha=0.1$.

A snapshot of the disk surface density at 4000 orbits (measured in units of the orbital period $P_{0}$ at $r_{0}$) is shown in Fig. 9. Strong density fluctuations are present throughout the disk as a result of the injection of external turbulence. We have confirmed that the density power spectrum follows a power-law scaling consistent with $n_{k}^{-5/3}$. The low mass ratio results in a weak density wave around the embedded object.

Under this circumstance, the inward migration rate is approximately 2-3 times faster than the classical type I migration rate as shown in the upper panel of Fig. 10, possibly due to turbulent density structures. During inward migration of the EMRI, an important feature is that the eccentricity can be excited as shown in Fig. 10. As shown in the lower panel of Fig. 10, the standard deviation of the radial oscillation from the mean radial motion $\delta r/r$, which is a measurement of the mean eccentricity of the EMRI $e$, is $\sim 0.01$. We have confirmed that such a radial oscillation (or the orbital eccentricity) does not evolve with time significantly after a few thousand orbits, indicating a convergence in time. The eccentricity excitation arises from the turbulent density structures, which generate the stochastic radial and azimuthal acceleration of the EMRI, as shown in Fig. 11. Based on the radial and azimuthal acceleration profile for the EMRI, we can obtain an effective eccentricity damping rate of the turbulent disk $\gamma_{e}\simeq 10^{-4}\Omega_{0}$ for this disk model using Eqs. (30) and (31). This damping rate for this particular model is on the same order of magnitude as the type I damping rate of the orbital eccentricity in the subthermal regime $\sim qh_{0}^{-4}\Sigma_{0}r_{0}^{2}/M_{\bullet}\Omega_{0}\simeq 5\times 10^{-4}\Omega_{0}$ [e.g., 65, 44, 66]. However, as we will show later, the parameter dependence of this damping rate differs from that of the classical Type I eccentricity damping rate.

II.2.2.4 Scaling Relation for Turbulent Disks.

The evolution of EMRI eccentricity can depend on various disk parameters, such as the turbulence strength $\gamma$, the disk aspect ratio $h_{0}$, the disk surface density $\Sigma_{0}$ and the EMRI mass ratio $q$. To comprehensively characterize the turbulence-induced eccentricity under general conditions, we conduct an extensive parameter survey across a range of model parameters.

The dependence of $\delta r$, the turbulence forces ($\delta F_{r}$ and $\delta\Gamma$), and the damping rate $\gamma_{e}$ on the disk turbulence strength $\gamma$ and the mass ratio $q$ are shown in Fig. 12. We have fixed the disk aspect ratio of $h_{0}=0.03$, and a surface density $\Sigma_{0}=10^{-4}\ M_{\bullet}/r_{0}^{2}$. The pow-law scaling for different quantities is shown as a dashed or dotted line in each panel. It can be seen that the relations $\delta F_{r}\propto q\gamma$, $\delta\Gamma\propto q\gamma$, $\gamma_{e}\propto\gamma^{-1}$ can roughly match the parameter dependence. The linear dependence of the turbulence force acting on the sBH results in a constant specific force ($f_{r}\equiv\delta F_{r}/q$, and $\tau_{\phi}=\delta\Gamma/q$) for the sBH. This behavior can be attributed to the fact that the turbulent potential is independent of the mass ratio $q$. The linear dependence of $\gamma$ is directly related to the scaling of the turbulence potential of Eq. (20). The scaling of the eccentricity damping rate $\gamma_{e}$ in the turbulence strength suggests a suppressed eccentricity damping for the stronger turbulence, which is probably due to the weaker coherent structures and/or Lindblad resonances in strongly turbulent environments. For the parameter explored here, we have $f_{r}\gg\tau_{\phi}/r_{0}$. Therefore, the radial fluctuation follows $\delta r\propto(f_{r}^{2}/\gamma_{e})^{1/2}\propto\gamma^{1.5}q^{0}$, which is roughly consistent with the scaling results shown in the upper left panel of Fig. 12. The deviation from linear scaling, especially the flattening of the scaling, at large $\gamma$ may arise from strong turbulence-induced disturbances in the disk, leading to a nonlinear response. To account for this deviation in the $\delta F_{r}-\gamma$ relation, we also tested an alternative scaling, $f_{r}\propto\gamma^{1/2}$, which appears to better capture the flattening of the trend at higher turbulence levels. The resulting scaling for $\delta r$, $\delta r\propto\gamma$, is shown as the dashed line in the upper left panel of Fig. 12. In the following, we will adopt both scalings in long-term N-body simulations to predict the orbital evolution of the EMRI in the vicinity of the MBH.

Now we further explore the scaling relation with respect to the disk aspect ratio. Here we have fixed a disk surface density $\Sigma_{0}=10^{-4}\ M_{\bullet}/r_{0}^{2}$, a turbulence strength $\gamma=10^{-4}$, and the EMRI mass ratio $q=5\times 10^{-6}$. The simulation results for different $h_{0}$ are shown in Fig. 13. We can see that $\delta F_{r}\propto h_{0}^{-1}$. There is a weak dependence of $\delta\Gamma$ on $h_{0}$, so we treat $\delta\Gamma$ as a constant with respect to $h_{0}$. For the damping rate $\gamma_{e}$, we identify a scaling relation $\gamma_{e}\propto h_{0}^{-4}$, which coincidentally resembles the classical type I scaling [e.g., 65, 44, 66]. Given these scalings, the radial fluctuation follows $\delta r$ as $\delta r\propto(f_{r}^{2}/\gamma_{e})^{1/2}\propto h_{0}$, which is consistent with the simulation results shown in the upper left panel of Fig. 13.

We also present the scaling of the time-averaged background force magnitude $|\langle F_{r}\rangle|$ and $|\langle\Gamma\rangle|$ in the right panels of Fig. 13. The “quasi-steady” (time-averaged) component of the radial force is comparable to or even larger than the stochastic one and shows almost no dependence on the disk aspect ratio $h_{0}$. In contrast, the “quasi-steady” component of the torque is generally weaker than the stochastic component and exhibits a power-law decay with $h_{0}^{-2}$, consistent with the scaling of the type I migration torque [67, 68].

Finally, we perform a set of additional simulations with different disk masses by varying $\Sigma_{0}$ from $10^{-7}M_{\bullet}/r_{0}^{2}$ to $10^{-4}M_{\bullet}/r_{0}^{2}$. The other model parameters are chosen as the turbulence strength $\gamma=10^{-4}$, a disk aspect ratio of $h_{0}=0.03$, and a mass ratio for EMRI $q=2\times 10^{-6}$. The simulation results are shown in Fig. 14. We find that all the quantities $\delta F_{r}$, $\delta\Gamma$, and $\delta r$ are linearly proportional to $\Sigma_{0}$. This implies that the damping rate $\gamma_{e}$ should be independent of $\Sigma_{0}$ over the three orders of magnitude explored, which is remarkably different from the classical type I scaling ($\propto\Sigma_{0}$) for the orbital eccentricity. The distinct scaling behavior of the turbulence-induced $\gamma_{e}$ highlights the need for further dedicated studies to uncover its physical origin.

Now we can summarize the full scaling of stochastic forces. For the stochastic radial force,

[TABLE]

where $s=0.5$ or $1.0$ based on the scaling of $\delta F_{r}$ on the turbulence strength $\gamma$ shown in Fig. 12. For the stochastic torque,

[TABLE]

and the damping rate is

[TABLE]

When the stochastic radial force dominates (e.g., for smaller $h_{0}$) the orbital eccentricity evolution, we have the equilibrium orbital perturbation, namely eccentricity as $\delta r/r\propto\Sigma_{0}h_{0}\gamma^{s+0.5}r^{1.25}$. Therefore, we can obtain

[TABLE]

Otherwise, for the $\tau_{\phi}$ dominating case (e.g., for larger $h_{0}$), $\delta r/r\propto\Sigma_{0}h_{0}^{2}\gamma^{1.5}r^{1.25}$, we have

[TABLE]

Note that these scalings depend solely on local disk conditions, neither on global radial profiles nor specific AGN disk models.

The above discussion assumes quasi-equilibrium in the evolution of eccentricities at each radius $r$, driven by the stochastic force. As the sBH inspirals towards the MBH, the gravitational wave radiation becomes increasingly significant. Consequently, there exists a critical radius $r_{\rm peak}$, within which the eccentricity damping timescale $1/\gamma_{\rm e,gw}$ by gravitational wave radiation becomes shorter than the turbulent disk equilibrium timescale $1/\gamma_{e}$ (also see Fig. 15). At radii smaller than this critical value (that is, $r\lesssim r_{\rm peak}$), the eccentricity evolution is predominantly governed by gravitational wave radiation, following a radius scaling of $r^{19/12}$ [41].

II.2.2.5 Longterm Evolution and Implications.

We implement the stochastic radial $f_{r}$ and torque $\tau_{\phi}$ with the damping rate $\gamma_{e}$ following the above scaling relation into N-body code Rebound [46]. We use the WHFast integrator [69] in Rebound to perform the long-term calculation. Migration and turbulence forces are implemented in Reboundx [70] to study the long-term evolution of EMRI. The “quasi-steady” component of the turbulence torque drives the inward migration of the sBH, whereas the stochastic components of the radial and azimuthal forces can excite the eccentricity of the EMRI. The eccentricity and semi-major axis damping due to gravitational wave radiation are also included in the N-body calculation.

As an illustrative example, we first consider a $\beta$-disk model, in which the disk surface density is higher in the inner region ($r\lesssim 100\ M_{\bullet}$) compared to a $\alpha$-disk model with the same disk parameter (refer to Fig. 3). As a result, the turbulence-driven eccentricity in the inner disk is higher. We choose the model parameter as $\alpha=0.1$, $\dot{M}_{\bullet}=0.5\dot{M}_{\rm Edd}$, $M_{\bullet}=4\times 10^{6}M_{\odot}$ and $q=5\times 10^{-6}$. The disk profiles are shown in Fig. 3. The scaling index $s$ for the radial turbulence force is set to $s=1.0$. Since the turbulence force is stochastic in nature, we simulate 100 Monte Carlo realizations with the same model parameters to explore the statistics properties of the EMRI.

The results for these N-body simulations are shown in Fig. 15. As the EMRI migrates inward, its orbital eccentricity increases as shown in the upper panel of Fig. 15, although there are some fluctuations due to the implementation of stochastic forces. The increase in eccentricity is attributed to the radial dependence of the turbulence force and damping rate, such that the equilibrium eccentricity varies with orbital radius. Based on the scaling shown in the previous Section, we find $e_{\rm sBH}\propto r^{-1.85}$, shown as the dashed line in the upper panel of Fig. 15. This radial dependence is in good agreement with our N-body simulation results. Within a radius of $r\sim 30\ M_{\bullet}$, the gravitational wave radiation becomes dominant and the eccentricity of the sBH decays with radius following $r^{19/12}$, as expected.

The lower panel of Fig. 15 shows the statistical distribution of the orbital eccentricity during the last inspiral stage of EMRI, for example, at $10\ M_{\bullet}$. We can see that the predicated EMRI eccentricity at $r=10\ M_{\bullet}$ is around $e_{\rm sBH}\simeq 2\times 10^{-4}$, which is detectable for future LISA/Taiji/TianQin observations [9, 10].

We have also explored the case where the power law index for the radial turbulence force is $s=0.5$. In this case, the radial dependence of the eccentricity is $e_{\rm sBH}\propto r^{-1.35}$ based on the scaling of the radial turbulence force. Due to the much steeper radial dependence of the azimuthal turbulence forces, they can also contribute comparably to the orbital eccentricity around $10\ M_{\bullet}$. As a result, the equilibrium eccentricity reaches a level comparable to $10^{-4}$ at $r\simeq 10\ M_{\bullet}$.

For an $\alpha$-disk model with the same model parameters, the disk scale height remains almost the same; however, the disk surface density scales as $\Sigma\simeq 10^{5}(r/100\ M_{\bullet})^{1.5}\ {\rm g\ cm^{-2}}$ in the inner region of the AGN disk ($r\lesssim 100M_{\bullet}$). This results in a significantly lower disk density at $10\ M_{\bullet}$, resulting in a much lower eccentricity of the EMRI $e_{\rm sBH}\lesssim 10^{-7}$ induced by disk turbulence. Such a low eccentricity lies below the detection limit of the LISA/Taiji/TianQin and could be undetectable.

To further explore the parameter regime in which the EMRI eccentricity is detectable by future observations, we show the EMRI eccentricity as a function of disk surface density $\Sigma$ and disk aspect ratio $h$ at $10\ M_{\bullet}$ in Fig. 16. Here we have adopted the typical disk viscosity of $\alpha=0.1$ (upper panel) and $\alpha=0.02$ (lower panel). The radial profile of a $\beta$-disk model is chosen to determine the radius at which $\gamma_{\rm e,gw}\simeq\gamma_{e}$. The eccentricity of the EMRI increases with both the density of the disk surface and the aspect ratio of the disk. For a disk aspect ratio of $h\simeq 0.1$ in the inner region of the AGN disk, the eccentricity of the EMRI becomes detectable when the density of the surface of the disk exceeds $\Sigma\gtrsim 5\times 10^{5}\ {\rm g\ cm^{-2}}$. Such a high surface density at $10\ M_{\bullet}$ is more likely to occur in a $\beta$-disk model than in a $\alpha$-disk model. A higher disk viscosity $\alpha$ (or stronger turbulence) leads to a greater eccentricity of the EMRI at given $\Sigma$ and $h$. This results from the positive scaling of the equilibrium eccentricity with turbulence strength, as shown in Eq. (35-36). However, it should be noted that for a given $\dot{M}_{\bullet}$, the surface density of a $\beta$-disk decreases with increasing $\alpha$ (see Fig. 3), leading to similar EMRI eccentricities in $\beta$-disk models with different viscosities $\alpha$.

In Fig. 16, we further map the locations of different $\beta$-disk models in the eccentricity contour. The trajectories for $\alpha=0.1$ (upper panel) and $\alpha=0.02$ (lower panel), corresponding to different values of $\dot{M}_{\bullet}$ in the $\beta$-disk model, are shown as colored lines. The three stars embedded in each line indicate the models with $\dot{M}_{\bullet}=0.02,0.1,0.5\dot{M}_{\rm Edd}$, respectively. As expected, a higher $\dot{M}_{\bullet}$ leads to a higher EMRI eccentricity, and the dependence of the eccentricity on the disk viscosity $\alpha$ at a fixed $\dot{M}_{\bullet}$ is weaker. These results suggest that precise measurements of EMRI eccentricities may allow us to constrain the AGN disk parameters and the underlying accretion model.

We note that Copparoni et al. [71] studied the evolution of EMRI using stochastic disk torques based on hydrodynamic simulations of [72]. However, these hydrodynamic simulations do not implement any MHD turbulence, which is different from the turbulent disk model considered in our work. The stochastic disk forces therein likely originate from the disk fluctuations induced by the perturber and may not share the same properties as the stochastic torques and forces arising from the intrinsic turbulent AGN disks.

II.2.2.6 Impact on waveform.

As sBH migrates inward in the accretion disk, the migration force generally modifies the long-term orbital evolution and, consequently, the gravitational waveform. Since the leading-order phase term of an EMRI scales as $1/q$, it is expected that a phase modulation of order $\mathcal{O}(1)$ radians is detectable. Therefore, a migration force that is at least a fraction $\mathcal{O}(q)$ of the gravitational radiation reaction should, in principle, be observable. This roughly translates to the requirement that the ratio between the gravitational wave damping timescale and the migration timescale be greater than or equal to $10^{-5}$-$10^{-4}$.

For one typical Monte Carlo realization adopted in Fig. 15, we show the semi-major axis evolution rate $1/\tau_{\rm a}$ ($\equiv\dot{r}/r$) of EMRI in Fig. 17. It is clearly seen that disk-driven migration dominates the evolution of the semi-major axis in the outer region, while gravitational wave radiation takes over in the inner region. Theoretically, the disk-driven migration rate of the sBH is

[TABLE]

where $\Omega_{\bullet}$ is the Keplerian frequency at $r=M_{\bullet}$. For a $\beta$ disk with $\dot{M}_{\bullet}=0.5\dot{M}_{\rm Edd}$, where $h_{0}\propto r^{-1}$ and $\Sigma_{0}\propto r^{-0.6}$, $\alpha=0.1$, we then have $\dot{r}/r\simeq 3\times 10^{-11}\Omega_{\bullet}(r/100M_{\bullet})^{1.9}$, which is fully consistent with the simulation results. Gravitational wave radiation ($\propto r^{-4}$) dominates at smaller radii. The transition between the two regimes occurs around $r\simeq 50\ M_{\bullet}$. Although the ratio of the disk-driven migration timescale $\tau_{\rm mig}$ to the gravitational wave inspiral timescale $\tau_{\rm gw}$ at $r\simeq 10\ M_{\bullet}$ is of the order of $\sim 10^{4}$, such a difference can still leave observable imprints on the EMRI waveform.

To explore the parameter space in which the impact of disk-driven migration on the EMRI waveform becomes detectable, we show the ratio $\tau_{\rm mig}/\tau_{\rm gw}$ in Fig. 18 for different values of disk surface density $\Sigma$ and aspect ratio $h$ at $10\ M_{\bullet}$. According to the type I and type II migration formulas, a lower disk aspect ratio $h$ and a higher surface density $\Sigma$ are expected to result in a smaller $\tau_{\rm mig}/\tau_{\rm gw}$, thus producing a stronger impact on the waveform, as illustrated in Fig. 18. This timescale ratio is generally independent of disk viscosity $\alpha$, as shown in Fig. 18, except in the very low regime $h$ where the opening of the gap is likely to occur, which is consistent with the type I/II migration rate. It should be noted that this mapping is independent of the disk’s radial structure or evolution history, since we are only concerned with the parameters of the disk at $10\ M_{\bullet}$.

Thus, we can further map the different AGN disk models in the contour $\tau_{\rm mig}/\tau_{\rm gw}$. Similarly to Fig. 16, we show the trajectories of both the $\alpha$ model and the $\beta$ model with different $\dot{M}_{\bullet}$ in Fig. 18. These trajectories depend on disk viscosity $\alpha$, and therefore occupy different locations in the contour $\tau_{\rm mig}/\tau_{\rm gw}$ for different values of $\alpha$. Unlike eccentricity detection by EMRI, a lower $\dot{M}_{\bullet}$ leads to a smaller $\tau_{\rm mig}/\tau_{\rm gw}$, and hence a stronger impact on the waveform. The three stars embedded in each colored line (representing the $\alpha$ model on the left and the $\beta$ model on the right) correspond to $\dot{M}_{\bullet}=0.02,\ 0.1,\ 0.5\dot{M}_{\rm Edd}$, respectively. For a disk viscosity of $\alpha=0.02$, the gap opening process could occur only for $\dot{M}_{\bullet}\lesssim 0.01\dot{M}_{\rm Edd}$ (that is, $h\lesssim 0.007$), in which case the migration proceeds at the type II rate [73, 74, 75]. The AGN disk with an even lower $\dot{M}_{\bullet}\lesssim 0.01\dot{M}_{\rm Edd}$ would accrete via an advection-dominated accretion flow [76], different from the $\alpha/\beta$ disk considered here. A higher $\alpha$ (e.g. $\alpha=0.1$ in the upper panel) would suppress the formation of the gap and, therefore, results in type I migration across the entire range of model parameters considered here. The transition from type I to type II migration can be seen at the lowest $h$ in the lower panel of Fig. 18. We find that for almost all values of $\dot{M}_{\bullet}$ explored in the $\beta$-disk model, the disk-driven migration has a strong enough effect to produce an observable waveform impact. For the $\alpha$-disk model, the waveform impact remains observable for $\dot{M}_{\bullet}\sim 0.02\dot{M}_{\rm Edd}$.

In addition to the migration effect, the stochastic force may also influence the secular evolution of the EMRI phase. A subtle issue arises from the fact that the corresponding phase modulation may contain a stochastic component, which may not be compatible with traditional matched filter analyses that rely on deterministic waveforms. It would be interesting to explore how turbulence-induced randomness can be incorporated into waveform modeling, and whether it can be measured with sufficient precision to constrain the properties of the turbulence. A potentially relevant discussion on the inference of stochastic forces can be found in [77].

III Inclination

Another important observable characterizing various EMRI formation channels is the inclination angle between the EMRI orbital angular momentum vector and the spin axis of the central MBH. Note that in Kerr spacetime, a geodesic generally does not conserve a three-dimensional angular momentum vector. However, the inclination angle $\iota$ can be defined using the relation $\cos\iota:=L_{z}/\sqrt{L_{z}^{2}+\mathcal{Q}}$, where $L_{z}$ is the component of the angular momentum along the spin axis, and $\mathcal{Q}$ is the Carter constant. This definition reduces to the conventional one in Newtonian mechanics for wide orbits. Moreover, it is straightforward to show that, during a generic inspiral process of an eccentric, inclined EMRI, $\iota$ generally changes by no more than a few degrees [78]. In other words, it is a rather robust observable that is insensitive to the time of measurement within the LISA frequency band.

The distribution of $\iota$ differ significantly between wet and dry formation channels. As discussed in Sec. III.1, the accretion to the MBH determines the initial orbital inclination of the wet EMRI, and both the Bardeen-Peterson timescale and the lifetime of the disk determine the final distribution. For dry EMRIs, it turns out that the critical factor for the distribution of $\iota$ arises from the critical angular momentum of the last stable orbit. This generally leads to an asymmetric distribution, as detailed in Sec. III.2.

III.1 Wet EMRI

Two major scenarios are commonly discussed to describe the accretion episodes of MBHs: coherent accretion and chaotic accretion. In coherent accretion, the orbital plane of the incoming gas remains roughly constant across episodes, so the MBH spin tends to align with the angular momentum of the disk. It is also natural to expect that the MBH spin becomes near-extremal after several e-folds of mass growth through coherent accretion. One subtlety is that MBH mergers, following galaxy mergers, can significantly alter the mass and spin of the remnant black hole. However, since such mergers are infrequent and the typical time between them is much longer than the Bardeen-Petterson (BP) alignment timescale, we do not consider this effect. In the chaotic scenario, the direction of accretion in each episode is assumed to be random. As a result, the MBH spin is likely misaligned with the normal of the accretion disk at the beginning of each episode, with the inclination gradually evolving due to the BP effect.

As a wet EMRI is generally embedded in the accretion disk during its migration, its initial orbital plane aligns with the disk plane. Therefore, in the coherent accretion scenario, we expect the inclination distribution $P(\cos\iota)$ to peak around $\cos\iota=1$, or $P(\cos\iota)\approx\delta(\cos\iota-1)$.

In the case of chaotic accretion, the angular momentum of the accretion disk is initially misaligned with the spin of the central MBH. This misalignment is gradually reduced due to Lense-Thirring (LT) precession [79], which arises from the frame-dragging effect in Kerr spacetime. The inner parts of the disk are forced to align with the spin of the MBH, with the characteristic warp radius given by [27, 80, 81]

[TABLE]

where $h$ is the aspect ratio of the disk and $\chi$ is the dimensionless spin of the central MBH. The warpped disk will exert a torque onto the center MBH and align or anti-align the BH spin with the angular momentum of the outer disk. The characteristic alignment or anti-alignment timescale over which this process occurs is known as the Bardeen-Petterson timescale $\tau_{\text{BP}}$[27]:

[TABLE]

which depends on the MBH mass $M_{\bullet}$, the dimensionless spin parameter $\chi$, the viscosity parameter $\alpha$, the MBH luminosity $L$ ($L_{\rm Edd}$ is the Eddington luminosity), and the radiative efficiency $\epsilon$.

Wet EMRIs embedded within the disk remain in the disk during their migration at large radii. The LT torque exerted on the EMRI tends to drive its orbit out of the disk plane, whereas the disk torque acts to realign off-plane orbits back into the plane. The LT torque dominates for $r\leq R_{\rm dec}$, where the decoupling radius $R_{\rm dec}$ -which is smaller than the radius where the LT torque overtakes the disk density wave torque - is given by [81]:

[TABLE]

According to the analysis in [78], even when a wet EMRI evolves within the decoupling radius, the inclination angle relative to the MBH spin changes by no more than a few degrees. Therefore, the EMRI inclination still closely traces the disk inclination at larger radii.

Assuming that the MBH is completely uncorrelated with the accretion disk at the beginning of each accretion episode, the distribution of the initial inclination angle $\iota_{0}\in(0,\pi)$ should be

[TABLE]

The evolution of the inclination angle between the center MBH spin and the orbital angular momentum of the outer disk depends on the ratio of the angular momentum $J_{\rm d}/J_{\bullet}$ [82], where $J_{\rm d}$ and $J_{\bullet}$ are the angular momentum of the disk and the angular momentum of the spin of MBH, respectively. If $\cos\iota_{0}>-J_{\rm d}/2J_{\bullet}$, the BH will eventually align with the outer disk, and if $\cos\iota_{0}<-J_{\rm d}/2J_{\bullet}$, the BH will finally anti-align with the outer disk.

In the $\alpha$-disk model, the angular momentum of the outer disk is significantly greater than that of the MBH, expressed as [26]

[TABLE]

in which $\epsilon$ is the radiative efficiency, $\dot{M}_{\bullet}$ is the MBH accretion rate.

For the typical MBH in this paper ( $M_{\bullet}=10^{5}-10^{7}M_{\odot}$), $J_{\rm d}\gtrsim 2J_{\bullet}$, which means that the BH will finally align with the outer disk for any initial inclination angle $\iota_{0}$. The time-dependent evolution of the inclination angle can be described as [82]

[TABLE]

where the $\tau$ is a coefficient and should have a relation with the alignment timescale $\tau_{\rm BP}$. To counteract the influence of the angular momentum of the disk and comparing the small angle case, we use $\tau=(1+J_{\rm d}/J_{\bullet})\tau_{\rm BP}$.

The solution of above equation for $J_{\rm d}/J_{\bullet}\gtrsim 2$ can be integrated directly to an implicit form

[TABLE]

where $C$ is the integral constant and can be obtained by setting $t=0$.

Starting from the initial distribution of the inclination angle described by Eq. 41, the time-dependent evolution of the inclination distribution of the disk can be obtained from

[TABLE]

which can be solved numerically.

In order to compute the distribution of $\iota$ of wet EMRIs, we also need to account for the time-dependent formation rate of wet EMRIs. According to the analysis in [36], wet-EMRI formation rate $\Gamma_{\rm M_{\bullet},EMRI}(t)$ typically peaks at the beginning of AGN episodes, and gradually declines due to the insufficient supplies of sBHs from the outer boundary of the nuclear star cluster. Assuming the disk lifetime to be $T_{\rm disk}$, the distribution of $\iota$ for wet EMRIs is then given by

[TABLE]

where the total number of wet EMRIs formed within the lifetime is

[TABLE]

Here, $\Gamma_{\rm M_{\bullet},EMRI}$ denotes the EMRI formation rate as well, with the subscript ${\rm M_{\bullet}}$ indicating its dependence on the MBH mass.

Notice that the time-dependent wet-EMRI formation rate based on the above FP calculations explicitly depends on the mass of the MBH, as does the associated $\iota$ distribution in Eq. (46).

A sample distribution is illustrated in Fig. 19, in which case the MBH mass is assumed to be $10^{6}M_{\odot}$, the dimensionless spin is $a=0.8$, and the disk properties are set to $\alpha=0.1,L=0.1L_{\text{Edd}},\epsilon=0.1$. We consider three different disk lifetimes: $T_{\text{disk}}=1\times 10^{6},\ 5\times 10^{6},\ 1\times 10^{7}\text{yr}$. Notice that according to Eq. (III.1), the BP timescale is approximately $1.7\times 10^{6}\text{yr}$.

For inclined disks, the BP mechanism is generally in operation, and the final state depends on both the disk lifetime and the BP timescale. In Fig. 19, we observe that when the disk lifetime is shorter than the BP timescale ($T_{\text{disk}}<\tau_{\text{BP}}$), the distribution of $\cos\theta$ remains approximately flat, indicating little alignment. However, when the disk lifetime exceeds the BP timescale ($T_{\text{disk}}>\tau_{\text{BP}}$), the inclination angle distribution becomes tilted, with peak corresponding to alignment ($\iota\approx 0$) with the MBH spin. As the disk lifetime increases, the tilted distribution becomes sharper. Therefore, the inclination distribution of wet EMRIs encodes important information about $\tau_{\rm BP}/T_{\rm disk}$.

For actual observations, we need to average over the mass function of MBHs to obtain the total inclination distribution of wet EMRIs. Similar to the treatment in Sec. II.1, we vary the power-law index of the MBH mass function from $\tilde{\gamma}=-0.3$ to $\tilde{\gamma}=0.3$, i.e.

[TABLE]

where the normalization constant $N$ is

[TABLE]

The range of mass integration is from $10^{5}M_{\odot}$ to $10^{7}M_{\odot}$, which is relevant for space-borne gravitational wave observation.

The resulting distribution function is shown in Fig. 20 for two disk lifetimes: $T_{\text{disk}}=10^{6}\ \text{yr}$ (left panel) and $T_{\text{disk}}=5\times 10^{6}\ \text{yr}$ (right panel). In the case of a shorter disk lifetime, the deviation in the mass function significantly affects the inclination angle distribution. However, for a longer disk lifetime, the deviation in the mass function only affects the region where $\cos\iota\lesssim 0.5$, due to the alignment effect.

When the lifetime of the disk is several times longer than the BP timescale, the alignment process is almost complete, and the wet EMRIs are more aligned with the MBH spin, as illustrated in the right panel.

We should note that the angle of inclination distribution of wet EMRIs in the case of chaotic accretion could be affected by the relationship between $J_{\rm d}$ and $J_{\bullet}$. If $J_{\rm d}\ll J_{\bullet}$, the critical angle that distinguishes alignment from anti-alignment of the MBH spin and the disk is $\arccos(0)=\pi/2$, which means that initially prograde cases will align with the MBH spin, while initially retrograde cases will anti-align with it. Consequently, in the chaotic accretion scenario, we may observe a bimodal inclination-angle distribution of wet EMRIs. Therefore, the distribution of the angle of inclination of wet EMRIs can not only reveal the BP effect in AGN disks but also provide information on the relationship between $J_{\rm d}$ and $J_{\bullet}$.

III.2 Dry EMRI

Assuming a spherically symmetric nuclear star cluster, the formation rate of dry EMRIs still depends on their inclination with respect to the spin of the MBH [83]. In general, for aligned orbits the angular momentum of the Last Stable Orbits (LSOs) is lower than that of misaligned orbits. Fig. 1 shows the loss cone boundary with $J_{\rm lc}=4M_{\bullet}$, but orbits with smaller inclination angles should have a lower loss cone boundary (in the $R$ direction) due to the lower LSO angular momentum. This implies that the critical energy $E_{\text{cr}}$ is lower, with a correspondingly larger integrated flux at the boundary $t_{\rm GW}=t_{J,\rm sca}$. Therefore, higher EMRI rates are expected for prograde orbits than for retrograde orbits.

The LSO radius and angular momentum for Kerr geodesics can be obtained following the method in [83], as a function of the black hole spin $a$ and the inclination angle $\iota$. They generally also depend on orbital eccentricity, but in practice, since the EMRIs formed in their early stages are highly eccentric, we essentially adopt the LSO values in the $e=1$ limit. In Fig. 21 we present a series of FP calculations following similar procedures in Sec. II.1, by choosing $\chi=0.8$ and the $\iota$-dependent LSO angular momentum $J_{\rm lc}(\cos\iota)$. We evolve all systems to 5 Gyrs. For more massive MBHs, the asymmetry of the “tilted” distributions becomes significantly larger than for lower-mass MBHs, which is because different MBH systems reach the equilibrium state at different times. During the evolution stage, the sBH component will gradually flow to the inner region from the outer region, reducing the differences in EMRI flux at the loss-cone boundary over time. In lower-mass MBH systems, equilibrium is reached relatively quickly, after which the systems evolve in a quasi-equilibrium state. Consequently, the EMRI flux at the loss-cone boundary shows smaller differences compared to a non-equilibrium state. In contrast, higher-mass MBH systems take much longer to reach equilibrium (e.g. $\mathcal{O}(1{\rm Gyrs})$ for heavier MBH and $\mathcal{O}(0.1{\rm Gyrs})$ for lighter MBH). As a result, their EMRI flux at the loss-cone boundary exhibits larger variations across different energies, directly affecting the EMRI rate. In fact, orbits with different inclination angles have different critical energies $E_{\rm cr}$, as discussed in the previous section. Therefore, the combined effect of evolution time and inclination angle results in varying degrees of tilt asymmetry in the inclination angle distribution for different MBH systems.

Similarly to the discussion in Sect. III.1, we also consider the average distribution of $\iota$ for dry EMRIs, with the power-law index of the MBH mass function ranging from $\tilde{\gamma}=-0.3$ to $\tilde{\gamma}=0.3$. The result is shown in Fig. 22. We find a stronger dependence of $\gamma$ at angles around $\iota\approx\pi/2$ and $\iota\approx 0$, although the asymmetry dependence on the $\cos\iota$ is still significant.

In reality, the spins of MBHs may vary between different systems, but EMRI observations should be able to measure them with great precision, ie $\mathcal{O}(10^{-5})$ [9]. The above analysis can be modified to account for the spin distribution in order to predict the $\iota$ dependence of dry EMRIs, which can then be compared with the observed population of events from LISA, Taiji, and TianQin [9, 10]. A possible source of deviation, even when considering uncertainties in $\tilde{\gamma}$, may arise from intrinsic asymmetry in the nuclear star cluster. For example, if most clusters tend to have stellar disks aligned with the spin of the MBH, the relative ratio between prograde and retrograde EMRIs would be further enhanced. In this sense, the prediction shown in Fig. 22, along with its spin-modified variants, can serve as baseline benchmarks to investigate the asymmetries of nuclear clusters.

From the analysis of the inclination angle distributions of dry and wet EMRIs, we find that both types exhibit tilted distributions due to different physical origins. For short disk lifetimes, the tilt amplitudes are comparable; however, for long disk lifetime, the asymmetry in wet EMRIs becomes significantly larger than in dry EMRIs.

IV sBH Mass

The population study of BBHs observed by the Advanced LIGO-Virgo-KAGRA (LVK) Collaboration has been one of the central topics in gravitational wave astronomy. The mass function of black holes within these binaries encodes critical information about the underlying formation channels and the physical properties of the formation mechanisms. Moreover, the mass function of sBHs in EMRI systems should be closely related to the mass spectrum of isolated sBHs within the nuclear cluster. In this section, we discuss the information that can be extracted from the comparison of the mass functions of stellar BBHs, dry EMRIs, and wet EMRIs.

The mass spectrum of isolated sBHs remains poorly understood. Aside from Galactic X-ray binaries, the only direct observational constraints on the mass distribution of sBHs come from gravitational wave detections by the LVK Collaboration. If the pairing mechanism that forms BBHs can be identified, it may be possible to infer the mass distribution of isolated single black holes from the observed primary mass distribution of BBHs. In [84], the authors investigated a range of possible single black hole formation channels: pulsational pair instability supernovae, chemically homogeneous evolution, Population III stars, globular cluster dynamics, and hierarchical mergers, to determine whether any of them can reproduce the observed BBH mass spectrum through pairing. They concluded that no single channel can fully account for the population properties currently observed by LVK.

Rather than assuming a specific formation channel, in this paper we adopt a phenomenological approach. We model the mass distribution of single black holes and apply a pairing function to construct binaries. By comparing the resulting primary mass distribution with LVK observations, we infer the underlying single-BH mass spectrum corresponding to the chosen parameterized pairing function. This inferred mass function of isolated black holes is then used as a benchmark in Sec.IV.2 and Sec.IV.1, to demonstrate how the sBH mass functions in wet and dry EMRIs can be predicted.

We assume a phenomenological “Power law + Gaussian peak” model for the distribution of primary masses in LVK BBHs, following [1, 2]

[TABLE]

in which $\mathcal{P}$ is the normalized power law distribution defined between $m_{\rm min}=4.59M_{\odot}$ and $m_{\rm max}=86.22M_{\odot}$ (from GWTC-2 [1]) with power-law index $\alpha_{\rm p}$, $\mathcal{G}$ is the normalized Gaussian distribution with mean $\mu_{m}$ and standard deviation $\sigma_{m}$, $N$ is the normalized factor, $S$ is the smoothing function that increases from zero to one over the interval $(m_{\rm min},\ m_{\rm min}+\delta_{m})$,

[TABLE]

with

[TABLE]

where the parameters are $\alpha_{\rm p}=3.5,\mu_{m}=34,\ \sigma_{m}=5.69,\ \lambda_{\rm peak}=0.038,\ \delta_{m}=4.82$ as consistent in [1, 2].

In addition, although the single/isolated BH mass distribution is also unknown, we apply a similar parameterization for the mass function with the “Power law + Gaussian peak” model:

[TABLE]

and with a parametrized pairing function as [84]

[TABLE]

Here $q$ is the mass ratio $m_{2}/m_{1}$ ($q<1$ by definition) and $\tilde{\beta}$ is the power-law index. If $\tilde{\beta}=0$, the distribution is flat and pairing of sBHs with different masses does not have a mass preference. If $\tilde{\beta}>0$, the pairing in the binary system is more likely toward the equal-mass end $q=1$, and if $\tilde{\beta}<0$, the masses are more likely to be asymmetric with a smaller $q$.

The joint mass distribution of the BBH system can be obtained with the single BH mass function and the pairing function:

[TABLE]

where $N_{\rm b}$ is the normalization constant. The distribution of the primary mass component $m_{1}$ is

[TABLE]

The factor of $2$ comes from the symmetry of $m_{1}$ and $m_{2}$. The underlying parameters $\{\tilde{\alpha}_{\rm p},\ \tilde{\mu}_{m},\ \tilde{\sigma}_{m},\ \tilde{\lambda}_{\rm peak},\ \tilde{\delta}_{m}\}$ are unknown but can be determined (subject to data and fitting uncertainties) by matching $P_{\rm marginal}(m_{1})$ with the LVK data. As a baseline for the following discussion, we shall also adopt the power law index $\tilde{\beta}=2.6$ as measured by fitting the BBH population data [84]. Other values within the uncertainty range of [84] are also possible, although it does not affect the qualitative signatures of the conclusions.

For example, for a $P_{\rm LVK}(m_{1})$ as shown in Fig. 23, and with $\tilde{\beta}=2.6$, we find that the optimal fitting parameters are $\{\tilde{\alpha}_{\rm p}=3.2,\tilde{\mu}_{m}=34,\tilde{\sigma}_{m}=5.65,\tilde{\lambda}_{\rm peak}=0.058,\tilde{\delta}_{m}=2.3\}$. The solid blue line shows the LVK primary mass distribution, while the green dashed line represents the fitted primary mass distribution obtained by combining the single BH mass distribution with a pairing function, demonstrating promising agreement with the observed data. The orange dotted line corresponds to the inferred isolated BH mass function under the “Power-law + Gaussian peak” model with optimal fitting parameters. It is important to note that both the primary mass distribution $P_{\rm LVK}(m_{1})$ and the power-law index $\tilde{\beta}$ have substantial uncertainties [1, 2, 84] due to limited data. As a result, the inferred isolated BH mass function also carries inherent uncertainties (the shaded region in Fig. 23) and should not be interpreted as the “true” distribution. However, it serves as a useful reference point for exploring potential features of the EMRI mass distributions discussed in Sec.IV.1 and Sec.IV.2. Similar to $P_{\rm LVK}(m_{1})$, the inferred single-BH mass function exhibits a prominent peak around $10\,M_{\odot}$, as well as a secondary peak near $35\,M_{\odot}$. Following GWTC-2 [1], we adopt the same upper mass cutoff of $m_{\rm max}=86.22\,M_{\odot}$ for the isolated single black hole mass distribution.

We also make one key assumption to facilitate our discussion of EMRI mass functions: we assume that the isolated sBH mass function responsible for producing LVK binaries is the same as that of isolated sBHs within nuclear star clusters. This assumption is not necessarily valid, as the typical environments in which stellar-origin BBHs form may differ significantly from those of nuclear star clusters. However, in the absence of additional observational constraints on the sBH mass function within nuclear clusters, we proceed with the isolated sBH mass function inferred from LVK observations to explore the EMRI mass spectrum in the following two subsections.

IV.1 Dry EMRI

The dry EMRI mass spectrum can be obtained by combining the isolated BH mass function with the mass-dependent dry EMRI formation rate. In Sec. II.1, the modeling of a nuclear star cluster using the Fokker-Planck equation assumes a stellar mass of $1\ M_{\odot}$ and an sBH mass of $10\ M_{\odot}$ (a “1+10” model). If the masses of sBHs follow a continuous distribution, a reasonable approximation since stars make up the majority of the population, then different mass components of sBHs evolve independently under an averaged potential. Therefore, we can evaluate the dry EMRI formation rate for a $xM_{\odot}$ sBH in a “1+$x$” model, and then weight this rate by the isolated sBH mass function to obtain the resulting dry EMRI mass spectrum.

The EMRI rates for various mass components are shown in Fig. 24, based on the $1+x$ model described in Sec. II.1. Heavier sBHs tend to have higher EMRI rates due to the effect of mass segregation. After the initial ramp-up stage, the rates gradually decline over time as a result of the insufficient supply of sBHs to the nuclear star cluster. In the inset of Fig. 24, we show the dry EMRI rate per object averaged over time within a canonical time $5$Gyrs. This dry EMRI rate acts as a weighting factor in the final dry EMRI mass distribution. By performing a time average, we obtain the final dry EMRI mass spectrum as follows:

[TABLE]

where $\Gamma(m,t)$ is per-object EMRI rate and $N_{\rm dry}$ is the normalization constant.

The corresponding dry EMRI mass function is shown in Fig. 25. Compared to the initial isolated mass spectrum, the dry EMRI mass function exhibits a more pronounced peak around $\sim 35M_{\odot}$, which is higher than the peak in the LVK primary mass function and even higher than the peak at $10M_{\odot}$. An intriguing implication is that secondary BHs of dry EMRIs are of heavy tailed distribution, therefore lead to more EMRI detections. With future observations from space-based detectors, the dry EMRI mass function, combined with refined theoretical modeling of the relative rates among different mass components, can be used to infer the mass function of isolated sBHs within nuclear star clusters. This, in turn, can be compared to the mass function derived from LVK observations, serving both as a test of the assumption made in Sec. IV and as a probe of environmental effects on binary formation.

IV.2 Wet EMRI

Isolated sBHs may be captured by the AGN disk, if present, and subsequently migrate toward the central MBH to form wet EMRIs. Therefore, similar to dry EMRIs, the mass function of wet EMRIs also critically depends on the mass function of isolated sBHs within the nuclear star cluster. However, additional factors, namely disk capture efficiency, accretion effects, and merger effects, also influence the final mass distribution. These additional factors will be addressed quantitatively in this section.

Starting from the mass function $P_{\rm single}(m)$ of isolated sBHs in the nuclear cluster, the number of sBHs captured by the accretion disk also depends on the mass of the sBH and the time-dependent per-object capture rate $\Gamma_{\rm cap}(m,t)$, similar to the scenario in dry EMRI. During the course of an entire AGN episode, the average mass distribution after disk capture depends on the disk lifetime, which is

[TABLE]

The per-object captured rate $\Gamma_{\rm cap}(m,t)$ may be determined by applying an FP simulation of the system with both the nuclear star cluster and the incorporated disk, i.e. see the detailed modeling in [36]. Given similar setups, we compute the capture rates shown in Fig. 26. Its inset displays the time-averaged capture rate as a function of mass for different disk lifetimes. The resulting mass functions $P_{\rm cap}$ for different disk lifetimes are presented in the upper left panel of Fig. 27, where a relative shortage of more massive sBHs can be found for scenarios with larger $T_{\rm disk}$. On the other hand, compared to the initial single mass function $P_{\rm single}(m)$, there is already a relative boost in the fraction of more massive sBHs, because inclination-angle diffusion is more efficient for more massive sBHs.

sBHs captured by the AGN disk generically migrate within the disk due to the back reaction from density-wave emission, absorption, and momentum transfer with the head wind [23], and gravitational pull from possible turbulent eddies of the disk. In particular, the head wind arises from gravitational capture of the surrounding disk material by the sBH, which tends to circularize within the Bondi radius, and form a supercritical accretion flow towards the sBH. In [85] we have adopted the “inflow-outflow” model [86] to describe the supercritical accretion flow and studied the change in mass and spin of sBH during the migration process. Different simulations have further confirmed the validity of the “inflow-outflow” model in the super-Eddington regime [87, 88].

Assuming the same sBH accretion model used in [85], and incorporating a refined treatment on type II migration (the “gap-opening” scenario) as discussed in detail in Appendix. C [89, 90], we trace the accretion rate of sBHs during their migration processes and obtain their final mass $m_{\rm f}$ as a function of the initial mass $m_{\rm i}$: $m_{\rm f}=m_{\rm f}(m_{\rm i})$, or vice versa $m_{\rm i}=m_{\rm i}(m_{\rm f})$. It is straightforward to see that the final mass distribution is given by

[TABLE]

As found in [85], accretion generally increases the mass of an sBH by $<25\%$ for $\alpha$-disk and $<70\%$ for $\beta$-disk. This is also consistent with the mass spectra of wet EMRIs shown in Figs. 27 and 28. For reference, we also show the accreted mass in the lower panel of Fig. 30.

The third critical factor that affects the mass function of wet EMRIs is the pairing and merger of sBHs within the disk. It has been suggested that sBH binaries may form and merge in AGNs, potentially contributing to some of the observed events in the GWTC catalog [91, 92, 93, 94, 95]. The pairing process arises mainly from the collision of circum-single disks from individual sBHs, which subsequently move sufficiently close to each other [49, 48, 96, 50]. However, the formation rate, the subsequent evolution, and the merger rate are still highly uncertain [91]. In this study, we therefore adopt a parametrized approach to (crudely) assume that each sBH has a constant probability $f_{\rm mer}$ of being captured in a binary and eventually merge with its companion.

Assuming that all isolated sBHs within the disk would form pairs with another sBH that eventually leads to merger, that is, $f_{\rm mer}=1$, and that we only consider 1G(generation)$+$1G scenarios, the resulting wet EMRI mass function is

[TABLE]

where we have assumed a flat pairing function for simplicity and have not incorporated accretion in this part of the merger description. If the merger fraction is not $100\%$, that is, $f_{\rm mer}\neq 1$, the corresponding mass function becomes a linear combination of $P_{\rm f=1}$ and $P_{\rm cap}$:

[TABLE]

As we can find in Fig. 28, $P_{\rm f=1}$ clearly shows displaced peaks at $m=20M_{\odot}$ and $m=45M_{\odot}$, corresponding to the merger scenarios with $10M_{\odot}+10M_{\odot}$ and $10M_{\odot}+35M_{\odot}$. The mass peak corresponding to the $35M_{\odot}+35M_{\odot}$ merger process is weak due to the large relative abundance of sBHs in the $10M_{\odot}$ peak. In the upper right panel of Fig. 27, we show results for different fractions of merger $f_{\rm mer}$. Depending on the fraction of merger $f_{\rm mer}$, the $m=20M_{\odot}$ peak may manifest itself as a bump around $20M_{\odot}$, and another peak may change from $35M_{\odot}$ to $45M_{\odot}$.

Notice that the formation of sBH binary may happen at any stage of the migration process, so that the individual sBH in the binary may have acquired a certain amount of mass before reaching the center MBH. In addition, accretion onto these sBHs likely continues during the inspiral evolution of the binary. As a result, it is likely more accurate to apply a mixed distribution between $P_{\rm mer}$ and $P_{\rm acc}$ when analyzing the final wet EMRI mass spectrum. In practice, computing $P_{\rm mer}$ starting from $P_{\rm cap}$ as the initial distribution or replacing $P_{\rm cap}$ in Eq. 60 by $P_{\rm acc}$ leads to final distributions with only minor differences, as shown in the bottom left panel of Fig. 27. Therefore, we shall not distinguish the sequence between accretion and merger in our analysis.

In Fig. 27 we present the mass distribution of wet EMRIs assuming a canonical model with $M_{\bullet}=4\times 10^{6}M_{\odot}$, $\dot{M}_{\bullet}=0.1\dot{M}_{\rm Edd}$, $f_{\rm mer}=0.1$, $T_{\rm disk}=10^{6}{\rm yr}$ and $\alpha$-disk with $\alpha=0.1$. In order to study various parameter dependences, we also vary $T_{\rm disk},f_{\rm mer}$ and disk models to compare the resulting mass functions. For example, in the upper left panel, we observe a relative decline in the peak around $35M_{\odot}$ for longer disk lifetimes, mainly due to the depletion and insufficient supply of more massive sBHs in the nuclear star cluster. In the upper right panel, the mass functions with respect to different $f_{\rm mer}$ are shown. As $f_{\rm mer}$ increases, the peak around $35M_{\odot}$ gradually moves to the right ($\sim 45M_{\odot}$), and the peak around $20M_{\odot}$ becomes more pronounced. In the bottom left panel, we compute the mass function assuming different sequence of merger and accretion effects. The resulting distributions are largely similar. Finally, in the bottom right panel, we also compare the predictions from two different disk models: the $\alpha$-disk model and the $\beta$-disk mode. In general, the $\beta$-disk model yields a more massive distribution, due to higher AGN disk density, therefore its higher accretion rate onto isolated sBHs in the disk. In general, the mass spectrum of wet EMRI is most sensitive to $f_{\rm mer}$, the merger fraction.

In addition to sBHs that were first formed in the stellar cluster and then captured onto the AGN disk, in situ formation of sBHs within the disk is also possible. In the outer regions of an AGN disk, where the disk fragments into dense clumps as a result of gravitational instability, stars can form. The mass function of stars formed in AGN disks may be top-heavy compared to the conventional Salpeter initial mass function [16]. Stars in AGN disks may not undergo chemical evolution because the fresh gas from the disk replenishes stellar cores faster than hydrogen is burned [97, 98]. If this is true, stars on AGN disks are effectively “immortal” [99] and therefore are not expected to collapse into compact objects and form EMRIs. Otherwise, massive stars in AGN disks may eventually collapse into compact stellar-mass objects, which then migrate toward the central MBH and become wet EMRIs [13, 16]. In this scenario, the mass functions of stars and sBHs embedded in AGN disks are highly uncertain because of the complex stellar environments and interactions between stars [100, 97, 98, 99].

V Conclusion

In this work, we have performed a first population study on three important observables of EMRIs: eccentricity, inclination, and sBH mass, all of which exhibit distinct distributions for wet and dry EMRI formation channels. Taking into account their evolutionary histories starting from highly eccentric inspirals, dry EMRIs are likely to retain eccentricities of $e\gtrsim 10^{-2}$ at $r_{\rm p}=10M_{\bullet}$. In contrast, the eccentricities of wet EMRIs are strongly damped by the accretion disk, but can be reexcited through multibody resonances and/or turbulent fluctuations within the disk. In both cases, the eccentricity is typically below $10^{-2}$ at $r_{\rm p}=10M_{\bullet}$, making it an effective parameter to distinguish between wet and dry formation channels. Furthermore, measurements of wet EMRI eccentricities can provide insight into key disk parameters, as illustrated in Eqs. (35, 36).

The distributions of the EMRI inclination angles also exhibit interesting features. For dry EMRIs, since the angular momentum of the LSO depends on the inclination, the EMRI formation rate is intrinsically $\iota$-dependent for spinning MBHs. The “tilted” distribution of $\iota$ thus serves as a direct test of the dry EMRI formation model, although it may also be influenced by potential deviations from spherical symmetry in the nuclear star cluster. For wet EMRIs, the $\iota$ distribution is typically aligned or tilted, depending on whether MBH accretion is predominantly coherent or chaotic. In the chaotic accretion scenario, the detailed shape of the resulting tilted distribution can be used to constrain the relative timescales of disk lifetimes and BP alignment.

The mass function of sBHs in wet EMRIs is likely different from that of sBHs in nuclear star clusters and dry EMRIs, primarily due to three key factors: capture efficiency, accretion, and mergers. The disk capture process is more efficient for more massive sBHs, enhancing their relative abundance at the high-mass end. Accretion generally increases the mass of the sBHs, changing the distribution toward larger values. In contrast, sBH binary formation and mergers introduce distinctive features into the mass spectrum. Depending on the fraction of merger $f_{\rm mer}$, new peaks may emerge in the resulting mass function.

Additional EMRI observables are not discussed in this work. For example, the spin of sBHs can be altered by accretion [85, 101], leading to a potentially different spin distribution for wet EMRIs compared to dry EMRIs. However, measuring the spin of sBHs is challenging, as recent studies incorporating sBH spin into EMRI waveforms suggest [28, 102], unless the eccentricity is sufficiently large. Therefore, it is unlikely that the spins of wet EMRIs can be accurately measured. In contrast, the spin of MBHs can be determined with extremely high precision, at the level of $\sim 10^{-5}$ [9]. Since the formation rate of wet EMRIs is largely insensitive to MBH spin, the observed spin distribution of MBHs in wet EMRIs should reflect their intrinsic spin distribution. On the other hand, the formation rate of dry EMRIs depends on the angular momentum of LSO, and therefore on the MBH spin (see the discussion in Sec. III.2). As a result, the MBH spin distributions inferred from wet and dry EMRIs may differ significantly. Similar argument applies to the mass function of MBH in EMRIs as well, as the formation efficiency of wet EMRIs and dry EMRIs may have a different dependence on $M_{\bullet}$, in addition to being influenced by the intrinsic mass function of MBHs.

In addition to the distributions of individual parameters, correlations between parameters may also exist for different EMRI formation channels. For example, in dry EMRIs, the “tilt” in the distribution of $\iota$ is expected to be more pronounced for MBHs with higher spin, implying a positive correlation between $\cos\iota$ and the MBH spin $a$ in their joint distribution. Moreover, if AGN accretion is intrinsically correlated with the mass of the central MBH, that is, if $\dot{M}_{\bullet}$ correlates with $M_{\bullet}$, then the eccentricity of wet EMRIs may also exhibit a non-trivial dependence on $M_{\bullet}$. Such multiparameter correlations and joint distributions remain largely unexplored but may encode valuable information about the underlying EMRI formation channels.

Last but not least, although we have focused only on wet and dry EMRIs, EMRIs formed via the Hill mechanism or the Kozai-Lidov mechanism in SMBH binaries may also constitute a sub-population in the event catalog of space-based detectors. It would be interesting to develop quantitative population models for these alternative formation channels and to investigate the distributions of their associated observables.

Acknowledgements.

VI Acknowledgments

The authors thank Yanbei Chen and Xinyu Li for interesting discussions. Y.P.L is supported in part by the Natural Science Foundation of China (Grants 12373070, and 12192223), the Natural Science Foundation of Shanghai (grant NO. 23ZR1473700). H.Y. is supported by the Natural Science Foundation of China (Grant 12573048). The calculations have made use of the High Performance Computing Resource in the Core Facility for Advanced Research Computing at Shanghai Astronomical Observatory.

Appendix A Post-Newtonian Coefficient

In the multibody resonance section, we consider the post-Newtonian force between the sBH and the center MBH. The force depends on their position and velocity in the form of [47],

[TABLE]

the relative position $\vec{r}=\vec{r}_{\text{sBH}}-\vec{r}_{\text{MBH}}$, $r=|\vec{r}|$, the relative velocity $\vec{v}=\vec{v}_{\text{sBH}}-\vec{v}_{\text{MBH}}$, the radial velocity $v_{r}=\vec{v}\cdot\frac{\vec{r}}{r}$, the total mass $M=m+M_{\bullet}$, the symmetric mass ratio $\eta=mM_{\bullet}/M^{2}$, and the PN coefficients are

[TABLE]

Appendix B Equilibrium Eccentricity in Post-Newtonian Analysis

In this paper, we study the resonance between two sBHs orbiting a central MBH, taking into account post-Newtonian corrections. In planetary systems, the case of two planets orbiting a central protostar is very common, and extensive analytical work has been done on such systems, particularly on the equations of motion for the orbital elements under resonance. Here, we include the post-Newtonian terms in equations of motion, which yields good agreement between the equilibrium eccentricity and the actual eccentricity, as shown in Fig. 4. The equations of motion for the orbital elements are formulated as [103]

[TABLE]

We use the labels ‘in’ and ‘out’ to represent the inner and outer objects, respectively. These equations of motion satisfy only the $j+1:j$ resonance, where ${t_{J},\ t_{e}}$ are the angular momentum loss timescale and eccentricity damping timescale, due to type I migration forces and gravitational radiation. Here, $a$ is the semi-major axis, $e$ is the eccentricity, $\mathcal{R}_{a}=a_{\rm in}/a_{\rm out}$ is the ratio of the semi-major axis, $n$ is the mean motion, $\varpi$ is the pericenter longitude, $\varphi$ is the resonance angle, $f$ is the resonance coefficient (which can be found in Table 8.1 of [104]), and $\mu$ is the mass ratio $m/M_{\bullet}$. We have added only the simple precession term from the 1PN effect to $\dot{\varpi}$.

If the orbital elements reach equilibrium, we can get $\dot{n}_{\rm in}/n_{\rm in}=\dot{n}_{\rm out}/n_{\rm out}$, $\dot{e}_{\rm in}=\dot{e}_{\rm out}=0$, $\dot{\varpi}_{\rm in}=\dot{\varpi}_{\rm out}$, and we can obtain the equilibrium eccentricity equation

[TABLE]

we have ignored the eccentricity term $(1-e^{2})$ in the 1PN precession term for the sake of simplicity in the calculation. Using these two equations, we can obtain the equilibrium eccentricity.

To study the influence of the post-Newtonian effect on the evolution of a resonant pair, we plot the $9:8$ resonance in Fig. 29, considering only the disk migration force and the 2.5PN gravitational radiation force. The figure clearly shows that the equilibrium eccentricity increases without the 1PN effect, while it decreases when the 1PN effect is included.

Appendix C Migration and Accretion

In the wet EMRI mass spectrum section, we mention that the sBHs captured by the AGN disk will migrate toward the central MBH and accrete gas from the disk. For the accretion model, we adopt the same supercritical accretion method as in [85]. A secondary disk forms around a sBH in the AGN disk. The inflow from the outer boundary of the secondary disk will not be completely accreted by the sBH; most of the inflow will escape as an outflow. In [85], the relation between the inflow from the outer boundary, $\dot{m}_{\rm in}$, and the actual accretion rate $\dot{m}_{\rm in,0}$ was numerically calculated from the “inflow-outflow” model and fitted as

[TABLE]

in which $\dot{m}_{\rm Edd}$ is the Eddington accretion rate, $\chi$ is dimensionless spin parameter of the sBH, $r_{\rm obd}$ is the inflow out boundary (see Eqs. (26-32) in [85] for more details about the out boundary).

During migration, an sBH will undergo type II or type-I migration, depending on whether it opens a gap. If the sBH is more massive, it may open a gap, which will significantly reduce the surface density around it.

We do not consider type II migration in Section II.2.1 because the post-Newtonian effect operates in the inner region of the AGN disk ($<100M_{\bullet}$), which corresponds to type I migration (without gaps) according to the following criteria. A gap opens if the mass ratio is greater than or equal to [73, 23, 74, 75]

[TABLE]

in which $q_{\rm crit}$ is the minimum mass ratio that can open a gap, $\alpha$ is the dimensionless disk parameter, $\beta=p_{\rm gas}/p_{\rm tot}$, $b=0$ for $\alpha$-disk and $b=1$ for $\beta$-disk, $h$ is the disk aspect ratio. In our model here, we use $\alpha=0.1$, the center MBH mass $M_{\bullet}=4\times 10^{6}M_{\odot}$

In the upper panel of Fig. 30, we show the minimum sBH mass that can open a gap at different radii for the $\alpha$- and $\beta$-disk. For simplicity, during the sBH migration, we consider that it first undergoes type II migration until it reaches the radius where it can no longer open a gap and then undergoes type I migration.

For type I migration, the specific torques that act on the sBH have three components: [23, 85]

(1) Type I migration torque

[TABLE]

in which $C_{\rm I}$ is the torque coefficient on the order of unity, $m$ is the sBH mass, $M_{\bullet}$ is the center MBH mass, $\Sigma$ is the surface density of the AGN disk, $\Omega$ is the gas Kepler angular velocity, $h$ is the aspect ratio. All disk quantities here are measured at the orbital location of the sBH.

(2) Gravitational wave torque

[TABLE]

which is Peter’s formula [41], assuming a circular orbit.

(3) Wind torque

[TABLE]

in which $\delta v_{\phi}$ is the relative bulk velocity in the $\phi$ direction between the local gas and sBH, $\dot{m}_{\rm wind}$ is the head wind strength which determines $\dot{m}_{\rm in}$ at the outside boundary.

The head wind can be obtained from Bondi accretion

[TABLE]

in which $\dot{m}_{\rm BHL}$ is the Bondi accretion rate, $r_{\rm BHL}$ is the Bondi radius, $H$ is the AGN disk height, $r_{\rm Hill}$ is the sBH Hill radius [85].

For type II migration, we only need to consider two components of specific torques acting on the sBH [23, 85]:

(1) Type II migration torque [105, 106]

[TABLE]

with $K^{\prime}=0.04q^{2}h^{-5}\alpha^{-1}$ and $\zeta=3\pi\alpha(1+K^{\prime})\times\left(\sqrt{2/\pi}(qh^{-3})^{-2}+3^{-1/3}(qh^{-3})^{-2/3}\right)$, in which $q=m/M_{\bullet}$ is the mass ratio.

(2) Gravitational wave torque, same formula with Eq. 70

The accretion rate of the sBH can be prescribed by combining Bondi accretion and Hill accretion [89, 90]

[TABLE]

with

[TABLE]

and $\dot{m}_{\rm std}=3\pi\alpha\Sigma_{0}H^{2}\Omega$ is the steady accretion rate of the AGN disk. The outer boundary of the sBH accretion for the sBH that opens the gap (type II migration) is the secondary disk size $\sim r_{\rm Hill}/3$ [107, 89].

In the lower panel of Fig. 30, we show the ratio of the increase in mass after accretion compared to the initial mass from the outside $r\sim 1000M_{\bullet}$ of the AGN disk to the inside $\sim 10M_{\bullet}$. The sBH in $\beta$-disk always accretes more mass due to the dense gas environment than in $\alpha$-disk.

Appendix D Comparison with [42]

In a recent work [42], Mancieri et al. studied the eccentricity distribution of dry EMRIs at plunge using Monte Carlo simulations and find a distribution peak at $e_{\rm pl}\approx 0.2$ (see Fig. 4 in [42]). We note that their results are broadly consistent with ours. In [42], the numerical evolution of the nuclear stellar cluster ends when the time reaches ten times the EMRI rate peak time – where the “peak time” refers to when the EMRI rate reach a maximal value – while in this work we fix the evolution time as 5 Gyr. In addition, a slightly different definition of GW timescale $t_{\rm GW}$ are used in [108] in defining boundaries of computational domain: “$(1-e)$ decay timescale” due to gravitational radiation (see Eq. 20 in [108]). Another difference lies in the way the semi-major axis $a$ and eccentricity $e$ are calculated from orbital energy $E$ and angular momentum $J$. In this work, we compute the pericenter $r_{p}$ and apocenter $r_{a}$ in the Newtonian potential of the central MBH plus nuclear cluster, then define semi-major axis as $a=(r_{a}+r_{p})/2$ and eccentricity as $e=(r_{a}-r_{p})/(r_{a}+r_{p})$. In Ref. [42], the authors use orbital parameter definitions in Schwarzschild spacetime framework $(E,L)\to(p,e)$ [11], where $p$ is the semi-latus rectum.

In Fig. 31, we show the boundaries obtained when adopting different gravitational-radiation timescales. The difference is only visible in the low-eccentricity regime, where a minor fraction of EMRIs are produced. For comparison with [42], we now use the same relativistic definitions of $a$, $e$ and calculate the eccentricity distribution of dry EMRIs at plunge, i.e., at $p=6+2e$ [11], and show the results in Fig. 32. For $M_{\bullet}=1\times 10^{5}M_{\odot}$, we display two time snapshots. The snapshot in $t_{\rm ref}\approx 2\times 10^{8}{\rm yr}$ corresponds to ten times the EMRI-rate peak time, and the dry EMRI eccentricity profile is similar to that in [42]. The snapshot at $t_{\rm ref}=5\times 10^{9}{\rm yr}$ shows the distribution after 5 Gyr of evolution, indicating that the eccentricity distribution also evolves over time, with the peak decreasing. In addition, we also show the figures for two other MBHs. For $M_{\bullet}=7\times 10^{5}M_{\odot}$, the peak EMRI rate time is approximately 5 Gyr, and the figure shows a distribution similar to that for $M_{\bullet}=1\times 10^{5}M_{\odot}$ at 0.2 Gyr. For $M_{\bullet}=4\times 10^{6}M_{\odot}$, 5 Gyr is less than ten times the peak time of the EMRI rate, and the averaged eccentricity distribution shows a higher peak at $e\simeq 0.3$.

We thereby confirm that our results are broadly consistent with those in [42]. We also note that the average distribution evolves with time because of the limit supply. However, the distributions are not fully consistent because the authors include gravitational wave radiation in their Monte Carlo simulations, which can naturally produce the so-called “cliff EMRI” [109, 108]. As a result, an additional small fraction of EMRIs of high eccentricities are found [42]. We plan to update our FP code accordingly in the future.

Appendix E Physical Quantities

In this appendix, we summarize the physical quantities used in this paper and provide their physical meanings in Table LABEL:tab:symbols. This compilation is intended to serve as a convenient reference for readers and to ensure consistency in notation throughout the main text.