Issue 
A&A
Volume 598, February 2017



Article Number  A71  
Number of page(s)  24  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/201629138  
Published online  02 February 2017 
The secular evolution of discrete quasiKeplerian systems
I. Kinetic theory of stellar clusters near black holes
^{1} Institut d’Astrophysique de Paris and UPMC, CNRS (UMR 7095), 98bis boulevard Arago, 75014 Paris, France
email: fouvry@iap.fr
^{2} Korea Institute of Advanced Studies (KIAS), 85 Hoegiro, Dongdaemungu, 02455 Seoul, Republic of Korea
^{3} Rudolf Peierls Centre for Theoretical Physics, University of Oxford, Keble Road, Oxford OX1 3RH, UK
Received: 17 June 2016
Accepted: 6 October 2016
We derive the kinetic equation that describes the secular evolution of a large set of particles orbiting a dominant massive object, such as stars bound to a supermassive black hole or a protoplanetary debris disc encircling a star. Because the particles move in a quasiKeplerian potential, their orbits can be approximated by ellipses whose orientations remain fixed over many dynamical times. The kinetic equation is obtained by simply averaging the BBGKY equations over the fast angle that describes motion along these ellipses. This socalled BalescuLenard equation describes selfconsistently the longterm evolution of the distribution of quasiKeplerian orbits around the central object: it models the diffusion and drift of their actions, induced through their mutual resonant interaction. Hence, it is the master equation that describes the secular effects of resonant relaxation. We show how it captures the phenonema of mass segregation and of the relativistic Schwarzschild barrier recently discovered in Nbody simulations.
Key words: galaxies: kinematics and dynamics / galaxies: nuclei / diffusion / gravitation
© ESO, 2017
1. Introduction
The stars in a stellar cluster surrounding a dominant supermassive black hole (BH) move in a quasiKeplerian potential. Their orbits are ellipses that maintain their spatial orientation for many orbital periods. So, for many purposes, the cluster can be thought of as a system of massive wires, in which the mass of each star is smeared out along the path traced by its quasiKeplerian orbit. The consequences of this idea were first developed by Rauch & Tremaine (1996), who showed that wirewire interactions greatly enhance the relaxation of the stars’ angular momenta when compared to conventional estimates that ignore the coherence of the stars’ orbits over many dynamical times and consider only uncorrelated twobody encounters. They named this phenomenon “resonant relaxation”, because such enhanced relaxation occurs more generally in any potential in which the threedimensional vector of stellar orbital frequencies Ω satisfies a commensurability condition of the form n·Ω ≃ 0 for some vector of integers n = (n_{1},n_{2},n_{3}).
Understanding the effects of these relaxation processes in galactic nuclei is important when predicting the rates of tidal disruptions of stars by black holes (e.g. Rauch & Tremaine 1996; Rauch & Ingalls 1998), predicting merging rates of binary supermassive black holes (e.g. Yu 2002) or of gravitational wave signatures from starBH interactions (e.g. Hopman & Alexander 2006; Merritt et al. 2011). Resonant relaxation provides as well a promising framework for explaining some of the puzzling features of the young stellar populations found at the centre of our own Galaxy (e.g. Kocsis & Tremaine 2011).
A first way to study the evolution of such star clusters is by direct Nbody simulations. Unfortunately, extracting physical insights from such simulations is challenging, because the complex dynamical processes involved are entangled and, more practically, the computational costs of running the simulations typically mean that one can run just a few realisations, each with relatively small N. For problems that focus on resonant relaxation phenomena, one can often do better by using Nwires codes (e.g. Kocsis & Tremaine 2015) in which individual stars are replaced by orbitaveraged Keplerian wires.
A complementary way of understanding these systems is by using the tools of kinetic theory. For plasmas, Balescu (1960) and Lenard (1960) have developed a rigorous kinetic theory that takes the most important collective effects into account. In this theory, the coupled evolution equations for the system’s onebody distribution function and its twobody correlation function are reduced to a single equation – the BalescuLenard equation – that describes the evolution of the onebody distribution function alone. See Chavanis (2010, 2013a,b), Fouvry et al. (2016a), for a review on the early development of kinetic theory for plasmas, stellar systems, and other longrange interacting systems. The original BalescuLenard formalism was developed for homogeneous plasmas. One way of generalising it to inhomogeneous selfgravitating systems, such as star clusters, was proposed by Gilbert (1968). Sridhar & Touma (2016a,b) have recently applied Gilbert’s methods to the secular evolution of a star cluster around a BH, which is the problem addressed by the present paper.
An alternative way of generalising the BalescuLenard formalism to inhomogeneous systems has been presented by Heyvaerts (2010) and Chavanis (2012), who reformulate the nonlinear kinetic equation in terms of the angleaction variables that are appropriate for spatially inhomogeneous multiperiodic systems. Fouvry et al. (2015a,b, 2016b) have applied this formalism to describe the secular response of tepid selfgravitating stellar discs. The resulting inhomogeneous BalescuLenard equation accounts for the selfinduced secular orbital diffusion of a selfgravitating system driven by the internal shot noise due to the finite number N of particles involved. In common with all results based on the BalescuLenard formalism, it is valid to order O(1 /N) in a formal expansion of the dynamics ordered by the small parameter 1 /N. Therefore, it describes the evolution of the system on timescales of about Nt_{d}, where t_{d} is the dynamical time. The secular interactions between particles need not be local in space: they need only correspond to gravitationally amplified longrange correlations via resonances.
In its original form, however, the BalescuLenard formalism assumes that resonances are localised in action space and not degenerate. Therefore, it must be reexamined before it can be applied to the degeneracies inherent to resonant systems. In this paper, we show how to account for these degeneracies in the case of a cluster of N particles orbiting a massive, possibly relativistic, central body. We first average the equations of motion over the fast angle associated with the orbital motion of the stars around the BH. Once such an averaging is carried out, it turns out that the general formalism of the inhomogenous BalescuLenard equation applies straightforwardly and yields the associated secular collisional equation. This equation captures the diffusion and drift of particles’ actions induced through their mutual resonant interaction at the frequency shifts present in addition to the mean Keplerian dynamics, for instance induced by the selfgravity of the cluster or relativistic effects. Hence it is well suited to describe the secular evolution of a large set of particles orbiting around a massive object, for instance to account for the longterm evolution of a disc or a sphere of (possibly relativistic) stars near a galactic centre, or a protoplanetary debris disc circling a star. As such, it captures the secular effects of a sequence of polarised wirewire interactions (corresponding to scalar or vector resonant relaxation) on the underlying orbital structure of the cluster.
The paper is organised as follows. Section 2 derives the BBGKY hierarchy of a system with a massive central body using canonical coordinates to account properly for the black hole’s motion. Section 3 introduces angleaction coordinates for such quasiKeplerian systems. Section 4 averages the corresponding dynamical equations over the fast angles of the Keplerian motion. Section 5 presents the degenerate one and multicomponent Keplerian BalescuLenard equations (Appendix B details both derivations, following the steps of Heyvaerts 2010, while also correcting for a minor issue in the multicomponent case). Section 6 discusses applications to the cases of razorthin axisymmetric and spherical clusters orbiting a massive central object, and compares our results to those of Sridhar & Touma (2017) and others, while Sect. 7 concludes. Appendix A outlines the relativistic precessions frequencies involved near a massive black hole, and Appendix C presents the stochastic counterpart of the Keplerian BalescuLenard equation.
2. The BBGKY hierarchy
Consider a system of N stars in motion about a central black hole of mass M_{•}, in which each star has mass μ. We assume that the total stellar mass M_{⋆} ≡ μN is small enough that the ratio (1)Let X_{•} be the location of the BH and X_{n} be the location of the nth star referred to an inertial frame. The Hamiltonian for the system is then given by (2)in which the canonical momenta are given by and . Here, U(  X ) corresponds to the interaction potential, that is U(  X  ) ≡ −G/  X  in the gravitational context. In Eq. (2), the first two terms correspond to the kinetic energy of the BH and the stars. The third term corresponds to the Keplerian potential of the BH, while the fourth term is associated with the pairwise interactions among stars. Finally, the third line of Eq. (2)accounts for the relativistic correction forces such as the Schwarzschild and LenseThirring precessions occurring in the vicinity of the BH (see Appendix A), where the normalisation prefactor μM_{⋆} was added for later convenience. For simplicity, we neglected any additional external perturbations, which could offset the system. This will be the subject of a future work.
Let us now rewrite the Hamiltonian from Eq. (2)as N decoupled Keplerian Hamiltonians plus perturbations. We follow Duncan et al. (1998) and carry out a canonical transformation to a new set of coordinates, the democratic heliocentric coordinates (x_{•},x_{1},...,x_{N}) defined as (3)where we have introduced the total mass of the system M_{tot} = M_{•} + M_{⋆}. In Eq. (3), x_{•} corresponds to the position of the system’s centre of mass and x_{i} to the locations of the stars in the frame centred on the BH. These relations have inversion (4)As obtained in Duncan et al. (1998), the associated canonical momenta (p_{•},p_{1},...,p_{N}) are (5)Within these new canonical coordinates, the Hamiltonian from Eq. (2)takes the form (6)which consists of N independent Keplerian Hamiltonians (first term of the first line) plus the twobody couplings among them (second term) plus additional kinetic terms (second line). The evolution of the total momentum p_{•} is given by . Without loss of generality, we may therefore assume that p_{•} = 0. The evolution of the barycentre position is then given by , and we therefore set x_{•} = 0. Introducing the notation , the Hamiltonian from Eq. (6)becomes (7)in which one of the kinetic terms in the second line has been transformed away.
In order to obtain a statistical description of the system, we now introduce its Nbody probability distribution function (PDF) P_{N}(Γ_{1},...,Γ_{N},t) defined so that P_{N}(Γ_{1},...,Γ_{N},t) dΓ_{1}...dΓ_{N} is at time t the probability of finding particle 1 within the volume element dΓ_{1} located at the phase space point Γ_{1} = (x_{1},v_{1}), particle 2 within dΓ_{2} of the phase space point Γ_{2} = (x_{2},v_{2}), and so on. We normalise P_{N} such that (8)It evolves according to Liouville’s equation (9)The dynamics of the individual particles are given by Hamilton’s equations, μdx_{i}/ dt = ∂H/∂v_{i} and μdv_{i}/ dt = −∂H/∂x_{i}, where the system’s Hamiltonian was obtained in Eq. (7). From P_{N}, we define reduced PDFs, (10)by integrating over the phase space locations of particles n + 1 to N. To obtain the evolution equation of any reduced PDF P_{n}, we integrate Liouville’s Eq. (9)over dΓ_{n + 1}...dΓ_{N} and use the fact that P_{N} and H are unchanged under permutations of their arguments. This leads to the general term of the BBGKY hierarchy (11)Here, we have written the force exerted by particle j on particle i as μℱ_{ij} = −μ∂U_{ij}/∂x_{i}, using the shorthand notation U_{ij} = U(x_{i}−x_{j}). The force exerted by the BH on particle i is denoted by M_{•}ℱ_{i0} = −M_{•}∂U_{i0}/∂x_{i} and the force associated with the relativistic corrections as M_{⋆}ℱ_{ir} = −M_{⋆}∂Φ_{rel}/∂x_{i}.
It is convenient to replace these PDFs by the reduced distribution functions (DFs) (12)in terms of which Eq. (11)can be rewritten as (13)To isolate the contributions to f_{n} that arise from correlations among particles, let us introduce the cluster representation of the DFs. We define the 2body correlation function g_{2} in terms of f_{1} and f_{2} via (14)Similarly, the 3body correlation function g_{3} is defined by (15)These correlation functions have simple dependence on the number of particles N. It is straightforward to check that the following normalisations hold: (16)As the individual mass scales like μ ~ 1 /N, one immediately has  f_{1}  ~ 1,  g_{2}  ~ 1 /N, and  g_{3}  ~ 1 /N^{2}. Using the decompositions from Eqs. (14)and (15), after some simple algebra, the first two equations of the BBGKY hierarchy from Eq. (13)become (17)and (18)where (1 ↔ 2) means that all preceding terms are written out again, but with indices 1 and 2 swapped.
We now use the scalings obtained in Eq. (16)to truncate Eqs. (17)and (18)at order 1 /N. Notice that the system includes two small parameters, namely 1 /N associated with the discreteness of the system and ε = M_{⋆}/M_{•} associated with the amplitude of the nonKeplerian components. As will be emphasised in the upcoming calculations, we will perform kinetic developments, where we only keep terms of the order ε and ε/N. In Eq. (17), all the terms are of order 1 /N or larger, and should therefore all be kept. In Eq. (18), the first four lines are of order 1 /N (except for the correction (ε/N)(v_{1} + v_{2})·∂g_{2}/∂x_{1} which may be neglected), while all the terms from the two last lines are of order 1 /N^{2} and may therefore be neglected. Notice that the first term on the fifth line of Eq. (18), which, while being of order 1 /N^{2}, can nevertheless get arbitrarily large as particles 1 and 2 get closer. This term accounts for strong collisions between two particles, which are not accounted for in the present formalism. In addition to these truncations, and in order to consider terms of order 1, let us finally introduce the system’s 1body DF F and its 2body autocorrelation function as (19)Moreover, in order to emphasise the various order of magnitude of the forces present in the problem, let us also rescale some of the quantities appearing in Eqs. (17)and (18). Let us first rescale the interaction potential using the mass of the central black hole, so as to have the relations (20)Similarly, the relativistic potential Φ_{r} = Φ_{rel} is also rescaled so that (21)Following these various truncations and renormalisations, Eq. (17)becomes (22)while Eq. (18)becomes (23)The next step of the calculation involves rewriting Eqs. (22)and (23)within appropriate angleaction coordinates allowing us to capture in a simple manner the dominant mean Keplerian motion due to the central BH. When considering Keplerian potentials, one has to deal with additional dynamical degeneracies between the orbital frequencies, which should be handled with care, as we will now detail.
3. Degenerate angleaction coordinates
In Eqs. (22)and (23), one can note the presence of an advection term v_{1}·∂/∂x_{1} + ℱ_{10}·∂/∂v_{1} associated with the Keplerian motion driven by the central black hole. The next step of the derivation is to introduce the appropriate angleaction coordinates (Goldstein 1950; Born 1960; Binney & Tremaine 2008) to simplify this integrable Keplerian motion. We therefore remap the physical coordinates (x,v) to the Keplerian angleaction ones (θ,J). Along the unperturbed Keplerian orbits, the actions J are conserved, while the angles θ are 2πperiodic, evolving with the frequency Ω_{Kep} defined as (24)where H_{Kep} is the Hamiltonian associated with the Keplerian motion due to the black hole. For 3D spherical potentials, the usual angles and actions (Binney & Tremaine 2008) are given by (25)where J_{r} is the radial action, L the magnitude of the angular momentum, and L_{z} its projection along the zaxis. The Keplerian Hamiltonian then becomes H_{Kep} = H_{Kep}(J_{r} + L). Another choice of angleaction coordinates in 3D is given by the Delaunay variables (Sridhar & Touma 1999; Binney & Tremaine 2008) defined as (26)In Eq. (26), (I = J_{r} + L,L,L_{z}) are the three actions of the system, while (w,g,h) are the associated angles. Here, w stands for the orbital phase or mean anomaly, g for the angle from the ascending node to the periapse, and h for the longitude of the ascending node. With these variables, one has H_{Kep} = H_{Kep}(I), so that the angle w advances at the frequency ẇ = Ω_{Kep} = ∂H_{Kep}/∂I, while the angles g and h are constant. The existence of these additional conserved quantities makes the Keplerian potential dynamically degenerate. This can have some crucial consequences on its longterm behaviour, as we will now detail.
To clarify the upcoming discussions we denote as d the dimension of the considered physical space, for instance d = 2 for a razorthin disc. In this space, we consider an integrable potential ψ and an associated angleaction mapping (x,v) → (θ,J). A potential is said to be degenerate if there exists n ∈ Z^{d} such that (27)where it is understood that the vector n is independent of J, so that the degeneracy is global. A given potential may have more than one such degeneracy, and we denote as k the degree of degeneracy of a potential, that is the number of linearly independent vectors n satisfying Eq. (27). For example, for the angleaction coordinates from Eq. (25), the frequencies and degeneracy vectors are given by (28)so that k = 2. Using the Delaunay angleaction coordinates from Eq. (26), one can similarly write (29)which also gives k = 2. The degree of degeneracy of the potential is independent of the chosen angleaction coordinates. The Delaunay variables from Eq. (26)appear as a simpler choice than the usual ones from Eq. (25), because of their simpler degeneracy vectors.
For a given degenerate potential, one can always remap the angleaction coordinates to get simpler degeneracies. Indeed, let us assume that in our initial angleaction coordinates (θ,J), we have at our disposal k degeneracy vectors n_{1}, ... , n_{k}. Thanks to a linear transformation, we may change coordinates (θ,J) → (θ′,J′), so that in the new coordinates the k new degeneracy vectors take the simple form , where e_{i} are the natural basis elements of Z^{d}. Following Morbidelli (2002), as the vectors n_{i} are assumed to be linearly independent, we may complete this family with d−k vectors n_{k + 1},...,n_{d} ∈ Z^{d} to have a basis over Q^{d}. We then define the transformation matrix of determinant 1 as (30)and the new angleaction coordinates (θ′,J′) are defined as (31)One can check that (θ′,J′) are indeed new angleaction coordinates, with J′ conserved and θ′ ∈ [0,2π]. Within these new coordinates, the k degeneracy vectors are immediately given by , that is the intrinsic frequencies satisfy for 1 ≤ i ≤ k. The degeneracies of the potential got simpler. In the upcoming calculations, we will always consider such simpler angleaction coordinates, and we introduce the notations (32)where θ^{s} and J^{s} respectively stand for the slow angles and actions, while θ^{f} and J^{f} stand for the fast angles and actions. Finally, we introduced ℰ as the vector of all the conserved quantities (for a Keplerian potential, this corresponds to a Keplerian elliptical wire). For a degenerate potential, the slow angles are the angles for which the associated frequencies are equal to 0, while these frequencies are nonzero for the fast angles. Let us finally define the degenerate angleaverage with respect to the fast angles as (33)We now use these various properties to rewrite Eqs. (22)and (23)using the angleaction coordinates appropriate for the Keplerian motion due to the central BH. In these coordinates, the Keplerian advection term becomes (34)A nice property of the average from Eq. (33), is that the collisionless advection term from Eq. (34)then naturally vanishes, so that one has (35)Finally, the mapping (x,v) → (θ,J) preserves the infinitesimal volumes so that dΓ = dxdv = dθdJ. In addition, it also preserves Poisson brackets, so that for two functions G_{1}(x,v), and G_{2}(x,v), one has (36)In order to shorten the notations, let us now introduce the rescaled selfconsistent potential Φ as (37)One can now rewrite Eq. (22)within these angleaction coordinates and it takes the form (38)where we have written and have introduced the notation (39)so that it corresponds to the Poisson bracket with respect to the variables 1. In Eq. (38), the terms of the second and third lines are associated with the additional kinetic terms appearing in the Hamiltonian from Eq. (7). As we will emphasise later on, once averaged over the fast Keplerian dynamics, these terms will be negligible at the order considered here. Similarly, one can straightforwardly rewrite Eq. (23)as (40)where the terms from the two last lines are associated with the additional kinetic terms from Eq. (7), and will become negligible once averaged over the fast Keplerian dynamics. The rewriting from Eq. (38)is particularly enlightening, since one can easily identify in its first line the three relevant timescales of the problem. These are: i) the dynamical timescale T_{Kep} = 1/Ω_{Kep} associated with the Keplerian advection term ; (ii) the secular collisionless timescale of evolution T_{sec} = ε^{1}T_{Kep} associated with the potential contributions ε [ Φ + Φ_{r}]; and finally (iii) the collisional timescale of relaxation T_{relax} = NT_{sec}, associated with the last term in the first line of Eq. (38).
4. Fast averaging the evolution equations
Starting from Eqs. (38)and (40), let us carry out an average over the degenerate angles as defined in Eq. (33). We recall that the main virtue of such an averaging is to naturally cancel out any contributions associated with the Keplerian advection term, as observed in Eq. (35). We start from Eq. (38)and multiply it by . In order to estimate the average of the various crossed terms in Eq. (38), let us assume that the DF of the system can be expanded as (41)where ϵ ≪ 1 is a small parameter of order 1 /N. This ansatz is the crucial assumption of the present derivation. Indeed, the BH’s domination on the dynamics strongly limits the efficiency of violent relaxation or phase mixing to allow for a rapid dissolution of any dependence on θ^{f}. Hence it is somewhat arbitrarily assumed here that this condition has been achieved, so that, for our purposes, the system starts in a phasedmixed state.
We now discuss in turn how the various terms appearing in Eq. (38)can be averaged with respect to the fast Keplerian angle. In the first Poisson bracket of Eq. (38), one should keep in mind that the selfconsistent potential Φ, introduced in Eq. (37), should be seen as a functional of F. As a consequence, this term takes the form (42)where the averaged selfconsistent potential was introduced as (43)In Eq. (43), for clarity, the notation was shortened for the selfconsistent potential as . The (doubly) averaged interaction potential is defined as (44)while the angleaveraged potential was also introduced as (45)where the prefactor 1/(2π)^{d−k}, was introduced for convenience. As emphasised in Eq. (42), one should note that at first order in ε and zeroth order in ϵ, the selfconsistent potential has to be computed while only considering the averaged system’s DF .
To deal with the second Poisson bracket of Eq. (38), the same double average as introduced in Eq. (44)should be performed on . As we did for Eq. (41), it is assumed that the 2body correlation can be developed as (46)At first order in ε and zeroth order in ϵ, the third term from Eq. (38)can immediately be rewritten as Finally, at first order in ε and zeroth order in ϵ, the terms from the two last lines of Eq. (38)will involve the quantities (47)The first identity comes from the fact that Keplerian orbits are closed, so that the mean displacement over one orbit is zero, while the second identity comes from the virial theorem. As these terms either vanish or do not depend on the slow coordinates θ^{s} and J^{s}, they will not contribute to the dynamics at the orders considered here once averaged over the fast angle. Therefore, keeping only terms of order ε and ε/N, one can finally rewrite Eq. (38)as (48)In Eq. (48), we note that all the functions appearing in the Poisson brackets only depend on . As a consequence, the Poisson brackets defined in Eq. (36)take the shortened form (49)so that only derivatives with respect to the slow coordinates appear. Let us finally introduce the rescaled time τ as (50)so that Eq. (48)becomes (51)One may use a similar angleaveraging procedure for the second equation of the BBGKY hierarchy. Indeed, multiplying Eq. (40)by , relying on the developments from Eqs. (41)and (46), and keeping only terms of order ε, Eq. (40)can finally be rewriten as (52)where one can note that all the additional kinetic terms of the two last lines of Eq. (40)vanish at the considered order, when averaged over the fast Keplerian angle.
Equations (51)and (52)are the main results of this section. They describe the coupled evolutions of the system’s averaged DF, and 2body correlation . A rewriting of the same pair of equations has recently been derived by Sridhar & Touma (2016a,b) using Gilbert’s method. At this stage, one could investigate at least four different dynamical regimes of evolution for the system:

I.
Considering Eq. (51), the Keplerian wires couldinitially be far from a quasistationary equilibrium, so that. One then expects that this outofequilibrium system will undergo a phase of violent relaxation (LyndenBell 1967), allowing it to rapidly reach a quasistationary equilibrium. We do not investigate this process here, but still rely on the assumption that the collisionless violent relaxation of the wires’ DF can be sufficiently efficient for the system to briefly reach a quasistationary stable state, which will then be followed by a much slower secular evolution, either collisionless or collisional.

II.
For a given DF of stationary wires, one could also investigate the possible existence of collisionless dynamical instabilities associated with the collisionless part of the evolution Eq. (51), namely . Such instabilities are not considered in the present paper, and we will assume, as will be emphasised in the upcoming derivations, that throughout its evolution the system always remains dynamically stable with respect to the collisionless dynamics. See for instance Tremaine (2005), Polyachenko et al. (2007), Jalali & Tremaine (2012) for examples of stability investigations in this context.

III.
Once it is assumed that the system has reached a quasistationary stable state, one can study the secular evolution of this system along quasistationary equilibria. Such a longterm evolution can first be induced by the presence of external stochastic perturbations. To capture such a secular collisionless evolution, one should neglect contributions from the collisional term in 1 /N in Eq. (51), and look for the longterm effects of stochastic perturbations. The formalism appropriate for such a secular collisionless stochastic forcing is similar to the one presented in Fouvry et al. (2015c) in the context of stellar discs. The specification of such externally forced secular dynamics to the case of dynamically degenerate systems is postponed to a future work.

IV.
During its secular evolution along quasistationary equilibria, the dynamics of an isolated system can also be driven by finiteN fluctuations. This amounts to neglecting the effects due to any external stochastic perturbations, and considering the contributions associated with the collisional term in 1 /N in Eq. (51). This requires to solve simultaneously the system of two coupled evolution Eqs. (51)and (52). This approach is presented in Sect. 5, where the analogs of the (bare, that is without collective effects) Landau equation and (dressed, that is with collective effects) BalescuLenard equation are derived in the context of degenerate dynamical systems, such as galactic nuclei. As will be emphasised later on, these diffusion equations, sourced by finiteN fluctuations capture the known mechanism of resonant relaxation (Rauch & Tremaine 1996). See BarOr & Alexander (2014) for a similar study of the effect of the finiteN stochastic internal forcing via the socalled ηformalism.
We note that one could also consider the secular evolution of a nonaxisymmetric set of eccentric orbits orbiting a black hole as an unperturbed collisionless equilibrium (corresponding to the expected configuration of the galactic centre of M 31 Tremaine 1995). The derivation of the associated BalescuLenard equation for such a configuration would first involve identifying new angleaction variables for the nonaxisymmetric configuration so as to satisfy II, and then extend the formalism accordingly. This will not be explored any further in this paper. Regarding item II, we expect that, depending on the relative mass of the considered cluster, there is a regime where the selfinduced orbital precession is significant, but the selfgravity of the wires is not strong enough to induce a collisionless instability. In this regime, accounting for the polarisation of the orbits becomes important in item III and IV. This motivates the rest of the paper.
5. The degenerate BalescuLenard equation
We now show how to obtain the closed kinetic equations – the degenerate BalescuLenard and Landau equations – when considering the 1 /N collisional contribution present in the evolution Eq. (51). It will be assumed that the system is isolated so that it experiences no external perturbations. Our aim is to obtain a closed kinetic equation involving only. To do so, we rely on the adiabatic approximation (or Bogoliubov’s ansatz) that the system secularly relaxes through a series of collisionless equilibria. In this context, collisionless equilibria are stationary (and stable) steady states of the collisionless advection component of Eq. (51). Therefore, it is assumed that throughout the secular evolution, one has (53)As already highlighted, it is expected that such collisionless equilibria are rapidly reached by the system (on a few T_{sec}), through an outofequilibrium mechanism related to violent relaxation. In addition, the symmetry of the system is expected to be such that the collisionless equilibria are of the form (54)so that, during its secular evolution, the system’s averaged DF does not have any slow angle dependence. Notice however that, despite the hypothesis from Eq. (54), the averaged autocorrelation evolving according to Eq. (52)still depends on the two slow angles and . We also assume that the symmetry of the system is such that (55)As we will see later on in Sects. 6.1 and 6.2, such symmetry is satisfied for instance for razorthin axisymmetric discs and 3D spherical clusters (see also Appendix A for the expression of the relativistic precession frequencies). Given Eqs. (54)and (55), the equilibrium condition from Eq. (53)is immediately satisfied. We introduce the precession frequencies Ω^{s} as (56)These frequencies correspond to the precession frequencies of the slow angles due to the joint contributions from the system’s selfconsistent potential and the relativistic corrections. Notice that they do not involve the Keplerian frequencies from Eq. (24)anymore and hence are not degenerate a priori. With them, one can for example easily rewrite the collisionless precession advection term from Eq. (52)as (57)where the precession frequencies associated with the slow angles come into play.
The two coupled evolution Eqs. (51)and (52)are now quasiidentical to the traditional coupled BBGKY equations considered in Heyvaerts (2010) to derive the inhomogeneous BalescuLenard equation for nondegenerate inhomogeneous systems. Various methods have been proposed in the literature to derive the closed kinetic equation satisfied by . Heyvaerts (2010) proposed a direct resolution of the BBGKY equations, based on Bogoliubov’s ansatz. Chavanis (2012) considered a rewriting of Eqs. (51)and (52)using the Klimontovich equation (Klimontovich 1967), and relied on a quasilinear approximation. Finally, in the limit where collective effects are not accounted for, Fouvry et al. (2016a) recently presented a new derivation of the relevant kinetic equation based on functional integrals.
In the present paper, the derivation proposed by Heyvaerts (2010) will be followed, by directly solving the two first averaged BBGKY Eqs. (51)and (52). The basic idea of this approach is to solve Eq. (52), so as to obtain the system’s autocorrelation as a functional of the system’s 1body DF . Injecting this expression in Eq. (51)yields finally a closed kinetic equation quadratic in . The detailed calculations required to derive the inhomogeneous degenerate BalescuLenard equation are presented in Appendix B.
5.1. The one component BalescuLenard equation
In its explicitly conservative form, the degenerate inhomogeneous BalescuLenard equation reads (58)In Eq. (58), we recall that d is the dimension of the physical space and k the number of degeneracies of the underlying zerothorder potential. The r.h.s. of Eq. (58)is the degenerate inhomogeneous BalescuLenard collision operator, which describes the secular diffusion induced by dressed finiteN fluctuations. It describes the distortion of Keplerian orbits as their actions diffuse through their selfinteraction. As expected, it vanishes in the large N limit. Notice the presence of the resonance condition operating on their precession frequencies encapsulated by the Dirac delta (using the shortened notation ), where , are integer vectors. In fact, Eq. (58)shows that the diffusion occurs along preferred discrete directions labelled by the resonance vectors . The integration over the dummy variable J_{2} scans action space for regions where the resonance condition is satisfied, and such resonant (possibly distant) encounters between orbits are the drivers of the collisional evolution. The resonance condition is illustrated in Fig. 1.
Fig. 1 Illustration of the resonance condition appearing in the degenerate inhomogeneous BalescuLenard Eq. (58). Topleft: a set of two resonant orbits precessing at the same frequency ω_{s}. Topright: in the rotating frame at frequency ω_{s} in which the two orbits are in resonance. Bottom: fluctuations of the system’s DF in action space caused by finiteN effects and showing overdensities for the blue and red orbits. The dashed line correspond to the critical resonant line in action space along which the resonance condition Ω^{s}(J) = ω_{s} is satisfied. The two set of orbits satisfy a resonance condition for their precession frequencies, and uncorrelated sequences of such interactions lead to a secular diffusion of the system’s orbital structure following Eq. (58). Such resonances are non local in the sense that the resonant orbits need not be close in action space nor in position space. As emphasised in Sect. 6.1 for axisymmetric razorthin discs, symmetry enforces , so that the two orbits are caught in the same resonance. 

Open with DEXTER 
Notice also that Eq. (58)involves the antisymmetric operator, , which when applied to weighs the relative number of pairwise resonant orbits caught in this resonant configuration. The quantities are the socalled dressed susceptibility coefficients: each distribution entering the r.h.s. of Eq. (58)is boosted by this susceptibility. These dressed coefficients include the effects of the gravitational wake induced by each wire, represented by the last term of Eq. (52); in constrast, bare susceptibility coefficients (introduced later) are obtained without taking this selfgravity into account. In order to solve Poisson’s nonlocal equation relating the DF’s perturbations and the induced potential perturbations, Kalnajs’ matrix method (Kalnajs 1976) can be used to implement a biorthonormal basis of potentials and densities ψ^{(p)} and ρ^{(p)} such that (59)where U stands for the rescaled interaction potential from Eq. (20). The dressed susceptibility coefficients appearing in Eq. (58)are then given by (60)where I is the identity matrix, and is the system’s averaged response matrix defined as (61)In Eq. (61), the averaged basis elements were defined following Eq. (33). Their Fourier transform with respect to the slow angles was also defined using the convention (62)The susceptibility coefficients from Eq. (60)quantify the polarisation cloud around each orbit, which triggers sequences of transient wakes (Julian & Toomre 1966; Toomre 1981). In the secular timeframe, these are assumed to be instantaneous, via the socalled Bogoliubov’s ansatz, as shown in Appendix B.1.
One can straightforwardly rewrite the BalescuLenard Eq. (58)as an anisotropic nonlinear diffusion equation, by introducing the appropriate drift and diffusion coefficients. Equation (58)then reads (63)where and are respectively the drift and diffusion coefficients associated with a given resonance . The secular dependence of these coefficients with the system’s averaged DF, , is not written out explicitly to simplify the notations but is a central feature of the present formalism. In Eq. (63), the drift coefficients and diffusion coefficients are given by (64)When collective effects are not accounted for (that is when the last term of Eq. (52)is neglected), the degenerate BalescuLenard Eq. (58)becomes the degenerate Landau equation (see Polyachenko & Shukhman 1982; Chavanis 2013b, for the nondegenerate case), which reads (65)Notice that this is just the previous BalescuLenard Eq. (58)with the dressed replaced by the bare susceptibility coefficients . The latter are related to the (partial) Fourier transform of the interaction potential (LyndenBell 1994; Pichon 1994; Chavanis 2013b) and read (66)so that the averaged interaction potential from Eq. (44)can be decomposed as (67)One should note that the kinetic Eqs. (58)and (65), while defined on the full action space J = (J^{s},J^{f}), do not allow for changes in the fast actions J^{f}. Indeed, if one defines the marginal DF, , as , Eqs. (58)and (65)immediately give (68)so that the collisional secular diffusion occurs only in the directions J^{f} = const.
5.2. Multiple components black hole environment
It is of prime importance to follow the joint longterm evolution of multiple types of stars or black holes orbiting a central supermassive black hole, as it will allow astronomers to capture their relative segregation, when the lighter black holes sink in towards the more massive one. In turn, this could allow us to predict the expected rate of mergers and accretion events.
As already emphasised in Heyvaerts (2010), Chavanis (2012), the BalescuLenard equation can also be written for a system involving multiple components (corresponding to say, a spectrum of stars and low mass black holes or debris of different masses orbiting the central object). The different components will be indexed by the letters “a” and “b”. The particles of the component “a” have a mass μ_{a} and follow the DF F^{a}. As briefly detailed in Appendix B.3 (which gives the details of all normalisations), the evolution of each DF is given by (69)where the dimensionless relative mass η_{a} = μ_{a}/M_{⋆} was introduced, and where is the total active mass of the system. In the multicomponent case, the dressed susceptibility coefficients are still given by Eq. (60). However, as expected, the response matrix now encompasses all the active components of the system which polarise so that In the limit where only one mass is considered, one has η_{a} = 1 /N_{a}, and the single mass BalescuLenard Eq. (58)is recovered. Equation (69)describes the evolution of the “a” population, and differs from Eq. (58)via the weight η_{a}, and the sum over “b” weighted by η_{b}. As in Eq. (63), one can introduce drift and diffusion coefficients to rewrite Eq. (69)as (70)where the drift and diffusion coefficients and depend on the position in action space J_{1}, the considered resonance , and the component “b” used as the underlying DF to estimate them. The drift coefficients and diffusion coefficients are given by (71)Equation (70)can finally be rewritten as (72)where the total drift and diffusion coefficients and are given by In Eq. (72), the total drift coefficients are multiplied by the dimensionless mass η_{a} of the considered component. This essentially captures the known process of segregation, when a spectrum of masses is involved, so that components with larger individual masses tend to narrower steady states. Indeed, the multicomponent BalescuLenard formalism captures the secular effect of multiple resonant (nonlocal) deflections of lighter particles by the more massive ones: the lighter population will drift towards larger radii, while the massive one will sink in. This can be seen for instance by seeking asymptotic stationary solutions to Eq. (72)by nulling the curly brace in its r.h.s.
5.3. Secular evolution increases Boltzmann entropy
Following closely the demonstration presented in Heyvaerts (2010), let us define the system’s entropy S(τ) as (73)Differentiating Eq. (73)once with respect to τ yields (74)Let us introduce the system’s diffusion flux, ℱ_{tot}(J_{1}), given by (75)with given by (76)such that Eq. (58)reads (77)Using integration by parts in Eq. (74)and ignoring boundary terms leads to (78)Given Eq. (75), Eq. (78)can be rewritten as with , and . This equation can symmetrised via the substitutions and J_{1} ↔ J_{2}, relying on the fact that , so that (79)As the entropy function satisfies s′′(x) = 1 /x (any double primitive of 1 /x would work too), the square braket of Eq. (79)can immediately be factored as (80)so that one finally gets dS/ dτ ≥ 0. This entropy increase corresponds to heat generation as the orbital structure of the cluster rearranges itself in a more eccentric configuration. The previous demonstration naturally extends for the multicomponent BalescuLenard Eq. (69). Indeed, defining the system’s total entropy S_{tot}, summed for all components, as (81)one can again show that for s′′(x) = 1 /x, one has dS_{tot}/ dτ ≥ 0, which does not necessarily imply that the entropy of each component increases.
6. Applications
Up to now we have considered the general framework of a system made of a finite number of particles orbiting a central massive object. We now examine in turn some more specific configurations of particles orbiting a black hole, and discuss how the results of the previous section can be further extended when considering specific geometries and physical secular processes, to highlight the wealth of possible implications one can draw from this framework. Detailed applications are postponed to followup papers.
6.1. Razorthin axisymmetric discs
Let us first specialise the degenerate BalescuLenard Eq. (58)to razorthin axisymmetric discs. For such systems, the dimension of the physical space is given by d = 2, while the number of dynamical degeneracies of the Keplerian dynamics is k = 1. Therefore, the resonance condition in Eq. (58)takes the simpler form of a 1D condition naively reading and the Delaunay angleaction variables from Eq. (26)become (82)Symmetries of the interaction potential lead to relationships among the susceptibility coefficients, which simplify the BalescuLenard equation. The rescaled interaction potential U_{12} from Eq. (20)takes the form (83)in which we introduce the usual polar coordinates (R,φ). Following Eqs. (3.28a) and (3.28b) of Binney & Tremaine (2008), the mapping from the physical polar coordinates to the Delaunay angleaction ones can be written as (84)where the semimajor axis a, eccentricity e, true anomaly f, and eccentric anomaly η are introduced as (85)Substituting Eq. (84)into Eq. (83), we immediately have that (86)As a consequence, the bare susceptibility coefficients from Eq. (66)for a razorthin disc are related to one another through (87)A similar result also holds for the dressed susceptibility coefficients from Eq. (60). Indeed, for any 2D razorthin system, one can assume the basis elements from Eq. (59)to be generically of the form (88)where ℓ^{p} and n^{p} are two integer indices, and are radial functions. Such a decomposition of the basis elements allows us to decouple the azimuthal and radial dependence of the basis elements. Noting that in the mapping from Eq. (84)only the azimuthal angle φ depends on the slow angle g, one obtains that the Fourier transformed basis elements satisfy (89)Substituting this into the expression (61)for the response matrix, we find that (90)Using Eqs. (89)and (90), the dressed susceptibility coefficients satisfy (91)just as the bare ones satisfy Eq. (87).
The symmetry corresponding to Eq. (91)allows us to get rid of the sum over the resonance index in the BalescuLenard Eq. (58), so that it becomes (92)in which we use the relation δ_{D}(αx) = δ_{D}(x)/  α  and introduce the (unique) total dressed susceptibility coefficient (93)Similarly, if we neglect selfgravity, then the symmetry in Eq. (87)applies and Eq. (92)becomes the associated Landau equation, in which the total dressed susceptibility coefficient is replaced by the bare one, (94)This Landau analog of Eq. (92)for razorthin axisymmetric discs with the bare susceptibility coefficients from Eq. (94)has already been derived in Sridhar & Touma (2017) via Gilbert’s equation.
The result of these simplifications is that the degenerate BalescuLenard Eq. (92)possesses a straightforward resonance condition in which resonant encounters can only occur between two orbits caught in the same resonance, as illustrated in Fig. 1. To compute the diffusion flux appearing in the r.h.s. of this equation, we employ the generic definition of the composition of a Dirac delta and a function (Hörmander 2003), which in a ddimensional setup takes the form (95)where g^{1}(0) = {x  g(x) = 0} is the hypersurface of dimension (d−1) defined by the constraint g(x) = 0, and dσ(x) is the surface measure on g^{1}(0). In our case, the resonance condition is given by the function (96)For a given value of J_{1}, and introducing ω = Ω^{s}(J_{1}), we define the critical resonant curve γ(ω) as (97)This curve corresponds to the set of all orbits which are in resonance with the precessing orbit of action J_{1}. Once this resonance line has been identified, the diffusion flux from Eq. (92)is straightforward to compute and reads (98)where to shorten the notations, we introduced the function G(J_{1},J_{2}) as (99)as well as the resonant contribution  ∇(Ω^{s}(J_{2}))  given by (100)We note that Eq. (98)is now a simple onedimensional integral involving a regular integrand.
In summary, because the quasistationary potentials and are known via Eqs. (43)and (A.8), one can compute the associated precession frequencies Ω^{s} (and their gradients). This allows for the determination of the critical resonant lines γ from Eq. (97). Following Eq. (98), it then only remains to integrate along these lines to determine the secular diffusion flux. Such an effective computation for razorthin discs in the Landau limit is postponed to a followup paper, as Eq. (43)involves a singular triple integral over wirewire interactions. Similarly, the study of the longterm evolution of quasistationary nonaxisymmetric razorthin discs (such as M 31) will also be the subject of a future work.
6.2. Spherical cluster around BH
We now turn to the application of the degenerate BalescuLenard Eq. (58)to spherically symmetric systems. The general procedure is the same as for the razorthin disc case, but now the dimension of the physical space is d = 3, while the number of Keplerian dynamical degeneracies is given by k = 2. The resonance condition in Eq. (58)becomes twodimensional and reads . In the 3D context, the Delaunay variables from Eq. (82)become (101)where g stands for the angle from the ascending node to the periapse, h for the longitude of the ascending node and w for the Keplerian orbital phase, that is the mean anomaly.
As was done in Eq. (92)for razorthin discs, let us now show how the 3D geometry allows us to further simplify the kinetic equation. Written in spherical coordinates (R,θ,φ), the rescaled interaction potential from Eq. (20)becomes (102)Following Eq. (5.20) from Merritt (2015), these (R,θ,φ) can be expressed as a function of the Delaunay angleaction variables, so that (103)where a, e, f and η were introduced in Eq. (85), and i is the orbit’s inclination, defined through cos(i) = L_{z}/L. Therefore the interaction potential of Eq. (102)and its angleaveraged version (Eq. (44)) have the symmetries (104)From this, it immediately follows that the bare susceptibility coefficients (Eq. (66)) are related to one another via (105)Here, we have written the resonance vectors as , so that the coefficient is the one associated with the slow angle h. A similar result holds for the dressed susceptibility coefficients defined in Eq. (60). Indeed, for any 3D system, the basis elements in Eq. (59)can be written as (106)where ℓ^{p}, m^{p} and n^{p} are three integer indices, are the usual spherical harmonics, and are radial functions. We note in the mappings from Eq. (103)that only the azimuthal angle φ depends on the slow angle h. Because the spherical harmonics are of the form , where are the associated Legendre polynomials, one immediately finds that the Fourier transformed basis elements satisfy (107)As a consequence, the expression (61)of the response matrix immediately gives (108)Equations (107)and (108)allow us to rewrite the dressed susceptibility coefficients as (109)showing that they are related to one another in the same way as the bare susceptibility coefficients of Eq. (105).
As a consequence, when considering a spherically symmetric system, one can simplify the resonance condition, and the BalescuLenard Eq. (58)becomes (110)To neglect collective effects in Eq. (110), one only has to make the substitution . Here, it is important to note that the 1.5PN relativistic precession frequencies obtained in Appendix A do depend on the action L_{z}, so that at this stage further simplifications of Eq. (110)are not possible. The computation of the diffusion flux in Eq. (110)proceeds as in Eq. (98)by identifying the critical surfaces of resonance. We do not detail these calculations here.
6.3. Relativistic barrier crossing in the vicinity of BHs
We now show how the degenerate BalescuLenard Eq. (58)naturally accounts for the presence of the “Schwarzschild barrier” encountered by stars diffusing towards the central BH. This Schwarzschild barrier was discovered by Merritt et al. (2011) in their simulations of spherically symmetric star clusters. Here, we show how it arises in the simpler case of a razorthin axisymmetric disc of stars around the BH, but the same fundamental idea applies to the 3D case. The secular collisional evolution of such a disc is governed by Eq. (92). The resonance condition in that equation is , in which the precession frequency Ω^{s}(J) of each of the two wires, defined in Eq. (56), is composed of two parts. The first is the contribution from the system’s selfconsistent Newtonian potential, (111)The second is the additional contribution from relativistic effects. We derive it in Appendix A. In the case of a razorthin disc, it reads (112)We now study how these precession frequencies depend on distance to the central BH. Following the timescale comparisons of Kocsis & Tremaine (2011), one expects the relativistic precession frequency to dominate close to BH (and, in fact, to diverge as the star gets closer to capture), while the selfconsistent one, , will be the largest for orbits in the vicinity of the considered disc. Such a behaviour is qualitatively illustrated in Fig. 2, where we represent the typical dependence of the precession frequencies as a function of the distance to the central BH.
Fig. 2 Illustration of the typical dependence of the precession frequencies and (Eqs. (111)and (112)) as a function of the distance to the central BH. The relativistic precession frequencies, diverge as the star gets closer to the central BH, while the selfconsistent ones are typically the largest for stars in the neighbourhood of the considered disc. The black dots give all the locations, whose precession frequency is equal to ω_{s} (illustrated by the dotted horizontal line). These positions are in resonance and will therefore have a nonvanishing contribution in the BalescuLenard Eq. (92). Because Eq. (92)involves the product of the system’s DF in the two resonating locations, the resonant coupling between the two outer points (which belong to the region where the disc dominates) will be much stronger than the couplings involving the inner point (which does not belong to core of the disc). As stars move inward, their precession frequencies increase up to a point where this prevents any resonant coupling with the disc’s region. This effectively stops the secular diffusion, and induces a diffusion barrier. 

Open with DEXTER 
Figure 2 shows that, for a given precession frequency ω_{s}, one can identify the actions J within the disc for which the resonance condition is satisfied.
Equation (92)involves the quadratic factor , which is the product of the system’s density at the two locations that are in resonance. As shown in Fig. 2, because the disc is only located in the outer regions of the BH, the resonant coupling between two locations within the disc will be much stronger than one involving a resonant location inside the inner edge of the disc, very close to the BH. Therefore, in Fig. 2, the coupling between the two outer black dots will be much larger than the couplings involving the inner dot. The situation becomes even worse if one wants to couple a region even closer to the BH, for which the precession frequency is too large to resonate with any part of the disc. In this situation, no efficient resonant couplings are possible and the secular diffusion is drastically suppressed. In short, the divergence of the relativistic precession frequencies in the neighbourhood of the BH means that stars whose orbits diffuse inwards closer to the BH experience a rise in their precession frequency, which prevents them from resonating anymore with the disc, strongly suppressing further inward diffusion. This is the socalled Schwarzschild barrier.
This explanation of the Schwarzschild barrier using the notion of resonant coupling is directly related to the explanation proposed in BarOr & Alexander (2014), which relies on the concept of adiabatic invariance. In their picture, a test star can undergo resonant relaxation only if the timescale associated with its relativistic precession is longer than the coherence time of the perturbations induced by the field stars and felt by the test star. Because the typical coherence time of the perturbations scales as the inverse of the typical precession frequency of the field stars (which lie within the cluster), the requirement for an efficient diffusion from the adiabatic invariance point of view is equivalent to the requirement from the point of view of the BalescuLenard resonance condition.
6.4. Solving the BalescuLenard equation by Monte Carlo sampling
This suppression of diffusion in the neigbourhood of the BH can also be illustrated by considering the orbitaveraged motion of individual wires. Equation (92)takes the form of a diffusion equation in action space, where one follows selfconsistently the evolution of the system’s DF. One could also be interested in describing the stochastic evolution of individual stellar wires, whose ensemble average is described by this diffusion equation. To do so, let us rewrite Eq. (92)as (113)Then, as detailed in Appendix C, one can write the corresponding Langevin equation that captures the dynamics of individual test wires. Here, we consider the case of a razorthin axisymmetric disc and denote as the position at time τ of a test wire in action space J = (L,I). Following Eq. (C.5), the dynamics of the test wire takes the form (114)in which the 1D Langevin coefficients that describe the diffusion of the wire in the Ldirection are given by (115)and the Langevin stochastic force Γ(τ) satisfies Eq. (C.6). As in Eq. (68), the individual fast action J^{f} = I is preserved during the wire’s evolution.
Equation (114)describes the diffusion of an individual test wire when embedded in the selfinduced noisy environment that is described by the drift and diffusion coefficients from Eq. (92). As such it could be used iteratively jointly with Eq. (93)– which depends on the sampled position of all orbits in action space via Eq. (61)– to effectively integrate Eq. (92)over cosmic time. This would simply involve discretising Eq. (114)in time as ℒ_{i + 1} = ℒ_{i} + (dℒ/dτ)_{i}Δτ, while sampling initial ℒ_{0}s to match the original distribution. Strikingly, Eq. (114)shares some similarity with the individual Hamilton’s equations associated with the Hamiltonian from Eq. (7), but the significant gain of the present work is to allow for individual timesteps, Δτ, which are orders of magnitudes larger than the original one required to solve for the trajectories of individual stars. It also deals seamlessly with postNewtonian orbit integration over the fast and slow angles.
A qualitative description of the dynamics of individual orbits from Eq. (114)is illustrated in Fig. 3.
Fig. 3 Illustration of the individual dynamics of stars in the (j,a) = (L/I,I^{2}/GM_{•}) space, as given by the Langevin Eq. (114). The grey region corresponds to the capture region, within which stars inevitably sink into the BH. As I is an invariant of the diffusion (see Eq. (114)), stars’ diffusion is onedimensional, conserves a, and occurs only in the jdirection. The background dotted curves illustrate the contour lines of the precession frequency given by the function (j,a) → Ω^{s}(j,a). As illustrated in Fig. 2, the precession frequencies get larger as particles approach the central BH, due to the contributions from the relativistic precession frequencies. The blue and red orbits precess at the same frequency ω_{s}, so that they belong to the same critical resonant line γ_{ωs}, which allows them to resonate one with another. As the precession frequencies diverge in the vicinity of the BH, such resonant couplings are significantly less likely as stars get closer to the BH, which effectively creates a diffusion barrier in action space, the socalled Schwarzschild barrier. 

Open with DEXTER 
Following the representations from BarOr & Alexander (2016), Fig. 3 represents the diffusion of stars in the (j,a) = (L/I,I^{2}/ (GM_{•})) space. As observed in Eq. (114), the fast action I of the stars is conserved during the diffusion, so that stars diffuse only in the jdirection along a = const. lines. When diffusing, individual particles may resonate with stars which precess at the same frequency, such as the blue and red particles in Fig. 3. However, as already illustrated in Fig. 2, the precession frequencies diverge as stars get closer to the BH. This increase in the precession frequencies will then forbid any resonant coupling between a star in this internal region and stars belonging the disc itself, where precession frequencies are much smaller. Resonances becoming impossible, the diffusion is stopped and stars cannot keep diffusing closer to the central BH. This suppression of the diffusion is the Schwarzschild barrier. A quantitative illustration of this damping of resonant couplings is postponed to a later paper, where we will effectively compute the precession frequencies in action space for a physically motivated razorthin axisymmetric disc.
6.5. Evolution of BH mass and spin
The calculations above successfully explain the existence of the socalled Schwarzschild barrier, which strongly suppresses the supply of tightly bound matter to the black hole. We note that the analysis of Merritt et al. (2011) suggests that, in practice, this suppression is probably tempered by simple twobody relaxation (not accounted for in the orbitaveraged approach followed in this paper), which provides an additional mechanism for transporting stars even closer to the BH, once resonant relaxation becomes inefficient. This mechanism was recently demonstrated in detail in BarOr & Alexander (2016), which showed that adiabatic invariance (in other words the damping of resonant relaxation) limits the effects of resonant relaxation to a region well away of the loss lines, so that the dynamics of accretion of stars by the BH is only very moderately affected by the presence of resonances. Nevertheless, one can calculate the rate at which stars are transported across any boundary in phase space within which resonant relaxation dominates, which is important for quantifying the growth rate of the central black hole. Consider then a fixed boundary in action space. From the divergence theorem, the flux of mass, dM/ dτ through that boundary, , due to secular diffusion is proportional to (116)where n is the exterior pointing normal vector. In Eq. (116), one can note that the contribution from a given resonance m^{s} takes the form of a preferential diffusion in the direction of m^{s}. This diffusion is therefore anisotropic because it is maximum for n ∝ m^{s} and equal to 0 for n·m^{s} = 0. We note that if a set of stars of various masses or black holes orbit the galactic centre, the net flux of each component can also be computed via Eq. (70)as This is likely to be of particular interest for predicting the distribution of heavy compact remnants, which, from equipartition arguments, are expected to sink more rapidly towards the centre. Similarly, the flux of angular momentum, dL/ dτ, can be computed, and is proportional to (117)and could contribute to either spinning up or down the central black hole, once the selfconsistent evolution of the black hole’s spin and loss of angular momentum via gravitational wave emission are taken care of. We note that if the disc is sufficiently selfgravitating, the diffusion in action space is likely to be dominated by a specific resonance (as was shown in Fouvry et al. 2015b).
7. Discussion and conclusion
Supermassive black holes absorb stars and debris whose orbits reach the losscone, the region of phase space corresponding to orbits on which they are either taken directly into the black hole or close enough to interact strongly with it. Such accretion affects the secular evolution of the SMBH’s mass and spin, which is of interest in understanding black hole demographics and AGN feedback (Volonteri et al. 2016). It also affects the matter that remains. For instance, the continuous loss of stars can resupply and reshape the central stellar distribution (e.g. Genzel et al. 2000). These dynamical processes have observable signatures, such as binary capture and hypervelocity star ejection (Hills 1988), the tidal heating and disruption of stars (Frank & Rees 1976), gravitational waves produced by inspiraling compact remnants (Abbott et al. 2016). All these signatures provide possible indirect evidence of the existence of the black hole and offer the opportunity of probing the theory of relativity in the strong field limit (Blanchet 2014). Understanding the dynamics of stars in the vicinity of supermassive black holes is in fact one of the prime goal of the new generation of interferometers such as Gravity (Jocou et al. 2014).
In this paper, we have specialised the recently developed kinetic theory of selfgravitating systems of N particles (Heyvaerts 2010; Chavanis 2012) to quasiKeplerian systems dominated by a massive central object, deriving the equation that governs the secular evolution of such systems to leading order in 1 /N. The selfconsistent dressed equations (Eq. (58)and its multicomponent and stochastic counterparts, Eqs. (69)and (C.5)respectively) account for the dynamical degeneracies in quasiKeplerian systems. Because purely Keplerian orbits do not precess, the dynamical evolution of such degenerate systems may differ significantly from that of fully selfgravitating systems, such as discs and spheroids. In particular, to a good approximation stars behave as if they were smeared out into orbitaveraged Keplerian wires and the evolution of the system modelled by following the dressed interactions among such wires. The coupling among these wires generates sequences of uncorrelated transient polarised density waves, which make the underlying stars’ orbits diffuse in phase space.
The quasiKeplerian BalescuLenard Eq. (58)is quadratic in the phaseaveraged distribution function and describes i) the selfgravity of the orbiting particles; ii) the discreteness of the cluster; iii) the resonances between such orbits; iv) a full spectrum of masses, via Eq. (69); and v) possible postNewtonian corrections, including relativistic precession induced by the rotation of the central black hole, if present. These last effects are encoded in the frequency shifts occurring in the resonance condition from the diffusion and drift coefficients. It is therefore the quasilinear selfconsistent master equation quantifying the effect of resonant relaxation. As such it provides a very rich framework to describe the evolution of galactic centres for cosmic times, or the secular evolution of debris discs – which is an interesting venue in the context of planet formation (e.g. Tremaine 1998).
A key step in the derivation of this equation is the phase averaging of the first two equations of the BBGKY hierarchy over the fast angles associated with the orbital motion of the bodies on their Keplerian orbits. In order to derive Eqs. (58)and (69), we assumed that the (spherical or coplanar) cluster was dynamically relaxed at every stage of its secular evolution. As the equations are averaged over the Keplerian fast angles, the corresponding actions are adiabatically preserved. Because of this phase average, the Keplerian BalescuLenard equation cannot capture mean motion resonances. Hence a limitation of the present formalism is that it is restricted to systems with a high degree of symmetry.
More generally, the averaging over fast angles means that traditional nonresonant twobody relaxation is not accounted for in the BalescuLenard equations we derive here. This is usually appropriate though, because the derivation of these equations ignores terms of order O(1 /N^{2}), which means that they are valid only on timescales ≲ Nt_{d}, where t_{d} is the dynamical timescale. Such timescales are typically expected to be much shorter than the nonresonant twobody relaxation time. When investigating specifically the vicinity of supermassive black holes, we found that the quasiKeplerian BalescuLenard equation captures naturally the presence of a Schwarzschild barrier, explains why it is not fully impermeable, and why it allows us to estimate for instance the mass and angular momentum fluxes of each component through its boundary. In its multicomponent formulation, the BalescuLenard equation also captures mass segregation and radial migration as entropy increases.
7.1. Comparison to other work
A number of other recent papers have tackled the dynamics of quasiKeplerian stellar systems. The closest to the present paper is the recent sequence of papers by Sridhar & Touma (2016a,b), who have already obtained equations equivalent to our Eqs. (51)and (52)following a different route starting from the approach of Gilbert (1968), which itself extended the work of Balescu (1960), Lenard (1960) from plasma physics. The “passive response” approximation they make in their analysis of razorthin axisymmetric discs (Sridhar & Touma 2017) corresponds to the Landau limit in which one uses the bare susceptibility coefficients from Eq. (94)in the BalescuLenard Eq. (92).
Another way of modelling such dynamics is by using some form of Monte Carlo approach in which the noise due to the discrete number of stars is treated as an externally imposed perturbation (e.g. Madigan et al. 2011; BarOr & Alexander 2014). This basic idea is very powerful, particularly if one wants to investigate additional perturbations that are genuinely external to the cluster. For example, the ηformalism introduced recently in BarOr & Alexander (2014) and implemented in detail in BarOr & Alexander (2016) is one such scheme. Imposing plausible constraints on the power spectrum of the discreteness noise, these papers recovered the location of the Schwarzschild barrier (explained in terms of adiabatic invariance), and investigated the role of 2body relaxation for the losscone problem. They showed that on longer timescales, 2body nonresonant relaxation completely erases the Schwarzschild barrier, and also argued that resonant relaxation is effective only in a restricted region of action space away from the losslines, so that its overall effect on plunge rates is small.
The BalescuLenard equation has a couple of important conceptual advantages over the ηformalism (and similar Monte Carlo schemes). First, the ηformalism requires assumptions about the statistical characteristics of the externally imposed discreteness noise felt by each wire. The BalescuLenard equation requires no such external input, because the system’s discreteness is described selfconsistently. Second, in the ηformalism the selfgravity of the response to the noise is difficult to account for. In the BalescuLenard equation, this full response is naturally present in the dressed susceptibility coefficients (Eq. (60)). Such collective effects can be crucial in systems close to marginal stability, where the associated polarisation can get very large (see e.g. Fouvry et al. 2015b, for an illustration in the case of razorthin stellar discs). We note that, just as in Monte Carlo schemes, the BalescuLenard approach also offers a natural way of including external potential fluctuations (see point III of Sect. 4).
At the heart of the ηformalism lies a distinction between “field” and “test” stars: the dynamics of the test stars are followed as they undergo the stochastic perturbations generated by the field stars. Such a split is also used in the restricted Nbody calculations recently presented in Hamers et al. (2014): the motion of each field star is followed along their precessing Keplerian orbits (with a precession induced by both relativistic effects and the system’s selfconsistent potential), but interactions among field stars are ignored. The test stars are then followed by direct integration of their motion in the timevarying potential due to the field stars: it does not rely on the averaging approximation. Such an approach is especially useful in order to get a better grasp of the typical stochastic perturbations generated by the cluster of field stars. Like the ηformalism, it ignores the interactions among field stars (and indeed among test stars) and there is no backinfluence of the test stars on the field ones.
7.2. Future work
The Langevin formulation of the BalescuLenard equation (Sect. 6.4 and Appendix C) combines the flexibility of Monte Carlo methods with a selfconsistent treatment of the dynamics. A subsequent improvement is offered by the possibility of adding the secondary effects of twobody relaxation and gravitationalwave losses to the resonant relaxation dynamics, on which the present paper was focused. Eventually, one could evolve jointly the BH and its environment. This would involve considering that the frequencies in Eq. (58)are time dependent, via the variation of the BH’s mass and spin as outlined in Sect. 6.3. In the context of razorthin axisymmetric discs, one would compute the drift and diffusion coefficients given by Eq. (92), following the steps of Fouvry et al. (2015b) transposed to the Keplerian framework. A few difficulties have to be overcome to perform such a computation. The first is the computation of the wirewire interaction potential (see e.g. Touma et al. 2009; Touma & Sridhar 2012) and its harmonic transform over the slow angles. Then, in order to account for the system’s selfgravity, one has to compute the system’s averaged response matrix from Eq. (61), which asks for the integration of a resonant function over action space, a daunting numerical task. Finally, on secular timescales, one has to deal with the selfconsistency of the diffusion, so that the system’s drift and diffusion coefficients should be updated along the diffusion. As was shown in the present paper, the net effect of Eq. (92)will be to induce diffusion along preferred ridges, whose properties are set by the distribution of stars within the cluster and their selfgravity. Depending on their starting point in action space, some orbits will be driven near the region where the black hole dominates diffusion. This will allow us to quantify for instance the relative importance of black hole spin on barrier crossing, and the efficiency at which a supermassive black hole is fed by its surrounding stellar cluster (as discussed in Sect. 6.3).
Acknowledgments
J.B.F. and C.P. thank the CNRS Inphyniti programme for funding and the KIAS for hospitality while this project was initiated. We thank PierreHenri Chavanis and Simon Prunet for comments. Special thanks to Walter Dehnen for his suggestion to consider democratic heliocentric coordinates. C.P. thanks Scott Tremaine for stimulating discussions on how to tackle this problem. This work has made use of the Horizon cluster, hosted by the Institut d’Astrophysique de Paris. Special thanks to Stephane Rouberol for customising the cluster for our purposes. This work is partially supported by the Spin(e) grants ANR13BS050005 of the French Agence Nationale de la Recherche (http://cosmicorigin.org).
References
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, Phys. Rev. Lett., 116, 061102 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Balescu, R. 1960, Phys. Fluids, 3, 52 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 BarOr, B., & Alexander, T. 2014, Class. Quant. Grav., 31, 244003 [NASA ADS] [CrossRef] [Google Scholar]
 BarOr, B., & Alexander, T. 2016, ApJ, 820, 129 [NASA ADS] [CrossRef] [Google Scholar]
 Binney, J., & Tremaine, S. 2008, in Galactic Dynamics, 2nd edn. (Princeton University Press) [Google Scholar]
 Blanchet, L. 2014, Liv. Rev. Relat., 17, 2 [CrossRef] [Google Scholar]
 Born, M. 1960, in The Mechanics of the Atom (F. Ungar Pub. Co.) [Google Scholar]
 Chavanis, P.H. 2010, J. Stat. Mech., 5, 19 [Google Scholar]
 Chavanis, P.H. 2012, Physica A, 391, 3680 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Chavanis, P.H. 2013a, Eur. Phys. J. Plus, 128, 126 [CrossRef] [EDP Sciences] [Google Scholar]
 Chavanis, P.H. 2013b, A&A, 556, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Duncan, M. J., Levison, H. F., & Lee, M. H. 1998, AJ, 116, 2067 [NASA ADS] [CrossRef] [Google Scholar]
 Fouvry, J.B., Pichon, C., & Chavanis, P.H. 2015a, A&A, 581, A139 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fouvry, J.B., Pichon, C., Magorrian, J., & Chavanis, P.H. 2015b, A&A, 584, A129 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fouvry, J.B., Pichon, C., & Prunet, S. 2015c, MNRAS, 449, 1967 [NASA ADS] [CrossRef] [Google Scholar]
 Fouvry, J.B., Chavanis, P.H., & Pichon, C. 2016a, Physica A, 459, 117 [NASA ADS] [CrossRef] [Google Scholar]
 Fouvry, J.B., Pichon, C., Chavanis, P.H., & Monk, L. 2016b, MNRAS, submitted [Google Scholar]
 Frank, J., & Rees, M. J. 1976, MNRAS, 176, 633 [NASA ADS] [CrossRef] [Google Scholar]
 Genzel, R., Pichon, C., Eckart, A., Gerhard, O. E., & Ott, T. 2000, MNRAS, 317, 348 [NASA ADS] [CrossRef] [Google Scholar]
 Gilbert, I. H. 1968, ApJ, 152, 1043 [NASA ADS] [CrossRef] [Google Scholar]
 Goldstein, H. 1950, in Classical mechanics (AddisonWesley) [Google Scholar]
 Hamers, A. S., Portegies Zwart, S. F., & Merritt, D. 2014, MNRAS, 443, 355 [NASA ADS] [CrossRef] [Google Scholar]
 Heyvaerts, J. 2010, MNRAS, 407, 355 [NASA ADS] [CrossRef] [Google Scholar]
 Hills, J. G. 1988, Nature, 331, 687 [NASA ADS] [CrossRef] [Google Scholar]
 Hopman, C., & Alexander, T. 2006, ApJ, 645, 1152 [NASA ADS] [CrossRef] [Google Scholar]
 Hörmander, L. 2003, The analysis of linear partial differential operators. I (SpringerVerlag) [Google Scholar]
 Jalali, M. A., & Tremaine, S. 2012, MNRAS, 421, 2368 [NASA ADS] [CrossRef] [Google Scholar]
 Jocou, L., Perraut, K., Moulin, T., et al. 2014, in Optical and Infrared Interferometry IV, 9146, 91461 [Google Scholar]
 Julian, W. H., & Toomre, A. 1966, ApJ, 146, 810 [NASA ADS] [CrossRef] [Google Scholar]
 Kalnajs, A. J. 1976, ApJ, 205, 745 [NASA ADS] [CrossRef] [Google Scholar]
 Klimontovich, I. 1967, in The statistical theory of nonequilibrium processes in a plasma (M.I.T. Press) [Google Scholar]
 Kocsis, B., & Tremaine, S. 2011, MNRAS, 412, 187 [NASA ADS] [CrossRef] [Google Scholar]
 Kocsis, B., & Tremaine, S. 2015, MNRAS, 448, 3265 [NASA ADS] [CrossRef] [Google Scholar]
 Lenard, A. 1960, Ann. Phys., 10, 390 [NASA ADS] [CrossRef] [Google Scholar]
 LyndenBell, D. 1967, MNRAS, 136, 101 [NASA ADS] [CrossRef] [Google Scholar]
 LyndenBell, D. 1994, in Lectures on stellar dynamics (Berlin: Springer Verlag) [Google Scholar]
 Madigan, A.M., Hopman, C., & Levin, Y. 2011, ApJ, 738, 99 [NASA ADS] [CrossRef] [Google Scholar]
 Merritt, D. 2015, in Astrophysical Black Holes (Springer) [Google Scholar]
 Merritt, D., Alexander, T., Mikkola, S., & Will, C. M. 2011, Phys. Rev. D, 84, 044024 [NASA ADS] [CrossRef] [Google Scholar]
 Morbidelli, A. 2002, in Modern Celestial Mechanics (Taylor & Francis) [Google Scholar]
 Pichon, C. 1994, Ph.D. Thesis, University of Cambridge [Google Scholar]
 Polyachenko, V. L., & Shukhman, I. G. 1982, Soviet Ast., 26, 140 [NASA ADS] [Google Scholar]
 Polyachenko, E. V., Polyachenko, V. L., & Shukhman, I. G. 2007, MNRAS, 379, 573 [NASA ADS] [CrossRef] [Google Scholar]
 Rauch, K. P., & Ingalls, B. 1998, MNRAS, 299, 1231 [NASA ADS] [CrossRef] [Google Scholar]
 Rauch, K. P., & Tremaine, S. 1996, New Astron., 1, 149 [NASA ADS] [CrossRef] [Google Scholar]
 Risken, H. 1996, in The FokkerPlanck Equation (Berlin, Heidelberg: Springer) [Google Scholar]
 Sridhar, S., & Touma, J. 1999, MNRAS, 303, 483 [NASA ADS] [CrossRef] [Google Scholar]
 Sridhar, S., & Touma, J. R. 2016a, MNRAS, 458, 4129 [NASA ADS] [CrossRef] [Google Scholar]
 Sridhar, S., & Touma, J. R. 2016b, MNRAS, 458, 4143 [NASA ADS] [CrossRef] [Google Scholar]
 Sridhar, S., & Touma, J. R. 2017, MNRAS, 465, 1856 [NASA ADS] [CrossRef] [Google Scholar]
 Toomre, A. 1981, in Structure and Evolution of Normal Galaxies, eds. S. M. Fall, & D. LyndenBell, 111 [Google Scholar]
 Touma, J. R., & Sridhar, S. 2012, MNRAS, 423, 2083 [NASA ADS] [CrossRef] [Google Scholar]
 Touma, J. R., Tremaine, S., & Kazandjian, M. V. 2009, MNRAS, 394, 1085 [NASA ADS] [CrossRef] [Google Scholar]
 Tremaine, S. 1995, AJ, 110, 628 [NASA ADS] [CrossRef] [Google Scholar]
 Tremaine, S. 1998, AJ, 116, 2015 [NASA ADS] [CrossRef] [Google Scholar]
 Tremaine, S. 2005, ApJ, 625, 143 [NASA ADS] [CrossRef] [Google Scholar]
 Volonteri, M., Dubois, Y., Pichon, C., & Devriendt, J. 2016, MNRAS, 460, 2979 [NASA ADS] [CrossRef] [Google Scholar]
 Yu, Q. 2002, MNRAS, 331, 935 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Relativistic precessions
Let us briefly detail the relativistic precessions encompassed in particular by the averaged potential correction present in Eqs. (51)and (52). In order to obtain explicit expressions for these corrections, the 3D Delaunay variables from Eq. (26)will be used. In addition, we assume here that the spin of the BH is aligned with the zdirection and introduce the BH’s spin parameter 0 ≤ s ≤ 1. To recover the expression of these precession frequencies, let us follow Merritt (2015). Equation (5.103) therein gives us that during one Keplerian orbit of duration T_{Kep} = 2π/ Ω_{Kep} = 2πI^{3}/ (GM_{•})^{2}, the 1PN Schwarzschild precession effect leads to a modification of the slow angle g given by (A.1)This is straightforwardly associated with a precession frequency given by (A.2)where the semimajor axis a and the eccentricity e satisfy a = I^{2}/ (GM_{•}) and 1−e^{2} = (L/I)^{2}. We also introduced the Hamiltonian as (A.3)Similarly, Eq. (5.118) of Merritt (2015) gives that the 1.5PN LenseThirring precession during one Keplerian orbit leads to a precession of the slow angle g given by (A.4)where it is assumed that the spin of the BH was aligned with the zdirection. This is immediately associated with a precession frequency given by (A.5)relying on the relation L_{z} = Lcos(i). Hence the Hamiltonian which accounts for the rotation of the BH reads (A.6)Such a Hamiltonian also induces relativistic precessions with respect to the second slow angle h associated with L_{z}. We do not detail here how these precessions are indeed correctly described by the Hamiltonian . Paying a careful attention to the normalisation prefactors used in Eqs. (2), (21), and (45), one finally gets the expression of the averaged 1PN and 1.5PN relativistic corrections appearing in Eqs. (51)and (52). These read (A.7)From this potential correction, following Eq. (56), one can immediately compute the associated precession frequencies with respect to the slow angles θ^{s}. They read (A.8)Note finally that gravitational waves emissions, along with the associated dissipations, are not considered here.
Appendix B: The degenerate collisional equation
For completeness, let us revisit the derivation of the BalescuLenard equation presented in Heyvaerts (2010) in this new quasiKeplerian regime. The starting point of this derivation is the two coupled averaged Eqs. (51)and (52), which involve the system’s averaged 1body DF , and its averaged 2body autocorrelation . The heart of the present derivation is the following: first one must solve the evolution Eq. (52), so as to obtain (Sect. B.1). Injecting this relation in Eq. (51), one obtains a closed kinetic equation involving only. Its simplification will be carried out in Section B.2. Section B.3 will present the specifics of the corresponding multicomponent derivation.
Appendix B.1: Solving for the autocorrelation
With the assumption of stationarity from Eq. (53)and the Bogoliubov’s ansatz from Eq. (54), one can rewrite Eq. (51)as (B.1)where the collision operator C is introduced as (B.2)relying on the fact that during the secular diffusion, is of the form , allowing us to perform an angleaverage with respect to . Similarly, relying on the definition of the precession frequencies from Eq. (56), Eq. (52)can be rewritten as (B.3)where in the r.h.s. the source term S_{2}(ℰ_{1},ℰ_{2},τ) obeys (B.4)Notice that Eq. (B.3)is linear in , symmetric in 1 and 2, and can therefore be solved by working out the Green’s function, , associated with the linear differential operator of the l.h.s. of Eq. (B.3). The solution for may therefore be written as (B.5)Injecting Eq. (B.5)into Eq. (B.3), one gets the propagation equation satisfied by . It reads (B.6)where we assumed that the source term S_{2}(t) was effectively turned on only for t ≥ 0, so that S_{2}(t< 0) = 0. In addition, the Green’s function has to satisfy the initial condition . When considering Eq. (B.6), it is worth noting that this propagation equation acts separately on the variables and (and the initial condition of is also separable). We may then solve Eq. (B.6)by factoring the 2body Green’s function as the product of two 1body Green’s function so that (B.7)where the 1body Green’s function satisfies the linearised 1body Vlasov equation, namely (B.8)with the initial condition . Heyvaerts (2010) interestingly notes that, if one were to account for contributions associated with strong collisions, such as in the fifth line of Eq. (18), the property of separability from Eq. (B.7)would not hold anymore. Because of causality, Eq. (B.8)only has to be solved for τ′ ≥ 0. To do so, we rely once again on Bogoliubov’s ansatz, and assume that the system’s 1body DF evolves on a slow secular collisional timescale T_{relax}, while the fluctuations evolve much faster on a secular collisionless timescale T_{sec}. As a consequence, in Eq. (B.8), which describes the evolution of fluctuations, we may assume to be frozen. Therefore, the correlations at a given time τ can be seen as functionals of evaluated at the very same time. To solve Eq. (B.8), we introduce the Laplace transfom following the convention (B.9)In Eq. (B.9), for the inverse Laplace transform, the Bromwich contour ℬ in the complex ωplane should pass above all the poles of the integrand, that is Im [ ω] should be large enough. The Laplace transform of Eq. (B.8)gives (B.10)where the source term on the r.h.s. comes from the initial condition. We now perform a Fourier transform with respect to the slow angles of Eq. (B.10), following the convention from Eq. (62). We multiply Eq. (B.10)by and get (B.11)Equation (B.11)introduced the bare susceptibility coefficients , associated with the Fourier transform of the interaction potential and defined in Eq. (66). Equation (B.11)can easily be rewritten as (B.12)At this stage, it is important to note that Eq. (B.12)takes the form of a Fredholm equation for the Green’s function , as it appears twice in the l.h.s., in particular under the form of an integral term. The trick to solve such an equation is to rely on Kalnaj’s matrix method (Kalnajs 1976), and introduce a basis of potential and densities (ψ^{(p)},ρ^{(p)}) presented in Eq. (59), thanks to which potential perturbations may be decomposed. Let us first decompose the rescaled interaction potential U from Eq. (20)on these basis elements. We consider the function x_{1} → U(  x_{1}−x_{2} ) and project it on the basis ψ^{(p)}(x_{1}). This takes the form U(  x_{1}−x_{2}  ) = ∑ _{p}u_{p}(x_{2}) ψ^{(p)}(x_{1}), where the coefficients u_{p}(x_{2}) are given by (B.13)As they were defined in Eq. (66)as the Fourier transform in angles of the averaged interaction potential , the bare susceptibility coefficients can immediately be written as (B.14)where the averaged Fouriertransformed basis elements were introduced in Eq. (62). In order to invert Eq. (B.12), we perform on the same operations than those acting on . This amounts to multiplying Eq. (B.12)by , so that it becomes (B.15)In order to further simplify Eq. (B.15), let us introduce the notations (B.16)Using the response matrix introduced in Eq. (61), Eq. (B.15)then becomes (B.17)Assuming that the system always remains dynamically stable, so that can be inverted (where I stands for the identity matrix), Eq. (B.17)leads to (B.18)Injecting this relation into Eq. (B.12), can be written as (B.19)where the dressed susceptibility coefficients were introduced in Eq. (60). Relying on the inverse Fourier transform in angles from Eq. (62), we finally obtain the expression of the 1body Green’s function as (B.20)
Appendix B.2: Simplifying the collision operator
Given the explicit calculation of the 1body Green’s function in Eq. (B.20), one may proceed to the evaluation of the collision operator from Eq. (B.2). Relying again on Bogoliubov’s ansatz in Eq. (B.5), we may perform the replacement . Given the factorisation of the Green’s function from Eq. (B.7), and the inverse Laplace transform from Eq. (B.9), the collision operator then takes the form (B.21)where the Laplace transformed 1body Green’s functions were introduced in Eq. (B.20). Let us then rewrite Eq. (B.21)simply as a function of the system’s 1body DF only. Integrating Eq. (B.21)with respect to , , , and , one gets (B.22)with the notations and . Using the explicit expression of the Fourier coefficients of the 1body Green’s function from Eq. (B.20), Eq. (B.22)becomes (B.23)with the shortened notations , as well as and . The rest of this section is devoted to simplifying Eq. (B.23), which still involves a triple integration over action space, two integrals over frequency space and one time integration. In the following, we will first integrate over two actions, then over time, and over the two remaining frequencies, the trickiest step.
Let us first deal with the integration and sum with respect to J_{2} and . It requires to evaluate (B.24)To obtain Eq. (B.24), we relied on the intrinsic definition of the dressed susceptibility coefficients (see Eq. (A.8) in Chavanis 2012) reading (B.25)which is straightforward to obtain given the basis decompositions of the susceptibility coefficients from Eqs. (60)and (B.14), and the definition of the response matrix from Eq. (61). Equation (B.23)then becomes (B.26)Next, the integration and sum with respect to J_{1}′ and are performed. These only act on the two last lines of Eq. (B.26). As previously, one relies on the intrinsic definition of the dressed susceptibility coefficients from Eq. (B.25). Two different contributions have to be dealt with: the first one is associated witht the gradient , and the second one with the gradient . The first contribution takes the form (B.27)Similarly, the second contribution takes the form (B.28)Let us now rewrite Eq. (B.26)while relying on the matrix method, that is by using the basis elements ψ^{(p)}. Within the basis, the bare and dressed susceptibility coefficients take the form of Eq. (B.14)and allow for a rewrite of Eq. (60)as (B.29)where the sums over the greek indices are implied. Following more closely Heyvaerts (2010), we introduced here the matrix , where the response matrix is given by Eq. (61). Finally, let us accordingly define the matrix H(ω) as (B.30)Combining the two contributions from Eqs. (B.27)and (B.28), and after some straightforward algebra, Eq. (B.26)becomes (B.31)Next the integrations with respect to τ′ and ω′ in Eq. (B.31)should be performed, which is technically demanding. This equation formally takes the form (B.32)The integration over τ′ is straightforward provided that ω + ω′ has a negative imaginary part. We therefore introduce p> 0 and perform the substitution ω + ω′ → ω + ω′−ip, so that the integration may be computed as (B.33)As the system is assumed to be linearly stable, the poles of the function ω′ → g(ω,ω′) are all in the lower half complex plane and the Bromwich contour ℬ′ has to pass above all these singularities. The only pole in ω′ which remains is then ω′ = −ω + ip and is located in the upper half plane. We carry the integral over ω′ using the residue theorem by closing the contour ℬ′ in the upper half complex plane – this is possible because the integrands decreases sufficiently fast at infinity like 1/  ω′  ^{2}. One therefore gets (B.34)We may now consider the integration with respect to ω in Eq. (B.31). First, we note that the fourth term of Eq. (B.31)vanishes when integrated upon ω. Indeed, by construction, the Bromwich contour ℬ has to pass above all the singularities of the functions of + ω. This contour may then be closed in the upper half complex plane, and, because it surrounds no singularities, gives a vanishing result for this term. Equation (B.31)may then be rewritten as (B.35)Let us now evaluate the term within brackets in the second term of Eq. (B.35). It reads (B.36)using the notation . When considering the limit p → 0, one should be careful with the fact ω = ω_{2} and ω = ω_{2} + ip are on opposite sides of the prescribed integration contour ℬ. Indeed, when lowering the integration contour to the real axis, the pole ω = ω_{2} remains below the contour, while the one in ω = ω_{2} + ip is above it. In this limit, the term in bracket in Eq. (B.36)becomes [1/(ω−ω_{2} + i0)−1/(ω−ω_{2}−i0)]. We may rely on Plemelj formula (B.37)where stands for Cauchy principal value. Equation (B.36)can then be evaluated and reads When lowering the contour ℬ to the real axis, one can also compute the integration with respect to ω of the first term in Eq. (B.35). Once again, the system being stable, the poles of are all located in the upper half plane, and there remains only one pole on the real axis in ω = ω_{1}. The contour ℬ is closed in the lower half plane and only encloses this second pole. Accounting for the direction of integration, the residue theorem gives a factor −2iπ and Eq. (B.35)then becomes (B.38)keeping track of the small positive imaginary part in the pole 1/(ω_{2}−ω_{1} + i0) associated with the fact that the contour ℬ passed above the pole ω = ω_{1}. Relying on the expression of the susceptibility coefficients from Eq. (B.29), Eq. (B.38)can immediately be rewritten as (B.39)where we made the change for the first term. We note that is real in Eq. (B.29).
Let us now rely on the fact that the collision term is also a real quantity. In Eq. (B.39), because of the prefactor “i”, we may restrict ourselves only to the imaginary part of the terms within brackets. The first term of Eq. (B.39)requires us to study (B.40)In order to compute the term within brackets, we rely on the identity (B.41)The inner term within parenthesis in Eq. (B.41)reads (B.42)where Plemelj formula was used once again. Combining Eqs. (B.40)and (B.42)yields (B.43)This contribution corresponds to the drift term in the BalescuLenard equation.
To evaluate the second term in Eq. (B.39), we make use of the relation , as demonstrated in note [83] of Chavanis (2012), while relying on Plemelj formula. This second term corresponds to the diffusion term in the BalescuLenard equation. All calculations are straightforward. Gathering these two contributions and keeping track of the signs of the various terms, we finally get the expression of the collision term as (B.44)This collision term, together with Eq. (B.1)finally yields the BalescuLenard Eq. (58). It now only involves the divergence of a flux corresponding to a simple integration over action space, and a physically motivated resonant condition and amplification factor, as discussed in the main text.
Appendix B.3: Multicomponent BalescuLenard derivation
Let us explain how one can adapt the formalisms presented in the main text to the situation where the system is composed of multiple components. The different components are indexed by the letters “a” and “b”. We assume that the component “a” is made of N_{a} particles of individual mass μ_{a}, and the total mass of this component is . When accounting for multiple components and placing ourselves within the democratic heliocentric coordinates from Eq. (3), the total Hamiltonian from Eq. (7)becomes (B.45)where stands for the position and velocity of the ith particle of component “a”. The various terms appearing in Eq. (B.45)are respectively the kinetic energy of the particles, the Keplerian potential due to the central BH, the selfgravity among a given component, the interaction between particles of different components, the relativistic potential corrections Φ_{r}, and finally the additional kinetic terms due to the change of coordinates from Eq. (3). One should pay attention to the normalisation of the component Φ_{r}. Indeed, we rewrite this potential as μ_{a}M_{⋆}Φ_{r}, where we introduce the total active mass of the system as . This allows for a rewriting similar to Eq. (7). The dynamics of individual particles is then given by the Hamilton’s equations associated with the Hamiltonian from Eq. (B.45). We now introduce the system’s total PDF , which gives the probability of finding at time t, the particle 1 of the component “a” at position with velocity , etc. As in Eq. (8), we normalise P_{tot} so that (B.46)Following Eq. (9), the dynamics of P_{tot} is governed by Liouville’s equation which becomes (B.47)We define the system’s reduced PDFs (see Eq. (10)) where one integrates P_{tot} over all particles, except n particles belonging respectively to the components a_{1}, ..., a_{n}. Our aim is now to write the two first equations of the associated BBGKY hierarchy. To get the evolution equation of , one proceeds as in Eq. (11), by integrating Eq. (B.47)over all particles except . In order to clarify the upcoming calculations, we will from now on neglect any contributions associated with the last additional kinetic terms from Eq. (B.45). Indeed, we justified in Eq. (47), that, because of the ansatz from Eqs. (41)and (46), once averaged over the fast Keplerian angle, these terms do not contribute the system’s dynamics at the considered order of our kinetic developments. Relying on the symmetry of P_{tot} with respect to interchanges of particles of the same component, one gets (B.48)In Eq. (B.48), we used the same notations as in Eq. (11), and introduced as ℱ_{1a0} the force exerted by the BH on particle 1^{a}, ℱ_{1ar} the force acting on particle 1^{a} associated with the relativistic corrections, and ℱ_{ij} the force between two particles. To obtain the second equation of the hierarchy, one may proceed similarly and integrate Eq. (B.47)for all particles, except 2. Two different cases should be considered, depending on whether one is considering or (with a ≠ b). Let us first consider the diffusion equation satisfied by . Integrating Eq. (B.47)with respect to all particles except and , one gets (B.49)Similarly, starting from Eq. (B.47), and integrating it with respect to and (for b ≠ a), one gets (B.50)As in Eq. (12), we now introduce the renormalised DFs , , and as (B.51)where we assumed that “a”, “b”, and “c” were associated with different components. With these new normalisations, Eq. (B.48)immediately becomes (B.52)where one should note that the sum over “b” runs for all components, which allows for a generic writing. Equations (B.49)and (B.50)can then be cast under the same generic form (B.53)Let us insist on the fact that Eq. (B.53)holds for both the cases where “a” and “b” are equal or different, and the sum on “c” runs for all components. As in Eqs. (14)and (15), one can now define the cluster representation of the DFs which, in this multicomponent context, reads (B.54)and (B.55)Following Eq. (16), we assume that scales like the inverse of the number of particles, while scales like the square of the inverse of the number of particles. Using the decompositions from Eqs. (B.54)and (B.55), and keeping only terms of order 1 /N_{a} or larger (where “a” runs over all the components), the first Eq. (B.52)of the BBGKY hierarchy becomes (B.56)while the second Eq. (B.53)becomes (B.57)Much like for Eq. (19), let us introduce the system’s 1body DF F^{a} and 2body autocorrelation as (B.58)noting the slightly different normalisations of , so as to ensure a symmetric rescaling with respect to “a” and “b”. We also follow Eqs. (20)and (21)to rescale the interaction potential as well as the relativistic corrections with the mass of the BH. Given these various renormalisations, Eq. (B.56)becomes (B.59)where the small parameter ε = M_{⋆}/M_{•} was introduced. Similarly, Eq. (B.57)becomes (B.60)where we introduced the small parameter η_{a} = μ_{a}/M_{⋆} of order 1 /N_{a}. Equations (B.59)and (B.60)are the direct analogs of Eqs. (22)and (23), when one considers a system with multiple components.
As was done in Sect. 3, one may now rewrite the two previous evolution equations within the appropriate angleaction coordinates for the BHinduced Keplerian motion. We perform a degenerate angleaverage as defined in Eq. (33), and assume that F^{a} and satisfy the ansatz from Eqs. (41)and (46). It is then straightforward to rewrite Eq. (B.59)as (B.61)where we used the rescaled time τ = (2π)^{d−k}εt from Eq. (50), with ε = M_{⋆}/M_{•}. Following Eq. (43), we also introduced the total averaged selfconsistent potential as (B.62)where the averaged potential is given by (B.63)In Eq. (B.63), the averaged interaction potential introduced in Eq. (44)was used. One can similarly rewrite Eq. (B.60)as (B.64)With the two coupled evolution Eqs. (B.61)and (B.64)while keeping track of the different mass prefactors, one can follow the path presented in the previous subsection to derive Eq. (69), the appropriate closed kinetic equation for .
Appendix C: From FokkerPlanck to Langevin
Following Risken & Frank (1996), let us briefly recall how one may obtain the stochastic Langevin equation describing the individual dynamics of a test particle starting from a FokkerPlanck equation describing the diffusion of the system’s DF as a whole. We start from the generic writing of the degenerate BalescuLenard equation from Eq. (63), written as an anisotropic selfconsistent FokkerPlanck equation. It reads (C.1)where, following the notations from Eq. (64), the drift vector A(J,τ) and diffusion tensor D(J,τ) are introduced as (C.2)One should keep in mind that in Eq. (63)the drift and diffusion coefficients also secularly depend on the system’s DF , but this was not written out to simplify the notations. Following the notations of Eq. (4.94a) in Risken & Frank (1996), we may immediately rewrite Eq. (C.1)as (C.3)where the first and secondorder diffusion coefficients read (C.4)One should pay attention to the fact that the diffusion of stars takes place in the full action domain J, while only gradients with respect to the slow actions J^{s} are present in Eq. (C.3). Of course, by enlarging the diffusion coefficients D^{(1)} and D^{(2)} with zero coefficients for all the adiabatically conserved fast actions J^{f}, it is straightforward to rewrite Eq. (C.3)as a diffusion equation in Jspace involving derivatives with respect to J.
Let us now focus on the individual dynamics of a given test particle. We denote as its position in action space at time τ. This test particle then undergoes a stochastic diffusion consistent with the averaged diffusion captured by the FokkerPlanck Eq. (C.3), namely a Langevin equation reading (C.5)where we introduced the Langevin vector and tensor h and g, as well as the stochastic Langevin forces Γ(τ), whose statistics satisfy (C.6)with I the identity matrix. Following Eq. (3.124) of Risken & Frank (1996), we may now express the Langevin coefficients from Eq. (C.5)as a function of the coefficients appearing in the FokkerPlanck Eq. (C.3). The secondorder diffusion tensor D^{(2)} being definite positive, we introduce as one of its square root, so that one has the component relations (C.7)Equation (C.7)therefore allows us to fully specify the properties of the diffusion of an individual particle as captured by the Langevin Eq. (C.5). Of course, selfconsistency requires that the diffusion coefficients D^{(1)} and D^{(2)}, and therefore the Langevin coefficients h and g should be updated as the system’s DF changes on secular timescales, as mentioned in the main text.
All Figures
Fig. 1 Illustration of the resonance condition appearing in the degenerate inhomogeneous BalescuLenard Eq. (58). Topleft: a set of two resonant orbits precessing at the same frequency ω_{s}. Topright: in the rotating frame at frequency ω_{s} in which the two orbits are in resonance. Bottom: fluctuations of the system’s DF in action space caused by finiteN effects and showing overdensities for the blue and red orbits. The dashed line correspond to the critical resonant line in action space along which the resonance condition Ω^{s}(J) = ω_{s} is satisfied. The two set of orbits satisfy a resonance condition for their precession frequencies, and uncorrelated sequences of such interactions lead to a secular diffusion of the system’s orbital structure following Eq. (58). Such resonances are non local in the sense that the resonant orbits need not be close in action space nor in position space. As emphasised in Sect. 6.1 for axisymmetric razorthin discs, symmetry enforces , so that the two orbits are caught in the same resonance. 

Open with DEXTER  
In the text 
Fig. 2 Illustration of the typical dependence of the precession frequencies and (Eqs. (111)and (112)) as a function of the distance to the central BH. The relativistic precession frequencies, diverge as the star gets closer to the central BH, while the selfconsistent ones are typically the largest for stars in the neighbourhood of the considered disc. The black dots give all the locations, whose precession frequency is equal to ω_{s} (illustrated by the dotted horizontal line). These positions are in resonance and will therefore have a nonvanishing contribution in the BalescuLenard Eq. (92). Because Eq. (92)involves the product of the system’s DF in the two resonating locations, the resonant coupling between the two outer points (which belong to the region where the disc dominates) will be much stronger than the couplings involving the inner point (which does not belong to core of the disc). As stars move inward, their precession frequencies increase up to a point where this prevents any resonant coupling with the disc’s region. This effectively stops the secular diffusion, and induces a diffusion barrier. 

Open with DEXTER  
In the text 
Fig. 3 Illustration of the individual dynamics of stars in the (j,a) = (L/I,I^{2}/GM_{•}) space, as given by the Langevin Eq. (114). The grey region corresponds to the capture region, within which stars inevitably sink into the BH. As I is an invariant of the diffusion (see Eq. (114)), stars’ diffusion is onedimensional, conserves a, and occurs only in the jdirection. The background dotted curves illustrate the contour lines of the precession frequency given by the function (j,a) → Ω^{s}(j,a). As illustrated in Fig. 2, the precession frequencies get larger as particles approach the central BH, due to the contributions from the relativistic precession frequencies. The blue and red orbits precess at the same frequency ω_{s}, so that they belong to the same critical resonant line γ_{ωs}, which allows them to resonate one with another. As the precession frequencies diverge in the vicinity of the BH, such resonant couplings are significantly less likely as stars get closer to the BH, which effectively creates a diffusion barrier in action space, the socalled Schwarzschild barrier. 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.