# Kinetic theory of spatially inhomogeneous stellar systems without collective effects

###### Key Words.:

Gravitation – Methods: analytical – Galaxies: star clusters: generalWe review and complete the kinetic theory of spatially inhomogeneous stellar systems when collective effects (dressing of the stars by their polarization cloud) are neglected. We start from the BBGKY hierarchy issued from the Liouville equation and consider an expansion in powers of in a proper thermodynamic limit. For , we obtain the Vlasov equation describing the evolution of collisionless stellar systems like elliptical galaxies. At the order , we obtain a kinetic equation describing the evolution of collisional stellar systems like globular clusters. This equation coincides with the generalized Landau equation derived by Kandrup (1981) from a more abstract projection operator formalism. This equation does not suffer logarithmic divergences at large scales since spatial inhomogeneity is explicitly taken into account. Making a local approximation, and introducing an upper cut-off at the Jeans length, it reduces to the Vlasov-Landau equation which is the standard kinetic equation of stellar systems. Our approach provides a simple and pedagogical derivation of these important equations from the BBGKY hierarchy which is more rigorous for systems with long-range interactions than the two-body encounters theory. Making an adiabatic approximation, we write the generalized Landau equation in angle-action variables and obtain a Landau-type kinetic equation that is valid for fully inhomogeneous stellar systems and is free of divergences at large scales. This equation is less general than the Lenard-Balescu-type kinetic equation recently derived by Heyvaerts (2010) since it neglects collective effects, but it is substantially simpler and could be useful as a first step. We discuss the evolution of the system as a whole and the relaxation of a test star in a bath of field stars. We derive the corresponding Fokker-Planck equation in angle-action variables and provide expressions for the diffusion coefficient and friction force.

## 1 Introduction

In its simplest description, a stellar system can be viewed as a collection of
classical point mass stars in Newtonian gravitational interaction (Binney and Tremaine 2008).
As understood early by Hénon (1964), self-gravitating systems experience
two successive types of relaxation: A rapid “collisionless” relaxation towards a
quasi stationary state (QSS) that is a virialized state in mechanical equilibrium
but not in thermodynamical equilibrium, followed by a slow “collisional”
relaxation. One might think that, due to the development of stellar encounters, the system will reach, for sufficiently large times, a statistical equilibrium state described by the Maxwell-Boltzmann distribution. However, it is well-known that unbounded stellar systems cannot be in strict statistical equilibrium^{1}

The first kinetic equation was written by Jeans (1915). Neglecting encounters between stars, he
described the dynamical evolution of stellar systems by the
collisionless Boltzmann equation coupled to the Poisson equation. This purely mean field description applies to large groups of stars such as elliptical galaxies whose ages are much less than the collisional relaxation time. A
similar equation was introduced by Vlasov (1938) in plasma physics to describe the collisionless evolution of a system of electric charges interacting by the Coulomb force. The collisionless Boltzmann equation coupled self-consistently to the Poisson equation is oftentimes called the Vlasov equation^{2}

The concept of collisionless relaxation was first understood by Hénon (1964) and King (1966). Lynden-Bell (1967) developed a statistical theory of this process and coined the term “violent relaxation”. In the collisionless regime, the system is described by the Vlasov-Poisson system. Starting from an unsteady or unstable initial condition, the Vlasov-Poisson system develops a complicated mixing process in phase space. The Vlasov equation, being time-reversible, never achieves a steady state but develops filaments at smaller and smaller scales. However, the coarse-grained distribution function usually achieves a steady state on a few dynamical times. Lynden-Bell (1967) tried to predict the QSS resulting from violent relaxation by developing a statistical mechanics of the Vlasov equation. He derived a distribution function formally equivalent to the Fermi-Dirac distribution (or to a superposition of Fermi-Dirac distributions). However, when coupled to the Poisson equation, these distributions have an infinite mass. Therefore, the Vlasov-Poisson system has no statistical equilibrium state (in the sense of Lynden-Bell). This is a clear evidence that violent relaxation is incomplete (Lynden-Bell 1967). Incomplete relaxation is due to inefficient mixing and non-ergodicity. In general, the fluctuations of the gravitational potential that drive the collisionless relaxation last only for a few dynamical times and die out before the system has mixed efficiently (Tremaine et al. 1986). Understanding the origin of incomplete relaxation, and developing models of incomplete violent relaxation to predict the structure of galaxies, is a very difficult problem (Arad & Johansson 2005). Some models of incomplete violent relaxation have been proposed based on different physical arguments (Bertin and Stiavelli 1984, Stiavelli and Bertin 1987, Hjorth and Madsen 1991, Chavanis 1998, Levin et al. 2008).

On longer timescales, stellar encounters
(sometimes referred to as “collisions” by an abuse of language) must
be taken into account. This description is particularly important for small groups of stars such as globular clusters whose ages are of the order of the collisional relaxation time. Chandrasekhar (1942,1943a,1943b) developed a kinetic
theory of stellar systems in order to determine the timescale of
collisional relaxation and the rate of escape of stars from globular
clusters^{3}^{4}^{5}

King (1960) noted that, if we were to describe the dynamical evolution of the cluster as a whole, the distribution of the field stars must evolve in time in a self-consistent manner so that the kinetic equation must be an integrodifferential equation. The kinetic equation obtained by King, from the results of Rosenbluth et al. (1957), is equivalent to the Landau equation, although written in a different form. It is interesting to note, for historical reasons, that none of the previous authors seemed to be aware of the work of Landau (1936) in plasma physics. There is, however, an important difference between stellar dynamics and plasma physics. Neutral plasmas are spatially homogeneous due to Debye shielding. By contrast, stellar systems are spatially inhomogeneous. The above-mentioned kinetic theories developed for an infinite homogeneous system can be applied to an inhomogeneous system only if we make a local approximation. In that case, the collision term is calculated as if the system were spatially homogeneous or as if the collisions could be treated as local. Then, the effect of spatial inhomogeneity is only retained in the advection (Vlasov) term which describes the evolution of the system due to mean-field effects. This leads to the Vlasov-Landau equation which is the standard kinetic equation of stellar dynamics. To our knowledge, this equation has been first written (in a different form), and studied, by Hénon (1961). Hénon also exploited the timescale separation between the dynamical time and the relaxation time to derive a simplified kinetic equation for , where is the individual energy of a star by unit of mass, called the orbit-averaged Fokker-Planck equation. In this approach, the distribution function , averaged over a short timescale, is a steady state of the Vlasov equation of the form which slowly evolves in time, on a long timescale, due to the development of “collisions” (i.e. correlations caused by finite effects or graininess). Hénon used this equation to obtain a more relevant value for the rate of evaporation from globular clusters, valid for inhomogeneous systems. Cohn (1980) numerically solved the orbit-averaged Fokker-Planck equation to describe the collisional evolution of star clusters. His treatment accounts both for the escape of high energy stars put forward by Spitzer (1940) and for the phenomenon of core collapse resulting from the gravothermal catastrophe discovered by Antonov (1962) and Lynden-Bell and Wood (1968) on the basis of thermodynamics and statistical mechanics. The local approximation, which is a crucial step in the kinetic theory, is supported by the stochastic approach of Chandrasekhar and von Neumann (1942) showing the preponderance of the nearest neighbor. However, this remains a simplifying assumption which is not easily controllable. In particular, as we have already indicated, the local approximation leads to a logarithmic divergence at large scales that is difficult to remove. This divergence would not have occurred if full account of spatial inhomogeneity had been given since the start.

The effect of spatial inhomogeneity was investigated by Severne
and Haggerty (1976), Parisot and Severne (1979), Kandrup
(1981), and Chavanis (2008). In particular, Kandrup (1981) derived
a generalized Landau equation from the Liouville equation by using
projection operator technics. This generalized Landau equation is interesting
because it takes into account effects of spatial inhomogeneity
which were neglected in previous approaches. Since the finite
extension of the system is properly accounted for, there is no
divergence at large scales^{6}

In addition, this equation is still approximate as it neglects
collective effects and considers binary collisions between “naked”
particles. As in any weakly coupled system, the particles engaged in
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. In the case of plasmas, collective effects are responsible for
Debye shielding and they are accounted for in the Lenard-Balescu
equation. They allow to eliminate the logarithmic divergence that
occurs at large scales in the Landau equation. For self-gravitating
systems, they lead to “anti-shielding” and are more difficult to
analyze because the system is spatially inhomogeneous^{7}

Recently, Heyvaerts (2010) derived from the BBGKY hierarchy issued from the Liouville equation a kinetic equation in angle-action variables that takes both spatial inhomogeneities and collective effects into account. To calculate the collective response, he used Fourier-Laplace transforms and introduced a bi-orthogonal basis of pairs of density-potential functions (Kalnajs 1971a). The kinetic equation derived by Heyvaerts is the counterpart for spatially inhomogeneous self-gravitating systems of the Lenard-Balescu equation for plasmas. Following his work, we showed that this equation could be obtained equivalently from the Klimontovich equation by making the so-called quasilinear approximation (Chavanis 2012). We also developed a test particle approach and derived the corresponding Fokker-Planck equation in angle-action variables, taking collective effects into account. This provides general expressions of the diffusion coefficient and friction force for spatially inhomogeneous stellar systems.

In a sense, these equations solve the problem of the kinetic theory of
stellar systems since they take into account both spatial
inhomogeneity and collective effects. However, their drawback is that
they are extremely complicated to solve (in addition of being
complicated to derive). In an attempt to reduce the complexity of the
problem, we shall derive in this paper a kinetic equation that is
valid for spatially inhomogeneous stellar systems but that neglects
collective effects. Collective effects may be less crucial in stellar
dynamics than in plasma physics. In plasma physics, they must be taken
into account in order to remove the divergence at large scales that
appears in the Landau equation. In the case of stellar systems, this
divergence is removed by the spatial inhomogeneity of the system, not
by collective effects. Actually, previous kinetic equations based on
the local approximation ignore collective effects and already give
satisfactory results. We shall therefore ignore collective effects and
derive a kinetic equation (in position-velocity and angle-action
variables) that is the counterpart for spatially inhomogeneous
self-gravitating systems of the Landau equation for plasmas. Our
approach has three main interests. First, the derivation of this
Landau-type kinetic equation is considerably simpler than the
derivation of the Lenard-Balescu-type kinetic equation, and it can be
done in the physical space without having to introduce Laplace-Fourier
transforms nor bi-orthogonal basis of pairs of density-potential
functions. This offers a more physical derivation of kinetic equations
of stellar systems that can be of interest for
astrophysicists. Secondly, our approach is sufficient to remove the
large-scale divergence that occurs in kinetic theories based on the
local approximation. It represents therefore a conceptual progress in
the kinetic theory of stellar systems. Finally, this equation is
simpler than the Lenard-Balescu-type kinetic equation derived by
Heyvaerts (2010), and it could be useful in a first step before
considering more complicated effects. Its drawback is to ignore
collective effects but this may not be crucial as we have explained
(as suggested by Weinberg’s work, collective effects become important
only when the system is close to instability). This Landau-type
equation was previously derived for systems with arbitrary long-range
interactions^{8}

The paper is organized as follows. In Section 2, we study the dynamical evolution of a spatially inhomogeneous stellar system as a whole. Starting from the BBGKY hierarchy issued from the Liouville equation, and neglecting collective effects, we derive a general kinetic equation valid at the order in a proper thermodynamic limit. For , it reduces to the Vlasov equation. At the order we recover the generalized Landau equation derived by Kandrup (1981) from a more abstract projection operator formalism. This equation is free of divergence at large scales since spatial inhomogeneity has been properly accounted for. Making a local approximation and introducing an upper cut-off at the Jeans length, we recover the standard Vlasov-Landau equation which is usually derived from a kinetic theory based on two-body encounters. Our approach provides an alternative derivation of this fundamental equation from the more rigorous Liouville equation. It has therefore some pedagogical interest. In Section 3, we study the relaxation of a test star in a steady distribution of field stars. We derive the corresponding Fokker-Planck equation and determine the expressions of the diffusion and friction coefficients. We emphasize the difference between the friction by polarization and the total friction (this difference may have been overlooked in previous works). For a thermal bath, we derive the Einstein relation between the diffusion and friction coefficients and obtain the explicit expression of the diffusion tensor. This returns the standard results obtained from the two-body encounters theory but, again, our presentation is different and offers an alternative derivation of these important results. For that reason, we give a short review of the basic formulae. In Section 4, we derive a Landau-type kinetic equation written in angle-action variables and discuss its main properties. This equation, which does not make the local approximation, applies to fully inhomogeneous stellar systems and is free of divergence at large scales. We also develop a test particle approach and derive the corresponding Fokker-Planck equation in angle-action variables. Explicit expressions are given for the diffusion tensor and friction force, and they are compared with previous expressions obtained in the literature.

## 2 Evolution of the system as a whole

### 2.1 The BBGKY hierarchy

We consider an isolated system of stars with identical mass in Newtonian interaction. Their dynamics is fully described by the Hamilton equations

(1) |

This Hamiltonian system conserves the energy , the mass , and the angular momentum . As recalled in the Introduction, stellar systems cannot reach a statistical equilibrium state in a strict sense. In order to understand their evolution, it is necessary to develop a kinetic theory.

We introduce the -body distribution function giving the probability density of finding, at time , the first star with position and velocity , the second star with position and velocity etc. Basically, the evolution of the -body distribution function is governed by the Liouville equation

(2) |

where

(3) |

is the gravitational force by unit of mass experienced by the -th star due to its interaction with the other stars. Here, denotes the exact gravitational potential produced by the discrete distribution of stars and denotes the exact force by unit of mass created by the -th star on the -th star. The Liouville equation (2), which is equivalent to the Hamilton equations (2.1), contains too much information to be exploitable. In practice, we are only interested in the evolution of the one-body distribution .

From the Liouville equation (2) we can construct the complete BBGKY hierarchy for the reduced distribution functions

(4) |

where the notation stands for . The generic term of this hierarchy reads

(5) |

This hierarchy of equations is not closed since the equation for the one-body distribution involves the two-body distribution , the equation for the two-body distribution involves the three-body distribution , and so on.

It is convenient to introduce a cluster representation of the distribution functions. Specifically, we can express the reduced distribution in terms of products of distribution functions of lower order plus a correlation function [see, e.g., Eqs. (8) and (2.2) below]. Considering the scaling of the terms in each equation of the BBGKY
hierarchy, we can see that there exist solutions of the whole BBGKY hierarchy such that the
correlation functions scale as in the proper thermodynamic limit defined in Appendix A. This implicitly assumes that the initial condition has
no correlation, or that the initial correlations respect this scaling^{9}

### 2.2 The truncation of the BBGKY hierarchy at the order

The first two equations of the BBGKY hierarchy are

(6) | |||||

(7) |

We decompose the two- and three-body distributions in the form

(8) |

where is the correlation function of order . Substituting Eqs. (8) and (2.2) in Eqs. (6) and (7), and simplifying some terms, we obtain

(10) | |||||

Equations (10) and (LABEL:trunc6) are exact for all but they are not closed. As explained previously, we shall close these equations at the order in the thermodynamic limit . In this limit and . On the other hand, and (see Appendix A).

The term in the l.h.s. of Eq. (10) is of order , and the term in the r.h.s. is of order . Let us now consider the terms in Eq. (LABEL:trunc6) one by one. The first four terms correspond to the Liouville equation. The Liouvillian describes the complete two-body problem, including the inertial motion, the interaction between the particles and the mean field produced by the other particles. The terms and are of order while the term is of order . Therefore, the interaction term can a priori^{10}

If we introduce the notations (distribution function) and (two-body correlation function), we get at the order :

(12) |

(13) |

We have introduced the mean force (by unit of mass) created on star by all the other stars

(14) |

and the fluctuating force (by unit of mass) created by star on star :

(15) |

Equations (12) and (2.2) are exact at the order . They form the right basis to develop the kinetic theory of stellar systems at this order of approximation. Since the collision term in the r.h.s. of Eq. (12) is of order , we expect that the relaxation time of stellar systems scales as where is the dynamical time. As we shall see, the discussion is more complicated due to the presence of logarithmic corrections in the relaxation time and the absence of strict statistical equilibrium state.

### 2.3 The limit : the Vlasov equation (collisionless regime)

In the limit , for a fixed interval of time (any), the correlations between stars can be neglected. Therefore, the mean field approximation becomes exact and the -body distribution function factorizes in one-body distribution functions

(16) |

Substituting the factorization (16) in the Liouville equation (2), and integrating over , , …, , we find that the smooth distribution function is the solution of the Vlasov equation

(17) |

This equation also results from Eq. (12) if we neglect the correlation function in the r.h.s. and replace by .

The Vlasov equation describes the collisionless evolution of stellar systems for times shorter than the relaxation time . In practice so that the domain of validity of the Vlasov equation is huge (see the end of Sec. 2.8). As recalled in the Introduction, the Vlasov-Poisson system develops a process of phase mixing and violent relaxation leading to a quasi-stationary state (QSS) on a very short timescale, of the order of a few dynamical times . Elliptical galaxies are in such QSSs. Lynden-Bell (1967) developed a statistical mechanics of the Vlasov equation in order to describe this process of violent relaxation and predict the QSS achieved by the system. Unfortunately, the predictions of his statistical theory are limited by the problem of incomplete relaxation. Kinetic theories of violent relaxation, which may account for incomplete relaxation, have been developed by Kadomtsev and Pogutse (1970), Severne and Luwel (1980), and Chavanis (1998,2008).

### 2.4 The order : the generalized Landau equation (collisional regime)

If we neglect strong collisions and collective effects, the first two equations (12) and (2.2) of the BBGKY hierarchy reduce to

(18) | |||||

(19) |

The first equation gives the evolution of the one-body distribution function. The l.h.s. corresponds to the (Vlasov) advection term. The r.h.s. takes into account correlations (finite effects, graininess, discreteness effects) between stars that develop due to their interactions. These correlations correspond to “collisions”.

Equation (2.4) may be viewed as a linear first order differential equation in time. It can be symbolically written as

(20) |

where is a mean field Liouvillian and is a source term expressible in terms of the one-body distribution. This equation can be solved by the method of characteristics. Introducing the Green function

(21) |

constructed with the mean field Liouvillian , we obtain

(22) |

where we have assumed that no correlation is present initially so that (if correlations are present initially, it can be shown that they are rapidly washed out). Substituting Eq. (2.4) in Eq. (18), we obtain

(23) |

In writing this equation, we have adopted a Lagrangian point of view. The coordinates appearing after the Greenian must be viewed as and . Therefore, in order to evaluate the integral (23), we must move the stars following the trajectories determined by the self-consistent mean field.

The kinetic equation (23) is valid at the order
so it describes the “collisional” evolution of the system
(ignoring collective effects) on a timescale of order . Equation (23) is a
non-Markovian integro-differential equation. It takes into account delocalizations in
space and time (i.e. spatial inhomogeneity and memory effects). Actually,
the Markovian approximation is justified
in the limit because the timescale
over which changes is long compared to the correlation
time over which the integrand in Eq. (23) has
significant support^{11}

(24) |

Similarly, we can compute the trajectories of the stars by assuming that the mean field is independent on and equal to its value at time so that and .

The structure of the kinetic equation (24) has a clear physical meaning. The l.h.s. corresponds to the Vlasov advection term due to mean field effects. The r.h.s. can be viewed as a collision operator taking finite effects into account. For , it vanishes and we recover the Vlasov equation. For finite , it describes the cumulative effect of binary collisions between stars. The collision operator is a sum of two terms: A diffusion term and a friction term. The coefficients of diffusion and friction are given by generalized Kubo formulae, i.e. they involve the time integral of the auto-correlation function of the fluctuating force. The kinetic equation (24) bears some resemblance with the Fokker-Planck equation. However, it is more complicated since it is an integro-differential equation, not a differential equation (see Section 3).

Equation (24) may be viewed as a generalized Landau equation. Since the spatial inhomogeneity and the finite extension of the system are properly taken into account, there is no divergence at large scales. There is, however, a logarithmic divergence at small scales which is due to the neglect of the interaction term in the Liouvillian (see Sec. 2.6). At large scales (i.e. for large impact parameters), this term can be neglected and the trajectories of the stars are essentially due to the mean field. Thus, we can make the weak coupling approximation leading to the Landau equation. This approximation describes weak collisions. However, at small scales (i.e. for small impact parameters), we cannot ignore the interaction term in the Liouvillian anymore and we have to solve the classical two-body problem. This is necessary to correctly describe strong collisions for which the trajectory of the particles deviates strongly from the mean field motion. When the mean field Greenian (constructed with ) is replaced by the total Greenian (constructed with ), taking into account the interaction term, the generalized Landau equation (24) is replaced by a more complicated equation which can be viewed as a generalized Boltzmann equation. This equation is free of divergence (since it takes both spatial inhomogeneity and strong collisions into account) but it is unnecessarily complicated because it does not exploit the dominance of weak collisions over (rare) strong collisions for the gravitational potential. Indeed, a star suffers a large number of weak distant encounters and very few close encounters. A better practical procedure is to use the generalized Landau equation (24) with a cut-off at small scales in order to take into account our inability to describe strong collisions by this approach.

The generalized Landau equation (24) was derived by Kandrup (1981) from the Liouville equation by using the projection operator formalism. It can also be derived from the Liouville equation by using the BBGKY hierarchy or from the Klimontovich equation by making a quasilinear approximation (Chavanis 2008).

### 2.5 The Vlasov-Landau equation

Self-gravitating systems are intrinsically spatially inhomogeneous. However, the collision operator at position can be simplified by making a local approximation and performing the integrals as if the system were spatially homogeneous with the density . This amounts to replacing by and by in Eq. (24). This local approximation is motivated by the work of Chandrasekhar and von Neumann (1942) who showed that the distribution of the gravitational force is a Lévy law (called the Holtzmark distribution) dominated by the contribution of the nearest neighbor. With this local approximation, Eq. (24) becomes

(25) |

where we have used . The Greenian corresponds to the free motion of the particles associated with the Liouvillian . Using Eqs. (124) and (125), the foregoing equation can be rewritten as

(26) |

where is expressed in terms of the Lagrangian coordinates. The integrals over and can be calculated explicitly (see Appendix C). We then find that the evolution of the distribution function is governed by the Vlasov-Landau equation

(27) |

where we have noted , , , and where represents the Fourier transform of the gravitational potential. Under this form, we see that the collisional evolution of a stellar system is due to a condition of resonance (with ) encapsulated in the -function. This -function expresses the conservation of energy.

The Vlasov-Landau equation can also be written as (see Appendix C):

(28) | |||||

(29) |

where

(30) |

is the Coulomb factor that has to be regularized with appropriate
cut-offs (see Section 2.6). The r.h.s. of Eq. (28) is
the original form of the collision operator given by Landau (1936) for
the Coulombian interaction^{12}^{13}

### 2.6 Heuristic regularization of the divergence

To obtain the Vlasov-Landau equation (28), we have made a local approximation. This amounts to calculating the collision operator at each point as if the system were spatially homogeneous. As a result of this homogeneity assumption, a logarithmic divergence appears at large scales in the Coulombian factor (30). In plasma physics, this divergence is cured by the Debye shielding. A charge is surrounded by a polarization cloud of opposite charges which reduces the range of the interaction. When collective effects are properly taken into account, as in the Lenard-Balescu equation, no divergence appears as large scales and the Debye length arises naturally. Heuristically, we can use the Landau equation and introduce an upper cut-off at the Debye length . For self-gravitating systems, there is no shielding and the divergence is cured by the finite extent of the system. The interaction between two stars is only limited by the size of the system. When spatial inhomogeneity is taken into account, as in the generalized Landau equation (24), no divergence occurs at large scales. Heuristically, we can use the Vlasov-Landau equation (28) and introduce an upper cut-off at the Jeans length which gives an estimate of the system size .

The Coulombian factor (30) also diverges at small scales. As explained previously, this is due to the neglect of strong collisions that produce important deflections. Indeed, for collisions with low impact parameter, the mean field approximation is clearly irrelevant and it is necessary to solve the two-body problem exactly. Heuristically, we can use the Landau equation and introduce a lower cut-off at the gravitational Landau length (the gravitational analogue of the Landau length in plasma physics) which corresponds to the impact parameter leading to a deflection at .

Introducing a large-scale cut-off at the Jeans length and a small-scale cut-off at the Landau length , and noting that , we find that the Coulombian factor scales as where is the number of stars in the
cluster^{14}

### 2.7 Properties of the Vlasov-Landau equation

The Vlasov-Landau equation conserves the mass and the energy . It also monotonically
increases the Boltzmann entropy in the sense that (-theorem). Due to the local approximation, the proof of these properties is the same^{15}

(31) |

where is the energy and the angular momentum. This distribution takes into account the escape of high energy stars and the anisotropy of the velocity distribution. It can be derived, under some approximations, from the Vlasov-Landau equation by using the condition that if the energy of the star is larger than the escape energy (Michie 1963, King 1965). The Michie-King distribution reduces to the isothermal distribution for low energies. In this sense, we can define a “relaxation” time for a stellar system. From the Vlasov-Landau equation (28), we find that the relaxation time scales as

(32) |

where is the root mean square (r.m.s.) velocity. Introducing the dynamical time , we obtain the scaling

(33) |

The fact that the ratio between the relaxation time and the dynamical time depends only on the number of stars and scales as was first noted by Chandrasekhar (1942).

A simple estimate of the evaporation time gives (Ambartsumian 1938, Spitzer 1940). More precise values have been obtained by studying the evaporation process in an artificially uniform cluster (Chandrasekhar 1943b, Spitzer and Härm 1958, Michie 1963, King 1965, Lemou and Chavanis 2010) or in a more realistic inhomogeneous cluster (Hénon 1961). Since , we can consider that the system relaxes towards a steady distribution of the form (31) on a timescale and that this distribution slowly evolves on a longer timescale as the stars escape^{16}

Even if the system were confined within a small box so as to prevent both the evaporation and the gravothermal catastrophe, there would be no statistical equilibrium state in a strict sense because there is no global entropy maximum (Antonov 1962). A configuration in which some subset of the particles are tightly bound together (e.g. a binary star), and in which the rest of the particles shares the energy thereby released, may have an arbitrary large entropy. However, such configurations, which require strong correlations, are generally reached very slowly (on a timescale much larger than ) due to encounters involving many particles. To describe these configurations, one would have to take high order correlations into account in the kinetic theory. These configurations may be relevant in systems with a small number of stars (Chabanol et al. 2000) but when is large the picture is different. On a timescale of the order of the one-body distribution function is expected to reach the Boltzmann distribution which is a local entropy maximum. This state is “metastable” but its lifetime is expected to be very large, scaling as , so that it is stable in practice (Chavanis 2006). In this sense, there exist “true” statistical equilibrium states for self-gravitating systems confined within a small box. However, we may argue that this situation is highly artificial.

### 2.8 Dynamical evolution of stellar systems: a short review

Using the kinetic theory, we can identify different phases in the dynamical evolution of stellar systems.

A self-gravitating system initially out-of-mechanical equilibrium undergoes a process of violent collisionless relaxation towards a virialized state. In this regime, the dynamical evolution of the cluster is described by the Vlasov-Poisson system. The phenomenology of violent relaxation has been described by Hénon (1964), King (1966), and Lynden-Bell (1967). Numerical simulations that start from cold and clumpy initial conditions generate a quasi stationary state (QSS) that fits the de Vaucouleurs law for the surface brightness of elliptical galaxies quite well (van Albada 1982). The inner core is almost isothermal (as predicted by Lynden-Bell 1967) while the velocity distribution in the envelope is radially anisotropic and the density profile decreases as . One success of Lynden-Bell’s statistical theory of violent relaxation is to explain the isothermal core of elliptical galaxies without recourse to “collisions”. By contrast, the structure of the halo cannot be explained by Lynden-Bell’s theory as it results from an incomplete relaxation. Models of incompletely relaxed stellar systems have been elaborated by Bertin and Stiavelli (1984), Stiavelli and Bertin (1987), and Hjorth and Madsen (1991). These theoretical models nicely reproduce the results of observations or numerical simulations (Londrillo et al. 1991, Trenti et al. 2005). In the simulations, the initial condition needs to be sufficiently clumpy and cold to generate enough mixing required for a successful application of the statistical theory of violent relaxation. Numerical simulations starting from homogeneous spheres (see, e.g., Roy and Perez 2004, Levin et al. 2008, Joyce et al. 2009) show little angular momentum mixing and lead to different results. In particular, they display a larger amount of mass loss (evaporation) than simulations starting from clumpy initial conditions. Clumps thus help the system to reach a “universal” final state from a variety of initial conditions, which can explain the similarity of the density profiles observed in elliptical galaxies.

On longer timescales, encounters between stars must be taken into account and the dynamical evolution of the cluster is governed by the Vlasov-Landau-Poisson system.
The first stage of the collisional evolution is driven by
evaporation. Due to a series of weak encounters, the energy of a
star can gradually increase until it reaches the local escape energy;
in that case, the star leaves the system.
Numerical simulations (Spitzer 1987, Binney and Tremaine 2008) show that
during this regime the system reaches a quasi-stationary state that
slowly evolves in amplitude due to evaporation as the system loses
mass and energy. This quasi stationary distribution function is
close to the Michie-King model (31). The system has a core-halo
structure. The core is isothermal while the stars in the outer halo
move in predominantly radial orbits. Therefore, the distribution
function in the halo is anisotropic. The density follows the
isothermal law in the central region (with a core of
almost uniform density) and decreases as in the
halo (for an isolated cluster with ). Due to evaporation, the halo expands while the core shrinks as required by energy conservation. During the evaporation process, the
central density increases permanently. At some point of the
evolution, the system undergoes an
instability related to the Antonov (1962) instability^{17}

At the present epoch, small groups of stars such as globular clusters (, , age , ) are in the collisional regime. They are either in quasistationary states described by the Michie-King model or experiencing core collapse. By contrast, large clusters of stars like elliptical galaxies (, , age , ) are still in the collisionless regime and their apparent organization is a result of an incomplete violent relaxation.

## 3 Test star in a thermal bath

### 3.1 The Fokker-Planck equation

We now consider the relaxation of a “test” star (tagged particle)
evolving in a steady distribution of “field” stars^{18}^{19}

(34) |

where is the static mean force created by the field stars with density . This equation does not present any divergence at large scales.

If we make a local approximation and use the Vlasov-Landau equation (26), we obtain

(35) |

Denoting the advection operator by , Eq. (35) can be written in the form of a Fokker-Planck equation

(36) |

involving a diffusion tensor

(37) |

and a friction force

(38) |

If we directly start from Eq. (27), which amounts to performing the integrals over and in the previous expressions, we obtain the Fokker-Planck equation

(39) |

with the diffusion and friction coefficients