{'fraction_non_alphanumeric': 0.041081021414379276, 'fraction_numerical': 0.025259400986091692, 'mean_word_length': 4.159168564920273, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Resonant Dynamical Friction in Nuclear Star Clusters: Rapid Alignment of an Intermediate-mass Black Hole with a Stellar Disk

We investigate the dynamical evolution of an intermediate-mass black hole (IMBH) in a nuclear star cluster hosting a supermassive black hole (SMBH) and both a spherical and a flattened disk-like distribution of stellar-mass objects. We use a direct N-body (phiGPU) and an orbit-averaged (N-ring) numerical integrator to simulate the orbital evolution of stars and the IMBH. We find that the IMBH's orbit gradually aligns with the stellar disk if their mutual initial inclination is less than 90 degree. If it is larger than 90 degree, i.e. counterrotating, the IMBH does not align. Initially, the rate of orbital reorientation increases linearly with the ratio of the mass of the IMBH over the SMBH mass and it is orders of magnitude faster than ordinary (i.e. Chandrasekhar) dynamical friction, particularly for high SMBH masses. The semimajor axes of the IMBH and the stars are approximately conserved. This suggests that the alignment is predominantly driven by orbit-averaged gravitational torques of the stars, a process which may be called resonant dynamical friction. The stellar disk is warped by the IMBH, and ultimately increases its thickness. This process may offer a test for the viability of IMBH candidates in the Galactic Center. Resonant dynamical friction is not limited to IMBHs; any object much more massive than disk particles may ultimately align with the disk. This may have implications for the formation and evolution of black hole disks in dense stellar systems and gravitational wave source populations for LIGO, VIRGO, KAGRA, and LISA.

The IMBHs are expected to sink deep into the Galactic Center due to dynamical friction (Levin et al. 2005;Mastrobuono-Battisti et al. 2014). In the vicinity of the central SMBH, the IMBHs settle on short-period (1−10 4 yr) approximately Keplerian orbits. Then the eccentricity and the orbital inclination change predominantly due to the orbit-averaged gravitational interactions of the stars, a process known as resonant relaxation (Rauch & Tremaine 1996). Since the mean spherical potential of the nuclear star cluster also drives rapid apsidal in-plane precession with 10 4 -10 5 yr period, this limits the resonant accumulation of torques that drive the eccentricity evolution. However, one component of the resonant torques that drive the reorientation of the orbital planes is not limited by apsidal precision; this process is known as vector resonant relaxation (VRR).
The timescale of VRR has been determined for singlemass and two-component clusters of stellar objects to be between 10 5 − 10 7 yr (Hopman & Alexander 2006;Eilon et al. 2009;Kocsis & Tremaine 2011Giral Martínez et al. 2020). This timescale is sufficiently short that the orbital inclinations may be expected to reach a quasi-stationary equilibrium distribution while the eccentricities and semimajor axes change much more slowly (Roupas et al. 2017;Takács & Kocsis 2018;Fouvry et al. 2019;Roupas 2020). However, entropy is expected to be maximized when the inclination distribution exhibits mass segregation, meaning that higher-mass objects are ultimately confined to smaller root-mean-square inclinations compared to lower-mass objects (Rauch & Tremaine 1996;Roupas et al. 2017;Szölgyén & Kocsis 2018;Gruzinov et al. 2020). The eccentricity of higher-mass objects may also become systematically different Gruzinov et al. 2020). The systematic change of the eccentricity of massive objects has been confirmed by numerical simulations (Levin et al. 2005;Alexander et al. 2007;Löckmann et al. 2009;Iwasawa et al. 2011;Sesana et al. 2011;Madigan & Levin 2012;Foote et al. 2020;Bonetti et al. 2020). Mass segregation in the vertical direction has also been identified in numerical simulations of initially strongly non-axisymmetric nuclear stellar disks that lack a spherical cusp of stars (Foote et al. 2020).
In this paper, we examine the dynamical mechanism which leads to the rapid reorientation of the orbital plane of a massive object (e.g. an IMBH) in response to a population of lower-mass stars orbiting a SMBH in a disk configuration, a process that may be called resonant dynamical friction (RDF) (Rauch & Tremaine 1996). In particular, we investigate how an IMBH settles into the midplane of a stellar disk, such as the clockwise disk around Sgr A in the Galactic Center (Bartko et al. 2009;Lu et al. 2009;Yelda et al. 2014;Gillessen et al. 2017), and study the response of the disk using numerical simulations.
We use two different numerical methods: (i) ϕGPU -a direct N-body simulator accelerated by graphical processing units (Berczik et al. 2011(Berczik et al. , 2013, and (ii) Nring -a secular N-body simulator that uses the orbitaland precession-period averaged interactions (Kocsis & Tremaine 2015).Both ϕGPU and N-ring account for the superposition of the interactions between all pairs of particles, but while ϕGPU simulates the instantaneous interaction between the point masses N-ring calculates the interaction between annuli covered during the orbital and precession period. The semimajor axes and eccentricities of different stars are free to change in ϕGPU while they are fixed by construction in Nring. The difference between the two methods allows us to identify the main dynamical mechanism driving the evolution and in particular to explore the contribution of Chandrasekhar's (ordinary nonresonant) dynamical friction (CDF) and resonant dynamical friction (RDF) to the process of inclination relaxation. We find that the N-ring simulations, which account for RDF but not CDF by construction, match the results of the ϕGPU simulations that include both CDF and RDF. The alignment process is much more rapid in the simulations than expected by a simple analytic estimate of CDF. We examine the process using numerical simulations with different initial IMBH orbital inclinations, IMBH masses, number of stars, and radial surface density profiles for the stellar disk. We measure the warp and the thickness of the disk, as well as the evolution of the IMBH's semimajor axis, eccentricity, and inclination. We construct a simple empirical analytic model for the alignment via RDF.

NONRESONANT DYNAMICAL FRICTION
We start with a simple analytic estimate of the alignment time of the IMBH due to Chandrasekhar's dynamical friction (CDF) with respect to the disk stars during disk crossings. The mean deceleration of the IMBH is given by Chandrasekhar's formula for a homogeneous medium as (Eq. 8.1a in Binney & Tremaine 2008): where ρ is the density of stars and ln(Λ) is the where v rel is the relative velocity of stars and the IMBH, h d is the thickness of the disk, and f (v) is the probability density function of disk star velocities in a small box around the crossing point of the IMBH. 1 The alignment time due to dynamical friction accounting for two disk crossings of duration t cross = h d /v IMBH,z per orbit of period t orb may be estimated as where v IMBH,z is the IMBH's velocity in the z-direction, a CDF,z = dv IMBH,z /dt and ι is the inclination angle of the IMBH with respect to the stellar disk.
In the limit of small IMBH orbital eccentricity and a cold thin disk of total mass m d and maximum radius of r d with density and ρ(r, z) = 0 otherwise. We define Cartesian coordinates (x, y, z) where the disk plane is spanned by {x, y} and the IMBH crosses the z = 0 midplane at x = 0, substitute v IMBH = (v K cos ι 0 , 0, v K sin ι 0 ) and assume that all stars have v = (v K , 0, 0) in the vicinity of the crossing point of the IMBH in the disk in Eq. (1) where v K = (Gm SMBH /r) 1/2 is the Keplerian velocity. Then the alignment time simplifies to where we have introduced the "local disk mass at r" as (5) The CDF alignment time is 43 Myr for ι = 45 • at r = 0.15 pc for our fiducial model with m SMBH = 10 6 M , m IMBH = 10 3 M , m d = 8, 191M , m d,loc = 2, 273M , γ = 0.75, and r d = 0.5 pc. For our alternative model with γ = 1.75 and all other parameters fixed, we have m d,loc = 1, 516M , so t align,CDF = 65 Myr. Note that the result depends on the radius and the stellar disk properties only through the orbital period and the local disk mass. It is otherwise independent of h d and the mass of individual disk stars, but it shows a strong inclination dependence as shown in Figure 12 below.
We also find a similar result by measuring the CDF alignment time in the numerical simulations by substituting the actual velocity distribution of stars ( Figure  14, see Appendix 5) around the IMBH within a cylinder around the crossing point of the IMBH in the disk (r = 0.15 pc) of height h d = 0.01pc (i.e. the local thickness of the stellar disk) and radius chosen arbitrarily to be r cyl = 0.05pc in Eq. (1) and find a similar result.
We note, however, that t align,CDF is an upper limit for the true alignment time for, several reasons. First, dynamical friction is enhanced significantly in an inhomogeneous cuspy density profile around SMBH (Arca-Sedda & Capuzzo-Dolcetta 2014b). Further, dynamical friction may be accelerated by orbit-averaged torques, which drives resonant relaxation (Madigan & Levin 2012). We confirm this estimate in detail below. For the system parameters listed above, the apsidal precession rate of the IMBH is 0.019 times the orbital angular velocity around the SMBH according to Eq. (A1) in Kocsis & Tremaine (2015). Thus, the apsidal precession period of the IMBH is t a,prec = 0.28 Myr, which is much larger than the orbital period, t orb = 5.4 kyr. However, it is much shorter than the initial nodal precession period of the IMBH and the stellar disk measured using N-ring to be t n,prec = 2.9 Myr, which approximately characterizes the RDF timescale. 2 Eq. (4) shows that t align,CDF ∝ ι 4 for small ι, implying that CDF may play an important role in aligning the orbit during the later phase of the evolution when the inclination is small (see Section 4.8 for further discussion). 3

NUMERICAL METHODS
We run simulations using ϕGPU (Berczik et al. 2011(Berczik et al. , 2013Li et al. 2012Li et al. , 2017) a direct-summation N -body code that uses the Hermite integration scheme with individual block time steps to solve the instantaneous equations of motion of point particles with gravitational softening. We also run N-ring (Kocsis & Tremaine 2015) with the same initial conditions. In contrast with ϕGPU, N-ring is a secular code that integrates the pairwise interactions of orbit-averaged stellar trajectories using a time-reversible symplectic parallel scheme using a multipole expansion. By construction, N-ring neglects variations in the semimajor axes and eccentricities, implying that it neglects two-body relaxation and scalar resonant relaxation, and thus simulates the effects of VRR only, while ϕGPU does not make any such approximations. We direct the readers to the references for a detailed discussion of these numerical methods. Nring neglects the contribution of multipoles with > 50 which amounts to softening with an angular separation of 1.5 • in angular momentum direction space, while ϕGPU applies softening with 10 −5 pc for stars and 10 −3 pc for the SMBH and IMBH. We use both codes to understand the relative contribution of VRR and two-body relaxation to the inclination relaxation.
First, we run 6 ϕGPU and N-ring simulations with different initial inclinations for the IMBH's orbit as specified below. Then, we explore the evolution for different disk surface density profiles and IMBH masses using N-ring to determine the main properties of resonant dynamical friction and to derive an empirical analytic model.

Initial conditions
Our models consist of a disk of N = 8191 equalmass 1 M stars together with an m IMBH = 10 3 M IMBH and a m SMBH = 10 6 M SMBH fixed at the center. In addition, we include a static Plummer-potential Φ P (r) = −Gm P / r 2 + r 2 P with m P = 10 5 M and r P = 0.2 pc, in order to mimic the gravity of the spherical stellar component in the galactic nucleus, which drives apsidal precession and quenches the coherent torques of eccentric orbits. 4 While this model is obviously far from a realistic representation of the Galactic Center, which has a Bahcall-Wolf cusp of 10 7 stars with a possible inner cavity, a population of S-stars, a warped and twisted stellar disk of 10 4 M stars, 10 4 black holes, a massive molecular torus, etc. (Genzel et al. 2010), it serves as a simplified model to gain an understanding of the orbital alignment of a massive point mass in an anisotropic system.
The stars are initially distributed in the simulation in an axisymmetric thin disk with height h d = 0.01 pc, which is uniform in cylindrical coordinates along the azimuth angle and the z-axis. The surface density in our fiducial model is chosen according to Eq. (3) with surface density ∝ r −0.75 between 0.03 − 0.5 pc from the SMBH. This matches the expectation for the cluster infall scenario (Berukoff & Hansen 2006), but it is much shallower than the observed surface density profile of the clockwise disk in the Galactic Center, which is proportional to r −2.1±0.4 (Lu et al. 2009), r −1.95±0.25 (Bartko et al. 2009), or r −1.4±0.2 (Bartko et al. 2010). We also explore a steeper density profile of r −1.75 separately below. The initial velocities of stars are perturbed with respect to the circular velocities 5 corresponding to the SMBH and the Plummer- is a randomly sampled value from the PDF of a Maxwell-Boltzmann distribution that has the mean at v c,i , and x i , y i are the coordinates of the ith star. We choose the direction of v c,i such that 95% and 5% of stars are orbiting respectively in the same sense and in the opposite sense as the IMBH in the disk. In comparison, 19 out of 90 WR/O-stars in the Galactic Center are consistent with belonging to a retrograde disk (Bartko et al. 2009) although they are also consistent with belonging to a spherical distribution (Lu et al. 2009;Yelda et al. 2014). This disk configuration is not in statistical equilibrium initially, thus we first let the system to evolve for a few million years until the disk reaches a mean-field equilibrium in which the inclination distributions of the inner, overlapping, and outer stars have mean and standard deviation cos 2 ι = 0.994, 0.996, 0.998 and σ(cos 2 ι) = 0.022, 0.036, 0.012, respectively; see Sec. 3.2 below. We show the relaxed distribution of orbital parameters in the Appendix; see Figure 13.
We use this relaxed disk model together with an IMBH parameterized by (m IMBH , a IMBH , e IMBH ) = (10 3 M , 0.15pc, 0.33) as the initial condition in both the ϕGPU and in the N-ring simulations. We run 6 simulations with ι 0 = 35 • , 45 • , 55 • , 65 • , 75 • , 135 • different orbital inclinations for the IMBH. The initial orbital parameters and masses of stars, as well as the IMBH are identical in the ϕGPU and the N-ring simulations. N-ring does not include the Plummer-potential component directly, but its influence is included in the assumption of rapid apsidal precession, as the interaction in N-ring is averaged over the orbital and the apsidal precession period. 5 Note that v c,i is the Keplerian circular velocity for a point mass of m SMBH +m P , which is somewhat larger than the circular velocity of the actual SMBH+Plummer potential by a factor (m SMBH + m P )/{m SMBH + m P [1 + (r P /r) 2 ] −3/2 } ≈ 1 + (m P /m SMBH ){1 − [1 + (r P /r) 2 ] −3/2 }, which is between 1 and 1.1.

Characterization of disk thickness and warp
To characterize the state of the disk in the simulation, we introduce the quadrupole tensor of the stellar angular momenta: where L i denotes the angular momentum vector of the i th star, and α and β respectively label the Cartesian vector components. The largest eigenvalue (denoted by q hereafter) quantifies the flatness, as it is given by where ι i is the orbital inclination of stars relative to the disk midplane, i.e., the angle subtended by L i and the principal eigenvector of Q αβ . 6 We characterize the root-mean-square warp angle or thickness of the disk by the quantity Note that ∆ι = 0 for a razor-thin flat disk, and ∆ι is independent of an overall tilt of the disk but it increases with the warp of a thin disk and with the thickness. In the limit of an isotropic distribution, or a razor thin disk warped by 180 • , ∆ι = cos −1 (3 −1/2 ) = 54.7 • .
We measure ∆ι for the inner, overlapping, and outer regions of the disk with respect to the IMBH, separately, by evaluating the sums over only the corresponding stars in Eq. (6). We define these regions depending on the peri-and apocenter (r p , r a ) of the stars with respect to the IMBH (r p,IMBH , r a,IMBH ) as • inner region: r a < r p,IMBH , • overlapping region: r p r a,IMBH or r p,IMBH r a , • outer region: r a,IMBH < r p .
The number of stars in the three regions are N in = 9.4% ± 0.64%, N overl = 18.1% ± 4.3%, and N out = 72.5% ± 4.0%, respectively, for our standard model. Here, the fluctuations are mostly due to variations in the IMBH eccentricity with time.

RESULTS
We first present the results of ϕGPU simulations for initial IMBH inclination ι 0 = 45 • . Figure 1 shows the scatter plot of orbital inclinations and semimajor axes of disk stars and the IMBH, cos ι i = L i · L tot /(|L i ||L tot |), for 6 representative snapshots as indicated in the panels. Here L tot is the total angular momentum of the system. We also provide an animated video showing the complete structural evolution of the disk as the IMBH sinks down and aligns with the disk. The stars are depicted by red, blue, or green points depending on their peri-and apocenter (r p , r a ) with respect to the IMBH (r p,IMBH , r a,IMBH ) (see definition below Eq. 8): inner orbits in red; outer orbits in blue; and overlapping orbits in green. The black filled circle indicates the IMBH. The figure shows that initially the disk is flat and thin and the IMBH is at ι 0 = 45 • . The inclination angles start to change visibly at 0.4 Myr. At 0.9 Myr, the inner regions develop a larger warp. At 1.8 Myr, the IMBH aligns the disk first and the disk is significantly warped in the inner and overlapping region. At 2.7 and 4.5 Myr, the IMBH settles into the midplane, and the disk warp is limited mostly to the inner region to within r < 0.08 pc, where the local disk mass Eq. (5) is smaller than the IMBH mass, but the disk thickness in the overlapping and the close by parts of the outer disk is much larger than initially. We examine the evolution of the orbital elements of the IMBH and the response of the disk in detail next. Figure 2 shows the semimajor-axis and eccentricity evolution of the IMBH. The IMBH orbit circularizes during the relaxation process as also found by Madigan & Levin (2012) for a corotating disk and by Bonetti et al. (2020). The mean eccentricity decays from 0.33 to 0.02 during 4.5 Myr. The semimajor axis is approximately constant during the evolution to within 6% until the IMBH aligns with the disk at 2.6 Myr, when the semimajor axis decreases suddenly by another 6% within 0.1 Myr then decreases gradually by 9% until 4.5 Myr.

Semimajor-axis and eccentricity evolution
The semimajor-axis distribution of stars shows a depletion around the IMBH's semimajor axis in the bottom left panel of Figure 13 in the Appendix. The scalar angular momentum of the IMBH varies by less than 6%.

Orbital alignment and warps
The top panel in Figure 3 shows the time evolution of the angle between L IMBH and L tot (black line), where L IMBH is the angular momentum vector of the IMBH. Colored curves show the angle between the respective principal angular momentum eigenvector 7 of the disk stars (i.e. that of Q αβ , see Eq. 6) and L tot separately for the inner, overlapping, and outer stars as defined below Eq. bottom panel characterizes the evolution of the warp and/or the disk thickness within the three regions. The top panel shows that the IMBH sinks into the plane of the disk in an accelerated way, and settles after around 2.6 Myr. The disk is initially flat, but the inner disk quickly becomes tilted as the orientation of the innermost part of the disk separates from that of the overlapping and the outer disk in the top panel. At around > 0.6 Myr, the overlapping part of the disk also decouples from the outer disk. The inner disk tilt peaks at ∼ 30 • at ∼ 2.5 Myr, while the overlapping stars' inclination peaks around ∼ 12 • at almost the same time.

8
The bottom panel of Figure 3 shows that the curvature or thickness of the disk increase in all three regions until the IMBH settles into the disk at 2.6 Myr, and this is more prominent in the inner and overlapping regions. The thickness/curvature continues to increase only in the inner and overlapping regions after the IMBH aligns with the disk. By 4.5 Myr, ∆ι is increased by respective factors of 7.4, 5.1, and 4.1 in the inner, overlapping, and outer regions.

Dependence on the initial IMBH inclination
Let us now compare the results of models with different initial inclinations. Figure 4 shows the evolution of the IMBH inclination in the top panel with respect to the disk (ι = cos −1 (L disk ·L IMBH )) for 6 different initial r a < r p,IMBH r p < r a,IMBH or r p,IMBH < r a r a,IMBH < r p ∆ι r a < r p,IMBH r p < r a,IMBH or r p,IMBH < r a r a,IMBH < r p Figure 3. Top: The time evolution of the orbital inclination of the IMBH and stars for the inner, overlapping, and outer region with respect to the total angular momentum, Ltot. Black curve shows the alignment angle of the IMBH, defined as the angle between the angular momentumLIMBH andLtot. We characterize the orientation of the stellar disk with the respective principal angular momentum eigenvectors (Eq. 6) of the stars in the given radial group, ζgroup. Red, green, and blue curves show the angle betweenζgroup andLtot for the stars in the inner, radially overlapping, and outer parts of the disk with respect to the IMBH. Bottom: The time evolution of the disk warp/thickness in the inner (red), overlapping (green), and outer region (blue) defined by Eq. (8).
inclinations, 35 • , 45 • , 55 • , 65 • , 75 • , and 135 • , from bottom to top. The bottom panel shows the evolution of ∆ι (Eq. 8). Solid and dashed lines show simulations with ϕGPU and N-ring, respectively, started from the same relaxed disk initial conditions (see Section 3.1). For ∆ι this is shown in separate subpanels for clarity, as the curves are intersecting.
The figure shows that the ϕGPU and N-ring results are in qualitative agreement. The alignment is somewhat more rapid in ϕGPU by a factor ∼ 10% (ι 0 = 35 • ), ∼ 5% (ι 0 = 45 • ), and ∼ 4% (ι 0 = 55 • ) than in the N-ring results. Regardless of the value ι 0 , the discrepancy is more significant once the alignment angle is less than ι 30 • . The difference is well-explained by the differences in the approximations; N-ring neglects twobody relaxation and SRR by assuming that the eccentricity and semimajor axes of all objects are conserved. As shown above by Eq. (4), energy exchange by close two-body encounters results in Chandrasekhar dynamical friction, which leads to an alignment timescale of 43 Myr at ι = 45 • , 9.4 Myr at 30 • and 1.9 Myr at 20 • . The total time needed for alignment is determined by the evolution at large ι if ι 0 > 30 • , where Chandrasekhar dynamical friction is insignificant. Note further that Figure 2 has shown that the eccentricity of the IMBH changes significantly over a Myr, which is a manifestation of scalar resonant relaxation. While this is also neglected in N-ring, VRR is relatively insensitive to eccentricity for e < 0.7 as shown in Figure 7 of Kocsis & Tremaine (2015). Thus, we conclude that the alignment is predominantly driven by VRR if ι 0 > 30 • , and in this case the total alignment time is well-modeled by N-ring.
The counter-rotating simulations do not lead to systematic deviations between the two methods, showing that VRR is dominant over scalar resonant relaxation and two-body relaxation in this case. Here, the IMBH's inclination angle and the disk exhibit small-amplitude periodic oscillations. In this case, the conservation of total angular momentum vector inhibits the IMBH's orbital flip. For counter-rotating orbits, two-body interactions between disk stars and the IMBH is extremely weak to drive angular momentum exchange between the IMBH and the disk. We ran the N-ring simulations for an extended timescale of 150 Myr for ι 0 = 105 • and 135 • , and found small-amplitude oscillations throughout the evolution.
The similarity of the results of the two numerical methods at large inclinations ι 0 30 • in Figure 4, shows that the IMBH's alignment is mainly driven by VRR. Rauch & Tremaine (1996) named this process resonant dynamical friction (RDF). Next, we explore this process by running N-ring simulations on longer timescales, different IMBH masses, and surface density profiles.
Note that α = 1 matches the value for CDF. For larger IMBH masses, the mass dependence of alignment rate is shallower, which is to be expected since in this case the local disk mass is comparable to or smaller than the IMBH mass, and angular momentum conservation inhibits the rapid relaxation process.
We find that the IMBH ultimately aligns with the stellar disk for m IMBH ≥ 250M for our standard The α exponent of the mass dependence of the alignment rate, dι/dt ∝ m α , measured as a function of inclination in a disk with surface density exponent γ = 0.75 (top) and 1.75 (bottom). We compare simulations with IMBH mass m1 = 250 M versus m2 = 500 M (blue) or m2 = 1000M (red) and plot ln(ι 1 /ι 2 )/ ln(m1/m2) where ι 1 = dι1/dt. The value of ι1 is the inclination of IMBH with mass m1 and similarly for ι 2 . All simulations were initialized with ι1 = ι2 = 45 • . In all cases, the mass dependence is close to linear, α = 1, initially.

Dependence on the disk mass, radius, SMBH mass
Due to the scale-free nature of the N-ring simulations, the simulation time in code units is converted to physical time by scaling with where r unit and m unit are the distance and mass units adopted in the simulation. Here, t orb (r) = 2π(Gm SMBH ) −1/2 r 3/2 is the orbital time. Note that the N-ring simulation in code units is independent of m SMBH ; this quantity enters the results only when converting the code units to physical units using Eq. (10). Indeed, N-ring follows orbit-averaged interactions among stellar objects; the influence of the SMBH is accounted for by the orbit averaging. Eq. (10) shows that, if the mass of the IMBH and the disk are both increased by the same factor κ, keeping the SMBH mass and all other parameters fixed, this is equivalent to changing m unit to κm unit and the evolution is identical if the time is scaled by the factor κ −1 . However, we have shown in Section 4.4 that the alignment rate scales initially approximately with m α IMBH where α ≈ 1 initially, and later by a somewhat different value for a range of IMBH and disk masses with m IMBH ≤ m d,loc . Based on dimensional analysis, 9 combining with Eq. (10) implies that if the RDF rate of alignment scales with the masses as dι/dt ∝ m α IMBH m β d,loc m δ SMBH then β = 1 − α and δ = 1. For α ≈ 1 this is independent of m d,loc . Note, however, that here ι is the relative angle between the IMBH and the disk, and the rate of reorientation of the IMBH only (with respect to an inertial frame) may be different if the disk angular momentum is comparable to or less than the IMBH's (see Section 4.8).
Furthermore, Eq. (10) shows that, if the semimajor axes of the IMBH and all disk stars are scaled by a factor κ, this results in an identical evolution if the time unit is increased according to the change in the orbital time, i.e., by a factor of κ 3/2 .

Dependence on the number of stars
To examine how the results depend on the total number of disk stars N , we compare simulations with approximately fixed total disk mass and IMBH mass, , and find that the respective discrepancies in the alignment time are less than 10% and 16%. These simulations assume γ = 0.75 and ι 0 = 45 • . We conclude that resonant dynamical friction in many cases does not depend on the mass of individual disk stars and their total number for fixed stellar surface density profile.
However, we found that the number of particles N significantly affects the long-term evolution for surface density profile γ = 1.75, high initial IMBH inclination, and IMBH masses less than the local disk mass. The time duration of the orbital inclination hang-up phase is increased for higher N and the oscillation amplitude is decreased in this phase; see further discussion in Sec. 4.7. 9 Indeed, denoting the dimensionless code units with overline: m IMBH = m IMBH /m unit ,m d,loc = m d,loc /m unit ,t = t/t unit , dι/dt = t unit dι/dt and given that the evolution in code units is independent of m SMBH , assuming that dι/dt ∝ m α IMBH m β d,loc for α and β real numbers, implies that dι/dt ∝ t −1 unit (m IMBH /m unit ) α (m d,loc /m unit ) β . Substituting t unit from Eq. (10) shows that β = 1 − α to ensure that dι/dt is independent of the assumed units.

4.7.
Dependence on the disk surface density profile, γ = 1.75 Up to this section, we mostly used stellar disk surface density profile ρ ∝ r −γ with γ = 0.75 in Eq. (3). Now using N-ring, we present the time evolution for a stellar disk with a steeper radial density profile with γ = 1.75, consistent with that in the Galactic Center (Bartko et al. 2009;Lu et al. 2009;Yelda et al. 2014). We keep all other parameters of the surface density profile (Eq. 3) and the IMBH the same as before. In this case, the local disk mass is 1.5× smaller at the IMBH semimajor axis, and the mass distribution is approximately uniform on a log scale m d,loc ∝ r 0.25 , while for γ = 0.75 most of the mass was at the outside m d,loc ∝ r 1.25 .
We initialize the IMBH (m IMBH = 1000M ) with ι 0 = 35 • , 45 • , 55 • , 65 • and 75 • , respectively, and follow the evolution of the system. Figures 7 and 8 show snapshots of the evolution of orbital inclinations vs. semimajor axes for 45 • and 75 • , respectively. In contrast to the case with γ = 0.75 where the inner disk gets warped rapidly, Figure 1, the inner disk remains flat for γ = 1.75 and ι 0 ≥ 55 • and the outer disk develops the warp starting from the outer edge. The warp propagates inward in time, at a rate that is smaller for larger ι 0 . If ι 0 ≥ 55 • , the IMBH orbital inclination exhibits a hang-up until the disk warp reaches close to the IMBH apoapsis, after which the IMBH quickly plunges into the disk. For highly inclined initial conditions, the orbital hang-up may be quite prolonged; see the yellow and purple curves in the top panel of Figure 9. Figure 9 shows the inclination of the IMBH (top) and the evolution of the disk warp/thickness (bottom) as a function of time for the 5 different cases. In comparison to the γ = 0.75 case (see Figure 4 and 11), here the local disk mass is 1.5× smaller, and the alignment occurs more slowly by ∼ 1.2× (ι 0 = 35 • ) and ∼ 1.6× (ι 0 = 45 • ). For larger ι 0 the alignment time is much longer due to the orbital hang-up by ∼ 4.3× (ι 0 = 55 • ), ∼ 8.5× (ι 0 = 65 • ), and ∼ 12.4× (ι 0 = 75 • ), respectively. In this case, the IMBH does not follow the same exponential decay as the empirical fit for γ = 0.75, Eq. (12). For simulations with ι 0 = 55 • , 65 • , 75 • , the IMBH inclination first decreases to an intermediate value of ∼ 45 • ) where the IMBH oscillates with a decreasing amplitude until structural changes in the disk allow the IMBH's direct and rapid infall.
Note that Chandrasekhar's dynamical friction has been neglected here, but that timescale is even longer than the alignment time found here (see the solid line and × symbols in Figure 12). CDF may possibly affect the evolution for thereby possibly shortening the orbital hang-up phase. A detailed study of the orbital hang-up phase requires direct N-body simulations over extended time periods, which are beyond the scope of this paper. Figure 10 shows how the number of disk stars, N , affects the orbital inclination hang-up phase for model (m IMBH , ι 0 , γ) = (1000M , 75 • , 1.75) for N = 2047, 4095, and 8191 and fixed total disk mass m d = 8191M . We find that the oscillation amplitude is increased around the intermediate equilibrium inclination of 45 • and the total time duration of the hang-up is decreased for smaller N . It is remarkable that, for systems that exhibit the orbital inclination hang-up phenomenon, the evolution remains sensitive to N even in the N 1000 limit. Indeed, Figure 10 shows that, while the initial rate of reoriention is independent of N for a fixed initial ι, the orbital hang-up phase may be significantly extended in time for very large N . We find that this is due to a discontinuity that forms in the inclination distribution in the disk in a narrow range of semimajor axis in the outer region, i.e. near a = 0.38 pc at 21 Myr for N = 8191 in Figure 8  is manifested in only a small number of stars, the evolution remains sensitive to N even for N = 8191. The final IMBH realignment takes place after the discontinuity propagates to overlapping radii. The discontinuity moves more slowly for a smoother disk, i.e. for larger N .

Analytic model
Following the arguments in Sec. 4.5, in the limit of a scale-free disk, and with an IMBH much larger than the mass of individual objects in the disk, we expect a scaling as 10 where α ≈ 1 initially as measured above (and it may be somewhat larger or smaller at later times) and f may depend on the IMBH inclination with respect to the actual state of the disk, which in turn may be expected to depend mostly on the initial value ι 0 , the ratio of IMBH mass to the local disk mass and the IMBH mass to the individual object mass in the disk. Note that, Similar to the top panel in Figure 9: the orbital inclination of the IMBH for ι0 = 75 • as a function of time relative to a stellar disk containing different numbers of disk stars N = 2047 (blue), 4095 (red), and 8191 (green), respectively. Other parameters are the same as in Figure 9 including the mass of the IMBH (1000M ), as well as the total mass and surface density profile of the disk: 8191M and r −1.75 . as long as ι ≈ ι 0 , the disk has not had time to change significantly and f is expected to depend mostly only on ι 0 . We determine this function by fitting to simulations next. We derive this function for the model with γ = 0.75 since this model is free from discontinuities and does not exhibit the orbital hang-up phenomenon, but we find that the resulting formula is a reasonable approximation to describe the initial rate of reorientation for the γ = 1.75 model as well. Figure 11 shows the results of N-ring for the evolution similar to the top panel of Figure 3 for the same surface density profile Eq. (3) with γ = 0.75, but on longer timescales, and solid lines show empirical fits to the evolution. The inclination angle decays at an expo-nentially growing rate in the simulations for ι 0 15 • , which may be fit for 35 • ≤ ι 0 ≤ 75 • initial inclination by where the parameters depend on ι 0 as These parameters have the following physical interpretation: ι max is the rate of alignment when the IMBH inclination reaches the edge of the warped disk, which is found to be approximately independent of the initial condition ι 0 ; ι min is the inclination when the IMBH reaches the edge of the warped disk, which increases with ι 0 to roughly the saturation value of ∆ι (Figure 3);and λ is related to the change of the alignment rate (ι 0 − ι max )/(ι 0 − ι min ). Note that these parameters are not universal for RDF. They may depend on the IMBH mass relative to the local disk mass, the radial dependence of the local disk mass, and the eccentricities. The particular scaling with the IMBH and SMBH mass and the Keplerian orbital time is motivated by Eq. (11). Based on Eq. (12), the rate of reorientation for RDF increases exponentially as which may be expressed in terms of the instantaneous inclination of the IMBH in the form of Eq. (11) as This model may be interpreted physically as follows. Initially, the alignment rate is very slow, as the stellar distribution and potential are nearly axisymmetric for which L z is conserved. This is analogous to the lack of dynamical friction in a homogeneous medium before any particle made a close encounter and before a density wake develops. The gravitational influence of the IMBH warps the disk, which in turn torques the IMBH's orbit to align toward the disk. The speed of reorientation grows exponentially in the simulations, which is typical for other types of instabilities in galactic dynamics (see, e.g. Sellwood 2013; Sellwood & Gerhard 2020). The instability here is driven by the IMBH, as the stellar disk was started from a stable configuration before deploying the IMBH. For larger ι 0 the disk has more time to respond and becomes more warped by the time the IMBH aligns with the disk. The rate of alignment is largest when the orbit first reaches the edge of the warped disk; here, ι = ι min and the first term vanishes in Eq. (17). Since this configuration is roughly independent of the initial ι 0 , in terms of the inclination between the IMBH and the local disk, the final rate of alignment is roughly universal for fixed (m d,loc , m IMBH , t orb , γ). However, the bottleneck for alignment is in the initial phase, which is highly sensitive to ι 0 .
Up to this point, we have discussed the relative angle between the IMBH and the disk ι = cos −1 (L IMBH · L disk ). Let us now determine the rate of reorientation of the IMBH only, i.e., relative to an inertial frame. Since L tot = L IMBH + L disk is conserved, we expect that (dι IMBH /dt)/(dι/dt) = [L disk /(L disk + L IMBH )] = m d,char /(m d,char + m IMBH ), where m d,char = m IMBH L disk /L IMBH is the "characteristic mass of the disk" at the IMBH. Substituting in Eq. (17) gives where Here, we substituted Eq. (3) for ρ with γ < 2 in the limit of a thin flat Keplerian disk and an IMBH on a circular orbit at r. Note that m d,char > m d as long as r ≤ r d . Thus, if the IMBH is much less massive than the disk, m IMBH m d , then µ IMBH ≈ m IMBH , as expected. Note that RDF is quite different from Chandrasekhar dynamical friction which is described by Eq. (4) as  (21), it is clear that RDF is more efficient than CDF by a factor of (m SMBH /m d,loc )/[1 + (m IMBH /m d,char )], i.e. approximately m SMBH /m d,loc for m IMBH m d < m d,char , and a factor that depends on inclination. For our standard disk model, this mass factor is 440. Thus, in comparison to CDF, the rate of alignment due to RDF is boosted most significantly for low-mass IMBHs (including heavy stellar-mass black holes), low-mass disks, and high-mass SMBHs. Note that for (m IMBH , m SMBH ) = (10 3 , 10 9 ) M and with all other parameters unchanged, the rate of reorientation due to Chandrasekhar's dynamical friction and resonant dynamical friction are respective factors of 10 6 and 10 3 slower. Thus, in this case t CDF /t RDF > 10 4 for ι 0 ∼ 45 • . Figure 12 compares the alignment time of an IMBH in N-ring for a disk of mass m d = 8, 191M and that due to Chandrasekhar's dynamical friction (Eq. 4) as a function of ι 0 . Dashed and solid black curves show the CDF (Eq. 4) for surface density profile exponent γ = 0.75 and 1.75, respectively, while the dotted black curve represents the empirical fit to RDF for γ = 0.75 (Eq. 17). Colored circles and crosses show the measured alignment times (up to the point when the relative inclination of the IMBH and the disk reaches ι = 15 • ) in the simulations with γ = 0.75 and 1.75, respectively. These correspond to values where the IMBH crosses 15 • in Figure 11 and 9 (top panel). The figure shows that RDF becomes dominant over CDF in aligning the IMBH beyond a critical inclination of ι 0 15 • .
The analytic model (Eq. 18) is successful at describing the initial rate of reorientation for both models with γ = 0.75 and 1.75, respectively. However, it fails to describe the total alignment time for systems that exhibit the inclination hang-up phenomenon (e.g. ι 0 ≥ 55 • for (N, m IMBH , γ) = (8191, 10 3 M , 1.75)). By comparing the disk structure for the different simulations shown in Figure 10, we find that systems that exhibit the orbital hang-up phase form a discontinuity in the inclination distribution in a narrow range of radii in the outer disk; see upper right panel of Figure 8 showing at 21 Myr in the region near a = 0.38 pc. The analytic model presented here describes the initial rate of reorientation before discontinuities may form.

DISCUSSION
We have examined the dynamical evolution of a nuclear star cluster around an SMBH with both a spherical and a stellar disk component, as well as an IMBH on an initially inclined orbit with respect to the disk. We ran numerical simulations with two codes, ϕGPU and Nring with the same initial conditions.
For an initial inclination of ι 0 = 45 • and for (m IMBH , m SMBH = (10 3 M , 10 6 M ), we found the eccentricity of the IMBH decreased from the initial value of 0.33 to 0.02 in 4.5 Myr in ϕGPU, while the eccentricity is fixed by construction in N-ring. This timescale for the IMBH's eccentricity decrease is much shorter than Chandrasekhar's dynamical friction timescale (Section 2), supporting previous findings by Madigan & Levin (2012) that resonant dynamical friction, which arises due to orbit-averaged torques, dominates over ordinary dynamical friction driven by hyperbolic encounters and decreases the eccentricity for a disk corotating with the IMBH. This is related to the process of scalar resonant relaxation.
To investigate the contribution of VRR to how the orbital plane of the IMBH aligns with the stellar disk, we compared ϕGPU direct N-body simulations with the orbit-averaged N-ring simulations that follow VRR and neglect two-body relaxation and scalar resonant relaxation. The results of the two types of methods were in approximate agreement with respect to the evolution of angular momentum vector directions (Figure 4). Since VRR is weakly sensitive to the orbital eccentricity for e < 0.7 if the spherical component of the cluster's potential drives rapid in-plane apsidal precession (Kocsis & Tremaine 2015), the change of eccentricity does not influence the VRR evolution strongly. We have confirmed with ϕGPU that the semimajor axis changes much more slowly during the alignment. The agreement between ϕGPU and N-ring for the evolution of the orbital inclination of the IMBH and the disk suggests that the alignment is mainly driven by resonant gravitational torques between the stellar disk and the IMBH.
We investigated the response of the stellar disk due to the interaction with an IMBH. In the case where disk surface density follows r −0.75 , we found that the inner region of the disk exhibits a strong warp, while the outer disk is initially not affected (Figure 1). Then, as the IMBH aligns with the midplane of the disk, the thickness of the disk increases in the radially overlapping regions and in the inner region. However, if the disk surface density is more cuspy and follows r −1.75 consistent with the clockwise disk of massive stars in the Galactic Center (Lu et al. 2009;Bartko et al. 2009Bartko et al. , 2010Yelda et al. 2014), the evolution is qualitatively different, since in this case the mass is nearly uniformly distributed on a logarithmic scale, and there is much more mass in the inner disk and less mass in the outer disk in comparison to the shallower surface density profile. The disk is then stiffer in the inner region and more prone to deformations in the outer regions (Figure 7 and 8). In case of high initial IMBH inclination (ι ≥ 45 • ) and if the fraction of IMBH mass to local disk mass is smaller than 0.7, then the IMBH alignment process is delayed and it takes place in three steps ( (ii) The orbital inclination settles and oscillates with a decreasing amplitude for an extended period of time In this phase, while the inner disk perturbations are small, the outer disk ultimately responds in an irreversible way by developing a strong discontinuous warp starting near the outer edge of the disk, which propagates inward.
(iii) Once the warp approaches the apoapsis of the IMBH orbit, the IMBH plunges rapidly into the disk. Ultimately, the inner disk is flat and aligned with the IMBH, and it becomes more warped and thicker in the outer regions.
We determined the alignment time of the IMBH with ϕGPU and N-ring for 6 different initial inclinations (ι 0 ) and found an analytic fit to the results. The evolution of the disk and the IMBH are in approximate agreement for ϕGPU and N-ring if ι 0 > 20 • . The alignment time is an exponential function of ι 0 if ι 0 < 90 • (Eq. 18 and Figure 12) and there is no alignment in the counter-rotating case, i.e., if ι 0 > 90 • . For the surface density profile of r −0.75 , the alignment of the IMBH with the stellar disk is between 9-28× faster for ι 0 = 35 • -75 • than the estimate for the Chandrasekhar type dynamical friction (Eq. 4). The orbital inclination alignment for ι > 20 • for component masses of (m SMBH , m disk , m IMBH ) = (10 6 M , 8 × 10 3 M , 10 3 M ) increases at an exponentially accelerating rate (Eq. 12). These findings are valid for the total alignment time for surface density r −0.75 for a wide range of masses and inclinations, and also for r −1.75 surface density as long as the initial inclination satisfies ι 0 45 • and the IMBH mass is not much smaller than the local disk mass. However, for surface density r −1.75 with large ι 0 and/or small IMBH mass, the exponential fit to the realignment (Eq. 12) is valid approximately only for phase (i), but the total alignment is delayed due to the orbital inclination hang-up phase (ii). Nevertheless, this is typically still much faster than the Chandrasekhar's dynamical friction ( Figure 12). The time duration of the inclination hang-up phase increases with the smoothness of the disk by increasing the number of disk particles for fixed total disk mass ( Figure 10). We also examined the dependence of the alignment rate on the mass of the IMBH and found that it is initially linearly proportional to mass, similarly to Chandrasekhar's dynamical friction, if the IMBH mass is much larger than the mass of disk particles (Figure 6).
At later times, as the disk gets significantly perturbed, the scaling with m IMBH can be either stronger or weaker, depending on the surface density exponent γ and m IMBH . Furthermore, we found that the resonant dynamical friction timescale is approximately proportional to m SMBH /m IMBH unlike Chandrasekhar's dynamical friction, which scales with m 2 SMBH /(m d m IMBH ). This implies that, while RDF is faster than CDF by a factor of 10-40 for our fiducial model with (m IMBH , m SMBH ) = (10 3 M , 10 6 M ), RDF is faster than CDF by a factor of more than ∼ 10 4 for (m IMBH , m SMBH ) = (10 3 M , 10 9 M ).
These findings may provide an explanation for the dynamical origin of the warped structures of the stellar disks in the Galactic Center (Bartko et al. 2009).
Our results indicate that an inclined massive perturber, like an IMBH, can strongly warp an initially flat stellar distribution, due to coherent gravitational torques. Depending on the mass of the IMBH and the surface density profile of the disk, it may cause a strongly perturbed inner or outer disk and a significantly thickened overlapping disk on a 1 − 100 Myr timescale after its formation or its arrival to the Galactic Center, see e.g., Figure 1.
Resonant dynamical friction may have applications beyond the relaxation of IMBHs examined in this paper. It may affect all objects in stellar clusters much more massive than the individual constituents of the disk, if present, including massive stars, stellar mass black holes (BHs), or the center of mass of massive bi-  Figure 13. Orbital parameter distributions of stars in the disk: the number density of semimajor axes and eccentricities, the distribution of inclinations as a function semimajor axis, and the surface density of the disk as a function of radius, respectively, from left to the right. First row shows the initial conditions at 0 Myr of the ϕGPU simulation, second row shows the parameters at 2.6 Myr, and the third row shows parameters at 4.5 Myr.
naries. Furthermore, it is also expected to operate in any type of disk with a high number of particles, including active galactic nucleus (AGN) accretion disks. Previously, it has been argued that stars and BHs crossing the disk on low-inclination orbits get captured by Chandrasekhar dynamical friction into the disk (Bartos et al. 2017;Panamarev et al. 2018;Tagawa et al. 2020). An interesting implication is that, if BHs settle into the disk, they interact dynamically and form BH-BH binaries efficiently, and frequent dynamical interactions and gas effects drive the BHs to merger, producing gravitational waves (GWs) detectable by LIGO, VIRGO, andKAGRA (McKernan et al. 2014, 2018;Bartos et al. 2017;Leigh et al. 2018;Yang et al. 2019;Tagawa et al. 2020Tagawa et al. , 2021Samsing et al. 2020). Mergers are also facilitated by Lidov-Kozai oscillations in anistropic systems (Heisler & Tremaine 1986;Petrovich & Antonini 2017;Hamilton & Rafikov 2019). The results in this paper show that resonant dynamical friction may accelerate the capture of objects in the accretion disks by a factor proportional to the SMBH mass over the local disk mass for large orbital inclinations. Pressure and viscosity in a gaseous disk do not inhibit the orbit-averaged torque from the IMBH, which leads to realignment and the warping of the disk (Bregman & Alexander 2012). Thus, RDF may efficiently catalyze the alignment of the orbital planes of BHs even in low-luminosity AGN or Seyfert galaxies with relatively small disk masses, which may not be possible for Chandrasekhar dynamical friction. In fact, this mechanism extends the scope of the "AGN merger channel" for GW source populations even beyond low-luminosity AGN and Seyfert galaxies, as it may organize BHs into disks also in nonactive galaxies with nuclear stellar disks.
Even more generally, resonant dynamical friction may be expected to operate in all systems whose mean-field potential admits action-angle variables where one of the precession frequencies is zero, e.g. for approximately spherical systems. RDF may accelerate the relaxation of massive subsystems including dark matter substructures in the spherical halo into galactic disks, which could then catalyze the formation of galactic spiral arms (D'Onghia et al. 2013), similarly to how an IMBH may excite spiral waves in a stellar disk in the Galactic center (Perets et al. 2018). It may also play a role in the formation and evolution of anisotropic groups of galaxies, including the disk of Milky Way satellites (Kroupa et al. 2005) or anisotropic galaxy clusters (Binney 1977).

ACKNOWLEDGMENTS
We gratefully thank Peter Berczik for providing us with the ϕGPU code and for his help in its use. We thank Scott Tremaine, Yohai Meiron, Hiromichi Tagawa, Nathan Leigh, Chris Hamilton, and John Magorrian for useful discussions. We are grateful to Mária Kolozsvári for help with logistics and administration related to the research. This work received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program under grant agreement No 638435 (GalNUC) and was supported by the Hungarian National Research, Development, and Innovation Office grant NKFIH KH-125675.Ákos Szölgyén was supported by theÚNKP-20-3 New National Excellence Program of the Ministry for Innovation and Technology from the source of the National Research, Development and Innovation Fund.

Distribution of orbital parameters
Here, we give the semimajor-axis, eccentricity, inclination, and mass distributions of this paper's fiducial relaxed stellar disk model that we introduced in Section 3.1. Figure 13 shows the number density of semimajor axes and eccentricities, the distribution of the cosine of the inclinations as a function semimajor axes, and the surface density as a function of radius (∝ r −0.75 ), respectively, in panels from the left to the right. The first row of panels represents the initial, relaxed disk model at 0 Myr when the IMBH is added to the simulation. The second and the third row show snapshots of ϕGPU simulation at 2.6 Myr (when the IMBH reaches the disk), and 4.5 Myr (end of the simu-lation). Here, the initial IMBH parameters are (m IMBH , a IMBH , e IMBH , ι 0 ) = (10 3 M , 0.15pc, 0.33, 45 • ).
The first column of panels show the evolution of disk stars' semimajor axes. As the IMBH's orbital inclination enters the disk at 2.6 Myr, the number density of stars around the IMBH's semimajor axis (0.15 pc) is strongly depleted. By 4.5 Myr, the number density decreases by a factor of ∼ 50% which suggests that the IMBH strongly scatters the disk stars out of their orbit, i.e., nonresonant dynamical friction becomes dominant over VRR.
The eccentricity density of stars does not vary much during the evolution, as shown by the second column of panels. However, the eccentricity distribution flattens as the peak decreases by 22% and the standard deviation increases by 50%. Note that, unlike Gualandris et al. (2012), we do not find the eccentric disk instability that forms a bimodal eccentricity distribution, but this may be due to the differences in our assumptions. The disk is initially axisymmetric in our simulations, and we account for the spherical potential of the stellar cusp, which drives differential apsidal precession and quenches scalar resonant relaxation.
The inclination of stars as a function of the semimajor axis shows a warp and an increase of the thickness of the disk in the third column of panels, which is also presented in Figure 1 and 3 in detail. Figure 14 shows the local 3D velocity density of stars in cylindrical coordinates at 0 Myr. The velocity distribution f (v r , v φ , v z ) is restricted to the stars in the neighborhood of the crossing point of the IMBH in the disk in a cylindrical box with radius r cyl = 0.05 pc and height h d = 0.01 pc.