Secular Dynamics around a Supermassive Black Hole via Multipole Expansion

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 Gauss, is to use the secular approximation of smearing out the $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 re-formulate the discretisation 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 $\ell_{\mathrm{max}}$, reducing the overall computational costs from $O(N^2K^2)$ to $O(NK\ell_{\mathrm{max}}^2)$. Finally, we show that the averaged dynamical system is equivalent to $2N$ 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.


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 get slowly distorted by the additional perturbations present in the system.
The steep potential well generated by the central supermassive 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, Bar-Or, & Chavanis 2019), which, owing to non-spherical 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; (v) the non-resonant relaxation time (Bahcall & Wolf 1976;Lightman & Shapiro 1977;Cohn & Kulsrud 1978;Bar-Or et al. 2013), during which localised two-body encounters between stars lead to the longterm relaxation of the stars' Keplerian energy, i.e., the wires' semi-major axes.
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 non-resonant relaxation in the nucleus of the Milky Way, direct numerical simulations of these dynamical systems re-main very challenging. These 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 6 stars has been presented (Panamarev et al. 2019), and even in this case most of the stars lie outside the central BH's sphere of influence. Conversely, simulations in the very relativistic regime are still limited to a small number of particles N∼10 2−3 , should they use direct N-body methods (Merritt et al. 2011) or effective ones (see, e.g., Madigan, Hopman, & Levin 2011;Hamers, Portegies Zwart, & Merritt 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 wires' orbital parameters. This is in particular at the heart of the Gauß method (Touma, Tremaine, & Kazandjian 2009) which provides an efficient algorithm to compute the force between two such wires. Should one be 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 wires' dynamics over their in-plane precession, replacing Keplerian wires with Keplerian annuli, offering new venues to perform numerical investigations of these long-term dynamics (Kocsis & Tremaine 2015).
The main benefit from 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, 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 2 max ), with N the total number of stars, K a parameter independent of N, and max the maximum harmonics considered in the multipole expansion. We also show that the system is equivalent to 2N classical spin vectors, and devise time integration schemes that exactly comply with this system's geometric constraints.
The present paper is organised 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 ap-proach of discretising averages over the Keplerian motions as sums over nodes, and shows how to utilize 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 comparing them with predictions from kinetic theory. In Section 7 we discuss limitations of the present approach and possible future improvements, and we conclude in Section 8.

THE ORBIT-AVERAGED HAMILTONIAN
We consider N stars with masses m i and positions R i , orbiting a supermassive BH of mass M • at location R • . The total Hamiltonian of this system is with the canonical momenta P • = M •Ṙ• and P i = m iṘi . In this equation, the potential contribution Φ i gr accounts for the (conservative) relativistic corrections induced by the central BH 1 , i.e., the Schwarzschild and Lense-Thirring precessions (see, e.g, Merritt 2013), whose detailed expressions are given in Appendix B.

Democratic coordinates
In order to emphasise the dominant influence of the central BH on the system's dynamics, we rewrite equation (1) using democratic coordinates centered on the BH (Duncan et al. 1998) and their canonical momenta where we introduced the total mass M tot = M • + M with M = N i=1 m i the total stellar mass. Thus, p • is just the total momentum, which we set to zero without loss of generality, so that p i = P i is the barycentric momentum of the ith star. Following this change of coordinates, the Hamiltonian (1) becomes (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 term become irrelevant constants.

Orbital Elements
In order to describe the Keplerian dynamics imposed by the central BH, we transform the (r, p) coordinates to the Delaunay variables (M, ω, Ω, Λ, L, L z ) (Binney & Tremaine 2008) for each star. In this notation, the dynamical angles M, ω, and Ω are, respectively, the mean anomaly, the argument of pericentre, and the longitude of the ascending node. The associated actions are where Λ is the circular angular momentum of an orbit with the same energy or same semi-major axis a, L the magnitude of the angular momentum vector, and L z its projection onto the z-axis, while e and I denote, respectively, the eccentricity and inclination of the orbit. Kepler's equation 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). Following this change of variables, equation (3) becomes with 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, Λ 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 If the BH is supermassive, i.e., if one has M • M , 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 Here and in the remainder of this paper, we use the notation · for an average over all unperturbed stellar orbits.
Since H is independent of the M i , the associated actions Λ 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 p 2 i and products of the orbit-averaged stellar momenta p i . In our barycentric frame, p i = 0 and p 2 i = m 2 i GM • /Λ i 2 by virtue of the virial theorem. The constant terms depending only on the Λ i do not induce any dynamics and can be omitted, so the orbit-averaged Hamiltonian (9) finally becomes

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 characterised by the five quantities ω, Ω, Λ, L, L z , of which Λ 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 r(ω, Ω, L, L z ), at fixed (M, Λ). Such 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 will now pursue.
Rather than integrating the motion w.r.t. the orbital elements, we keep track of the wire dynamics through the dimensionless vectors where a hat denotes a unit vector as usual. Here, h is the angular momentum scaled to the circular angular momentum at the same energy and e is the eccentricity vector, which points in the direction of pericentre and has magnitude |e| = e. While only the four orbital elements ω, Ω, L, L z evolve, the two vectors h, e have six dynamical variables in total. The two associated degeneracies are captured by the two identities h · e = 0, h 2 + e 2 = 1.
In terms of these vectors, position and momentum are 3.1. Orbit-averaged rates of change Using equations (11) one can compute the time derivatives of the vectors h and e from the derivatives of the canonical variables r and p. The combined force from all other wires at position r i on wire i is where · \i denotes an average over all orbits except i and F(r) ≡ −∂ϕ/∂r with ϕ(r) ≡ −|r| −1 . At each value of the mean anomaly M i along the orbit, the force F i induces the local changes 2 In order to obtain the rates of changeḣ i andė 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 obtaiṅ Here, we used the shorthand F ij ≡ F(r i − r j ) and the averages are over all orbits including i.

Keeping all conservation laws
The expressions (16) satisfy the constraints (12) for any conservative field F ij ). 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 numerical treatment. The resulting numerical average will inevitably carry a small error with the consequence that the 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, since the conditions (12) are essential for the interpretation of h i and e i as a description of the wires. It is therefore important to construct numerical expressions forḣ i andė i that, despite their discretisation errors, keep the constraints (12) valid to machine precision.
In order to achieve that goal, we note that equations (16) are equivalent to Milankovitch's (1939) relations Rosengren & Scheeres 2014) Here, the averaged Hamiltonian H is expressed as a function of the vectors h i , e i (and the constants Λ i ) from all wires. These relations are equivalent to the canonical equations of motion, but much more useful. First, the vectors h and e are well-defined for all orbits, even those with zero eccentricity, unit eccentricity, and zero inclination. Second, they have the beautiful property that the identities (12) and the conservation of total energy are explicitly satisfied for any form of H , i.e., These properties suggest an alternative, fully conservative approach to numerically computing the ratesḣ i andė i : instead of following established practice of discretising the orbit averages in equations (16), we first discretise the orbit average for H and subsequently use Milankovitch's relations (17) to obtainḣ i andė i 3 , as detailed in Section 4.1.

Reformulation as a spin system
Following Klein (1924), we introduce the vectors then the identities (12) become Thus, the system of N wires is fully described by 2N independent unit vectors b = b i± , in other words it is a classical spin system. The dynamics of these vectors is simplẏ b =ḣ ±ė.
Since the vectors b remain on the unit sphere, this equation can be expressed as a precessioṅ The precession vectors Ω are not uniquely determined through equation (21), since a component parallel to b does not affectḃ (one would need to knowb to construct such a component). The most conservative choice for the precession vectors-in the sense that |Ω| 2 is minimized-is therefore for which Ω · b = 0, so that b moves along a great circle for constant Ω.
Milankovitch's equations (17), when re-expressed in terms of the vectors b, reaḋ These have the standard form for the equations of motion of classical spin systems. Comparing with equation (22), one might identify 2Λ −1 ∂ H /∂b with Ω. However, the gradient of the Hamiltonian is not unique, since owing to the constraints (20), H is only determined up to additive terms of the form (b 2 − 1)H with an arbitrary functionH. Such terms alter the gradient, but have no effect onḃ, 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 coming 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 a N = 10 4 stellar cluster mimicking Sgr A*. We refer to Appendices D, E, and F for the detailed description of these setups.

CALCULATING THE RATES OF CHANGE
Our computed rates of changeḣ i andė 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ḣ andė from H gr , namely the Schwarzschild and Lense-Thirring relativistic precessions.

Discretised rates of change
In order to numerically calculate the orbit averages, we approximate them by discrete midpoint sums over K positions along each orbit 4 : for any function f (M). In other words, each eccentric wire of mass m i is replaced by K nodes of masses µ ik . Placing the nodes' mean anomalies M ik uniformly (Gürkan & 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 large eccentricity, and hence produces a slow convergence with K. Instead, we place the nodes uniformly in eccentric anomaly: with ∆E = 2π/K. With this choice, the nodes are placed symmetrically w.r.t. pericentre but there is no node exactly at that position. To complete the discretisation (25), we specify the node masses With these specifications, the discretised potential energy of node k on wire i becomes such that the total interaction energy between the wires is computed as This expression depends on the h i and e i through the node positions r ik but also through their masses µ 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 H resulting from this dependence through µ ik induce via Milankovitch's relations (17) the rates of change (see Ap- whereμ i ≡ m i /K is the mean node mass on the wire. These rates of change do not appear in the traditional equations (16) based on forces.
The rates of change induced by the dependence of H on h i and e i through r ik are (see Appendix A) Here, is the gravitational force acting on node k of wire i and generated by the nodes of all other wires. Upon close inspectioṅ which is identical to the discretisation of equation (16a) and implies conservation of total angular momentum L tot = i Λ i h i , as is expected from Noether's theorem (since our discretisation of the wire averages remains invariant under spatial rotations). Also note that the discretisation of equation (16b) clearly differs from equation (31b). The advantage of equation (31b) is that it satisfiesḣ i · h i +ė i · e i = 0 anḋ h i · e i +ė i · h i = 0 exactly, as required.
We note from equations (30) and (31) that the computation ofḣ i andė i due to interactions with N wires requires O(N 2 K 2 ) computations. As shown below, this can be reduced to O(NK 2 max ) when approximating F(r ik − r jl ) by its expansion in spherical harmonics up to order max . To this end, it is advantageous to split 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 f all ik and f self ik must compute the exact same wire self-gravity such that when computing the difference f ik = f all ik − f self ik these erroneous nonphysical contributions exactly cancel. This is important, since 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.

Multipole Expansion
In order to accelerate the computation of the O(N 2 K 2 ) pairwise node interactions, we utilise the expansion of the Newtonian interaction kernel ϕ(r) = −|r| −1 in spherical harmonics Here, we have defined real-valued upper and lower solid spherical harmonics, U m (r) and T m (r) respectively, which are described explicitly in Appendix C. Note that there exist essentially two flavours 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); (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 min /r +1 max , as in equation (35), here absorbed into the definitions of U m and T m . Our goal here is to follow this second avenue and tailor it to the case of Keplerian wires. The expansion (35) is isotropic, even when truncated at order 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 ensures, by design, the conservation of the system's total linear momentum.
Inserting the expansion (35) into equation (34a), we have with an analogous expression for f self ik . The important difference between this expression and equation (34a) is that the dependence on the positions r ik and r jl has been factorized. Hence, the gradient terms depending on r ik can be taken outside the sum over nodes after splitting it into an inner and outer part, giving where P m (r) = j,l: r jl <r Q m (r) = j,l: r jl >r 1 sort all nodes in radius r α . 2 set f all α = 0 for all nodes α 3 set P m = 0 4 for all nodes α in order of increasing radius r α : 5 set Q m = 0 6 for all nodes α in order of decreasing radius r α : Table 1. The algorithm to compute f all α for all NK nodes. The forces on node α due to all other nodes at smaller and large radii, respectively, are computed in steps 4.1 and 6.1 from the inner and outer multipoles P m and Q m , which in turn are accumulated in steps 4.2 and 6.2.
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 ( max , K), via equation (10), which can be computed efficiently using the same algorithm as for the wire forces. In the limit max → +∞, i.e. in the absence of any multipole expansion, equation (10) provides us with the full total energy, E direct tot , which we can also compute via a direct O(N 2 K 2 ) sum over all pairwise node-node couplings.

Computing the wire forces
We begin by sorting the nodes into ascending order of radius 5 . Since the nodes on each wire are already sorted when created, this costs only O(NK lnN) 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 α = {i, k}, which we call the node index. Once the nodes are sorted in radius, the multipoles P m and Q m can be calculated by increasing and decreasing recurrence, respectively. The algorithm for the computation of f all α for all nodes is summarized in Kozai-Lidov 10 1 10 2 10 3 10 4 K 10 -8 10 -6 10 -4 10 -2 Cluster Figure 1. Illustration of the relative errors inḣ andė as a function of the number of nodes K for fixed multipole truncation order max , for Kozai-Lidov oscillations (Appendix E) and an N = 10 4 stellar cluster (Appendix F). We used 500 independent realisations 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.
In Figure 1, we illustrate the dependence of the relative errors in bothḣ andė as a function of the number of nodes per wire K.
For the Kozai-Lidov system, for a fixed value of 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 that case, the double orbit-average integral from equation (9) naturally splits into the product of two 1D integrals, which the multipole algorithm from equation (38) catches. Importantly, the integrands of each of these 1D integrals are 2π-periodic w.r.t. the eccentric anomaly. In that case, the midpoint sampling from equation (26) ensures an exponential convergence w.r.t. 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 ≤ e ≤ 0.99, see Appendix F -radial overlaps are ubiquitous and the exponential convergence w.r.t. K does not hold anymore. More precisely, for a fixed value of max , the error decreases as K −2.5 forḣ and as K −1.5 forė. We believe that this difference in the rates of convergence stems from the fact that the discontinuity in Cluster Figure 2. Illustration of the relative error in the truncated total energy, E tot ( max , K), as a function of the number of nodes K for fixed multipole truncation order max , w.r.t. to calculations with K max = 100 for the top panel and K max = 20 480 for the bottom panel. In both panels, the colored regions correspond to the 16% and 84% levels over 500 independent realisations. equation (36) is only in the radial force, which does not enter the equation forḣ.
We reach similar conclusions in Figure 2, where we present the error in the truncated total energy, E tot ( max , K) as a function of K. For the radially non-overlapping Kozai-Lidov simulations, it converges exponentially w.r.t. 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ḣ andė as a function of max for a fixed number of nodes K. For the Kozai-Lidov simulations, we find that the error decreases exponentially as e −1.6 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 min /r max ) ≤ η , with η the largest distance ratio that occurs as the inner and outer stars run through their orbits. For the present setup, one readily finds η = 1/5, see Appendix E. As a consequence, when truncated at order max , the typical error is expected to scale like η max e −1.6 max , matching the numerical measurement from Figure 3.
In the same figure, we note that the convergence w.r.t. max is much slower for the N = 10 4 clusters. This is a conse-  quence of the very slow convergence of the Legendre expansion (35) in any regime where the wires overlap in radius. We also note thatė is more accurate thanḣ. This is likely due to the analytical contribution from the Schwarzschild precession (see equation B6), which only affectsė and somewhat reduces the relative errors. Let us finally emphasise 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 w.r.t. max of the relative error between the truncated total energy, E tot ( max , K), and the full total energy, E direct tot . For the Kozai-Lidov simulations, we recover an exponential convergence w.r.t. max , as e −1.6 max . For the N = 10 4 clusters, where numerous orbits are radially overlapping, we empirically find that E tot ( max , K) converges roughly like −1.5 max w.r.t. 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 non-coplanar orbits, which mostly compose the present N = 10 4 clusters, a given harmonic, , is found to give a con- Cluster Figure 4. Illustration of the relative error in the truncated total energy, E tot ( max , K), as a function of the multipole truncation order max for a fixed number of nodes (K = 100 for Kozai-Lidov and K = 20480 for the cluster), using the same convention as in Figure 2. The full total energy, E direct tot , was computed by direct summation of equation (7) using respectively K = 20480 for the Kozai-Lidov simulations and K = 1 000 for the clusters.
tribution to the total truncated energy scaling like −2.5 . Once summing over all the harmonics 0 ≤ ≤ max , this ultimately leads to a relative error in E tot ( max , K) scaling like −1.5 max . Finally, in Figure 5, we illustrate the joint dependence of the relative error in the truncated total energy, E tot ( max , K), w.r.t. the full total energy, E direct tot , as a function of max and K for the cluster simulations. For the present system, we find that errors saturate for K 2 max . Because the computational complexity of the discretised harmonic expansion scales like O(K 2 max ), for a fixed value of max , ensuring converged estimations of E tot ( max , K) requires O( 3 max ) operations.

TIME INTEGRATION
In our application, the state of the system is described by the N sets y = {b + , b − } i . Given a timestep τ, an integration scheme is a mapping y n → y n+1 , where subscripts denote different time slices. Standard Runge-Kutta integration of the three-dimensional vectors b generally leads to violations of the constraints |b| = 1. Instead, we present here integration schemes that by construction exactly comply with these constraints. For convenience we define the velocity field B(y) =ḣ ±ė, so the equations of motions readḃ = B(y).

Explicit Lie group methods
10 -8 10 -7 10 -6 10 -5 2 10 1 10 2 10 1 10 2 10 3 ℓmax K K = 2 ℓmax K ℓ max 2 = cst. Figure 5. Illustration of the median relative error in the truncated total energy, E tot ( max , K), w.r.t. the full total energy, E direct tot , as a function of max and K for the cluster simulations, using the same parameters as in Figure 4. Errors are empirically found to saturate for K 2 max (black line), while the computational complexity scales like O(K 2 max ) (gray lines).
For a fixed value of the precession vector Ω = b × B (see equation 23) and an initial condition b 0 , equation (22) can be integrated exactly for a duration t, via Rodrigues' rotation formula, to the new location with |b(t)| = |b 0 |. In order to exactly preserve the constraints |b| = 1, explicit integration schemes analogous to Runge-Kutta schemes can be devised via Lie group methods (see, e.g., §IV.8 in Hairer et al. 2006) by concatenating appropriately 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 which imposes that U(t) evolves according tȯ with (23). Here, dφ −1 U is the inverse of the differential of the rotation map (see Theorem 3 in Munthe-Kaas 1999). It generically reads with B k the Bernoulli numbers and ad k U the k-th power of the adjoint operator. In the present context, this operator simply becomes 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φ −1 U appropriately truncated (Theorem IV.8.5 in Hairer et al. 2006). For a secondorder scheme, one can use dφ −1 U • Ω = Ω, from which one constructs the two-stage explicit midpoint rule (coined MK2). Starting from an initial state b n , it proceeds via Naturally the rotations are performed over all elements of y simultaneously. When constructing higher order methods, one must account for the fact that the rotations (39) do not commute, i.e. better approximations of dφ −1 U have to be used. Here, we implement a four-stage scheme (coined MK4) based on the classical fourth-order Runge-Kutta scheme. It reads where dφ −1 U is truncated at second order in U so that equation (42) becomes The schemes (44) and (45) are (i) explicit, (ii) conserve |b| = 1 exactly 6 , (iii) require two or four computations of the derivatives, respectively, and (iv) are, respectively, second-and fourth-order accurate. Because they are direct translations of usual Runge-Kutta methods, it is straightforward to design Munthe-Kaas schemes of higher order.

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 6 We systematically perform the re-normalisation b ← b/|b| after every evaluation of equation (39) to prevent a drift of |b| from round-off errors.

0
1×10 Figure 6. Illustration in the Kozai-Lidov simulations of the evolution of the inner star's eccentricity, e in , as well as its inclination w.r.t. the outer star,ĥ in ·ĥ 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.
The scheme from equation (47) is (i) second-order accurate; (ii) implicit; (iii) exactly conserves the constraints from equation (20)  The MK2, MK4 and MD2 schemes conserve the constraints for the vectors h, e to machine precision, which as far as we know none of the previous studies (e.g., Touma et al. 2009;Tremaine et al. 2009;Hamers & Portegies Zwart 2016) have managed.

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 3-body problem. Both methods predict similar oscillations of the inner star 8 . On average, the direct integration used an integration timestep of the order τ 10 −4 , while for the multipole integration, benefiting from  Figure 7. Illustration of the relative error in the total energy E tot w.r.t. its initial value after a fixed time T = 1.28, as a function of the considered integration timestep, τ. The colored regions correspond to the 16% and 84% levels over 500 independent realisations of the "Pair" Hamiltonian. We compute the relative errors by comparison with the initial total energy. See Appendix D for details.
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 in E tot after a finite time, as a function of the integration timestep τ. We recover that MK2 and MD2 are second-order accurate with a finite-time error converging like O(τ 2 ), while MK4 is fourth-order accurate.
In Figure 8, for the same Hamiltonian, we investigate 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; (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 ( max , K) and L tot for the N = 10 4 clusters. This figure exhibits similar asymptotic trends as in 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 # forces /(τ/τ dyn ), with # forces the average number of force evaluations per integration timestep and τ dyn the shortest dynamical time of the clusters at hand (see equation F31). While the explicit schemes always perform a constant number of force evaluations per timestep, the symplectic scheme typically requires from 5 evaluations (smallest τ) up to 12 (largest τ). As one reduces τ, the error in E tot ( max , K) keeps shrinking, illustrating that our discretisation scheme indeed corresponds to a Hamiltonian system (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 tot ( max , K), while the converse holds regarding L tot .  Figure 8. Illustration of the relative error in the total energy E tot (top) and total angular momentum L tot (bottom) w.r.t. their initial values, as a function of the number of integration steps, t/τ, using the same convention as in Figure 7. Here, we used τ = 10 −2 as our integration timestep.

SCALAR RESONANT RELAXATION
Using the method described above, we may now use our simulations of a 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 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 = {[∆h(T )] 2 }/T , where { · } stands for an ensemble average, and we used T =25 kyr, 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 mul- 10 -8 10 -7 |ΔLtot|/|Ltot| Cluster MD2 MK2 MK4 Figure 9. Illustration of the relative error in the truncated total energy E tot ( max , K) (top) and total angular momentum L tot w.r.t. their initial values, as a function of the number of integration steps, t/τ, for the N = 10 4 clusters, using the same convention as in Figure 8. Here, we used max = 10, K = 100, and τ = 10 −2 τ dyn as our integration timestep, see Appendix F. tipole method. Astrophysical applications to more diverse physical systems will be the subject of future work.

DISCUSSION
First, we note that the computation of the rates of change (ḣ,ė) for all the nodes has a computational complexity in O(NK 2 max ). This is the main benefit of the present method which, at the cost of (truncated) harmonic expansions and discretised orbit averages, has a computational complexity scaling linearly in the total number of wires.
We also note that the computations of {ḣ,ė} all and {ḣ,ė} self are completely independent (and of similar computational difficulty), so that they can be easily be performed in parallel. Similarly, for each of these computations, the recurrences over P m and Q m , i.e. steps 4 and 6 of the algorithm given in Table 1, are also independent one from another (and of similar computational difficulty), so that they can also be performed in parallel. As a consequence, one can easily benefit from a parallelisation of the previous algorithm over four cores, to compute in parallel {ḣ,ė} all and {ḣ,ė} self due to nodes at smaller as well as larger radii.
Further parallelisation requires the parallel computation of the prefix sums from equations (38) and is much more complex, but not impossible. We plan to investigate the appli-   Figure 10. Error-cost relation for the cluster simulations. Simulations were performed with (i) the symplectic scheme (MD2; star symbols) or the explicit ones (MK2 and MK4; circle symbols); with (ii) different integration timesteps, τ/τ dyn = 2 k × 10 −3 (0 ≤ k ≤ 7; different colors); (iii) up to different finite times, T/τ dyn = 0.128 × 2 −k (0 ≤ k ≤ 5): the larger T , the thicker the line. The cost of a given simulation is estimated via # forces /(τ/τ dyn ), with # forces the average number of force evaluations per integration timestep, e.g., two for MK2. Median errors in E tot ( max , K) (top) and L tot (bottom) are estimated over 128 realisations, with max = 10 and K = 100. cation 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 self , produces no torques and hence only affects each wire's apsidal precession rates,ė ·k/e. In addition, these terms scale trivially with semi-major axis and hence depend non-trivially only on eccentricity. An alternative possibility for correcting for these self-gravity contributions is therefore to use the scaling in a and a onedimensional 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 costs. Since our present implementation shares these costs between 4 CPUs, such a scheme would reduce this to 2 CPUs but not speed-up the computation overall. However, in an implementation featuring the parallel scan (see above) this would reduce the wall-clock time at 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 wires' orbital parameters, in particular 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 → 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 θ = 1, where the maximum errors hardly depend on max . For the errors to decrease significantly and hence offer a better efficiency at a given accuracy, one needs θ < 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 ), 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 in the number of wires. More precisely, its computational complexity scales like O(N 2 K 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 parallelisation; (iii) because it relies on the computation of pairwise interactions, Gauß' method can also be further parallelised over pairs of stars, contrary to the present multipole method; (iv) Gauß' method is straightforward to apply in the case of softened interactions ). One of the main drawbacks of Gauß' method is that its computational difficulty gets prohibitively large as N gets larger, 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 tot , scaling like t. It would be worthwhile to improve the algorithm to always ensure an unbiased error growth in √ t, following Brouwer's (1937) law, while still complying (exactly) with the constraints from equation (20). Similarly, one could investigate further the design of higher order structure-preserving integration schemes (see, e.g., Hairer et al. 2006).

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 the system's orbital constraints. We used this method to recover the predicted statistical properties of resonant relaxation of stellar eccentricities. This algo-rithm has two main advantages, namely relying on a explicit orbit average over the fast Keplerian motion to reduce the range of dynamical times in the system, and offering, through the multipole expansion, a computational complexity scaling linearly with the total number of particles. The code will be made available on request.
The present work is only a first step towards the development of faster integration methods for the dynamics of galactic nuclei, here seen as an archetype 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 inplane precession frequency for very eccentric wires. Similarly, the system's dynamical range is also enlarged by the ranges in semi-major 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 parallelisation, 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 anonymous referee for a very detailed and insightful report. JBF 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. ST acknowledges the support of the Natural Sciences and Engineering Research Council of Canada (NSERC), funding reference number RGPIN-2020-03885. BB 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 DISCRETISED EQUATIONS OF MOTION
Here, we detail the derivation of equations (30) and (31) from Milankovitch's relations (17) and the discretised Hamiltonian H as given in equation (29). We first consider the terms induced by the dependencies of H on h i and e i through the node masses µ ik . From equation (27), withμ i ≡ m i /K, and therefore from equation (29 Milankovitch's equations (17) allow us then to obtain equations (30). Next, we derive the terms induced by the dependencies of H on h i and e i through the node positions r ik . For a fixed vector z, equation (13a) gives us the derivative of the node positions via 9 Differentiating equation (29) and inserting these relations, we get Here, the 1PN contribution corresponds to the Schwarzschild precession that drives an in-plane (prograde) precession of the wire's pericentre, without any change in the wire's orbital plane. It readṡ with Ω 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's spin and leads to an orbital precession of the wire. It readṡ with |S| ≤ 1 the normalised spin angular momentum of the central BH.

C. SPHERICAL HARMONICS
In this Appendix, we spell out explicitly our convention and our implementation for the calculation of the solid spherical harmonics, following definitions very similar to the ones of Dehnen (2014).
The (complex) surface spherical harmonics are defined as with P m the associated Legendre functions. The (complex) solid spherical harmonics are defined as In order to work only with real-valued quantities, we define the associated real-valued solid harmonics as and Let us note that this definition differs from the definition in Eqs. (58a) and (58b) of Dehnen (2014), where here we added a factor √ 2 to the expressions for m 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 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 will now detail. The important property of all the coming recurrences is that, for a given location r, they allow for the computation of U m (r), T m (r) and their derivatives with complexity O( 2 max ).

C.1. Computing U m
The upper real-valued solid harmonics satisfy the initialisation condition They satisfy some boundary recurrence relations for m = ± . For = 1, it reads with r = (x, y, z), while for ≥ 2, it reads Finally, inside the boundary of the ( , m)-domain, they satisfy the generic recurrence relation which can also be applied to the cases |m| = − 1, if one follows the convention U ±( −1) −2 = 0.

C.2. Computing T m
The lower real-valued solid harmonics satisfy the initialisation condition They satisfy some boundary recurrence relations for m = ± . For = 1, it reads while for ≥ 2, it reads Finally, inside the boundary of the ( , m)-domain, they satisfy the generic recurrence relation which can also be applied to the cases |m| = − 1, if one follows the convention T ±( −1) The gradients of the upper spherical harmonics, ∇U m , can be determined using our previous computation of U m . For m = 0, we have For m = 1, we use while for m ≥ 2, this relation becomes (C23) Similarly, the gradients of the lower spherical harmonics, ∇T m , are computed using the predetermined values of T m . For m = 0, we have For m = 1, we use while for m ≥ 2, this relation becomes (C26)

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 where a, b are constant coefficients and we work in dimensionless units Λ 1 = Λ 2 = 1. Following the Milankovitch equations (17), the equations of motion for the two wires reaḋ In Figures 7 and 8, we considered 500 independent realisations where (i) the coefficients a,b were drawn uniformly in [0,1]; (ii) the wires' initial orientationsĥ,ê were drawn uniformly on the unit sphere; (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 max ) or nodes sampling (via K).

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 out = 10, semi-major axis a out = 10 and initial eccentricity e out = 0.5, while for the inner star we take m in = 1, a in = 1 and e in = 0.01. At the initial time, we set the orientation of the outer orbit to (ĥ out ,ê out ) = (ê z ,ê y ). Finally, for the inner star, we imposeĥ in = cos I inêz − sin I inêy andê in = cos ω inêx + sin ω inĥin ×ê x , with the initial inclination I in = 60 • and pericentre phase ω in = 90 • . To perform the multipole integration, we use K = 100 nodes, max = 10 for the harmonic truncation, and τ = 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 realisations of the system where we drew the initial inclinations and pericentre phase of the inner star uniformly within ±5 • around their fiducial values, I in and ω in .
In order to check the validity of our multipole algorithm, we compare it with direct integrations of the associated 3-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 two stars' initial mean anomalies were picked at random. We considered a total of 100 different realisations, and found that e in andĥ in ·ĥ out only differed by 10 −4 at most over these realisations. As such, in Figure 6, we could safely limit ourselves to only representing the median values obtained over the available realisations. Finally, we used the mappings from equations (2b) and (11) to infer the stars' orbital elements from their individual locations and velocities.

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 × 10 6 M . We assume that the initial eccentricity distribution of the stars follows a thermal distribution, f e ∝ e, within the range e min ≤ e ≤ e max , with (e min , e max ) = (0, 0.99). For the initial distribution in the stars' semi-major axes, we assume that that it follows a (truncated) power-law distribution of the form f E (E) ∝ |E| p , with E = E Kep the Keplerian energy, and the power index p < 3 2 . With such a choice, the number of stars per unit a is given by with γ = p + 3 2 . For our fiducial model, we assume that the stellar distribution follows γ = 7 4 , the slope of the expected equilibrium (Bahcall & Wolf 1977). We also introduced N 0 = g(γ)N(< a 0 ) with g(γ)=2 −γ √ π Γ(1+γ)/Γ(γ− 1 2 ). Here, N(< a 0 ) is the number of stars physically within a radius a 0 , which should not be confused with N 0 , the number of stars with a semi-major axis smaller than a 0 . For the fiducial model, we assume that the total stellar mass 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 . It finally remains to carefully pick the value of m . To do so, we fix ourselves a range of interest in semi-major axis, namely a min ≤ a ≤ a max , with (a min , a max ) = (5, 100) in milliparsecs, to which the sampling of the particles' initial semi-major axes is limited. The fiducial simulations are performed with a total number of stars, N = 10 4 . Following equation (F29), this imposes a self-consistent relation between N and m , namely Such a choice ensures that despite the limited range in semi-major axis, the stellar self-consistent potential remains close to the full expected one. For the fiducial model with N = 10 4 , we find m 8.63 M . In all these simulations, the integration timestep gets dictated by the timescale associated with the (fast) in-plane relativistic precessions of the wires with small semi-major axes and large eccentricities. As such, we fix our integration timestep to be τ = 10 −2 τ dyn 4.6 yr with τ dyn = 1 Ω gr (a min , e max ) .
where Ω gr = (3GM • Ω Kep (a))/(c 2 a(1−e 2 )) follows from equation (B6). With such a choice, for the (fast) wire with (a, e) = (a min , e max ), the entire precession of the pericentre's 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 τ direct = 10 −2 /Ω Kep (a min ), multipole runs would use timesteps that are τ/τ direct ∼ 200 times larger than those of the direct integration. Integrating for a single timestep has a complexity scaling like O(NK 2 max ) for the multipole method and O(N 2 ) for direct integration. For the present values N = 10 4 , 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 parallelise each simulation over four cores. We truncate the harmonic expansion at max = 10, and use K = 100 nodes for each wire. We performed 500 realisations of the fiducial system, each of them integrated up to t max = 25kyr.
With these choices of parameters, each realisation 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 the a-bins from Figure 11. For each a-bin, we split the wires in 20 logarithmic bins in h for 10 −1 ≤ h ≤ 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 ∆h=h(t max ) − h(0), and use the fractional moments (Bar-Or et al. 2013) with { · } the ensemble average over realisations. In practice, we used δ=10 −3 , and finally set D hh ={∆h) 2 }/t max . We emphasise 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 that prediction, we considered a maximum harmonic number given by 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 min ≤ a ≤ a max and e min ≤ e ≤ e 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.