Resonant Dynamical Friction Around a Super-Massive Black Hole: Analytical Description

Dynamical friction is a process whereby collective gravitation interactions between a massive body and a large number of smaller ones (whose total mass dominates) lead to a net effect on the motion of the former, in a manner mimicking friction (Chandrasekhar, 1943; Binney & Tremaine, 2008). The usual approach (e.g. Lynden-Bell & Kalnajs 1972; Tremaine & Weinberg 1984) to modelling dynamical friction treats the potential of the large mass as a perturbation to the whole system’s Hamilton’s equations of motion, which are solved perturbatively. At first order, the perturber scatters smaller objects from one orbit (in the background potential) to another, and so creates an over-density behind it; the force exerted on the perturber is then a consequence of the net gravitational pull of this wake.

Dynamical friction may also arise as a consequence of orbit-averaged interactions: suppose that the entire system consists of particles orbiting in the potential $\varphi_{0}$ of some object, which is so massive that it is entirely unaffected by them. Every particle primarily follows the orbits of this potential, but secular effects may appear due to their interaction with other particles (Arnold et al., 2006; Morbidelli, 2002), which may be accounted for by averaging the correction to the potential, $\varphi_{1}$, due to the self-gravity of the particles, over the orbits of $\varphi_{0}$. If one of the particles is much heavier than the others, these orbit-averaged equations describe a dynamical-friction-like force (e.g. Tremaine & Weinberg 1984). Indeed, such a phenomenon was recently observed numerically in the context of an intermediate-mass black hole (IMBH) entering a disc of stars around a super-massive black hole (SMBH) by Szölgyén et al. (2021), which was termed ‘resonant dynamical friction’. This process pertains to the alignment of the inclination of the IMBH’s orbit with the disc, and the main purpose of this paper is to offer an analytical explanation of this phenomenon. A possibly-related process occurs in rotating clusters (Szölgyén & Kocsis, 2018; Gruzinov et al., 2020).

In the inner-most regions of the Galactic centre, there is a disc of young massive stars (Levin & Beloborodov, 2003; Paumard et al., 2006; Gallego-Cano et al., 2018; Ali et al., 2020; Schödel, R. et al., 2020; von Fellenberg et al., 2022) orbiting around the SMBH Sgr A*, whose mass is $M_{\bullet}\approx 4\times 10^{6}~{}M_{\odot}$ (Schödel et al., 2002; Eisenhauer et al., 2005; Ghez et al., 2005; Ghez et al., 2008; GRAVITY Collaboration et al., 2018; Event Horizon Telescope Collaboration et al., 2022). The observed distribution of the arguments of the ascending node of these stars exhibits unusual features (Gallego-Cano et al., 2018; Ali et al., 2020; von Fellenberg et al., 2022), which has been interpreted as evidence for two distinct discs (e.g. Ali et al. 2020), but the existence of a second separate disc may be moot (von Fellenberg et al., 2022). Many phenomena can lead to features in the distribution of the arguments of the ascending node, $\Omega$, for example in the context of the Solar system: evidence for the hypothesised existence of Planet 9 is related, in part, to a non-uniform distribution of $\Omega$ (see, e.g. Batygin & Brown 2016). In this paper, we show that a massive perturber aligns the stars’ arguments of ascending node with a $90^{\circ}$ offset with respect to its own, and this configuration in turn causes resonant dynamical friction driving the alignment of the perturber’s orbital inclination.

Numerical simulations by Szölgyén et al. (2021) showed that an IMBH of mass $m_{\rm p}=\textrm{a few}\times 1000~{}M_{\odot}$ on an inclined orbit with respect to a stellar disc around an SMBH, warps the disc and the relative inclination decreases until the IMBH becomes embedded in the disc. The time-scale for this process was found to be similar to the time-scale of resonant relaxation (Rauch & Tremaine, 1996), but much shorter than the time-scale of inclination change due to Chandrasekhar dynamical friction.¹¹1Generally, ‘resonance’ refers to the relation between the fundamental frequencies of motion of different bodies in the system (Rauch & Tremaine, 1996). If the mean-field potential admits action-angle variables, resonance occurs between all bodies for which the angles change respectively with frequencies $(\Omega_{1},\Omega_{2},\Omega_{3})$ such that they satisfy a resonance condition: $n_{1}\Omega_{1}+n_{2}\Omega_{2}+n_{3}\Omega_{3}=0$, with integer $(n_{1},n_{2},n_{3})$ prefactors, not all zero. In particular, for vector resonant relaxation, the mean field potential is generated by the SMBH and the spherical extended stellar mass distribution and/or general relatistic corrections, which drive elliptic motion and apsidal precession, respectively, while the argument of node is fixed, $\Omega_{3}=0$. The resonance condition holds with $n_{1}=n_{2}=0$ and $n_{3}$ arbitrary for all bodies in the system, leading to a rapid coherent change in the $z$-component of the angular momentum vectors (Rauch & Tremaine, 1996). Physically, the orbits in the spherical potential cover punctured discs which interact coherently until the angular momentum vectors reorient. The coherent accumulation of torques due to this global resonance drives an accelerated evolution relative to the energy diffusion driven by stochastic incoherent two-body encounters. This study also compared an $N$-ring code and $N$-body simulations, and found a very good agreement between the two approaches; thus, the underlying physical mechanism is probably the same as what gives rise to resonant relaxation.

While it is unknown whether IMBHs exist in the relevant mass range in the Universe, if they do, the Galactic centre is a relatively likely place to find them (Yu & Tremaine, 2003; Goodman & Tan, 2004; Gualandris & Merritt, 2009; Gualandris et al., 2010; Kocsis et al., 2011; Kocsis et al., 2012; McKernan et al., 2014; Arca-Sedda & Capuzzo-Dolcetta, 2019; Naoz et al., 2020; GRAVITY Collaboration et al., 2020).

This paper explores the orbit-averaged gravitational interactions between $N$ stars lying initially in a flat disc and a perturber, which may be thought of as an IMBH, on an inclined orbit. We start by describing the problem and the orbit-averaged equations of motion in §2, we then solve these equations in §3 – which are singular, when studied in the context of perturbation theory – and derive solutions for the arguments of the ascending node of the disc stars and the perturber. We find that they align to have a phase-difference of $90^{\circ}$, and then show how resonant dynamical friction arises from these equations. Most of the calculations are done in appendices, and §3 summarises their conclusions. Then, we perform a set of $N$-body simulations to test our model, which are described in §4, and whose results are presented in §5. We discuss our results in §6 and summarise in §7.

Suppose that one has a thin disc of $N\gg 1$ stars of masses $m_{n}$ ($n\in\{1,\ldots,N\}$) of total mass $M_{\rm d}\equiv\sum_{n=1}^{N}m_{n}$, surrounding a SMBH whose mass is $M_{\bullet}$, and a particle of mass $m_{\rm p}$, the “perturber”, is placed on an inclined orbit. We assume that the following mass-hierarchy holds

and that the stars have small eccentricities and mutual inclinations, so the Hamiltonian is expanded to second order in them, but not in the eccentricity and the inclination of the perturber, which can be large. In the case we study here, the perturber moves under the influence of the entire disc, while each of the disc stars is affected primarily by the perturber, but also by the other stars (treated as a perturbation to the former). This is different from the case of pairwise interactions (cf. Kocsis & Tremaine 2015), because the perturber is affected by the collective gravity of all disc stars.

Moreover, since vector resonant relaxation (VRR) is the dominant perturbation with respect to the Keplerian orbit around the central SMBH (Kocsis & Tremaine, 2011, 2015; Fouvry et al., 2019), the Hamiltonian (i.e. the disturbing function) is orbit-averaged both over the mean anomalies and over the arguments of periapsis. This is justified by the fact that the semi-major axes and eccentricities are approximately conserved as two-body relaxation (and consequently Chandrasekhar dynamical friction) and scalar resonant relaxation take place on much longer time-scales, which may therefore be safely ignored. Hence the leading order Keplerian term in the Hamiltonian $-\sum_{i=1}^{n}Gm_{i}/(2a_{i})$ is a constant which may be omitted when solving for the evolution of the inclinations and the arguments of ascending node, which specify the directions of the angular momentum vectors.

The Hamiltonian governing the system is

where $\mathscr{H}_{\rm LL}$ is the Laplace-Lagrange Hamiltonian given by equation (32) of Kocsis & Tremaine (2011) which describes the orbit-averaged gravitational interaction between a thin disc of stars on nearly circular and nearly co-planar Keplerian orbits and which reduces to a system of harmonic oscillators, as given explicitly in equation (26). This approximation is justified by our assumption that the inclinations and eccentricities of the disc stars are small, and corrections to $\mathscr{H}_{\rm LL}$ are third order in the inclinations or eccentricities. $\mathscr{H}_{\rm p}$ describes the orbit-averaged interaction between the stars and the perturber (Kocsis & Tremaine, 2015; Roupas, 2020)

where $a_{i}$ are the semi-major axes, $\alpha_{pn}=\frac{\min\left\{a_{\rm p},a_{n}\right\}}{\max\left\{a_{\rm p},a_{n}\right\}}$, $P_{\ell}$ is the $\ell$-th Legendre polynomial, the dimensionless coefficient $s_{pnl}$ is a function of $\alpha_{pn}$,²²2Explicitly, $s_{pn2}$ is $$s_{pn2}\equiv\iint_{0}^{\pi}\frac{\mathrm{d}x\mathrm{d}y}{\pi^{2}}~{}\frac{\left[\min\left\{1+e_{\rm in}\cos x,\frac{1+e_{\rm out}\cos y}{\alpha_{pn}}\right\}\right]^{3}}{\left[\max\left\{\alpha_{pn}(1+e_{\rm out}\cos x,1+e_{\rm in}\cos y\right\}\right]^{2}},$$ where $e_{\rm in}$ is the eccentricity of the particle with the lesser semi-major axis, and similarly for $e_{\rm out}$. $e_{\rm p}$, $e_{n}$, and the multipole index $\ell$ only, where $s_{pn\ell}=1$ for circular orbits, and $\theta_{pn}$ is the angle between the angular momentum vector of the perturber and that of star $n$, viz.

Since $\mathscr{H}$ is already doubly-averaged over two of the three angle variables, the semi-major axes and eccentricities of all bodies are constant, and the only dynamical variables are their inclinations $\left\{i_{n}\right\}$ and arguments of ascending node $\left\{\Omega_{n}\right\}$. Up to constant coefficients, $\cos i_{n}$ and $\Omega_{n}$ are canonical momenta and position variables, respectively.

In this paper we perform a perturbative expansion in the parameter $\varepsilon\equiv m_{\rm p}/M_{\rm d}\ll 1$ keeping $M_{\rm d}$ fixed, but the main result of this paper concerning the arguments of the ascending node will be valid even for $\varepsilon=1$, albeit for short times. For two functions $f(\varepsilon)$ and $g(\varepsilon)$, we write $f(\varepsilon)=O(g(\varepsilon))$ if $\left|\frac{f(\varepsilon)}{g(\varepsilon)}\right|$ becomes bounded as $\varepsilon\to 0$, and $f(\varepsilon)=\textrm{ord}\left(g(\varepsilon)\right)$ if $\displaystyle\lim_{\varepsilon\to 0}\left|\frac{f(\varepsilon)}{g(\varepsilon)}\right|$ is non-zero.

The expansion, however, does not simply amount to the naïve one of setting $\mathscr{H}_{\rm p}=\textrm{ord}\left(\varepsilon\right)$ and $\mathscr{H}_{\rm LL}=\textrm{ord}\left(1\right)$, because, as we shall see, the case $\varepsilon=0$ is singular, and therefore it calls for a special type of expansion. In fact, we will show that $\mathscr{H}_{\rm p}$ dominates the dynamics of the disc particles for most of the relevant time; in the notation of the introduction, one would have $\varphi_{0}\equiv\mathscr{H}_{\rm LL}$ and $\varphi_{1}\equiv\mathscr{H}_{\rm p}$. This sets the phenomenon of resonant dynamical friction apart from other dynamical friction phenomena (e.g. Tremaine & Weinberg 1984; Ostriker 1999; Binney & Tremaine 2008; Magorrian 2021; Banik & van den Bosch 2021; Desjacques et al. 2022; Dootson & Magorrian 2022), where usually $\varphi_{1}\ll\varphi_{0}$ is treated as a perturbation, which acts only to scatter medium particles from one orbit in $\varphi_{0}$ to another, or trap them around resonant orbits in $\varphi_{0}$ (where dynamical friction arises from the back-reaction of this scattering). Here, quite the opposite occurs, and, in fact, if one did attempt to perform a naïve expansion, one would have found zero dynamical friction force acting on the perturber, to $\textrm{ord}\left(\varepsilon^{2}\right)$.

One can use Lagrange’s equations of motion, which yield (Murray & Dermott, 2000)

where the coefficient $\gamma_{\rm p}$ is defined as

where $\mu_{\rm p}\equiv m_{\rm p}M_{\bullet}/(m_{\rm p}+M_{\bullet})\approx m_{\rm p}$, and the frequency $n_{\rm p}$ is $\sqrt{G(M_{\bullet}+m_{\rm p})/a_{\rm p}^{3}}$. Likewise, one may define $\gamma_{n}$ for any of the disc particles, by replacing the index p by $n$.

Thanks to the mass hierarchy, Eq. (1), one might expect $\mathscr{H}_{\rm p}$ to be truly a perturbation relative to the disc’s self-interaction Hamiltonian $\mathscr{H}_{\rm LL}$, governing particles $1,\ldots,N$; this is accurate, however, only if the disc is not thin. Here, though, the thinness of the disc implies that it is in fact $\mathscr{H}_{\rm p}$ which dominates both the dynamics of the perturber and the disc, because $\mathscr{H}_{\rm LL}\propto~{}\left[\textrm{disc thickness}\right]^{2}$ for a thin disc. The explicit equations of motion are given in appendix A, and are then solved perturbatively there. Additionally, we provide a test case for our results in appendix B, where $i_{\rm p}$ is also taken to be small. There the equations of motion may be solved analytically, and the solutions are used to verify the results here.

We refer the readers to the appendices for details, and state the main results here: At $\textrm{ord}\left(\varepsilon^{0}\right)$, the disc-perturber interaction in $\mathscr{H}_{\rm p}$ induces a nodal precession of the perturber’s argument of the ascending node, $\Omega_{\rm p}$ with a frequency $\nu_{\rm p}\propto n_{\rm p}\frac{M_{\rm d}}{M_{\bullet}}$. This is non-zero even in the limit $\varepsilon\to 0$.

At the next order, the thinness of the disc introduces a very short time-scale, much shorter than $\nu_{\rm p}^{-1}$, over which the arguments of the ascending node of the disc stars align with $\Omega_{\rm p}$, with a phase difference of $90^{\circ}$. While this phenomenon appears at $\textrm{ord}\left(\varepsilon\right)$, it is much faster, and can be derived correctly only by explicitly accounting for this additional short time-scale. The short time-scale, $b_{{\rm p}n}$ is defined by

The equation of motion for $\Omega_{n}$ is

this equation has two time-scales: $\nu_{\rm p}=\dot{\Omega}_{\rm p}$, and $b_{{\rm p}n}$. Even though $b_{{\rm p}n}\propto\varepsilon$, it is still the case that $\left|b_{{\rm p}n}\right|/\left|\nu_{\rm p}\right|\gg 1$ because the disc is so thin; we show in appendix A.3 that for stars with $\left|b_{{\rm p}n}\right|>\left|\nu_{\rm p}\right|$, the solution to this equation yields that $\Omega_{n}$ quickly aligns such that

the alignment occurs after a time $\sim b_{{\rm p}n}^{-1}\ll\nu_{\rm p}^{-1}$, i.e. much faster than one nodal precession period. This is a striking phenomenon: the arguments of the ascending node of the stars in the disc with $\nu_{\textrm{p}}/b_{\textrm{p}n}\leq 1$ align themselves with that of the perturber, with a $90^{\circ}$ phase difference!

Stars for which $\left|\nu_{\rm p}/b_{\textrm{p}n}\right|>1$ behave quite differently: for them, the relevant limit of their governing differential equation, equation (40), is the oscillatory one, which implies that $\Omega_{n}(t)$ simply oscillates about $\Omega_{n}(0)$, regardless of $\Omega_{\rm p}$; no alignment occurs.

With these solutions for the arguments of the ascending node, one can solve the equations of motion for the inclinations, $i_{n}$ and $i_{\rm p}$; we do so in appendix C, and find that upon substituting the aligned $\Omega_{n}$ and $\Omega_{\rm p}$ the aligned stars in the disc exert a coherent torque on the perturber, which leads to a net decrease in its inclination. We show in appendix C that³³3Recall that $b_{{\rm p}n}$ depends on time via $i_{\rm p}$.

When substituting the solution (69) for $i_{n}(t)$, we find that at early times, $\mathrm{d}i_{\rm p}/\mathrm{d}t\propto-t/\sin i_{\rm p}$; the coefficient must, by dimensional analysis, have units of frequency squared, and we define the dynamical friction time-scale to be the reciprocal of the root of that coefficient, i.e., we define

at small $t$. This definition ensures that $\tau_{\rm RDF}$ gives a prediction for the time-scale over which $i_{\rm p}$ decreases significantly. We obtain

where

is the orbital period of the perturber and the proportionality constant is of order unity. We derive the proportionality coefficient in appendix C, but for a power-law surface density density profile $\Sigma\sim r^{-\beta}$ let us state the result here (see equation (73), and equation (72) for a general circular surface density profile):

where we have introduced the local disc mass at $r=a_{p}$

Note the mass scaling $\tau_{\rm RDF}\propto m_{\rm p}^{-1/2}$ but $di_{\rm p}/dt\propto m_{\rm p}$ (Eq. 11) as expected for dynamical friction.⁴⁴4These results are in agreement with the numerical results in figures 5 and 6 of Szölgyén et al. (2021) (cf. figure 8 below). In contrast, for Chandrasekhar dynamical friction, (Szölgyén et al., 2021)

This implies that the total alignment timescale for resonant dynamical friction is reduced by a factor of order $(m_{\rm p}M_{\rm d,loc})^{1/2}/M_{\bullet}$ in comparison.

We expect equation (11) to be valid under the following sufficient conditions: $\varepsilon\ll 1$ and the disc is thin relative to $\varepsilon$ (i.e., the perturber has not yet entered the disc); in other words, we need

We expect, however, that it will remain approximately accurate until the perturber enters the disc. The reason for that is, that the time-scale for the disc’s thickening (cf. equation (69)) is comparable to $\tau_{\rm RDF}$. The node alignment takes place on a time-scale $b_{{\rm p}n}(0)^{-1}$, which is shorter than $\tau_{\rm RDF}$ by a factor of $\sim\sin[i_{n}(0)]/\sqrt{\varepsilon}$.⁵⁵5In this case, our having defined $\tau_{\rm RDF}$ in the early time limit does not pose a problem. The coefficient of $t$ changes with time as more and more stars fall out of the nodal alignment, by equation (11). Hence, the total time until the perturber enters the disc is just the inverse of the right-hand-side of equation (12), where $t$ is evaluated as the time when a considerable fraction of the disc stars stop being aligned – which is just, again, $\tau_{\rm RDF}$.

To test these results, we perform a set of numerical simulations, which we now describe. We use a modified version of the direct $N$-body code phi-GRAPE (Harfst et al., 2007) which benefits from the 4th order Hermite integration scheme using block time-steps. Despite the original design for the GRAPE cards the code is adapted for modern GPUs and is widely used in various astrophysical simulations (see e.g. Li et al. 2019; Shukirgaliyev et al. 2021).

The simulations performed here fall into two categories: those with a disc and a perturber around the SMBH only, which are closer to the analytical model described above, and some more realistic ones which also include a ‘live’ spherical component. In addition, we also vary the initial eccentricity and inclination distributions of the disc, as well as the initial inclination of the perturber, and its mass (from $\varepsilon=1/8$ to $\varepsilon=1$). Details of the numerical set-up are described in appendix D; here, it suffices to show in figure 1 a plot of the various time-scales involved in the process, for the numerical set-up we consider. The parameters specified in appendix D and table 1 imply that in figure 1 $t_{\rm orb}=524$ years, and $\tau_{\rm RDF}=9347t_{\rm orb}\approx 4.9\times 10^{6}$ years. The expected synchronisation time-scale of the nodes is, on the other hand, much smaller than this $\tau_{\rm RDF}$, as shown in figure 1, for all stars that satisfy $\nu_{\rm p}<b_{{\rm p}n}$ initially.

As one expects the predictions of §3 to be valid until the perturber enters the disc, let us start by comparing them with those of the numerical simulations described above, when the perturber is still far from the disc, starting with those on the first row of table 1 – a thin, circular disc, with $\varepsilon=1/8$. A scatter-plot of the inclinations of all the stars in the disc, versus their semi-major axes, is shown in figure 2 – once derived from the analytical model of §3, and once from the numerical simulations. Similarly, we show a similar plot of the argument of the ascending node in figure 3. Following Szölgyén et al. (2021), we split the stars to 3 regions depending on the peri- and apocentres of the stars ($r_{p,*}$, $r_{a,*}$) with respect to the perturber ($r_{p,\textrm{p}}$, $r_{a,\textrm{p}}$): the inner region where the stellar orbits are within the pericentre of the perturber ($r_{a,*}<r_{p,\textrm{p}}$), the middle region where orbits of stars overlap with the perturber’s orbit ($r_{p,*}\leq r_{a,\textrm{p}}$ and $r_{a,*}\geq r_{p,\textrm{p}}$) and the outer region where the stellar orbits lie outside of the apocentre of the perturber ($r_{p,*}>r_{a,\textrm{p}}$). For the readers’ convenience, we also show a plot of the evolution of the average argument of ascending node of the disc within the regions as a function of time, compared with the inclination in figure 4. Indeed, figure 4 shows that the analytical predictions of $\dot{\Omega}_{\rm p}=\textrm{const}$ and $\Omega_{n}-\Omega_{\rm p}\approx\pi/2$ are satisfied by the numerical simulation not only at $t=1~{}\textrm{Myr}$, but up to about $4.5~{}\textrm{Myr}$, when that they persist until the perturber enters the disc.⁶⁶6Or, equivalently, until the disc is broken up by it – a situation which doesn’t occur for these models, but cf. figure 7.

Both the model of §3 and the simulation exhibit a sharp transition between an alignment of $\Omega_{n}$ with $\Omega_{\rm p}+\pi/2$, and an essentially uniform distribution of $\Omega_{n}$. According to §3 and appendix A.3, the transition occurs when $\left|\nu_{\rm p}\right|=\left|b_{\textrm{p}n}\right|$, i.e. when (recall that $M_{\bullet}\gg m_{\rm p},m_{n}$ for all $n$)

As $\nu_{\rm p}=\textrm{ord}\left(1\right)$, and $\frac{m_{\rm p}}{M_{\rm d}\sin i_{n}}\gg 1$, the transition can only occur at an $a_{n}>a_{\rm p}$, as is evident from figure 3.

A main difference in figures 2 and 3 between the analytical model and the simulations – the larger spread in values of orbital inclinations and longitudes of the ascending node – is caused by effects which are neglected in the analytical model. These may include two-body interactions or scalar resonant relaxation, but a detailed investigation of these is beyond the scope of this work.

One can see that the analytical treatment both captures the essential features of the simulations, and that the simulations do indeed exhibit an alignment of the nodes, as predicted by equation (10).

Next, we test the predictions of §3 (and appendix C) for $i_{\rm p}$: we show, in figure 5, the solution of equation (11), for the same numerical set-up as in figure 2, and compare it with the simulation. While the two do not completely agree, they do so approximately, and one can see that the treatment of §3 offers a way to understand the alignment of the inclination of the perturber with the disc.

We now proceed to show that the phenomena of §3 are not peculiar to the specific initial conditions considered above, that is to those of the first row of table 1. We vary both the mass-ratio $\varepsilon$ between the perturber and the disc, the latter’s thickness and the treatment of the spherical part of the system; we also change the perturber’s eccentricity and finally its initial inclination. The results are summarised in figures 6 and 7.

The first alteration we consider is turning the spherical component of the system into a ‘live’ sphere. The spherical component was a pure monopole, Plummer potential, above, and here, by sampling $N_{s}$ particles from a spherically symmetric power-law distribution, and allowing them to evolve and interact with all other particles, the ‘sphere’ could have non-zero higher multipoles, and a non-zero net angular momentum. This is the case with all the simulations shown in figures 6 and 7. From the first row of figure 6 we see that in this case the node alignment between the disc and the perturber still persists, especially with the stars that have a similar semi-major axis to $a_{\rm p}$, until the perturber coalesces with the disc. The effect is weakened, because the ‘sphere’ has a non-vanishing quadrupole moment by chance, which competes with the influence of the perturber. Upon increasing $m_{\rm p}$, the time it takes the perturber to reach the disc is decreased, but one can still see an alignment of the nodes before that. As the sphere becomes too massive (right panel of the penultimate row of the figures) the effect becomes somewhat less pronounced, and likewise for a thicker disc (bottom right panel).

Another noteworthy property is that, as one expects, if the perturber starts inside the disc (bottom right panel), there is no alignment, because then $\mathscr{H}_{\rm p}$ is far from dominating the dynamics.

The alignment time-scale in equation (13) is proportional to $m_{\rm p}^{-1/2}$. In figure 8 we test this dependence on the peturber’s mass, by running the same simulation as in the first row of table 1, but with varying IMBH mass. The total time it takes the perturber – starting at $i_{\rm p}(0)=45^{\circ}$ – to get to $i_{\rm p}=15^{\circ}$ is seen to follow the theoretical prediction closely.

Having shown that the perturber induces a rapid alignment of the $\left\{\Omega_{n}\right\}$ with its argument of the ascending node, let us inquire what the most general conditions under which one can expect this phenomenon to arise are; that is, are the requirements that $m_{\rm p}\ll M_{\rm d}$, and that the disc stars’ inclinations be much smaller than $\varepsilon$ really necessary? Does the phenomenon persist when including higher multipoles? What if the perturber is not a single particle, but its contribution to $\mathscr{H}$ is replaced by the shot noise fluctuations of a spherical stellar distribution, or a general $\mathscr{H}_{\rm p}$?

In the most general setting, the full Hamiltonian may be decomposed as

where $\mathscr{H}_{\rm d}$ governs the self-interaction of the disc (which we assume is $\mathscr{H}_{\rm LL}$), $\mathscr{H}_{s}$ describes the interaction of the perturber with itself, and $\mathscr{H}_{\rm p}$ includes all the interactions between the disc particles and the perturber(s). In the double orbit-averaged limit, the most general form for $\mathscr{H}_{\rm p}$ we consider here is

where the index $n$ runs over the disc particles, $p$ refers to the perturber(s), and $\ell\in\mathbb{N}$ is the multipole index, starting at quadrupole, $\ell=2$. The angle $\theta_{pn}$ is the angle between the angular momentum vector of perturbing particle $p$ and that of disc particle $n$. It is this angle which contains the entire dependence on $\Omega_{n}$, $i_{n}$, $\Omega_{p}$ and $i_{p}$. Let $\beta$ denote one of these four angles; we will restrict ourselves to the case where

where now $\theta_{n}$ is the angle between the angular momentum direction of particle $n$, and the total angular momentum of the perturbers $\hat{\mathbf{J}}_{\rm pert}$, i.e., we restrict ourselves to the case where the entire dependence on $p$ is in $\frac{\partial\theta_{pn}}{\partial\beta}$. This approximation is valid when the disc is thin. We may now write

in complete analogy with the inclination and argument of the ascending node of the single perturber above. By examining the way we derived equation (10), we see that it is necessary that $\dot{\Omega}_{\rm pert}\neq 0$, as well as that $\mathscr{H}_{\rm p}$ dominate the dynamics of the disc particles. This is possible for the thin disc limit considered above. If these two conditions hold, the derivation can be repeated, and a rapid alignment occurs. For the single perturber case, as $\varepsilon$ increases, both of these conditions are still met, but the disc thickens faster, over a time-scale $A_{n}$ (which is explicitly given in appendix C), and we know that the alignment persists only until the disc’s thickness reaches $\textrm{ord}\left(\varepsilon\right)$. Consequently, as the alignment occurs over a time-scale $b_{\textrm{p}n}^{-1}$, and the thickening requires time $A_{n}$, one must have

as a necessary condition for $\Omega_{n}$ to align itself with $\Omega_{\rm p}$ in accordance with (10), i.e. precisely that the disc be thin, i.e. that inequality (18) be satisfied. These are the conditions under which such a phenomenon occurs, and it persists as long as they are met.⁷⁷7For $\varepsilon\sim 1$, these conditions could be met, but the disc’s thickness would grow too much over the course of its evolution, and hence the assumption of $\mathscr{H}_{\rm d}=\mathscr{H}_{\rm LL}$ would have to be modified at late times. Furthermore, it is also possible that for the outer extent of the disc, one would not have $\mathscr{H}_{\rm p}\gg\mathscr{H}_{\rm LL}$, and indeed, whether or not this alignment occurs varies from star to star.

Very recently Levin (2022) studied a rotating spherical cluster and a (rotating) disc, which quickly align together around an SMBH; this alignment was derived from a statistical mechanical view-point using the fluctuation-dissipation theorem. If the perturber were replaced by a rotating cluster, then the discussion in the previous paragraph would apply. Indeed, the mass-scaling of the resonant friction time-scale derived by that work is the same as that in equation (13).

This time-dependence is much faster that the ordinary $\tau_{\rm CDF}\sim M_{\bullet}^{2}/(m_{\rm p}M_{\rm d})\times t_{\rm orb}/\ln\Lambda$ one would find for Chandrasekhar dynamical friction. This implies that $\tau_{\rm RDF}\sim\sqrt{\tau_{\rm CDF}\,t_{\rm orb}\ln\Lambda}$ – generally much shorter than $\tau_{\rm CDF}$. For our setting, $M_{\bullet}=4\times 10^{6}M_{\odot}$, $M_{\rm d}=2000M_{\odot}$, $m_{\rm p}\approx 250M_{\odot}$, $t_{\rm orb}=524$ years (for the same density-profile as in §4 and in figure 1), we obtained $\tau_{\rm RDF}\approx 4.9~{}\textrm{Myr}$, while $\tau_{\rm CDF}\approx 1.7~{}\textrm{Gyr}$ for $\ln\Lambda=10$. For smaller Coulomb logarithms, the latter becomes already of the order of the age of the Universe. For even more massive super-massive black holes, the difference is even more extreme: for instance for $m_{\rm p}=1000M_{\odot}$, $M_{\bullet}=10^{9}M_{\odot}$, and $M_{\rm d}=10^{4}M_{\odot}$, the orbital time is multiplied by $\sqrt{10^{9}/(4\times 10^{6})}$, and we find $\tau_{\rm RDF}\approx 4.3~{}\textrm{Gyr}$ while $\tau_{\rm CDF}\approx 8.3\times 10^{13}~{}\textrm{years}\gg H_{0}^{-1}$. In other words, an alignment of the inclination due to resonant dynamical friction is possible in systems, where ordinary, Chandrasekhar, dynamical friction would take much too long to do so.

In this paper we developed an analytical model, based on resonant relaxation and singular perturbation theory, which captures the essential features of resonant dynamical friction. We showed that the singular nature of the equations of motion implies, in the case of a thin disc, a rapid alignment of the arguments of the ascending node of the perturber and the disc particles, which then gives rise to a coherent torque, which aligns the perturber with the disc. We showed that the predictions of this model mesh well with results of $N$-body simulations, and found that the node alignment occurs for a wide range of initial conditions. In particular, one may safely conclude that the perturber would not re-orient to the disc if the disc particles were treated as test-particles, for it is their gravity that sources an $\textrm{ord}\left(1\right)$ value of $\dot{\Omega}_{\rm p}$, but on the other hand, the perturber’s contribution to the potential does dominate their motion, at least initially.

The instantaneous alignment timescale $(\mathrm{d}i_{\rm p}/\mathrm{d}t)^{-1}$ is proportional to $m_{\rm p}^{-1}$, but the total alignment time-scale as a function of SMBH, IMBH, and local disc mass $(M_{\bullet},m_{\rm p},M_{\rm d,loc})$ is proportional to $M_{\bullet}/(m_{\rm p}M_{\rm d,loc})^{1/2}$ – much faster than Chandrasekhar dynamical friction timescale which scales as $M_{\bullet}^{2}/(m_{\rm p}M_{\rm d,loc})$. We expect some of the results of this analysis to extend to more general perturbers of the disc, such as a ‘live’ spherical component instead of a single point-particle perturber, or the case of two stellar discs, or to certain planetary systems.

We wish to thank the anonymous referee for reviewing our paper insightfully. We are grateful to Yuri Levin and Mor Rozner for helpful comments on the manuscript. Y. B. G. is grateful for the kind hospitality of the Rudolf Peierls Centre for Theoretical Physics at the University of Oxford and to St. Hugh’s College, Oxford, where some work on this research was done. Y. B. G. is supported by the Adams Fellowship Programme of the Israeli Academy of Sciences and Humanities. T. P. and B. K. received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 Programme for Research and Innovation ERC-2014-STG under grant agreement No. 638435 (GalNUC), and were also supported by the Science and Technology Facilities Council (STFC) Grant Number ST/W000903/1. T. P. acknowledges the support of the Science Committee of the Ministry of Science and Higher Education of the Republic of Kazakhstan (Grant No. AP14870501).

The data underlying this article will be shared on reasonable request to the corresponding author.