Formation Rate of Extreme Mass Ratio Inspirals in Active Galactic Nuclei
FFormation Rate of Extreme Mass Ratio Inspirals in Active Galactic Nucleus
Zhen Pan ∗ and Huan Yang
1, 2, † Perimeter Institute for Theoretical Physics, Ontario, N2L 2Y5, Canada University of Guelph, Guelph, Ontario N1G 2W1, Canada
Extreme Mass Ratio Inspirals (EMRIs) are important sources for space-borne gravitational wave detectors,such as LISA (Laser Interferometer Space Antenna) and TianQin. Previous EMRI rate studies have focusedon the “loss cone” scenario, where stellar-mass black holes (sBHs) are scattered into highly eccentric orbitsnear the central massive black hole (MBH) via multi-body interaction. In this work, we calculate the rate ofEMRIs of an alternative formation channel: EMRI formation assisted by the accretion flow around accretingmassive black holes. In this scenario, sBHs and stars on inclined orbits are captured by the accretion disk, andthen subsequently migrate towards the MBH, under the influence of density wave generation and head wind.By solving the Fokker-Planck equation incorporating both sBH-sBH / sBH-star scatterings and sBH / star-diskinteractions, we find that an accretion disk usually boosts the EMRI formation rate per individual MBH by O (10 − ) compared with the canonical “loss cone” formation channel. Taking into account that the fractionof active galactic nucleus (AGNs) is ∼ O (10 − − − ), where the MBHs are expected to be rapidly accreting,we expect EMRI formation assisted by AGN disks to be an important channel for all EMRIs observed by space-borne gravitational wave detectors. These two channels also predict distinct distributions of EMRI eccentricitiesand orbit inclinations, which can be tested by future gravitational wave observations. I. INTRODUCTION
With the Laser Interferometer Space Antenna (LISA) andTianQin planned for launch in early 2030s [1, 2], the mHzband will be available for gravitational-wave (GW) observa-tion. One primary target source of space-borne GW detectorsis extreme mass ratio inspiral (EMRI), which usually consistsof a stellar-mass compact object, e.g. a black hole or a neutronstar, and a massive black hole (MBH). The stellar-mass objectmay stay in the LISA band for years and complete 10 ∼ circles around the MBH before their final mergers [3]. Be-cause of such large number of cycles, small modification ofthe metric of the EMRI system [4] and possible environmen-tal / astrophysical e ff ects [5–7] may accumulate over the dura-tion of the waveform and generate detectable phase shift. Witha population of events, the distributions of masses and spins ofthe host MBHs may be measured as a way to infer the growthhistory of galactic centre MBHs [8, 9].One important problem related to EMRIs is to evaluate theirevent rate. The “canonical” EMRI formation is expected to bea stellar-mass black hole (sBH) captured by a MBH via multi-body scatterings in the core of a galaxy [10], which has beenthe main assumption for previous rate calculations [11, 12].Given the MBH mass and the initial distributions of surround-ing stars and sBHs, the EMRI rate per MBH can be obtainedby solving the Fokker-Planck equation or by N-body simu-lations [13, 14]. In addition to the generic rate per MBH, theLISA detectable EMRI rate also depends on the mass functionof MBHs at di ff erent redshifts, the fraction of MBH living instar clusters, and the number fraction of sBHs in star clusters.Taking account of all these uncertainties with semi-analyticmodels, Babak et al. [12] forecasted that there will be tens tothousands EMRIs detected by LISA per year, and in a simi- ∗ [email protected] † [email protected] lar analysis, Fan et al. [15] forecasted a slightly lower EMRIdetection rate by TianQin.In this work, we consider another possible EMRI formationchannel, where a sBH in the core of a galaxy is captured bythe accretion disk around an accreting MBH. As extensivelystudied in the context of star-disk-satellite systems [16–19],a planet within inclined orbit with repesct to the disk excitesdensity waves that drive the planet’s inward migration, cir-cularize the planet’s orbit and drive the planet orbit towardthe disk plane. In a MBH-disk-sBH system, similar processesalso work and assist the EMRI formation. A sBH initiallyresides on a inclined orbit crossing the accretion disk gener-ally moves towards lower incliantion orbits and eventually getcaptured by the accretion disk. sBHs within the disk generallyinteract with the disk through wind e ff ects, density wave gen-ration and dynamic friction. As we can find in the discussionsin Sec. III, sBHs migrate within the disk under these e ff ects.The timescale of sBH capture and migration apparently de-pends on the disk profile. In particular, certain disk models[20, 21] predict local density maxima at distance ∼ O (10 ) −O (10 ) times the MBH size away. As a result, the migrationtorque due to density wave generation changes sign at theselocations where the migrating object are trapped. This trap-ping mechanism has been extensively discussed in recent liter-ature [22–32] as a possile way to generate hierarchical stellar-mass black hole mergers that are observable by ground-basedgravitational-wave detectors. However, we notice that thesestudies have overlooked an important component of the diskforce: the head wind, which orginates from the accretion ofdisk materials onto the moving sBH [33]. With the wind in-fluence included, we find that the total torque is alway posi-tive for realistic disk parameters, so that the disk trap is un-likely giving rise to hierarchical stellar-mass binary mergers,nor does it stop the sBHs migrating towards the MBH to be-come EMRIs.We incorprate all the relevant disk-star / sBH interactionsinto a Fokker-Planck code, and then compute the EMRI for-mation rate with varying MBH mass and disk profiles. The a r X i v : . [ a s t r o - ph . H E ] F e b initial distribution of the star cluster is specified according toTremaine’s stellar cluster model [34]. In the limit of zero diske ff ects, th code reproduces known results from previous “loss-cone” calculations. With the presence of an “realistic” accre-tion disk (Sec. IV), we find that the disk-assisted EMRI for-mation is generically O (10 − ) times faster than the “losscone” mechanism for the same MBH. Taking into account thefraction of AGNs observed O (10 − − − ) [35, 36], we con-clude that disk-assisted EMRI formation may be an importantor even dominant channel for LISA EMRI observation.Interestingly, disk-assisted EMRIs tend to have low eccen-tricities ( e (cid:39)
0) and low inclinations ( ι (cid:46) .
1) when they enterthe LISA band. These distributions are very di ff erent fromthe ones predicted by “loss cone” formation (0–0.2 and 0– π/
2, respectively). This means that it is possible to seperateout these two channels with a population of events. With therate inferred for each channel from observations, one can thenfurther constrain the distribution of stars around MBHs andAGN / disk physics. In particular, EMRIs within AGNs mayproduce both gravitational wave and electromagnetic signalsfor multi-messenger observations [33, 37].The paper is organized as follows. In Section II, we firstbriefly review the canonical EMRI formation channel via losscone and numerically calculate the EMRI rate for a fiducialMBH + star / sBH cluster system. In Section III, we introducea few commonly used AGN accretion disk models, exploredi ff erent interactions between sBHs / stars with accretion disks,and discuss the existence problem of migration traps in AGNdisks. In Section IV, we incorporate the sBH / star-disk inter-actions into the Fokker-Planck equation and numerically cal-culate the accretion disk-assisted EMRI rate. Summary anddiscussion are given in Section V. Some numerical details areplaced in Appendix A and B.Throughout this paper, we use geometrical units G = c = II. REVIEW OF EMRI FORMATION VIA LOSS CONE
In this section, we briefly review how stars fall into a MBHvia the loss cone mechanism, following Refs. [38–40], andthen compute the EMRI formation rate for a given MBH-stellar cluster system by numerically solving the Fokker-Planck equation. Many techical details discussed are also use-ful for the rate calculation with the presence of AGN disks.
A. Basics
Let us consider a stellar mass BH (sBH) orbiting around aMBH with mass M • , which locates in the center of a galaxybeing surrounded by a stellar cluster with velocity dispersion σ . Assume the sBH is on an eccentric orbit with eccentricity e and semi-major axis length a , its specific orbit energy andspecific orbital angular momentum are E : = φ ( r ) − v = M • a , J = M • (cid:114) − e E , (1) where φ ( r ) = M • / r is the (positive) gravitational potential.For later convenience, we also define R ≡ J / J c ( E ), where J c ( E ) is the specific orbital angular momentum of a sBH withspecific energy E on a circular orbit. For the point-mass grav-itational potential, we have R = EJ / M • = − e . ThesBH gradually spirals inward as GW emission takes away en-ergy and angular momentum on a timescale t gw . On the otherhand, stars and sBHs in the cluster continously scatter the sBHvia mutual gravitational interactions, changing its orbit angu-lar momentum on a relaxation timescale t J . As a result, theinspiral is susceptible to scatterings if t gw > t J and the sBHwill either be scattered to an wider orbit or plunges into theMBH. Therefore a stable EMRI form only if the GW dissi-pation dominates over scatterings, i.e., the orbit must be tightand highly eccentric ( e →
1) to enable e ffi cient GW emission, t gw < t J .In the phase space, there is an region bounded by theenergy-dependent angular momentum, J lc ( E ), where stellarmass BHs initially populating the region promptly fall intothe MBH within one orbital period [41]. This part of phasespace is usually referred as the “loss cone”. The overall infallrate is set by the rate of di ff usion / relaxation processes whichdrive sBHs to the loss cone by successive two-body scatter-ings, and the condition of infall is written as P ( E ) < t J [42],where P ( E ) is the orbital period. For comparison, the con-dition of stable inspirals t gw < t J is much stronger, becausethe orbital period P ( E ) is usually much shorter than the GWdissipation timescale t gw . For the problem we are discussing,the relevant orbits are nearly zero-energy (with a (cid:29) M • and E (cid:39) J lc , bh ( E (cid:39) = M • . (2)The above discussion equally applies to stars around a MBH,except the star loss cone is determined by tidal disruption and J lc , star ( E ) is slightly larger than J lc , bh ( E ) [43]. For numericalconvenience, we simply take J lc , star ( E (cid:39) = M • . B. EMRI rate via loss cone
1. Initial condition
For a given MBH, to accurately compute the EMRI ratevia the loss-cone mechanism, we need to know the distribu-tion functions of the surrounding stars and sBHs, f i ( t , (cid:126) x , (cid:126) v )( i = star , bh). As argued in Refs. [44, 45], these distributionfunctions are approximately functions of the action variables: f i ≈ f i ( t , E , R ). In alignment with previous studies [13, 14],we use Tremaine’s MBH + stellar cluster model [34, 46] as theinitial condition for the Fokker-Planck evolution. Assumingthere are two components in the stellar cluster: light starswith mass m star and heavy sBH with mass m bh , and the totalstar / sBH mass in the cluster are M star and M bh , respectively,the number densities of stars and sBHs in the Tremaine’s clus-ter model are given by n star ( r ) = M star m star − γ π r a r γ ( r + r a ) − γ , n bh ( r ) = δ × n star ( r ) , (3)with r a being the density transition radius, γ being the den-sity scaling power index, and δ being the number fraction ofsBHs. From Eq. (3), we see n i ( r ) ∼ r − γ ( i = star, bh) for r (cid:28) r a , n i ( r ) ∼ r − − γ/ for r = r a , and n i ( r ) ∼ r − for r (cid:29) r a .Di ff erent combinations of model parameters γ and r a producerich cluster profiles. For example, the Galactic nuclear stellarcluster is approximately described by the Tremaine’s modelwith γ = . r a = r h = M • /σ , where the star densityprofile is n star ( r ) ∼ r − . within the influence radius r h and be-comes as steep as n star ( r ) ∼ r − at a distance a few times largerthan r h [47, 48]From the density profiles (3), one can obtain the magnitudeof the gravitational potential φ ( r ) = M • r + M star + M bh r a − γ − (cid:32) rr + r a (cid:33) − γ , (4)where the second term is contributed by the stars and sBHs.In the case that the initial distribution functions f i ( t = , E , R )only depend on the energy E , they are related to the positionspace number density by [34] f i ( t = , E , R ) = √ π ) ddE (cid:90) E dn i d φ d φ √ E − φ , (5)where n i ( r ) has been written as an implicit function of φ ( r ).In the more general case, f i depends on both E and R . Inorder to invert the distribution function f i ( E , R ) to find thenumber density n i ( r ), we first list the properties of star orbitsin given potential field φ ( r ) [45]. From the energy definition E = φ − v /
2, we have2( φ − E ) = v = J r + v r , (6)where v r is the radial velocity. For a circular orbit of energy E , its orbit radius and angular momentum J c are determinedby J c ( E ) = − r c φ (cid:48) ( r c ) , φ ( r c ) − E ) = J c r c . (7)For a general non-circular orbit with parameters ( E , R ), itsturning points (apsis / periapsis) r ± are determined by2( φ ( r ± ) − E ) = J r ± , (8)and its orbit period P ( E , R ) is defined as P ( E , R ) = (cid:90) r + r − drv r . (9) Defining the number density in the ( E , R ) phase space as N i ( E , R ) dEdR : = (cid:82) r + r − d rd v f i ( E , R ), we have [44, 45] N i ( E , R ) = π P ( E , R ) J c ( E ) f i ( E , R ): = C ( E , R ) f i ( E , R ) . (10)With these listed properties, one can show that the position-space number density n i ( r ) is related to the distribution func-tion f i ( E , R ) by [45] n i ( r ) = π r (cid:90) φ ( r )0 dEJ c ( E ) (cid:90) R max dRv r f i ( E , R ) , (11)where R max ( r , E ) = r ( φ ( r ) − E ) / J c ( E ), and v r ( r , E , R ) = φ − E ) − J / r = ( R max − R ) J c ( E ) / r . In the case of isotropicdistribution f i = f i ( E ), the above equation simplifies as [49] n i ( r ) = π (cid:90) φ ( r )0 dE (cid:112) φ ( r ) − E ) f i ( E ) . (12)
2. Fokker-Planck equation
Given initial distributions of stars and sBHs, f i ( t = , E , R ),their evolution is governed by the orbit-averaged Fokker-Planck equation [45] C ∂ f ∂ t = − ∂∂ E F E − ∂∂ R F R , (13)with C the weight function defined in Eq. (10) and F E , R theflux in the E / R direction: − F E = D EE ∂ f ∂ E + D ER ∂ f ∂ R + D E f , − F R = D RR ∂ f ∂ R + D ER ∂ f ∂ E + D R f , (14)where the di ff usion coe ffi cients {D EE , D ER , D RR } i and theadvection coe ffi cients {D E , D R } i are functions of f i ( t , E , R )[45, 50] and we detail their calculation in Appendix A. In par-ticular, the local relaxation timescale of the system is approx-imately [44, 50–52] t rlx ( r ) = . σ (cid:80) i n i ( r ) m i ln Λ , (15)where the Coulomb’s logarithm ln Λ weakly depends the totalnumber of stars within the influence radius and we take ln Λ =
10 in this work.We aim to evolve f i ( t , E , R ) according to Eq. (13) and sub-ject to following boundary conditions [53]. On the E → f i ( t , E , R ) | E → = f i ( t = , E , R ) | E → , (16)i.e., the distributions far away from the central MBH barelyevolve due to its long relaxation timescale. On the R = R direction should vanish for bothstars and sBHs, F R | R → = . (17) log ( E ) R log ( f star ) log ( E ) R log ( f bh ) E f E f ( t ) f star ( t i ) f bh ( t i ) f star ( t f ) f bh ( t f ) FIG. 1. Final distribution functions f star ( t f , E , R ) (left panel), f bh ( t f , E , R ) (middle panel). Inital and final R -integrated distribution functions(right panel). All the distribution functions are shown in units of 10 pc − / (2 πσ ) / and energies E are shown in units of σ . On the loss cone boundary R = R lc ( E ) : = J / J c ( E ), the fluxin the R direction has been derived in Ref. [44] as − F R C (cid:39) (cid:18) D RR R (cid:19) R → f ( R )ln( R / R lc ) + F ( y lc ) , (18)where D RR : = D RR / C , y lc : = R lc / [( D RR / R ) R → P ], P = P ( E , R ) is the orbital period, R is any small R in the rangeof R lc ≤ R (cid:28) F ( y lc ) (cid:39) / y lc for y lc (cid:46) F ( y lc ) (cid:39) . y − / for y lc (cid:38)
1. For y lc (cid:29) R = R lc , Eq. (18) issimplified as f ( R lc ) (cid:39)
3. EMRI rate with the loss-cone mechanism
We consider a fiducial model of MBH + star / sBH clusterwith M • = × M (cid:12) , and two components in the cluster: lightstars m star = M (cid:12) , and heavy BHs with mass m bh = M (cid:12) . Weassume that stellar velocity dispersion σ follows the M • − σ relation [54, 55] M • = . × M (cid:12) (cid:32) σ / s (cid:33) . . (19)The initial star / sBH distributions are specified according tothe Tremaine’s model outlined in Section. II B 1 with the to-tal star mass M star = M • , the density transition radius r a = r h = M • /σ , the density power index γ = . δ = − . For reference, theinitial total number of stars within the influence sphere is N star ( r < r h ) = . × in this stellar cluster.To numerically solve the Eq. (13), we introduce a new di-mensionless variable Z = ln(1 + E /σ ) following Ref. [44]and implement an uniform 128 ×
128 grid on the ( R , Z ) space.From the initial distribution functions f i ( t = , E , R ), we com-pute all the di ff usion coe ffi cients and the advection coe ffi -cients, with which we evolve a discretized version of Eq. (13)with time step δ t = . C j , k min . ∆ R D R , ( ∆ R ) D RR , ∆ Z D E dEdZ , ( ∆ Z ) D EE (cid:32) dEdZ (cid:33) i , j , k , where i = { star , sBH } , j and k are the grid indices. We updatethe coe ffi cients according to the new distribution functions ev-ery 100 steps, and we stop the simulation at t f = × × f i ( t f , E , R ) in thefirst two panels and the R -integrated distribution functions¯ f i ( t , E ) : = (cid:90) f i ( t , E , R ) dR , (20)in the third panel, where ¯ f star ( t f , E ) ∼ E and ¯ f bh ( t f , E ) ∼ E . for E /σ (cid:38)
10. The sBH density profile ¯ f bh here turns outto be shallower than that obtained from solving 1-d Fokker-Planck equations [13]. The steeper profile from 1-d calcula-tion is expected because sBHs’ leaving the cluster and fallinginto the MBH via the loss cone cannot be incorporated in the1-d Fokker-Planck equation, and consequently sBHs whichare supposed to fall into the MBH, are instead accumulatedin the large E regime. The phase space distribution func-tion being ¯ f ( E ) ∼ E ξ is approximately equivalent to theposition space density being n ( r ) ∼ r − (3 / + ξ ) [42], As a re-sult, the star / sBH number density profile in the final state are n star ( r ) ∼ r − . and n bh ( r ) ∼ r − . for r (cid:46) . r h . We can alsofind that the density profile of the more massive, sBH com-ponent is steeper than that of the star component, a phenon-menon known as mass segregation, which has been shown toenhance the EMRI rate [13, 14, 56].In fact, mass segregation not only occurss in the E - dimen-sion, but also in the R -dimension. In Fig. 2, we show the finaldistributions of the two components at E = E gw [see Eq. (21)]in the R -direction, where the sBH component is more concen-trated on circular orbits. As a result, the mass segregation inthe R -direction is expected to mildly reduce the EMRI rate,which is mainly contributed by sBHs on highly eccentric or-bits.With the distribution function f bh ( t , E , R ) (and all the dif-fusion and advection coe ffi cients) obtained, we are ready tocalculate the EMRI rate for the loss cone mechanism. As dis-cussed in Refs. [14, 39], the EMRI condition t gw < t J is ap-proximately formulated as a < r gw = . r h . Therefore the R f i ( E g w , R ) f i ( E g w , R = ) f star f bh FIG. 2. Mass segregation in the R -direction, where the heaviercomponent is more concentrated on circular orbits. t [Gyr] X ( t )[ G y r ] F E ( t )| E = E gw emri ( t ) dN bh ( E > E gw )/ dt FIG. 3. The time dependence of the inflow rate ¯ F E | E = E gw , the EMRIrate Γ emri and the sBH number growth rate dN bh ( E > E gw ) / dt . EMRI rate per MBH via loss cone is given by Γ emri = (cid:90) E > E gw (cid:126) F · d (cid:126) l , (21)where E gw = M • / (2 r gw ), (cid:126) F = ( F E , F R ), and d (cid:126) l = ( dE , dR ) isthe line element along the boundary of the loss cone. Accord-ing to the flux conservation in the steady state, the EMRI rateshould be equal to the inflow rate ¯ F E at E = E gw ¯ F E | E = E gw = (cid:90) F E ( E , R ) | E = E gw dR . (22)In Fig. 3, we plot the time dependence of three di ff erent rates:¯ F E | E = E gw , Γ emri , and the sBH number growth rate dN bh ( E > E gw ) / dt . We find that the EMRI rate has reached an quasi-steady state at t ∼ Γ emri (cid:39)
210 Gyr − . In ad-dition, the inflow rate ¯ F E | E = E gw initially increases upto t ∼ F E | E = E gw (cid:39) Γ emri at t f = N bh | E > E gw increases from its initialvalue (cid:39)
10 to the final value (cid:39) N bh ( E > E gw ) starts to slowly decrease because theconsumption rate via loss cone Γ emri becomes larger than theinflow supply rate ¯ F E | E = E gw . e m r i [ G y r ] = , = 1.5= 2 , = 1.5= , = 1.8 M M / M e m r i [ G y r ] FIG. 4. Upper panel:the dependence of the characteristic EMRIrate per MBH ˆ Γ emri on the MBH mass M • for di ff erent initial clustermodels, where ( δ = δ = − , γ = .
5) is our fiducial model, ( δ = δ , γ = .
5) is a similar cluster model with higher sBH fraction,and ( δ = δ , γ = .
8) is a Galactic nuclear stellar cluster like model.Lower panel: the average rate ¯ Γ emri . To explore the dependence of EMRI rate on the MBH mass,we have performed additional 4 simulations similar to thefiducial model case, except with di ff erent MBH mass M • . Toquantify the EMRI rate via loss cone, we define the averagerate ¯ Γ emri : = t f (cid:90) t f Γ emri ( t ) dt , (23)and a characteristic rateˆ Γ emri : = Γ emri ( t ) | when N bh ( t , E > E gw ) maximizes , (24)i.e., the EMRI rate when the total sBH number N bh ( E > E gw )is maximal in the range of t ∈ [0 , t f ]. For the fiducial modelwe see that ˆ Γ emri = Γ emri ( t = t f ) (see Fig. 3), while ˆ Γ emri turnsout to be the EMRI rate at some earlier time Γ emri ( t < t f ) forcases with lighter MBHs, for which the relaxation timescalesare shorter and the peak N sBH ( E > E gw ) comes earlier. In theupper panel of Fig. 4, we show the characteristic rate ˆ Γ emri asa function of M • . We find that the MBH mass dependence ismild, which can be approximated as ∝ M − . • . This scalingcan be qualitatively understood as follows [39]:ˆ Γ emri ∝ N bh ( E > E gw ) t rlx ( r gw ) ∝ σ M • , (25)where we have used N bh ( E > E gw ) ∝ M • , t rlx ( r gw ) ∝ σ / n star ( r gw ) ∝ σ r / N star ( E > E gw ) ∝ σ r / M • , and r gw ∝ M • /σ . In combination with the M • − σ relation[Eq. (19)], we have ˆ Γ emri ∝ M − . • , which is close to the scal-ing fitted from our numerical results.To explore the e ff ect of the initial cluster profile, we also runtwo sets of simulations with two di ff erent initial stellar clustermodels: one with higher sBH fraction ( δ = × − , γ = . δ = − , γ = . N star ( r < r h ) = . × . Wefind the dependence of ˆ Γ emri on M • can be approximately fittedby the same scaling ˆ Γ emri ∝ M − . • in both cases. As expected,both higher sBH fraction and steeper density profile enhancethe characteristic rate ˆ Γ emri . In more details, the EMRI rate in-creases by less than a factor of 2 when we double the sBHfraction, because the sBH-sBH coupling becomes strongerwhich tends to flatten the density profile of sBHs. The EMRIrate ˆ Γ emri in the stellar cluster with a steeper density profileis higher by a factor of ∼
3, which is mainly contributed byhigher sBH numbers and shorter relaxation timescale (higherdensity of stars) within the influence sphere. The normalizedEMRI rate ˆ Γ emri / N ( r < r h ) in the two stellar clusters di ff ersby only ∼ Γ emri is more relevant,whose dependence on the MBH mass M • and on the initialcluster model are shown in the lower panel of Fig. 4. Wesee no simple scaling between the average rate and the MBHmass for any initial cluster model. The overall behavior isthe ratio ¯ Γ emri / ˆ Γ emri is lower for lower MBH mass, while theratio depends more on the initial condition for MBHs on thehigh mass end. For a MBH on the low mass end, sBHs in thevicinity of the MBH are rapidly depleted and the average rateis limited by the number of sBHs available. For a MBH on thehigh mass end, the relaxation timescale of the system is longand the two rates are di ff erent by no more than a factor of 2.In Ref. [14], stars and sBHs in the cluster are evolved fol-lowing 1-d ( E -direction) Fokker-Planck equations, in whichno star / sBH loss is takend account of. Therefore the num-bers of stars and sBHs in the cluster are conserved and thereexists a steady state, base on which the EMRI rate was calcu-lated assuming the standard logarithmic distribution in the R -direction and the relaxation timescale (15) in the R -direction[39]. Though there is no strictly steady state in our approachbecause the cluster continually loses stars / sBHs via the losscone. The characteristic rate ˆ Γ emri is roughtly comparable withthe steady-state EMRI rate in Ref. [14], in terms of either mag-nitude or the scaling with M • . III. SBH-ACCRETION DISK INTERACTIONS
About 1% low-redshift ( z (cid:46)
1) galaxies and as high as 10%high-redshift (1 (cid:46) z (cid:46)
3) galaxies are active [35, 36] in whichMBHs are expected to be rapidly accreting gas in a disk con-figuration. The interaction with disk could completely reshapethe distribution of sBHs. In this section, we will first discussmodels of AGN accretion disks, then introduce relevant disk-sBH interactions, including density waves and wind. Otherpossible interactions, e.g., dynamic friction [57, 58] and heat-ing torque [59, 60], are negligible as we will explain in Sec-tion III C.
A. AGN disk models α/β -disk:
With the α -viscosity prescription [61] and thethin disk assumption, the 1-d disk structure is governed by thefollowing equations [20]: σ SB T ff = π ˙ M • Ω , (26a) T = (cid:32) τ + + τ (cid:33) T ff , (26b) τ = κ Σ , (26c) β b H Ω Σ = ˙ M • Ω πα , (26d) p rad = τ σ SB T ff , (26e) p gas = ρ kTm H , (26f) β = p gas p rad + p gas (26g) Σ = ρ H , (26h) c s = H Ω = (cid:112) p tot /ρ , (26i) κ = κ ( ρ, T mid ) , (26j)with T mid being the middle plane temperature, T e ff beingthe e ff ective radiation temperature, τ being the disk opaticaldepth, H = rh being the scale height of the disk, c s being thelocal sound speed and κ being the gas opaticity [62, 63]. Theparameter b can be either 0 or 1 depending on whether the vis-cosity is proportional to the toal pressure ( α -disk) or the gaspressure ( β -disk).In the outer region, the disk will be prone to the disk self-gravity if the Toomre’s stability parameter Q = c s Ω π Σ (cid:39) Ω πρ (27)is less than unity. Following Ref. [20], we assume some ex-ternal feeding back mechanism heats the disk and maintainsa minimum value of the Toomre’s parameter as Q min (cid:39)
1. Inouter parts where Q = Q min , Eq. (26a) is replaced by Eq. (27).Given the value of accretion rate ˙ M • and the gasopacity function κ ( ρ, T mid ), all the disk variables { T e ff , T mid , τ, Σ , ρ, H , p rad , p gas , β, c s , κ } can be numericallysolved as functions of the angular velocity Ω , which isapproximated by the Kepler angular velocity Ω = (cid:112) − φ (cid:48) ( r ) / r .In this work, we assume the potential φ ( r ) [Eq.(4)] of thefiducial MBH + cluster model.As two fiducial disk models, we caculate the structure ofan α -disk and an β -disk both with M • = × M (cid:12) , α = .
1, ˙ M • = . M Edd • , where the Eddington accretion rate isrelated to the Eddington luminosity by ˙ M Edd • : = L Edd • / .
1. Weplot the fiducial α/β -disk structure in left panels of Fig. 5: thedisk surface density Σ , the middle plane temperature T mid , thedisk optical depth τ and the disk aspect ratio h as functionsof radius r . The two disks only di ff er within radius ∼ M • ,beyond which the gas pressure dominates over the radiationpressure, so that the di ff erence in the viscosity prescriptionsof the two disks is negligible. In the inner region, the radiationpressure dominates, so that the viscosity in the α -disk is largerthan in the β -disk, which results in a larger radial gas velocityand a lower gas surface density. TQM disk:
In the TQM disk model [21], the disk an-gular momentum is assumed to be carried away by globaltorques instead of local viscosity, and the gas inflow velocityis parameterized as a constant fraction of local sound speed: v gas , r = Xc s . In outer parts of the disk, star formation is as-sumed to heat the disk and maintain its stability against diskself gravity. In addition to the gas pressure and the radiationpressure, a turbulence pressure driven by supernova explosionin the disk is also incorporated. The 1-d disk structure is gov-erned by following equations: σ SB T ff = π ˙ M • Ω + (cid:15) (cid:63) ˙ Σ (cid:63) , (28a) T = (cid:32) τ + + τ (cid:33) T ff , (28b) τ = κ Σ , (28c)˙ M ( r ) = Xc s (2 π r Σ ) , (28d) c s = H Ω = (cid:112) p tot /ρ , (28e) p gas = ρ kT mid m H , (28f) p rad = τ σ SB T ff , (28g) p tb = (cid:15) (cid:63) ˙ Σ (cid:63) , (28h) Σ = ρ H , (28i)˙ M ( r ) = ˙ M • + (cid:90) rr min π r ˙ Σ (cid:63) dr , (28j) κ = κ ( ρ, T mid ) . (28k)In inner parts where Q >
1, the star formation ceases ( ˙ Σ (cid:63) = M ≡ ˙ M • and theturbulence pressure p tb vanishes. In outer parts, the Toomre’sstability parameter is assumed to be Q = ρ = Ω π . (29)In Fig. 6, we show the structure of an example TQM diskwith M • = × M (cid:12) , ˙ M • = . M Edd • , X = . (cid:15) (cid:63) = − .A salient feature in the disk is a opacity gap at r ∼ M • ,inside which the disk is optically thin τ <
1. In Ref. [21],a sharp density increase was found on the inner edge of theopacity gap, while the density increase in our solution is rathermild. This di ff erence is traced back to di ff erent equations ofradiation pressure assumed: in Ref. [21], p rad = σ SB T wasassumed, which should hold only in the optically thick regimeand break down inside the opacity gap, while our equation ofradiation pressure [Eq. (28g)] is more general.The α -viscosity prescription is consistent with the turbu-lence viscosity driven by magnetorotational instability in in-ner parts of accretion disks where the gas is fully ionized[64–66]. In outer parts, the physical mechanisms of the an-gular momentum transport and the external heating processes(in addition to the disk viscosity heating) maintaining thedisk stability are still open issues. In this work, we followRef. [20] to consider α/β disks assuming α -viscosity prescrip-tion throughout the disk and certain implicit heating processin outer parts of the disk. In consistent with Ref. [21], wealso consider TQM disks where a more e ffi cient angular mo-mentum transport mechanism is assumed. In addition, starformation in outer parts of AGN disk is explictly taken intoaccount as the external heating process. To our best knowl-edge, we expect α/β disk models to be a closer description toinner parts of AGN disks in nature, while it is not clear whichdisk model works better or whether any of them accuratelydescribes the nature in outer regions. B. Density waves
As extensively studied in the context of star-disk-satellitesystems [16–19], a planet excites density waves consisting ofthree components: regular density waves excited excited bythe circular motion of the planet, eccentricity waves excitedby the non-circular motion and bending waves excited by themotion normal to the disk. The regular density waves exerta negative torque on the planet and drive an inward migra-tion (commonly called type-I migration) on a timescale t mig , I ;the eccentricity waves work to damp the eccentricity e of theplanet orbits on a timescale t wav and the bending waves workto drive the planet onto the disk on a same timescale t wav . Sim-ilar processes should also work in the MBH-disk-sBH system,with torque arising from density waves [18, 19]˙ J mig , I = C I m bh M Σ M r Ω h , (30)where M = M ( < r ) is the total mass within radius r , C I = − . + d ln Σ / d ln r + . d ln T mid / d ln r [67]. In some spe-cial cases we will show later, the surface density Σ ( r ) is a fastincreasing function of radius r , and the type-I torque becomespositive. The corresponding migration timescale t mig , I and ec-centricity / inclination damping timescale t wav are t mig , I = J | ˙ J mig , I | = r Ω | ˙ J mig , I | ∼ Mm bh M Σ r h Ω , t wav = Mm bh M Σ r h Ω , (31)where in the “ ∼ ” sign we have take | C I | = ff usion term in theFokker-Planck equation as long as it is better quantified.A gap in the disk might be opened if the sBH is so massivethat its tidal torque moves gas away faster than the viscosityreplenishing rate. The gap width can be estimated as [33, 71–73] ∆ (cid:39) . r Ω ν m M / r , (32)where ν = ˙ M / (3 π Σ ) is the kinetic viscosity coe ffi cient. Thegap opening requires H < ∆ , r Hill < ∆ , (33)where H is the disk thickness and r Hill = ( m bh / M ) / r is theHill radius of the sBH inside which the tidal field of the sBHdominates. As long as a gap opens, type-I migration turnso ff and the sBH is subject to type-II migration. FollowingRef. [74], the type-II torque on the sBH can be estimated as˙ J mig , II = − π r Σ m bh r Ω | v gas , r | , (34)where v gas , r = − ˙ M / (2 π r Σ ) is the gas inflow velocity. Thecorresponding timescale of type-II migrationis defined as t mig , II = r Ω | ˙ J mig , II | . (35) The above analysis equally applies to stars except with a lowermass m star .In the three fiducial disk models (Figs. 5 and 6), the gapopening condition (33) is not satisfied. In the upper row ofFig. 7, we show the structure of an comparison α -disk withlow viscosity α = .
01 where the gap opening condition issatisfied if a sBH orbits around the MBH in the range of ∼ (10 , ) M • . C. Wind
For a sBH embedded in the gas disk, its gravitational at-traction inluences the surrounding gas meterials, so that theytend to flow towards to the sBH. If the disk is not rotatingand the sBH has no relative motion with repsect to the disk,these gases should flow towards the sBH in a nearly-sphericalmanner. On the other hand, if the disk is rotating and the sBHhas nonzero velocity relative to nearby materials, the accretioncannot be spherical. In addition, the accreting materials gen-erally carry nonzero angular momentum relative to the sBH,so that they tend to circularize and form certain local disk orbuldge profile to organize the accretion flow. Depending onthe heating processes and magnetic fields, a major part of cap-tured materials may be carried away in the form of outflowand only the remaining part is accreted [75, 76]. Becaues ofthe circularization process, it is reasonable to expect that theoutflow materials carry minimal net momentum with respectto the sBH.As a result, the “head wind” with respect to the sBH arecaptured at places where the sBH gravity becomes important,and the momentum carried by the wind eventually transfers tothe sBH. The amount of gas captured per unit time ˙ m gas can beestimated according to the Bondi-Hoyle-Lyttleton (BHL) ac-cetion rate (plus some disk environment corrections [see Ref.33, for a detailed analysis])˙ m BHLgas m bh = πρ m bh ( v + c s ) / , (36)within Bondi radius r BHL = m bh / ( v + c s ), where v rel is therelative velocity between the sBH and the local gas, v = ( δ v φ + δ v dr ) + δ v r , with δ v φ , δ v r being the relative bulk ve-locity in the φ , r direction, and δ v dr being the relative velocitycoming from the di ff erential rotation of the gas [33]: δ v φ = − γ hc s (37a) δ v r = | v gas , r − v bh , r | = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − ˙ M π r Σ − ˙ JdJ / dr (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (37b) δ v dr = (cid:18) m bh M (cid:19) / h − c s , (37c)where γ = d ln ρ/ d ln r and ˙ J is the change rate of the specificangular momentum of the sBH (or equivalently the specific r / M T mid h r / M J gw J wind + J mig,I J mig,I r / M t gw t wind t mig,I t wav r / M T mid h r / M J gw J wind + J mig,I J mig,I r / M t gw t wind t mig,I t wav FIG. 5. Fiducial α/β -disk in the upper / lower row. Left panel: disk structure with Σ [g / cm ] the surface density, T [eV] the middle planetemperature, τ the disk optical depth and h : = H / r the disk aspect ratio. Middle panel: the torques exerted on the sBH from GW emission˙ J gw [38], disk wind ˙ J wind [40] and density waves ˙ J mig , I [30], where ˙ J mig , I changes its sign close to the local density maxima. Right panel:the corresponding timescales on which di ff erent torque change the sBH orbital angular momentum by order of unity t i : = J / | ˙ J i | . with i = gw / wind / mig , I. For comparison, the typical disk life span ∼ (10 , ) yrs is plotted as a horizontal gray band. r / M T mid h r / M J gw J wind + J mig,I J mig,I r / M t gw t wind t mig,I t wav FIG. 6. Same to Fig. 5 except for a fiducial TQM disk. torque exerted on the sBH) due to sBH-disk interactions andGW emission, i.e., ˙ J = ˙ J mig , I , II + ˙ J wind + ˙ J gw , where˙ J gw = − m bh M (cid:18) Mr (cid:19) / (38)is the angular momentum loss rate due to GW emission (as-suming a circular orbit), J wind is the loss rate due to the windcapture [Eq. (40)], and ˙ J mig , I , II = ˙ J mig , I or ˙ J mig , II depends onwhich type of migration is operating.In the disk environment, the above accretion rate is subjectto corrections from limited disk thickness H and finite Hill radius r Hill of the sBH [33]˙ m gas = ˙ m BHLgas × min . { , H / r BHL , r Hill / r BHL } . (39)As a result, the specific torque exerted on the sBH from thewind is written as ˙ J idwind = − r δ v φ ˙ m gas m bh , (40)where the upper script “id” is used to denote quantities of in-disk (id) sBHs. For sBHs with inclined orbits (so that part0 r / M T mid h r / M J gw J wind + J mig,I,II J mig,I,II r / M t gw t wind t mig,I,II t wav r / M T mid h r / M J gw J wind + J mig,I J mig,I r / M t gw t wind t mig,I t wav FIG. 7. Upper row:same to Fig. 5 except for a comparison α -disk with low viscosity α = .
01, where the gap opening condition is satisfiedwhen an orbiting sBH is located in the range of ∼ (10 , ) M • (vertical gray band in the middle panel).Lower row:same to Fig. 5 except for a comparison α -disk with M • = × M (cid:12) and ˙ M • = . M Edd • , where the type-I torque becomes positivein two disconnected regions. of their orbits are outside of the disk (od)), when they hit theaccretion disk, the relative velocity v rel ∼ r Ω (cid:29) c s is usuallymuch greater than that of the in-disk sBHs. As a result, thewind capture radius is much smaller which greatly reducesthe wind e ff ect, so that we simply take ˙ J odwind =
0. In summary,migration timescales of in-disk sBHs and those outside are t bh , idmig = J | ˙ J mig , I , II + ˙ J gw + ˙ J wind | , t bh , odmig = J | ˙ J mig , I + ˙ J gw | , (41)where we take ˙ J wind = ff er-ent radii ( r sBH ± H ) is supersonic, and the dynamic friction withrespect to gas in this region may be important. In this case, thesupersonic relative velocity is mostly due to the shear (di ff er-ential rotation) of the accretion flow instead of the local pres-sure gradient, and the standard dynamic friction is automati-cally incorporated in the migration torque [58]. For inclinedsBHs, dynamic friction is still weaker than the e ff ect of den-sity waves : consider a sBH on an inclined orbit, penetratingthe gas disk with relative velocity v rel ∼ r Ω , the dynamic fric-tion (per unit mass) on the sBH is f df ∼ G m bh ρ/ v [57, 58],and the timescale for the dynamic friction to change the sBH’sorbit is t df ∼ ( r Ω ) / ( f df H Ω ), which is ∼ t mig , I h − (cid:29) t mig , I . Another possible contribution which we do not include hereis the heating torque [59] arising from the asymmetric dis-tribution of low-density gas around the sBH due to the ac-cretion heating and the shear of disk flow. As estimated in[60], the heating torque might be comparable with the type-Imigration torque assuming a thermal feedback of Eddingtonluminosity, but neither gravity nor dynamical feedback fromthe sBH. In fact, the local gas distribution is sensitive to boththe sBH gravity and the gas outflow which carries away theangular momentum of the sBH as explained in the beginningof the section. A simple estimate shows that the disk flowshear timescale 1 / ∆Ω ∼ / ( λ c d Ω / dr ) > / ( Ω h ) for the char-acteristic size λ c < H of the thermal feedback is much longerthan the dynamical timescale λ c / c s < / Ω on the same lengthscale. Therefore the local gas distribution should be more sen-sitive to the dynamical processes of gas inflow and outflow,and the heating torque should be much weaker than the esti-mate assuming neither sBH gravity nor dynamical feedback.Equating the in-disk migration timescale with the disk lifetime T disk defines a critical radius within which sBH can suc-cessfully migrate into the central MBH, i.e., the critical radius r crit ( T disk ) is defined by (cid:90) r crit r min t bh , idmig ( r ) d ln r = T disk . (42)In right panels of Fig. 5, we plot the timescales of di ff erentprocesses: t mig , I , t wav [Eq. (31)], t gw : = J / | ˙ J gw | [Eq. (38)], and1 t wind : = J / | ˙ J wind | [Eq. (40)].Di ff erent from sBHs, the gas accretion onto stars is morecomplicated considering that both the radiation heating andsolar wind are supposed to alter the local gas environment, byheating up and blowing away surrounding gas. In the presenceof an strong isotropic outflow from stars, the star-gas interac-tion could be completely di ff erent from the classical Bondiaccretion [77, 78]. Here we simply take ˙ J starwind = t star , idmig (cid:39) t star , odmig = m bh m star t bh , odmig . (43) D. Migration traps
In middel panels of Fig. 5, we plot the torques exerted onthe sBH ˙ J gw , ˙ J mig , I and ˙ J wind for α and β type of disks, respec-tively. For the α disk, we find that the torque J mig , I changesits sign around the local density maxima r ∼ M • , which isknown as the migration trap [79, 80]. However, the migrationtrap is not present for the β -disk, simply because there is nosign change in ˙ J mig , I as the surface density decreases mono-tonically with r . In addition, although there is a sign changein ˙ J mig , I for the α -disk, ˙ J wind and ˙ J gw dominate in the regionwhere ˙ J mig , I is positive. As a result, the combined torque neverchanges sign and there is no migration trap in the α -disk ei-ther.To compare with previous studies about migration traps[e.g., 80], we also calculate the disk structure of an compari-son α -disk with M • = × M (cid:12) , ˙ M • = . M Edd and α = . ∼ M • and ∼ M • ) where ˙ J mig , I changes its sign fromnegative to positive in the decreasing r direction. These tworadii are called migration traps by previous studies, and therehave been extensive studies on the consquences of migrationtraps in AGN disks accumulating compact objects [22–30].As shown in Fig. 7, the migration traps are supposed to beovercome by two counteracting processes: GW emission andwind. We have explored the parameter space α ∈ (0 . , . M • ∈ (0 . , .
5) ˙ M Edd • , M • ∈ (10 , ) M (cid:12) , where no migra-tion trap is found in either α -disks or β -disks.In the example TQM disk model of Ref. [21], a salient fea-ture is the presence of a opacity gap and consequently a sharpdensity increase on its inner edge, where the type-I migrationtorque changes sign according to Eq. (30) and has been in-terpreted as a possible location of migration trap [e.g., 80].As mentioned in Section III A, the sharp density increase onthe edge of the opacity gap is in fact resulted by the improperequation of radiation pressure. With a more general equationof radiation pressure, we find the density increase is muchmilder and there is no sign change in the type-I migrationtorque. We also explored the parameter space X ∈ (0 . , . M • ∈ (0 . , .
5) ˙ M Edd • , M • ∈ (10 , ) M (cid:12) , where no migra-tion trap is found in TQM disks.To summarize, we find no migration trap in the three AGNdisk models in a large parameter space we considered. In α -disks, there are locations where the type-I migration torque changes sign, but the total torque is always negative becauseof the negative torque from head wind and GW emission. In β -disks, there is no sign change in the type-I migration torquebecause of the monotonical density and temperature profiles.In TQM disks, there is no sign change in the type-I migrationtorque either as explained above [81].In previous studies of hierarchical binary black hole (BBH)mergers in migration traps of AGN disks, the existence of mi-gration traps was established on a fiducial α -disk model inRef. [20] and a fiducial TQM disk model in Ref. [21]. Asshown above, the migration traps no longer stand after tak-ing account of the head wind and / or using a more reason-able equation of radiation pressure. Therefore, the analysishere raises concerns about the feasibility of hierarchical BBHformation channel in AGN disks. (see [e.g., 82, 83] for theimpact of migration traps to general BBH mergers in AGNdisks).According to the three disk models considered in this work,we do not expect any migration trap, but it does not excludethe possibility that these disk models are not good approxima-tions to the AGN accretion disks in nature. If a migration trapindeed exists in an AGN disk, sBHs would be trapped and hi-erarchical BBH mergers would consequently happen until theremnant BH is so massive that it opens a gap and tears downthe trap. The critical BH mass for gap opening is sensitiveto the trap location and the local disk structure [33]. For ex-ample, in the scenario considered in [32] to associate possibleAGN flares with BBH merger in the disk, the trap is assumedto be ∼ M • away from the MBH with mass M • ∼ M (cid:12) and the the critical BH mass is ∼ O (10 ) M (cid:12) assuming an α -disk with α = . M • = . M Edd • . Thismeans for typical 10 M (cid:12) sBHs, O (10)-times mergers are ex-pected before a gap opens up and the final BH starts its typeII migration. IV. EMRI FORMATION ASSISTED BY AGN ACCRETIONDISKSA. Fokker-Planck equation
As shown in the right panels of Figs. 5, 6, and in Eq. (15),the timescale of orbit eccentricity / inclination decay t wav , asdriven by density wave geneartion, is much shorter than thesBH migration timescale t odmig and the stellar cluster relaxationtimescale t rlx . As a result, one may naively expect that allstars / sBHs are captured onto the disk on the shortest timescale t wav , which turns out to be incorrect. In fact, due to the densedistribution of scatters (dominated by stars) in the disk, a largeportion of stars / sBHs captured by the disk will be scatteredback into the cluster. As demonstrated in previous studies ofstar-disk interactions [84–86], a local equilibrium is built be-tween the net rate of stars captured onto the disk and the rateof inward migration within the disk. Assuming a net fraction µ star of in-cluster stars are captured on to the disk and migrateinward in a timescale t star , idmig , the in-cluster star loss rate due to2disk capture can be formulated as (cid:32) ∂ f star ∂ t (cid:33) cap = − µ star f star t star , idmig . (44)From another aspect, the loss rate should be proportional tothe ratio of two timescales, the inclination damping timescale t starwav and the timescale in which a star is scattered by in-disk scatters (mainly stars), where the latter is inversely pro-portional to the local density of in-disk scatters n star , id , i.e., (cid:16) ∂ f star ∂ t (cid:17) cap ∝ / ( t starwav n star , id ). In the same way, the in-clustersBH loss rate due to disk capture depends on the inclinationdamping timescale t bhwav and the local density of in-disk scat-ters (mainly stars) n star , id . As a result, we obtain1 f bh (cid:32) ∂ f bh ∂ t (cid:33) cap = t starwav t bhwav f star (cid:32) ∂ f star ∂ t (cid:33) cap = − µ star m bh m star t star , idmig , (45)The net fraction µ star should fall in the range ∼ ( h , h being the disk aspect ratio. The exact value of µ star depends onthe detailed balance between the rate of stars captured to thedisk, the fraction of which scattered away from the disk, andthe inward migration rate of in-disk scatters, which requireseperate numerical studies to determine.With the above analysis, we rewrite the orbit-averagedFokker-Planck equation (13) for sBHs / stars in the cluster (ex-cluding in-disk sBHs / stars) as C ∂ f ∂ t = − ∂∂ E F E − ∂∂ R F R + S , (46)where F E , F R are defined in Eq. (14), with the advection co-e ffi cients modified as D E , bh → D E , bh − C Et bh , odmig , D R , bh → D R , bh − C − Rt bh , odwav , D E , star → D E , star − C Et star , odmig , D R , star → D R , star − C − Rt star , odwav , (47)where the corrections are due to interactions with the accre-tion disk and GW emission. The (negative) source term arisesfrom stars / sBHs capture onto the accretion disk, S bh = − µ cap C f bh t star , idmig , S star = − µ cap m star m bh C f star t star , idmig . (48)where we have defined µ cap = µ star m bh / m star , which we expectto be in the range of ∼ ( h , m bh / m star . In this paper, we treat µ cap as a free parameter and take two representative numbers µ cap = µ cap = . T disk , only sBHs within some criticalradius r crit ( T disk ) have enough time to migrate within the disk T disk [yr] E c r i t FIG. 8. The dependence of E crit [ σ ] on the disk life time T disk [yr]for the fiducial α -disk shown in Fig. 5. to reach the central MBH within T disk . Therefore, the EMRIrate is formulated as Γ emri ( t ; T disk ) (cid:39) (cid:90) (cid:90) E > E crit − S bh ( E , R ) dEdR , (49)where E crit ( T disk ) : = φ ( r crit ( T disk )) / f i ( t , E , R ) usingEq. (46) for a period of time T disk . On the low energy E → / star-diskinteractions and two body scatterings are slow, therefore f i ( t , E , R ) | E → = f i ( t = , E , R ) | E → . (50)On the R = F R for both sBHsand stars still applies F R | R → = . (51)On the R = R lc boundary, we again impose the vanishing fluxcondition F R | R = R lc = . (52)This boundary condition is di ff erent from the one imposed inthe no-disk case, simply because the fast eccentricity dampingby density waves drives stars / sBHs away from the loss cone. B. Numerical method
The numerical method for solving Eq. (46) is the sameas the one used for solving Eq. (13) in Section II. The onlyextra numerical subtlety is due to the scale separation: thetimescales of migration and orbit eccentricity decay t mig , t wav are smaller than the cluster relaxation timescale t rlx . As a re-sult, the advection coe ffi cient will be much larger than the dif-fusion coe ffi cient: |D E / E | (cid:29) D EE / E and |D R | (cid:29) D RR . Toavoid numerical di ffi culties for resolving the large scale sep-arations, we choose to regularize the advection coe ffi cients in3 log ( E ) R log ( f star ) log ( E ) R log ( f bh ) E f star ( t = 0) f star ( t = 10 yr) f star ( t = 10 yr) f bh ( t = 0) f bh ( t = 10 yr) f bh ( t = 10 yr) log ( E ) R log ( f star ) log ( E ) R log ( f bh ) E f star ( t = 0) f star ( t = 10 yr) f star ( t = 10 yr) f bh ( t = 0) f bh ( t = 10 yr) f bh ( t = 10 yr) FIG. 9. We show the distribution functions for the case of fast disk capture ( µ cap =
1) in the upper row: f star ( t = yr , E , R ) (left panel), f bh ( t = yr , E , R ) (middle panel), the time dependence of EMRI rate Γ emri , where f i are shown in units of 10 pc − / (2 πσ ) / and E is shownin units of σ . The counterparts of the slow disk capture ( µ cap = .
1) case are in the lower row.
Eq. (47) as follows: D E , bh → D E , bh − C Et bh , odmig + (cid:15) T , D R , bh → D R , bh − η (1 − R ) D , D E , star → D E , star − C Et bh , odmig + (cid:15) T m star m bh , D R , star → D R , star − η (1 − R ) D m star m bh , (53)where (cid:15) is a small number ensuring a numerically resolvablescale separation, T =
10 Gyr, η is a large number ensuringa large scale separation between the regularized advection co-e ffi cient D R and the di ff usion coe ffi cient D RR , and D ( E ) isdefined as the maximal value of D RR , bh ( E , R ) for given E .We choose (cid:15) = − and η = as default regulariza-tion parameters. In Appendix B, we will show that di ff erentchoices of these two regularization parameters do not a ff ectthe EMRI rate assisted by accreton disks as long as (cid:15) is su ffi -ciently small. C. Comparison with loss-cone rate
In order to compare with the canonical EMRI rate asso-ciated with the loss cone mechanism, we consider a fiducialmodel assuming the same MBH + sBH / star cluster model andparameters as those considered in Section II and shown inFig. 5. We take the final state of Eq. (13) (Fig. 1) as the initial condition of Eq. (46), and evolve the equation for a period oftime T disk .In Fig. 9, we show the sBH / star distribution functions f i ( t , E , R ). In combination with the initial condition (Fig. 1),we see that sBHs are migrating toward larger R and larger E driven by the density waves. In Fig. 10, we show the timedependence of the disk assisted EMRI rates for two di ff erentcases: µ cap = µ cap = .
1. For both cases, we find that thedependence of Γ emri ( t ; T disk ) on T disk is relatively weak (the Γ emri ( t ; T disk ) curves for di ff erent T disk roughtly overlap witheach other). This is because of the weak dependence of E crit on T disk (Fig. 8). For short time t , Γ emri mainly depends onthe initial distributions f i ( t = , E , R ) and is proportional tothe parameter µ cap . For long time t , we expect an equilibriumbetween sBHs migrating inward from E < E crit and sBHs cap-tured onto the disk in the region E > E crit , where the EMRIrate is determined by the sBH supply rate irrespective of µ cap .For the case of µ cap =
1, we find that the initial EMRI rate perMBH is Γ emri ( t = T disk ) ∼ × Gyr − (irrespective of disklife time); after a rapid initial settling (see Fig. 11), Γ emri then( t (cid:46)
20 Myr) slowly increases as more sBHs migrate inwardfrom E < E crit than those captured onto the disk; after t ∼ Γ emri slowly decreases due tothe decreasing supply arising from the decreasing number ofsBHs around E ∼ E crit . For µ cap = .
1, the evolution track issimilar except starting with a lower initial EMRI rate.In both examples we see that the disk-assisted rate perMBH is much higher ( O (10 − )) than the loss-cone rate,4 T disk [Myr] e m r i [ G y r ] cap = 1 cap = 0.1 T disk = 10 yr T disk = 10 yr T disk = 10 yrNo Disk FIG. 10. The time dependence of disk assisted EMRI rate per MBH Γ emri ( t ; T disk ) for di ff erent disk life times T disk and di ff erent disk cap-ture parameters µ cap . The dashed line is the average EMRI rate vialoss cone of the fiducial model. t [Myr] e m r i [ G y r ] M = 10 MM = 4 × 10 MM = 10 MM = 4 × 10 MM = 10 M FIG. 11. The disk-assisted EMRI rate Γ emri ( t ; T disk = yrs) fordi ff erent MBH masses M • , where take µ cap = which is due to much more e ffi cient capture and transportmechanisms by the disk. It turns out that this observation isgenerally true as we vary the parameters for MBH, disk andstar cluster. In the following subsection we explore the depen-dence of disk-assisted EMRI rate on various model parame-ters. D. EMRI rate for di ff erent models In this subsection, we investigate the dependence of thedisk-assisted EMRI rate on di ff erent parameters: the MBHmass M • , the accretion rate ˙ M • and the α parameter of α -diskmodels. We also explore the EMRI rate for the TQM diskmodel and di ff erent cluster initial condition. We take the fidu-cial model M • = × M (cid:12) , ˙ M • = . M Edd • , α = . , µ cap = M • , we initialize the distributions of surrounding stars and sBHs according to the Tremaine’scluster model outlined in Section II B 1, and evolve the distri-butions using the Fokker-Planck equation (13) to t = t = yr. In Fig. 11, we show the EMRIrate Γ emri ( t ; T disk = yrs) for each M • (the results of di ff er-ent T disk and µ cap can be easily inferred from reading Fig. 10).The initial EMRI rates are ∼ O (10 ) Gyr − , except for thelow mass M • = M (cid:12) case, for which the high EMRI ratevia loss cone has consumed sBHs close to the central MBH.Their subsequent evolution basically follow the descriptiongiven in the Sec. IV C: an increasing phase where the supplyfrom inward migration dominates and then a decreasing phasewhere the supply-consumption equilibrium is built. The rel-evant timescale is longer for larger M • . As a result, Γ emri isroughtly a constant in 10 yrs for M • = M (cid:12) , while Γ emri changes by orders of magnitude for low MBH massses. Inconstrast with the loss-cone EMRI channel, the average disk-assisted EMRI rate (for long disk life time) Γ emri increaseswith the MBH mass M • , because the capacity of sBH reserv-ior ( ∝ M • ) is larger for more massive MBHs.We conclude this section by collecting the average disk-assisted EMRI rate (cid:104) Γ emri ( T disk ) (cid:105) : = T disk (cid:90) T disk Γ emri ( t ; T disk ) dt (54)of the representative models in Table I, and briefly outline itsparameter dependence as follows. Parameter µ cap . The average EMRI rate (cid:104) Γ emri ( T disk (cid:105) ) isshould be proportional to the parameter µ cap for T disk →
0, while is independent of µ cap for long T disk as explainedabove. For the fiducial model with M • = × M (cid:12) , (cid:104) Γ emri ( T disk = yrs) (cid:105) ∼ O (10 ) Gyr − (irrespective ofthe parameter µ cap and the α -disk model used), which ishigher than the average loss-cone EMRI rate ¯ Γ emri by a fac-tor of ∼ O (10 ) (Fig. 4), and (cid:104) Γ emri ( T disk = yrs (cid:105) ) (cid:38) O (10 µ cap ) Gyr − is higher by a factor of (cid:38) µ cap . MBH mass M • . The average EMRI rate (cid:104) Γ emri ( T disk ) (cid:105) with long T disk is usually higher for a heavier MBH be-cause more sBHs are available in the stellar cluster ( ∝ M • )while (cid:104) Γ emri ( T disk ) (cid:105) with short T disk is more initial conditiondependent. In a model same to the fiducial model exceptwith a lighter MBH M • = M (cid:12) , (cid:104) Γ emri ( T disk = yrs) (cid:105) ∼O (10 ) Gyr − (irrespective of µ cap and the disk model used),which is higher than the average loss-cone EMRI rate ¯ Γ emri by a factor of ∼
15, and (cid:104) Γ emri ( T disk = yrs (cid:105) ) ∼ (2 − × Gyr − (depending on the value of µ cap ) is higher by a fac-tor of (30 − Disk model.
The dependence of disk-assisted EMRI rateon the α viscosity parameter and on the accretion rate ˙ M • areweak, which can be understood from their impact on the mi-gration timescales of stars and sBHs. Accretion disks withsmaller α parameters / higher accretion rates are thicker (higherdisk aspect ratio h ) and therefore the migration timescalesarising from density waves are longer [Eq. (31)], so that theEMRI rate is lower for a disk with a smaller α and a higher ˙ M .Note that this trend does not sustain all the way to the regime5 TABLE I. Average EMRI rate per MBH (cid:104) Γ emri ( T disk ) (cid:105) [Gyr − ] for di ff erent models. In the 1st column are the parameters of initial stellar clusterprofiles, in the 2nd / M • and the parameter µ cap , and in the 4th column is the α parameter in (default) α -disksor X parameter in TQM disks, and in the 5th column is the MBH accretion rate. In a few models with TQM disks, the EMRI rates are nearlyzero for short disk lifetime T disk because the migration timescale in inner parts of the disk is longer than T disk , consequently almost no sBHsuccessfully migrates into the MBH within T disk .( γ, δ ) M • / M (cid:12) µ cap α or X ˙ M • / ˙ M Edd • (cid:104) Γ emri ( T disk = yrs) (cid:105) (cid:104) Γ emri ( T disk = yrs) (cid:105) (cid:104) Γ emri ( T disk = yrs) (cid:105) (1 . , . × − − . × . × . × × . × . × . × × . × . × . × × . × . × . × × . × . × . × × − − − . × . × . × × . × . × . × × . × . × . × × . × . × . × × . × . × . × × − × − . × . × . × − − . × . × . × × − − . × . × . × − − . × . × . × × − ) TQM − ∼ ∼ . × × ∼ . × . × × ∼ . × . × × . × . × . × × . × . × . × (1 . , . × − − . × . × . × × . × . × . × × . × . × . × × . × . × . × × . × . × . × (1 . , . × − − . × . × . × × . × . × . × × . × . × . × × . × . × . × × . × . × . × of extremely low accretion rates (cid:46) α ˙ M Edd • [87], where thethin accretion disk description breaks down. But (cid:104) Γ emri ( T disk (cid:105) )is sensitive to which disk model is assumed ( α -disk or TQMdisk). Due to the large di ff erence in the two disk models, thesBH migration timescales t bh , idmig are quite di ff erent in these twodisks, and consequently the critical radii r crit ( T disk ) [Eq. (49)]in these two disks di ff er significantly for short T disk , where r crit ( T disk = yrs) ∼ M • for α disks and r crit ( T disk = yrs) ∼ M • for TQM disks (see Figs. 5 and 6). This ex-plains why (cid:104) Γ emri ( T disk = yrs) (cid:105) is much smaller for TQMdisks. However, we notice that the α -viscosity prescriptionis more favored by the current knowledge of turbulence vis-cosities driven by magnetorotational instability in inner partsof AGN disks. The sharp di ff erence in (cid:104) Γ emri ( T disk = yrs) (cid:105) for TQM disks is likely an artifact of too e ffi cient angular mo-mentum transport assumed in TQM disks. Intial density profile of stellar cluster.
The linear depen-dence of (cid:104) Γ emri ( T disk ) (cid:105) on the number fraction of sBHs δ is nat-ural for long T disk , while the dependence is more complicatefor short T disk because the sBH fraction has been substantiallychanged in the pre-disk evolution. The dependence on the ini-tial stellar density profile (parameterized by the power index γ ) varies for di ff erent MBH masses. For the Galactic nuclear stellar cluster like model ( γ = . (cid:104) Γ emri ( T disk ) (cid:105) is higher in the case of heavy MBHs (irrespective of T disk ). Forlighter MBHs, the dependence reverses because more sBHshave been depleted via the loss cone in the pre-disk phaseand therefore less sBHs are available for disk capture. In thismodel with a steeper intial stellar density profile, the averageEMRI rate (cid:104) Γ emri ( T disk ) (cid:105) changes only by a factor of few (cid:46) (cid:104) Γ emri ( T disk = yrs) (cid:105) is higher than the average loss-conerate by O (10 − ) for M • = × M (cid:12) and by O (10 − )for M • = M (cid:12) (irrespective of disk models, initial stellarprofiles and the value of parameter µ cap ). In the case of shortdisk lifetime T disk = yrs, (cid:104) Γ emri ( T disk ) (cid:105) depends on thedisk model used and the value of parameter µ cap : (cid:104) Γ emri ( T disk ) (cid:105) is higher than the loss-cone rate by O (10 µ cap ) for M • = × M (cid:12) and is higher by O (10 µ cap ) for M • = M (cid:12) as-suming an α/β -disk, while (cid:104) Γ emri ( T disk ) (cid:105) is much smaller as-suming a TQM disk.6 V. SUMMARY AND DISCUSSIONA. Model uncertainties
There are several caveats in this analysis that possibly a ff ectthe estimate for disk-assisted EMRI rates. Disk model . In this work, we have used the α , β andTQM model to describe the profiles of accretion disks aroundMBHs. The di ff erence between α -disk model and β -diskmodel mainly comes from the prescription of modelling diskviscosities, where the β -disk viscosity prescription was intro-duced to avoid the α -disk thermally instability in the radia-tion pressure dominated region [88, 89]. We also note thatthese two models yiels similar disk structures for larger dis-tance r (cid:38) M • . On the other hand, the TQM model was de-veloped for consistently modelling the star formation in AGNdisks and the disk structure in large distance. In partiuclar,we point out that in Ref. [21], an equation of radiation pres-sure which holds only in optically thick regime was used inboth optically thick and thin regimes, which gives rises to anartificially large density variation on the edge of the opacitygap. Although these models are the state-of-the-art tools fordescribng the AGN accretion disks, there is still large room forimprovement before accurately describing the reality, becauseof various simplification / approximations taken in the modelsand the complexity of disk conditions / states in nature. There-fore the real disk profiles may or may not be accurately de-scribed by these models, which is a possible source of uncer-tainty in the rate analysis. From α -disks and TQM disks, wefind di ff erent disk models a ff ect the disk-assisted EMRI ratebecause of the di ff ference in the migration timescales whichmildly a ff ect the EMRI rate for long disk lifetime but changethe EMRI rate hugely for short disk lifetime ∼ yrs becausethe migration timescale is longer in inner parts of TQM disks(Fig. 6) and only nearby sBHs ( r (cid:46) M • ) can migrate intothe MBH within the disk lifetime. Based on the current under-standing of turbulence viscosities driven by magnetorotationalinstability in fully ionized accretion disks, the α -viscosity isa good approximation. Therefore the sharp di ff erence in thedisk-assisted EMRI rate for short disk lifetime in TQM diskmodels is likely the consequence of the artificial angular mo-mentum transport assumed in TQM disks. Repeating AGNs . Disk lifetime characterizes the time dura-tion of the active phase of AGNs. However, it is possible thatbefore the current active phase, there are already a sequence ofactive phases of AGN, with various lifetimes [90, 91]. Theseactive periods may introduce significant change in star clus-ter distributions due to disk-star / sBH interactions. In fact, ifwe neglect the evolution of the star cluster distribution duringthe “quiet” periods between those active periods, the presenceof previous active periods e ff ectively extends the disk lifetimefrom T disk to T disk + t in Eq. (49) and shifts t to t + t in Fig. 11,where t is the summation of the lifetimes of all previous ac-tive cycles. In other words, the integration upper and lowerlimit may also need to be shifted by t in Eq. (54). In addition,the distribution of stars and sBHs may still evolve during thequiet phase, which further complicates the picture. To fullyaccount for these e ff ects, we will need (from observations) in- formation about the fraction of active cycles v.s. quite cyclesfor AGNs and the total duration of AGN outside which there isno more active phase. In this work, we use T disk ∈ (10 − )yrs as examples, while the total duration of AGN active phases(e ff ective disk lifetime) should be longer ∼ (10 − ) yrs ac-cording to Soltan’s argument [92]. Initial condition of stellar clusters.
We initialized the stel-lar cluster following the Tremaine’s cluster model, and ex-plored the dependence of both the loss-cone EMRI rate andthe disk-assisted EMRI rate on the sBH fraction δ and the den-sity profile (parameterized by γ ). We find the loss-cone EMRIrate dependence on δ is shallower than linear scaling and thedisk-assisted EMRI rate dependence is linear (for long disklifetime). Di ff erent initial density profiles a ff ect the loss-coneEMRI rate by changing the total number of stars / sBHs withinthe influence sphere N star , bh ( r < r h ). After a few Gyrs, thedistributions within the influence sphere have reached a lo-cal equilibrium and the details of initial distributions has beenmostly erased except the total number of stars / sBHs. With theaccretion disk turned on, the disk-assisted EMRI rate againroughly only depends on the total number instead of othererased details of the initial distributions. Therefore we do notexpect much uncertainty in the EMRI rate estimation arisingfrom unknowns in the initial condition of stellar clusters ex-cept the total number of stars within the influence N star ( r < r h )which can be inferred from the MBH mass as ∼ M • / M (cid:12) . Torque for inclined orbits . Based on the studies for plane-tary systems, a point mass moving along an inclined orbit withrespect to a disk excites density waves of various kinds thatmodify the orbit period, eccentricity and inclination in time[16–19]. However, we notice that these studies mainly focuson low-inclinnation orbit. For highly inclined orbits, whilethe qualitative density wave generation and propagation pic-ture should still apply, the actual torque may deviate from theformulas derived or fitted for low-inclination orbits. The mi-gration speed of sBHs is proportional to the magnitude of themigration torque. The disk-assisted EMRI rate for short disklifetime is determined by the capture rate of sBHs within thecritical radius r crit , so that it is insenstive to the torque; for longdisk lifetime the rate is determined by the migration supply-capture consumption equilibrium, so that it should be propor-tional to the torque magnitude. The uncertainty in the torquemagnitude should proportionally propagate to the EMRI rateestimation. B. Application and future work
With the EMRI rate computed for di ff erent models, thenext natural step is to predict the corresponding event rate forspace-borne detectors such as LISA and TianQin, based onmass distribution of MBH, star cluster distribution, disk pa-rameters and detector sensitivity. We will leave this part asfuture work. It is however evident from the rates listed in Ta-ble I, the loss-cone rate (Fig. 4) and the AGN fraction [35, 36],that the disk-assisted EMRIs should be a good fraction of allEMRIs detected by LISA and TianQin.It is then important to explore how to distinguish disk-7assisted EMRIs from loss-cone EMRIs within future obser-vations. Based on the analysis in [12], the eccentricity of loss-cone EMRIs ranges from 0 to 0 . . ff ectively zero eccentricity considering that theeccentricity damping timescale is much shorter than the mi-gration timescale, and the inclination should be confined bythe disk thickness ι (cid:46) h .As a good fraction of EMRIs detected by LISA shouldcome from systems with AGN, it is possible that the elec-tromagnetic emission from some of these AGNs can be ob-served. This brings up the opportunity for multi-messengeranalysis for these EMRIs. According to Ref. [9], a fractionof low-redshift ( z (cid:46) .
3) EMRIs can be traced back to theirhost galaxies with LISA observations alone, and host galax-ies of ∼
50% EMRIs in low-redshift ( z (cid:46) .
5) AGNs canbe identified with LISA observations alone considering thelower density of AGNs. If the host galaxy of such EMRIcan be identified, the distance measurement from gravitaitonalwave observable and redshift measurement from optical ob-servables should allow accurate determination on the Hub-ble’s constant. On the other hand, for those distant EMRIswithout host galaxy identification, one may still be able tomeasure the Hubble’s constant using all AGNs in the errorvolume with the same statistical method introduced in [93].In addition, for certain disk profiles the EMRI waveform maybe significantly modified, so that certain disk properties areable to be constrained with GW observations [7, 33]. Theseinformation can be further compared with electromagnetic ob-servables from the AGN to help reveal the unknowns aboutaccretion physics.Lastly, in [6] we observed that a pair of sBHs embedded inthe acretion disk may be locked into mean motion resonanceand then migrate together towards the MBH. The resonancebreaks when the pair is close to the MBH, at which stage theinner EMRI should be a ff ected by the gravitational force fromthe outer sBH, so that the gravitaitonal waveform should becorrespondingly satisfied [5]. Such resonance locking for apair or a chain of objects has been discussed previously forplanetary systems [94], which is also interesting to explore inthis disk-assisted EMRI scenario. ACKNOWLEDGMENTS
Z.P. and H.Y. thank B´eatrice Bonga for instructive discus-sions during early stage of this work. Z.P. and H.Y. also thankXinyu Li, Cole Miller, Neal Dalal and Barry McKernan forvery helpful discussions and comments. Z. P. and H. Y. aresupported by the Natural Sciences and Engineering ResearchCouncil of Canada and in part by Perimeter Institute for Theo-retical Physics. Research at Perimeter Institute is supported inpart by the Government of Canada through the Department ofInnovation, Science and Economic Development Canada andby the Province of Ontario through the Ministry of Collegesand Universities.
Appendix A: Di ff usion and advection coe ffi cients in theFokker-Planck equation (13) Following Ref. [50], we extend the calculation of the dif-fusion and the advection coe ffi cients of a single-componentcluster in Ref. [38, 44, 45] to our two-component (stars andsBHs) case. We first define a few auxiliary functions: F ( i )0 ( E , r ) = (4 π ) m i ln Λ (cid:90) E −∞ dE (cid:48) ¯ f i ( E (cid:48) ) , F ( i )1 ( E , r ) = (4 π ) m i ln Λ (cid:90) φ ( r ) E dE (cid:48) (cid:32) φ − E (cid:48) φ − E (cid:33) / ¯ f i ( E (cid:48) ) , F ( i )2 ( E , r ) = (4 π ) m i ln Λ (cid:90) φ ( r ) E dE (cid:48) (cid:32) φ − E (cid:48) φ − E (cid:33) / ¯ f i ( E (cid:48) ) , (A1)where i = { star , bh } , ln Λ the Coulomb’s logarithm which takeas ln Λ =
10, and ¯ f i ( E ) : = (cid:90) f ( E , R ) dR . (A2)With these auxiliary functions, the coe ffi cients are written as D ( i ) EE = π J c (cid:90) r + r − drv r v ( F ( i )0 + F ( i )2 ) + ( i ↔ j ) , D ( i ) E = − π J c (cid:90) r + r − drv r v F ( i )1 + m i m j × ( i ↔ j ) , D ( i ) ER = π J (cid:90) r + r − drv r (cid:32) v v c − (cid:33) ( F ( i )0 + F ( i )2 ) + ( i ↔ j ) , D ( i ) RR = π R (cid:90) r + r − drv r (cid:40) r v v t (cid:32) v v c − (cid:33) + v r F ( i )0 + r v v r F ( i )1 + r v v t (cid:32) v v c − (cid:33) − v r F ( i )2 (cid:41) + ( i ↔ j ) , D ( i ) R = − π Rr c (cid:90) r + r − drv r (cid:32) − v c v (cid:33) F ( i )1 + m i m j × ( i ↔ j ) , (A3)where j = { star , bh } , i (cid:44) j , and v t = J / r is the tangentialvelocity. Appendix B: Sanity check for the regularization algorithm
As discussed in Section IV B, the scale separations betweendi ff erent timescales can be as large as 10 orders of magni-tude. Therefore we need to the regularize the advection coef-ficients ensuring the scale separations are numerically resolv-able. There are two parameters our regularization algorithm[Eq. (53)]: (cid:15) and η , where (cid:15) is a small number determiningan constant floor of migrate timescale as (cid:15) T , and η is a largenumber determining a numerically resolvable scale separationof D RR and D R . In the main text, we choose (cid:15) = − and η = as our default values of the two regularization param-eters. In Fig. 12, we show the time dependence of EMRI rateof the fiducial fast-disk-capture model with di ff erent (cid:15) and η ,8 t [Myr] e m r i [ G y r ] = 10 , = 10 = 2 × 10 , = 2 × 10 FIG. 12. The EMRI rate Γ ( t ; T disk = yrs) of the fiducial modelis independent of the values of regularization parameters (cid:15) and η weused in the main text. where we see the EMRI rate has no dependence on the twoparameters. [1] J. Baker, J. Bellovary, P. L. Bender, E. Berti, R. Caldwell,J. Camp, J. W. Conklin, N. Cornish, C. Cutler, R. DeRosa,M. Eracleous, E. C. Ferrara, S. Francis, M. Hewitson,K. Holley-Bockelmann, A. Hornschemeier, C. Hogan, B. Ka-mai, B. J. Kelly, J. Shapiro Key, S. L. Larson, J. Livas, S. Man-thripragada, K. McKenzie, S. T. McWilliams, G. Mueller,P. Natarajan, K. Numata, N. Rioux, S. R. Sankar, J. Schnittman,D. Shoemaker, D. Shoemaker, J. Slutsky, R. Spero, R. Steb-bins, I. Thorpe, M. Vallisneri, B. Ware, P. Wass, A. Yu,and J. Ziemer, arXiv e-prints , arXiv:1907.06482 (2019),arXiv:1907.06482 [astro-ph.IM].[2] J. Mei, Y.-Z. Bai, J. Bao, E. Barausse, L. Cai, E. Canuto,B. Cao, W.-M. Chen, Y. Chen, Y.-W. Ding, H.-Z. Duan,H. Fan, W.-F. Feng, H. Fu, Q. Gao, T. Gao, Y. Gong, X. Gou,C.-Z. Gu, D.-F. Gu, Z.-Q. He, M. Hendry, W. Hong, X.-C.Hu, Y.-M. Hu, Y. Hu, S.-J. Huang, X.-Q. Huang, Q. Jiang,Y.-Z. Jiang, Y. Jiang, Z. Jiang, H.-M. Jin, V. Korol, H.-Y.Li, M. Li, M. Li, P. Li, R. Li, Y. Li, Z. Li, Z. Li, Z.-X. Li,Y.-R. Liang, Z.-C. Liang, F.-J. Liao, Q. Liu, S. Liu, Y.-C. Liu,L. Liu, P.-B. Liu, X. Liu, Y. Liu, X.-F. Lu, Y. Lu, Z.-H. Lu,Y. Luo, Z.-C. Luo, V. Milyukov, M. Ming, X. Pi, C. Qin, S.-B.Qu, A. Sesana, C. Shao, C. Shi, W. Su, D.-Y. Tan, Y. Tan,Z. Tan, L.-C. Tu, B. Wang, C.-R. Wang, F. Wang, G.-F. Wang,H. Wang, J. Wang, L. Wang, P. Wang, X. Wang, Y. Wang, Y.-F.Wang, R. Wei, S.-C. Wu, C.-Y. Xiao, X.-S. Xu, C. Xue, F.-C.Yang, L. Yang, M.-L. Yang, S.-Q. Yang, B. Ye, H.-C. Yeh,S. Yu, D. Zhai, C. Zhang, H. Zhang, J.-d. Zhang, J. Zhang,L. Zhang, X. Zhang, X. Zhang, H. Zhou, M.-Y. Zhou, Z.-B.Zhou, D.-D. Zhu, T.-G. Zi, and J. Luo, Progress of Theoreticaland Experimental Physics (2020), 10.1093 / ptep / ptaa114,ptaa114, https: // academic.oup.com / ptep / advance-article-pdf / doi / / ptep / ptaa114 / / ptaa114.pdf.[3] P. Amaro-Seoane, H. Audley, S. Babak, J. Baker, E. Ba-rausse, P. Bender, E. Berti, P. Binetruy, M. Born, D. Bortoluzzi,J. Camp, C. Caprini, V. Cardoso, M. Colpi, J. Conklin, N. Cor-nish, C. Cutler, K. Danzmann, R. Dolesi, L. Ferraioli, V. Fer-roni, E. Fitzsimons, J. Gair, L. Gesa Bote, D. Giardini, F. Gib-ert, C. Grimani, H. Halloin, G. Heinzel, T. Hertog, M. He-witson, K. Holley-Bockelmann, D. Hollington, M. Hueller,H. Inchauspe, P. Jetzer, N. Karnesis, C. Killow, A. Klein,B. Klipstein, N. Korsakova, S. L. Larson, J. Livas, I. Lloro, N. Man, D. Mance, J. Martino, I. Mateos, K. McKenzie, S. T.McWilliams, C. Miller, G. Mueller, G. Nardini, G. Nelemans,M. Nofrarias, A. Petiteau, P. Pivato, E. Plagnol, E. Porter,J. Reiche, D. Robertson, N. Robertson, E. Rossi, G. Russano,B. Schutz, A. Sesana, D. Shoemaker, J. Slutsky, C. F. Sopuerta,T. Sumner, N. Tamanini, I. Thorpe, M. Troebs, M. Vallisneri,A. Vecchio, D. Vetrugno, S. Vitale, M. Volonteri, G. Wanner,H. Ward, P. Wass, W. Weber, J. Ziemer, and P. Zweifel, arXive-prints , arXiv:1702.00786 (2017), arXiv:1702.00786 [astro-ph.IM].[4] S. Babak, J. Gair, A. Sesana, E. Barausse, C. F. Sopuerta, C. P.Berry, E. Berti, P. Amaro-Seoane, A. Petiteau, and A. Klein,Phys. Rev. D , 103012 (2017), arXiv:1703.09722 [gr-qc].[5] B. Bonga, H. Yang, and S. A. Hughes, Phys. Rev. Lett. ,101103 (2019), arXiv:1905.00030 [gr-qc].[6] H. Yang, B. Bonga, Z. Peng, and G. Li, Phys. Rev. D ,124056 (2019), arXiv:1910.07337 [gr-qc].[7] E. Barausse, V. Cardoso, and P. Pani, Phys. Rev. D , 104059(2014), arXiv:1404.7149 [gr-qc].[8] E. Berti and M. Volonteri, Astrophys. J. , 822 (2008),arXiv:0802.0025 [astro-ph].[9] Z. Pan and H. Yang, Astrophys. J. , 163 (2020),arXiv:2007.03783 [astro-ph.CO].[10] P. Amaro-Seoane, Living Reviews in Relativity , 4 (2018),arXiv:1205.5240 [astro-ph.CO].[11] J. R. Gair, S. Babak, A. Sesana, P. Amaro-Seoane, E. Barausse,C. P. Berry, E. Berti, and C. Sopuerta, J. Phys. Conf. Ser. ,012021 (2017), arXiv:1704.00009 [astro-ph.GA].[12] S. Babak, J. Gair, A. Sesana, E. Barausse, C. F. Sopuerta,C. P. L. Berry, E. Berti, P. Amaro-Seoane, A. Petiteau, andA. Klein, Phys. Rev. D , 103012 (2017), arXiv:1703.09722[gr-qc].[13] M. Preto and P. Amaro-Seoane, Astrophys.J.Lett. , L42(2010), arXiv:0910.3206 [astro-ph.GA].[14] P. Amaro-Seoane and M. Preto, Classical and Quantum Gravity , 094017 (2011), arXiv:1010.5781 [astro-ph.CO].[15] H.-M. Fan, Y.-M. Hu, E. Barausse, A. Sesana, J.-d. Zhang,X. Zhang, T.-G. Zi, and J. Mei, Phys. Rev. D , 063016(2020), arXiv:2005.08212 [astro-ph.HE].[16] P. Goldreich and S. Tremaine, Astrophys. J. , 857 (1979).[17] P. Goldreich and S. Tremaine, Astrophys. J. , 425 (1980). [18] H. Tanaka, T. Takeuchi, and W. R. Ward, Astrophys. J. ,1257 (2002).[19] H. Tanaka and W. R. Ward, Astrophys. J. , 388 (2004).[20] E. Sirko and J. Goodman, MNRAS , 501 (2003),arXiv:astro-ph / , 167 (2005), arXiv:astro-ph / , 460 (2012), arXiv:1206.2309 [astro-ph.GA].[23] B. McKernan, K. E. S. Ford, B. Kocsis, W. Lyra, and L. M.Winter, MNRAS , 900 (2014), arXiv:1403.6433 [astro-ph.GA].[24] N. C. Stone, B. D. Metzger, and Z. Haiman, MNRAS , 946(2017), arXiv:1602.04226 [astro-ph.GA].[25] I. Bartos, B. Kocsis, Z. Haiman, and S. M´arka, Astrophys. J. , 165 (2017), arXiv:1602.03831 [astro-ph.HE].[26] B. McKernan, K. E. S. Ford, J. Bellovary, N. W. C. Leigh,Z. Haiman, B. Kocsis, W. Lyra, M. M. Mac Low, B. Metzger,M. O’Dowd, S. Endlich, and D. J. Rosen, Astrophys. J. , 66(2018), arXiv:1702.07818 [astro-ph.HE].[27] Y. Yang, I. Bartos, V. Gayathri, K. E. S. Ford, Z. Haiman,S. Klimenko, B. Kocsis, S. M´arka, Z. M´arka, B. McKernan,and R. O’Shaughnessy, Phys. Rev. Lett. , 181101 (2019),arXiv:1906.09281 [astro-ph.HE].[28] Y. Yang, I. Bartos, Z. Haiman, B. Kocsis, Z. M´arka, N. C. Stone,and S. M´arka, Astrophys. J. , 122 (2019), arXiv:1903.01405[astro-ph.HE].[29] A. Secunda, J. Bellovary, M.-M. Mac Low, K. E. S. Ford,B. McKernan, N. W. C. Leigh, W. Lyra, and Z. S´andor, As-trophys. J. , 85 (2019), arXiv:1807.02859 [astro-ph.HE].[30] A. Secunda, J. Bellovary, M.-M. Mac Low, K. E. S. Ford,B. McKernan, N. W. C. Leigh, W. Lyra, Z. S´andor, andJ. I. Adorno, Astrophys. J. , 133 (2020), arXiv:2004.11936[astro-ph.HE].[31] B. McKernan, K. E. S. Ford, and R. O’Shaughnessy, MNRAS , 4088 (2020), arXiv:2002.00046 [astro-ph.HE].[32] M. J. Graham et al. , Phys. Rev. Lett. , 251102 (2020),arXiv:2006.14122 [astro-ph.HE].[33] B. Kocsis, N. Yunes, and A. Loeb, Phys. Rev. D , 024032(2011), arXiv:1104.2322 [astro-ph.GA].[34] S. Tremaine, D. O. Richstone, Y.-I. Byun, A. Dressler, S. M.Faber, C. Grillmair, J. Kormendy, and T. R. Lauer, Astronomi-cal Journal , 634 (1994), arXiv:astro-ph / , 1309 (2009),arXiv:0901.1109 [astro-ph.CO].[36] M. Macuga, P. Martini, E. D. Miller, M. Brodwin, M. Hayashi,T. Kodama, Y. Koyama, R. A. Overzier, R. Shimakawa, K.-i. Tadaki, and I. Tanaka, Astrophys. J. , 54 (2019),arXiv:1805.06569 [astro-ph.GA].[37] S. McGee, A. Sesana, and A. Vecchio, Nature Astronomy ,26 (2020), arXiv:1811.00050 [astro-ph.HE].[38] S. L. Shapiro and A. B. Marchant, Astrophys. J. , 603(1978).[39] C. Hopman and T. Alexander, Astrophys. J. , 362 (2005),arXiv:astro-ph / , 129 (2016),arXiv:1508.01390 [astro-ph.GA].[41] C. Cutler, D. Kennefick, and E. Poisson, Phys. Rev. D , 3816(1994).[42] A. P. Lightman and S. L. Shapiro, Astrophys. J. , 244(1977). [43] A star orbiting around a MBH will be disrupted as long as itsperiapsis is r p (cid:46) r star ( M • / m star ) / , where r star ∼ km isthe typical star radius. As a result, we find J lc , star ( E (cid:39) = M • (cid:112) a (1 − e ) / M • = M • (cid:112) r p (1 + e ) / M • (cid:39) M • (cid:112) r p / M • , wherewe have used Eq. (1) in the first equality, the relation r p = a (1 − e ) in the second and e (cid:39) , 1087 (1978).[45] H. Cohn, Astrophys. J. , 1036 (1979).[46] W. Dehnen, MNRAS , 250 (1993).[47] R. Launhardt, R. Zylka, and P. G. Mezger, A&A , 112(2002), arXiv:astro-ph / , A47(2014), arXiv:1403.6657 [astro-ph.GA].[49] D. F. Cherno ff and M. D. Weinberg, Astrophys. J. , 121(1990).[50] J. Binney and S. Tremaine, Galactic dynamics (1987).[51] J. Spitzer, Lyman and M. H. Hart, Astrophys. J. , 399(1971).[52] J. N. Bahcall and R. A. Wolf, Astrophys. J. , 214 (1976).[53] In principle, we should evolve both the distributions f i and thepotential field φ ( r ) self-consistently. For the problem we arediscussing, the potential field φ ( r ) barely changes [13] becauseduring the evolution time range the distributions evolve mainlywithin the influence radius where the potential field is domi-nated by the MBH.[54] S. Tremaine, K. Gebhardt, R. Bender, G. Bower, A. Dressler,S. M. Faber, A. V. Filippenko, R. Green, C. Grillmair, L. C.Ho, J. Kormendy, T. R. Lauer, J. Magorrian, J. Pinkney,and D. Richstone, Astrophys. J. , 740 (2002), arXiv:astro-ph / , 198 (2009),arXiv:0903.4897 [astro-ph.GA].[56] T. Alexander and C. Hopman, Astrophys. J. , 1861 (2009),arXiv:0808.3150 [astro-ph].[57] S. Chandrasekhar, Astrophys. J. , 255 (1943).[58] E. C. Ostriker, Astrophys. J. , 252 (1999), arXiv:astro-ph / , 4204 (2017), arXiv:1708.09807[astro-ph.EP].[60] A. M. Hankla, Y.-F. Jiang, and P. J. Armitage, Astrophys. J. , 50 (2020), arXiv:2005.03785 [astro-ph.EP].[61] N. I. Shakura and R. A. Sunyaev, A&A , 33 (1973).[62] D. R. Alexander and J. W. Ferguson, Astrophys. J. , 879(1994).[63] C. A. Iglesias and F. J. Rogers, Astrophys. J. , 943 (1996).[64] S. A. Balbus and J. F. Hawley, Astrophys. J. , 214 (1991).[65] S. A. Balbus and J. F. Hawley, Reviews of Modern Physics ,1 (1998).[66] R. G. Martin, C. J. Nixon, J. E. Pringle, and M. Livio, NewAstronomy , 7 (2019), arXiv:1901.01580 [astro-ph.HE].[67] S. J. Paardekooper, C. Baruteau, A. Crida, and W. Kley, MN-RAS , 1950 (2010), arXiv:0909.4552 [astro-ph.EP].[68] R. P. Nelson, A&A , 1067 (2005), arXiv:astro-ph / ,1413 (2006), arXiv:astro-ph / , 1233 (2009), arXiv:0907.1897 [astro-ph.EP].[71] D. N. C. Lin and J. Papaloizou, Astrophys. J. , 846 (1986). [72] G. Bryden, X. Chen, D. N. C. Lin, R. P. Nelson, and J. C. B.Papaloizou, Astrophys. J. , 344 (1999).[73] A. Crida, A. Morbidelli, and F. Masset, Icarus , 587 (2006),arXiv:astro-ph / , 758 (1995), arXiv:astro-ph / , 79 (2014), arXiv:1306.1871 [astro-ph.HE].[76] J. C. McKinney, A. Tchekhovskoy, A. Sadowski, andR. Narayan, MNRAS , 3177 (2014), arXiv:1312.6127[astro-ph.CO].[77] A. Gruzinov, Y. Levin, and C. D. Matzner, MNRAS , 2755(2020), arXiv:1906.01186 [astro-ph.HE].[78] X. Li, P. Chang, Y. Levin, C. D. Matzner, and P. J. Armitage,MNRAS , 2327 (2020), arXiv:1912.06864 [astro-ph.HE].[79] W. Lyra, S.-J. Paardekooper, and M.-M. Mac Low, Astro-phys.J.Lett. , L68 (2010), arXiv:1003.0925 [astro-ph.EP].[80] J. M. Bellovary, M.-M. Mac Low, B. McKernan, and K. E. S.Ford, Astrophys.J.Lett. , L17 (2016), arXiv:1511.00005[astro-ph.GA].[81] In fact, Dittmann and Miller [95] also noted that the migrationtraps in TQM disks no longer stand if a more updated opacityis used in solving the disk structure.[82] N. W. C. Leigh, A. M. Geller, B. McKernan, K. E. S. Ford,M. M. Mac Low, J. Bellovary, Z. Haiman, W. Lyra, J. Samsing,M. O’Dowd, B. Kocsis, and S. Endlich, MNRAS , 5672 (2018), arXiv:1711.10494 [astro-ph.GA].[83] B. McKernan, K. E. S. Ford, R. O’Shaugnessy, andD. Wysocki, MNRAS , 1203 (2020), arXiv:1907.04356[astro-ph.HE].[84] E. Y. Vilkoviskij and B. Czerny, A&A , 804 (2002),arXiv:astro-ph / , 240 (2016),arXiv:1604.05309 [astro-ph.GA].[86] T. Panamarev, B. Shukirgaliyev, Y. Meiron, P. Berczik, A. Just,R. Spurzem, C. Omarov, and E. Vilkoviskij, MNRAS ,4224 (2018), arXiv:1802.03027 [astro-ph.GA].[87] M. J. Rees, M. C. Begelman, R. D. Blandford, and E. S. Phin-ney, Nature (London) , 17 (1982).[88] A. P. Lightman and D. M. Eardley, Astrophys.J.Lett. , L1(1974).[89] T. Piran, Astrophys. J. , 652 (1978).[90] A. King and C. Nixon, MNRAS , L46 (2015),arXiv:1507.05960 [astro-ph.HE].[91] K. Schawinski, M. Koss, S. Berney, and L. F. Sartori, MNRAS , 2517 (2015), arXiv:1505.06733 [astro-ph.GA].[92] A. Soltan, MNRAS , 115 (1982).[93] B. F. Schutz, Nature , 310 (1986).[94] S. M. Mills, D. C. Fabrycky, C. Migaszewski, E. B. Ford, E. Pe-tigura, and H. Isaacson, Nature , 509 (2016).[95] A. J. Dittmann and M. C. Miller, MNRAS493