J Heyvaerts

A Balescu-Lenard type kinetic equation for the collisional evolution of stable self-gravitating systems

Keywords:

A kinetic equation for the collisional evolution of stable, bound, self gravitating and slowly relaxing systems is established, which is valid when the number of constituents is very large. It accounts for the detailed dynamics and self consistent dressing by collective gravitational interaction of the colliding particles, for the system's inhomogeneity and for different constituent's masses. It describes the coupled evolution of collisionally interacting populations, such as stars in a thick disk and the molecular clouds off which they scatter. The kinetic equation derives from the BBGKY hierarchy in the limit of weak, but non-vanishing, binary correlations, an approximation which is well justified for large stellar systems. The evolution of the one-body distribution function is described in action angle space. The collective response is calculated using a biorthogonal basis of pairs of density-potential functions. The collision operators are expressed in terms of the collective response function allowed by the existing distribution functions at any given time and involve particles in resonant motion. These equations are shown to satisfy an H-theorem. Because of the inhomogeneous character of the system, the relaxation causes the potential as well as the orbits of the particles to secularly evolve. The changing orbits also cause the angle Fourier coefficients of the basis potentials to change with time. We derive the set of equations which describes this coupled evolution of distribution functions, potential and basis Fourier coefficients for spherically symmetric systems. In the homogeneous limit, which sacrifices the description of the evolution of the spatial structure of the system but retains the effect of collective gravitational dressing, the kinetic equation reduces to a form similar to the Balescu-Lenard equation of plasma physics.

INTRODUCTION AND MOTIVATION

The description of collisional relaxation in a self-gravitating system usually rests on a Fokker-Planck equation in which the diffusion and braking coefficients are calculated in the local approximation, taking the finite dimension of the system into account by limiting the impact parameter of the collisions to a length of order of the system's size [START_REF] Chandrasekhar | Principles of stellar dynamics[END_REF][START_REF] Chandrasekhar | [END_REF][START_REF] Binney | Galactic Dynamics[END_REF][START_REF] Spitzer | Dynamical evolution of globular clusters[END_REF]. Although characteristic relaxation times may be somewhat overestimated by this approximation due to the neglect of collective self-gravitational effects (Weinberg 1993), such a kinetic equation may provide in practice a reasonable description of the collisional relaxation of gravitationally bound systems. It nevertheless rests on assumptions which, from a principle point of view, are unsatisfactory because the motion of particles during the collision is regarded as rectilinear and uniform and the system's inhomogeneity, which is basically the reason why collisions with an infinite impact parameter do not occur, is treated by way of an ill-defined cutoff. Moreover, the collective response of the system is not taken into account, since the Fokker-Planck collision term only considers binary collisions between naked particles. A self-gravitating medium, unlike an electrical plasma, does not respond to the presence in it of a particle by screening its interaction potential with other particles. As a result, even distant particles effectively interact, while in electrical, globally neutral, plasmas, the effective interaction distance is limited to the Debye length. In a self-gravitating system, the distance between interacting particles is only limited by the system's inhomogeneity. The spatial structure of the system matters as well as the details of the particle orbits.

The consistent inclusion of collective screening effects in a kinetic equation for electrically interacting weakly coupled particles has been one of the major theoretical achievements in plasma physics when [START_REF] Peyrissac | Equations cinétiques des fluides et des plasmas[END_REF] and Lenard (1960) could derive an equation surpassing in consistency the simple Fokker-Planck equation [START_REF] Spitzer | Dynamical evolution of globular clusters[END_REF]. It is the aim of this paper to derive a similar equation for self-gravitating systems. The task is slightly more difficult because the screening of the electrical interaction at the, usually small, Debye length allows, in electrically interacting systems, to take the homogeneous and uniform motion limits. These limits cannot be taken in a self-gravitating system. We overcome this difficulty by expressing the kinetic equation in action angle space rather than in position momentum space. This is possible when the Hamiltonian corresponding to the average potential U (r) of the system is integrable. It is nevertheless uneasy in general to toggle from one to the other space, although this is certainly possible for spherically symmetric potentials, for flat systems (which may however be unstable) and for special thick disk potentials. Numerical methods could be used to achieve the necessary transformation (Pichon & Cannon 1997;McMillan & Binney 2008). As an illustrative example, we shall give special attention to spherically symmetric potentials, expanding their kinetic equation into a system which almost entirely avoids any calculation in the position-momentum space. The system's inhomogeneity requires that solutions to the Poisson equation are easily found for any inhomogeneous mass distributions. This is achieved by projecting on a biorthogonal basis of pairs of density-potential functions.

Many astrophysical systems which have evolved to a quasi-stationary collisionless equilibrium still keep evolving on time scales longer than the dynamical time as a result of gravitational noise induced by their own constituents or by external ones. We disregard external perturbators, which we define as unbound to the system, although, as did Weinberg (2001b), these could be treated, if numerous and frequent enough, as a given, non-evolving, population providing a source of gravitational noise for other populations. Loosely bound satellites or remote star populations are regarded as internal to the system. This is possible because our set of kinetic equations allows to simultaneously follow different mass populations. Dwarf satellite galaxies could be regarded for example as one such mass population. Globular clusters, dwarf galaxies, disk galaxies and their haloes are examples of bound systems still evolving as a result of internal noise caused by particle discreteness. Such systems are the object of our study. As in any weakly coupled system, the particles suffering collisions are dressed by the polarization clouds caused by their own influence on other particles. Collisions between dressed particles have quantitatively different outcomes than collisions between naked ones (Weinberg 1998). This may reflect in significant differences in calculated effective relaxation times and braking or diffusion coefficients, especially when the system, though stable, is not too far from instability (Weinberg 1993). It is therefore useful to account for collective dressing when calculating such processes as secular thick disk evolution, mass segregation in galaxies or in star clusters, or the damping by dynamical friction of galactic populations on high energy orbits. For simplicity, the kinetic equations to be derived below assume that the system is stationary on a dynamical time scale. They thus cannot address questions in which the distribution in angle variable matters, such as the dissolution of freshly accreted satellites, although a simple extension of the theory could. Since however our equations describe the coupled evolution of all populations present in the system, they are well suited to study, for example, the simultaneous evolution by dynamical friction and diffusion of a stellar population and the population of molecular clouds off which these stars scatter.

The collective response of a self gravitating system to the presence of a perturbing body has been considered by a number of authors, analytically [START_REF] Weinberg | [END_REF](Weinberg , 1995;;Murali & Tremaine 1998;Saha & Jog 2006) or numerically (Thielheim & Wolff 1984;[START_REF] Gnedin | [END_REF]. Sometimes, the reaction of this perturbation on the perturbing body itself is calculated, as did [START_REF] Kalnajs | IAU Colloq. 10, Gravitational N-body problem[END_REF], who computed the drag on a large body moving in an homogeneous medium, taking the collective response of this medium into account, and Tremaine & [START_REF] Weinberg | [END_REF], who considered the global, self-consistent, perturbation caused by a satellite or a barred structure in a spherically symmetric system and its reaction on the perturbator object by the effect of dynamical friction. The secular evolution of the system in response to such perturbations has been considered by Weinberg (2001a), who considered general types of perturbations on a galaxy, and by Pichon & Aubert (2006) who considered perturbations caused by the cosmological environment on dark matter haloes. This evolution is of course in principle observable in N-body simulations, which however have their own difficulties in calculating the long term evolution of such systems [START_REF] Balescu | Equilibrium and non equilibrium statistical mechanics[END_REF]. A number of authors (Murali 1999;Weinberg 2001a;Pichon & Aubert 2006) have studied the collective perturbations caused in a massive spherical galactic halo by its environment. They could calculate the response of this system by resorting to a representation of the particle's motion in action and angle variables, a method first used by Kalnajs (1977). We follow them on this road. They also made good use of a basis of biorthogonal pairs of density-potential functions. Weinberg (1993) first derived a kinetic equation for the collisional relaxation of a self gravitating system along these lines. His equation accounts for the self-consistent gravitational dressing of the particles, but is otherwise simplified, the geometry supposedly being that of an homogeneously filled periodic cube. The inhomogeneous nature of the system should be described more accurately, still accounting for collective gravitational dressing effects. This is specifically the aim of this paper. Chavanis (2007) presented a Gravitational Balescu-Lenard equation 3 similar approach to ours for one-dimensional systems, the constituents of which interact by a general long range force. In this paper we further elaborate in section 8 on the structural evolution of the inhomogeneous system and on the secular evolution of the orbits.

CUTTING THE BBGKY HIERARCHY

Reduction of the hierarchy to a kinetic equation

The Liouville equation for the N-body distribution function of a system of interacting particles can be translated into a hierearchy of equations, the BBGKY hierarchy, for the reduced 1-body, 2-body, 3-body etc .. distribution functions (Balescu 1963;[START_REF] Binney | Galactic Dynamics[END_REF]. The equation for the 1-body distribution function also involves the 2-body distribution, the equation for the 2-body distribution involves the 3-body distribution and so on. The kinetic equation being meant to be an autonomous equation for the 1-body distribution f1(r, p, t), its derivation necessarily involves some approximation allowing to cut this hierarchy. This is usually done at the level of the equation of evolution of the 2-body distribution function, reducing it to a relation between the 2-body and the 1-body distribution functions. The simplified equation for the 2-body distribution f2(r1, p1, r2, p2, t) is then solved in terms of the 1-body distribution f1(r1, p1, t) and the result, once introduced in the first equation of the hierarchy, provides the desired kinetic equation for f1.

Plasma physics knows of two such successful approximations: rare and short range interactions, allowing to ignore 3-body collisional effects on the evolution of the 2-body distribution function, leading to the Boltzmann equation [START_REF] Uhlenbeck | Lectures in statistical mechanics[END_REF] and weakly coupled, collective, systems in which the 3-body correlations may be neglected and the 2-body correlations considered weak, leading to the Balescu-Lenard equation [START_REF] Peyrissac | Equations cinétiques des fluides et des plasmas[END_REF]Lenard 1960). The weak correlation approximation is valid when the number of particles in the effective interaction sphere, the Debye sphere, is large. This approximation is also valid for self-gravitating systems with a large number N of simultaneously interacting particles. The coupling in this case is indeed weak, the ratio of the average interaction energy to the average kinetic energy scaling as N -2/3 . This provides a solid basis for the derivation of a kinetic equation. The larger N , the more valid the approximation is. For systems with a very large number of bodies, the resulting kinetic equation is almost exact, but for the description of strong collisions.

The constituents of the system are considered to be point-like objects of different masses, which we refer to as particles. They need not all be stars, but could be other entities as well, such as molecular clouds, bound clusters, a population of satellites or lumps of dark matter in the halo of a galaxy. The kinetic equations to be derived below are valid as long as most collisions are weak, which implies that the collisional evolution time of any type of particles remains long compared to the dynamical time. We assume that the masses of the consituents come in a finite set. Each mass group is labeled by a lower case latin letter.

Notations

An efficient and concise notation is needed. Some weakly relevant variables, such as time, will often be omitted from the list of arguments of some functions. The subscripts 1 or 2 on one-or two-body distributions or correlation functions will also be omitted, the number of arguments indicating the number of bodies involved. The 1-body distribution function of particles of species a (that is, of mass ma) is denoted by f a , the 2-body distribution function of a pair of particles of species a and b (where a and b may be equal or different) is f ab and the corresponding 2-body correlation function is g ab = f ab -f a f b . The space and momentum integral of a 1-body distribution function is the total number of particles of the considered species. Similarly, the space and momentum integral of 2-body distribution functions is the total number of pairs of the considered species. When a = b, pairs should be regarded as ordered entities. The position and momentum (r1, p1) of a particle is simply noted 1, for brevity. The three angle and three action variables of this particle similarly form a pair of vectors (w1, J1). The same shorthand notation, 1, is used where the context commands. The notation d1 represents either d 3 r1d 3 p1 or d 3 w1d 3 J1. These phase space volume elements are equal because both sets of variables are canonical. The velocity of particle 1 is v1. The gradient with respect to a vectorial variable u, like r, p, w or J, is denoted by ∇u. The derivative with respect to time is noted ∂t.

G being Newton's constant, the gravitational force suffered by a particle of species a with dynamical variables 1 (that is at r1 with momentum p1) from a particle of species b with dynamical variables 2 is:

F ab (1, 2) = Gmam b r2 -r1 | r2 -r1 | 3 • (1)
We ignore any external force, be it tidal or exerted by some closeby external body. The collective gravitational force F 0 a (1) exerted at r1 on a particle of species a is the 1-body and species average of F ab (1, 2):

F 0 a (1) = b d2 F ab (1, 2) f b (2) • (2)
The gravitational potential U (r1) from which this force derives is:

U (r1) = - b d2 Gm b | r2 -r1 | f b (2) • (3)

Weak correlations in terms of one-body distributions

The first equation of the BBGKY hierarchy can be written:

∂tf a (1) + v1 • ∇r 1 f a (1) + F 0 a (1) • ∇p 1 f a (1) = - b d2 F ab (1, 2) • ∇p 1 g ab (1, 2) • (4)
Neglecting 3-body correlations, the second equation of the BBGKY hierarchy can be written as:

∂tg ab (1, 2) + v1 •∇r 1 + v2 •∇r 2 g ab (1, 2) + F 0 a (1)•∇p 1 + F 0 b (2)•∇p 2 g ab (1, 2) + c d3 g bc (2, 3) Fac(1, 3)•∇p 1 f a (1) + c d3 g ac (1, 3) F bc (2, 3)•∇p 2 f b (2) = F ab (1, 2) • (∇p 2 -∇p 1 )f a (1)f b (2)•(5)
Equation ( 5) is linear in the correlation function and has on its right hand side a source term S ab (1, 2, t) which is a functional of the 1-body distribution functions, namely:

S ab (1, 2, t) = F ab (1, 2) • (∇p 2 -∇p 1 )f a (1)f b (2) • (6)
The solution for g ab (1, 2, t) can be found in terms of the sources S by working out the Green's function, or propagator, of the operator on the left hand side of equation ( 5). This Green's function is a matrix in particle species space, G ab pq (1, 2, 1 ′ , 2 ′ , τ ), in terms of which the correlation function can be expressed as:

g ab (1, 2, t) = p,q ∞ 0 dτ d1 ′ d2 ′ G ab pq (1, 2, 1 ′ , 2 ′ , τ ) S pq (1 ′ , 2 ′ , t -τ ) • (7)
Equation ( 7) expresses the correlation function as a functional g ab (1, 2; f ) of the 1-body distributions. Once the 2-body propagator has been found, the solution (7) for g ab (1, 2) may be substituted on the right hand side of equation ( 4), which then depends explicitly, and only, on the 1-body distributions. We call it the collision operator C a (f ) for species a:

C a (f ) = - b d2 F ab (1, 2) • ∇p 1 g ab (1, 2; f ) • (8)
The initial value of the 2-body propagator is:

G ab pq (1, 2, 1 ′ , 2 ′ , 0) = δ a p δ b q δ(1 -1 ′ ) δ(2 -2 ′ ) , (9)
where δ(1 -1 ′ ) is a Dirac function and δ a p a Kronecker symbol. By substituting equation (7) in equation ( 5), it can be shown that the 2-body propagator can be factored into the product of two 1-body propagators:

G ab pq (1, 2, 1 ′ , 2 ′ , τ ) = G a p (1, 1 ′ , τ ) G b q (2, 2 ′ , τ ) • (10)
Had we considered strong interactions as well, the correlation function g ab (1, 2) would not have been negligible compared to f a (1)f b (2) and the right hand side term of equation ( 5) would have been changed by the substitution of f a f b +g ab to f a f b . The 2-body propagators would in this case not factor as in equation ( 10). In the weak correlation approximation considered here, the 1-body propagators G a p (1, 1 ′ , τ ) satisfy the linearized Vlasov equations:

∂τ G a p (1, 1 ′ , τ ) + (v1 •∇r 1 + F 0 a (1)•∇p 1 ) G a p (1, 1 ′ , τ ) + c d2 G c p (2, 1 ′ , τ ) Fac(1, 2) • ∇p 1 f a (1) = 0 , (11)
with initial condition G a p (1, 1 ′ , 0) = δ a p δ(1 -1 ′ ). The solution of equation ( 11) has to be found for τ ≥ 0 only, because of causality. According to Bogoliubov's synchronisation hypothesis (Bogoliubov 1946), the 1-body distribution functions can be regarded as constant in equations ( 7) and ( 11) because they evolve on the relaxation time scale, which is much longer than the time required for the correlation function to reach an equilibrium, given the present value of the 1-body distributions. The correlations at a given time t then are functionals of the one particle distribution functions at this very same time.

Equation ( 11) can be solved by means of a Laplace transform with respect to the time lapse τ . The Laplace transform f (ω) of a function of time f (t) depends on a complex argument ω. The transformation and its inverse are defined by:

f (ω) = ∞ 0 f (t) e iωt dt and f (t) = 1 2π B f (ω) e -iωt dω • (12)
The direct transform is convergent only when the imaginary part of ω exceeds some ordinate of convergence, above which the function f (ω) is regular. Below it, it is defined by analytical continuation. The Bromwich contour B which appears in the inverse transformation runs parallel to the real axis from -∞ to +∞ above all singularities of f (ω). Equation ( 11) is Laplace-transformed into:

-iω G a p (1, 1 ′ , ω) + (v1 •∇r 1 + F 0 a (1)• ∇p 1 ) G a p (1, 1 ′ , ω) + c d2 G c p (2, 1 ′ , ω) Fac(1, 2) • ∇p 1 f a (1) = δ a p δ(1 -1 ′ ) • (13)

PARTICLE MOTIONS AND BASIS FUNCTIONS IN ANGLE AND ACTION VARIABLES

Angle and action variables

The particle motions in the self gravitational field are complex in general. This precludes a direct solution of equation ( 11) by integration along unperturbed trajectories. It is preferable to change the position and momentum variables for a set of canonical angle and action variables [START_REF] Goldstein | Classical Mechanics[END_REF]. So doing, the description of the motion becomes simple, all the complexity being embodied in the relation between position and momentum variables and angle and action variables. By definition, the Hamiltonian H in angle and action variables depends only on the three actions J1, J2, J3, which we regard as the three components of an action vector J. The three actions are constants of the motion and the three angles w1, w2, w3 which similarly form the components of an angle vector w, vary linearly in time. The angular frequency of the angle wi is Ωi = ∂H/∂Ji. The frequencies Ωi form the components of a frequency vector Ω which depends on J. For brevity, we use shorthand notations, such as:

Ω1 ≡ Ω(J1) , Ω ′ 1 ≡ Ω(J ′ 1 ) • (14)
The derivative following the motion is (v

• ∇r + F 0 • ∇p).
The actions being first integrals, this operator translates in angle and action variables into (dw/dt) • ∇w, that is, v • ∇r + F 0 • ∇p = Ω • ∇w The last, collective, term of the left hand side of equation ( 13) must be expressed in action and angle variables. It is of the frequently met general form:

λ d2 M (2) r2 -r1 | r2 -r1 | 3 • ∇p 1 N (1) • (15)
The force in equation ( 15) can be expressed in terms of a "potential" φ such that:

d2 M (2) r2 -r1 | r2 -r1 | 3 = -∇r 1 φ(r1) • (16)
This potential and the "mass distribution" D from which it derives, depend on the function M (2) only. They are defined by:

φ(r1) = -d2 M (2) | r2 -r1 | , D(r2) = d 3 p2 M (2) • (17) Since φ(1) is independent of p1, ∇r 1 φ(1)•∇p 1 N (1) is the Poisson bracket { φ(1), N (1)
}. This bracket being invariant on a change of canonical variables, the expression (15) can be written as:

d2 M (2) r2 -r1 | r2 -r1 | 3 • ∇p 1 N (1) = -∇r 1 φ • ∇p 1 N (1) = -∇w 1 φ • ∇J 1 N (1) -∇J 1 φ • ∇w 1 N (1) • (18)
All functions depend periodically, with period 2π, on the angles, with respect to which a discrete Fourier transform can be made. All components of the associated wave vector k are relative integers. The transform of any function f (w, J) and the inverse transform are defined by:

f (w, J) = k f k (J) e i k•w and f k (J) = d 3 w 8π 3 f (w, J) e -i k•w • (19)
Each integral in the second term of equation ( 19) is over the 2π period of one of the components of w. The transform of the Dirac function δ(w) is 1/8π 3 and the transform of unity is δ(k), where δ is here a triple Kronecker symbol. The position r of a particle is a function of its angle and action variables, w and J. The simple Fourier transforms with respect to the angles w1 of ψ α (r1) and G a p (1, 1 ′ , ω) and the double Fourier transform of the propagator with respect to angles w1 and w ′ 1 are:

ψ α (1) ↔ ψ α k 1 (J1) G a p (1, 1 ′ , ω) ↔ G ap k 1 (J1, 1 ′ , ω) G a p (1, 1 ′ , ω) ↔ G ap k 1 k ′ 1 (J1, J ′ 1 , ω) • (20)

Biorthogonal density-potential bases

A basis of biorthogonal density-potential pairs is effective in calculating the potential φ(1) defined by eq. ( 17). Many such bases have been proposed [START_REF] Kalnajs | [END_REF]Clutton-Brock 1972, 1973;[START_REF] Kalnajs | [END_REF]Aoki & Iye 1978;Aoki, Noguchi & Iye 1979;Saha 1991;Hernquist & Ostriker 1992;Robijn & Earn 1996;Brown & Papaloizou 1998;Rhamati & Jallali 2009). A basis element is labeled by a greek letter. The dummy index rule is used for these basis indices. Let D α (r) and ψ α (r) be the density and the potential of the element α of the basis. The potential ψ α derives from the density distribution D α and is related to it by:

ψ α (r) = -d 3 r ′ D α (r ′ ) | r ′ -r | • ( 21
)
The basis is biorthogonal and normalized, such that:

d 3 r D α (r) (ψ β (r)) * = -δ α β • (22)
The symbol on the right of equation ( 22) is a generalized Kronecker. The minus sign results from the fact that when α = β the left hand side of equation ( 22) necessarily is negative. The functions to be expanded on the basis being real, the complex conjugates of D α and ψ α ,

D α(r) ≡ (D α (r)) * and ψ α(r) ≡ (ψ α (r)) * , (23)
also form an element α of the basis, which in general is different from α. The variable r being a length and the Kronecker δ in ( 22) being dimensionless, equations ( 21) and ( 22) imply that D α and ψ α have dimensions L -5/2 and L -1/2 respectively. Any density distribution D(r) and its associated potential φ(r) can be expanded on the basis as:

D(r) = aα D α (r) φ(r) = aαψ α (r) • (24)
The basis functions ψ α (r) are not real in general, which implies that ψ α -k = (ψ α k ) * . The notation ψ α * k denotes the complex conjugate of ψ α k . The notation ψ α k is adopted for the k-Fourier transform of the function ψ α(r) ≡ (ψ α (r)) * . In general,

ψ α k = ψ α * k .
Complex conjugation implies however that:

ψ α k = (ψ α -k ) * • (25)
The coefficients of the expansions (24) can be calculated by using the biorthogonality relation ( 22) and expressed in angle and action variables by using the density-potential basis and angle Fourier coefficients. In particular, the coefficient aα associated with the density field of equation ( 17) is:

aα = -d2 M (2) (ψ α (2)) * = -8π 3 k d 3 J2 M k (J2) ψ α * k (J2) • (26)
The expression ( 15) is transformed in angle and action variables by using the density-potential basis and angle Fourier coefficients into:

λ d2 M (2) r2 -r1 | r2 -r1 | 3 • ∇p 1 N (1) = -λ k 1 aα e ik 1 •w 1 ψ α k 1 (1) ik1 •∇J 1 N (1) -(∇J 1 ψ α k 1 (1))•∇w 1 N (1) , (27)
where the expansion coefficients aα are given in terms of the function M (2) by equation ( 26). Other expressions of aα will be established below for the case when M (2) is the propagator, as in equation ( 13).

THE LINEARIZED VLASOV PROPAGATOR

In a relaxing system, a collisionless equilibrium is supposedly reached on a time scale shorter than the relaxation time, so that the system is stationary on the dynamical time scale. This means that the distributions f a (1) really are functions f a (J1) of the actions only. The third, "collective", term on the left of equation ( 13) is of the form displayed in equation ( 15). The corresponding factor λ and functions N (1) and M (2) particularize in this case to λ = Gmamc, N (1) = f a (1) and

M (2) = G c p (2, 1 ′ , ω).
For these functions N (1) and M (2), the coefficients aα of equation ( 26) are:

a cp α (1 ′ , ω) = -8π 3 k d 3 J G cp k (J, 1 ′ , ω) ψ α * k (J) • (28)
Species-cumulative coefficients A p α are defined by:

A p α (1 ′ , ω) = a ma a ap α (1 ′ , ω) • (29)
Equation ( 13) for the 1-body propagators is Fourier transformed with respect to w1 using equation ( 27), which gives:

G ap k 1 (J1, 1 ′ , ω) = i 8π 3 δ a p δ(J1 -J ′ 1 ) e -i k 1 •w ′ 1 ω -k1 • Ω1 -Gma (k1 •∇J 1 f a (1)) ω -k1 • Ω1 ψ β k 1 (J1) A p β (1 ′ , ω) • (30)
The coefficients A (eq. ( 29)) can be expressed in terms of the Fourier transform of the propagators by using equation ( 28):

A p α (1 ′ , ω) = -8π 3 c k d 3 J mc ψ α * k (J) G cp k (J, 1 ′ , ω) • (31)
Operating on equation (30) as on the function G in equation ( 31), a linear system is obtained for the species-cumulative coefficients A. It can be written:

ε αβ (ω) A p β (1 ′ , ω) = σ p α (1 ′ , ω) , (32)
σ p α (1 ′ , ω) = -i mp k 1 ψ α * k 1 (J ′ 1 ) e -i k 1 •w ′ 1 ω -k1 •Ω ′ 1 , (33)
ε αβ (ω) = δ αβ - a k 1 d 3 J1 8π 3 Gm 2 a ψ α * k 1 (1) ψ β k 1 (1) k1 •∇J 1 f a (1) ω -k1 •Ω1 • (34)
The solution of equation ( 32), obtained by inverting the matrix ε αβ (ω), is then introduced in equation ( 30), giving the Fourier and Laplace transform of the propagator. So doing, a function D appears in the solution, which is defined by:

1 D k 1 k ′ 1 (J1, J ′ 1 , ω) = ψ α k 1 (J1) ε -1 (ω) αβ ψ β * k ′ 1 (J ′ 1 ) • (35)
Performing the inverse Fourier and Laplace transforms of equation ( 30), the 1-body propagator itself is eventually found:

G a p (1, 1 ′ , τ ) = B dω 2π e -i ωτ k 1 k ′ 1 i e i(k 1 •w 1 -k ′ 1 •w ′ 1 ) 8π 3 (ω -k1 • Ω1) δ a p δ(k1-k ′ 1 ) δ(J1 -J ′ 1 ) + 8π 3 Gmamp (k1 •∇J 1 f a (1)) (ω -k ′ 1 •Ω ′ 1 ) D k 1 k ′ 1 (J1, J ′ 1 , ω) • (36)
5 THE KINETIC EQUATION

Explicit writing of the kinetic equation

The correlation function is obtained from the solution (36) for the 1-body propagator by using equations ( 10) and ( 7). The kinetic equation and its collision operator are then given by equations ( 4) and ( 8). Thanks to the Bogoliubov synchronisation hypothesis, this equation is local in time, because the source term S pq (1 ′ , 2 ′ , t -τ ) in equation ( 7) can be regarded as independent of τ and equal to its value at τ = 0. The collision operator for the evolution of the distribution function of species a, C a (f ), is defined by equation ( 8) and can be written as:

C a (f ) = - p, q ∞ 0 dτ d1 ′ d2 ′ B dω 2π B ′ dω ′ 2π e -i (ω+ω ′ )τ b d2 • • • • • • (F ab (1, 2) • ∇p 1 ) G a p (1, 1 ′ , ω) G b q (2, 2 ′ , ω ′ ) (Fpq(1 ′ , 2 ′ ) • (∇ p ′ 2 -∇ p ′ 1 )) f p (1 ′ , t)f q (2 ′ , t) • (37)
The somewhat lengthy transformations that must be performed to express this equation in terms of the angle and action variables, of the density-potential basis and of its angle Fourier transforms are described in appendix A. They eventually yield the following final form of the kinetic equations:

∂tf a (J1) = b k 1 k 2 d 3 J2 8π 4 G 2 m 2 a m 2 b k1 •∇J 1 δ(k1 •Ω1 -k2 •Ω2) |D k 1 k 2 (J1, J2, k1 •Ω1)| 2 (k1 •∇J 1 -k2 •∇J 2 ) f a (J1)f b (J2) , ( 38
)
where D is defined by equation ( 35) and the response matrix elements ε αβ needed to determine D are expressed in terms of the 1-body distribution functions by equation (34). No convective term Ω1•∇w 1 f a (1) appears on the left hand side of eq.( 38) because in a slowly relaxing system the distribution functions f a (1) are meant to depend only on the actions.

Physical content of the kinetic equation

Equation ( 38) describes the relaxation of the distribution functions caused by the, supposedly weak, noise created by the discreteness of the particles accompanied by their associated gravitational polarization cloud (Weinberg 1998;Rostoker & Rosenbluth 1960). This is shown by working out the Fokker-Planck equation for the evolution of actions of the particles in this random field. The potential of a mass m2 with action-angle variables J2, w2 on a particle 1 with action-angle variables J1, w1 is:

Ũ2(1, t) = - k 1 k 2 Gm2 exp(i(k1 •w1 -k2 •w2)) D k 1 k 2 (J1, J2, k2 •Ω2) • (39)
The fluctuating part of the potential created by the discreteness of the dressed particles is the sum over all particles 2 and all non-vanishing k1, k2 of potentials like (39). The rate of change of the action J1 of a particle 1 in this fluctuating field is:

J1 = 2 k 1 =0 k 2 =0 Gm1m2 ik1 exp(i(k1 •w1 -k2 •w2) D k 1 k 2 (J1, J2, k2 •Ω2) • (40)
The braking and diffusion coefficients of the corresponding Fokker-Planck equation are obtained from, respectively, the first and second moments of the random change ∆J1 suffered by the particle 1 in a time ∆t. The averaging is performed on the values of the angle variables of particles 2 and on their action distribution functions. Equation ( 38) is recovered that way. In the calculation of the braking coefficient, small departures from uniform angular motion should be accounted for, as shown by [START_REF] Ecker | Theory of Fully ionized plasmas[END_REF] in a similar context. The Fokker-Planck form of equation ( 38), although equivalent to it, looks more complex than equation ( 38) itself because the braking coefficient involves the derivative of a Dirac distribution.

Accounting for strong collisions

Equation ( 38) and its quasi-homogeneous limit, equation ( 47), both result from a weak collision theory. Strong collisions involving substancial deviation of at least one of the colliding particles are not adequately described. This inappropriate description of the rare strong collisions can be fixed by limiting the range of impact parameters to values larger than some critical limit b ab cr which depends on the masses of the colliding species. This critical impact parameter for particles of species a and b is such that the typical kinetic energy in their relative motion be equal to their interaction energy, that is:

GM R mam b ma + m b = Gmam b b ab cr • ( 41
)
Here M is the total system's mass and R a typical global size of it. Were this cut to be omitted, the expressions of the coefficients in equations ( 38) and (47) would diverge logarithmically at large wavenumbers, where the response function ε approaches unity. This divergence results from the neglect of large deviations in strong collisions. A physically sound result is obtained by limiting the K integration in equation ( 47) to the domain | K |< K ab cr where:

K ab cr = 2π b ab cr = 2π R M ma + m b • (42)
Similarly the summations on the angle Fourier variables ki (i = 1 or 2) in equation ( 38) should be limited, in the term associated to species a and b, to values such that the physical wavenumbers along the quasi-intersecting orbits be smaller than K ab cr . This modulus of the physical wavenumber can be crudely related to the dimensionless angle wavenumber by K = k/R, where R is a typical global size of the system and k the modulus of the angle Fourier variable. Thus, the summation on k1 and k2 in equation ( 38) should be limited to wave vectors, the modulus of which is bounded by:

| ki |< 2π M ma + m b • (43)
When solving equation ( 38), the secular evolution of the response matrix ε, the system's collective potential U (r, t) and the Fourier transform coefficients ψ α k (J) should be followed in time together with the 1-body distributions. We return to this, in the case of spherical potentials, in section 8. Prior to that, let us discuss various limits and approximate forms of equation ( 38) and show that, as it should, it satisfies an H-theorem. The irreversibility stemms from the fact that information is lost when the real issues of collisions are replaced in the equation by average ones, in particular by averaging over the angles of the colliding particles.

6 LIMITING CASES

Homogeneous limit

Although the limit of an homogeneous medium cannot be rigourously taken for a self gravitating system, it is nevertheless possible to assume local homogeneity at the price of artificially limiting the interaction distance between particles by cutting it at some characteristic size of the system. So doing, the effects of the collective dressing of the particles are still retained, albeit less precisely, but the effects of the structure of the system are only sketchily accounted for.

In this limit, the system is regarded as homogeneously filling a large cubix box of side L, on the surface of which periodic boundary conditions apply. This is the geometry considered by Weinberg (1993). Due to the assumed homogeneity, the collective force F 0 vanishes and the unperturbed motion is rectilinear and uniform, whatever the state of relaxation of the system. The action variables are then proportional to the components of the momentum p and the angle variables are proportional to the components of the position r. Since the angles must be variables of period 2π, the angle vector must be w = 2πr/L, which implies that the action vector is J = Lp/2π. The angle Fourier vector is k = LK/2π where K is the usual wave vector of Fourier transforms with respect to position. The frequency vector Ω is 2πv/L, so that k

• Ω = K • v.
The density-potential basis consists of functions proportional to complex exponentials, like exp(i K•r). A given element of the basis, α say, is characterized by its wave vector K. This can be accounted for in the notation by writing this wave vector as Kα, the corresponding angle wave vector being noted kα. The density function and the potential of the element α of the basis are both proportional to exp (i Kα •r). Their normalization factor must be such that the biorthogonality relation ( 22) be satisfied, the density D α (r) and the potential ψ α (r) being related by equation ( 21). These constraints result in:

D α (r) = | KαL | 2 √ π L 5/2 e i Kα•r ψ α (r) = - 2 √ π L 1/2 e i Kα•r | KαL | • (44)
The ψ α k 's are the Fourier transforms of ψ α (r) with respect to the angles w, namely:

ψ α k = - 2 √ π L 1/2 δ(kα -k) | KαL | , (45)
where δ(kαk) is a triple Kronecker symbol. In this case the ψ α k 's do not depend on the actions and remain fixed while the relaxation proceeds. The response matrix ε, calculated from its definition (34), is diagonal, its element αα being given, for ω in the upper half complex plane, by:

ε αα (ω) = 1 - q 4π G m 2 q | Kα | 2 d 3 p Kα •∇pf q (p) ω -Kα •v • (46)
For real ω, a +i0 should be added to the singular denominator. Since α enters this relation by its wave vector Kα, ε αα (ω) can be regarded as a function ε(Kα, ω), or, more generally, as a function of a wave vector K and of a frequency ω. Because of the diagonality of the response matrix ε and the simplicity of equation ( 45), the writing of the kinetic equation ( 38) simplifies to:

∂tf a (p) = b d 3 p ′ ∇p • Q ab (p, p ′ ) • (∇p -∇ p ′ )f a (p)f b (p ′ ) , (47)
where the tensor Q ab is defined by:

Q ab (p, p ′ ) = 2G 2 m 2 a m 2 b d 3 K KK K 4 δ(K•(v -v ′ )) | ε(K, K•v) | 2 • (48)
Equations ( 47)-( 48) are identical to equation (29) of Weinberg (1993) when the quasi homogeneity of the system and associated absence of collective and external forces are accounted for. Equation ( 47) can be written explicitly as a Fokker-Planck equation in the form:

∂tf a (p) = -∇p • Aaf a (p) + 1 2 ∇p∇p : Baf a (p) , (49)
where the momentum drag and diffusion coefficients are:

Aa(p) = b d 3 p ′ f b (p ′ ) (∇p -∇ p ′ ) • Q ab (p, p ′ ) , (50)
Ba(p) = 2 b d 3 p ′ f b (p ′ ) Q ab (p, p ′ ) • (51)
For electrical instead of gravitational interactions, the gravitational constant G should be replaced by 1/4πǫo in MKSA units, ǫo being the dielectric permittivity of vacuum. The electric force between like charges being repulsive instead of attractive, the minus sign before the second term of equation ( 46) should be changed to a positive sign and the masses replaced by the charges of the particles. Equation ( 47) then reduces to the Balescu-Lenard equation for homogeneous and multispecies plasmas [START_REF] Peyrissac | Equations cinétiques des fluides et des plasmas[END_REF]. It implicitly accounts for the screening effect, which is embodied in the dielectric function. The integral on wave vector space in equation ( 48) then need not be cut at small wave vectors because | ε(K, K• v) | diverges as | K | approaches zero. For self gravitational systems the small | K | limit is unphysical, due to the absence of screening. The distance between interacting particles is limited in this case by the inhomogeneity of the system, a feature which is lost in the local approximation. If one were to insist on the quasi-homogeneous approximation, the integration over wave vectors in equation ( 48) would have to be artificially limited from below to some minimum modulus Kmin ∼ 2π/R, where R is a characteristic size of the system. Little would then be gained over a more traditional Fokker-Planck approximation, but for the fact that equation (47) still accounts for the collective dressing of the colliding particles.

Non-collective homogeneous limit

When these collective effects are themselves neglected, which amounts to take ε = 1 in equation ( 48), the usual local Fokker-Planck equation ( 49) is recovered, with braking and diffusion coefficients given by expressions ( 50) and ( 51), ε now supposedly being equal to unity. As above, the integral on wavevectors in equation ( 48) should limited to a lower cutoff at | K |= Kmin, to account for the finite size of the system, and to an upper cutoff | K |= Kmax, to account for strong collisions (section 5.3). The coulomb logarithm is ln Λ, where Λ = Kmax/Kmin. When ε equals unity, the integration over wave vectors in equation ( 48) can easily be performed. The result, which involves the relative velocity of the colliding particles g = vv ′ , is:

Aa(p) = -4πG 2 ln(Λ) b mam b (ma + m b ) d 3 p ′ f b (p ′ ) g g 3 , (52)
Ba(p) = +4πG 2 ln(Λ) b m 2 a m 2 b d 3 p ′ f b (p ′ ) I g 2 -gg g 3 , ( 53
)
where I is the unit second rank tensor. This is identical to the Fokker-Planck equation presented, for example, in [START_REF] Binney | Galactic Dynamics[END_REF], equations (8A.10). The collective dressing becomes important when | ε | -2 in equation ( 48) largely differs from unity. As shown by Weinberg (1993), this happens when the system is not far from being unstable, for example when its size becomes of order of the Jeans length, the complex zeroes of ε(K, ω) lying close, but below, the real axis.

AN H THEOREM

Equation ( 38) satisfies an H theorem which states that the statistical entropy:

s(t) = - a d 3 w1 d 3 J1 f a (w1, J1, t) ln (f a (w1, J1, t)) , (54)
increases with time. Because the relaxing distribution functions depend on actions only, the integral over angles reduces to a mere multiplication by a factor 8π 3 , so that:

ds(t) dt = -8π 3 a d 3 J1 1 + ln (f a (J1, t)) ∂tf a (J1, t) • (55)
The time derivative of f a is given by equation ( 38) which can be symmetrized by substituting to the first operator k1•∇J 1 the operator k1•∇J 1 -k2•∇J 2 . The contribution associated with the added operator k2•∇J 2 vanishes on integration over J2. This can be seen by using the flux divergence theorem in action space, recognizing that the surface integral over the boundary of the physical J2 domain vanishes. Indeed, the expression on the right of the first operator k1•∇J 1 in equation ( 38) represents, up to its sign, the flux in action space at J1 caused by collisions with particles having action J2 or the flux at J2 caused by collisions with particles of action J1. The physical domain is limited in action space by a boundary at a finite distance and extends to infinity in certain directions. The flux through the boundary at finite distance vanishes because the action vector of no particle can evolve through this boundary from the physical to the unphysical domain. The flux at infinity vanishes because f b (J2) decreases fast enough. This justifies the above-suggested substitution. The expression of ∂tf a (1) given by equation ( 38), modified as described, when inserted in equation ( 55), gives the following expression for ds/dt:

ds(t) dt = -64π 7 a b k 1 k 2 d 3 J1 d 3 J2 G 2 m 2 a m 2 b 1 + ln(f a (J1)) • • • • • • (k1 •∇J 1 -k2 •∇J 2 ) δ(k1 •Ω1 -k2 •Ω2) |D k 1 k 2 (J1, J2, k1 •Ω1)| -2 (k1 •∇J 1 -k2 •∇J 2 ) f a (J1)f b (J2) • ( 56
)
This expression is further symmetrized by combining it with the equivalent expression obtained by exchanging species indices a and b, actions J1 and J2 and Fourier variables k1 and k2. The resulting expression is then integrated by parts, using the flux divergence theorem in either J1 or J2 space. As explained above, the contribution of the integral on the boundary of the action domain or at infinity vanishes. We are eventually left with the positive expression:

ds(t) dt = +32π 7 a b k 1 k 2 d 3 J1 d 3 J2 G 2 m 2 a m 2 b δ(k1 •Ω1 -k2 •Ω2) (k1 •∇J 1 -k2 •∇J 2 ) f a (J1)f b (J2) 2 |D k 1 k 2 (J1, J2, k1 •Ω1)| 2 f a (J1)f b (J2) • ( 57
)
This establishes that the statistical entropy of the system is monotonically increasing. Since the entropy of a self gravitating system of a given total mass and energy is not bounded from above [START_REF] Binney | Galactic Dynamics[END_REF], the increase of the statistical entropy does not lead, as in homogeneous gases or plasmas, to a state of thermodynamic equilibrium. When the system becomes sufficiently centrally condensed, a gravothermal instability develops [START_REF] Hénon | [END_REF]Antonov 1962;Lynden-Bell & Wood 1968). The quenching of this instability by the formation of binaries is not described by equation ( 38), because the formation of binary systems results from triple collisions (that is, from third order correlations) and is of the strong interaction type.

EVOLUTION OF A SPHERICAL POTENTIAL AND BASIS FOURIER COEFFICIENTS

When the relaxation proceeds, the distribution functions f a (J, t) evolve according to equations ( 38). This causes a slow secular change in the average potential U (r, t) and in the response matrix ε(ω) (equation ( 34)). The basis potential functions ψ α (r) are time-independent, but their Fourier transforms with respect to angles w are not because they depend on the orbits of the particles, which slowly change with the potential. The Fourier coefficients ψ α k (J) are then indirectly related to the slowly evolving potential U (r, t).

Thus, equation ( 38) is not an autonomous equation for the distribution functions f a (J, t). The response matrix ε(ω) of equation ( 34) is a functional of those, which also depends on the angle Fourier transforms ψ α k (J, t) as does the quantity D present in equations ( 35) and ( 38). The kinetic equation ( 38) must then be completed with equations describing the evolution in time of the average potential U (r, t) and of the angle Fourier coefficients ψ α k (J, t) of the basis potentials. This aspect of the system's evolution is not considered by Chavanis (2007). In this section, the time t will be restaured, though only where necessary, in the list of arguments of functions. The potential U (r, t) derives from the mass density:

D(r, t) = a ma d 3 p f a (r, p, t) • (58)
The corresponding gravitational potential is obtained from its expansion on the biorthogonal density-potential basis by equation ( 24). From appendix B, it is found that its partial time-derivative is:

∂tU (r, t) = - a 8π 3 Gma ψ α (r) d 3 J ∂t f a (J, t) (ψ α 0 (J, t)) * • (59)
Equation ( 59), as well as equations ( 34), ( 35) and ( 38), call for an equation for the time-evolution of the coefficients ψ α k (J, t):

ψ α k (J, t) = d 3 w e -ik•w ψ α (r(w, J)) • (60)
An explicit expression for these coefficients when the potential is spherical is derived in appendix B, which also gives a summary of angle and action variables for a particle moving in a spherical potential. In this case, the coefficients ψ α k (J, t) vanish when the wave vector k has non-vanishing k2 or k3 components. The non-vanishing coefficients depend only on the radial k1 component, hereafter noted k. For conciseness, the variables J, which are mere parameters, are omitted wherever possible. After some calculations, we find that:

ψ α k (t) = 8π 2 Ω1(t) r A (t) r P (t) cos W k (r, t) ψ α (r) dr 2(E(t) -U (r, t)) -J 2 2 /r 2 , where W k (r, t) = r r P (t) k Ω1(t) dr ′ 2(E(t) -U (r ′ , t)) -J 2 2 /r ′ 2 • (61)
When the potential changes, the radial distances rP and rA of the periapse and apoapse change accordingly: the bounds of the integrals in equations ( 61), (B5) and (B2) are time-dependent. These integrals are singular, though convergent, the periapse and apoapse distances being the zeroes of the square root denominator:

q(r, t) = 2(E(t) -U (r, t)) -J 2 2 /r 2 • (62)
These zeroes are simple when the orbit is not circular and merge into a double zero when it is. This latter situation can be dealt with by a limit process, in which simple zeroes rP and rA are made to converge to eachother. We then assume that rP and rA are simple zeroes. An index P or A denotes the value of a function of r at rP or rA respectively. Integrals like:

I(t) = r A (t) r P (t) dr m(r, t) q(r, t) , (63)
or similar ones can be expressed in terms of a variable ξ, the values of which remain constant at the changing periapse and apoapse. This variable is defined by:

r = rP (t) + ξ (rA(t) -rP (t)) • (64)
To each of these two types of radial variables, r or ξ, a time variable, t or τ , can be associated, it being meant that t ≡ τ . This introduces two types of time derivatives: ∂t, which is at constant r, and ∂τ , which is at constant ξ. Ordinary time derivatives with respect to τ and t are identical and are denoted by a dot. Partial time derivatives with respect to τ and t differ and are related by:

∂τ = ∂t + ṙA r -rP rA -rP + ṙP rA -r rA -rP ∂r • (65)
The partial derivatives with respect to r and ξ are simply proportional:

∂ ξ = (rA -rP ) ∂r.
At the periapse or apoapse the function q(r, t) vanishes, whatever τ . Hence, ∂τ q = 0 at these points. Differentiating the equation q(r, t) = 0, ṙP and ṙA are found:

ṙP = - ∂tq(rP , t) ∂rq(rP , t) ṙA = - ∂tq ( rA, t) ∂rq(rA, t) • (66)
The partial derivatives of q(r, t) (equation ( 62)) are:

∂tq(r, t) = 2( Ė -∂tU (r, t)) , ∂rq(r, t) = 2 J 2 2 r 3 -∂rU (r, t) • (67)
It can be checked from equations ( 66) and ( 65) that (∂τ q)(rP , t) = ∂τ q(rA, t) = 0. Equation ( 65) implies that:

∂τ r = ṙP rA -r rA -rP + ṙA r -rP rA -rP , ∂τ q = ∂tq + ṙP rA -r rA -rP + ṙA r -rP rA -rP ∂rq • (68)
The time-derivative of I(t) (equation ( 63)) is found by changing the variable r to ξ:

İ = ṙA -ṙP rA -rP I + r A (t) r P (t) m(r, t) dr q(r, t) ∂τ m m(r, t) - 1 2 ∂τ q q(r, t) • (69)
It is important to note that the last term in the parenthesis of the integral in equation ( 69) is regular since the numerator, ∂τ q, vanishes at rP and rA, where q(r, t) does. The right hand side of equation ( 69) then consists of convergent integrals. When this method is used to calculate the time derivatives of E(t) and Ω1(t) from equations (B2) and (B5), the following results are obtained:

Ė = Ω1 π r A (t) r P (t) dr ∂tU (r, t) q(r) , Ω1 = -Ω1 ṙA -ṙP rA -rP + Ω 2 1 2π r A (t) r P (t)
dr q(r, t) ∂τ q q(r, t)

•

The same method is used to calculate ∂τ W k :

∂τ W k (r, t) = Ω1 Ω1 W k (r, t) + ṙA -ṙP rA -rP W k (r, t) - k1Ω1 2 r r P (t) dr ′ q(r ′ , t) ∂τ q(r ′ , t) q(r ′ , t) • (71)
The angle Fourier coefficient ψ α k (J, t) is given by equation ( 61), which is of a form similar to equation ( 63). Using the general result (69), the time derivative of ψ α k (J, t) is found to be:

ψα k (J, t) = Ω1 Ω1 ψ α k + ṙA -ṙP rA -rP ψ α k -8π 2 Ω1 r A (t) r P (t) ψ α (r) dr q(r, t) sinW k (r, t) ∂τ W k -cosW k (r, t) ψ ′ α (r) ∂τ r ψ α (r) - ∂τ q 2q(r, t) • (72)
Equation ( 72) describes the time-evolution of the Fourier coefficients ψ α k (J, t). The auxiliary τ -derivatives which enter this equation are given by equations ( 71), ( 70), ( 68), ( 66) and ( 67). The other quantities entering equation ( 72) are expressed in terms of the potential by equations ( 61), ( 62), (B4) and (B5). All these relations eventually link the time-evolution of ψ α k to the time-evolution of U (r, t), which is itself described by equation (59). Equations ( 38), ( 59) and (72) form the system of coupled equations for the distribution functions f a (J, t), the average potential U (r, t) and the angle Fourier coefficients ψ α k (J, t) that we have been seeking for in this section. This system involves the auxiliary equations mentioned above, as well as equations ( 34)-(35).

CONCLUSIONS

Kinetic equations for the collisional evolution of the constituents of self-gravitating inhomogeneous systems have been derived. These equations (38) surpass in consistency the usual Fokker-Planck equations ( 49) -( 52) -( 53). The latter are unsatisfactory from a principle point of view, being local and non-collective.

By contrast, the proposed equations fully account for the system's inhomogeneity and for the collective gravitational dressing of the colliding particles.

Equations ( 38) describe the evolution of distribution functions in action and angle space, which is possible when the hamiltonian associated with the average potential is integrable.

Physically, these equations describe the evolution of the distribution functions in action space as a result of the weak gravitational noise caused by the discreteness of the particles, dressed with the polarization clouds that their own gravity induces in the system. This gravitational polarization is accounted for in equation ( 38) in a manner that is fully consistent with the distribution functions, as they are at the moment.

Properties of the kinetic equation

Equation ( 38) is the sum of a second order derivative term with respect to actions and of a first order one. It therefore basically is of the Fokker-Planck type, although it is definitely simpler in the form of expression (38). The diffusion coefficient involved depends on the 1-body distributions themselves, in particular through the factor | D | -2 which represents the effect of the dressing of the colliding particles by the gravitational polarization induced around them by their own influence.

Unlike in electrical plasmas, the polarization dressing in self-gravitational systems does not cause any screening of the interaction, which remains effective even between distant particles. The mutual distance of such particles is limited only by the finite size of the system. Were the gravitational influence of particles on their surrounding to be neglected, the response matrix ε (equation ( 34)) would reduce to unity and the coefficients of the corresponding Fokker-Planck kinetic equation would simply be averages by the distribution functions of functions of velocity, as in equations ( 52) -( 53).

It is apparent from the developments of appendix A, which lead to equation ( 38), that the k component in angle Fourier space of the gravitational polarization response given to a particle has frequency ω = k•Ω. This means that the polarization cloud which accompanies a particle forms a structure in angle space which vary as w -Ωt: it corotates in angle with that particle.

The presence of the Dirac function δ(k1 • Ω1 -k2 • Ω2) in equation ( 38) indicates that particles interact resonantly. This certainly is an important physical property of remote interactions, for which the components of the angle wave vectors k1 and k2 must be small. For closer encounters, the modulus of these wave vectors is larger and the resonance condition k1 • Ω1 = k2 • Ω2 becomes less selective, being more easily satisfied.

The correlation function has been calculated on the basis of a linearized theory, which is justified by the weakness of the average interactions in this many-body system. This means that the trajectories of the particles during the collision are regarded as being the unperturbed trajectories. Similarly, the gravitational polarization cloud around any one of the colliding particles is calculated as if the partner in the collision were not present: equation ( 38) is still a weak collision approximation. A cutoff at small impact parameters is therefore needed to account for the rare strong collisions. Equation (38) takes full account of the inhomogeneity of the system, which is embodied in the dependence of the distribution functions on the actions J's. It requires no artificial cutoff at large impact parameters. The details of the trajectories followed by the particles in the present gravitational potential are also fully accounted for, being implicit in the relations which link the angle and action variables to the position and momentum ones. These relations depend on the actual global gravitational potential of the system, which slowly evolves in time together with the distribution functions.

The density-potential basis functions ψ α (r) are choosen at the beginning of the calculation once and for all, but their angle Fourier transforms ψ α k (J), which depend on the actual trajectories of the particles, change with time because the trajectory of a particle of given actions slowly evolves with the general potential of the system as the relaxation proceeds. As long as it suffers no collision, a given particle keeps its vector J fixed because the actions are adiabatic invariants. Collisions, however, cause a secular evolution of the functions f a (J), which is exactly what equation (38) describes.

The description of particle motions is made simple by the use of action and angle variables. Their complexity is embodied in the supposedly known relation between position and momentum variables and action and angle variables. The usefulness of equation ( 38) is therefore limited to systems for which this relation can be established, either analytically or, possibly, numerically (Pichon & Cannon 1997;McMillan & Binney 2008).

While the relaxation proceeds, the gravitational potential and the orbits of the particles evolve. As a result, the kinetic equation must be completed by evolution equations for the potential and for other relevant quantities. Section 8 establishes, for spherical systems, this set of coupled equations. Weinberg M.D., 2001b, MNRAS, 328, 321 APPENDIX A: DERIVATION OF THE KINETIC EQUATION Equation (37) must be expressed in terms of angle and action variables, using the adopted density-potential basis. The integration on the dynamical state of the particle of species b should be carried out first. The integral over the variables 2 in the second line of eq.( 37) is similar to equation ( 15), with λ = Gmam b , M (2) = G b q (2, 2 ′ , ω ′ ) and N (1) = G a p (1, 1 ′ , ω) and is thus expressed in angle and action variables by equation ( 27). Since M (2) is as in equation ( 13), the aα coefficients are those of equation ( 28), the species indices being now b and q instead of c and p and the parameters being 2 ′ , ω ′ instead of 1 ′ , ω. This leads to the following change in equation ( 37

): b d2 (F ab (1, 2) • ∇p 1 ) G a p (1, 1 ′ , ω) G b q (2, 2 ′ , ω ′ ) = i Gmamq k 2 k ′ 2 ε -1 (ω ′ ) αβ ψ β * k ′ 2 (J ′ 2 ) e i(k 2 •w 1 -k ′ 2 •w ′ 2 ) ω ′ -k ′ 2 • Ω ′ 2 ψ α k 2 (J1) ik2 • ∇J 1 -(∇J 1 ψ α k 2 (J1)) • ∇w 1 G a p (1, 1 ′ , ω) • (A1)
Note that, as a general rule, operators act on everything on their right, up to the end of the expression or to a closing delimiter.

Using the relation (A1), the collision operator (37) can be written as:

C a (f ) = - p, q ∞ 0 dτ d1 ′ d2 ′ B dω 2π B ′ dω ′ 2π e -i (ω+ω ′ )τ k 2 k ′ 2 i Gmamq (ε -1 (ω ′ )) αβ ψ β * k ′ 2 (J ′ 2 ) e i(k 2 •w 1 -k ′ 2 •w ′ 2 ) ω ′ -k ′ 2 • Ω ′ 2 ψ α k 2 (J1) ik2 • ∇J 1 -(∇J 1 ψ α k 2 (J1)) • ∇w 1 G a p (1, 1 ′ , ω) Fpq(1 ′ , 2 ′ ) • ∇ p ′ 2 -∇ p ′ 1 f p (1 ′ )f q (2 ′ ) • (A2)
The 1-body propagator G a p (1, 1 ′ , ω) is then Fourier-expanded with respect to both angles w1 and w ′ 1 according to equation ( 19) and this expansion is inserted in equation (A2). It then appears that C a (f ) depends on w1 as exp(i(k1 + k2) • w1). Since during relaxation f a (1) remains a function of J1 only, it is possible to average over w1 without loss of information, which brings a Kronecker factor δ(k1 + k2), such that k1 = -k2 = k. The angle-averaged form of the collision operator is:

C a (f ) = - p, q k ∞ 0 dτ d1 ′ d2 ′ B dω 2π B ′ dω ′ 2π e -i (ω+ω ′ )τ k•∇J 1 ψ α -k (1)(ε -1 (ω ′ )) αβ   k ′ 2 Gmamq ψ β * k ′ 2 (2 ′ ) ω ′ -k ′ 2 • Ω ′ 2 e -i k ′ 2 •w ′ 2     k ′ 1 e i k ′ 1 •w ′ 1 G ap kk ′ 1 (1,1 ′ , ω) Fpq(1 ′ ,2 ′ )•(∇ p ′ 2 -∇ p ′ 1 )f p (1 ′ )f q (2 ′ )   • (A3)
The following calculations are somewhat similar to those carried out for an homogeneous plasma by [START_REF] Ichimaru | Basic principles of plasma physics[END_REF]. The expression (A3) can be split into two parts, one, C a 1 (f ), being associated with the operator ∇ p ′ 1 in the last parenthesis and the other, C a 2 (f ), being associated with the operator ∇ p ′ 2 , so that:

C a (f ) = C a 1 (f ) + C a 2 (f ) • (A4) The terms Fpq(1 ′ ,2 ′ ) •∇ p ′ 2 f q (2 ′ ) and Fpq(1 ′ ,2 ′ ) •∇ p ′ 1 f p (1 ′
), multiplied by other functions of 1 ′ and 2 ′ respectively, are subject to an integration over these variables. The structure of these expressions being similar to equation ( 15), they are expressed as in equation ( 18), noting that f p (1 ′ ) and f q (2 ′ ) do not depend on angles. The coefficients aα of the development on the density-potential basis (equation ( 24)) which enter equation ( 18) are calculated from equation ( 26), with appropriate M functions. Integration over angles w ′ 1 or w ′ 2 can then easily be carried out. We are left with:

C a 1 (f ) = -i p q ∞ 0 dτ B dω 2π B ′ dω ′ 2π e -i (ω+ω ′ )τ k (8π 3 G) 2 mampm 2 q (k•∇J 1 ) ψ α -k (1)(ε -1 (ω ′ )) αβ   k ′ 1 d 3 J ′ 1 G ap kk ′ 1 (1,1 ′ , ω)ψ γ -k ′ 1 (1 ′ ) k ′ 1 • ∇ J ′ 1 f p (1 ′ )     k ′ 2 d 3 J ′ 2 ψ β * k ′ 2 (2 ′ )ψ γ * -k ′ 2 (2 ′ ) f q (2 ′ ) ω ′ -k ′ 2 • Ω ′ 2   , (A5)
c 0000 RAS, MNRAS 000, 000-000

C a 2 (f ) = +i p q ∞ 0 dτ B dω 2π B ′ dω ′ 2π e -i (ω+ω ′ )τ k (8π 3 G) 2 mampm 2 q (k•∇J 1 ) ψ α -k (1)(ε -1 (ω ′ )) αβ   k ′ 1 d 3 J ′ 1 G ap kk ′ 1 (1,1 ′ , ω)ψ γ * k ′ 1 (1 ′ )f p (1 ′ )     k ′ 2 d 3 J ′ 2 ψ β * k ′ 2 (2 ′ )ψ γ k ′ 2 (2 ′ ) k ′ 2 •∇ J ′ 2 f q (2 ′ ) ω ′ -k ′ 2 • Ω ′ 2   • (A6)
The integrations over τ and ω ′ which appear in equations (A5) -(A6) are of the general form

h(ω) = ∞ 0 dτ B ′ dω ′ 2π e -i (ω+ω ′ )τ f (ω) g(ω ′ ) • (A7)
The integration over τ is regular, and straightforward, when ω + ω ′ has a negative imaginary part. Otherwise the result must be obtained by analytical continuation. This means that, whatever ω:

h(ω) = B ′ dω ′ 2π -i ω + ω ′ f (ω) g(ω ′ ) • (A8)
The contour B ′ passes above all singularities of g(ω ′ ). For a stable system these are all below or on the real axis. When ω is low enough in the lower half complex plane C -for -ω to be above B ′ , the integration on ω ′ can be carried out by closing the contour B ′ at infinity in the upper half complex plane C + , using the theorem of residues at the unique singularity in the closed up contour, which is at ω ′ = -ω. The closing of B ′ in C + is possible because in the present case (see equations (A5 -A6)) g(ω ′ )/(ω + ω ′ ) decreases at infinity as | ω ′ | -2 , which means that for such ω's, h(ω) = f (ω) g(-ω). Analytical continuation extends this result to other ω's. The two parts of the collision operator then reduce to:

C a 1 (f ) = -i B dω 2π k (8π 3 G) 2 ma(k•∇J 1 ) ψ α -k (1) (ε -1 (-ω)) αβ   p mp k ′ 1 d 3 J ′ 1 G ap kk ′ 1 (1,1 ′ , ω) ψ γ -k ′ 1 (1 ′ ) k ′ 1 •∇ J ′ 1 f p (1 ′ )     q m 2 q k ′ 2 d 3 J ′ 2 ψ β * k ′ 2 (2 ′ )ψ γ * -k ′ 2 (2 ′ )f q (2 ′ ) -ω -k ′ 2 • Ω ′ 2   , (A9) C a 2 (f ) = +i B dω 2π k (8π 3 G) 2 ma(k•∇J 1 ) ψ α -k (1) (ε -1 (-ω)) αβ   p mp k ′ 1 d 3 J ′ 1 G ap kk ′ 1 (1,1 ′ , ω)ψ γ * k ′ 1 (1 ′ )f p (1 ′ )     q m 2 q k ′ 2 d 3 J ′ 2 ψ γ k ′ 2 (2 ′ )ψ β * k ′ 2 (2 ′ ) k ′ 2 •∇ J ′ 2 f q (2 ′ ) -ω -k ′ 2 • Ω ′ 2   • (A10)
The last parenthesis in the second line of equation ( A10) is (δ βγ -ε βγ (-ω))/(8π 3 G). Similarly, the last parenthesis in eq. ( A9) is H βγ (-ω)/8π 3 , where the matrix H αβ (ω) is defined by:

H αβ (ω) = 8π 3 q m 2 q k ′ d 3 J ′ ψ α * k ′ (J ′ ) ψ β * -k ′ (J ′ ) f q (J ′ ) ω -k ′ • Ω(J ′ ) • (A11)
The double Fourier transform G ap kk ′ 1 of the propagator which is present in the first parentheses of equations ( A9) -(A10) can be read from equation (36):

G ap k k ′ 1 (1, 1 ′ , ω) = i (ω -k • Ω1) 1 8π 3 δ a p δ(k+k ′ 1 ) δ(J1 -J ′ 1 ) + Gmpma (k•∇J 1 f a (1)) (ω + k ′ 1 •Ω ′ 1 ) ψ λ k (J1) ε -1 (ω) λµ ψ µ * -k ′ 1 (J ′ 1 ) • (A12)
Using this, the first parenthesis of the second line of equation (A9), V γ a 1 (1, ω), can be written as:

V γ a 1 (1, ω) = - i 8π 3 ma (k •∇J 1 f a (1)) ω -k • Ω1 ψ λ k (1) ε -1 (ω) λγ • (A13)
The first parenthesis of the second line of equation (A10), V γ a 2 (1, ω), is similarly calculated and expressed in terms of the matrix H defined by equation (A11):

V γ a 2 (1, ω) = i 8π 3 maf a (1) ω -k • Ω1 ψ γ * -k (1) + i 8π 3 Gma (k •∇J 1 f a (1)) ω -k • Ω1 ψ λ k (1) ε -1 (ω) λµ H µγ (+ω) • (A14)
Inserting equation (A13) in equation (A9) we get:

c 0000 RAS, MNRAS 000, 000-000

Gravitational Balescu-Lenard equation 17

C a 1 (f ) = - B dω 2π k G 2 m 2 a (k • ∇J 1 ) ψ α -k (1) ε -1 (-ω) αβ ψ λ k (1) ε -1 (+ω) λγ H βγ (-ω) k • ∇J 1 f a (1) ω -k • Ω1 • (A15)
Inserting equation (A14) in equation (A10) we get:

C a 2 (f ) = - B dω 2π k G 2 m 2 a (k • ∇J 1 ) ψ α -k (1) 1 ω -k • Ω1 (ε -1 (-ω)) αγ -δ αγ 1 G f a (1) ψ γ * -k (1) + (k •∇J 1 f a (1)) ψ λ k (1)(ε -1 (+ω)) λµ H µγ (+ω) • (A16)
Gathering equations ( A15) and (A16), the following expression is obtained for the collision operator:

C a (f ) = - B dω 2π k G 2 m 2 a (k•∇J 1 ) ψ α -k (1) 1 ω -k • Ω1 + (ε -1 (-ω)) αβ H βγ (-ω) ψ λ k (1) (ε -1 (+ω)) λγ (k •∇J 1 f a (1)) + (ε -1 (-ω)) αγ ψ λ k (1) (ε -1 (+ω)) λµ H µγ (+ω) (k •∇J 1 f a (1)) + (ε -1 (-ω)) αγ -δ αγ ψ γ * -k (1) 1 G f a (1) -(k •∇J 1 f a (1)) ψ λ k (1) (ε -1 (+ω)) λµ H µα (+ω) • (A17)
The last term on the third line vanishes on integration over ω. Indeed, the Bromwich contour must pass over all singularities of the function f (ω) in eq. (A7), that is, in equation ( A17), over all singularities of functions of +ω. The contour B can be closed at infinity in the upper complex plane, which gives, for the fourth term of the square bracket, a vanishing result. The two terms in the second line of equation (A17) can be associated, yielding the following expression for C a (f ):

C a (f ) = - B dω 2π k 1 G 2 m 2 a (k1 •∇J 1 ) ψ α -k 1 (1) 1 ω -k1 • Ω1 ε -1 (-ω)) αγ -δ αγ ψ γ * -k 1 (1) f a (1) G - B dω 2π k 1
G 2 m 2 a (k1•∇J 1 ) ψ α -k 1 (1) (k1 •∇J 1 f a (1)) ω -k1 • Ω1 (ε -1 (-ω)) αβ ψ λ k (1) (ε -1 (+ω)) λγ H βγ (-ω) + H γβ (+ω)

• (A18)

To evaluate the second line of equation (A18), the integration contour B may be lowered to the real axis. Rigourously, ω pertains the upper complex half plane and can only be consider real in a limit sense when the contour B descends to the real axis. Real singularities at ω = k • Ω must therefore be avoided by the contour by skirting them from above, so that:

1 ω -k • Ω → 1 ω -k • Ω + i0 = P ω -k • Ω -iπδ(ω -k • Ω) , ( A19
)
where P is the principle value distribution. Conversely, when ω descends to the real axis, -ω rises to it from below, so that:

1 -ω + k • Ω → (-1) ω -k • Ω -i0 = - P ω -k • Ω -iπδ(ω -k • Ω) • (A20)
The sum H βγ (-ω) + H γβ (+ω) calculated in this limit is:

H βγ (-ω) + H γβ (+ω) = -16iπ 4 q k ′ 1 d 3 J ′ 1 m 2 q ψ γ * k ′ 1 (1 ′ )ψ β * -k ′ 1 (1 ′ ) δ(ω -k ′ 1 • Ω ′ 1 )f q (1 ′ ) • (A21)
Thanks to the Dirac function in equation (A21), the second line of equation ( A18) is easily integrated over ω. Where conciseness demands, we note:

ω1 = k1•Ω1 , ω ′ 1 = k ′ 1 •Ω ′ 1 •
The first line of equation (A18) can be disposed of by closing the ω integration contour in the lower half complex plane, which is possible because the integrand declines fast enough at infinity. The system being supposedly stable, all the singularities of (ε -1 (-ω)) αγ are in the upper half plane. The contour then encloses only the real singularity at ω = k1 • Ω1 and its sense brings a factor -2iπ when using the residue theorem. The expression of the collision operator then becomes:

C a (f ) = i Gm 2 a k 1 (k1 •∇J 1 ) ψ α
-k 1 (1) (ε -1 (-ω1)) αβ -δ αβ ψ β * -k 1 (1)fa(1)

+ q k ′ 1 d 3 J ′ 1 8π 3 Gm 2 q ψ α -k 1 (1) ε -1 (-ω ′ 1 ) αβ ψ β * -k ′ 1 (1 ′ ) ψ λ k 1 (1) ε -1 (+ω ′ 1 ) λγ ψ γ * k ′ 1 (1 ′ ) k1•∇J 1 (f a (1)f q (1 ′ )) ω ′ 1 -ω1 + i 0 • (A22)

c 0000 RAS

c 0000 RAS, MNRAS 000, 000-000

Acknowledgements

I am indebted to C. Pichon and D. Aubert for awaking my interest for this problem. Their work (Pichon & Aubert 2006) introduced me to the tools which make it possible to extend the methods of plasma physics to these systems. I thank the referee, C. Pichon, for many suggestions which helped to improve the first version of this paper. I also thank the Observatoire Astronomique and the Université de Strasbourg for accepting me as professor emeritus.

This expression is then symmetrized. Half of the term on the second line of eq. ( A22) is added to half of the same expression, modified by changing k1 into -k1 and k ′ 1 into -k ′ 1 . This leaves it invariant, except for the last denominator, which is changed into -(ω ′ 1 -ω1 -i 0). The half difference brings a contribution -iπδ(ω ′ 1 -ω1). The term on the first line of equation (A22) may be similarly symmetrized. When changing k1 to -k1, the argument of the inverse response function changes sign. The change of the response function when the sign of its real frequency argument, ωr say, is modified may be found by noting that its real and imaginary parts, ε ′ αβ and ε" αβ , are:

The conjugation relations (25) can be used to show that:

where the basis element α associated with α is defined by equation ( 23). Using this, the expression T , defined by:

which is present in the first line of equation ( A22) is symmetrized to:

As expected (Nelson & Tremaine 1999), this expression involves the antihermitian part of the matrix ε -1 , which may be expressed in terms of the antihermitian part of the matrix ε as:

The matrix ε † -ε is calculated from equations (A23) -(A24):

The term T in equation (A27) can then be written as:

and inserted in the first line of equation (A22). The square modulus factor in equation (A30 35)). The second line of equation (A22) can be treated similarly. From equation (35), one of the parentheses is

) -1 and the other is its complex conjugate, which can be shown by using the conjugation relation ( 25). When all these symmetrizations and substitutions are made, the collision operator is finally written as: The motion of a particle in a spherically symmetric potential is best described in spherical coordinates r, θ, ϕ, the variable r being the distance to the center, θ the colatitude measured from the pole associated with the coordinate polar axis z and ϕ the azimuth measured from an arbitrarily defined origin in the equatorial plane. Let U (r) be the gravitational potential, an increasing function of r approaching 0 at infinity, which is provisionnally treated as constant in time. The fact that U (r) actually slowly evolves as the relaxation proceeds is addressed in section 8. Without loss of generality, the particle may be assumed to be of unit mass. A dot indicating time derivative, the conjugate momenta to r, θ, ϕ are:

c 0000 RAS, MNRAS 000, 000-000

Gravitational Balescu-Lenard equation 19

Figure B1. Angular parameters associated with the projection M of the particle on the unit sphere. O is the origin of the spherical coordinates, N the origin of the colatitudes and V the ascending node of the orbit. The orbital plane is OVM.

In a constant potential, the energy E of a particle is a first integral, as is the vectorial angular momentum L, that is, its modulus L, its projection Lz on the polar axis and the direction of its projection onto the equatorial plane. The angle and action variables are deduced from the variables r, θ, ϕ, pr, p θ , pϕ by a canonical transformation, the generating function of which is the solution to the Hamilton-Jacobi equation. [START_REF] Goldstein | Classical Mechanics[END_REF] shows how to construct angle and action variables w1, w2, w3, J1, J2, J3 in the case of a newtonian potential. Similar results for a general spherical potential are also well known. They can be found, for example, in Tremaine & [START_REF] Weinberg | [END_REF] or in Saha (1991). A summary is presented in this appendix.

The orbit in a spherically symmetric potential being plane, the periods of the azimuthal and latitudinal motions are equal. This introduces some freedom in defining the actions, which can be taken advantage of to impose that one of the angles, w3 say, be a first integral, associated with the direction of the equatorial projection of the angular momentum. The origin of the constant angle w3 can be chosen to coincide with the origin of the azimuths and the origin of the radial angle variable w1 may be placed at some fiducial periapse. The angle and action variables then are given by the following expressions:

P represents the position of a particle passing at the point M = r, θ, ϕ, with momenta pr, p θ , pϕ when it reaches a fiducial periapse of its orbit. For given position and momenta, the action and angle variables in equations (B2)-(B3) depend on the potential U (r). The radii rP and rA are the distances to the origin of the periapses and the apoapses of the orbit of a particle with actions J = (J1, J2, J3). They are given by the equation:

and they depend on E and J2, that is, on J1 and J2 but not on J3. There being many different periapses, the fiducial one must be defined not only by its spatial location, but also by the time at which the particle passes there. The sign ± in equations (B3) should be taken as + when the particle visits the fiducial periapse P before it passes at M andotherwise. Ω1 and Ω2 are the pulsations of the radial and latitudinal motions respectively (equations (B5)). The angle ψ is the azimuth of the present particle's position in the orbital plane, measured from the ascending node (figure B1). The angle w2, which varies linearly in time, is the mean angular motion of the particle in the plane of its orbit. The constant angle w3 is the azimuth of the ascending node in the equatorial plane. The angles w1 and w2 are expressed as radial integrals following the sense of the particle's motion, whence the presence of an absolute value of the differential element in equations (B3). The boundaries of integrals on r ′ have not been written as rP and r because, depending on the relative position of the particle and the fiducial periapse, the integral may be extended over several successive senses of the radial oscillations. The ratio J3/J2 is the cosine of an inclination angle β (figure B1) defined by cos β = J3/J2. The latitude of the particle oscillates between ±β. The frequency Ω3 vanishes and the frequencies Ω1 and Ω2 are given by:

They are both positive. The angle variables w1 and w2 then increase linearly with time with the frequencies Ω1 and Ω2, changing by 2π in a complete, respectively radial and latitudinal, oscillation. Equations ( B2)-(B3) give the angle and action variables in terms of the position and momentum variables. These relations may be inverted to give the latter in terms of the former. This however involves the inversion of the implicit relation (B2) to obtain E as a function of J1 and J2 and of the first of equations (B3) to obtain r as a function of w1, J1, J2.

B2 Basis Fourier coefficients for spherical potentials

The basis expansion coefficients which correspond to the density distribution (58) are obtained from equations ( 17) and ( 26). When, as here, the distribution functions do not depend on angles, their integrals over angles in equation ( 26) are simply proportional to the k = 0 Fourier coefficient of ψ α * , or equivalently of ψ α (equation ( 23)), which is the complex conjugate of ψ α 0 (J, t) (equation ( 25)). This coefficient depends on time, due to the slow variation of the orbits. We find that:

From equations ( 21) and ( 24), the gravitational potential is U (r, t) = G aα(t) ψ α (r), α being a index. Its partial time derivative is given by equation ( 59). We also need some explicit expression for the angle Fourier coefficients ψ α k (J, t) We also need some explicit expression for the angle Fourier coefficients ψ α k (J, t) of the basis potentials (equation ( 60)). The relation of the position r to the angle and action variables w and J depends on the potential U (r, t), and thus on t. One could think of evaluating ψ α k (J, t) for a given potential U (r, t) by just calculating the integral over angles in equation ( 60). The position vector r would then have to be expressed in terms of the angle vector w, for given actions. This cannot be done explicitly in general, since the relations (B2)-(B3) would have to be inverted. For spherical potentials, it is easier to change the variables of integration w1, w2, w3 in equation ( 60) for position-type variables r, ψ, w3 (figure B1). For given actions J, these variables are related by the equations (B3) which can also be written, with the notations explained above, as w1 = W1(J, M (r), t) and w2 = ψ + W2(J, M (r), t), where:

The jacobian of the transformation from w1, w2, w3 to r, ψ, w3 is | dW1/dr |. For a spherical potential, the basis potential functions can be chosen to depend only on the radial distance r. Equation (60) then becomes:

The integrations over the angles ψ and w3 reduce to 4π 2 δ 0 k 2 δ 0 k 3 where the δ's are Kronecker symbols. Thus the coefficients ψ α k differ from zero only when the k2 and k3 components vanish. The radial integration is over a complete oscillatory cycle of the variable r, from rP to rA and back. The coefficients ψ α k (J, t) depend on the potential U (r, t) and on the k1 component of k, hereafter simply noted k. Equation (B8) then reduces to:

| dr | ψ α (r) e -ikW 1 (M (r),t)

2(E -U (r, t)) -J 2 2 /r 2 • (B9)

The cycle integral over r in equation (B9) can be separated into an ascending part, in which r increases from rP to rA, and a descending part in which it decreases from rA to rP . Let W + 1 (r) be the value of W1(M (r), t) during the ascending part and W - 1 (r) its value during the descending part. W1(M (r), t) is a monotonically increasing function along the oscillation. Its value at the apoapse is π. Equation (B7) shows that π -W + 1 (r) = W - 1 (r) -π. Defining W k (r, t) = kW + 1 (r, t), equation (B9) is turned into equation (61).

Accepted by MNRAS, 2010, April 22 ;Received 2010 April 13 ;in original form 2010 February 25