<figure><img src="image_1.png" /></figure>

HAL
open science

# Dynamical traction and black hole orbital migration

Christian M Boily, Thibaut L François, Jonathan Freundlich, Francoise

Combes, Anne-Laure Melchior, Y. Hénin

<figure><img src="image_2.png" /></figure>

► To cite this version:

Christian M Boily, Thibaut L François, Jonathan Freundlich, Francoise Combes, Anne-Laure Melchior, et al.. Dynamical traction and black hole orbital migration. Astronomy & Astrophysics - A&A, 2025, 701, pp.A10. 10.1051/0004-6361/202453104. hal-05237438

HAL Id: hal-05237438

https://hal.science/hal-05237438v1

Submitted on 3 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Dynamical traction and black hole orbital migration

## I. Angular momentum transfer and a fragmentation-driven instability*

C. M. Boily<fnref n="1" />,**, D, T. L. François<fnref n="1" />,<fnref n="2" />, J. Freundlich<fnref n="1" />, F. Combes<fnref n="3" />,<fnref n="4" />, A.-L. Melchior<fnref n="3" />, and Y. Hénin<fnref n="1" />,<fnref n="5" />

Received 21 November 2024 / Accepted 7 July 2025

## ABSTRACT

*Context.* Observations with the James Webb Space Telescope of massive galaxies at redshift $z$ up to ~14, many of which have quasar activity, require careful accounting of the orbital migration of seed black holes to the heart of the host galaxy on timescales of 300 Myr or less.

*Aims.* We investigate the circumstances that allow a black hole to remain at the barycentre of a system when the equilibrium galactic stellar core is anisotropic.

*Methods.* We developed a Fokker-Planck treatment to analyse the migration of a massive black hole. The analysis focuses on ex-
changes of a black hole’s orbital angular momentum with stars. Furthermore, we used a set of N-body calculations to study the
response of stellar orbits drawn from a Miyamoto-Nagai disc embedded in a larger, isotropic isochrone (Hénon) background poten-
tial.

*Results.* When the black hole has little angular momentum initially but orbits in a sea of stars drawn from an odd $f[E, L_z]$ velocity distribution function, a wake in the stellar density sets in that pulls on the black hole and transfers angular momentum to it. We call this transfer ‘dynamical traction’, in contrast with the more familiar Chandrasekhar dynamical friction. We argue that this phenomenon takes place whenever the kinetic energy drawn from $f[E, L_z]$ has an excess of streaming motion over its (isotropic) velocity dispersion. We illustrate this process for a black hole orbiting in a dynamically warm disc with no sub-structures. We then show that for a dynamically cold disc, the outcome depends on both the orbit of the black hole and that of the stellar sub-structures stemming from a Jeans instability. When the stellar clumps have much binding energy, a black hole may scatter off of them after they are formed. In the process, the black hole may be dislodged from the centre and migrate outward due to dynamical traction. When the stellar clumps are less bound, they may still migrate to the centre, where they either dissolve or merge with the black hole. The final configuration is similar to a nuclear star cluster that may yet be moving at ~10 km s with regard to the barycentre of the system.

*Conclusions.* The angular momentum transferred to a black hole by dynamical traction delays the migration to the galactic centre by several hundred million years. The efficiency of angular momentum transfer is a strong function of the fragmented (cold) state of the stellar space density. In a dynamically cold environment, a black hole is removed from the central region through a two-stage orbital migration instability. A criterion against this instability is proposed in the form of a threshold in isotropic velocity dispersion compared to streaming motion (i.e. angular momentum). For a black hole to settle at the heart of a galaxy on timescales of ~300 Myr or less, a large fraction of the system angular momentum must be dissipated or, alternatively, the black hole must grow in situ in an isotropic environment devoid of massive sub-structures.

**Key words.** galaxies: active – galaxies: bulges – galaxies: formation – galaxies: kinematics and dynamics – quasars: general – quasars: supermassive black holes

## 1. Introduction

In this article, we investigate the orbital evolution of a black hole
moving freely at the heart of a galaxy in the later stages of a
major merger. The motivation for this work stems from recent
observations by the James Webb Space Telescope (JWST) of
galaxies in the mass range of $10^9-10^{10}$ solar-masses at cosmo-
logical redshifts up to $z \approx 12$ (Naidu et al. 2022; Labbé et al.
2023; Finkelstein et al. 2023). The short $\le 500$ Myr timescales
implied challenge standard models of the formation of galax-
ies and super-massive black holes (SMBH; see e.g. Mo et al.

2010; Boylan-Kolchin 2023; Matthee et al. 2024; Barro et al. 2024) because they point to a high level of efficiency of the formation processes. For example, the rapid growth of SMBH on a timescale of a few hundred million years by accretion or coalescence suggests that seed BHs of $10^4 M_{\odot}$ or more form already at the epoch of recombination (Bogdán et al. 2024). In these scenarios, the growth of a single BH by gas accretion at the Eddington limit or by stochastic shredding and loss-cone capture of main-sequence stars is made more likely by placing the BH in the pit of the (local) gravitational well, where gas and stars are accreted on shorter timescales (e.g. Mo et al. 2010; Schawinski et al. 2015; Pfister et al. 2021; see also the review by Volonteri 2012). However, this basic picture bypasses

<footnote n="1">¹ Observatoire astronomique de Strasbourg, Université de Strasbourg, CNRS UMR 7550, 10 rue de l'Université, Strasbourg F-67000, France</footnote>

<footnote n="2">² School of Mathematics and Physics, University of Surrey, Guildford, Surrey GU2 7XH, UK</footnote>

<footnote n="3">³ LERMA, Sorbonne Université, Observatoire de Paris, Université PSL, CNRS, F-75014 Paris, France</footnote>

<footnote n="4">⁴ Collège de France, 11 Place Marcelin Berthelot, 75005 Paris, France</footnote>

<footnote n="5">⁵ Télécom Physique Strasbourg, Pôle API, 300 Bd Sébastien Brant, F-67412 Illkirch Graffenstaden Cedex, France</footnote>

<footnote>* Dedicated to the memory of our colleague T. Padmanabhan.</footnote>

<footnote>** Corresponding author: christian.boily@astro.unistra.fr</footnote>
the inevitable mass build-up of the host galaxy through multiple mergers. At the same time as a BH accretes from its immediate surroundings, it drifts in space in a time-dependent galactic potential. Several recent studies have included BH wandering in their dynamical modelling (e.g. Bellovary et al. 2019; Pfister et al. 2019; Lapiner et al. 2021). Bright cluster galaxies in the Illustris-TNG cosmological simulations point to SMBHs wandering over a timescale of order 1 Gyr (Chu et al. 2023). The expectation that an SMBH will come to settle at the barycentre of the host galaxy rests on the theory of dissipation of orbital energy by dynamical friction (Chandrasekhar 1943; Binney & Tremaine 2008) and from observations of local Universe AGN galaxies, where it almost invariably sits at the photometric centre of the galactic bulge. The Milky Way's SMBH, which is of a relatively low mass ($\approx 4 \times 10^6 M_\odot$; see e.g. Chatzopoulos et al. 2015), is a prime example of this configuration. Yet the recent detection of an off-centred SMBH at $z \approx 7.15$ by Übler et al. (2024) puts this expectation in the spotlight and is cause for reassessment of the conditions that enable an SMBH to stay at rest at the heart of a young galaxy. We investigate how an SMBH can remain at the barycentre of a cored galaxy when it is surrounded by a disc of stars (i.e. stars supported by angular momentum). This configuration is designed to represent the final stage of an in-falling SMBH carried by a major (dry) galaxy merger in an advanced stage with much orbital angular momentum still present in the surrounding stars (see for example Chapon et al. 2013). Thus, this setup complements earlier work exploring an SMBH wandering in an isotropic background of stars (de Vita et al. 2018; Merritt 2013). Our focus is on the orbital evolution of an SMBH accreted by a young galaxy with an anisotropic flow of stars with net rotation (a proto-galactic disc). We view the SMBH as an external $m = 1$ perturbation mode acting on a system in near-equilibrium. The expectation for an accretion event carrying little angular momentum is that it will pick up angular momentum from background stars (as the perturber has a slow pattern speed; see Tremaine & Weinberg 1984), and this will slow down its migration to the centre. We focus on the case of an SMBH orbiting in the same plane as the proto-galactic disc, and our goal is to quantify its orbital migration under gravity alone (no gas dissipation). Since gravitational dynamics are scale-free, hereafter we use the acronym BH as a generic term.

Our study should be put in the perspective of earlier works. In a system with periodic orbits and a flat (cored) stellar mass density, all orbits have periods in the ratio $\approx 1:1$. The breakdown of resonant motion between the perturber and the background stars near the system barycentre leads to core stalling (e.g. Read et al. 2006; Boily et al. 2008; Cole et al. 2012; Kaur & Sridhar 2018; Banik & van den Bosch 2022) when dynamical friction shuts off entirely (see also Petts et al. 2015; Banik & van den Bosch 2021; Chiba & Kataria 2024, for orbit-based analyses). This leaves the orbital energy and angular momentum of a massive perturber unchanged. Kalnajs et al. (1972) showed that dynamical friction in a self-gravitating disc of uniform angular speed $\Omega$ vanishes everywhere. His analysis of global modes of perturbation extends to individual orbits when the disc is fully phase mixed. The perturbation of stellar orbits of any energy by a massive BH leads to positive and negative torques acting on the BH, as the stars respond to the perturbation (the BH was dubbed a catalyst of energy exchange between stars in Boily et al. 2008). Recently, Banik & van den Bosch (2022) and Kaur & Stone (2022) showed that the exact cancellation of gravitational torques depends on details of the phase-space distribution function, $f$. The general conclusion is

that dynamical friction will be much reduced in any realistic near-constant density environment with the BH on a circular orbit.

The case of a rotating system with a disc-like ordered motion is different because of the strongly anisotropic velocity field. The response of the stars to an in-falling BH leads to the growth of density patterns, as first shown by Julian & Toomre (1966), where a first-order perturbation analysis developed exponential growth in the sheared-box approximation. In a galactic disc in equilibrium, spiral density waves cause mechanical torques that redistribute angular momentum throughout the disc (e.g. Lynden-Bell & Kalnajs 1972; Tremaine & Weinberg 1984; Hopkins & Quataert 2011; Anglés-Alcázar et al. 2015). Tremaine & Weinberg (1984) performed a general analysis based on resonant orbits and argued that friction originates from global resonant modes exerting a net torque on a massive body. The response of the stars remains correlated in phase-space to the perturber's orbit, which is then said to be dressed by these response modes (see Weinberg 1998; Fouvry et al. 2015; and especially Lau & Binney 2021 for a detailed analysis). A strong response by the stars leads to polarised orbits. The self-gravity of the perturbed stars triggers non-linear effects (modes couple to each other) that are difficult to track analytically. Indeed we show with numerical examples how a fragmented (dynamically cool) stellar disc has a much different impact on the BH than a warm one.

We review the basics of dynamical friction in Sect. 2 before developing our perturbation analysis based on a Fokker-Planck treatment in Sect. 3. Having obtained some analytical direction, we then move on to present a series of restricted N-body models that detail the time evolution of the BH's perturbed orbit (Sect. 4). We show in Sect. 5 that the orbital migration of the BH triggers the growth of $m = 1$ and even $m = 2$ (bar) modes in the stars. For such modes, further evolution of the orbit through gravitational torques becomes possible. A brief quantitative comparison of the numerical results with theoretical expectations is presented in Sect. 6. We discuss the implication of the perturbed configurations for the short-term accretion history of the BH in Sect. 7.

## 2. Dynamical friction

The main driver of BH orbital evolution is the dynamical friction exerted by low-mass field stars (and gas) on the BH. The density wake forming behind a BH of mass $M_\bullet$ moving at $v_\bullet$ causes a frictional acceleration of the BH (Chandrasekhar 1943; Binney & Tremaine 2008, Sect. 8.1):

$$ \frac{d\mathbf{v}_\bullet}{dt} = \langle \Delta\mathbf{v}_\bullet \rangle_\parallel = -\frac{\Gamma M_\bullet}{\pi m_\star} \int d^3\mathbf{v} f(\mathbf{v}) \frac{\mathbf{v}_\bullet - \mathbf{v}}{|\mathbf{v}_\bullet - \mathbf{v}|^3}, \quad (1) $$

where $M_\bullet$, $v_\bullet$ are the BH's mass and velocity; $\mathbf{v}$ the background stars velocity of phase-space density $f(\mathbf{v})$, and the integral is performed over all bound stars (orbital energy $E < 0$). The Coulomb term reads $\Gamma = 4\pi G^2 m_\star^2 \ln \Lambda$ with

$$ \Lambda = \frac{p_{\text{max}}}{p_{90^{\circ}}} \approx 0.3N, $$

where $p$ is the impact parameter evaluated for a hyperbolic stellar orbit in the star-BH barycentre, and $\Lambda$ is the ratio of the maximum value to the value leading to a 90° deflection angle (e.g. Spitzer 1987; Heggie & Hut 2003). The last relation to the number of stars $N$ applies for a system of bounding volume $V$ and uniform mean mass density, $\rho_\star = Nm_\star/V$.
<figure><img src="image_4.png" /><figcaption>Fig. 1. Sketch of gravitational accelerations for two different situations. The panels show a massive body (blue dot) orbiting among stars, and the system's barycentre at the bottom is the focus of the global gravitational acceleration $g_{cm}$ (open circle). The mean field acceleration gives rise to a velocity component $\delta v_\bullet$ along the vertical y-axis. In panel (a), a massive perturber moves at some velocity $v_\bullet$ in an isotropic field of stars; the result of the polarisation of the stellar orbits is a net dynamical friction force $a_\parallel$, parallel but in the opposite sense to $v_\bullet$. In panel (b), the massive body is ejected from the origin at a velocity $\delta v_\bullet$ pointing upwards. It orbits in a strongly anisotropic (rotating) stellar field and slows down due to the steady negative acceleration $g_{cm}$. However, dynamical friction now contributes two terms, one $a_\parallel$ as before, and a second one, $a_\perp$, orthogonal to $\delta v_\bullet$. That non-zero component due to the streaming stars causes a steady transfer of angular momentum to the perturber. We refer to this second component as dynamical traction acting on the perturber.</figcaption></figure>

As a consequence of Eq. (1), when the field stars are distributed isotropically in velocity space (i.e. when $f[v] = f[v]$ is a function of the norm $v$ alone) the net deceleration runs parallel to $v_\bullet$ as depicted on Fig. 1(a). The steady loss of kinetic energy brings the massive body at the barycentre of the galaxy, the only stable point in equilibrium for a body at rest. The time interval for this to take place is estimated from

$$ \Delta t_{\text{fric}} = \frac{6}{5 \ln \Lambda} \frac{M(<r)}{M_{\bullet}} \times \sqrt[2]{\frac{3\pi}{G\rho_{\star}}} = \frac{6}{5 \ln \Lambda} \frac{M(<r)}{M_{\bullet}} t_{\text{dyn}}, \quad (2) $$

with the background stars of density $\rho_\star$ adding to a mass inside the orbital radius, $M(<r)$, much larger than the massive body, $M_\bullet$. We also introduced a definition of the dynamical time, $t_{dyn}$. (In numerical applications, the ratio $\Delta t_{fric}/t_{dyn} \approx 3$; see Sect. 4.)

The type of kinetic energy diffusion that results from (1) is drawn from summing over non-resonant star-BH hyperbolic encounters. The unperturbed stellar trajectories are straight lines and vanish to infinity after the encounter with the perturber. In reality, the host galaxy spawns a finite volume, and stars with a similar orbital energy visit each other and the perturber BH repeatedly. This causes resonant coupling between the orbits when their period ratio is commensurate (Binney & Tremaine 2008). Such a resonant response leads to efficient energy exchanges and the notion of buoyancy and core-stalling (Read et al. 2006; Banik & van den Bosch 2022; Kaur & Sridhar 2018). We revisit this subject briefly in Sect. 7. In the next paragraph Sect. 3, we outline the coupling between the BH orbit and the background stars based on the non-resonant Fokker-Planck treatment.

### 3. Position of the problem

To achieve our goal, several simplifying assumptions are required. The first is that a well-defined galactic rotation pattern is in place. We also presume that the BH has already bound to the

host system (no dissipative terms, no gas) and remains embedded in the rotating flow of stars. Also, the BH is a small fraction of the total mass (less than 1% of the stellar disc, but ≈0.1% of the system, in line with observational trends; Kormendy & Ho 2013) and much more massive than the background stars. The problem then is to identify a small set of initial conditions for the BH's orbit that are idealised yet useful to our analysis.

#### 3.1. Orders of magnitude calculation

Let a BH sit at rest at the origin of coordinates surrounded by a spheroid of stars drawn from a distribution function $f(E, L_z)$, where $E$ is the energy and $L_z$ the z-component of the angular momentum. For simplicity, we lump together the undressed BH mass and the stars strongly bound to it (when stellar orbits evolve adiabatically in the BH potential, e.g., Binney & Tremaine 2008; Spitzer 1987, Sect. 2.3). The stars give rise to an axi-symmetric potential $\Phi(R, z)$ and the whole system is in equilibrium (à la Tremaine et al. 1994). In the context of a young galaxy in an advanced stage of formation, some sub-structures of significant mass may still induce time-dependant perturbations in the global potential. We set the BH in motion through an $m = 1$ perturbation mode, which we model as a hyperbolic (shock) encounter with a body of mass $M_p$, impact parameter $p$, passing by the BH at velocity $V_p$ parallel to the z-axis (Spitzer 1987, Sect. 5.2). The velocity imparted to the BH is then

$$ \|\delta v_{\bullet}\| \approx \frac{2GM_p}{pV_p}, \quad (3) $$

with a velocity vector, $\delta v_\bullet$, that we align with the y-axis of a Cartesian frame of reference (see Fig. 1[b]). When the (dynamically hot) spheroid of stars spawns a uniform density $\rho_\star$, the orbital motion is flagged with a reference angular frequency $\omega = \sqrt{4\pi G\rho_\star/3}$, and the amplitude of the BH orbit can be
estimated from the solution for a harmonic oscillator,

$$ \delta\mathbf{x}_{\bullet} \approx \frac{\delta\mathbf{v}_{\bullet}}{\omega}. \qquad (4) $$

The impact parameter, *p*, and the relative velocity, $V_p$, in Eq. (3) are free for us to choose. However, the mass, $M_p$, can be constrained by requiring that the perturbing body survives the encounter. For an extended body (such as a clump of gas or a star cluster), this will be achieved if $M_p$ exceeds the local Jeans mass (Binney & Tremaine 2008, Sect. 3):

$$ M_p \ge \frac{4\pi}{3} \rho_{\star} \left(\frac{\lambda_J}{2}\right)^3 = \frac{\pi^{\frac{5}{2}}}{6} \left(\frac{\sigma_{\star}}{\sqrt{G\rho_{\star}}}\right)^3 \rho_{\star} \equiv M_J. \qquad (5) $$

In this equation the length $\lambda_J$ is the Jeans wavelength and $\sigma_\star$ is the stars' one-dimensional velocity dispersion. A massive BH will generate a strong tidal field in a volume where it contributes much of the mass; furthermore, the impact parameter, *p*, should be greater than the size of the perturbing body, $M_p$; otherwise, it could overlap with the BH at the closest approach. For these reasons, in numerical applications we set

$$ p \ge \max \left\{ \frac{\lambda_J}{2}, l_{\bullet} \right\}, \qquad (6a) $$

$$ l_{\bullet} = \left[ \frac{3}{4\pi} \frac{M_{\bullet}}{\rho_{\star}} \right]^{\frac{1}{3}} \qquad (6b) $$

for self-consistency. The length $l_\bullet$ is sometimes referred to as the BH radius of influence (as in Merritt 2004 up to a factor 2; $l_\bullet$ should not to be confused with the event horizon of relativistic gravitation).

We considered a naked BH mass but included all stars with a binding energy (to the BH) exceeding the energy perturbation applied to the BH itself. For stars in a quasi-Keplerian orbit around the BH, this rough approximation implies that the BH mass includes stars confined to a Bahcall & Wolf (1976) cusp where $\rho_\star \propto r^{-7/4}$. We discuss the implications of this simplification in the closing section (Sect. 8). Below we derive expectations from a Fokker-Planck treatment of the diffusion of kinetic energy. Mathematical details can be found in Appendix A.

## 3.2. Sample distribution function

A Fokker-Planck treatment requires a choice of the phase-space distribution function $f(\mathbf{x}, \mathbf{v})$. We applied a Dirac $\delta$-function to select loop orbits otherwise drawn randomly from $f$. The distribution function (d.f.) takes the form

$$ f(\mathbf{x}, \mathbf{v})|_{\mathbf{x}} \equiv f_{\mathbf{x}}(\mathbf{v}_{\star}) = f_{\perp}(\mathbf{v}_{\perp}) \mathcal{D}(v_{\phi} - R\Omega), \qquad (7) $$

where $R$, $\phi$, and $z$ are cylindrical coordinates; $\Omega(R, z = 0)$ is the local angular frequency; and $\mathcal{D}(v)$ is the one-dimensional Dirac operator (with units of $v^{-1}$; see e.g. Fukushige & Heggie 2000; Renaud et al. 2011, for further applications of the Dirac operator). The $R, z$ components of the velocity field are made to be irrotational and isotropic by fixing

$$ f_{\perp}(\mathbf{v}_{\perp}) = \frac{\rho(R, z)/m_{\star}}{2\pi\sigma_{\star}^2} \exp\left(-\frac{v_{\perp}^2}{2\sigma_{\star}^2}\right). \qquad (8) $$

We set the azimuthal velocity component of the BH, $v_\phi$, equal to the parallel component, $v_\parallel$, in Eqs. (A.3a-c), so

$$ v_{\phi} \approx v_{\parallel} $$

throughout. This will prove sufficiently accurate to recover the evolution with the restricted N-body numerical models of Sect. 4.

With the definition

$$ \mathcal{K}^2 \equiv (R\Omega - v_\phi)^2 + v_\perp^2 \qquad (9) $$

and the help of Eqs. (A.3), (A.6) and the relations (7) and (8), we found

$$ \langle \Delta v_{\phi} \rangle = \Gamma \frac{M_{\bullet} v_{\phi}}{m_{\star} v} f_{\perp}(0) [\mp 2\pi + xe^{x^2} \operatorname{erfc}(x)] \times \frac{d\mathcal{K}}{dv}, \qquad (10) $$

where $x \equiv \mathcal{K}\sqrt{\pi/2\sigma_\star^2}$ and $\operatorname{erfc}(x)$ is the complementary error function. The key feature of Eq. (10) is that the contribution to the variation in kinetic energy $v_\phi \langle \Delta v_\phi \rangle$ changes sign with $R\Omega - v_\phi$ (see Eq. [A.7]). It is the exact same feature of gain or loss of kinetic energy depending whether the orbit trails or leads a resonant density pattern, as found in earlier studies (Tremaine & Weinberg 1984; Rauch & Tremaine 1996; Weinberg 1989; Kocsis & Tremaine 2015). We do not treat the growth of resonant features here, as we focus only on the response of the background stars and its impact on the BH's orbit.

## 3.3. The transition to dynamical traction

Figure 2 plots the rates of specific kinetic energy variation entering the right-hand side of Eq. (A.7) as function of the azimuthal velocity component, $v_\phi$. For the reference calculations displayed on Fig. 2(a), the Coulomb potential was set to $\ln \Lambda = 21$ ($N = 10^5$) and $m_\star = 1.6 \times 10^4 \text{ M}_\odot$, $M_\bullet = 1 \times 10^7 \text{ M}_\odot$; the rotational velocity $R\Omega \approx 70 \text{ km s}^{-1}$ at cylindrical radius $R = 1.5 \text{ kpc}$ fixes the angular frequency. Two sets of results are shown, obtained for a stellar velocity dispersion $\sigma_\star = 10 \text{ km s}^{-1}$ (in blue) and $20 \text{ km s}^{-1}$ (in black).

The high sensitivity to $\sigma_\star$ is clearly visible for the quadratic terms $\langle (\Delta v_\perp)^2 \rangle$ and $\langle (\Delta v_\parallel)^2 \rangle$. With $\sigma_\star = 10 \text{ km s}^{-1}$, both curves rise steeply at low $v_\phi$, and reach a maximum value well off-scale on the figure. Thus at low $v_\phi$, the gain or loss of kinetic energy is completely dominated by random scattering. The expectation, therefore, is that unless the BH is on an orbit with a significant azimuthal velocity, it is unlikely to acquire further streaming motion. Instead, a gain in radial and vertical velocity would set it on a new orbit of increased eccentricity (no gain in angular momentum) at a different energy, $E$. The turnover from a gain in radial- to one in azimuthal kinetic energy occurs at $v_\phi \approx 45 \text{ km s}^{-1}$, when the rate $v_\parallel \langle (\Delta v_\parallel)^2 \rangle \approx 20(\text{km s}^{-1})^2/\text{Myr}$ reaches a much larger amplitude than the quadratic terms. By contrast, when the velocity dispersion $\sigma_\star = 20 \text{ km s}^{-1}$, the rates $\langle (\Delta v_\perp)^2 \rangle$ and $\langle (\Delta v_\parallel)^2 \rangle$ peak at a much lower value of $v_\phi \approx 6 \text{ km s}^{-1}$, and decrease rapidly thereafter. The gain in azimuthal kinetic energy becomes dominant at $v_\phi \approx 20 \text{ km s}^{-1}$ and larger values, which translates in a net gain in azimuthal motion from that point on-ward. This behaviour is expected because a warmer stellar population will be harder to polarise, and so the effect of random scattering off of the BH orbit is much reduced (less effective gravitational focusing). The streaming flow remains the same, however, so the response in azimuth is much less affected by variations in $\sigma_\star$. We say that the gain in angular momentum is triggered by a dynamical traction transition<fnref n="1" /> because for some

<footnote n="1">¹ We have not shown that this transition takes place through a topological critical point in solution space; hence, we do not refer to this dynamical process as a bifurcation. We thank the referee for bringing this to our attention.</footnote>
<figure><img src="image_5.png" /><figcaption>Fig. 2. Rates of specific kinetic energy diffusion, $\langle \Delta E \rangle$, as a function of the azimuthal velocity, $v_\phi$. The rates are evaluated from Eqs. (A.7) and (A.3) for two values of the stellar velocity dispersion, $\sigma_\star \cdot 20 \text{ km s}^{-1}$ (in black or grey); and $10 \text{ km s}^{-1}$ (in blue or skyblue). A) Setting $M_\bullet = 1 \times 10^7 \text{ M}_\odot$, with $M_\bullet/m_\star \approx 10^3$, the figure graphs each quadratic component and the first-order parallel one ($v_\parallel \langle \Delta v_\parallel \rangle$, dark solid lines; we note that it has the same sign as $R\Omega - v_\phi$). B) The same but for a BH with $M_\bullet = 1 \times 10^8 \text{ M}_\odot$ ($M_\bullet/m_\star \text{tar} \approx 10^4$). The much higher inertia of the BH leads to a ten-fold increase in the parallel term but only subtle differences in the diffusion (quadratic) ones. The light and dotted-blue curves peak off scale in either case. The (central) angular velocity $\Omega = 0.05/\text{Myr}$ is such that $R\Omega \approx 70 \text{ km s}^{-1}$ at $R = 1.5 \text{ kpc}$, a good match to the circular velocity at that radius.</figcaption></figure>

critical value of $v_\phi$, the orbit of the BH gains angular momentum systematically, transforming its orbit from a low-$L_z$ box orbit to a high-$L_z$ loop. The onset of this seemingly irreversible trend is a sensitive function of the velocity field of the background stars. We cover this point more fully with applications later in Sect. 3.5.

Another important parameter is the choice of BH mass. On Fig. 2(b), we graph the results for the same 1D velocity dispersions $\sigma_\star$, for a BH of mass $M_\bullet = 1 \times 10^8 \text{ M}_\odot$. The ten-fold increase in $M_\bullet$ means that the rate $v_\parallel \langle \Delta v_\parallel \rangle$ also increases ten times (see the multiplicative factor in Eq. [A.7]). The most noticeable effect is a slight shift to the left of the curves for $\langle (\Delta v_\perp)^2 \rangle$ and $\langle (\Delta v_\parallel)^2 \rangle$, which both cross the run of $v_\parallel \langle \Delta v_\parallel \rangle$ at lower values of $v_\phi$. For example, for the case with $\sigma_\star = 20 \text{ km s}^{-1}$ the turnover to a regime where $v_\parallel \langle \Delta v_\parallel \rangle$ is dominant takes place at $v_\phi \approx 12 \text{ km s}^{-1}$ (it was $\approx 20 \text{ km s}^{-1}$ on Fig. 2[a] for the lower-mass BH). This suggests that the more massive BH may yet gain angular momentum from a more slowly rotating stellar flow (all other quantities being unchanged). To back up this intuition, we have checked that reducing the angular frequency $\Omega$ from $\approx 7.0 \times 10^{-1} \text{ Myr}^{-1}$ to $3.5 \times 10^{-1} \text{ Myr}^{-1}$ brings the crossover value to $v_\phi \approx 20 \text{ km s}^{-1}$, when $\sigma_\star = 10 \text{ km s}^{-1}$; and to $v_\phi \lesssim 5 \text{ km s}^{-1}$, when $\sigma_\star = 20 \text{ km s}^{-1}$.<fnref n="2" /> The crossover value goes down by a factor of more than two from the reference calculation, highlighting the strong link with the rotation curve. The calculations with the reduced angular speed $\Omega$ still had a circular velocity $\approx 35 \text{ km s}^{-1}$ at $R = 3/2 \text{ kpc}$, significantly larger than the two values for the velocity dispersion $\sigma_\star$ that we tested. A reduced $\Omega$ implies longer timescales for the onset of the dynamical traction transition (i.e. when the net gain in angular momentum will exceeds the loss to dynamical friction, and transform the orbit of the BH to a loop orbit).

There is no gain or loss in azimuthal energy from dynamical traction when $v_\phi = R\Omega$ (co-moving BH) because $v_\parallel \langle \Delta v_\parallel \rangle = 0$ there. In the event of the BH running ahead of the stars, the rate $v_\parallel \langle \Delta v_\parallel \rangle < 0$ changes sign. If and when that is the case, the BH loses kinetic energy and angular momentum, when it was gaining before. This is the signature of an 1:1 resonant trapping at co-rotation. That result was obtained by an argument of symmetry, rather than a proper linear perturbation treatment of the

stellar response. No such analysis was needed here because we focused on the perturber (the BH) and not the system as a whole.

### 3.4. Time evolution

A self-consistent treatment of the orbital evolution of the BH with diffusion coefficients would require a full 2D Fokker-Planck integration. The configuration we are treating is far for equilibrium and will evolve on a dynamical timescale. For that reason, we make more progress with a set of restricted N-body integrations (see Sect. 4). Here we focus on the linear and diffusion terms of Eq. (A.7) to compute the increase in kinetic energy of each parallel and orthogonal degrees of freedom by treating their independent evolution. With this rough treatment we can estimate quantitatively the timescale over which the dynamical traction becomes dominant and the BH acquires significant angular momentum, $L_z$.

The specific kinetic energy $E_k$ of the BH is a sum over quadratic components of the velocity $v^2 = v_\parallel^2 + v_\perp^2$. We split $E_k = E_\parallel + E_\perp$ such that

$$E_{\parallel} = \frac{1}{2}v_{\parallel}^{2} + \frac{1}{2}\sigma_{\parallel}^{2}, \qquad (11)$$

$$E_{\perp} = \frac{1}{2}\sigma_{\perp}^{2}. \qquad (12)$$

We computed e.g. $E_\perp$ from Eq. (A.7) with a simple summation over time; likewise for the two components of $E_\parallel$. For instance,

$$\sigma_{\perp} = \sqrt{2E_{\perp}} = \left( 2 \int_{0}^{t} \langle (\Delta v_{\perp})^{2} \rangle dt \right)^{\frac{1}{2}}.$$

By assuming that the diffusion terms always have the same (positive) sign, we optimise the growth of kinetic energy orthogonally to the rotation flow. Since the flow of rotating stars singles out the torque acting on the BH, we argue that this approach slows down the transfer of angular momentum to the BH because the BH would take longer to develop a velocity vector that runs parallel, or nearly parallel, to the flow of stars.

To avoid a singularity in the diffusion coefficients, we set up the BH with a small velocity vector such that $v_\parallel = v_\perp \approx 3 \text{ km s}^{-1}$ initially at a cylindrical radius $R = 1.5 \text{ kpc}$; we also

<footnote n="2">² These results were obtained with $M_\bullet = 1 \times 10^7 \text{ M}_\odot$.</footnote>
<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Amplitude of the BH velocity components as a function of time. The figure graphs the velocity dispersion for each parallel and orthogonal component (denoted as $\sigma_{\parallel \text{or } \perp}$ in the legend) and the (first-order integrated) streaming velocity ($v_{\parallel}$, black solid line). The period of a circular orbit at radius $R = 1$ kpc (1.5 kpc) is $\approx 61$ Myr (95 Myr).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Critical line of stability for various configurations. The vertical axis is the BH mass $M_{\bullet}$ plotted against azimuthal velocity $v_{\phi}$. At fixed $M_{\bullet}$, the minimum azimuthal velocity $v_{\phi}$ to reach such that the dynamical traction sets in is found by crossing one of the black curves, from left, to right. Four curves are plotted, each corresponding to a different value of the stellar velocity dispersion $\sigma_{\star}$ as indicated in the legend. The blue dot marks the estimated coordinates of the Milky Way's BH. The red curve highlights an asymptotic trend of slope $\approx -2$ found in the limit of low $v_{\phi}$.</figcaption></figure>

set $\sigma_{\star} = 15 \text{ km s}^{-1}$ in Eq. (8) to boost the growth of velocity components (dispersion and stream; see Fig. 2). The results of the integrals over time are displayed on Fig. 3 for a BH of mass $M_{\bullet} = 1 \times 10^7 \text{ M}_{\odot}$, with a ratio $M_{\bullet}/m_{\star} = 10^3$. The initial rapid growth of the velocity dispersion components $\sigma_{\parallel \text{or } \perp}$ is clearly visible, but we note how it quickly levels off as the BH velocity increases (large relative velocities diminish the efficiency of energy diffusion). On the other hand, the component parallel to the stellar flow, i.e. $v_{\parallel} = v_{\phi}$, reaches a larger amplitude than either of the two dispersions $\sigma_{\parallel \text{or } \perp}$ after $\approx 270$ Myr of evolution. The square streaming velocity, $v_{\parallel}^2$, of the BH soon exceeds the sum of the two squared dispersion components $\sigma_{\perp}^2 + \sigma_{\parallel}^2$ at times $t \approx 300$ Myr (see Sect. 3.5 below). In other words, there is more specific kinetic energy in the streaming motion (carrying angular momentum) after 300 Myr of evolution. This is close to three full periods of $\approx 95$ Myr at a distance $R = 1.5$ kpc and with a rotation flow $R\Omega = 100 \text{ km s}^{-1}$. Not surprisingly, perhaps, this coincides exactly with the estimated time to dissipate kinetic energy in isotropic systems given by Eq. (2). But where one effect drains out orbital energy, the traction mode opposes this trend.

### 3.5. The onset of the transition to dynamical traction

As seen from Fig. 3, at some point in time the BH acquires a velocity vector such that its main component is $v_{\phi} \approx v_{\parallel}$ so that it moves almost (though not exactly) in the direction of the flow of stars. Its $z$-component angular momentum $L_z = Rv_{\phi}$ also increases with time. We stress that the time integral that we performed did not take into account any motion in $R$ or $z$. Despite this caveat, it is interesting to predict when the growth of angular momentum through the transition to dynamical traction must set in and transform the BH orbit into a loop orbit. We define this turning point as the value of $v_{\phi}$ such that the rate $v_{\parallel} \langle \Delta v_{\parallel} \rangle$ exceeds the sum of the rates due to the diffusion terms (see Fig. 2). We find that turning point for $M_{\bullet}$ in the bracket $[10^5, 10^9] \text{ M}_{\odot}$ and plot the results for four values of the stellar velocity dispersion $\sigma_{\star}$. We look at a few cases more closely below.

#### 3.5.1. Example for an $M_{\bullet} = 10^7 \text{ M}_{\odot}$ BH

The outcome for that case is shown on Fig. 4. Consider the case of a BH with $M_{\bullet} = 10^7 \text{ M}_{\odot}$. At low values of $v_{\phi}$, the BH gains and loses velocity components isotropically from the background stars. The gain and loss of such increments acquired stochastically mean that we do not expect a net increase in angular momentum, and the configuration is deemed stable. If and when the BH has acquired a significant azimuthal component, however, the situation changes. When $\sigma_{\star} \approx 70 \text{ km s}^{-1}$ and $v_{\phi} \ge 7 \text{ km s}^{-1}$, the positive transfer rate $v_{\parallel} \langle \Delta v_{\parallel} \rangle$ noticeably exceeds the other terms, and hence $v_{\phi}$ will continue to grow. Thus as the $M_{\bullet} = 10^7 \text{ M}_{\odot}$ BH crosses the dotted line in Fig. 4, the system switches from a stable to an unstable configuration. The physical reason for this is that the BH continues to induce a response in the stellar density as described originally by Chandrasekhar (1943), whereas the weak response coming from the quadratic diffusion coefficients only becomes weaker as the motion of the BH aligns with that of the flow of stars. (To first order one has $\langle \Delta v_{\parallel} \rangle = 0$ by symmetry, which is exact if we set $v_{\phi}$ strictly equal to $v_{\parallel}$.)

When the stellar velocity dispersion, $\sigma_{\star}$, is reduced, a BH triggers a response from a wider range of stellar orbits. Consequently the stochastic evolution of the BH's orbit persists until a larger azimuthal velocity component is reached and the transition to a systematic increase in $L_z$ takes place. As an example, when $\sigma_{\star}$ is reduced from 70 to 30 km s$^{-1}$, we find that the transition occurs at $v_{\phi} \approx 35 \text{ km s}^{-1}$ (solid line on Fig. 4). If we fix $\sigma_{\star} = 50 \text{ km s}^{-1}$, the transition point shifts back to $v_{\phi} \approx 15 \text{ km s}^{-1}$. A similar interpretation follows from increasing the mass $M_{\bullet}$ at constant $\sigma_{\star}$. A more massive BH will trigger a strong linear response from the flow of stars and will thus
acquire angular momentum steadily even if it is surrounded by stars of a relatively high velocity dispersion. In contrast, if the stellar velocity distribution function is very cold (low $\sigma_*$), a BH would trigger a strong response out to a large radius, in all directions. In that situation, it would not undergo a steady growth of its angular momentum unless it is given a large azimuthal velocity, for example, from an external perturbation. The shift over to a critical value of $v_\phi$ is sketched with an arrow on Fig. 4.

### 3.5.2. The case of the Milky Way BH

We have indicated with a large blue dot on Fig. 4 the position of the Milky Way’s BH at $M_\bullet = 4 \times 10^6$ M$_\odot$ and with a velocity in azimuth $\approx 2$ km s$^{-1}$, an upper limit for the velocity of Sgt A* in the Galactic frame of reference (Reid & Brunthaler 2004). The situation is complicated by the presence of a nuclear cluster of stars (NSC) around the BH (Neumayer et al. 2020 for a review). Chatzopoulos et al. (2015) estimate the half-light radius of the NSC at $\approx 7$ pc, where the velocity dispersion is on the order of 35 km s$^{-1}$. Their modelling of the NSC suggests significant rotation, with an isotropic rotator model (Binney & Tremaine 2008) of angular frequency$^3$ $\Omega_{nsc} \approx 30$ km s$^{-1}/3.9$ pc. Despite the high NSC angular velocity, the mean azimuthal velocity $v_\phi \approx 0.8 \sigma_*$ $\approx 28$ km s$^{-1}$ for the stars inside the half-light radius (Neumayer et al. 2020). The solid line on Fig. 4 gives a critical azimuthal velocity $v_\phi \approx 40$ km s$^{-1}$ for a stellar velocity dispersion $\sigma_* \approx 35$ km s$^{-1}$. This is significantly larger than the measured rotation for the stars. In conclusion, the BH would not suffer a dynamical traction transition even if it had acquired a $v_\phi$ of the same amplitude as that of the stars. The ratio $v_\phi/\sigma_* \approx 0.8$ implies that the NSC is nearly isotropic in velocity space, so the BH would sink back to the origin quickly (Chatzopoulos et al. 2015).

That line of reasoning follows from experience with fully isotropic velocity distribution functions, $f(E)$. For that case a massive body always loses momentum (Binney & Tremaine 2008). Consider the estimated mass of the NSC’s stellar population at $\approx 9 \times 10^6$ M$_\odot$, which exceeds that of the BH itself only by a factor $\times 2$. If we apply an $m = 1$ perturbation to the BH, its motion would trigger a strong response from the NSC stars and bring back the BH to the barycentre on a short mass-segregation timescale given by

$$t_{\text{ms}} \propto t_r \approx \frac{m_*}{M_*} \times \frac{N_* t_{\text{dyn}}}{\ln N_*} \sim \frac{2}{\ln N_*} t_{\text{dyn}}, \quad (13)$$

with $t_r$ denoting the two-body relaxation time and $M_* \approx N_* m_*/2$. The NSC dynamical time is

$$t_{\text{dyn}} \sim \frac{r_h}{\sqrt{v_*^2}} \approx 0.2 \, \text{Myr}$$

with a half-light radius, $r_h$, of $\approx 7$ pc and a root-mean-square velocity of $\approx 30$ km s$^{-1}$. (See e.g., Merritt (2013), Šubr & Haas (2014) or Panamarev & Kocsis (2022) for further details and discussion on the collisional dynamics of the NSC.)

On the other hand, the impact of the nuclear stellar disc (NSD) on a scale of $\sim 150$ pc of an estimated mass $\sim 10^9$ M$_\odot$ is much more likely to trigger a sustained transition to dynamical traction (Launhardt et al. 2002; Fritz et al. 2016;

Nogueras-Lara et al. 2020; Nieuwmuister et al. 2024; (but see Zoccali et al. 2024) for a critical assessment). The dynamical timescale $t_{dyn} \sim 1/\sqrt{G\rho}$ equals approximately 1 Myr. With the rotation velocity of $\approx 50$ km s$^{-1}$ at radius $r = 100$ pc from the Gaia survey (Gaia Collaboration 2023), we computed a period of $\approx 12$ Myr for a circular orbit at that radius. The growth rate of azimuthal velocity displayed in Fig. 2 over several $\times 10$ Myr would be significant. We note that the NSD has been the site of recent star formation events over times of 30 to 100 Myr (Nogueras-Lara et al. 2020), while some low-metallicity stars may be as young as 100 Myr within the central-most NSC (scale of 10 pc; see von Fellenberg et al. 2022; Chen et al. 2023). These are indirect indications of past perturbation on the right timescale for the full growth of the transition to dynamical traction. We do not attempt to model the central region of the Milky Way in this contribution. Instead we present a first numerical exploration of the onset of the transition in Sects. 4 and 5.

## 4. Restricted N-body calculations

Recently, Kaur & Sridhar (2018) and Banik & van den Bosch (2022) studied the exchange of mechanical energy between stars and an orbiting BH analytically and with a set of numerical integrations (see also Boily et al. 2008)). These studies show how the response of stars either drains energy from or invests energy in the BH’s orbit depending on the character of the stellar orbits and their parameters such as their phase, inclination angle, and so on. We view these analyses as a strong hint that different families of stellar orbits couple to a different degree with an BH, and in particular that anisotropy in the stellar velocity field is a key ingredient to the BH’s orbital evolution. Here we borrow from Boily et al. (2008), who noted that co-planar stellar orbits couple much more effectively to the BH than those at high inclination angle (see their Sect. 5.3 and Fig. 14b). This allows for a simplification of the numerical treatment, which we lay out below.

### 4.1. Initial conditions

We first divided the stellar population into spherical and disc components, whereby the dynamically cooler disc component is made up of stellar orbits with an angular momentum orthogonal to the orbital plane of the BH. We froze the other (dynamically hot) stellar orbits completely. The system to evolve is made of two potentials. For the frozen potential, we picked the spherical isochrone model of Hénon:

$$\Phi_{\text{Is}}(r) = - \frac{GM_{\text{Is}}}{a + \sqrt{a^2 + r^2}}, \quad (14)$$

where $a$ is a radial scale length and $M_{\text{Is}}$ is the total mass. For the live potential, we opted for an axi-symmetric Miyamoto-Nagai disc in cylindrical coordinates, $R, \theta, z$ (see Sect. 2 of Binney & Tremaine 2008),

$$\Phi_D(R, z) = - \frac{GM_D}{\sqrt{(a + \mathcal{Z})^2 + R^2}}, \quad (15)$$

where we defined $\mathcal{Z}^2 = z^2 + b^2$; $M_D$ is the total disc mass; and $a, b$ are radial and vertical scale lengths, respectively. We opted to fix $a = 1.5$ kpc in both potentials and the ratio $b/a = 1/6$ for the disc potential.

These potentials were combined, and their joint stellar phase-space distribution function was solved in action space with the

<footnote>$^3$ Notice how $\Omega_{nsc} \approx 10^4$ km s$^{-1}$/kpc $\approx 20$ times the angular rotation frequency at galactic scales: Gaia data suggests $\Omega \approx 50$ km s$^{-1}$/kpc for the innermost 100 pc of the Milky Way (Gaia Collaboration 2023; see also Argüelles & Collazo 2023).</footnote>
Agama package (Vasiliev 2019). A restricted N-body realisation of the stellar disc was obtained from a Monte Carlo sampling of the joint distribution function. We noted how the response of centrophobic loop orbits subjected to the motion of a massive perturber is the strongest when they are co-planar. The root-mean square variation in mechanical energy for loop orbits reaches ≈80% of the peak co-planar value whenever their orbital angular momentum remains within ±≈18° of the polar vector of the perturber’s orbital plane (see Fig. 14b of Boily et al. 2008, it becomes progressively weaker for larger inclination angles). The solid angle spawned around the pole is close to 5% of the whole sphere. An isotropic distribution function would yield orbits with angular momentum vectors covering the full sphere uniformly. The disc potential (Eq. (15)) of aspect ratio b/a = 1/6 gives rise to near-circular motion within a wedge of opening angle = arctan 1/6 ≈ 9.5° for close to 85% of the orbits (the distribution function becomes progressively more isotropic as we approach the centre, within a cylindrical radius R < a; it is very nearly isotropic inside a small region bound by R ≈ b = 250 pc). This sits comfortably within the solid angle for which we have estimated the orbits will respond the strongest to the perturbation. With these parameters we therefore focus on discy (nearly co-planar) orbits. We could have extracted co-planar orbits directly from the spherical isochrone distribution. Instead, we recall our motivation for more anisotropic systems in the context of a galaxy merger, or the end-state of the formation of a disc galaxy. For that reason, we chose to boost the sample of disc orbits by fixing the mass parameters in Eqs. (14) and (15) in the ratio

$$ \frac{M_D}{M_{Is}} = \frac{1}{6} \approx 0.167. $$

This increased the fraction of the disc orbits by about three times what we would have drawn from an isotropic d.f. in spherical symmetry.

A physical mass of $M_{Is} = 1 \times 10^{10} M_{\odot}$ was chosen to match the mass of the central kiloparsec of a typical galaxy. The central atomic hydrogen density of the isochrone sphere $\rho_{Is} \approx 7.16 \text{ cm}^{-3}$ is only slightly lower than that for the disc, for which $\rho_D \approx 8.6 \text{ cm}^{-3}$. These values are roughly twice the estimated value of $3.9 \text{ cm}^{-3}$ for the Milky Way at the solar radius (Bland-Hawthorn & Gerhard 2016).

## 4.2. Numerical integration

We performed numerical experiments using the Amuse suite of community codes to integrate the equations of motion (Pelupessy et al. 2013; Portegies Zwart & McMillan 2018). The python package Amuse v.5.2 (November 2023) offers an interface to treat each component of a system as a set endowed with its own methods and attributes. We opted to couple a live set of stars to a subset of BHs and test stars using the bridge methodology (Fujii et al. 2007; Portegies Zwart et al. 2013). A bridge easily allows one to couple a Barnes-Hut Tree integrator (Barnes & Hut 1986) and the direct-summation 4th-order accurate Hermite integrator PH4 (Portegies Zwart & McMillan 2018) to integrate the disc stars and BH orbits, respectively. Both integrators are ‘bridged’ to a frozen background potential that stands in for the extended galactic centre. The BHTREE code has a setup with critical opening angle $\Theta_c$ which we sat to 0.5 (down from the default 0.75), and a critical number of near-neighbours for multipolar expansion which we boosted, from the default 12, to 16 and in some cases, 32. This tactic meant that more stars were integrated with a 2nd-order accurate leapfrog symplectic scheme, whereas distant stars were integrated

with a quadrupolar expansion of the stars’ potential. A large softening length $l = 125$ pc was used for the two-body star-star interaction. We found this to be the largest value that still gives a sensible resolution of the vertical structure of the disc (scale height $b = 250$ pc) while also allowing any Jeans-type mode of instability to be resolved. The stars-BH interactions were treated with a much smaller softening kernel where $l = 10$ pc. The more compact kernel allowed for larger deflection angles of the stars and a accurate rendering of the BH potential (the effective size of the stellar mass elements ~40 pc, see below). The parameters used for the models and integrators are listed in Table 1.

Calculations were done by sampling the disc d.f. with $N = 100000$ stars, so each star has a mass $m_{\star} = 0.16 \times 10^3 M_{\odot}$. The BH has a reference mass $M_{\odot} = 1.25 \times 10^7 M_{\odot}$ for a mass ratio $M_{\odot}/m_{\star} = 781$. This takes us clearly in the galactic dynamics domain with each ‘star’ the mass equivalent of a rich open cluster of (real) stars. The disc density $\rho_D \approx 8.6 H/cc$ implies a mass inside the volume spawned by $l = 125$ pc of $1.67 \times 10^6 M_{\odot}$, or about 100 of our model stars. We argue this makes all galactic sub-structures well resolved on a scale ~100 pc and above. The physical mass of the BH exceeds the stellar mass (by a factor ≈7.5) in a spherical volume spawned by $l$; hence, the influence of the BH is correctly resolved on scales exceeding ~100 pc. The dynamics of the stars is also well recovered on such scales, because their individual mass is recovered from a spherical volume of radius ≈ 45 pc. Finally, we designed a global units converter to ensure that the Amuse integrators and the Agama module work hand in hand.

A series of test cases are presented in Appendix C. We checked that the parametric setup of the numerical integrators allow circular orbits to be recovered in the frozen potentials to a high accuracy. For a live disc component, we note that the total system angular momentum is conserved to about 0.4% accuracy for a runtime of $t = 1.5$ Gyr, which is sufficient to monitor the transfer of the orbital angular momentum from the stars to the BH, with three significant digits. In Appendix D, we address the question of the stability of the configurations with a frozen background component.

# 5. Results

We proceed with the problem of an in-coming BH as outlined in Sect. 3. To keep the initial conditions close to the reference configuration of Sect. 4.2, we use the same disc/halo configuration but shift the BH to a radius $a = 1.5$ kpc, at the edge of the core of approximately constant density (recall that the same radial scale length was used for the disc and isochrone halo). No angular momentum is given to the BH initially. The initial configuration is similar to what is depicted on Fig. 1(b), after an $m = 1$ mode dislodged the BH from the origin. The analysis leading to Fig. 3 did not factor in the orbital evolution (especially the radial distance, $R$) which must occur on a timescale ~100 Myr for the chosen radial orbit. This zero-momentum BH orbit is an idealised case of a scenario where the BH is being accreted by the system.

## 5.1. The case of a warm disc

When they form in a cold disc, density waves can carry angular momentum in the same way as the gravitational focusing leading to dynamical friction. To minimise the impact of such waves, we discuss first the case of a warm disc (Table 1) which remains
<table><caption>**Table 1.** Default parameters of the potentials and Amuse numerical integrators.</caption><thead><tr><th colspan="4">Isochrone (Hénon) potential</th><th colspan="4">Miyamoto-Nagai disc</th></tr><tr><th>M<sub>ls</sub></th><th>a</th><th>σ<sub>r</sub></th><th>Comment</th><th>M<sub>D</sub></th><th>a</th><th>b</th><th>{σ<sub>R</sub>, σ<sub>z</sub>}</th><th>Comment</th></tr></thead><tbody><tr><td>M<sub>☉</sub><br>1 × 10<sup>10</sup></td><td>kpc<br>3/2</td><td>km s<sup>-1</sup><br>22.7</td><td>frozen, infinite extent</td><td>M<sub>☉</sub><br>1.6 × 10<sup>9</sup></td><td>kpc<br>3/2</td><td>kpc<br>1/4</td><td>km s<sup>-1</sup><br>{4.2, 4.6}</td><td>warm:<br>σ<sub>R</sub> or z ≈ 12 km s<sup>-1</sup></td></tr><tr><th colspan="4">Barnes-Hut integrator</th><th colspan="4">Hermite (direct) integrator</th></tr><tr><th>Name</th><th>N<sub>2b</sub></th><th>θ<sub>c</sub></th><th>l<sub>★</sub> [pc]</th><th>Name</th><th>N<sub>2b</sub></th><th>l<sub>★</sub> [pc]</th><th>Mode</th><th>Accuracy</th></tr><tr><td>BHTREE</td><td>16 or 32</td><td>0.50</td><td>125</td><td>PH4</td><td>2</td><td>10</td><td>openmp</td><td>4th order</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Case when a BH is let go from $R = 1.5$ kpc in a (dynamically) warm stellar disc (Table 1). In panel (a), the run of $L_z$ is shown as a black solid curve, and radial and azimuthal velocity components are shown in red. In panel (b), the cylindrical radius $R$ and the amplitude of $v_★$ are displayed, each with their own colour-coded vertical axis. In the early stage, dynamical friction makes the orbit shrink; however, the perturber BH acquires angular momentum throughout (panel [a], black curve). Both $R$ and $v_★$ increase systematically from $t \approx 500$ Myr onward. The dotted line at the bottom is the $z$-coordinate bound to $\pm 15$ pc. Compare with the test case displayed on Fig. C.2.</figcaption></figure>

axially symmetric throughout the integration time (barring the response to the BH's motion).

### 5.1.1. Radial orbit

For a BH set on a radial orbit in a warm disc, the orbital evolution develops significant dynamical friction early on as the BH falls from rest to the system's barycentre. The run over time is graphed on Fig. 5 for the radius, velocity components, and $L_z$. One can clearly see how both the radial velocity $v_r$ and cylindrical radius $R$ diminish during the first ~500 Myr of evolution. This corresponds to about 4 periods of the circular orbit at the initial radius, in good agreement with the timescale given by Eq. (2). But where the dissipative friction drops in amplitude with the velocity $v_★$ (the Stokes limit), the dynamical traction exerted by the stream of stars remains steady. In fact, on Fig. 5(a), the black curve mapping the angular momentum $L_z$ increases constantly, as does the azimuthal velocity component (solid curve in red). If the effect of dynamical friction is to kill off the radial velocity component $v_r$, one sees that dynamical traction, on the contrary, boosts the orthogonal $v_φ$ component steadily. Hence, the BH never reaches the centre, and in fact it begins to migrate outward as it gains angular momentum (Fig. 5[a] and [b]).

### 5.1.2. Circular orbit

To determine whether the BH's orbital evolution is correctly recovered from the Fokker-Planck treatment of Sect. A, we restarted the numerical integration but now with the BH initially set on a circular orbit (same stellar disc and setup otherwise). The expectation from Eq. (A.3) and Fig. 2 is that the BH should remain trapped at constant radius, since there is in principle no net gain or loss of kinetic energy extracted from locally co-moving stars.

The results are displayed on Fig. 6 with the same colour-coding and line styles as on Fig. 5. The expected circular orbit shows significant evolution over time, for several reasons. Firstly, large deflection-angle stars-BH encounters (collisions) lead to small but noticeable time-evolution of the BH's radius and angular momentum during the first ~200 Myr. This is close to two orbital periods. Collisions take place because the epicyclic motion stemming from the stellar velocity ellipsoid brings stars on and out of reach of the BH radius of influence, and this happens over a few revolutions about the centre. Secondly, the (global) reflex response of the stars to the BH's perturbation also contributes to destabilise the BH's orbit on the same timescale<fnref n="4" />.

<footnote n="4">⁴ It would have been simple to set up a quiet start to nullify that effect, for example by growing a BH adiabatically, but that would have been contrary to the context of an accretion event viewed here as an external perturbation.</footnote>
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Similar to Fig. 5 but now the BH started at $R = 1.5$ kpc on a circular orbit. We note the sharp response of the axi-symmetric disc to the BH, which triggers radial motion. The disc adjusts to the black hole perturbation and settles progressively over several orbital times to a new equilibrium. The BH's orbit, however, does not migrate significantly over the full integration time. Panel (a): Run of $L_z$ (black solid curve). Radial and azimuthal velocity components are shown in red. Panel (b): Cylindrical radius $R$ (in black) and amplitude of $v_\bullet$ (in turquoise). The dotted line at the bottom is the z-coordinate bound to $\pm 17$ pc.</figcaption></figure>

Further orbital evolution takes place over a longer timescale of ~500 Myr as the BH radius varies with time in the bracket 1.3–1.7 kpc (Fig. 6[b]). We have not explored the reasons for this behaviour in any detail. However, we note that small radial displacements will be amplified if and when dynamical traction (or, friction) takes place since it increases (or decreases) the orbital angular momentum and thus boosts the orbit’s eccentricity. Also worth noting is that the timescale of 500 Myr compares well with the time of ~300 Myr needed to transition to an azimuthally dominant increase of kinetic energy (see Fig. 2). Our conclusion is that exact balance between dynamical friction and traction is not enough to maintain a BH on a circular orbit due to the response of the stars, first locally through strong scattering and then globally by the relaxation to a new equilibrium configuration of the system as a whole. Despite these caveats, we find no systematic and continued radial migration for a BH starting with a high-$L_z$ orbit, in contrast with the case of initially having only radial motion.

The results for these example orbits give confidence that the basic mechanics of friction/traction and transfer of angular momentum is well recovered by the numerical treatment. We also find that the amplitude of the motion and time-evolution stay well above the noise level graphed on Fig. C.2. The higher velocity dispersion of disc stars, of total amplitude $\approx 25$ km s$^{-1}$, exceeds the circular velocity in the BH’s potential at the linear resolution $l = 125$ pc, for which we computed

$$v_{\text{circ}} = \sqrt{\frac{GM_{\bullet}}{l}} \approx 18 \text{ km s}^{-1}.$$

The relative velocity between stars and the BH orbiting at the same radius (so sharing the same velocity in azimuth) may yet be high enough to prevent a significant number of stars from binding with the BH. Visual inspection of density maps suggests that only a few stars bind up with the BH. That situation will change when the BH orbits in a dynamically cooler disc, as we argue in the next sub-section.

## 5.2. The case of a dynamically cool disc

We now turn to the more complex case of a cool disc, when the vertical stellar velocity dispersion is reduced to $\sigma_z \approx 4.5$ km s$^{-1}$;

the disc effective scale height $\approx 100$ pc is lower by more than a factor of two compared with the warm disc (when $h \approx 240$ pc). Self-gravitating Jeans fragments and large-scale spirals are triggered by the BH’s gravity, which we view as an $m = 1$ perturbation mode. In the remainder of this section, the discussion follows the reference configuration for a dynamically cooler disc listed in Table 1 unless stated otherwise.

### 5.2.1. Radial BH orbit

Let the BH start from rest at $R = 1.5$ kpc initially, as done before for the warm disc (Sect. 5.1.1). The early evolution is staged on Figs. 7 and 8, where we find the BH initially falling to the system’s barycentre due to the mean gravitational field. A full orbit is completed in $\approx 123$ Myr. These first few orbital times show the BH triggering large-amplitude density waves, and these last for several hundred Myrs, seen both in density maps and in the velocity field (see Fig. 8). In the time interval running up to $t \approx 400$ Myrs, the BH acquires significant angular momentum, as shown on Fig. 9(a), where we plot the z-component of $L = M_{\bullet} R_{\bullet} \times v_{\bullet}$ along with the radial and azimuthal components of $v_{\bullet}$. At $t = 200$ Myr, we compute $L_z \approx 120 \times 10^6 M_{\odot} \text{ kpc/km s}^{-1}$ or $\approx 20 \times$ the peak value obtained in the test calculation with initially $v_{\bullet} = 0$ (Figs. C.1 and C.2). Therefore this early evolution is not the result of numerical inaccuracies or approximations made. Furthermore, a shift in the orbital phase of the BH has no effect (Fig. 10).

Despite the strong coupling and efficient transfer of angular momentum between stars and BH, dynamical friction soon takes over and drives the rapid decrease in radius, bringing the BH systematically closer to the centre (Fig. 9[b]). Specifically on Fig. 8, bottom frame when $t = 452$ Myr, we find the BH within $r_{\bullet} \approx 52$ pc of the centre, well below our spatial resolution limit of 125 pc. At that time, it is moving at a velocity of $\approx 11.3$ km s$^{-1}$ in the direction of the large banana-shaped clump directly above it to the north. The phase-space coordinates of the BH are consistent with harmonic motion in the frozen background potential as depicted on Fig. C.1. However, the kiloparsec-size over-density breaks the axi-symmetry of the background potential and exerts a net pull on the BH. We wondered what impact that would have on its orbital evolution.

To investigate this, we first draw an outline of the over-dense structure with a surface density threshold of $>0.06$, or four stan-
<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Evolution of the system when the BH is released from rest at a radius of $r = 1.5$ kpc (from top left to bottom right). The black arrow indicates the velocity vector. Large density fluctuations soon develop in the inner ~1 kpc where the BH potential becomes dominant. After $t \approx 100$ Myr of evolution, strong density fluctuations appear, including large-scale voids (deep blue). The BH used in the calculation had a mass $M_{\bullet} = 1.25 \times 10^7$ M$_{\odot}$.</figcaption></figure>

dard deviations above the mean<fnref n="5" />. We display the selection with three density contour levels on Fig. 8. Stars that fall inside that area were picked up by estimating the density at the position of each star using the same Gaussian kernel as for the gridded maps of the disc. The set of stars selected together make up a mass of $32 \times 10^6$ M$_{\odot}$ (2002 member stars), or close to $2.6 \times$ the mass of the BH; we refer to this set of stars as s452. To inspect the internal cohesion of s452, we computed its total mechanical energy $E = E_k + W$ and found $E \approx 8046.6$ (km/s)$^2$ M$_{\odot} > 0$ with a ratio $E_k/|W| \gtrsim 1.77$ to 1.88, depending on whether the was trimmed of a subset of the most distant stars to its barycentre, or not. The positivity of the mechanical energy is robust against details of the selection procedure, and implies that the clump is not bound. An unbound over-density should be interpreted as a wake caused by dynamical friction triggered by the massive perturber (Chandrasekhar 1943; Binney & Tremaine 2008). To assess the likelihood that a subset of s459's member stars are not,

by themselves, a bound sub-system, we re-computed the energy $E$ of each star in the rest frame of the selected clumped stars taken in isolation. We identified stars that had a negative energy, and found out that these make up 62% of the 2002 stars selected. We then removed all the stars with a positive energy and recomputed the potential of the new stellar clump with (now) only 843 stars: that system had a total positive energy, with a kinetic- to gravitational energy ratio very close to $E_k/|W| = 1.40$. The conclusion we draw from this analysis is that self-gravity plays only a minor role in the interaction of the clump with the BH, because unbound stars contribute a too-important fraction of the clump's gravitational potential.

s452 member stars disperse progressively over a timescale of ~200 Myr, close to two orbital periods, somewhat slowed down by their own gravity. We ask whether the initial configuration of s452 may efficiently drag the BH out of the origin. Switching to the rest-frame of the BH, we find the clump moving away on a north, north-west trajectory at a velocity of $\approx 9.8$ km s$^{-1}$ (the velocity vector is displayed on Fig. 8). This configuration matches very much the $m = 1$ mode of perturbation discussed

<footnote n="5">⁵ In computational units, the surface density $\Sigma$ falls in the range [2.68·10⁻⁴, 9.01·10⁻²], with a mean value 1.25(6)·10⁻² and standard deviation $\sigma_{\Sigma} = 1.12(4) \cdot 10^{-2}$.</footnote>
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Transition to a fragmented morphology for a galactic disc perturbed by a massive BH. Left-hand panels: Density map of stars with peaks indicated in red. Right-hand panels: Velocity field. A black dot indicates the BH of $1.25 \times 10^7 M_{\odot}$ (or, ≈1% the mass of the stars), with its velocity vector also in black. The configurations are displayed at three different times: 144 (top), 308 (middle), and 452 Myr (bottom row) of evolution. Contour levels in yellow (left-hand panel) and grey (right-hand panel) identify the member stars of the clump formed at $t = 452$ Myr. The BH sits at ≈52 pc of the coordinates centre. The relative velocity between the clump's centre of mass motion and the BH is indicated with an arrow on both panels, pointing in a north, north-west direction (west is to the left).</figcaption></figure>

in Sect. 3, here on a radial scale of ≈1 kpc. The distance of the centre of mass of the clump to the BH is ≈340 pc, which gives rise to a monopole acceleration of $GM_{\text{clump}}/||\mathbf{R}_{\bullet} - \mathbf{R}_{\text{clump}}||^2 \approx 0.00131 \text{ kpc/Myr}^2$. An acceleration of this amplitude maintained for $\Delta t \approx 40$ Myr (the dynamical timescale) would displace the BH by

$$ \Delta x \approx \frac{GM_{\text{clump}}}{||\mathbf{R}_{\bullet} - \mathbf{R}_{\text{clump}}||^2} \frac{\Delta t^2}{2} \approx 1 \text{ kpc.} $$

That this expectation is not fully realised is seen on Figs. 9(b), where the BH orbital radius increases to a much lesser extent over ≈50 Myr of evolution. The steady, nearly linear rise of the
<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Case of a BH falling from rest on a radial orbit. a) The solid black curve graphs the angular momentum accrued over time (left-hand axis with black labels); the right-hand axis is the scale for each velocity component (red curves). b) The cylindrical radius (in black, left-hand axis) and velocity (right-hand axis, in turquoise) are shown as a function of time. The dotted line at the bottom is the z coordinate bound to ±12 pc.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Similar to Fig. 9 but for when the BH is launched from the origin at $v_y = 80 \text{ km s}^{-1}$ and comes to rest at a radius of $R \approx 1.35 \text{ kpc}$. The run of all quantities are remarkably similar, but for a shift in time of ≈200 Myr (≈ two orbital times) when $L_z$ peaks. We also note that the phase of dynamical friction does not bring the BH close to the origin, as its apocentre distance always exceeds 250 pc. The gain in angular momentum seen in panel (a; black curve) begins before 500 Myr.</figcaption></figure>

BH's radial distance and velocity in the interval $t = 452$ Myr to ≈ 1 Gyr seen on Fig. 9(b) instead suggests that the acceleration field applied to the BH does not remain steady, but fluctuates in amplitude over time<fnref n="6" />. In other words, the clump of stars exerting a pull on the BH certainly does not preserve enough of its mass and morphology to have a significant impact on its long-term orbital evolution. Because the size of the clump compares well with the Jeans length $\lambda_J \sim 670$ pc of Eq. (D.1), we initially concluded, wrongly, that the self-gravity of the clump enhanced significantly the amplitude and duration of the traction on the BH. The dissolution of s452 argues against that.

### 5.2.2. The core region, fragments

The BH perturbs the core region significantly as soon as its orbital radius ≈ 200 pc or lower, since inside that volume the

mass $M_\bullet$ compares with the total gravitational mass of stellar and frozen components. We found the velocity dispersion near the centre increasing over time to more than $20 \text{ km s}^{-1}$ (up to a peak value of $25.5 \text{ km s}^{-1}$ around the time $t \approx 452$ Myr), distorting both the radial flow and rotation curve. A high velocity dispersion facilitates the coupling between the BH and the streaming stars (e.g. Fig. 2), because it reduces the amplitude of Brownian motion by the perturber. After the BH has reached deep in the core region, the axi-symmetry of the stellar flow is broken: we argued in Sect. 5.2.1 that a single impulse by a massive clump of stars gives enough linear momentum to shift the BH away from the centre in a region where the stellar stream becomes dominant again. The angular momentum is transferred to the BH progressively thereafter, driven by dynamical traction from the streaming stars. It is worth noting that the BH closely follows the mean disc rotation from the time $t = 750$ Myr onward. The BH has essentially the same specific angular momentum as the surrounding stars. After a closer look at Fig. 9(b), we noted that it shows an orbital radius oscillating with significant amplitude, and the oscillation increases as time passes. The BH is on an eccentric orbit, and the implication is that it rotates faster and then slower

<footnote n="6">⁶ We checked that the bend in the angular momentum $L_z = \mathbf{R} \times \mathbf{v} \propto R v_\phi$ grows more slowly than $t^{5/2}$ in the interval $t$: [450, 1500] Myr, consistent with a quasi-linear relation $\propto t^\gamma$, $\gamma < 1.25$ for both $v_\phi$ and $R$.</footnote>
(in angular speed) than the stars as it shifts from the perigee to the apogee of its galactic orbit. It is therefore trapped around the radius of co-rotation with the surrounding stars (see Fig. 2 and Appendix A).

### 5.3. Dynamically cold discs and BH orbits with angular momentum

We turn to a more likely situation where the BH starts off on an orbit that carries angular momentum. The case of the circular orbit is the most illustrative. The radial and azimuthal velocity dispersion of the stars mean that the BH falls victim of the asymmetric drift (Binney & Tremaine 2008, Sect. 4.4, p. 326), implying a faster BH velocity in azimuth than the mean for the surrounding stars. The faster BH motion causes a net negative drag (from the stars) and a loss of angular momentum. In practice the stellar velocity dispersion is so low that the asymmetric drift causes a difference of only ~1 km s⁻¹ with circular motion, much lower still than the 1D dispersion of ≈4 km s⁻¹ (see Tab. 1). We expect the asymmetric drift to have a negligible impact here.

On the other hand, a reduced stellar velocity dispersion boosts the local gravitational focusing of the stars by the BH. This phenomenon enhances the number of stars that bind with the BH. At the spatial resolution $r = l = 125$ pc, the circular velocity in the BH's gravitational field,

$$v_c = \left(\frac{GM\bullet}{l}\right)^{\frac{1}{2}} \approx 18 \text{ km s}^{-1},$$

is much higher than the 3D velocity dispersion (i.e. ≈7 km s⁻¹) of the unperturbed stellar disc component. A significant number of stars in the annulus $R = 1.5 \pm 0.125$ kpc will therefore interact strongly with the BH (the annulus covers close to 2.8 % of the disc's surface but contains about 5% of the disc's mass; the BH makes up ≈0.8% of the disc's mass).

#### 5.3.1. Evolution from a circular orbit

The time evolution is mapped out on Fig. 11. It is noticeable on the frames running from $t = 40$ Myr to 308 Myr that the BH attracts stars and "cleans up" its orbit so that a trail of low-density (a gap) forms. The bridges of higher density that form on either side of its orbit lead to the efficient angular momentum transport and orbital evolution, as in proto-planetary circumstellar disc dynamics (so-called Type-II migration; e.g., Kley & Nelson 2012; Papalozou 2021). An importance difference is that no ram pressure exists in the disc material, and the spiral features are never close to axial symmetry. Despite these caveats, it is clear that gravitational torques acting on the BH's orbit, rather than the net effect of dynamical traction / friction, drive this rapid inward migration: a wake of stars slowing down the BH would show up as a trailing over-density on Fig. 11, in the frames $t = 61$ Myr to 308 Myr. Instead we find an area devoid of stars in front and behind the BH. Compare this situation with the case of the radial orbit (Fig. 7; see also Fig. 8): clearly for that case, an over-dense wake trails the motion of the BH. The rapid evolution of the disc response over a few periods causes the clumps of disc material to interact and merge with each other, and they eventually form dense substructures with a roundish shape. This situation is depicted in frame $t = 915$ Myr when the BH is within 20 pc of the centre. The BH is virtually at rest at that time, since its

velocity $v\bullet \approx 9$ km s⁻¹ lies within the numerical uncertainty for a BH starting at the centre of coordinates with no motion (Fig. C.2).

It is important to note that the fragmented mass density at $t = 915$ Myr breaks the axi-symmetry of the system. A direct consequence of this is that the BH does not remain still, but is pulled out of the origin by one large clump, or a few small ones. The energy and angular momentum gained by the BH means that by $t = 1553$ Myr, it is now orbiting at a radius $R \approx 300$ pc; it has acquired significant angular momentum which persists until at least $t = 1750$ Myr, when the disc material defines a bar-like potential on a scale of ~500 pc (the BH radial orbit is then at a radius ≈710 pc). The run of radius, velocity components and angular momentum is plotted on Fig. 12. It is striking how the angular momentum drops to essentially zero at $t \approx 915$ Myr, and rises systematically thereafter. The same comment applies to its orbital radius, seen on Fig. 12(b). The fact that the BH's orbit does not circularise exactly (small but persistent oscillations in radius $R$) may be one reason why symmetry breaking is enhanced near the origin, i.e. on a scale matching the BH's radius of influence, or $r \approx 300$ pc.

#### 5.3.2. Eccentric orbits

To underscore the singular character of the (initially) circular BH motion, we ran two more cases with the BH starting on mildly eccentric orbits. On Fig. 13, we display the runs of orbital angular momentum $L_z$ (panel [a]) and radius (panel [b]) for the initially circular orbit discussed above ($e = 0$), along with two orbits that begin with lower $L_z$. The first one has $L_z \rightarrow 5L_z/6$ leading to an orbit of small eccentricity ($e \approx 0.55$); the second one has $L_z$ reduced by a second factor $5/6$ so $L_z \rightarrow 0.7L_z$ compared to circular motion, and a larger orbital eccentricity ($e \approx 0.67$). We estimated the eccentricities from the ratio of radii at peri- and apo-centres.

To understand the trends in $L_z$ and radii over time, it is useful to think in terms of epicyclic motion. In this regime, the angular momentum is conserved as the object moves in and out radially. We may assume that this is true of the unperturbed stellar orbits. Here, the BH may dominate the gravitational potential locally and polarises stellar orbits of lower and higher $L_z$, inside and outside of its orbit, respectively. When the BH starts off on a nearly circular orbit, its own epicycle does not generate significant radial motion compared with the radial stellar velocity dispersion $\sigma_R \approx 4.3$ km s⁻¹. Thus in the limit of small $e$, the BH act as a gravitational focal point. It is able to open a gap along its orbit, when it becomes a catalyst for angular momentum transfer, from the stars inside its orbit, to the stars outside it. This is a very similar situation to the instability studied by (Julian & Toomre 1966) for a thin, cold fluid of stars. On the contrary, when the BH is launched from a low-angular momentum orbit, of higher eccentricity $e$, it is unable to clear a gap in the stellar population because its own significant radial motion effectively boosts its velocity relative to the stars. This in turn makes the BH a less effective gravitational attractor. The BH orbit labelled large $e \approx 0.67$ on Fig. 13(b) gives a mean-amplitude (rms) radial velocity of $v_R \approx 6.8$ km s⁻¹ which is in fact larger than the stellar radial velocity dispersion, at least initially. Hence the BH stops being a point of gravitational focusing, and fewer stars respond strongly to its gravitational pull. A gap does not form, and the BH orbit preserves its angular momentum more efficiently.
<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Evolution of the BH set on a circular orbit initially of radius $r = 1.5$ kpc (top left to bottom right). The black arrow indicates the velocity vector. Large spiral waves and Jeans fragments develop quickly due to the lower stellar velocity dispersion. The BH binds with stars and opens a gap in its wake, which leads to a steady loss of angular momentum and inward migration. This is most visible in the frames with time running from $t = 41$ Myr to 915 Myr. At $t = 915$ Myr the BH sits at $r \approx 20$ pc with a rest-frame velocity $\approx 9$ km s$^{-1}$. Non-axi-symmetric clumps of stars dislodge the BH which acquires significant angular momentum by dynamical traction in the later stages of evolution (on a timescale $\sim 500$ Myr). Calculations performed with a BH of mass $M_{\bullet} = 1.25 \times 10^{7}$ M$_{\odot}$.</figcaption></figure>

## 5.4. Transition from gravitational torques to dynamical friction

Our study of BH orbital evolution points to a more complex picture than for warm discs. The time evolution of a BH orbit initially with zero angular momentum is first driven by dynamical friction in competition with dynamical traction (see Sect. 5.1.1). This happens both in a dynamically warm- or cold stellar disc, and can lead to a transition, from a contracting orbit losing angular momentum, to an expanding orbit gaining angular momentum. When the disc is cold and the BH starts on a near-circular orbit, the mechanics at play is not that of dynamical friction because a massive body will suck in stars that are close in phase-space: the formation of a gap of under-dense regions stretching in front and behind the massive body is the signature of this process (Fig. 11). We have argued, based on two cases of BH orbits that start on mildly eccentric orbits, that gap formation shuts off if the BH is on a sufficiently eccentric orbit, and the interplay of dynamical friction / traction kicks in. This transition to an evolution of the orbital momentum driven by gravitational focusing

(friction or traction) occurs at a critical orbital angular moment for fixed values of $M_{\bullet}$ and stellar velocity dispersion $\sigma_{\star}$ (the radial component mostly controlling the Toomre stability parameter $Q$; see Eq. (D.2)). We have not identified a relation between $Q, L_z$, and $M_{\bullet}$ that quantifies this transition. We anticipate that this point is worthy of a study all by itself, to which we hope to return in a future contribution. In the next section, we seek to justify in more quantitative details the application of the Fokker-Planck treatment to cold systems.

## 6. Comparison with theory

Our analysis of kinetic energy diffusion of Appendix A was developed with a smooth homogeneous background in equilibrium in mind. The surface density maps of our numerical cold discs models, indeed images of high-$z$ galaxies (e.g., the GOODS survey, Elmegreen et al. (2007)), on the contrary, generally point to a complex, fragmented mass distribution in the
<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Runs of angular momentum, radius, and velocity components as a function of time for a BH started on a circular orbit. a) The solid black curve graphs the angular momentum (left-hand axis with black labels). The right-hand axis is the scale for each velocity component (red curves). b) The cylindrical radius (in black, left-hand axis) is graphed together with the norm of the velocity (right-hand axis, in turquoise). The dotted line at the bottom is the *z* coordinate bound to ±10 pc. The radius and velocity components reach a minimum at *t* = 915 Myr.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 13.** Runs of angular momentum and radius for three BH orbits: circular; an elliptical orbit with a small eccentricity, *e* ≈ 0.55; and a third case with a large eccentricity, *e* ≈ 0.67. a) Run of the angular momentum for each case as indicated in the legend. The solid black curve graphs *L*<sub>*z*</sub> for the circular initial orbit displayed on Fig. 12. b) Cylindrical radius as a function of time for the same three cases. We note how the BH orbit with the more elliptical orbit actually evolves at a near-constant mean radius.</figcaption></figure>

core region visited by the BH. We limit our exploration to the
case when the BH starts on a radial orbit in a dynamically
cold disc (Table 1). We want to assess whether the basic pic-
ture of repeated gain and loss of azimuthal velocity is indeed at
work in the later stages of evolution, when the BH has acquired
angular momentum from the streaming of stars (see Fig. 2 and
Appendix A).

To that end, we graph on Fig. 14 the runs of specific angu-
lar momentum Lz and azimuthal velocity vφ as a function of the
cylindrical radius R for the BH orbit displayed on Figs. 7 and
8. The rotation curve for the model galaxy is displayed in solid
orange, and two phases of evolution are displayed with differ-
ent colours: prior to Lz reaching a maximum (in dark blue); and
at later stages when Lz decreases episodically (in red). The two
phases are shifted for clarity (the light-blue dotted lines con-
nect the end-point of one phase, to the start of the next), with
large light-blue arrows indicating the flow of time. The output
data was sampled at intervals of 10 Myr, from t = 1008 Myr to
t = 1800 Myr, with the transition through a peak Lz value at
t = 1522 Myr.

We discuss first the early phase, from 1008 to 1522 Myr (dark blue curves). Fig. 14(a) graphs continued increase in Lz with time, despite significant fluctuations (positive or negative) in radial migration. Let us focus on a time slice around the coordinates marked with the circled number❶. Near the symbol❶, the steady increase in Lz is evident as the radius increases from ≈1.2 to 1.5 kpc; notice how Lz (in blue) sits above the expected value (orange curve) at R ≈ 1.2 kpc. Consider now the evolution of vφ over the same interval, where we find again that the BH has excess azimuthal velocity at R = 1.2 kpc, sitting above the orange curve. Its larger azimuthal velocity means that it jumps ahead of the stars, which exert a negative acceleration on it. Therefore, vφ decreases in absolute terms as we move forward in time toward a larger radius R. The evolution at constant Lz would make vφ drop and fall on a path parallel to the black straight line, which is clearly not the case (Fig. 14[b]). In fact, the azimuthal velocity still decreases such that it crosses the orange curve, and soon the BH has a significant deficit of azimuthal velocity (dot-dash blue line, near❶ on the figure). The situation is now reversed, with the BH lagging behind the stellar flow
<figure><img src="image_17.png" /><figcaption>**Fig. 14.** A) Specific angular momentum, $L_z$, as a function of the cylindrical radius, $R$, during the orbital migration of the BH. The large sky-blue arrows indicate the flow of time, $t$. The orange solid lines trace the expected values for the model galaxy. The curve in blue (red) shows the coordinates before (after) $L_z$ has reached a maximum. The black lines indicate evolution at a constant $L_z$. The circles numbered 1 and 2 single out specific phases of evolution. B) Run of azimuthal velocity, $v_\phi$, against radius, $R$. The axes have been shifted for clarity. The colour-coding for the other curves and symbols is as for (A).</figcaption></figure>

at that radius: the dynamical traction exerted by the stars now pulls it forward and it soon catches up and overtakes the stellar flow. This see-saw pattern repeats itself up to $R \approx 1.8$ kpc which is near the point where the rotation curves flattens off; we find for our model galaxy that $dv_\phi/dR < 0$ when $R \ge 2$ kpc, off-scale on Fig. 14(b).

The second time slice of interest is marked with the circled number ②. In this later stage, the orbit evolves differently, with episodic losses of angular momentum, followed by time intervals of near-constant $L_z$. This is illustrated most clearly near the points ② on Fig. 14(a), when the horizontal black line (of constant angular momentum) traces well the orbit of the BH; in (b), the case of the azimuthal velocity, the black curve as a slope matching that of the data, in red. It is interesting to observe that the time intervals of near-constant $L_z$ are close to one full orbital time, or $\approx 100$ Myr. We suspected that this could be an indication of a bar-like structure or some other non-axisymmetric density fluctuations developing inside the BH's orbit. But the much increased radial stellar velocity dispersion, now exceeding $\sigma_R \approx 18$ km s$^{-1}$ close to $R = 0$ kpc (it started at $4.2$ km s$^{-1}$) works against bar- or any substructure formation. A casual inspection of density maps, not shown here, instead suggests that at the BH's orbit of $R \approx 1.7$ kpc where $\sigma_R \approx 8$ km s$^{-1}$ is much lower, the BH is still able to open a gap (albeit of low contrast) in the stellar distribution. In other words, the drop in angular momentum is most likely caused by gravitational torques, which here again work to lower its orbital angular momentum. We have not investigated the duration of such torques or their net impact quantitatively but suspect that they account for episodes of rapid angular momentum loss. By contrast, a steady decline in $L_z$ would be attributable to dynamical friction. This interpretation implies that the time-evolution of the BH's orbit is coupled to both the inner and outer disc stellar populations through their velocity dispersion, i.e. their local stability to Jeans modes of fragmentations.

A future study with spectral analysis should provide a more definitive diagnostic of angular momentum transport in these systems.

A related question which sheds light on the analysis is whether the time derivative $\langle \Delta \mathbf{v}_\parallel \rangle$ obtained numerically matches the values computed from Eqs. (A.3) and (A.4). We tailor these equations to the numerical setup by fixing $N_s = 1 \times 10^5$ and $m_\star = 1.6 \times 10^4$ M$_\odot$. We chose to round up the ratio $M_\star/m_\star$ from 731, as for the curves displayed on Fig. 2, to $\approx 10^3$. This simplifies the algebra and is of little consequence to the outcome. We also set $\sigma_\perp = \sqrt{2}\sigma_\star$ for the two-dimensional velocity dispersion. All quantities can be expressed in units of (km s$^{-1}$)$^2$/Myr with a single multiplicative constant. In the case of $\sigma_\star = 10$ km s$^{-1}$, we read on Fig. 2 a value of

$$v_\parallel\langle\Delta\mathbf{v}_\parallel\rangle \approx 40 \ (\text{km s}^{-1})^2/\text{Myr}$$

when $v_\parallel \approx 60$ km s$^{-1}$. This velocity component can also be read off Fig. 9(b) in the time interval $>1000$ Myr. If we isolate for $\langle \Delta \mathbf{v}_\parallel \rangle$ from the above relation, we find that $\langle \Delta \mathbf{v}_\parallel \rangle \approx 0.67$ km s$^{-1}$/Myr. A naive time-integration over the sampling interval $dt = 10$ Myr leads to $\Delta v_\parallel \approx 6.6$ km s$^{-1}$. Inspection of the output logs from the numerical integration gives us typical time-variations $\Delta v_\phi \approx 4$ to $5$ km s$^{-1}$, over the same time interval. In the case when $\sigma_\star = 20$ km s$^{-1}$, the same exercise leads to $\langle \Delta \mathbf{v}_\parallel \rangle \approx 1/6$ km s$^{-1}$/Myr and $\Delta v_\parallel \approx 1.66$ km s$^{-1}$ (same sampling rate). Thus we conclude that the analytical rates in Eq. (A.3) are of the right amplitude, but they fall short when the dispersion, $\sigma_\star$, increases significantly above $10$ km s$^{-1}$. Inspecting the runs for times $t \ge 329$ Myr we found a steady radial decrease in $\sigma_\star$ in the interval $R = 0$ to $2$ kpc, when $\sigma_\star$ falls from $\ge 20$ km s$^{-1}$, to $\approx 10$ km s$^{-1}$. In other words, the variations in $\Delta v_\phi$ become larger as we move outwards since locally $\sigma_\star$ is reduced. This comforts the conclusion that dynamical traction accounts quantitatively for much of the outward orbital evolution of the BH.
<figure><img src="image_18.png" /><figcaption>**Fig. 15.** Two examples of the interaction between the central BH and stellar sub-structures. In (A, top row), when the background spherical Isochrone potential is partially frozen, the disc stars form clumps of higher mass and binding energy. Such structures may dislodge the BH. This two-stage instability then grows as the BH extracts angular momentum from the flow of rotating stars around the centre. In (B, bottom row) the situation for a self-gravitating disc plus isochrone potential means that the clumps are much less massive and of a lower binding energy. In that case, a clump will dissolve in the strong tidal field of the BH. The centre of coordinates is indicated with a red dot. (See text for more details).</figcaption></figure>

## 7. Discussion

The problem of a massive perturber evolving in cored galactic potentials was revisited here with a focus on a background of stellar orbits with non-zero streaming motion drawn from a simple distribution function. A local treatment of gravitational focusing derived from Fokker-Planck diffusion coefficients was performed. The effects of resonant motion was left out of the analysis. A consequence of anisotropic gravitational focusing, we argued, is that the steady streaming of stars leads to an effective traction on the BH after a few orbital periods (see Appendix A and Fig. 3). Dynamical traction competes with the more familiar Chandrasekhar friction and we have illustrated the transition from one regime to the other for a BH on a radial orbit in a hot disc (Fig. 6). This is the expected evolution for low-angular momentum BH orbits.

Theoretical expectations are largely fulfilled in configurations where the streaming stars are stable against Jeans fragmentation modes (dynamically hot discs). This is not so with cold discs, where we found that stellar sub-structures are key players in the orbital evolution of the BH. Indeed stellar clumps shut out dynamical traction almost entirely, so a BH may reach the galactic centre on a timescale of ~300 Myr (compare Figs. 6 and 9). In this idealised framework, the formation of stellar clumps

is triggered by the response of the stellar orbits to the migrating BH (no clumps form when the BH is at rest at the origin: Fig. C.1). Thus when the BH reaches the centre, the stellar mass distribution is far from axial symmetry. The configuration will be unstable if and when sizeable clumps orbit near the centre. We spoke of a two-stage instability whereby an *m* = 1 mode pulls the BH out of the centre, at which point dynamical traction boosts its angular momentum and outward migration ensues. Hence the degree of binding energy of these clumps (their internal cohesion, not resolved properly in our numerical integrations) is crucial to the stability of the configuration. We address this specific point below.

### 7.1. Fragmentation and self-gravity

Fragments that develop through Toomre-unstable 'cold' stellar orbits in the disc triggered by the BH perturbation mimic observations of young galaxies undergoing active star formation. The extended and spherical (isochrone) component made up of 'hot' stellar orbits couples more weakly with the BH; we froze the spherical component for that reason in the reference numerical setup. A drawback is that the stellar clumps themselves do not exchange mechanical energy with the Isochrone potential (by
polarising orbits from that component, for instance). In that way
the binding energy of the clumps is enhanced artificially, and we
expect that relaxing the constraint of a frozen component will
change the outcome of the numerical integration significantly.
From hereon, we use the name ‘bulge’ as shorthand for spheri-
cal isochrone stellar component.

To inspect the impact of a live bulge, we performed two
calculations, the first with a fraction of live bulge orbits such
that they match in number (and mass) those of the live disc, or
100k mass elements each; and a second run with a fully self-
gravitating bulge and disc. The second run was performed with
the Amuse GPU-enhanced tree code BONSAI with a total of
8.4 × 10⁶ mass elements. The model parameters and bulge/disc
mass ratios are those of the standard model of Table 1. The ker-
nel sizes for numerical integration was reduced to l = 32 pc to
account better for the dynamics of the clumps.

Fig. 15 graphs for each case the distribution in space of a
subset of stars that form a bound sub-structure at a time when
the BH has reached the centre. The run shown on the top row
(a) for a partially self-gravitating bulge yielded stellar clumps
of strong binding energy. The run shows how the clumps merge
together and form a structure which far outweighs the BH (by
4× its mass). The large structure of ≈3318 mass elements seen at
the end of the sequence hovers close to the centre of coordinates
indicated with a red dot. At the end of the time sequence, we
find that 90.3% of the mass elements identified as belonging to
the clump are bound together. This is a prime example of the
important rôle that the net clump binding energy can have on the
stability of the BH’s position in the galaxy.

Compare these results with the second run displayed on
Fig. 15(b), where now the whole system is self-gravitating. Now
some 2887 mass elements add up to ≈3.8 × 10⁶ M☉ (1/3 the mass
of the BH). The important difference is that now only about 62%
of the mass of the clump has negative binding energy. As soon
as it approaches the centre the tidal field of the BH dissolves
the structure which forms a loose spiral pattern at the end of
the run. The BH remains at the heart of the system, and may
acrete stars and (possibly) form a nuclear star cluster. No insta-
bility will develop that can remove the BH from the centre unless
more massive sub-structures form in the course of time: nothing
of the sort happened over the course of 1.5 Gyr of evolution for
that model.

These calculations tell us that if the anisotropic stellar stream
is embedded in a large, much more massive hot stellar compo-
nent of nearly isotropic velocity field (like the current bulge),
then fragmentation modes that are triggered by a migrating BH
are unlikely to develop strong biding energy and create a situa-
tion whereby the massive compact object can be ejected from the
centre. The BH migrates to the centre on a timescale ~300 Myr
as anticipated from dynamical friction. We have not determined
the minimum mass ratio of isotropic- to anisotropic (disc-like)
components that will lead to this configuration: the calculations
that we performed suggest that this ratio has to be at least 2:1. It
would be interesting to extract the degree of anisotropy (essen-
tially, rotation) from calculations of galaxy formation simula-
tions that resolve the central kiloparsec adequately to provide a
better set of constraints on the initial configuration to use for a
refined study of BH migration.

7.2. *Cool and warm discs*

During the early stages of evolution of a low-$L_z$ orbit, the BH
undergoes dynamical friction. Two cases were considered, one
where the disc is dynamically warm, the other, cool. No sub-

structures form in a warm disc, and the evolution is very smooth
(Fig. 5).

The situation is much different in cooler environments, when
the disc fragments and density waves carry angular momentum
outward. The BH may well sink to the centre, but be ejected by a
symmetry-breaking configuration of stellar clumps which it trigger-
gered on the way in. When that happens. the BH acquires angular
momentum by dynamical traction from the surrounding stream-
ing stars. We have checked that the BH’s azimuthal velocity $v_φ$
increases at a rate well matched with Fokker-Planck energy dif-
fusion coefficients in that situation (Fig. 14). This is important
because these coefficients are obtained from local measurements
of rotation velocity, dispersion and surface brightness. Therefore
the dynamical state of an off-centred BH in an external galaxy
could possibly be gauged directly from a well-resolved young
galaxy with a rotating bulge (or highly anisotropic bulge). Sev-
eral such galaxies have long been observed at redshifts z ~ 1
to 3 (e.g., Elmegreen et al. 2007). The recent JWST detection of
an off-centre AGN at z ~ 7.3, when the Universe is ≈700 Myr
old (Übler et al. 2024), offers a great opportunity to tests and
understand better BH dynamics in the early stages of the forma-
tion of a galaxy. The environment of a BH in a young galaxy
should include gas dynamics and the formation of stars for a
realistic diagnostic of the internal dynamics at the core of young
galaxies. We stressed already how the binding energy of stel-
lar sub-structures, surely determined by dissipation processes at
formation, will impact on the stability of the BH at the centre of
galaxies. Analyses of the impact of gas on the orbit of a massive
perturber include, e.g., Ostriker (1999), Chapon et al. (2013) for
radiative gas; and the recent work by Lescaudron et al. (2023)
which includes turbulence. These last authors find a direct link
between the amplitude of the turbulent power-spectrum and the
orbital evolution of the BH. The density fluctuations brought by
turbulent modes parallel the role played here by stellar clumps.
There is a bridge to be built between such studies, because a
star-forming region will proceed from gas cooling in turbulent
configurations, and form large stellar aggregates (Renaud et al.
2021). It is intriguing that the BH discovered in the GA-NIFS
survey (Übler et al. 2024) of an estimated mass ≈5 × 10⁷ M☉ is
offset from the host galaxy’s photometric centre by R ≈ 1 kpc in
projection (see their Fig. 2; the broad-line region that is clearly
mis-aligned). Should the physical scales be confirmed, the BH
displacement would be on the same scale as the one we found
in our models. More detailed modelling is required to assess the
likelihood of an event such as this, yet such observations point
to a physically plausible configuration.

8. Conclusions

We have presented analytical arguments and showed with
numerical N-body calculations that the mechanics of dynami-
cal traction (see Sects. 3.3 and 5.4) work hand in hand with
the classic picture of orbital decay of a massive BH derived
from dynamical friction. They represent the same phenomena of
gravitational focusing but have opposite effects. An anisotropic
velocity field drawn from a d.f of the form $f(E, L_z)$ may lead to
a net gain of angular momentum by a massive BH and outward
migration. This is true for a dynamically warm stellar velocity
field. The competition between dynamical friction and traction
causes a continued transition in the solution space when one
effect comes to dominate over time.

A fragmented, dynamically cold configuration helps bring
the perturber BH to the barycentre faster, but these very frag-
ments make it difficult to keep it there. We identified a two-stage
instability that results in outward migration of a massive BH by dynamical traction. No dissipative terms were included in our analysis, and the scale-free nature of gravitational dynamics means that all quantities of time, mass, and length can be rescaled to fit observational constraints.

Our work was focused on a galactic system of total mass =1.16 × 10¹⁰ M☉ in equilibrium with a characteristic length a = 1.5 kpc. The orbital outward migration of an M☉ = 1.25 × 10⁷ M☉ BH in this system may reach ~1 kpc over a timescale of ~750 Myr to 1 Gyr. It is noteworthy that the recent AGN identified by the JWST at redshift z ≈ 7.3 shows an offset from its photometric centre by the same length in projection (Übler et al. 2024). Ongoing and future work with fully self-gravitating calculations should shed light on the dynamics of young galaxies by resolving the core region of these fascinating systems better.

**Acknowledgements.** We thank the anonymous referee for a thoughtful report that led to a much improved manuscript. We acknowledge G. Monari, B. Famaey and F. Renaud for their comments and suggestions about this work. K. Kraljic also pointed out several recent papers bearing on black hole dynamics in galaxies. The research reported here made extensive use of the community projects Agama and Amuse. We thank the developers for making their software available (and documented). Part of this work was supported by the French-Ukrainian DNIPRO research grant number Ts 17/2–1 awarded in 2022. Carbon footprints: this project made use of 1680 hours of organically grown intelligence (OGI), equivalent to ~30 kWh of energy (Balasubramanian 2021), together with 6000 CPUh of computational resources expended at the HPC centre of the University of Strasbourg (grants g2023a97-c and -g) during the academic year 2023-24. We thank M. Ringenbach and D. Brusson for support.

## References

Anglés-Alcázar, D., Özel, F., Davé, R., et al. 2015, ApJ, 800, 127

Argüelles, C. R., & Collazo, S. 2023, Universe, 9, 372

Bahcall, J. N., & Wolf, R. A. 1976, ApJ, 209, 214

Balasubramanian, V. 2021, Proc. Natl. Acad. Sci., 32, 118

Banik, U., & van den Bosch, F. C. 2021, ApJ, 912, 43

Banik, U., & van den Bosch, F. C. 2022, ApJ, 926, 215

Barnes, J., & Hut, P. 1986, Nature, 324, 446

Barro, G., Pérez-González, P. G., Kocevski, D. D., et al. 2024, ApJ, 963, 128

Bellovary, J. M., Cleary, C. E., Munshi, F., et al. 2019, MNRAS, 482, 2913

Binney, J., & Tremaine, S. 2008, Galactic Dynamics, 2nd edn. (Princeton: The University Press)

Bland-Hawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529

Bogdán, Á., Goulding, A. D., Natarajan, P., et al. 2024, Nat. Astron., 8, 126

Boily, C. M., Padmanabhan, T., & Paiement, A. 2008, MNRAS, 383, 1619

Boylan-Kolchin, M. 2023, Nat. Astron., 7, 731

Chandrasekhar, S. 1943, ApJ, 97, 255

Chapon, D., Mayer, L., & Teyssier, R. 2013, MNRAS, 429, 3114

Chatzopoulos, S., Fritz, T. K., Gerhard, O., et al. 2015, MNRAS, 447, 948

Chen, Z., Do, T., Ghez, A. M., et al. 2023, ApJ, 944, 79

Chiba, R., & Kataria, S. K. 2024, MNRAS, 528, 4115

Chu, A., Boldrini, P., & Silk, J. 2023, MNRAS, 522, 948

Cole, D. R., Dehnen, W., Read, J. I., & Wilkinson, M. I. 2012, MNRAS, 426, 601

de Vita, R., Trenti, M., & MacLeod, M. 2018, MNRAS, 475, 1574

Elmegreen, D. M., Elmegreen, B. G., Ravindranath, S., & Coe, D. A. 2007, ApJ, 658, 763

Finkelstein, S. L., Bagley, M. B., Ferguson, H. C., et al. 2023, ApJ, 946, L13

Fouvry, J.-B., Binney, J., & Pichon, C. 2015, ApJ, 806, 117

Fritz, T. K., Chatzopoulos, S., Gerhard, O., et al. 2016, ApJ, 821, 44

Fujii, M., Iwasawa, M., Funato, Y., & Makino, J. 2007, PASJ, 59, 1095

Fukushige, T., & Heggie, D. C. 2000, MNRAS, 318, 753

Gaia Collaboration (Drimmel, R., et al.) 2023, A&A, 674, A37

Heggie, D., & Hut, P. 2003, The Gravitational Million-Body Problem: A Multidisciplinary Approach to Star Cluster Dynamics (IOP Publishing Ltd)

Hopkins, P. F., & Quataert, E. 2011, MNRAS, 415, 1027

Julian, W. H., & Toomre, A. 1966, ApJ, 146, 810

Kalnajs, A. J. 1972, in IAU Colloq. 10: Gravitational N-Body Problem, ed. M. Lecar, Astrophys. Space Sci. Lib., 31, 13

Kaur, K., & Sridhar, S. 2018, ApJ, 868, 134

Kaur, K., & Stone, N. C. 2022, MNRAS, 515, 407

Kley, W., & Nelson, R. P. 2012, ARA&A, 50, 211

Kocsis, B., & Tremaine, S. 2015, MNRAS, 448, 3265

Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511

Labbé, I., van Dokkum, P., Nelson, E., et al. 2023, Nature, 616, 266

Lapiner, S., Dekel, A., & Dubois, Y. 2021, MNRAS, 505, 172

Lau, J. Y., & Binney, J. 2021, MNRAS, 507, 2241

Launhardt, R., Zylka, R., & Mezger, P. G. 2002, A&A, 384, 112

Lescaudron, S., Dubois, Y., Beckmann, R. S., & Volonteri, M. 2023, A&A, 674, A217

Lynden-Bell, D., & Kalnajs, A. J. 1972, MNRAS, 157, 1

Matthee, J., Naidu, R. P., Brammer, G., et al. 2024, ApJ, 963, 129

Merritt, D. 2004, Phys. Rev. Lett., 92, 201304

Merritt, D. 2013, Dynamics and Evolution of Galactic Nuclei, 1st edn. (Princeton, NJ: Princeton Univ. Press)

Merritt, D., Berczik, P., & Laun, F. 2007, AJ, 133, 553

Mo, H., van den Bosch, F. C., & White, S. 2010, Galaxy Formation and Evolution, 1st edn. (Princeton, NJ: Princeton Univ. Press)

Naidu, R. P., Oesch, P. A., van Dokkum, P., et al. 2022, ApJ, 940, L14

Neumayer, N., Seth, A., & Böker, T. 2020, A&A Rev., 28, 4

Nieuwmuaster, N., Schultheis, M., Sormani, M., et al. 2024, A&A, 685, A93

Nogueras-Lara, F., Schödel, R., Gallego-Calvente, A. T., et al. 2020, Nat. Astron., 4, 377

Ostriker, E. C. 1999, ApJ, 513, 252

Panamarev, T., & Kocsis, B. 2022, MNRAS, 517, 6205

Papaloizou, J. C. B. 2021, Celest. Mech. Dyn. Astron., 133, 30

Pelupessy, F. I., van Elteren, A., de Vries, N., et al. 2013, A&A, 557, A84

Petts, J. A., Gualandris, A., & Read, J. I. 2015, MNRAS, 454, 3778

Pfister, H., Volonteri, M., Dubois, Y., Dotti, M., & Colpi, M. 2019, MNRAS, 486, 101

Pfister, H., Dai, J. L., Volonteri, M., et al. 2021, MNRAS, 500, 3944

Portegies Zwart, S., & McMillan, S. 2018, Astrophysical Recipes: The art of AMUSE, IOP Astronomy/AAS (Bristol, UK: IOP Publishing)

Portegies Zwart, S., McMillan, S. L. W., van Elteren, E., Pelupessy, I., & de Vries, N. 2013, Comput. Phys. Commun., 184, 456

Rauch, K. P., & Tremaine, S. 1996, New Astron., 1, 149

Read, J. I., Goerdt, T., Moore, B., et al. 2006, MNRAS, 373, 1451

Reid, M. J., & Brunthaler, A. 2004, ApJ, 616, 872

Renaud, F., Gieles, M., & Boily, C. M. 2011, MNRAS, 418, 759

Renaud, F., Romeo, A. B., & Agertz, O. 2021, MNRAS, 508, 352

Schawinski, K., Koss, M., Berney, S., & Sartori, L. F. 2015, MNRAS, 451, 2517

Spitzer, L. 1987, Dynamical Evolution of Globular Clusters (Princeton University Press)

Šubr, L., & Haas, J. 2014, ApJ, 786, 121

Toomre, A. 1964, ApJ, 139, 1217

Tremaine, S., & Weinberg, M. D. 1984, MNRAS, 209, 729

Tremaine, S., Richstone, D. O., Byun, Y.-I., et al. 1994, AJ, 107, 634

Übler, H., Maiolino, R., Pérez-González, P. G., et al. 2024, MNRAS, 531, 355

Vasiliev, E. 2019, MNRAS, 482, 1525

Volonteri, M. 2012, Science, 337, 544

von Fellenberg, S. D., Gillessen, S., Stadler, J., et al. 2022, ApJ, 932, L6

Weinberg, M. D. 1989, MNRAS, 239, 549

Weinberg, M. D. 1998, MNRAS, 299, 499

Zoccali, M., Rojas-Arriagada, A., Valenti, E., et al. 2024, A&A, 684, A214
## Appendix A: Fokker-Planck approach

When a massive body moves through space at velocity $v$, multiple interactions with background stars lead to the diffusion of kinetic energy by repeated small-angle scattering. This is put on a quantitative footing by computing the velocity diffusion coefficients as functions of the local stellar density $\rho_*$ and velocity distribution function $f_x(v)$. The scattering events are summed up by integrating over the whole of space. When stars are distributed isotropically about the orbital trajectory of the massive body, pair-wise deflections cancel out and the net velocity kicks orthogonal to the trajectory add up to zero (to first order). For that reason, it is useful to split the net rate of change<fnref n="7" /> $\langle\Delta v\rangle$ in parallel and orthogonal components relative to the direction of motion of the body,

$$ \langle\Delta v\rangle = \langle\Delta v_\parallel\rangle + \langle\Delta v_\perp\rangle, \quad (A.1) $$

so that the velocity vector is incremented as

$$ v \rightarrow v + \Delta v = v + \langle\Delta v_\parallel\rangle \delta t + \langle\Delta v_\perp\rangle \delta t $$

over the time interval $\delta t$. The kinetic energy per unit mass of the body, $E = v \cdot v/2$, becomes

$$ E' = \frac{1}{2} (v + \Delta v)^2 = \frac{1}{2} v \cdot v + v \cdot \Delta v + \frac{1}{2} (\Delta v_\parallel)^2 + \frac{1}{2} (\Delta v_\perp)^2 $$

so that the rate of change $\langle\Delta E\rangle = \langle E' - E\rangle$ may be expressed as

$$ \langle\Delta E\rangle = v_{\parallel} \cdot \langle\Delta v_{\parallel}\rangle + \frac{1}{2} \langle(\Delta v_{\parallel})^2\rangle + \frac{1}{2} \langle(\Delta v_{\perp})^2\rangle. \quad (A.2) $$

In this equation, the 1st-order orthogonal term $v_\perp \cdot \langle\Delta v_\perp\rangle$ cancels by an argument of symmetry. The diffusion coefficients in (A.2) are computed from the (local) stellar velocity distribution function $f_x(v_\star)$ in integral form,

$$ \langle\Delta v_\parallel\rangle = \Gamma \left(1 + \frac{M_\bullet}{m_\star}\right) \frac{d}{dv} h(v|f_x(v_\star)) \quad (A.3a) $$

$$ \frac{1}{2} \langle (\Delta v_\parallel)^2 \rangle = \frac{1}{2} \Gamma \frac{d^2}{dv^2} g(v | f_x(v_\star)) \quad (A.3b) $$

$$ \frac{1}{2} \langle (\Delta v_\perp)^2 \rangle = \frac{1}{v} \Gamma \frac{d}{dv} g(v | f_x(v_\star)) \quad (A.3c) $$

where $\Gamma = 4\pi G^2 m_\star^2 \ln \Lambda$ accounts for the Coulomb term ($\Lambda$); integrals over the stellar distribution function $f_x(v_\star)$ define the Rosenbluth potentials $g, h$ as (Binney & Tremaine 2008, §7.4.4, Eqs. 7.83)

$$ \begin{aligned} h(v | f_x(v_\star)) &\equiv \int \frac{f_x(v_\star)}{\|v - v_\star\|} \mathrm{d}^3 v_\star \\ g(v | f_x(v_\star)) &\equiv \int f_x(v_\star) \|v - v_\star\| \mathrm{d}^3 v_\star . \end{aligned} \quad (A.4) $$

To ease the evaluation of the integrals (A.4), define $\mathcal{K}(v_\phi, v_\perp)$ as

$$ \mathcal{K}^2 \equiv (R\Omega - v_\phi)^2 + v_\perp^2 \quad (A.5) $$

and pick the sign of $\mathcal{K}$ to be that of $R\Omega - v_\phi$ (positive for a fast pattern speed / slow azimuthal velocity). From this function we obtain the derivative

$$ \frac{d\mathcal{K}}{dv} = -\frac{1}{\mathcal{K}} \left( R\Omega \frac{v_\phi}{v} - v \right) \quad (A.6) $$

where all quantities are scalars, $v$ being the norm of the velocity, and the chain rule may be applied to the Rosenbluth potentials (A.4) to compute derivatives which appear in (A.3), i.e. through the operator

$$ \frac{d}{dv} = \frac{d\mathcal{K}}{dv} \times \frac{d}{d\mathcal{K}}. $$

For instance, the second order derivative in (A.3b) becomes

$$ \frac{d^2}{dv^2} g = \frac{d}{dv} \left[ \frac{d\mathcal{K}}{dv} \frac{d}{d\mathcal{K}} \right] g = \frac{d^2 g}{d\mathcal{K}^2} \left( \frac{d\mathcal{K}}{dv} \right)^2 + \frac{dg}{d\mathcal{K}} \frac{d^2 \mathcal{K}}{dv^2}. $$

The first- and second-order derivatives of the Rosenbluth potential with respect to $\mathcal{K}$ are given in Appendix B for completeness.

The Fokker-Planck equations consists in integrating the diffusion coefficients to determine the time-evolution of the distribution function in equilibrium. The application to a single massive body reduces the Fokker-Planck treatment to a response equation where the first-order term in (A.2) dominates over the quadratic ones: this leads to the classic dynamical friction evolution with $dE/dt < 0$ due to a force term parallel but in the opposite sense to the motion of the body (see Fig. 1a and e.g. Binney & Tremaine 2008). The loss of kinetic energy and inspiral of the massive body appears to be common to any distribution function giving rise to an isotropic velocity field.

### A.1. Anisotropic configurations, dynamical traction

The situation is different if we focus on orbits that are not isotropic in space. A second look at Fig. 1 shows that the relative velocity with the background stars would remain the same if the massive body was at rest, and all the stars were flowing down the velocity vector, $-v_\bullet$. From that point of view, the dynamical friction becomes instead dynamical traction. Dynamical traction would drag a BH in any direction defined by the bulk background flow. Contrary to a dissipative frictional term, however, the traction on a body increases its kinetic energy and can be non-zero even if it is at rest in an inertial frame. For this to happen, the velocity field must be anisotropic (as on Fig. 1[b]).

This suggest that we embed a BH in a background potential made of a hot, isotropic component, and a second, cool component, made of stars drawn from a family of orbits carrying a net positive angular momentum. The situation is sketched on Fig. 1b. To simplify the discussion, we plot only stars rotating about the system barycentre and not the hot component. We suppose that the BH has come to rest at a distance $\approx \delta x_\bullet$ given by (4)). The bulk motion is now in azimuth and we may rewrite (A.2) in cylindrical coordinates as

$$ \langle\Delta E\rangle = v_\phi \cdot \langle\Delta v_\phi\rangle + \frac{1}{2} \langle(\Delta v_\phi)^2\rangle + \frac{1}{2} \langle(\Delta v_\perp)^2\rangle. \quad (A.7) $$

The azimuthal acceleration $\dot{v}_\phi = \langle\Delta v_\phi\rangle$ gives rise to a torque

$$ \mathbf{T} = \mathbf{r} \times \mathbf{M}_\bullet \dot{\mathbf{v}} = \|\delta \mathbf{x}_\bullet \cdot \langle\Delta \mathbf{v}_\phi\rangle\| \mathbf{M}_\bullet \hat{\mathbf{z}} $$

where we have combined the radial and z component of the velocity in the orthogonal component $v_\perp$. This net torque will transfer angular momentum, from the stars, to the BH. It is then a simple matter to estimate the magnitude of the angular momentum $L_z$ accrued by the BH as

$$ L_z \approx \|\mathbf{T}\| \delta t = \|\delta \mathbf{x}_\bullet \cdot \langle\Delta \mathbf{v}_\phi\rangle\| \mathbf{M}_\bullet \frac{2\pi}{\omega}. \quad (A.8) $$

<footnote n="7">⁷ We follow the syntax of (Spitzer 1987) and replace the time derivative operator $d/dt$ with angled-brackets $\langle\dots\rangle$.</footnote>
In deriving $L_z$, we assumed the torque to be constant over one full period of the stars. They are bound together by their self-gravity and the potential generated by the hot component, which we took to be spherically symmetric and uniform (no aggregates, fragments, or radial profile).

## Appendix B: Computing the diffusion coefficients

For completeness we give here an outline of the integrals leading to (A.3). Starting with the distribution function (7), we compute the Rosenbluth potential $g(v|f_x(v_\star))$ in (A.4) as a function of the norm $v$ of the velocity of the massive body by expanding the vector dot products. As a simplification, we extend the upper bound of integration in velocity space to infinity. Inserting our choice of d.f. (7), the integral transforms to

$$ 
\begin{aligned}
g[v|f_x(v_\star)] &= \int d^3v_* f_x(v_\star) ||v - v_*|| \\
&= \int \pi dv_{\star,\perp}^2 dv_{\star,\phi} f_{\perp}(v_{\star,\perp}) \mathcal{D}(R\Omega - v_{\star,\phi}) \times \\
&\qquad \sqrt{[v - v_*] \cdot [v - v_*]} \quad (B.1) \\
&= \int 2\pi v_{\star,\perp} dv_{\star,\perp} f_{\perp}(v_{\star,\perp}) \left[ (R\Omega - v_\phi)^2 + v_{\perp}^2 + v_{\star,\perp}^2 \right]^{1/2}
\end{aligned}
$$

and the term under the radical defines the space-dependent parameter $\mathcal{K}(R, z)$ as

$$ \mathcal{K}^2(R, z = 0) = (R\Omega - v_\phi)^2 + v_{\perp}^2. $$

The integral (B.1) is easily computed for the Gaussian form (8) taken for $f_\perp$. The simplification of a Dirac operator implies that the azimuthal motion of the background stars is largely dominated by the streaming flow and we neglect the energy invested in random motion for that degree of freedom. The same calculation leading to (B.1) can be repeated to compute the first Rosenbluth potential $h(v|f_x(v_\star))$. We find

$$ 
\begin{aligned}
h[v|f_x(v_\star)] &= \int_0^\infty \frac{\pi dv_{\star,\perp}^2 f_\perp(v_{\star,\perp})}{[\mathcal{K}^2 + v_{\star,\perp}^2]^{1/2}} \\
&= -2\pi \sqrt{\mathcal{K}^2} f_\perp(0) + \int_0^\infty \frac{\pi v_{\star,\perp}}{\sigma_*^2} dv_{\star,\perp} f_\perp(v_{\star,\perp}) \\
&\qquad \times \sqrt{\mathcal{K}^2 + v_{\star,\perp}^2} \quad (B.2)
\end{aligned}
$$

where the last integral is once more easily integrated by substituting (8). We note that $\sqrt{\mathcal{K}^2}$ yields the absolute value $|\mathcal{K}|$; the expression can be used in symbolic calculations to avoid conditional expressions.

The task of computing the diffusion coefficients reduces to computing the first- and second-order derivatives of $h$ and $g$ with respect to $\mathcal{K}$, and then use the relation (A.6). After some straightforward algebra we find (dropping the dependence on d.f. $f_x(v_\star)$)

$$ 
\begin{aligned}
\frac{d}{d\mathcal{K}}h(v) &= \mp 2\pi f_{\perp}(0) + \frac{\mathcal{K}}{\sqrt{2\sigma}} f_{\perp}(0) \left[ \sqrt{\pi} e^{\mathcal{K}^2/2\sigma_*^2} \text{erfc} \left( \frac{\mathcal{K}}{\sqrt{2\sigma_*}} \right) \right] \\
\frac{d}{d\mathcal{K}}g(v) &= \frac{2}{3}\pi f_{\perp}(0) \left\{ -3\mathcal{K}^2 + 2\mathcal{K} e^{\mathcal{K}^2/2\sigma_*^2} \sqrt{2\sigma^2\pi} \right\} \quad (B.3)
\end{aligned}
$$

$$ \frac{d^2}{d\mathcal{K}^2} g(v) = \frac{2}{3}\pi f_{\perp}(0) \left\{ -6\mathcal{K} + 2\sqrt{2\sigma_*^2\pi} e^{\mathcal{K}^2/2\sigma_*^2} \left( 1 + \frac{\mathcal{K}^2}{\sigma^2} \right) \right\} $$

and recall that $f_\perp(0) = \rho(R, z)/(2\pi m_* \sigma_*^2)$. The $\mp$ multiplicative factor in the first of Eq. (B.3) carries the sign of the derivative $d\sqrt{\mathcal{K}^2}/d\mathcal{K}$. We note that ERFC is the complementary error function.

## Appendix C: Validation tests

We ran tests with an axisymmetric configuration consisting of the isochrone potential, to which we added the Miyamoto-Nagai stellar disc sampled with the Agama-generated self-consistent distribution functions. The BH was not included in the generation of the system's distribution function, but rather was put in afterwards, at rest, at the centre of coordinates as an $m = 0$ perturbation. Our reasoning here is that if the numerical procedure failed, it would produce a de facto growing instability that would persist. With a stable numerical setup, we expect a ring-like pattern to develop quickly before being washed out over a few dynamical times due to shearing and phase-mixing.

The results are displayed on Figs. C.1 and C.2. On the first figure, we graph the surface density of the stars from $t = 0$ up to 1008 Myr of evolution, corresponding to roughly 10 orbital periods at $R = 1$ kpc. Large-scale spiral patterns can be seen to develop and be present at all stages of evolution. The inner region runs differently, with a ring ($m = 0$) mode soon appearing after one period (frame $t = 102$ Myr) before vanishing at time $t = 308$ Myr and beyond. This is the expected evolution given the initial configuration. In the later stages of evolution, we find non-axisymmetric density fluctuations on a scale of ~1 kpc, but hardly any in the core region. The free motion of the BH triggers fluctuations in the potential which remain on an amplitude coherent with the numerical resolution of the calculation. For example, the BH remains confined to its radius of influence (~180 pc; see Fig. C.2[b]) and essentially bound by the numerical resolution, indicated with a black arrow on the figure. For reference, we give the amplitude of the linear velocity of harmonic motion on a scale of the BH radius of influence as a red arrow: the residual motion of the BH which builds up over ~1.5 Gyr of evolution is therefore in agreement with analysis and sets a reference for the interpretation of our results.

### C.1. Angular momentum conservation

A crucial aspect for our study of angular momentum transfer is the conservation of angular momentum by the stars, and by the BH itself. The time evolution of the BH's angular momentum and velocity vector are displayed on Fig. C.2(a). The $z$-component of the angular momentum is shown in black (left-hand vertical axis). It appears to evolve somewhat stochastically, and crucially begins to change sign after $t = 750$ Myr; the same applies to the radial and azimuthal velocity components, which oscillate about 0 km s$^{-1}$ (the statistical expectancy). Overall the total system angular momentum was conserved to a relative error of 0.1% after $t = 1000$ Myr of evolution, and to 0.4% after $t = 1500$ Myr: the numerical drift gives an important clue to what will come later when the BH is set in motion. This global behaviour is similar to Brownian motion due to the graininess of the potential (Merritt et al. 2007). This is confirmed by monitoring the (3D) radius and velocity over time (Fig. C.2[b]). We performed a second set of tests (no shown here) with the Barnes-Hut integrator replaced with the direct-summation Hermite code PH4 to ensure that the
<figure><img src="image_19.png" /><figcaption>Fig. C.1. Time evolution of a galactic disc with a BH at rest at the centre of coordinates. The black dot indicates the BH of mass = 1.25 × 10⁷ M☉. A white dot indicates a test star which remains on a steady orbit. The time runs from left to right as indicated. Notice how an over-density in the shape of a ring develops early on but is then wiped out by differential rotation. After ≈ 1 Gyr of integration time, the BH remains close to the barycentre with a velocity of ~7 km s⁻¹ in amplitude.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. C.2. Time evolution of the velocity, position and angular momentum of the BH as function of time. The BH initially sits at rest at the system's barycentre. a) The solid black curve graphs the angular momentum, which shows an irregular pattern of small amplitude; the scale's upper bound of ≈ 0.5 is the value expected for a radial orbit with $r = 0.125$ kpc at velocity $v_φ = 6.6$ km s⁻¹; the right-hand axis in red is the scale for individual velocity components; b) The radius (in black, left-hand axis) and velocity $|v_·|$ (in turquoise, right-hand axis) as function of time. The arrows indicate the spatial resolution = 0.125 kpc of the calculation, in black; and the expected velocity for a BH orbiting in the frozen background potential inside its influence radius, here ≈ 0.180 kpc, in turquoise. The black dotted line is the z coordinate. Significant yet resolved motion develops only after a runtime ~1.5 Gyr.</figcaption></figure>

trends are not enhanced or affected by the choice of integration method. The dispersion in the results obtained with different integrators for these reference configurations are standards against which to identify a genuine physical trend in the calculations.

## Appendix D: Stability analysis

We are interested in the stability of the stellar disc before a perturbation (the BH) is introduced. Details of the calculations and definitions can be found in Binney & Tremaine (2008).

### D.1. Jeans modes, Toomre's Q

The scale at which self-gravity becomes significant is given by the Jeans length,

$$ \lambda_J = \sqrt{\pi \sigma_*^2 / G \rho_0} \qquad (D.1) $$

where $σ_*$ is the three-dimensional stellar velocity dispersion, and $ρ_0$ the central density. We compute for $σ_* ≈ 7.3$ km s⁻¹ a

length $λ_J ≈ 478$ pc from (D.1), significantly larger than the vertical scale height $b = 250$ pc used for the initial configuration. Therefore, we do not expect significant evolution of the vertical disc structure.

A stellar distribution is stable (in a local sense) against all axisymmetric radial modes of fragmentation when (see §6.2 of Binney & Tremaine 2008)

$$ Q \equiv \frac{\sigma_{RK}}{3.36 G \Sigma} > 1. \qquad (D.2) $$

This Toomre (1964) Q parameter combines the radial epicyclic frequency $κ$, surface density $Σ$ as well as the radial velocity dispersion, $σ_R$. The epicyclic frequency is easily computed from (15) and (14) after summing over the potentials ((see Binney & Tremaine 2008, Eqs. 3.79)). This leads to a mean value

$$ \bar{Q} \approx 0.87, $$

with a trend of decreasing Q with increasing radius; in fact, Q peaks at $Q ≈ 1.2$ near the origin, and comes down to $Q ≈ 0.45$
<figure><img src="image_21.png" /><figcaption>Fig. C.3. The figure graphs the critical length $l_{\star}$ of Eq. (D.4) as a function of the cylindrical radius, $R$. The Jeans lengths $\lambda_J$ are shown in blue. The dotted line has $\lambda_J$ = constant, the central value. The solid blue line recomputes the Jeans length using the local density $\rho(R, 0)$. The dot-dash blue line is the situation when both the disc and halo are live components. Notice how $l_{\star} \approx \lambda_J$ around the origin. When $l_{\star} < \lambda_J$, the differential rotation contributes to erode a fragmentation mode.</figcaption></figure>

when $R = 5$ kpc. Therefore we expect over-dense patterns to form more easily outside the region $R = a$ while no or only a few Jeans-type fragments should develop in the nearly harmonic core. We recall that $Q \approx 1$ is a threshold value, and hence the self-gravity of large fluctuations may yet develop to some degree (see Renaud et al. 2021 for a critical assessment of $Q$).

## D.2. Differential rotation and self-gravity

The disc stars are not fully self-gravitating so the connection between velocity dispersion, angular velocity and surface density is not as in Toomre (1964)'s analysis leading to (D.2). To quantify the importance of the external potential, we compare the scale $l_{\star}$ at which differential rotation is important, to the Jeans length $\lambda_J$. We fix $z = 0$ throughout.

Two points separated radially by a length $l_{\star}$ have a relative azimuthal velocity to first-order in $l_{\star}/R$ given by

$$ \Delta v_{\phi} = l_{\star} \left( \Omega + R \frac{d\Omega}{dR} \right). \qquad (D.3) $$

We now focus on the situation where an over-dense clump of mass $M_J$ forms which spawns a volume of diameter $= 2\lambda_J$. A star approaching $M_J$ will not be bound to it if $E = -GM_J/\lambda_J + \frac{1}{2}(\Delta v_{\phi})^2 > 0$, or

$$ l_{\star} > \sqrt{\frac{2\pi G\rho}{3\Delta^2}} \lambda_J, \qquad (D.4) $$

where we have made use of the effective frequency<fnref n="8" />

$$ \Delta \equiv \left| \frac{\kappa^2}{2\Omega} - \Omega \right|. $$

Note that the density $\rho$ in (D.4) is that of the live (disc) component alone; and that $\Delta = \Omega$ in the limit when $d\Omega/dR = 0$

(solid body rotation). When $l_{\star} \le \lambda_J$, the shear dissolves a Jeans-unstable fragment and stabilises the disc against local collapse; otherwise when $l_{\star} > \lambda_J$, the differential rotation cannot stop the Jeans fragmentation mode from growing.

The transition from one regime to the other is displayed on Fig. C.3 for the reference configuration (see §4.1). The black curve on Fig. C.3 is the solution for $l_{\star}$ from (D.4). Inspection of the figure shows that shearing takes place on a scale significantly smaller than the local Jeans length $\lambda_J$ = constant in the region inside radius $R = 0.82$ kpc (dotted blue line). This holds true throughout the region $R = a = 1.5$ kpc if we use the local value of $\lambda_J$ (solid blue curve; the cross-over to $l_{\star} > \lambda_J$ occurs at $R = 1.55(3)$ kpc).

In conclusion, we find that the Toomre criterion points to a marginally unstable core region with $\bar{Q} \sim 0.87 < 1$ on average. Because the stellar orbits are not fully self-gravitating, we argue that differential rotation helps curtail the growth of fragmentation modes in the central region: they should not form there unless triggered by a perturbation. By comparison, the tidal field of a background potential almost always contributes to slow down fragmentation modes; they have a limited impact when the density gradients are weak, which is the case in the core region (see Binney & Tremaine 2008, §5); and e.g. Renaud et al. 2011 for a general tensor treatment of tidal fields). On the larger scale, however, we computed a low value $Q \approx 0.45$ and so from the analysis we can anticipate the formation of large spiral features outside the core length $a$. This is borne out by the outcome of numerical integrations shown on Fig. C.1.

<footnote n="8">⁸ Note that $\Delta = 0$ here does not have the usual meaning of a resonance since the ratio $\kappa:\Omega = \sqrt{2}:1$ is not rational.</footnote>