# Secular Dynamics around a Supermassive black hole via Multipole Expansion

Jean-Baptiste Fouvry, Walter Dehnen, Scott Tremaine, Ben Bar-Or

► To cite this version:

Jean-Baptiste Fouvry, Walter Dehnen, Scott Tremaine, Ben Bar-Or. Secular Dynamics around a Supermassive black hole via Multipole Expansion. The Astrophysical Journal, 2022, 931, 10.3847/1538-4357/ac602e. insu-03748270

HAL Id: insu-03748270

https://insu.hal.science/insu-03748270

Submitted on 9 Aug 2022

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Secular Dynamics around a Supermassive black hole via Multipole Expansion

Jean-Baptiste Fouvry<sup>1,2 ID</sup>, Walter Dehnen<sup>3,4,5</sup>, Scott Tremaine<sup>2,6 ID</sup>, and Ben Bar-Or<sup>2 ID</sup>

<sup>1</sup> Institut d'Astrophysique de Paris, UMR 7095, CNRS, Sorbonne Université, 98 bis Boulevard Arago, F-75014 Paris, France

<sup>2</sup> Institute for Advanced Study, Princeton, NJ, 08540, USA

<sup>3</sup> Astronomisches Rechen-Institut, Zentrum für Astronomie der Universität Heidelberg, Mönchhofstraße 12-14, D-69120 Heidelberg, Germany

<sup>4</sup> Universitäts-Sternwarte der Ludwig-Maximilians-Universität, Scheinerstraße 1, D-81679 München, Germany

<sup>5</sup> School of Physics and Astronomy, University of Leicester, University Road, Leicester LE1 7RH, UK

<sup>6</sup> Canadian Institute for Theoretical Astrophysics, University of Toronto, 60 St. George Street, Toronto, ON M5S 3H8, Canada

Received 2020 November 3; revised 2022 March 1; accepted 2022 March 21; published 2022 May 18

## Abstract

In galactic nuclei, the gravitational potential is dominated by the central supermassive black hole, so stars follow quasi-Keplerian orbits. These orbits are distorted by gravitational forces from other stars, leading to long-term orbital relaxation. The direct numerical study of these processes is challenging because the fast orbital motion imposed by the central black hole requires very small timesteps. An alternative approach, pioneered by Gauß, is to use the secular approximation of smearing out *N* stars over their Keplerian orbits, using *K* nodes along each orbit. In this study, we propose three novel improvements to this method. First, we reformulate the discretization of the rates of change of the variables describing the orbital states to ensure that all conservation laws are exactly satisfied. Second, we replace the pairwise sum over nodes by a multipole expansion up to order *l*<sub>max</sub>, reducing the overall computational cost from *O*(N<sup>2</sup>K<sup>2</sup>) to *O*(NKl<sub>max</sub><sup>2</sup>). Finally, we show that the averaged dynamical system is equivalent to 2*N* interacting unit spin vectors and provide two time integrators: a second-order symplectic scheme, and a fourth-order Lie-group Runge-Kutta method, both of which are straightforward to generalize to higher order. These new simulations recover the diffusion coefficients of stellar eccentricities obtained through analytical calculations of the secular dynamics.

*Unified Astronomy Thesaurus concepts: Galactic center (565); Galaxy nuclei (609); Gravitation (661); Celestial mechanics (211)*

## 1. Introduction

Supermassive black holes (BHs) are ubiquitous in external galaxies (Kormendy & Ho 2013), where their active feedback plays a critical role in regulating galaxy formation through cosmic time (Heckman & Best 2014). Yet, the details of their diet and their impact on the stellar cluster that surrounds them (the galactic nucleus) remain open and challenging questions. Indeed, galactic nuclei are among the densest stellar systems in the universe. Despite the high stellar density, the gravitational potential in galactic nuclei is dominated by the central supermassive BH. As a result, stars follow quasi-Keplerian orbits, which are slowly distorted by the additional perturbations present in the system.

The steep potential well generated by the central super-massive BH implies the existence of a wide range of dynamical timescales in the system, and the evolution of the stellar cluster involves numerous dynamical processes acting on radically different timescales (Rauch & Tremaine 1996; Hopman & Alexander 2006; Merritt 2013; Alexander 2017). These successively include (i) the dynamical time associated with the fast Keplerian motion; on timescales longer than this, the stellar orbits can be regarded as eccentric massive wires; (ii) the in-plane precession time of the Keplerian wires generated by the relativistic corrections from the BH and the stellar mean potential; (iii) the vector resonant relaxation time (see, e.g., Kocsis & Tremaine 2015; Fouvry et al. 2019), which, owing to nonspherical stellar fluctuations and the relativistic corrections

induced by a spinning BH, leads to the reshuffling of the orientations of the orbital planes of the wires; (iv) the scalar resonant relaxation time (see, e.g., Rauch & Tremaine 1996; Bar-Or & Alexander 2016; Sridhar & Touma 2016; Bar-Or & Fouvry 2018), during which resonant torques between the precessing wires lead to a diffusion of the wires' eccentricities; and (v) the nonresonant relaxation time (Bahcall & Wolf 1976; Lightman & Shapiro 1977; Cohn & Kulsrud 1978; Bar-Or et al. 2013), during which localized two-body encounters between stars lead to the long-term relaxation of the stars' Keplerian energy, i.e., the semimajor axes of the wires.

As a result of this wide range of dynamical times, from a few years for the fast Keplerian motion (even a few minutes for stars near the event horizon) up to a Hubble time for nonresonant relaxation in the nucleus of the Milky Way, direct numerical simulations of these dynamical systems remain very challenging. They were first performed with grid methods (see, e.g., Jacobs & Sellwood 2001; Kazandjian & Touma 2013); it is only recently that an effective direct simulation of a galactic nucleus with *N* = 10<sup>6</sup> stars has been presented (Panamarev et al. 2019), and even in this case, most of the stars lie outside the sphere of influence of the central BH. Conversely, simulations in the very relativistic regime are still limited to a small number of particles *N* ~ 10<sup>2-3</sup> whether they use direct *N*-body methods (Merritt et al. 2011) or effective methods (see, e.g., Madigan et al. 2011; Hamers et al. 2014).

To circumvent these intrinsic difficulties, one has to resort to additional assumptions. Traditionally, the secular approximation smears out the stars along their underlying fast Keplerian motion, in other words, replaces stars with Keplerian wires. Describing the dynamics of the stellar cluster amounts then to describing the long-term evolution of the orbital parameters of the wires. This is in particular at the heart of the

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

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
Gauß method (Touma et al. 2009), which provides an efficient algorithm for computing the force between two such wires. If one is interested in the process of vector resonant relaxation, i.e., the relaxation of the orientations of orbital planes, this same approach can be further leveraged to also average the wire dynamics over their in-plane precession, replacing Keplerian wires with Keplerian annuli, which offers new venues for performing numerical investigations of these long-term dynamics (Kocsis & Tremaine 2015).

The main benefit of these approaches is that any explicit average over the fast orbital motion offers a reduction of the range of timescales in the system. These methods can then explore longer timescales that are beyond the reach of naive direct methods. Yet, most of these approaches suffer from relying on the computation of all the individual forces between objects, i.e., these methods come with a numerical complexity scaling like $O(N^2)$, with N the total number of stars.

In the present paper, we show how a multipole expansion (Hénon 1964; Aarseth 1967; Henon 1973; van Albada & van Gorkom 1977; Fry & Peebles 1980; Villumsen 1982; White 1983; McGlynn 1984; Meiron et al. 2014; Dehnen 2014) yields a numerical scheme that can integrate the secular dynamics of Keplerian wires with a complexity scaling like $O(NK\ell_{\text{max}}^2)$, with N the total number of stars, K a parameter independent of N, and $\ell_{\text{max}}$ the maximum harmonics considered in the multipole expansion. We also show that the system is equivalent to 2N classical spin vectors, and we devise time-integration schemes that exactly comply with this system's geometric constraints.

The present paper is organized as follows. In Section 2 we describe the Hamiltonian and obtain its orbit average to account for the domination of the central BH. In Section 3 we derive the equations of motion for the orbital elements of the Keplerian wires. Section 4 details our numerical approach of discretizing averages over the Keplerian motions as sums over nodes, and shows how to develop a multipole expansion of the Newtonian pairwise interactions. In Section 5 we describe numerical time-integration schemes appropriate for the system of N Keplerian wires. In Section 6 we illustrate these new numerical methods by measuring the diffusion coefficients of stellar eccentricities and by comparing them with predictions from kinetic theory. In Section 7 we discuss the limitations of the present approach and possible future improvements, and we conclude in Section 8.

## 2. The Orbit-averaged Hamiltonian

We consider N stars with masses $m_i$ and positions $\mathbf{R}_i$, orbiting a supermassive BH of mass $M_\bullet$ at location $\mathbf{R}_\bullet$. The total Hamiltonian of this system is

$$H = \frac{\mathbf{P}_\bullet^2}{2M_\bullet} + \sum_{i=1}^{N} \left[ \frac{\mathbf{P}_i^2}{2m_i} - \frac{GM_\bullet m_i}{|\mathbf{R}_i - \mathbf{R}_\bullet|} \right] + \sum_{i=1}^{N} m_i \Phi_{\text{GR}}^i - \sum_{i<j}^{N} \frac{Gm_i m_j}{|\mathbf{R}_i - \mathbf{R}_j|}, \quad (1)$$

with the canonical momenta $\mathbf{P}_\bullet = M_\bullet \dot{\mathbf{R}}_\bullet$ and $\mathbf{P}_i = m_i \dot{\mathbf{R}}_i$. In this equation, the potential contribution $\Phi_{\text{GR}}^i$ accounts for the (conservative) relativistic corrections induced by the central BH,<fnref n="7" /> i.e., the Schwarzschild and Lense-Thirring precessions (see, e.g,

Merritt 2013), whose detailed expressions are given in Appendix B.

### 2.1. Democratic Coordinates

In order to emphasize the dominant influence of the central BH on the system dynamics, we rewrite Equation (1) using democratic coordinates centered on the BH (Duncan et al. 1998) and their canonical momenta,

$$\mathbf{r}_\bullet = \frac{1}{M_{\text{tot}}} \left[ \mathbf{M}_\bullet \mathbf{R}_\bullet + \sum_{i=1}^{N} m_i \mathbf{R}_i \right], \quad \mathbf{r}_i = \mathbf{R}_i - \mathbf{R}_\bullet, \quad (2a)$$

$$\mathbf{p}_\bullet = \mathbf{P}_\bullet + \sum_{i=1}^{N} \mathbf{P}_i, \quad \mathbf{p}_i = \mathbf{P}_i - \frac{m_i}{M_{\text{tot}}} \left[ \mathbf{P}_\bullet + \sum_{i=1}^{N} \mathbf{P}_i \right], \quad (2b)$$

where we introduced the total mass $M_{\text{tot}} = M_\bullet + M_\bullet$, with $M_\bullet = \sum_{i=1}^N m_i$ the total stellar mass. Thus, $\mathbf{p}_\bullet$ is just the total momentum, which we set to zero without loss of generality, so that $\mathbf{p}_i = \mathbf{P}_i$ is the barycentric momentum of the ith star. Following this change of coordinates, the Hamiltonian (1) becomes

$$H = \sum_{i=1}^{N} \left[ \frac{\mathbf{p}_i^2}{2m_i} - \frac{GM_\bullet m_i}{|\mathbf{r}_i|} \right] + \sum_{i=1}^{N} m_i \Phi_{\text{GR}}^i - \sum_{i<j}^{N} \frac{Gm_i m_j}{|\mathbf{r}_i - \mathbf{r}_j|} + \frac{1}{2M_\bullet} \left[ \sum_{i=1}^{N} \mathbf{p}_i \right]^2. \quad (3)$$

The first term is the sum of N independent Kepler Hamiltonians, the second term is associated with the relativistic corrections to stellar motion induced by the central BH, the third term captures the pairwise interactions between the stars, and finally, the last term is the kinetic energy of the central BH. As we describe below, after averaging over the fast Keplerian motions, the first and last terms become irrelevant constants.

### 2.2. Orbital Elements

In order to describe the Keplerian dynamics imposed by the central BH, we transform the $(\mathbf{r}, \mathbf{p})$ coordinates to the Delaunay variables $(M, \omega, \Omega, \Lambda, L, L_z)$ (Binney & Tremaine 2008) for each star. In this notation, the dynamical angles $M, \omega$, and $\Omega$ are the mean anomaly, the argument of pericenter, and the longitude of the ascending node, respectively. The associated actions are

$$\Lambda = m \sqrt{GM_\bullet a}, \quad L = \Lambda \sqrt{1 - e^2}, \quad L_z = L \cos I, \quad (4)$$

where $\Lambda$ is the circular angular momentum of an orbit with the same energy or same semimajor axis $a$, $L$ is the magnitude of the angular momentum vector, and $L_z$ is its projection onto the z-axis, while $e$ and $I$ denote the eccentricity and inclination of the orbit, respectively. Kepler's equation

$$M = E - e \sin E \quad (5)$$

introduces the eccentric anomaly $E$, which relates to the orbital radius, i.e., the distance from the central BH, via $r = a(1 - e \cos E)$.

<footnote n="7">⁷ The notation $\Phi_{\text{GR}}^i$ is somewhat misleading, as the relativistic changes given in Appendix B cannot rigorously be derived from a potential. We nonetheless use this notation for convenience.</footnote>
Following this change of variables, Equation (3) becomes

$$H = -\sum_{i=1}^{N} \frac{m_i^3 (GM_i)^2}{2\Lambda_i^2} + H_{GR} + H_{\star} + \frac{1}{2M_{\bullet}} \left[ \sum_{i=1}^{N} \mathbf{p}_i \right]^2 . \quad (6)$$

with

$$H_{GR} = \sum_{i=1}^{N} m_i \Phi_{GR}^{i}, \quad H_{\star} = -\sum_{i<j}^{N} \frac{Gm_i m_j}{|\mathbf{r}_i - \mathbf{r}_j|} . \quad (7)$$

## 2.3. Orbit Average

The first term in the Hamiltonian (6) describes N Keplerian orbits around the BH, and only depends on a single set of actions, $\Lambda_i$. As a result, under this Hamiltonian, all variables but $M_i$ are conserved, and the mean anomaly M evolves with the Keplerian orbital frequency,

$$\Omega_{\text{Kep}}(a) = \sqrt{\frac{GM_{\bullet}}{a^3}} . \quad (8)$$

If the BH is supermassive, i.e., if one has $M_{\bullet} \gg M_{\star}$, the dynamics is dominated by this fast Keplerian motion and the evolution of all other variables is much slower. Therefore, we average the stellar dynamics over these fast orbital motions through the so-called secular approximation to obtain the orbit-averaged Hamiltonian,

$$\langle H \rangle (\{\omega, \Omega, \Lambda, L, L_z\}_i) \equiv \oint \frac{dM_1}{2\pi} \dots \oint \frac{dM_N}{2\pi} H. \quad (9)$$

Here and in the remainder of this paper, we use the notation $\langle \cdot \rangle$ for an average over all unperturbed stellar orbits.

Since $\langle H \rangle$ is independent of $M_i$, the associated actions $\Lambda_i$ are conserved and the first term of Equation (6) averages to a constant. Upon expanding, the last term of Equation (6) consists of the orbit-averaged stellar kinetic energies $\langle \mathbf{p}_i^2 \rangle$ and products of the orbit-averaged stellar momenta $\langle \mathbf{p}_i \rangle$. In our barycentric frame, $\langle \mathbf{p}_i \rangle = 0$ and $\langle \mathbf{p}_i^2 \rangle = (m_i^2 GM_{\bullet}/\Lambda_i)^2$ by virtue of the virial theorem. The constant terms depending only on $\Lambda_i$ do not induce any dynamics and can be omitted, so the orbit-averaged Hamiltonian (9) finally becomes

$$\langle H \rangle = \langle H_{GR} \rangle + \langle H_{\star} \rangle . \quad (10)$$

## 3. The Equations of Motion

The orbit-averaged Hamiltonian (10) describes the dynamics of N gravitationally coupled Keplerian wires subject to relativistic precession. Each wire is characterized by the five quantities $\omega, \Omega, \Lambda, L$, and $L_z$, of which $\Lambda$ is conserved through the orbit-averaged dynamics. Of course, the evolution equations for the four other coordinates can be obtained from Hamilton's canonical equations of motion, which involves obtaining the corresponding derivatives of $\mathbf{r}(\omega, \Omega, L, L_z)$ at fixed $(M, \Lambda)$. These calculations are rather involved and can become degenerate, e.g., at $I=0$ (equatorial orbits), $e=0$ (circular orbits), or $e=1$ (radial orbits). An equivalent alternative is to work directly with the forces acting on the wires, as we pursue below.

Rather than integrating the motion with respect to the orbital elements, we keep track of the wire dynamics through the

dimensionless vectors

$$\mathbf{h} = \frac{\mathbf{r} \times \mathbf{p}}{\Lambda}, \quad \mathbf{e} = \frac{\mathbf{p} \times (\mathbf{r} \times \mathbf{p})}{m^2 GM_{\bullet}} - \hat{\mathbf{r}}, \quad (11)$$

where a hat denotes a unit vector, as usual. Here, $\mathbf{h}$ is the angular momentum scaled to the circular angular momentum at the same energy, and $\mathbf{e}$ is the eccentricity vector, which points in the direction of the pericenter and has the magnitude $|\mathbf{e}| = e$.

While only the four orbital elements $\omega, \Omega, L$, and $L_z$ evolve, the two vectors $\mathbf{h}$ and $\mathbf{e}$ have six dynamical variables in total. The two associated degeneracies are captured by the two identities

$$\mathbf{h} \cdot \mathbf{e} = 0, \quad \mathbf{h}^2 + \mathbf{e}^2 = 1. \quad (12)$$

In terms of these vectors, position and momentum are

$$\mathbf{r} = a(\cos E - e)\hat{\mathbf{e}} + a \sin E \mathbf{h} \times \hat{\mathbf{e}}, \quad (13a)$$

$$\mathbf{p} = -\frac{\Lambda}{r} \sin E \hat{\mathbf{e}} + \frac{\Lambda}{r} \cos E \mathbf{h} \times \hat{\mathbf{e}}. \quad (13b)$$

### 3.1. Orbit-averaged Rates of Change

Using Equation (11), one can compute the time derivatives of the vectors $\mathbf{h}$ and $\mathbf{e}$ from the derivatives of the canonical variables $\mathbf{r}$ and $\mathbf{p}$. The combined force from all other wires at position $\mathbf{r}_i$ on wire $i$ is

$$\mathcal{F}_i(\mathbf{r}_i) = - \frac{\partial \langle H_{\star} \rangle_{\setminus i}}{\partial \mathbf{r}_i} = \sum_{j \neq i} Gm_i m_j \langle \mathbf{F}(\mathbf{r}_i - \mathbf{r}_j) \rangle_{\setminus i}, \quad (14)$$

where $\langle \cdot \rangle_{\setminus i}$ denotes an average over all orbits except for $i$ and $\mathbf{F}(\mathbf{r}) = -\partial\varphi/\partial\mathbf{r}$ with $\varphi(\mathbf{r}) = -|\mathbf{r}|^{-1}$. At each value of the mean anomaly $M_i$ along the orbit, the force $\mathcal{F}_i$ induces the local changes$^8$

$$\dot{\mathbf{h}}_{i,\text{local}}(M_i) = \frac{1}{\Lambda_i} \mathbf{r}_i \times \mathcal{F}_i - \frac{\dot{\Lambda}_i}{\Lambda_i} \mathbf{h}_i, \quad (15a)$$

$$\dot{\mathbf{e}}_{i,\text{local}}(M_i) = \frac{a_i}{\Lambda_i^2} [\mathcal{F}_i \times (\mathbf{r}_i \times \mathbf{p}_i) + \mathbf{p}_i \times (\mathbf{r}_i \times \mathcal{F}_i)]. \quad (15b)$$

In order to obtain the rates of change $\dot{\mathbf{h}}_i$ and $\dot{\mathbf{e}}_i$, these local changes are averaged over the unperturbed Keplerian orbit $i$. The last term in Equation (15a) averages to zero for conservative forces, and we obtain

$$\dot{\mathbf{h}}_i = \frac{1}{\Lambda_i} \sum_{j \neq i} Gm_j m_i \langle \mathbf{r}_i \times \mathbf{F}_{ij} \rangle, \quad (16a)$$

$$\dot{\mathbf{e}}_i = \frac{a_i}{\Lambda_i^2} \sum_{j \neq i} Gm_j m_i \langle \mathbf{F}_{ij} \times (\mathbf{r}_i \times \mathbf{p}_i) + \mathbf{p}_i \times (\mathbf{r}_i \times \mathbf{F}_{ij}) \rangle. \quad (16b)$$

Here, we used the shorthand $\mathbf{F}_{ij} \equiv \mathbf{F}(\mathbf{r}_i - \mathbf{r}_j)$ and the averages are over all orbits including $i$.

### 3.2. Keeping all Conservation Laws

Expressions (16) satisfy constraints (12) for any conservative field $\mathbf{F}_{ij}$ (Touma et al. 2009). Gauß showed that the average over orbit $j$ could be expressed in closed form. However, the double average over two interacting Keplerian wires $i$ and $j$ cannot be expressed in closed form, but rather requires

<footnote>$^8$ The term $\hat{\mathbf{r}}$ in Equation (11) does not contribute to the dynamics.</footnote>
numerical treatment. The resulting numerical average will inevitably carry a small error, with the consequence that constraints (12) are no longer exactly honoured. While in general, small numerical errors in conserved quantities are not a serious problem, this particular situation is awkward because conditions (12) are essential for the interpretation of $\mathbf{h}_i$ and $\mathbf{e}_i$ as a description of the wires. It is therefore important to construct numerical expressions for $\dot{\mathbf{h}}_i$ and $\dot{\mathbf{e}}_i$ that, despite their discretization errors, keep constraints (12) valid to machine precision.

In order to achieve this goal, we note that Equation (16) is equivalent to Milankovitch's (Milankovitch 1939) relations (Tremaine et al. 2009; Rosengren & Scheeres 2014),

$$ \dot{\mathbf{h}}_i = -\frac{1}{\Lambda_i} \left( \mathbf{h}_i \times \frac{\partial \langle H_* \rangle}{\partial \mathbf{h}_i} + \mathbf{e}_i \times \frac{\partial \langle H_* \rangle}{\partial \mathbf{e}_i} \right), \quad (17a) $$

$$ \dot{\mathbf{e}}_i = -\frac{1}{\Lambda_i} \left( \mathbf{h}_i \times \frac{\partial \langle H_* \rangle}{\partial \mathbf{e}_i} + \mathbf{e}_i \times \frac{\partial \langle H_* \rangle}{\partial \mathbf{h}_i} \right). \quad (17b) $$

Here, the averaged Hamiltonian $\langle H_* \rangle$ is expressed as a function of the vectors $\mathbf{h}_i$, $\mathbf{e}_i$ (and the constants $\Lambda_i$) from all wires.

These relations are equivalent to the canonical equations of motion, but they are much more useful. First, the vectors $\mathbf{h}$ and $\mathbf{e}$ are well defined for all orbits, even those with zero eccentricity, unit eccentricity, and zero inclination. Second, they have the beautiful property that identities (12) and the conservation of total energy are explicitly satisfied for any form of $\langle H_* \rangle$, i.e.,

$$ \frac{1}{2} \frac{d(\mathbf{h}_i^2 + \mathbf{e}_i^2)}{dt} = \dot{\mathbf{h}}_i \cdot \mathbf{h}_i + \dot{\mathbf{e}}_i \cdot \mathbf{e}_i = 0, \quad (18a) $$

$$ \frac{d(\mathbf{h}_i \cdot \mathbf{e}_i)}{dt} = \dot{\mathbf{h}}_i \cdot \mathbf{e}_i + \dot{\mathbf{e}}_i \cdot \mathbf{h}_i = 0, \quad (18b) $$

$$ \frac{d\langle H_* \rangle}{dt} = \sum_i \frac{\partial \langle H_* \rangle}{\partial \mathbf{h}_i} \cdot \dot{\mathbf{h}}_i + \frac{\partial \langle H_* \rangle}{\partial \mathbf{e}_i} \cdot \dot{\mathbf{e}}_i = 0. \quad (18c) $$

These properties suggest an alternative, fully conservative approach to numerically computing the rates $\dot{\mathbf{h}}_i$ and $\dot{\mathbf{e}}_i$: instead of following established practice of discretizing the orbit averages in Equation (16), we first discretize the orbit average for $\langle H_* \rangle$ and subsequently use Milankovitch's relations (17) to obtain $\dot{\mathbf{h}}_i$ and $\dot{\mathbf{e}}_i$,<fnref n="9" /> as detailed in Section 4.1.

### 3.3. Reformulation as a Spin System

Following Klein (1924), we introduce the vectors

$$ \mathbf{b}_{i+} = \mathbf{h}_i + \mathbf{e}_i \quad \text{and} \quad \mathbf{b}_{i-} = \mathbf{h}_i - \mathbf{e}_i; \quad (19) $$

then identities (12) become

$$ |\mathbf{b}_{i\pm}| = 1. \quad (20) $$

Thus, the system of $N$ wires is fully described by $2N$ independent unit vectors $\mathbf{b} = \mathbf{b}_{i\pm}$, in other words, it is a classical spin system. The dynamics of these vectors is simply

$$ \dot{\mathbf{b}} = \dot{\mathbf{h}} \pm \dot{\mathbf{e}}. \quad (21) $$

Since the vectors $\mathbf{b}$ remain on the unit sphere, this equation can be expressed as a precession,

$$ \dot{\mathbf{b}} = \Omega \times \mathbf{b}. \quad (22) $$

The precession vectors $\Omega$ are not uniquely determined through Equation (21) because a component parallel to $\mathbf{b}$ does not affect $\dot{\mathbf{b}}$ (one would need to know $\ddot{\mathbf{b}}$ to construct such a component). The most conservative choice for the precession vectors—in the sense that $|\Omega|^2$ is minimized—is therefore

$$ \Omega = \mathbf{b} \times \dot{\mathbf{b}}, \quad (23) $$

for which $\Omega \cdot \mathbf{b} = 0$, so that $\mathbf{b}$ moves along a great circle for constant $\Omega$.

Milankovitch's Equation (17), when re-expressed in terms of the vectors $\mathbf{b}$, reads

$$ \dot{\mathbf{b}} = \frac{2}{\Lambda} \frac{\partial \langle H_* \rangle}{\partial \mathbf{b}} \times \mathbf{b}. \quad (24) $$

These have the standard form for the equations of motion of classical spin systems. Compared with Equation (22), one might identify $2\Lambda^{-1}\partial\langle H_* \rangle/\partial\mathbf{b}$ with $\Omega$. However, the gradient of the Hamiltonian is not unique because owing to constraints (20), $\langle H_* \rangle$ is only determined up to additive terms of the form $(\mathbf{b}^2 - 1)\tilde{H}$ with an arbitrary function $\tilde{H}$. Such terms alter the gradient, but have no effect on $\dot{\mathbf{b}}$, reflecting the gauge invariance of Equation (24).

We are now set to compute the rates of change for each wire (Section 4), and to perform the time-integration of their evolution (Section 5). Throughout the next sections, we will test our algorithm in three systems, namely a simple analytical “pair” Hamiltonian of two oscillating wires, Kozai-Lidov oscillations of $N=2$ stars (whose orbits do not overlap radially) around a central BH, and an $N=10^4$ stellar cluster mimicking Sgr A*. We refer to Appendices D, E, and F for the detailed description of these setups.

## 4. Calculating the Rates of Change

Our computed rates of change $\dot{\mathbf{h}}_i$ and $\dot{\mathbf{e}}_i$ differ from the continuous forms (16) or (17) first by estimating the averaging integrals by discrete sums (detailed in Section 4.1), and second by approximating the forces using spherical harmonics (detailed in Section 4.2). We also refer to Appendix B for explicit expressions of the contributions to $\dot{\mathbf{h}}$ and $\dot{\mathbf{e}}$ from $\langle H_{GR} \rangle$, namely the Schwarzschild and Lense-Thirring relativistic precessions.

### 4.1. Discretized Rates of Change

In order to numerically calculate the orbit averages, we approximate them by discrete midpoint sums over $K$ positions along each orbit,<fnref n="10" />

$$ m_i \oint \frac{dM_i}{2\pi} f(M_i) \rightarrow \sum_{k=1}^{K} \mu_{ik} f(M_{ik}), \quad (25) $$

for any function $f(M)$. In other words, each eccentric wire of mass $m_i$ is replaced by $K$ nodes of masses $\mu_{ik}$. Placing the nodes' mean anomalies $M_{ik}$ uniformly (Gürkan &

<footnote n="9">⁹ This method for deriving discretized equations of motion is analogous to conservative formulations of smoothed particle hydrodynamics (SPH), where instead of discretizing the Euler equation directly, one discretizes the fluid Lagrangian, such that the Euler-Lagrange equations of motion give rise to a discretized form of the Euler equation that honors all conservation laws (Springel & Hernquist 2002).</footnote>

<footnote n="10">¹⁰ Adopting a different value of $K$ for each wire, e.g., depending on eccentricity, is a straightforward extension.</footnote>
Hopman 2007) is a bad idea, not only because it requires solving Kepler’s Equation (5) for each node, but also because it only poorly samples the pericentric passage, in particular for high eccentricity, and hence produces a slow convergence with K. Instead, we place the nodes uniformly in eccentric anomaly,

$$M_{ik} = E_k - e_i \sin E_k, \quad E_k = \left(k - \frac{1}{2}\right) \Delta E \quad (26)$$

with $\Delta E = 2\pi/K$. With this choice, the nodes are placed symmetrically with respect to the pericenter, but there is no node at exactly this position. To complete the discretization (25), we specify the node masses,

$$\mu_{ik} = \frac{m_i}{2\pi} \frac{dM_i}{dE} \Delta E = \frac{m_i}{K} (1 - e_i \cos E_k). \quad (27)$$

With these specifications, the discretized potential energy of node k on wire i becomes

$$\phi_{ik} = \mu_{ik} \psi_{ik}, \quad \psi_{ik} = \sum_{j \neq i}^{N} \sum_{l=1}^{K} G \mu_{jl} \varphi(\mathbf{r}_{ik} - \mathbf{r}_{jl}), \quad (28)$$

such that the total interaction energy between the wires is computed as

$$\langle H_{\star} \rangle = \frac{1}{2} \sum_{i=1}^{N} \sum_{k=1}^{K} \phi_{ik}. \quad (29)$$

This expression depends on $\mathbf{h}_i$ and $\mathbf{e}_i$ through the node positions $\mathbf{r}_{ik}$, but also through their masses $\mu_{ik}$. The latter depend on the eccentricities $e_i$ because of our sampling of the nodes in eccentric rather than mean anomaly. The derivatives of $\langle H_{\star} \rangle$ resulting from this dependence through $\mu_{ik}$ induce via Milankovitch’s relations (17) the rates of change (see Appendix A)

$$\dot{\mathbf{h}}_i|_{\mu} = 0, \quad (30a)$$

$$\dot{\mathbf{e}}_i|_{\mu} = \frac{1}{\Lambda_i} \mathbf{h}_i \times \dot{\mathbf{e}}_i \sum_{k=1}^{K} \bar{\mu}_i \psi_{ik} \cos E_k, \quad (30b)$$

where $\bar{\mu}_i \equiv m_i/K$ is the mean node mass on the wire. These rates of change do not appear in the traditional Equation (16) based on forces.

The rates of change induced by the dependence of $\langle H_{\star} \rangle$ on $\mathbf{h}_i$ and $\mathbf{e}_i$ through $\mathbf{r}_{ik}$ are (see Appendix A)

$$\dot{\mathbf{h}}_i|_{\mathbf{r}} = \frac{a_i}{\Lambda_i} [\hat{\mathbf{e}}_i \times (\mathbf{C}_i - e_i \mathbf{F}_i) + (\mathbf{h}_i \times \hat{\mathbf{e}}_i) \times \mathbf{S}_i], \quad (31a)$$

$$\dot{\mathbf{e}}_i|_{\mathbf{r}} = \frac{a_i}{\Lambda_i} \left[ \mathbf{F}_i \times \mathbf{h}_i - e_i \mathbf{S}_i + \frac{1}{e_i} [\mathbf{S}_i \cdot \hat{\mathbf{e}}_i - \mathbf{C}_i \cdot (\mathbf{h}_i \times \hat{\mathbf{e}}_i)] \hat{\mathbf{e}}_i \right] \quad (31b)$$

with $\mathbf{F}_i \equiv \sum_j \mathbf{d}_{ik}$, $\mathbf{C}_i \equiv \sum_k \mathbf{f}_{ik} \cos E_k$, and $\mathbf{S}_i \equiv \sum_k \mathbf{f}_{ik} \sin E_k$. Here,

$$\mathbf{f}_{ik} = -\frac{\partial \phi_{ik}}{\partial \mathbf{r}_{ik}} = \sum_{j \neq i}^{N} \sum_{l=1}^{K} G \mu_{ik} \mu_{jl} \mathbf{F}(\mathbf{r}_{ik} - \mathbf{r}_{jl}) \quad (32)$$

is the gravitational force acting on node k of wire i and generated by the nodes of all other wires. Upon close

inspection,

$$\dot{\mathbf{h}}_i|_{\mathbf{r}} = \frac{1}{\Lambda_i} \sum_{k=1}^{K} \mathbf{r}_{ik} \times \mathbf{f}_{ik}, \quad (33)$$

which is identical to the discretization of Equation (16a) and implies conservation of total angular momentum $\mathbf{L}_{\text{tot}} = \sum_i \Lambda_i \mathbf{h}_i$, as expected from Noether’s theorem (because our discretization of the wire averages remains invariant under spatial rotations). Also note that the discretization of Equation (16b) clearly differs from Equation (31b). The advantage of Equation (31b) is that it satisfies $\dot{\mathbf{h}}_i \cdot \mathbf{h}_i + \dot{\mathbf{e}}_i \cdot \mathbf{e}_i = 0$ and $\dot{\mathbf{h}}_i \cdot \mathbf{e}_i + \dot{\mathbf{e}}_i \cdot \mathbf{h}_i = 0$ exactly, as required.

We note from Equations (30) and (31) that the computation of $\dot{\mathbf{h}}_i$ and $\dot{\mathbf{e}}_i$ due to interactions with N wires requires $O(N^2K^2)$ computations. As shown below, this can be reduced to $O(NKl_{\text{max}}^2)$ when approximating $\mathbf{F}(\mathbf{r}_{ik} - \mathbf{r}_{jl})$ by its expansion in spherical harmonics up to order $l_{\text{max}}$. To this end, it is advantageous to split $\mathbf{f}_{ik} = \mathbf{f}_{ik}^{\text{all}} - \mathbf{f}_{ik}^{\text{self}}$ into the contributions

$$\mathbf{f}_{ik}^{\text{all}} = \sum_{\{j,l\} \neq \{i,k\}}^{N,K} G \mu_{ik} \mu_{jl} \mathbf{F}(\mathbf{r}_{ik} - \mathbf{r}_{jl}), \quad (34a)$$

$$\mathbf{f}_{ik}^{\text{self}} = \sum_{l \neq k}^{K} G \mu_{ik} \mu_{il} \mathbf{F}(\mathbf{r}_{ik} - \mathbf{r}_{il}), \quad (34b)$$

due to all of the other nodes on all of the wires, and due to all of the other nodes on the same wire, respectively.

We stress that our algorithms for $\mathbf{f}_{ik}^{\text{all}}$ and $\mathbf{f}_{ik}^{\text{self}}$ must compute the exact same wire self-gravity, such that when computing the difference $\mathbf{f}_{ik} = \mathbf{f}_{ik}^{\text{all}} - \mathbf{f}_{ik}^{\text{self}}$, these erroneous nonphysical contributions cancel exactly. This is important because the wire self-gravity can be substantial—it diverges logarithmically with the number K of nodes in the case of exact (unsoftened) Newtonian gravity. Similarly, our algorithms must also ensure that a single Keplerian wire is a stationary configuration of the orbit-averaged dynamics.

## 4.2. Multipole Expansion

In order to accelerate the computation of the $O(N^2K^2)$ pairwise node interactions, we use the expansion of the Newtonian interaction kernel $\varphi(\mathbf{r}) = -|\mathbf{r}|^{-1}$ in spherical harmonics,

$$\varphi(\mathbf{r}_i - \mathbf{r}_j) = - \sum_{\ell=0}^{\infty} \sum_{m=-\ell}^{\ell} \begin{cases} U_{\ell}^{m}(\mathbf{r}_i) T_{\ell}^{m}(\mathbf{r}_j) & \text{if } \mathbf{r}_i < \mathbf{r}_j, \\ T_{\ell}^{m}(\mathbf{r}_i) U_{\ell}^{m}(\mathbf{r}_j) & \text{if } \mathbf{r}_i > \mathbf{r}_j. \end{cases} \quad (35)$$

Here, we have defined real-valued upper and lower solid spherical harmonics, $U_{\ell}^{m}(\mathbf{r})$ and $T_{\ell}^{m}(\mathbf{r})$, respectively, which are described explicitly in Appendix C. Note that there exist essentially two flavors of codes based on spherical harmonic expansions: (i) codes in which the radial forces are evaluated using a basis-function expansion (see, e.g., Hernquist & Ostriker 1992; Saha 1993), and (ii) codes in which the radial forces associated with a given spherical harmonic are evaluated exactly (see, e.g., Villumsen 1982) using formulas involving the usual factor $r_{\text{min}}^{\ell}/r_{\text{max}}^{\ell+1}$, as in Equation (35), here absorbed into the definitions of $U_{\ell}^{m}$ and $T_{\ell}^{m}$. Our goal here is to follow this second avenue and tailor it to the case of Keplerian wires. Expansion (35) is isotropic, even when truncated at order $l_{\text{max}}$, and hence, by virtue of Noether’s theorem, will conserve the
total angular momentum $L_{tot}$ of the system. However, the expansion is not invariant under translations, and therefore the total linear momentum is usually not conserved in spherical-harmonic codes (Binney & Tremaine 2008, §2.9.4). Fortunately, as already highlighted in Equation (9), Keplerian wires do not induce any force on the central object, which by design ensures the conservation of the total linear momentum of the system.

Inserting expansion (35) into Equation (34a), we have

$$f_{ik}^{\text{all}} = \sum_{j,l} \sum_{\ell,m} G \mu_{ik} \mu_{jl} \begin{cases} U_{\ell}^{m}(\mathbf{r}_{jl}) \nabla T_{\ell}^{m}(\mathbf{r}_{ik}) & \text{if } r_{jl} < r_{ik} \\ T_{\ell}^{m}(\mathbf{r}_{jl}) \nabla U_{\ell}^{m}(\mathbf{r}_{ik}) & \text{if } r_{jl} > r_{ik} \end{cases}, \quad (36)$$

with an analogous expression for $f_{ik}^{\text{self}}$. The important difference between this expression and Equation (34a) is that the dependence on positions $\mathbf{r}_{ik}$ and $\mathbf{r}_{jl}$ has been factorized. Hence, the gradient terms depending on $\mathbf{r}_{ik}$ can be taken outside the sum over nodes after splitting it into an inner and outer part, giving

$$f_{ik}^{\text{all}} = \sum_{\ell,m} \mu_{ik} P_{\ell}^{m}(\mathbf{r}_{ik}) \nabla T_{\ell}^{m}(\mathbf{r}_{ik}) + \sum_{\ell,m} \mu_{ik} Q_{\ell}^{m}(\mathbf{r}_{ik}) \nabla U_{\ell}^{m}(\mathbf{r}_{ik}), \quad (37)$$

where

$$P_{\ell}^{m}(\mathbf{r}) = \sum_{j,l: r_{jl} < r} G \mu_{jl} U_{\ell}^{m}(\mathbf{r}_{jl}), \quad (38a)$$

$$Q_{\ell}^{m}(\mathbf{r}) = \sum_{j,l: r_{jl} > r} G \mu_{jl} T_{\ell}^{m}(\mathbf{r}_{jl}) \quad (38b)$$

are the multipoles of the distribution of all nodes inside and outside of radius $r$, respectively.

We note that the expansion from Equation (35) may also be used to define a cluster's truncated total energy, $E_{tot}(\ell_{max}, K)$ via Equation (10), which can be computed efficiently using the same algorithm as for the wire forces. In the limit $\ell_{max} \to +\infty$, i.e., in the absence of any multipole expansion, Equation (10) provides us with the full total energy, $E_{tot}^{\text{direct}}$, which we can also compute via a direct $O(N^2K^2)$ sum over all pairwise node-node couplings.

### 4.3. Computing the Wire Forces

We begin by sorting the nodes into ascending order of radius.<fnref n="11" /> Because the nodes on each wire are already sorted when created, this costs only $O(NK \ln N)$ rather than $O(NK \ln(NK))$ operations and is completely subdominant in the total operation budget of our code.

For brevity, we denote the pair of indices $i$ labeling the wire and $k$ labeling the node on the wire by $\alpha = \{i, k\}$, which we call the node index. Once the nodes are sorted in radius, the multipoles $P_{\ell}^{m}$ and $Q_{\ell}^{m}$ can be calculated by increasing and decreasing recurrence, respectively. The algorithm for the computation of $f_{\alpha}^{\text{all}}$ for all nodes is summarized in Table 1.

The total operation count of this algorithm is $O(NK\ell_{max}^2)$, with $\ell_{max}$ the maximum harmonic order considered. The same algorithm, but restricted to the nodes from only one wire,

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Illustration of the relative errors in $\dot{h}$ and $\dot{e}$ as a function of the number of nodes $K$ for fixed multipole truncation order $\ell_{max}$ for Kozai-Lidov oscillations (Appendix E) and an $N = 10^4$ stellar cluster (Appendix F). We used 500 independent realizations of the Kozai-Lidov system in the top panel and a single cluster for the bottom panel. The colored regions correspond to the 16% and 84% levels. We compute the relative errors by comparison to calculations with $K = 100$ for the top panel and $K = 20,480$ for the bottom panel.</figcaption></figure>

<table><caption>**Table 1**
The Algorithm for Computing $f_{\alpha}^{\text{all}}$ for All NK Nodes</caption><tr><td>1</td><td>sort all nodes in radius $r_{\alpha}$</td></tr><tr><td>2</td><td>set $f_{\alpha}^{\text{all}} = 0$ for all nodes $\alpha$</td></tr><tr><td>3</td><td>set $P_{\ell}^{m} = 0$</td></tr><tr><td>4</td><td>for all nodes $\alpha$ in order of increasing radius $r_{\alpha}$:</td></tr><tr><td>4.1</td><td>$f_{\alpha}^{\text{all}} \leftarrow f_{\alpha}^{\text{all}} + \mu_{\alpha} \sum_{\ell,m} P_{\ell}^{m} \nabla T_{\ell}^{m}(\mathbf{r}_{\alpha})$</td></tr><tr><td>4.2</td><td>$P_{\ell}^{m} \leftarrow P_{\ell}^{m} + G \mu_{\alpha} U_{\ell}^{m}(\mathbf{r}_{\alpha})$</td></tr><tr><td>5</td><td>set $Q_{\ell}^{m} = 0$</td></tr><tr><td>6</td><td>for all nodes $\alpha$ in order of decreasing radius $r_{\alpha}$:</td></tr><tr><td>6.1</td><td>$f_{\alpha}^{\text{self}} \leftarrow f_{\alpha}^{\text{self}} + \mu_{\alpha} \sum_{\ell,m} Q_{\ell}^{m} \nabla U_{\ell}^{m}(\mathbf{r}_{\alpha})$</td></tr><tr><td>6.2</td><td>$Q_{\ell}^{m} \leftarrow Q_{\ell}^{m} + G \mu_{\alpha} T_{\ell}^{m}(\mathbf{r}_{\alpha})$</td></tr></table>

**Note.** The forces on node $\alpha$ due to all other nodes at smaller and large radii are computed in steps 4.1 and 6.1 from the inner and outer multipoles $P_{\ell}^{m}$ and $Q_{\ell}^{m}$, respectively, which in turn are accumulated in steps 4.2 and 6.2.

calculates $f_{\alpha}^{\text{self}}$ for that wire and again requires $O(NK\ell_{max}^2)$ operations for all wires.

### 4.4. Convergence of the Force Calculations

In Figure 1 we illustrate the dependence of the relative errors on both $\dot{h}$ and $\dot{e}$ as a function of the number of nodes per wire $K$.

For the Kozai-Lidov system, for a fixed value of $\ell_{max}$, we find that the error decreases exponentially as $e^{-K}$. This is a direct consequence of the absence of any radial overlap between the two stellar orbits at play. In this case, the double orbit-averaged integral from Equation (9) naturally splits into the product of two 1D integrals, which the multipole algorithm

<footnote n="11">¹¹ For the self-gravity to be correctly removed, it is mandatory for the radius sortings used to compute $f_{ik}^{\text{all}}$ and $f_{ik}^{\text{self}}$ to be consistent.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Illustration of the relative error in the truncated total energy, $E_{\text{tot}}(\ell_{\text{max}}, K)$, as a function of the number of nodes $K$ for fixed multipole truncation order $\ell_{\text{max}}$ with respect to calculations with $K_{\text{max}} = 100$ for the top panel and $K_{\text{max}} = 20,480$ for the bottom panel. In both panels, the colored regions correspond to the 16% and 84% levels over 500 independent realizations.</figcaption></figure>

from Equation (38) catches. Importantly, the integrands of each of these 1D integrals are $2\pi$-periodic with respect to the eccentric anomaly. In this case, the midpoint sampling from Equation (26) ensures an exponential convergence with respect to the number of nodes $K$ (Trefethen & Weideman 2014).

For the $N = 10^4$ clusters, given that the orbits exhibit a wide range of eccentricities—$0 \le e \le 0.99$, see Appendix F—radial overlaps are ubiquitous and the exponential convergence with respect to $K$ no longer holds. More precisely, for a fixed value of $\ell_{\text{max}}$, the error decreases as $K^{-2.5}$ for $\dot{h}$ and as $K^{-1.5}$ for $\dot{e}$. We believe that this difference in the rates of convergence stems from the fact that the discontinuity in Equation (36) is only in the radial force, which does not enter the equation for $\dot{h}$.

We reach similar conclusions in Figure 2, where we present the error in the truncated total energy, $E_{\text{tot}}(\ell_{\text{max}}, K)$, as a function of $K$. For the radially nonoverlapping Kozai–Lidov simulations, it converges exponentially with respect to $K$, while for the $N = 10^4$ clusters, it is found to converge like $K^{-2.5}$, likely a consequence of the absence of any discontinuities in the expansion from Equation (35).

In Figure 3 we investigate the relative errors in $\dot{h}$ and $\dot{e}$ as a function of $\ell_{\text{max}}$ for a fixed number of nodes $K$. For the Kozai–Lidov simulations, we find that the error decreases exponentially as $e^{-1.6 \ell_{\text{max}}}$. Such a rapid convergence stems from the absence of any radial overlap between the two stellar orbits. Indeed, glancing back at Equation (35), the radial dependence of a typical term of the Legendre expansion is of the form $(r_{\text{min}}/r_{\text{max}})^\ell \le \eta^\ell$, with $\eta < 1$ the distance ratio closest to unity that occurs as the inner and outer stars run through their orbits. For the present setup, one readily finds $\eta = 1/5$, see Appendix E. As a consequence, when truncated at order $\ell_{\text{max}}$, the typical error is expected to scale like $\eta^{\ell_{\text{max}}} \simeq e^{-1.6 \ell_{\text{max}}}$, matching the numerical measurement from Figure 3.

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Illustration of the relative error in $\dot{h}$ and $\dot{e}$ as a function of $\ell_{\text{max}}$ for a fixed number of nodes ($K = 100$ for Kozai–Lidov and $K = 20,480$ for the cluster), using the same conventions as in Figure 1. Relative errors are computed by comparison to force calculations with $\ell_{\text{max}} = 20$ for Kozai–Lidov and $\ell_{\text{max}} = 100$ for the cluster.</figcaption></figure>

In the same figure, we note that the convergence with respect to $\ell_{\text{max}}$ is much slower for the $N = 10^4$ clusters. This is a consequence of the very slow convergence of the Legendre expansion (35) in any regime where the wires overlap in radius. We also note that $\dot{e}$ is more accurate than $\dot{h}$. This is likely due to the analytical contribution from the Schwarzschild precession (see Equation (B2)), which only affects $\dot{e}$ and somewhat reduces the relative errors. Let us finally emphasize that the truncation of the multipole approach effectively softens the interaction potential, thereby avoiding the singularities that would otherwise occur when wires cross. As such, this method is well suited to studying scalar and vector-resonant relaxation (Rauch & Tremaine 1996), which are dominated by large-scale effects.

We investigate the same convergence in Figure 4 by considering the dependence with respect to $\ell_{\text{max}}$ of the relative error between the truncated total energy, $E_{\text{tot}}(\ell_{\text{max}}, K)$, and the full total energy, $E_{\text{tot}}^{\text{direct}}$. For the Kozai–Lidov simulations, we recover an exponential convergence with respect to $\ell_{\text{max}}$, as $e^{-1.6 \ell_{\text{max}}}$. For the $N = 10^4$ clusters, where numerous orbits are radially overlapping, we empirically find that $E_{\text{tot}}(\ell_{\text{max}}, K)$ converges roughly like $\ell_{\text{max}}^{-1.5}$ with respect to the multipole truncation order. This is compatible with the asymptotic scalings predicted in Appendix B5 of Kocsis & Tremaine (2015) in the context of vector resonant relaxation. Indeed, for radially overlapping noncoplanar orbits, which mostly compose the present $N = 10^4$ clusters, a given harmonic, $\ell$, is found to contribute to the total truncated energy scaling like $\ell^{-2.5}$. Once summing over all the harmonics $0 \le \ell \le \ell_{\text{max}}$, this ultimately leads to a relative error in $E_{\text{tot}}(\ell_{\text{max}}, K)$ scaling like $\ell_{\text{max}}^{-1.5}$.

Finally, in Figure 5 we illustrate the joint dependence of the relative error on the truncated total energy, $E_{\text{tot}}(\ell_{\text{max}}, K)$, with respect to the full total energy, $E_{\text{tot}}^{\text{direct}}$, as a function of $\ell_{\text{max}}$ and
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Illustration of the relative error in the truncated total energy, $E_{tot}(\ell_{max}, K)$, as a function of the multipole truncation order $\ell_{max}$ for a fixed number of nodes ($K=100$ for Kozai–Lidov and $K=20,480$ for the cluster), using the same convention as in Figure 2. The full total energy, $E_{tot}^{direct}$, was computed by direct summation of Equation (7) using $K=20,480$ for the Kozai–Lidov simulations and $K=1000$ for the clusters, respectively.</figcaption></figure>

$K$ for the cluster simulations. For the present system, we find that errors saturate for $K \simeq 2 \ell_{max}$. Because the computational complexity of the discretized harmonic expansion scales like $O(K\ell_{max}^2)$, for a fixed value of $\ell_{max}$, ensuring converged estimates of $E_{tot}(\ell_{max}, K)$ requires ($O\ell_{max}^3$) operations.

## 5. Time Integration

In our application, the state of the system is described by the $N$ sets $y = \{\mathbf{b}_+, \mathbf{b}_-\}_{i}$. Given a timestep $\tau$, an integration scheme is a mapping $y_n \mapsto y_{n+1}$, where subscripts denote different time slices. Standard Runge–Kutta integration of the three-dimensional vectors $\mathbf{b}$ generally leads to violations of the constraints $|\mathbf{b}|=1$. Instead, we present here integration schemes that by construction exactly comply with these constraints. For convenience, we define the velocity field $\mathbf{B}(y) = \dot{\mathbf{b}} \pm \dot{\mathbf{e}}$, so the equations of motions read $\dot{\mathbf{b}} = \mathbf{B}(y)$.

### 5.1. Explicit Lie Group Methods

For a fixed value of the precession vector $\Omega = \mathbf{b} \times \mathbf{B}$ (see Equation (23)) and an initial condition $\mathbf{b}_0$, Equation (22) can be integrated exactly for a duration $t$, via Rodrigues’ rotation formula, to the new location

$$ \begin{aligned} \mathbf{b}(t) &= \cos(\Omega t)\mathbf{b}_0 + \sin(\Omega t)\hat{\Omega} \times \mathbf{b}_0 \\ &\quad + [1 - \cos(\Omega t)](\mathbf{b}_0 \cdot \hat{\Omega})\hat{\Omega} \\ &= \phi[t\Omega] \circ \mathbf{b}_0 \end{aligned} \quad (39) $$

with $|\mathbf{b}(t)| = |\mathbf{b}_0|$. In order to exactly preserve the constraints $|\mathbf{b}|=1$, explicit integration schemes analogous to Runge–Kutta schemes can be devised via Lie group methods (see, e.g., Section IV.8 in Hairer et al. (2006)) by appropriately

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Illustration of the median relative error in the truncated total energy, $E_{tot}(\ell_{max}, K)$, with respect to the full total energy, $E_{tot}^{direct}$, as a function of $\ell_{max}$ and $K$ for the cluster simulations, using the same parameters as in Figure 4. Errors are empirically found to saturate for $K \simeq 2 \ell_{max}$ (black line), while the computational complexity scales like $O(K\ell_{max}^2)$ (gray lines).</figcaption></figure>

concatenating rotations using precession vectors obtained at various intermediate stages. In practice, we use Munthe-Kaas (1999) integrators for their simplicity. To do so, one rewrites the dynamics as

$$ \mathbf{b}(t) = \phi[U(t)] \circ \mathbf{b}_0, \quad (40) $$

which imposes that $U(t)$ evolves according to

$$ \dot{U} = d\phi_U^{-1} \circ \{\Omega(\phi[U] \circ \mathbf{b}_0)\}, \quad (41) $$

with $\Omega(\mathbf{b}) = \mathbf{b} \times \mathbf{B}(\mathbf{b})$ following Equation (23). Here, $d\phi_U^{-1}$ is the inverse of the differential of the rotation map (see Theorem 3 in Munthe-Kaas 1999). It generically reads

$$ d\phi_U^{-1} \circ \Omega = \sum_{k=0}^{+\infty} \frac{B_k}{k!} ad_U^k \circ \Omega, \quad (42) $$

with $B_k$ the Bernoulli numbers and $ad_U^k$ the $k$-th power of the adjoint operator. In the present context, this operator simply becomes

$$ ad_U^0 \circ \Omega = \Omega; \quad ad_U^k \circ \Omega = U \times (ad_U^{k-1} \circ \Omega). \quad (43) $$

A Munthe-Kaas scheme proceeds then by using a classical Runge–Kutta scheme applied to $U(t)$ and “correcting” the intermediate precession vectors with $d\phi_U^{-1}$ appropriately truncated (Theorem IV.8.5 in Hairer et al. 2006). For a second-order scheme, one can use $d\phi_U^{-1} \circ \Omega = \Omega$, from which one constructs the two-stage explicit midpoint rule (coined MK2). Starting from an initial state $\mathbf{b}_n$, it proceeds via

$$ \begin{aligned} \Omega_1 &= \Omega(\mathbf{b}_n), \\ \mathbf{b}_2 &= \phi\left[\frac{\tau}{2}\Omega_1\right] \circ \mathbf{b}_n, \\ \Omega_2 &= \Omega(\mathbf{b}_2), \\ \mathbf{b}_{n+1} &= \phi[\tau\Omega_2] \circ \mathbf{b}_n. \end{aligned} \quad (44) $$

Naturally, the rotations are performed over all elements of $y$ simultaneously.

When constructing higher-order methods, one must account for the fact that rotations (39) do not commute, i.e., better approximations of $d\phi_U^{-1}$ have to be used. Here, we implement a
four-stage scheme (coined MK4) based on the classical fourth-order Runge–Kutta scheme. It reads

$$
\begin{align*}
\tilde{\Omega}_1 &= \Omega(b_n), \\
b_2 &= \phi\left[\frac{\tau}{2}\tilde{\Omega}_1\right] \circ b_n, \\
\tilde{\Omega}_2 &= d\phi_{\frac{1}{2}\tilde{\Omega}_1}^{-1} \circ \Omega(b_2), \\
b_3 &= \phi\left[\frac{\tau}{2}\tilde{\Omega}_2\right] \circ b_n, \\
\tilde{\Omega}_3 &= d\phi_{\frac{1}{2}\tilde{\Omega}_2}^{-1} \circ \Omega(b_3), \\
b_4 &= \phi[\tau\tilde{\Omega}_3] \circ b_n, \\
\tilde{\Omega}_4 &= d\phi_{\tau\tilde{\Omega}_3}^{-1} \circ \Omega(b_4), \\
b_{n+1} &= \phi\left[\tau\left(\frac{1}{6}\tilde{\Omega}_1 + \frac{1}{3}\tilde{\Omega}_2 + \frac{1}{3}\tilde{\Omega}_3 + \frac{1}{6}\tilde{\Omega}_4\right)\right] \circ b_n, \tag{45}
\end{align*}
$$

where $d\phi_U^{-1}$ is truncated at second order in $U$ so that Equation (42) becomes

$$d\phi_U^{-1} \circ \Omega = \Omega - \frac{1}{2}U \times \Omega + \frac{1}{12}U \times (U \times \Omega). \quad (46)$$

Schemes (44) and (45) are (i) explicit, (ii) conserve $|b| = 1$ exactly,<fnref n="12" /> (iii) require two or four computations of the derivatives, respectively, and (iv) are second- and fourth-order accurate, respectively. Because they are direct translations of the usual Runge–Kutta methods, it is straightforward to design Munthe-Kaas schemes of higher order.

## 5.2. Symplectic Scheme

Following McLachlan et al. (2014), we also consider an integrator relying on the spherical midpoint method (coined MD2). This is based on the implicit relation

$$\frac{\mathbf{b}_{n+1} - \mathbf{b}_n}{\tau} = \mathbf{B} \left( \frac{\mathbf{b}_{n+1} + \mathbf{b}_n}{|\mathbf{b}_{n+1} + \mathbf{b}_n|} \right). \quad (47)$$

This relation can be solved via the fixed-point iteration

$$\mathbf{b}_{n+1}^{(0)} = \mathbf{b}_n; \quad \mathbf{b}_{n+1}^{(k+1)} = \mathbf{b}_n + \tau \mathbf{B} \left( \frac{\mathbf{b}_n + \mathbf{b}_{n+1}^{(k)}}{|\mathbf{b}_n + \mathbf{b}_{n+1}^{(k)}|} \right), \quad (48)$$

terminating when $||\mathbf{b}_{n+1}^{(k+1)} - \mathbf{b}_{n+1}^{(k)}||_\infty \le \epsilon$, for a given tolerance $\epsilon$. Naturally, these iterations are performed for all the elements of $\mathbf{b}$ simultaneously. In practice, we impose $\epsilon = 2.2 \times 10^{-15}$, which typically requires 5–12 iterations.

The scheme from Equation (47) is (i) second-order accurate, (ii) implicit, (iii) exactly conserves the constraints from Equation (20),<fnref n="13" /> (iv) symplectic (McLachlan et al. 2014), and (iv) conserves $L_{tot} = \sum_i \Lambda_i h_i$ because it is a global linear invariant. Because MD2 is a symmetric scheme, one can use symmetric composition techniques (see, e.g., Section V.3.2 in Hairer et al. 2006) to devise higher-order symplectic integrators.

The MK2, MK4, and MD2 schemes conserve the constraints for the vectors $\mathbf{h}$ and $\mathbf{e}$ to machine precision, which as far as we know, none of the previous studies (e.g., Touma et al. 2009;

<figure><img src="image_9.png" /><figcaption>Figure 6. Illustration in the Kozai–Lidov simulations of the evolution of the inner star’s eccentricity, $e_{in}$, as well as its inclination with respect to the outer star, $\hat{h}_{in} \cdot \hat{h}_{out}$. Full lines correspond to the multipole integration (with the MK4 scheme from Equation (45)), while dashed lines were obtained by direct integration using the IAS15 integrator (Rein & Spiegel 2015). See Appendix E for details.</figcaption></figure>

Tremaine et al. 2009; Hamers & Portegies Zwart 2016) have managed.

## 5.3. Convergence of the Time Integrations

As a first check of the sanity of the present algorithm, we compare in Figure 6 the time-evolution of the Kozai–Lidov system as predicted by the multipole approach and a direct integration of the associated three-body problem.

Both methods predict similar oscillations of the inner star.<fnref n="14" /> On average, the direct integration used an integration timestep of about $\tau \approx 10^{-4}$, while for the multipole integration, benefiting from its explicit orbit average, we used $\tau = 10$, i.e., a significantly larger timestep.

In order to further assess the performance of the integration schemes, we now consider the simple two-body pair Hamiltonian from Appendix D. In Figure 7 we illustrate the dependence of the relative error on $E_{tot}$ after a finite time as a function of the integration timestep $\tau$. We recover that MK2 and MD2 are second-order accurate with a finite-time error converging like $O(\tau^2)$, while MK4 is fourth-order accurate.

In Figure 8 we investigate for the same Hamiltonian the long-time trends for the errors in $E_{tot}$ and $L_{tot}$. As expected, the explicit schemes show errors growing like $t$. There are two main improvements in the symplectic scheme: (i) the error in $E_{tot}$ is bounded on long timescales, and (ii) the error in $L_{tot}$ only grows via the accumulation of round-off errors. Similarly, in Figure 9, we investigate the errors in $E_{tot}(\ell_{max}, K)$ and $L_{tot}$ for the $N = 10^4$ clusters. This figure exhibits similar asymptotic trends as Figure 8.

Finally, in Figure 10, we illustrate the error-cost relation in the cluster simulations. The cost of a given simulation is estimated via $\langle \#\text{forces} \rangle / (\tau/\tau_{dyn})$, with $\langle \#\text{forces} \rangle$ the average number of force evaluations per integration timestep, and $\tau_{dyn}$ the shortest dynamical time of the clusters at hand (see Equation (F3)). While the explicit schemes always perform a constant number of force evaluations per timestep, the symplectic scheme typically requires from 5 evaluations (smallest $\tau$) up to 12 (largest $\tau$). As one reduces $\tau$, the error in $E_{tot}(\ell_{max}, K)$ keeps shrinking, illustrating that our discretization scheme indeed corresponds to a Hamiltonian system

<footnote n="12">¹² We systematically perform the renormalization $\mathbf{b} \leftarrow \mathbf{b}/|\mathbf{b}|$ after every evaluation of Equation (39) to prevent a drift of $|\mathbf{b}|$ from round-off errors.</footnote>

<footnote n="13">¹³ The outcome of the fixed-point search in Equation (48) is systematically renormalized via $\mathbf{b} \leftarrow \mathbf{b}/|\mathbf{b}|$ to prevent a drift of $|\mathbf{b}|$ from round-off errors.</footnote>

<footnote n="14">¹⁴ When represented in Figure 6, the MK2, MK4, and MD2 integrations are indistinguishable.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Illustration of the relative error in the total energy $E_{\text{tot}}$ with respect to its initial value after a fixed time $T = 1.28$ as a function of the considered integration timestep, $\tau$. The colored regions correspond to the 16% and 84% levels over 500 independent realizations of the pair Hamiltonian. We compute the relative errors by comparison with the initial total energy. See Appendix D for details.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Illustration of the relative error in the total energy $E_{\text{tot}}$ (top) and total angular momentum $L_{\text{tot}}$ (bottom) with respect to their initial values as a function of the number of integration steps, $t/\tau$, using the same convention as in Figure 7. Here, we used $\tau = 10^{-2}$ as our integration timestep.</figcaption></figure>

(Equation (18c)). Finally, as already hinted in Figure 9, for the present rather short-duration integrations, at a given cost, the explicit schemes outperforms the symplectic one regarding the conservation of $E_{\text{tot}}(\ell_{\text{max}}, K)$, while the converse holds regarding $L_{\text{tot}}$.

## 6. Scalar Resonant Relaxation

Using the method described above, we may now use our simulations of an $N = 10^4$ cluster to model the dynamics around the supermassive BH SgrA*. We use these simulations to study scalar resonant relaxation (Rauch & Tremaine 1996), the process through which Keplerian wires can relax their eccentricities through secular interactions with one another. We

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Illustration of the relative error in the truncated total energy $E_{\text{tot}}(\ell_{\text{max}}, K)$ (top) and total angular momentum $L_{\text{tot}}$ (bottom) with respect to their initial values as a function of the number of integration steps, $t/\tau$, for the $N = 10^4$ clusters, using the same convention as in Figure 8. Here, we used $\ell_{\text{max}} = 10$, $K = 100$, and $\tau = 10^{-2} \tau_{\text{dyn}}$ as our integration timestep, see Appendix F.</figcaption></figure>

measure the associated diffusion coefficients, and compare them with the theoretical predictions put forward by Bar-Or & Fouvry (2018). All the details for these measurements and predictions are spelled out in Appendix F.

In Figure 11 we present the main result of this comparison, the (finite-time) diffusion coefficients $D_{hh} = \{[\Delta h(T)]^2\}/T$, where $\{\cdot\}$ stands for an ensemble average, and we used $T = 25$ kyr, $\ell_{\text{max}} = 10$, and $K = 100$. This figure shows a good agreement between the diffusion coefficients inferred from the numerical simulations and from kinetic theory. It is likely that the remaining mismatch stems from the difficulty of measuring/predicting secular diffusion coefficients on such a short finite time, $T$; see Equation (24) in Bar-Or & Fouvry (2018). Figure 11 further illustrates the sanity of the present algorithm, which reproduced the intricate resonant dynamical interactions of Keplerian wires in a galactic nucleus. Of course, the present application is only a first illustration of the physical mechanisms that can be investigated with our multipole method. Astrophysical applications to more diverse physical systems will be the subject of future work.

## 7. Discussion

First, we note that the computation of the rates of change $(\dot{h}, \dot{e})$ for all the nodes has a computational complexity in $O(NK\ell_{\text{max}}^2)$. This is the main benefit of the present method, which, at the cost of (truncated) harmonic expansions and discretized orbit averages, has a computational complexity that scales linearly in the total number of wires.

We also note that the computations of $\{\dot{h}, \dot{e}\}_{\text{all}}$ and $\{\dot{h}, \dot{e}\}_{\text{self}}$ are completely independent (and of similar computational difficulty), so that they can be easily be performed in parallel.
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Error-cost relation for the cluster simulations. Simulations were performed with (i) the symplectic scheme (MD2; stars) or the explicit ones (MK2 and MK4; circles), with (ii) different integration timesteps, $\tau/\tau_{dyn} = 2^k \times 10^{-3}$ ($0 \le k \le 7$; different colors), and (iii) up to different finite times, $T/\tau_{dyn} = 0.128 \times 2^{-k}$ ($0 \le k \le 5$); the longer T, the thicker the line. The cost of a given simulation is estimated via $\langle \#\text{forces} \rangle / (\tau/\tau_{dyn})$, with $\langle \#\text{forces} \rangle$ the average number of force evaluations per integration timestep, e.g., two for MK2. Median errors in $E_{\text{tot}}(\ell_{\text{max}}, K)$ (top) and $L_{\text{tot}}$ (bottom) are estimated over 128 realizations, with $\ell_{\text{max}} = 10$ and $K = 100$.</figcaption></figure>

Similarly, for each of these computations, the recurrences over $P_{\ell}^{m}$ and $Q_{\ell}^{m}$, i.e., steps 4 and 6 of the algorithm given in Table 1, are also independent from one another (and of similar computational difficulty), so that they can also be performed in parallel. As a consequence, one can easily benefit from a parallelization of the previous algorithm over four cores to compute $\{\hat{h}, \hat{e}\}_{\text{all}}$ and $\{\hat{h}, \hat{e}\}_{\text{self}}$ in parallel due to nodes at smaller as well as larger radii.

Further parallelization requires the parallel computation of the prefix sums from Equation (38) and is much more complex, but not impossible. We plan to investigate the application of the parallel scan algorithm by Ladner & Fischer (1980) in a future study.

Since the spherical harmonic expansion is isotropic, each wire's self-gravity, $f^{\text{self}}$, produces no torques and hence only affects each wire's apsidal precession rates, $\dot{e} \cdot \hat{k}/e$. In addition, these terms scale trivially with semimajor axis and hence depend nontrivially only on eccentricity. An alternative possibility to correct for these self-gravity contributions is therefore to use the scaling in *a* and a one-dimensional interpolation in *e* from a table computed in advance for the self-gravity-induced apsidal precession at a given *a*. In practice, this could essentially halve the computational cost. Since our present implementation shares this cost among 4 CPUs, a scheme like this would reduce this to 2 CPUs, but would not speed up the computation overall. However, in an

implementation featuring the parallel scan (see above), this would reduce the wall-clock time at a fixed number of CPUs.

A limitation of the present approach is the difficulty of using different timesteps for different wires. Indeed, in order to compute the prefix sums (38), it is essential to scroll through the entire population of nodes. Similarly, it would be of interest to investigate adapting the number of nodes *K* as a function of the orbital parameters of the wires, in particular, as a function of eccentricity. These aspects deserve further analyses.

Our algorithm has problems with nodes at identical or very similar radii. Indeed, the spherical harmonic expansion (35) fails to converge in the limit $r_i \to r_j$, even if $|r_i - r_j| > 0$, i.e., even if the nodes are spatially well separated. This problem can be solved by reverting to direct summation of the forces between nodes with near-identical radii: this only requires a slight modification of our algorithm. Such an approach also calls for force softening in order to avoid artificially large forces between nodes that are much closer to one another than to their respective wire-companions. This observation was already made in Figure 4 of Dehnen (2014). In essence, our present setup has a limiting opening angle $\theta = 1$, where the maximum errors hardly depend on $\ell_{\text{max}}$. For the errors to decrease significantly and hence offer a better efficiency at a given accuracy, one needs $\theta < 1$, which can be achieved by using direct summation over nodes with similar radii. We reserve the development and testing of this more general approach for a future publication. In the present work, we merely avoid situations with nodes at identical radii (except for nodes on the same wire, where it does not matter) by not having two wires with exactly the same (a,e).

The present multipole method differs from an alternative approach based on Gauß' method (Touma et al. 2009), which does not rely on a harmonic expansion. As a result, Gauß' method must unavoidably compute all the pairwise interactions between the wires, i.e., it has a computational complexity that scales quadratically with the number of wires. More precisely, its computational complexity scales like $O(N^2K_G)$, with $K_G$ the typical number of sampling points used to perform the second orbit average, which cannot be performed analytically. Gauß' method has several strengths: (i) it is very well suited to perform time integration of orbit-averaged problems with a small number of wires, (ii) because the computation of the orbit-averaged force between two wires is a numerically expensive calculation, i.e., it has a computational complexity in $O(K_G)$, the method can also significantly benefit from parallelization, (iii) because it relies on the computation of pairwise interactions, Gauß' method can also be further parallelized over pairs of stars, in contrast to the present multipole method, and (iv) Gauß' method is straightforward to apply in the case of softened interactions (Touma et al. 2009). One of the main drawbacks of Gauß' method is that its computational difficulty becomes prohibitively large as $N$ increases, an aspect that the multipole method alleviates by design.

In Figure 8 we noted that the accumulation of round-off errors in the symplectic scheme led to a biased growth of the error in $L_{\text{tot}}$, scaling like $t$. It would be worthwhile to improve the algorithm to always ensure an unbiased error growth in $\sqrt{t}$, following Brouwer's law (Brouwer 1937), while still complying (exactly) with the constraints from Equation (20). Similarly, one could investigate further the design of higher-order
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Illustration of the (finite-time) diffusion coefficients $D_{hh} = \{[\Delta h(T)]^2\}/T$, with $h = \sqrt{1 - e^2}$ in a galactic nucleus mimicking SgrA*. Wires are split into different bins of semimajor axes (a) in milliparsecs, represented with different colors. The dashed curves correspond to the theoretical predictions following the formalism of Bar-Or & Fouvry (2018). The thick curves correspond to the numerical measurements performed in the multipole simulations and their associated bootstrap errors. See Appendix F for details.</figcaption></figure>

structure-preserving integration schemes (see, e.g., Hairer et al. 2006).

## 8. Conclusion

In this paper, we described the use of a multipole expansion to perform orbit-averaged simulations of stellar dynamics in galactic nuclei. We also presented integration schemes complying exactly with these systems' orbital constraints. We used this method to recover the predicted statistical properties of the resonant relaxation of stellar eccentricities. This algorithm has two main advantages: it relies on an explicit orbit average over the fast Keplerian motion to reduce the range of dynamical times in the system, and through the multipole expansion, it offers a computational complexity that scales linearly with the total number of particles. The code will be made available on request.

The present work is only a first step toward the development of faster integration methods for the dynamics of galactic nuclei, here seen as an prototype of dynamically degenerate systems, i.e., systems exhibiting a global resonance condition in their orbital frequencies. As highlighted in the main text, there are three main limitations to the present method. First, even after the orbit average over the fast Keplerian motion induced by the central BH, the range of dynamical times in a typical nucleus still remains very significant. This is partly a consequence of the divergence of the relativistic in-plane precession frequency for very eccentric wires. Similarly, the system's dynamical range is also enlarged by the ranges in semimajor axes and individual masses. For these reasons, even while using the present method, long time-integrations still remain challenging. A second limitation of the present method is associated with its difficult parallelization, in particular the prefix sums, as well as with the difficulty of tailoring it to multi-timestep schemes. A third limitation stems from the current lack of any force softening and direct summations between nearby nodes. Such improvements are expected to increase the numerical accuracy and stability of the scheme. All these aspects deserve further investigation.

Finally, when augmented by a second orbit average over the in-plane precession, the present methods can also be used to perform fast simulations of the dynamics of massive annuli in

galactic nuclei, which undergo vector resonant relaxation (Kocsis & Tremaine 2015; Szölgyén & Kocsis 2018).

We thank the referee, Jihad Touma, for a very detailed and insightful report that led to major improvements in the paper. J. B.F. acknowledges support from Program number HST-HF2-51374, which was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. S.T. acknowledges the support of the Natural Sciences and Engineering Research Council of Canada (NSERC), funding reference number RGPIN-2020-03885. B.B. is supported by the Martin A. and Helen Chooljian Membership at the Institute for Advanced Study. This work is partially supported by grant Segal ANR-19-CE31-0017 of the French Agence Nationale de la Recherche, and by the Idex Sorbonne Université. This work has made use of the Infinity Cluster hosted by the Institut d'Astrophysique de Paris. We thank Stéphane Rouberol for running this cluster smoothly.

## Appendix A
### Derivation of the Discretized Equations of Motion

Here, we detail the derivation of Equations (30) and (31) from Milankovitch's relations (17) and the discretized Hamiltonian $\langle H_* \rangle$ as given in Equation (29).

We first consider the terms induced by the dependences of $\langle H_* \rangle$ on $h_i$ and $e_i$ through the node masses $\mu_{ik}$. From Equation (27),

$$ \frac{\partial \mu_{ik}}{\partial h_i} = 0, \quad \frac{\partial \mu_{ik}}{\partial e_i} = -\bar{\mu}_i \cos E_k \hat{e}_i. \qquad (A1) $$

Here, $\bar{\mu}_i \equiv m_i/K$, and therefore from Equation (29),

$$ \left. \frac{\partial \langle H_* \rangle}{\partial h_i} \right|_{\mu} = 0, \quad \left. \frac{\partial \langle H_* \rangle}{\partial e_i} \right|_{\mu} = - \hat{e}_i \sum_k \bar{\mu}_i \psi_{ik} \cos E_k. \qquad (A2) $$

Milankovitch's Equation (17) allows us then to obtain Equation (30).

Next, we derive the terms induced by the dependences of $\langle H_* \rangle$ on $h_i$ and $e_i$ through the node positions $r_{ik}$. For a fixed
vector $\mathbf{z}$, Equation (13a) gives us the derivative of the node positions via<fnref n="15" />

$$ 
\begin{aligned}
\frac{\partial(\mathbf{z} \cdot \mathbf{r}_{ik})}{\partial \mathbf{h}_i} &= a_i \sin E_k \hat{\mathbf{e}}_i \times \mathbf{z}, \\
\frac{\partial(\mathbf{z} \cdot \mathbf{r}_{ik})}{\partial \mathbf{e}_i} &= -a_i \mathbf{z} + \frac{a_i}{e_i} \cos E_k (\mathbf{z} - \mathbf{z} \cdot \hat{\mathbf{e}}_i \hat{\mathbf{e}}_i) \\
&\quad - \frac{a_i}{e_i} \sin E_k [(\mathbf{h}_i \times \mathbf{z}) - (\mathbf{h}_i \times \mathbf{z}) \cdot \hat{\mathbf{e}}_i \hat{\mathbf{e}}_i].
\end{aligned}
\quad (A3) $$

Differentiating Equation (29) and inserting these relations, we obtain

$$ 
\begin{aligned}
\left. \frac{\partial \langle H_* \rangle}{\partial \mathbf{h}_i} \right|_{\mathbf{r}} &= -a_i \hat{\mathbf{e}}_i \times \mathbf{S}_i, \\
\left. \frac{\partial \langle H_* \rangle}{\partial \mathbf{e}_i} \right|_{\mathbf{r}} &= a_i \mathbf{F}_i + \frac{a_i}{e_i} [(\mathbf{h}_i \times \mathbf{S}_i - \mathbf{C}_i) \\
&\quad - (\mathbf{h}_i \times \mathbf{S}_i - \mathbf{C}_i) \cdot \hat{\mathbf{e}}_i \hat{\mathbf{e}}_i],
\end{aligned}
\quad (A4) $$

with $\mathbf{F}_i \equiv \sum_k \mathbf{f}_{ik}$, $\mathbf{C}_i \equiv \sum_k \mathbf{f}_{ik} \cos E_k$, and $\mathbf{S}_i \equiv \sum_k \mathbf{f}_{ik} \sin E_k$, and $\mathbf{f}_{ik}$ is the gravitational force on node $k$ of wire $i$, defined by Equation (32). With these, Milankovitch's Equation (17) gives Equation (31).

## Appendix B
### The Relativistic Precessions

For the sake of completeness, in this appendix, we explicitly spell out the 1PN and 1.5PN relativistic precessions induced by the central BH. Following Equations (B4) and (B7) of Touma et al. (2009), the respective precession velocities are given by

$$ \dot{\mathbf{h}}_{\text{GR}} = \dot{\mathbf{h}}_{\text{GR}}^{\text{1PN}} + \dot{\mathbf{h}}_{\text{GR}}^{\text{1.5PN}}, \quad \dot{\mathbf{e}}_{\text{GR}} = \dot{\mathbf{e}}_{\text{GR}}^{\text{1PN}} + \dot{\mathbf{e}}_{\text{GR}}^{\text{1.5PN}}. \quad (\text{B1}) $$

Here, the 1PN contribution corresponds to the Schwarzschild precession that drives an in-plane (prograde) precession of the wire pericenter, without any change in the orbital plane of the wire. It reads

$$ \dot{\mathbf{h}}_{\text{GR}}^{\text{1PN}} = 0, \quad \dot{\mathbf{e}}_{\text{GR}}^{\text{1PN}} = \frac{3GM_*\Omega_{\text{Kep}}(a)}{c^2a(1-e^2)^{3/2}} \mathbf{h} \times \mathbf{e}, \quad (\text{B2}) $$

with $\Omega_{\text{Kep}}(a)$ the Keplerian frequency given by Equation (8). The 1.5PN precession is the Lense-Thirring precession. It is sourced by the central BH spin and leads to an orbital precession of the wire. It reads

$$ 
\begin{aligned}
\dot{\mathbf{h}}_{\text{GR}}^{\text{1.5PN}} &= \frac{2(GM_*)^{3/2}\Omega_{\text{Kep}}(a)}{c^3 a^{3/2} (1 - e^2)^{3/2}} \mathbf{S} \times \mathbf{h}, \\
\dot{\mathbf{e}}_{\text{GR}}^{\text{1.5PN}} &= \frac{2(GM_*)^{3/2}\Omega_{\text{Kep}}(a)}{c^3 a^{3/2} (1 - e^2)^{3/2}} \left[ \mathbf{S} - \frac{3\mathbf{h}(\mathbf{h} \cdot \mathbf{S})}{1 - e^2} \right] \times \mathbf{e},
\end{aligned}
\quad (\text{B3}) $$

with $|\mathbf{S}| \le 1$ the normalized spin angular momentum of the central BH.

## Appendix C
### Spherical Harmonics

In this appendix, we explicitly spell out our convention and our implementation for the calculation of the solid spherical harmonics, following definitions very similar to those of Dehnen (2014).

The (complex) surface spherical harmonics are defined as

$$ Y_{\ell}^{m}(\hat{\mathbf{r}}) = (-1)^{m} \sqrt{\frac{(\ell - m)!}{(\ell + m)!}} P_{\ell}^{m}(\cos(\theta)) e^{im\phi}, \quad (\text{C1}) $$

with $P_{\ell}^{m}$ the associated Legendre functions. The (complex) solid spherical harmonics are defined as

$$ 
\begin{aligned}
\Upsilon_{\ell}^{m}(\mathbf{r}) &= [(\ell - m)! (\ell + m)!]^{-1/2} r^{\ell} Y_{\ell}^{m}(\hat{\mathbf{r}}), \\
\Theta_{\ell}^{m}(\mathbf{r}) &= [(\ell - m)! (\ell + m)!]^{1/2} r^{-\ell-1} Y_{\ell}^{m}(\hat{\mathbf{r}}).
\end{aligned}
\quad (\text{C2}) $$

In order to work only with real-valued quantities, we define the associated real-valued solid harmonics as

$$ U_{\ell}^{m} = \begin{cases} \sqrt{2} \operatorname{Im}(\Upsilon_{\ell}^{|m|}) = \frac{1}{\sqrt{2}} i & [\Upsilon_{\ell}^{-m} - \Upsilon_{\ell}^{-m*}] \quad \text{if } m < 0, \\ \operatorname{Re}(\Upsilon_{\ell}^{|m|}) = \frac{1}{2} & [\Upsilon_{\ell}^{m} + \Upsilon_{\ell}^{m*}] \quad \text{if } m = 0, \\ \sqrt{2} \operatorname{Re}(\Upsilon_{\ell}^{|m|}) = \frac{1}{\sqrt{2}} & [\Upsilon_{\ell}^{m} + \Upsilon_{\ell}^{m*}] \quad \text{if } m > 0, \end{cases} \quad (\text{C3}) $$

and

$$ T_{\ell}^{m} = \begin{cases} \sqrt{2} \operatorname{Im}(\Theta_{\ell}^{|m|}) = \frac{1}{\sqrt{2}} i & [\Theta_{\ell}^{-m} - \Theta_{\ell}^{-m*}] \quad \text{if } m < 0, \\ \operatorname{Re}(\Theta_{\ell}^{|m|}) = \frac{1}{2} & [\Theta_{\ell}^{m} + \Theta_{\ell}^{m*}] \quad \text{if } m = 0, \\ \sqrt{2} \operatorname{Re}(\Theta_{\ell}^{|m|}) = \frac{1}{\sqrt{2}} & [\Theta_{\ell}^{m} + \Theta_{\ell}^{m*}] \quad \text{if } m > 0, \end{cases} \quad (\text{C4}) $$

where we note that this definition differs from the definition in Equations (58a) and (58b) of Dehnen (2014), where here we added a factor $\sqrt{2}$ to the expressions for $m \ne 0$. With the present convention, the Legendre expansion is unchanged when using the real-valued harmonics, i.e., one has for $r_1 < r_2$

$$ 
\begin{aligned}
\frac{1}{|\mathbf{r}_1 - \mathbf{r}_2|} &= \sum_{\ell,m} \Upsilon_{\ell}^{m}(\mathbf{r}_1) \Theta_{\ell}^{m*}(\mathbf{r}_2) \\
&= \sum_{\ell,m} U_{\ell}^{m}(\mathbf{r}_1) T_{\ell}^{m}(\mathbf{r}_2).
\end{aligned}
\quad (\text{C5}) $$

However, compared to Dehnen (2014), this prefactor will slightly affect some of the recurrence relations used to compute the real-valued harmonics and their gradients, as we now detail. The important property of all the coming recurrences is that for a given location $\mathbf{r}$, they allow for the computation of $U_{\ell}^{m}(\mathbf{r})$, $T_{\ell}^{m}(\mathbf{r})$ and their derivatives with complexity $O(\ell_{\max}^2)$.

<footnote n="15">¹⁵ Since $e = |\mathbf{e}| = \sqrt{1 - \mathbf{h}^2}$, the eccentricity dependence can equivalently be considered a dependence on $\mathbf{h}$. Indeed, Milankovitch's Equation (17) benefits from a gauge freedom in regard to the constraints from Equation (12) (see the discussion after Equation (24) or Appendix A in Tremaine et al. 2009), so that the final expressions for $\mathbf{h}$ and $\dot{\mathbf{e}}$ are not affected by this choice.</footnote>
### C.1. Computing $U_{\ell}^m$

The upper real-valued solid harmonics satisfy the initialization condition

$$U_0^0 = 1. \qquad (C6)$$

They satisfy some boundary recurrence relations for $m = \pm \ell$. For $\ell = 1$, it reads

$$\begin{aligned} U_1^{-1} &= \frac{1}{2\ell} [y \sqrt{2} U_0^0 + x \times 0], \\ U_1^1 &= \frac{1}{2\ell} [x \sqrt{2} U_0^0 - y \times 0], \end{aligned} \qquad (C7)$$

with $r = (x, y, z)$, while for $\ell \ge 2$, it reads

$$\begin{aligned} U_{\ell}^{-\ell} &= \frac{1}{2\ell} [y U_{\ell-1}^{\ell-1} + x U_{\ell-1}^{-(\ell-1)}], \\ U_{\ell}^{\ell} &= \frac{1}{2\ell} [x U_{\ell-1}^{\ell-1} - y U_{\ell-1}^{-(\ell-1)}]. \end{aligned} \qquad (C8)$$

Finally, inside the boundary of the $(\ell, m)$-domain, they satisfy the generic recurrence relation

$$U_{\ell}^{m} = \frac{2\ell - 1}{\ell^2 - m^2} z U_{\ell-1}^{m} - \frac{1}{\ell^2 - m^2} r^2 U_{\ell-2}^{m}, \qquad (C9)$$

which can also be applied to the cases $|m| = \ell - 1$, if one follows the convention $U_{\ell-2}^{\pm(\ell-1)} = 0$.

### C.2. Computing $T_{\ell}^m$

The lower real-valued solid harmonics satisfy the initialization condition

$$T_0^0 = \frac{1}{r}. \qquad (C10)$$

They satisfy some boundary recurrence relations for $m = \pm \ell$. For $\ell = 1$, it reads

$$\begin{aligned} T_1^{-1} &= (2\ell - 1) \left[ \frac{y}{r^2} \sqrt{2} T_0^0 + \frac{x}{r^2} \times 0 \right], \\ T_1^1 &= (2\ell - 1) \left[ \frac{x}{r^2} \sqrt{2} T_0^0 - \frac{y}{r^2} \times 0 \right], \end{aligned} \qquad (C11)$$

while for $\ell \ge 2$, it reads

$$\begin{aligned} T_{\ell}^{-\ell} &= (2\ell - 1) \left[ \frac{y}{r^2} T_{\ell-1}^{\ell-1} + \frac{x}{r^2} T_{\ell-1}^{-(\ell-1)} \right], \\ T_{\ell}^{\ell} &= (2\ell - 1) \left[ \frac{x}{r^2} T_{\ell-1}^{\ell-1} - \frac{y}{r^2} T_{\ell-1}^{-(\ell-1)} \right]. \end{aligned} \qquad (C12)$$

Finally, inside the boundary of the $(\ell, m)$-domain, they satisfy the generic recurrence relation

$$\begin{aligned} T_{\ell}^{m} &= (2\ell - 1) \frac{z}{r^2} T_{\ell-1}^{m} \\ &\quad - [(\ell - 1)^2 - m^2] \frac{1}{r^2} T_{\ell-2}^{m}, \end{aligned} \qquad (C13)$$

which can also be applied to the cases $|m| = \ell - 1$ if one follows the convention $T_{\ell-2}^{\pm(\ell-1)} = 0$.

### C.3. Computing $\nabla U_{\ell}^m$

The gradients of the upper spherical harmonics, $\nabla U_{\ell}^m$, can be determined using our previous computation of $U_{\ell}^m$. For $m = 0$, we have

$$\nabla U_{\ell}^{0} = \begin{pmatrix} -\frac{1}{\sqrt{2}} U_{\ell-1}^{1} \\ -\frac{1}{\sqrt{2}} U_{\ell-1}^{-1} \\ U_{\ell-1}^{0} \end{pmatrix}. \qquad (C14)$$

For $m = 1$, we use

$$\begin{aligned} \nabla U_{\ell}^{-1} &= \begin{pmatrix} -\frac{1}{2} U_{\ell-1}^{-2} \\ \frac{1}{2} [\sqrt{2} U_{\ell-1}^{0} + U_{\ell-1}^{2}] \\ U_{\ell-1}^{-1} \end{pmatrix}; \\ \nabla U_{\ell}^{1} &= \begin{pmatrix} \frac{1}{2} [\sqrt{2} U_{\ell-1}^{0} - U_{\ell-1}^{2}] \\ -\frac{1}{2} U_{\ell-1}^{-2} \\ U_{\ell-1}^{1} \end{pmatrix}, \end{aligned} \qquad (C15)$$

while for $m \ge 2$, this relation becomes

$$\begin{aligned} \nabla U_{\ell}^{-m} &= \begin{pmatrix} \frac{1}{2} [U_{\ell-1}^{-(m-1)} - U_{\ell-1}^{-(m+1)}] \\ \frac{1}{2} [U_{\ell-1}^{m-1} + U_{\ell-1}^{m+1}] \\ U_{\ell-1}^{-m} \end{pmatrix}; \\ \nabla U_{\ell}^{m} &= \begin{pmatrix} \frac{1}{2} [U_{\ell-1}^{m-1} - U_{\ell-1}^{m+1}] \\ \frac{1}{2} [-U_{\ell-1}^{-(m-1)} - U_{\ell-1}^{-(m+1)}] \\ U_{\ell-1}^{m} \end{pmatrix}. \end{aligned} \qquad (C16)$$

### C.4. Computing $\nabla T_{\ell}^m$

Similarly, the gradients of the lower spherical harmonics, $\nabla T_{\ell}^m$, are computed using the predetermined values of $T_{\ell}^m$. For $m = 0$, we have

$$\nabla T_{\ell}^{0} = \begin{pmatrix} -\frac{1}{\sqrt{2}} T_{\ell+1}^{1} \\ -\frac{1}{\sqrt{2}} T_{\ell+1}^{-1} \\ -T_{\ell+1}^{0} \end{pmatrix}. \qquad (C17)$$
For $m=1$, we use

$$ 
\begin{aligned}
\nabla T_{\ell}^{-1} &= \begin{pmatrix} -\frac{1}{2} T_{\ell+1}^{-2} \\ \frac{1}{2} [\sqrt{2} T_{\ell+1}^{0} + T_{\ell+1}^{2}] \\ -T_{\ell+1}^{-1} \end{pmatrix}; \\
\nabla T_{\ell}^{1} &= \begin{pmatrix} \frac{1}{2} [\sqrt{2} T_{\ell+1}^{0} - T_{\ell+1}^{2}] \\ -\frac{1}{2} T_{\ell+1}^{-2} \\ -T_{\ell+1}^{1} \end{pmatrix},
\end{aligned} 
\quad (C18) $$

while for $m \ge 2$, this relation becomes

$$ 
\begin{aligned}
\nabla T_{\ell}^{-m} &= \begin{pmatrix} \frac{1}{2} [T_{\ell+1}^{-(m-1)} - T_{\ell+1}^{-(m+1)}] \\ \frac{1}{2} [T_{\ell+1}^{m-1} + T_{\ell+1}^{m+1}] \\ -T_{\ell+1}^{-m} \end{pmatrix}; \\
\nabla T_{\ell}^{m} &= \begin{pmatrix} \frac{1}{2} [T_{\ell+1}^{m-1} - T_{\ell+1}^{m+1}] \\ \frac{1}{2} [-T_{\ell+1}^{-(m-1)} - T_{\ell+1}^{-(m+1)}] \\ -T_{\ell+1}^{m} \end{pmatrix}.
\end{aligned} 
\quad (C19) $$

## Appendix D
Pair Hamiltonian

In this appendix, we present a simple analytical two-body pair Hamiltonian used to assess the performance of our integration schemes on long timescales. The Hamiltonian is

$$ H = -a \mathbf{h}_1 \cdot \mathbf{h}_2 - b \mathbf{e}_1 \cdot \mathbf{e}_2, \quad (D1) $$

where $a, b$ are constant coefficients, and we work in dimensionless units $\Lambda_1 = \Lambda_2 = 1$. Following the Milankovitch Equation (17), the equations of motion for the two wires read

$$ 
\begin{aligned}
\dot{\mathbf{h}}_1 &= a \mathbf{h}_1 \times \mathbf{h}_2 + b \mathbf{e}_1 \times \mathbf{e}_2, \\
\dot{\mathbf{e}}_1 &= b \mathbf{h}_1 \times \mathbf{e}_2 + a \mathbf{e}_1 \times \mathbf{h}_2,
\end{aligned} 
\quad (D2a) $$

$$ 
\begin{aligned}
\dot{\mathbf{h}}_2 &= -a \mathbf{h}_1 \times \mathbf{h}_2 - b \mathbf{e}_1 \times \mathbf{e}_2, \\
\dot{\mathbf{e}}_2 &= -a \mathbf{h}_1 \times \mathbf{e}_2 - b \mathbf{e}_1 \times \mathbf{h}_2.
\end{aligned} 
\quad (D2b) $$

In Figures 7 and 8 we considered 500 independent realizations where (i) the coefficients $a, b$ were drawn uniformly in [0, 1], (ii) the initial orientations $\hat{\mathbf{h}}, \hat{\mathbf{e}}$ of the wires were drawn uniformly on the unit sphere, and (iii) the initial eccentricity, $e$, was drawn uniformly in [0, 1]. Because this pair Hamiltonian does not involve any numerical orbit average, it allows for a simple test of the integration schemes without any noise stemming from the harmonic truncation (via $\ell_{\max}$) or node sampling (via $K$).

## Appendix E
Kozai-Lidov Oscillations

In this appendix, we detail the parameters of the Kozai-Lidov simulations presented throughout the main text. We work in dimensionless units, setting in particular $G=1$. The mass of the central BH is taken to be $M_* = 10^6$. The stellar

cluster is composed of $N=2$ stars. The outer star is of mass $m_{\text{out}} = 10$, semimajor axis $a_{\text{out}} = 10$, and initial eccentricity $e_{\text{out}} = 0.5$, while for the inner star, we take $m_{\text{in}} = 1$, $a_{\text{in}} = 1$, and $e_{\text{in}} = 0.01$. At the initial time, we set the orientation of the outer orbit to $(\hat{\mathbf{h}}_{\text{out}}, \hat{\mathbf{e}}_{\text{out}}) = (\hat{\mathbf{e}}_z, \hat{\mathbf{e}}_y)$. Finally, for the inner star, we impose $\hat{\mathbf{h}}_{\text{in}} = \cos I_{\text{in}} \hat{\mathbf{e}}_z - \sin I_{\text{in}} \hat{\mathbf{e}}_y$ and $\hat{\mathbf{e}}_{\text{in}} = \cos \omega_{\text{in}} \hat{\mathbf{e}}_x + \sin \omega_{\text{in}} \hat{\mathbf{h}}_{\text{in}} \times \hat{\mathbf{e}}_x$, with the initial inclination $I_{\text{in}} = 60^\circ$ and pericenter phase $\omega_{\text{in}} = 90^\circ$. To perform the multipole integration, we use $K=100$ nodes, $\ell_{\max} = 10$ for the harmonic truncation, and $\tau = 10$ for the integration timestep. We did not account for any relativistic precessions. Finally, in Figures 1–4 to estimate typical dispersions, we performed 500 independent realizations of the system where we drew the initial inclinations and pericenter phase of the inner star uniformly within $\pm 5^\circ$ around their fiducial values, $I_{\text{in}}$ and $\omega_{\text{in}}$.

In order to check the validity of our multipole algorithm, we compare it with direct integrations of the associated three-body problem, as presented in Figure 6. These direct integrations were performed using the IAS15 integrator (Rein & Spiegel 2015), available through REBOUND (Rein & Liu 2012), which automatically adjusts the integration timestep. For the direct integrations, the initial mean anomalies of the two stars were picked at random. We considered a total of 100 different realizations, and found that $e_{\text{in}}$ and $\hat{\mathbf{h}}_{\text{in}} \cdot \hat{\mathbf{h}}_{\text{out}}$ only differed by $10^{-4}$ at most over these realizations. As such, in Figure 6, we could safely limit ourselves to only representing the median values obtained over the available realizations. Finally, we used the mappings from Equations (2b) and (11) to infer the orbital elements of the stars from their individual locations and velocities.

## Appendix F
A Typical Stellar Cluster

In this appendix, we briefly detail the parameters of the fiducial model considered throughout the numerical applications. Mimicking SgrA* (Gillessen et al. 2017), we take the mass of the central BH to be $M_* = 4 \times 10^6 M_{\odot}$. We assume that the initial eccentricity distribution of the stars follows a thermal distribution, $f_e \propto e$, within the range $e_{\min} \le e \le e_{\max}$, with $(e_{\min}, e_{\max}) = (0, 0.99)$. For the initial distribution in the semimajor axes of the stars, we assume that it follows a (truncated) power-law distribution of the form $f_e(E) \propto |E|^p$, with $E = E_{\text{Kep}}$ the Keplerian energy, and the power index $p < \frac{3}{2}$. With this choice, the number of stars per unit $a$ is given by

$$ N(a) = (3 - \gamma) \frac{N_0}{a_0} \left( \frac{a}{a_0} \right)^{2-\gamma}, \quad (\text{F1}) $$

with $\gamma = p + \frac{3}{2}$. For our fiducial model, we assume that the stellar distribution follows $\gamma = \frac{7}{4}$, the slope of the expected equilibrium (Bahcall & Wolf 1977). We also introduced $N_0 = g(\gamma)N(<a_0)$ with $g(\gamma) = 2^{-\gamma}\sqrt{\pi} \Gamma(1+\gamma)/\Gamma(\gamma-\frac{1}{2})$. Here, $N(<a_0)$ is the number of stars that physically lie within a radius $a_0$, which should not be confused with $N_0$, the number of stars with a semimajor axis smaller than $a_0$.

For the fiducial model, we assume that the total stellar mass that is physically enclosed within the radius of influence $r_h = 2$ pc is $M_*$. Assuming that all the stars have the same individual mass, $m_*$, this implies that we make the choice
$a_0 = r_h = 2$ pc and $N(<a_0) = M_*/m_\star$. It finally remains to carefully pick the value of $m_\star$. To do so, we fixed a range of interest in semimajor axis, namely $a_{min} \le a \le a_{max}$, with $(a_{min}, a_{max}) = (5, 100)$ in milliparsecs, to which the sampling of the initial semimajor axes of the particles is limited. The fiducial simulations are performed with a total number of stars, $N=10^4$. Following Equation (F1), this imposes a self-consistent relation between $N$ and $m_\star$, namely

$$m_{\star} = \frac{M_{\bullet}}{N} g(\gamma) \left[ \left( \frac{a_{\max}}{a_0} \right)^{3-\gamma} - \left( \frac{a_{\min}}{a_0} \right)^{3-\gamma} \right]. \quad (\text{F2})$$

This choice ensures that despite the limited range in semimajor axis, the stellar self-consistent potential remains close to the full expected one. For the fiducial model with $N=10^4$, we find $m_\star \simeq 8.63 M_\odot$.

In all these simulations, the integration timestep is dictated by the timescale associated with the (fast) in-plane relativistic precessions of the wires with small semimajor axes and high eccentricities. As such, we fix our integration timestep to be

$$\tau = 10^{-2} \tau_{\text{dyn}} \simeq 4.6 \text{ yr}; \tau_{\text{dyn}} = \frac{1}{\Omega_{\text{GR}}(a_{\text{min}}, e_{\text{max}})}. \quad (\text{F3})$$

Here $\Omega_{\text{GR}} = (3GM_*\Omega_{\text{Kep}}(a))/(c^2a(1-e^2))$ follows from Equation (B2). With this choice, for the (fast) wire with $(a, e) = (a_{\text{min}}, e_{\text{max}})$, the entire precession of the pericenter phase requires ~600 integration timesteps. Finally, we neglected the spin of the central BH, so that the sole relativistic precession was the 1PN Schwarzschild precession, as dictated by Appendix B.

It is informative to (crudely) compare the performance of the multipole runs with tentative direct N-body integrations. Fixing the direct integration timestep to $\tau_{\text{direct}} = 10^{-2}/\Omega_{\text{Kep}}(a_{\text{min}})$, multipole runs would use timesteps that are $\tau/\tau_{\text{direct}} \sim 200$ times larger than those of the direct integration. Integrating for a single timestep has a complexity scaling like $O(NKl_{\text{max}}^2)$ for the multipole method and $O(N^2)$ for direct integration. For the present values $N=10^4$, $l_{\text{max}} = 10$, and $K=100$, we find therefore that the multipole runs would typically be ~200 times less computationally intensive than direct N-body integrations.

As discussed in Section 7, we parallelize each simulation over four cores. We truncate the harmonic expansion at $l_{\text{max}} = 10$, and use $K=100$ nodes for each wire. We performed 500 realizations of the fiducial system, each of them integrated up to $t_{\text{max}} = 25$ kyr. With these choices of parameters, each realization took about 9 hr of computation for the MK4 scheme (Equation (45)) on four cores. Figure 11 presents measurements from the MK4 integrations.

To measure the diffusion coefficients presented in Figure 11, we proceed as follows. First, as $a$ is conserved for the secular dynamics, we split the wires according to the $a$-bins from Figure 11. For each $a$-bin, we split the wires into 20 logarithmic bins in $h$ for $10^{-1} \le h \le 1$, according to $h(t=0)$. To suppress the pollution stemming from the tails associated with large angular momentum changes, for each $a$-bin and $h$-bin, we compute $\Delta h = h(t_{\text{max}}) - h(0)$, and use the fractional moments (Bar-Or et al. 2013)

$$\{ (\Delta h)^2 \} = \lim_{\delta \to 0} \frac{\{ |\Delta h|^\delta \}^{2/\delta}}{2 \left( \Gamma \left[ \frac{1+\delta}{2} \right] / \sqrt{\pi} \right)^{2/\delta}}, \quad (\text{F4})$$

with {$\cdot$} the ensemble average over all realizations. In practice, we used $\delta = 10^{-3}$, and finally set $D_{hh} = \{\Delta h\}^2/t_{\text{max}}$. We emphasize that we did not perform any linear fits in time and only computed finite-time diffusion coefficients. Finally, to estimate the associated errors, we proceed by bootstrap resamplings, and in Figure 11 we represent the 16% and 84% error levels.

In the same Figure 11, we also represent the theoretical predictions for the finite-time diffusion coefficients computed following Bar-Or & Fouvry (2018) and the associated code scrrpy. In this prediction, we considered a maximum harmonic number given by $l_{\text{max}} = 10$. In order to match the exact setup considered here, we made two additional modifications to scrrpy compared to Bar-Or & Fouvry (2018). First, we restricted the range of $a$ and $e$ of the underlying stellar cluster to $a_{\text{min}} \le a \le a_{\text{max}}$ and $e_{\text{min}} \le e \le e_{\text{max}}$. Second, rather than having an exact resonance condition on the precession frequencies, we computed finite-time diffusion coefficients, following Equation (24) of Bar-Or & Fouvry (2018). This last modification is an important contributor to the agreement seen in Figure 11.

## ORCID iDs

Jean-Baptiste Fouvry https://orcid.org/0000-0002-0030-371X

Scott Tremaine https://orcid.org/0000-0002-0278-7180

Ben Bar-Or https://orcid.org/0000-0002-8927-4571

## References

Aarseth, S. 1967, in Les Nouvelles Méthodes de la Dynamique Stellaire, Vol. 47 (Paris: Editions du Centre National de la Recherche Scientifique)

Alexander, T. 2017, ARA&A, 55, 17

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

Bahcall, J. N., & Wolf, R. A. 1977, ApJ, 216, 883

Bar-Or, B., & Alexander, T. 2016, ApJ, 820, 129

Bar-Or, B., & Fouvry, J.-B. 2018, ApJL, 860, L23

Bar-Or, B., Kupi, G., & Alexander, T. 2013, ApJ, 764, 52

Binney, J., & Tremaine, S. 2008, Galactic Dynamics (2nd ed.; Princeton, NJ: Princeton Univ. Press)

Brouwer, D. 1937, AJ, 46, 149

Cohn, H., & Kulsrud, R. M. 1978, ApJ, 226, 1087

Dehnen, W. 2014, CompAC, 1, 1

Duncan, M. J., Levison, H. F., & Lee, M. H. 1998, AJ, 116, 2067

Fouvry, J.-B., Bar-Or, B., & Chavanis, P.-H. 2019, ApJ, 883, 161

Fry, J. N., & Peebles, P. J. E. 1980, ApJ, 236, 343

Gillessen, S., Plewa, P. M., Eisenhauer, F., et al. 2017, ApJ, 837, 30

Gürkan, M. A., & Hopman, C. 2007, MNRAS, 379, 1083

Hairer, E., Lubich, C., & Wanner, G. 2006, Geometric numerical integration (2nd ed.; Berlin: Springer)

Hamers, A. S., & Portegies Zwart, S. F. 2016, MNRAS, 459, 2827

Hamers, A. S., Portegies Zwart, S. F., & Merritt, D. 2014, MNRAS, 443, 355

Heckman, T. M., & Best, P. N. 2014, ARA&A, 52, 589

Hénon, M. 1964, Annales d'Astrophysique, 27, 83

Henon, M. 1973, A&A, 24, 229

Hernquist, L., & Ostriker, J. P. 1992, ApJ, 386, 375

Hopman, C., & Alexander, T. 2006, J. Phys.: Conf. Ser., 54, 321

Jacobs, V., & Sellwood, J. A. 2001, ApJL, 555, L25

Kazandjian, M. V., & Touma, J. R. 2013, MNRAS, 430, 2732

Klein, O. 1924, ZPhy, 22, 109

Szölgyén, Á., & Kocsis, B. 2018, PhRvL, 121, 101101

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

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

Ladner, R. E., & Fischer, M. J. 1980, J. ACM, 27, 831

Lightman, A. P., & Shapiro, S. L. 1977, ApJ, 211, 244

Madigan, A.-M., Hopman, C., & Levin, Y. 2011, ApJ, 738, 99

McGlynn, T. A. 1984, ApJ, 281, 13

McLachlan, R. I., Modin, K., & Verdier, O. 2014, PhRvE, 89, 061301

Meiron, Y., Li, B., Holley-Bockelmann, K., & Spurzem, R. 2014, ApJ, 792, 98
Merritt, D. 2013, Dynamics and Evolution of Galactic Nuclei (Princeton, NJ: Princeton Univ. Press)

Merritt, D., Alexander, T., Mikkola, S., & Will, C. M. 2011, PhRvD, 84, 044024

Milankovitch, M. 1939, Bull. Serb. Acad. Sci. Math. Nat. A, 6, 1

Munthe-Kaas, H. 1999, ApNM, 29, 115

Panamarev, T., Just, A., Spurzem, R., et al. 2019, MNRAS, 484, 3279

Rauch, K. P., & Tremaine, S. 1996, NewA, 1, 149

Rein, H., & Liu, S. F. 2012, A&A, 537, A128

Rein, H., & Spiegel, D. S. 2015, MNRAS, 446, 1424

Rosengren, A. J., & Scheeres, D. J. 2014, CeMDA, 118, 197

Saha, P. 1993, MNRAS, 262, 1062

Springel, V., & Hernquist, L. 2002, MNRAS, 333, 649

Sridhar, S., & Touma, J. R. 2016, MNRAS, 458, 4143

Touma, J. R., Tremaine, S., & Kazandjian, M. V. 2009, MNRAS, 394, 1085

Trefethen, L. N., & Weideman, J. A. C. 2014, SIAM Rev., 56, 385

Tremaine, S., Touma, J., & Namouni, F. 2009, AJ, 137, 3706

van Albada, T. S., & van Gorkom, J. H. 1977, A&A, 54, 121

Villumsen, J. V. 1982, MNRAS, 199, 493

White, S. D. M. 1983, ApJ, 274, 53