Free Access
Volume 598, February 2017
Article Number A71
Number of page(s) 24
Section Astrophysical processes
Published online 02 February 2017

© ESO, 2017

1. Introduction

The stars in a stellar cluster surrounding a dominant supermassive black hole (BH) move in a quasi-Keplerian 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 quasi-Keplerian orbit. The consequences of this idea were first developed by Rauch & Tremaine (1996), who showed that wire-wire 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 two-body encounters. They named this phenomenon “resonant relaxation”, because such enhanced relaxation occurs more generally in any potential in which the three-dimensional vector of stellar orbital frequencies Ω satisfies a commensurability condition of the form n·Ω ≃ 0 for some vector of integers n = (n1,n2,n3).

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 star-BH 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 N-body 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 N-wires codes (e.g. Kocsis & Tremaine 2015) in which individual stars are replaced by orbit-averaged 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 one-body distribution function and its two-body correlation function are reduced to a single equation – the Balescu-Lenard equation – that describes the evolution of the one-body 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 long-range interacting systems. The original Balescu-Lenard formalism was developed for homogeneous plasmas. One way of generalising it to inhomogeneous self-gravitating 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 Balescu-Lenard formalism to inhomogeneous systems has been presented by Heyvaerts (2010) and Chavanis (2012), who reformulate the non-linear kinetic equation in terms of the angle-action variables that are appropriate for spatially inhomogeneous multi-periodic systems. Fouvry et al. (2015a,b, 2016b) have applied this formalism to describe the secular response of tepid self-gravitating stellar discs. The resulting inhomogeneous Balescu-Lenard equation accounts for the self-induced secular orbital diffusion of a self-gravitating system driven by the internal shot noise due to the finite number N of particles involved. In common with all results based on the Balescu-Lenard 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 Ntd, where td is the dynamical time. The secular interactions between particles need not be local in space: they need only correspond to gravitationally amplified long-range correlations via resonances.

In its original form, however, the Balescu-Lenard formalism assumes that resonances are localised in action space and not degenerate. Therefore, it must be re-examined 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 Balescu-Lenard 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 self-gravity 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 long-term evolution of a disc or a sphere of (possibly relativistic) stars near a galactic centre, or a proto-planetary debris disc circling a star. As such, it captures the secular effects of a sequence of polarised wire-wire 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 angle-action coordinates for such quasi-Keplerian systems. Section 4 averages the corresponding dynamical equations over the fast angles of the Keplerian motion. Section 5 presents the degenerate one and multi-component Keplerian Balescu-Lenard equations (Appendix B details both derivations, following the steps of Heyvaerts 2010, while also correcting for a minor issue in the multi-component case). Section 6 discusses applications to the cases of razor-thin 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 Balescu-Lenard 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 Xn 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 Lense-Thirring 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,x1,...,xN) defined as (3)where we have introduced the total mass of the system Mtot = M + M. In Eq. (3), x corresponds to the position of the system’s centre of mass and xi 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,p1,...,pN) 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 two-body 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 N-body probability distribution function (PDF) PN1,...N,t) defined so that PN1,...N,t) dΓ1...N is at time t the probability of finding particle 1 within the volume element 1 located at the phase space point Γ1 = (x1,v1), particle 2 within 2 of the phase space point Γ2 = (x2,v2), and so on. We normalise PN such that (8)It evolves according to Liouville’s equation (9)The dynamics of the individual particles are given by Hamilton’s equations, μdxi/ dt = ∂H/vi and μdvi/ dt = −∂H/xi, where the system’s Hamiltonian was obtained in Eq. (7). From PN, 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 Pn, we integrate Liouville’s Eq. (9)over n + 1...N and use the fact that PN 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 = −μ∂Uij/xi, using the shorthand notation Uij = U(|xixj|). The force exerted by the BH on particle i is denoted by Mi0 = −MUi0/xi and the force associated with the relativistic corrections as Mir = −MΦrel/xi.

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 fn that arise from correlations among particles, let us introduce the cluster representation of the DFs. We define the 2-body correlation function g2 in terms of f1 and f2 via (14)Similarly, the 3-body correlation function g3 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 | f1 | ~ 1, | g2 | ~ 1 /N, and | g3 | ~ 1 /N2. 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 non-Keplerian 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)(v1 + v2g2/x1 which may be neglected), while all the terms from the two last lines are of order 1 /N2 and may therefore be neglected. Notice that the first term on the fifth line of Eq. (18), which, while being of order 1 /N2, 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 1-body DF F and its 2-body 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 angle-action 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 angle-action coordinates

In Eqs. (22)and (23), one can note the presence of an advection term v1·/x1 + 10·/v1 associated with the Keplerian motion driven by the central black hole. The next step of the derivation is to introduce the appropriate angle-action 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 angle-action 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 HKep 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 Jr is the radial action, L the magnitude of the angular momentum, and Lz its projection along the z-axis. The Keplerian Hamiltonian then becomes HKep = HKep(Jr + L). Another choice of angle-action coordinates in 3D is given by the Delaunay variables (Sridhar & Touma 1999; Binney & Tremaine 2008) defined as (26)In Eq. (26), (I = Jr + L,L,Lz) 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 HKep = HKep(I), so that the angle w advances at the frequency = ΩKep = HKep/∂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 long-term 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 razor-thin disc. In this space, we consider an integrable potential ψ and an associated angle-action mapping (x,v) → (θ,J). A potential is said to be degenerate if there exists n ∈ Zd 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 angle-action coordinates from Eq. (25), the frequencies and degeneracy vectors are given by (28)so that k = 2. Using the Delaunay angle-action 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 angle-action 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 angle-action coordinates to get simpler degeneracies. Indeed, let us assume that in our initial angle-action coordinates (θ,J), we have at our disposal k degeneracy vectors n1, ... , nk. 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 ei are the natural basis elements of Zd. Following Morbidelli (2002), as the vectors ni are assumed to be linearly independent, we may complete this family with dk vectors nk + 1,...,nd ∈ Zd to have a basis over Qd. We then define the transformation matrix of determinant 1 as (30)and the new angle-action coordinates (θ′,J′) are defined as (31)One can check that (θ′,J′) are indeed new angle-action 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 ≤ ik. The degeneracies of the potential got simpler. In the upcoming calculations, we will always consider such simpler angle-action coordinates, and we introduce the notations (32)where θs and Js respectively stand for the slow angles and actions, while θf and Jf 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 non-zero for the fast angles. Let us finally define the degenerate angle-average with respect to the fast angles as (33)We now use these various properties to rewrite Eqs. (22)and (23)using the angle-action 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 G1(x,v), and G2(x,v), one has (36)In order to shorten the notations, let us now introduce the rescaled self-consistent potential Φ as (37)One can now rewrite Eq. (22)within these angle-action 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 TKep = 1/ΩKep associated with the Keplerian advection term ; (ii) the secular collisionless timescale of evolution Tsec = ε-1TKep associated with the potential contributions ε [ Φ + Φr]; and finally (iii) the collisional timescale of relaxation Trelax = NTsec, 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 phased-mixed 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 self-consistent 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 self-consistent potential was introduced as (43)In Eq. (43), for clarity, the notation was shortened for the self-consistent potential as . The (doubly) averaged interaction potential is defined as (44)while the angle-averaged potential was also introduced as (45)where the prefactor 1/(2π)dk, was introduced for convenience. As emphasised in Eq. (42), one should note that at first order in ε and zeroth order in ϵ, the self-consistent 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 2-body 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 Js, 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 angle-averaging 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 2-body 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 quasi-stationary equilibrium, so that. One then expects that this out-of-equilibrium system will undergo a phase of violent relaxation (Lynden-Bell 1967), allowing it to rapidly reach a quasi-stationary 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 quasi-stationary 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 quasi-stationary stable state, one can study the secular evolution of this system along quasi-stationary equilibria. Such a long-term 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 long-term 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 quasi-stationary equilibria, the dynamics of an isolated system can also be driven by finite-N 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) Balescu-Lenard 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 finite-N fluctuations capture the known mechanism of resonant relaxation (Rauch & Tremaine 1996). See Bar-Or & Alexander (2014) for a similar study of the effect of the finite-N stochastic internal forcing via the so-called η-formalism.

We note that one could also consider the secular evolution of a non-axisymmetric 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 Balescu-Lenard equation for such a configuration would first involve identifying new angle-action variables for the non-axisymmetric 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 self-induced orbital precession is significant, but the self-gravity 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 Balescu-Lenard equation

We now show how to obtain the closed kinetic equations – the degenerate Balescu-Lenard 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 Tsec), through an out-of-equilibrium 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 razor-thin 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 self-consistent 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 quasi-identical to the traditional coupled BBGKY equations considered in Heyvaerts (2010) to derive the inhomogeneous Balescu-Lenard equation for non-degenerate 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 quasi-linear 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 1-body DF . Injecting this expression in Eq. (51)yields finally a closed kinetic equation quadratic in . The detailed calculations required to derive the inhomogeneous degenerate Balescu-Lenard equation are presented in Appendix B.

5.1. The one component Balescu-Lenard equation

In its explicitly conservative form, the degenerate inhomogeneous Balescu-Lenard 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 zeroth-order potential. The r.h.s. of Eq. (58)is the degenerate inhomogeneous Balescu-Lenard collision operator, which describes the secular diffusion induced by dressed finite-N fluctuations. It describes the distortion of Keplerian orbits as their actions diffuse through their self-interaction. 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 J2 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.

thumbnail Fig. 1

Illustration of the resonance condition appearing in the degenerate inhomogeneous Balescu-Lenard Eq. (58). Top-left: a set of two resonant orbits precessing at the same frequency ωs. Top-right: 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 finite-N 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 razor-thin 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 so-called 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 self-gravity into account. In order to solve Poisson’s non-local 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 so-called Bogoliubov’s ansatz, as shown in Appendix B.1.

One can straightforwardly rewrite the Balescu-Lenard Eq. (58)as an anisotropic non-linear 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 Balescu-Lenard Eq. (58)becomes the degenerate Landau equation (see Polyachenko & Shukhman 1982; Chavanis 2013b, for the non-degenerate case), which reads (65)Notice that this is just the previous Balescu-Lenard Eq. (58)with the dressed replaced by the bare susceptibility coefficients . The latter are related to the (partial) Fourier transform of the interaction potential (Lynden-Bell 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 = (Js,Jf), do not allow for changes in the fast actions Jf. 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 Jf = const.

5.2. Multiple components black hole environment

It is of prime importance to follow the joint long-term 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 Balescu-Lenard 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 Fa. 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 multi-component 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 /Na, and the single mass Balescu-Lenard 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 J1, 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 multi-component Balescu-Lenard formalism captures the secular effect of multiple resonant (non-local) 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(J1), 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 J1J2, 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 multi-component Balescu-Lenard Eq. (69). Indeed, defining the system’s total entropy Stot, summed for all components, as (81)one can again show that for s′′(x) = 1 /x, one has dStot/ 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 follow-up papers.

6.1. Razor-thin axisymmetric discs

Let us first specialise the degenerate Balescu-Lenard Eq. (58)to razor-thin 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 angle-action variables from Eq. (26)become (82)Symmetries of the interaction potential lead to relationships among the susceptibility coefficients, which simplify the Balescu-Lenard equation. The rescaled interaction potential U12 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 angle-action ones can be written as (84)where the semi-major 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 razor-thin 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 razor-thin system, one can assume the basis elements from Eq. (59)to be generically of the form (88)where p and np 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 Balescu-Lenard 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 self-gravity, 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 razor-thin 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 Balescu-Lenard 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 d-dimensional 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 J1, and introducing ω = Ωs(J1), 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 J1. 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(J1,J2) as (99)as well as the resonant contribution | ∇(Ωs(J2)) | given by (100)We note that Eq. (98)is now a simple one-dimensional integral involving a regular integrand.

In summary, because the quasi-stationary 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 razor-thin discs in the Landau limit is postponed to a follow-up paper, as Eq. (43)involves a singular triple integral over wire-wire interactions. Similarly, the study of the long-term evolution of quasi-stationary non-axisymmetric razor-thin 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 Balescu-Lenard Eq. (58)to spherically symmetric systems. The general procedure is the same as for the razor-thin 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 two-dimensional 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 razor-thin 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 angle-action 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) = Lz/L. Therefore the interaction potential of Eq. (102)and its angle-averaged 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, mp and np 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 Balescu-Lenard 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 Lz, 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 Balescu-Lenard 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 razor-thin 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 self-consistent Newtonian potential, (111)The second is the additional contribution from relativistic effects. We derive it in Appendix A. In the case of a razor-thin 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 self-consistent 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.

thumbnail 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 self-consistent 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 non-vanishing contribution in the Balescu-Lenard 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 so-called Schwarzschild barrier.

This explanation of the Schwarzschild barrier using the notion of resonant coupling is directly related to the explanation proposed in Bar-Or & 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 Balescu-Lenard resonance condition.

6.4. Solving the Balescu-Lenard equation by Monte Carlo sampling

This suppression of diffusion in the neigbourhood of the BH can also be illustrated by considering the orbit-averaged motion of individual wires. Equation (92)takes the form of a diffusion equation in action space, where one follows self-consistently 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 razor-thin 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 L-direction are given by (115)and the Langevin stochastic force Γ(τ) satisfies Eq. (C.6). As in Eq. (68), the individual fast action Jf = I is preserved during the wire’s evolution.

Equation (114)describes the diffusion of an individual test wire when embedded in the self-induced 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 0s 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 post-Newtonian 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.

thumbnail Fig. 3

Illustration of the individual dynamics of stars in the (j,a) = (L/I,I2/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 one-dimensional, conserves a, and occurs only in the j-direction. 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 so-called Schwarzschild barrier.

Open with DEXTER

Following the representations from Bar-Or & Alexander (2016), Fig. 3 represents the diffusion of stars in the (j,a) = (L/I,I2/ (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 j-direction 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 razor-thin axisymmetric disc.

6.5. Evolution of BH mass and spin

The calculations above successfully explain the existence of the so-called 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 two-body relaxation (not accounted for in the orbit-averaged 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 Bar-Or & 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 ms takes the form of a preferential diffusion in the direction of ms. This diffusion is therefore anisotropic because it is maximum for nms and equal to 0 for n·ms = 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 self-consistent 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 self-gravitating, 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 loss-cone, 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 hyper-velocity 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 self-gravitating systems of N particles (Heyvaerts 2010; Chavanis 2012) to quasi-Keplerian 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 self-consistent dressed equations (Eq. (58)and its multi-component and stochastic counterparts, Eqs. (69)and (C.5)respectively) account for the dynamical degeneracies in quasi-Keplerian systems. Because purely Keplerian orbits do not precess, the dynamical evolution of such degenerate systems may differ significantly from that of fully self-gravitating systems, such as discs and spheroids. In particular, to a good approximation stars behave as if they were smeared out into orbit-averaged 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 quasi-Keplerian Balescu-Lenard Eq. (58)is quadratic in the phase-averaged distribution function and describes i) the self-gravity 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 post-Newtonian 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 quasi-linear self-consistent 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 Balescu-Lenard 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 non-resonant two-body relaxation is not accounted for in the Balescu-Lenard equations we derive here. This is usually appropriate though, because the derivation of these equations ignores terms of order O(1 /N2), which means that they are valid only on timescales Ntd, where td is the dynamical timescale. Such timescales are typically expected to be much shorter than the non-resonant two-body relaxation time. When investigating specifically the vicinity of supermassive black holes, we found that the quasi-Keplerian Balescu-Lenard 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 multi-component formulation, the Balescu-Lenard 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 quasi-Keplerian 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 razor-thin axisymmetric discs (Sridhar & Touma 2017) corresponds to the Landau limit in which one uses the bare susceptibility coefficients from Eq. (94)in the Balescu-Lenard 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; Bar-Or & 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 Bar-Or & Alexander (2014) and implemented in detail in Bar-Or & 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 2-body relaxation for the loss-cone problem. They showed that on longer timescales, 2-body non-resonant 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 loss-lines, so that its overall effect on plunge rates is small.

The Balescu-Lenard 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 Balescu-Lenard equation requires no such external input, because the system’s discreteness is described self-consistently. Second, in the η-formalism the self-gravity of the response to the noise is difficult to account for. In the Balescu-Lenard 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 razor-thin stellar discs). We note that, just as in Monte Carlo schemes, the Balescu-Lenard 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 N-body 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 self-consistent potential), but interactions among field stars are ignored. The test stars are then followed by direct integration of their motion in the time-varying 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 back-influence of the test stars on the field ones.

7.2. Future work

The Langevin formulation of the Balescu-Lenard equation (Sect. 6.4 and Appendix C) combines the flexibility of Monte Carlo methods with a self-consistent treatment of the dynamics. A subsequent improvement is offered by the possibility of adding the secondary effects of two-body relaxation and gravitational-wave 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 razor-thin 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 wire-wire 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 self-gravity, 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 self-consistency 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 self-gravity. 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).


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 Pierre-Henri 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 ANR-13-BS05-0005 of the French Agence Nationale de la Recherche (


  1. Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, Phys. Rev. Lett., 116, 061102 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  2. Balescu, R. 1960, Phys. Fluids, 3, 52 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  3. Bar-Or, B., & Alexander, T. 2014, Class. Quant. Grav., 31, 244003 [NASA ADS] [CrossRef] [Google Scholar]
  4. Bar-Or, B., & Alexander, T. 2016, ApJ, 820, 129 [NASA ADS] [CrossRef] [Google Scholar]
  5. Binney, J., & Tremaine, S. 2008, in Galactic Dynamics, 2nd edn. (Princeton University Press) [Google Scholar]
  6. Blanchet, L. 2014, Liv. Rev. Relat., 17, 2 [NASA ADS] [CrossRef] [Google Scholar]
  7. Born, M. 1960, in The Mechanics of the Atom (F. Ungar Pub. Co.) [Google Scholar]
  8. Chavanis, P.-H. 2010, J. Stat. Mech., 5, 19 [Google Scholar]
  9. Chavanis, P.-H. 2012, Physica A, 391, 3680 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  10. Chavanis, P.-H. 2013a, Eur. Phys. J. Plus, 128, 126 [CrossRef] [EDP Sciences] [Google Scholar]
  11. Chavanis, P.-H. 2013b, A&A, 556, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Duncan, M. J., Levison, H. F., & Lee, M. H. 1998, AJ, 116, 2067 [NASA ADS] [CrossRef] [Google Scholar]
  13. Fouvry, J.-B., Pichon, C., & Chavanis, P.-H. 2015a, A&A, 581, A139 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Fouvry, J.-B., Pichon, C., Magorrian, J., & Chavanis, P.-H. 2015b, A&A, 584, A129 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Fouvry, J.-B., Pichon, C., & Prunet, S. 2015c, MNRAS, 449, 1967 [NASA ADS] [CrossRef] [Google Scholar]
  16. Fouvry, J.-B., Chavanis, P.-H., & Pichon, C. 2016a, Physica A, 459, 117 [NASA ADS] [CrossRef] [Google Scholar]
  17. Fouvry, J.-B., Pichon, C., Chavanis, P.-H., & Monk, L. 2016b, MNRAS, submitted [Google Scholar]
  18. Frank, J., & Rees, M. J. 1976, MNRAS, 176, 633 [NASA ADS] [CrossRef] [Google Scholar]
  19. Genzel, R., Pichon, C., Eckart, A., Gerhard, O. E., & Ott, T. 2000, MNRAS, 317, 348 [NASA ADS] [CrossRef] [Google Scholar]
  20. Gilbert, I. H. 1968, ApJ, 152, 1043 [NASA ADS] [CrossRef] [Google Scholar]
  21. Goldstein, H. 1950, in Classical mechanics (Addison-Wesley) [Google Scholar]
  22. Hamers, A. S., Portegies Zwart, S. F., & Merritt, D. 2014, MNRAS, 443, 355 [NASA ADS] [CrossRef] [Google Scholar]
  23. Heyvaerts, J. 2010, MNRAS, 407, 355 [NASA ADS] [CrossRef] [Google Scholar]
  24. Hills, J. G. 1988, Nature, 331, 687 [NASA ADS] [CrossRef] [Google Scholar]
  25. Hopman, C., & Alexander, T. 2006, ApJ, 645, 1152 [NASA ADS] [CrossRef] [Google Scholar]
  26. Hörmander, L. 2003, The analysis of linear partial differential operators. I (Springer-Verlag) [Google Scholar]
  27. Jalali, M. A., & Tremaine, S. 2012, MNRAS, 421, 2368 [NASA ADS] [CrossRef] [Google Scholar]
  28. Jocou, L., Perraut, K., Moulin, T., et al. 2014, in Optical and Infrared Interferometry IV, 9146, 91461 [Google Scholar]
  29. Julian, W. H., & Toomre, A. 1966, ApJ, 146, 810 [NASA ADS] [CrossRef] [Google Scholar]
  30. Kalnajs, A. J. 1976, ApJ, 205, 745 [NASA ADS] [CrossRef] [Google Scholar]
  31. Klimontovich, I. 1967, in The statistical theory of non-equilibrium processes in a plasma (M.I.T. Press) [Google Scholar]
  32. Kocsis, B., & Tremaine, S. 2011, MNRAS, 412, 187 [NASA ADS] [CrossRef] [Google Scholar]
  33. Kocsis, B., & Tremaine, S. 2015, MNRAS, 448, 3265 [NASA ADS] [CrossRef] [Google Scholar]
  34. Lenard, A. 1960, Ann. Phys., 10, 390 [NASA ADS] [CrossRef] [Google Scholar]
  35. Lynden-Bell, D. 1967, MNRAS, 136, 101 [NASA ADS] [CrossRef] [Google Scholar]
  36. Lynden-Bell, D. 1994, in Lectures on stellar dynamics (Berlin: Springer Verlag) [Google Scholar]
  37. Madigan, A.-M., Hopman, C., & Levin, Y. 2011, ApJ, 738, 99 [NASA ADS] [CrossRef] [Google Scholar]
  38. Merritt, D. 2015, in Astrophysical Black Holes (Springer) [Google Scholar]
  39. Merritt, D., Alexander, T., Mikkola, S., & Will, C. M. 2011, Phys. Rev. D, 84, 044024 [NASA ADS] [CrossRef] [Google Scholar]
  40. Morbidelli, A. 2002, in Modern Celestial Mechanics (Taylor & Francis) [Google Scholar]
  41. Pichon, C. 1994, Ph.D. Thesis, University of Cambridge [Google Scholar]
  42. Polyachenko, V. L., & Shukhman, I. G. 1982, Soviet Ast., 26, 140 [NASA ADS] [Google Scholar]
  43. Polyachenko, E. V., Polyachenko, V. L., & Shukhman, I. G. 2007, MNRAS, 379, 573 [NASA ADS] [CrossRef] [Google Scholar]
  44. Rauch, K. P., & Ingalls, B. 1998, MNRAS, 299, 1231 [NASA ADS] [CrossRef] [Google Scholar]
  45. Rauch, K. P., & Tremaine, S. 1996, New Astron., 1, 149 [NASA ADS] [CrossRef] [Google Scholar]
  46. Risken, H. 1996, in The Fokker-Planck Equation (Berlin, Heidelberg: Springer) [Google Scholar]
  47. Sridhar, S., & Touma, J. 1999, MNRAS, 303, 483 [NASA ADS] [CrossRef] [Google Scholar]
  48. Sridhar, S., & Touma, J. R. 2016a, MNRAS, 458, 4129 [NASA ADS] [CrossRef] [Google Scholar]
  49. Sridhar, S., & Touma, J. R. 2016b, MNRAS, 458, 4143 [NASA ADS] [CrossRef] [Google Scholar]
  50. Sridhar, S., & Touma, J. R. 2017, MNRAS, 465, 1856 [NASA ADS] [CrossRef] [Google Scholar]
  51. Toomre, A. 1981, in Structure and Evolution of Normal Galaxies, eds. S. M. Fall, & D. Lynden-Bell, 111 [Google Scholar]
  52. Touma, J. R., & Sridhar, S. 2012, MNRAS, 423, 2083 [NASA ADS] [CrossRef] [Google Scholar]
  53. Touma, J. R., Tremaine, S., & Kazandjian, M. V. 2009, MNRAS, 394, 1085 [NASA ADS] [CrossRef] [Google Scholar]
  54. Tremaine, S. 1995, AJ, 110, 628 [NASA ADS] [CrossRef] [Google Scholar]
  55. Tremaine, S. 1998, AJ, 116, 2015 [NASA ADS] [CrossRef] [Google Scholar]
  56. Tremaine, S. 2005, ApJ, 625, 143 [NASA ADS] [CrossRef] [Google Scholar]
  57. Volonteri, M., Dubois, Y., Pichon, C., & Devriendt, J. 2016, MNRAS, 460, 2979 [NASA ADS] [CrossRef] [Google Scholar]
  58. 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 z-direction 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 TKep = 2π/ ΩKep = 2πI3/ (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 semi-major axis a and the eccentricity e satisfy a = I2/ (GM) and 1−e2 = (L/I)2. We also introduced the Hamiltonian as (A.3)Similarly, Eq. (5.118) of Merritt (2015) gives that the 1.5PN Lense-Thirring 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 z-direction. This is immediately associated with a precession frequency given by (A.5)relying on the relation Lz = 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 Lz. 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 Balescu-Lenard equation presented in Heyvaerts (2010) in this new quasi-Keplerian regime. The starting point of this derivation is the two coupled averaged Eqs. (51)and (52), which involve the system’s averaged 1-body DF , and its averaged 2-body 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 multi-component 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 angle-average 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 S2(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 S2(t) was effectively turned on only for t ≥ 0, so that S2(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 2-body Green’s function as the product of two 1-body Green’s function so that (B.7)where the 1-body Green’s function satisfies the linearised 1-body 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 1-body DF evolves on a slow secular collisional timescale Trelax, while the fluctuations evolve much faster on a secular collisionless timescale Tsec. 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 x1U( | x1x2 |) and project it on the basis ψ(p)(x1). This takes the form U( | x1x2 | ) = ∑ pup(x2) ψ(p)(x1), where the coefficients up(x2) 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 Fourier-transformed 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 1-body Green’s function as (B.20)

Appendix B.2: Simplifying the collision operator

Given the explicit calculation of the 1-body 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 1-body Green’s functions were introduced in Eq. (B.20). Let us then rewrite Eq. (B.21)simply as a function of the system’s 1-body 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 1-body 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 J2 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 J1 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 Balescu-Lenard 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 Balescu-Lenard 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 Balescu-Lenard 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: Multi-component Balescu-Lenard 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 Na 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 self-gravity 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 μaMΦ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 Ptot so that (B.46)Following Eq. (9), the dynamics of Ptot is governed by Liouville’s equation which becomes (B.47)We define the system’s reduced PDFs (see Eq. (10)) where one integrates Ptot over all particles, except n particles belonging respectively to the components a1, ..., an. 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 Ptot 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 1a, 1ar the force acting on particle 1a 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 multi-component 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 /Na 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 1-body DF Fa and 2-body 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 /Na. 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 angle-action coordinates for the BH-induced Keplerian motion. We perform a degenerate angle-average as defined in Eq. (33), and assume that Fa 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π)dkεt from Eq. (50), with ε = M/M. Following Eq. (43), we also introduced the total averaged self-consistent 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 Fokker-Planck 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 Fokker-Planck equation describing the diffusion of the system’s DF as a whole. We start from the generic writing of the degenerate Balescu-Lenard equation from Eq. (63), written as an anisotropic self-consistent Fokker-Planck 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 second-order 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 Js 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 Jf, it is straightforward to rewrite Eq. (C.3)as a diffusion equation in J-space 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 Fokker-Planck 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 Fokker-Planck Eq. (C.3). The second-order 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, self-consistency 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

thumbnail Fig. 1

Illustration of the resonance condition appearing in the degenerate inhomogeneous Balescu-Lenard Eq. (58). Top-left: a set of two resonant orbits precessing at the same frequency ωs. Top-right: 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 finite-N 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 razor-thin discs, symmetry enforces , so that the two orbits are caught in the same resonance.

Open with DEXTER
In the text
thumbnail 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 self-consistent 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 non-vanishing contribution in the Balescu-Lenard 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
thumbnail Fig. 3

Illustration of the individual dynamics of stars in the (j,a) = (L/I,I2/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 one-dimensional, conserves a, and occurs only in the j-direction. 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 so-called Schwarzschild barrier.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.