A The thermodynamic limit

Kinetic theory of spatially inhomogeneous stellar systems without collective effects

Key Words.:
Gravitation – Methods: analytical – Galaxies: star clusters: general

We 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 equilibrium1 due to the permanent escape of high energy stars (evaporation) and the gravothermal catastrophe (core collapse). Therefore, the statistical mechanics of stellar systems is essentially an out-of-equilibrium problem which must be approached through kinetic theories.

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 equation2, or the Vlasov-Poisson system.

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 clusters3. To simplify the kinetic theory, he considered an infinite homogeneous system of stars. He started from the general Fokker-Planck equation and determined the diffusion coefficient and the friction force (first and second moments of the velocity increments) by considering the mean effect of a succession of two-body encounters4. However, his approach leads to a logarithmic divergence at large scales that is more difficult to remove in stellar dynamics than in plasma physics because of the absence of Debye shielding for the gravitational force5. Chandrasekhar and von Neumann (1942) developed a completely stochastic formalism of gravitational fluctuations and showed that the fluctuations of the gravitational force are given by the Holtzmark distribution (a particular Lévy law) in which the nearest neighbor plays a prevalent role. From these results, they argued that the logarithmic divergence has to be cut-off at the interparticle distance (see also Jeans 1929). However, since the interparticle distance is smaller than the Debye length, the same arguments should also apply in plasma physics, which is not the case. Therefore, the conclusions of Chandrasekhar and von Neumann are usually taken with circumspection. In particular, Cohen et al. (1950) argue that the logarithmic divergence should be cut-off at the Jeans length which gives an estimate of the system’s size. Indeed, while in 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 size. Therefore, the Jeans length is the gravitational analogue of the Debye length. These kinetic theories lead to a collisional relaxation time scaling as , where is the dynamical time and the number of stars in the system. Chandrasekhar (1949) also developed a Brownian theory of stellar dynamics and showed that, from a qualitative point of view, the results of kinetic theory can be understood very simply in that framework. In particular, he showed that a dynamical friction is necessary to maintain the Maxwell-Boltzmann distribution of statistical equilibrium and that the coefficients of friction and diffusion are related to each other by an Einstein relation which is a manifestation of the fluctuation-dissipation theorem. This relation is confirmed by his more precise kinetic theory based on two-body encounters. It is important to emphasize, however, that Chandrasekhar did not derive the kinetic equation for the evolution of the system as a whole. Indeed, he considered the Brownian motion of a test star in a fixed distribution of field stars (bath) and derived the corresponding Fokker-Planck equation. This equation has been used by Chandrasekhar (1943b), Spitzer and Härm (1958), Michie (1963), King (1965), and more recently Lemou and Chavanis (2010) to study the evaporation of stars from globular clusters in a simple setting.

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 scales6. Furthermore, this approach clearly show which approximations are needed in order to recover the traditional Landau equation. Unfortunately, the generalized Landau equation remains extremely complicated for practical applications.

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 inhomogeneous7. If we consider a finite homogeneous system, and take collective effects into account, one finds a severe divergence of the diffusion coefficient when the size of the domain reaches the Jeans scale (Weinberg 1993). This divergence, which is related to the Jeans instability, does not occur in a stable spatially inhomogeneous stellar system. Some authors like Thorne (1968), Miller (1968), Gilbert (1968,1970), and Lerche (1971) attempted to take collective effects and spatial inhomogeneity into account. However, they obtained very complicated kinetic equations that have not found application until now. They managed, however, to show that collective effects are equivalent to increasing the effective mass of the stars, hence diminishing the relaxation time. Since, on the other hand, the effect of spatial inhomogeneity is to increase the relaxation time (Parisot and Severne 1979), the two effects act in opposite direction and may balance each other.

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 interactions8 in various dimensions of space (see Chavanis 2007,2008,2010) but we think that it is important to discuss these results in the specific case of self-gravitating systems with complements and amplification.

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 scaling9. If this scaling is satisfied, we can consider an expansion of the BBGKY hierarchy in terms of the small parameter . This is similar to the expansion of the BBGKY hierarchy is plasma physics in terms of the small parameter where represents the number of charges in the Debye sphere (Balescu 2000). However, in plasma physics, the system is spatially homogeneous (due to Debye shielding which restricts the range of interaction) while, for stellar systems, spatial inhomogeneity must be taken into account. This brings additional terms in the BBGKY hierarchy that are absent in plasma physics.

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 priori10 be neglected in the Liouvillian. This corresponds to the weak coupling approximation where only the mean field term is retained. The fifth term in Eq. (LABEL:trunc6) is a source term expressible in terms of the one-body distribution; it is of order . If we consider only the mean field Liouvillian and the source term , as we shall do in this paper, we can obtain a kinetic equation for stellar systems that is the counterpart of the Landau equation in plasma physics. The sixth term is of order and it corresponds to collective effects (i.e. dressing of the particles by the polarization cloud). In plasma physics, this term leads to the Lenard-Balescu equation. It takes into account dynamical screening and regularizes the divergence at large scales that appears in the Landau equation. In the case of stellar systems, there is no large-scale divergence because of the spatial inhomogeneity of the system. Therefore, collective effects are less crucial in the kinetic theory of stellar systems than in plasma physics. However, this term has been properly taken into account by Heyvaerts (2010) who obtained a kinetic equation of stellar systems that is the counterpart of the Lenard-Balescu equation in plasma physics. The last two terms are of the order and they will be neglected. In particular, the three-body correlation function , of order , can be neglected at the order . In this way, the hierarchy of equations is closed and a kinetic equation involving only two-body encounters can be obtained.

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 support11. Therefore, we can compute the correlation function (2.4) by assuming that the distribution function is “frozen” at time . This corresponds to the Bogoliubov ansatz in plasma physics. If we replace and by and in Eq. (23) and extend the time integral to , we obtain

(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 interaction12. It applies to weakly coupled plasmas. We note that the potential of interaction only appears in the constant which merely determines the relaxation time. The structure of the Landau equation is independent on the potential. The Landau equation was originally derived from the Boltzmann equation in the limit of weak deflections (Landau 1936)13. In the case of plasmas, the system is spatially homogeneous and the advection term is absent in Eq. (28). In the case of stellar systems, when we make the local approximation, the spatial inhomogeneity of the system is only retained in the advection term of Eq. (28). This is why this kinetic equation is referred to as the Vlasov-Landau equation. This is the fundamental kinetic equation of stellar systems.

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 cluster14.

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 same15 as for the spatially homogeneous Landau equation (Balescu 2000). Because of these properties, we might expect that a stellar system will relax towards the Boltzmann distribution which maximizes the entropy at fixed mass and energy. However, we know that there is no maximum entropy state for an unbounded self-gravitating system (the Boltzmann distribution has infinite mass). Therefore, the Vlasov-Landau equation does not relax towards a steady state and the entropy does not reach a stationary value. Actually, the entropy increases permanently as the system evaporates. But since evaporation is a slow process, the system may achieve a quasistationary state that is close to the Boltzmann distribution. A typical quasistationary distribution is the Michie-King model

(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 escape16. The characteristic time in which the system’s stars evaporate is . The evaporation is one reason for the evolution of stellar systems. However, as demonstrated by Antonov (1962) and Lynden-Bell and Wood (1968), stellar systems may evolve more rapidly as a result of the gravothermal catastrophe. In that case, the Michie-King distribution changes significatively due to core collapse. This evolution has been described by Cohn (1980), and it leads ultimately to the formation of a binary star surrounded by a hot halo. Such a configuration can have arbitrarily large entropy. Cohn (1980) finds that the entropy increases permanently during core collapse, confirming that the Vlasov-Landau equation has no equilibrium state.

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) instability17 and the gravothermal catastrophe sets in (Lynden-Bell and Wood 1968). This instability is due to the negative specific heat of the inner system that evolves by losing energy and thereby growing hotter. The energy lost is transferred outward by stellar encounters. Hence the temperature always decreases outward, and the center continually loses energy, shrinks, and heats up. This leads to core collapse. Mathematically speaking, core collapse would generate a finite time singularity. When the evolution is modeled by the orbit-averaged-Fokker-Planck equation, Cohn (1980) finds that the collapse is self-similar, that the central density becomes infinite in a finite time, and that the density behaves as . The invariant profile found by Cohn differs from the Michie-King distribution (for which ) beyond a radius of about . Larson (1970) and Lynden-Bell and Eggleton (1980) find similar results by modeling the evolution of the system by fluid equations. Alternatively, when the evolution is modeled by the original Vlasov-Landau equation, Lancellotti and Kiessling (2001) argue that the density should behave as in the final stage of the collapse. In all cases, the authors find a singular density profile at the collapse time that is integrable (or diverges logarithmically) at . This means that the core contains very little mass. In reality, if we come back to the -body system, there is no singularity and core collapse is arrested by the formation of binary stars due to three-body collisions. These binaries can release sufficient energy to stop the collapse and even drive a re-expansion of the cluster in a post-collapse regime (Inagaki and Lynden-Bell 1983). Then, in principle, a series of gravothermal oscillations should follow (Bettwieser and Sugimoto 1984).

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” stars18. Due to the encounters with the field stars, the test star has a stochastic motion. We call the probability density of finding the test star at position with velocity at time . The evolution of can be obtained from the generalized Landau equation (24) by considering that the distribution function of the field stars is fixed. Therefore, we replace by and by where is the steady distribution of the field stars19. This procedure transforms the integro-differential equation (24) into the differential equation

(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