Dressed diffusion and friction coefficients in inhomogeneous multicomponent self-gravitating systems

TL;DR

This paper derives comprehensive expressions for diffusion and friction coefficients in inhomogeneous, multicomponent self-gravitating systems, linking kinetic equations and providing tools for modeling their secular evolution.

Contribution

It introduces a self-consistent formalism for diffusion and friction in multicomponent systems, connecting Fokker-Planck, Balescu-Lenard, and Landau equations, with practical Langevin equations.

Findings

Derived general expressions for diffusion and friction coefficients.

Demonstrated the equivalence of Fokker-Planck, Balescu-Lenard, and Landau equations.

Provided Langevin equations as practical tools for simulations.

Abstract

General self-consistent expressions for the coefficients of diffusion and dynamical friction in a stable, bound, multicomponent self-gravitating and inhomogeneous system are derived. They account for the detailed dynamics of the colliding particles and their self-consistent dressing by collective gravitational interactions. The associated Fokker-Planck equation is shown to be fully consistent with the corresponding inhomogeneous Balescu-Lenard equation and, in the weak self-gravitating limit, to the inhomogeneous Landau equation. Hence it provides an alternative derivation to both and demonstrates their equivalence. The corresponding stochastic Langevin equations are presented: they can be a practical alternative to numerically solving the inhomogeneous Fokker-Planck and Balescu-Lenard equations. The present formalism allows for a self-consistent description of the secular evolution of different populations covering a spectrum of masses, with a proper accounting of the induced secular mass segregation, which should be of interest to various astrophysical contexts, from galactic centers to protostellar discs.

Equations

$$F_{\mathrm{tot}}^{\mathrm{a}}(\mathbf{x},\mathbf{v},t)=F^{\mathrm{a}}(\mathbf{x})+\delta F^{\mathrm{a}}(\mathbf{x},t);U(\mathbf{x},t)=U_0(\mathbf{x})+\delta U_{\mathrm{t}}(\mathbf{x},t)+\delta U_{\mathrm{p}}(\mathbf{x},t),$$

$$\delta F^{\mathrm{a}}\ll F^{\mathrm{a}};\delta U_{\mathrm{t}},\delta U_{\mathrm{p}}\ll U_0.$$

$$H_{\mathrm{a}}=\frac{\mu_{\mathrm{a}}}{2}\bm{v}^{2}\!+\!\mu_{\mathrm{a}}\big{[}U_{0}(\bm{x})\!+\!\delta U_{\mathrm{t}}(\bm{x},t)\!+\!\delta U_{\mathrm{p}}(\bm{x},t)\big{]}=H_{0}\!+\!\delta H_{\mathrm{t}}\!+\!\delta H_{\mathrm{p}}\,.$$

$$\mathrm{d}\Gamma =\mathrm{d}\mathbf{x}\mathrm{d}\mathbf{v}=\mathrm{d}\mathbf{\theta }\mathrm{d}\mathbf{J}.$$

$${\psi }^{(\alpha )}(\mathbf{x})=\int \mathrm{d}{\mathbf{x}}^{\prime }{\rho }^{(\alpha )}({\mathbf{x}}^{\prime })u(|\mathbf{x}-{\mathbf{x}}^{\prime }|);\int \mathrm{d}\mathbf{x}{\psi }^{(\alpha )}(\mathbf{x}){\rho }^{(\beta )\ast }(\mathbf{x})=-{\delta }_{\alpha }^{\beta },$$

$${\partial }_t{F}_{\mathrm{tot}}^{\mathrm{a}}+{\dot{\mathbf{x}}}_{\mathrm{a}}\cdot \frac{\partial F_{\mathrm{tot}}^{\mathrm{a}}}{\partial {\mathbf{x}}_{\mathrm{a}}}+{\dot{\mathbf{v}}}_{\mathrm{a}}\cdot \frac{\partial F_{\mathrm{tot}}^{\mathrm{a}}}{\partial {\mathbf{v}}_{\mathrm{a}}}=0,$$

$$H_0={\mu }_{\mathrm{a}}\mathbf{\Omega }\cdot \mathbf{J};\delta H_{\mathrm{t}}={\mu }_{\mathrm{a}}\delta U_{\mathrm{t}}({\mathbf{x}}_{\mathrm{a}}(\mathbf{\theta },\mathbf{J}),t);\delta H_{\mathrm{p}}={\mu }_{\mathrm{a}}\delta U_{\mathrm{p}}({\mathbf{x}}_{\mathrm{a}}(\mathbf{\theta },\mathbf{J}),t).$$

$${\mu }_{\mathrm{a}}{\partial }_t{F}^{\mathrm{a}}+{\mu }_{\mathrm{a}}{\partial }_t\delta F^{\mathrm{a}}+{\partial }_{\mathbf{\theta }}(F^{\mathrm{a}}+\delta F^{\mathrm{a}})\cdot {\partial }_{\mathbf{J}}(H_0+\delta H_{\mathrm{t}}+\delta H_{\mathrm{p}})-{\partial }_{\mathbf{J}}(F^{\mathrm{a}}+\delta F^{\mathrm{a}})\cdot {\partial }_{\mathbf{\theta }}(H_0+\delta H_{\mathrm{t}}+\delta H_{\mathrm{p}})=0.$$

$${\partial }_t\delta F^{\mathrm{a}}+\mathbf{\Omega }\cdot {\partial }_{\mathbf{\theta }}\delta F^{\mathrm{a}}-{\partial }_{\mathbf{J}}{F}^{\mathrm{a}}\cdot {\partial }_{\mathbf{\theta }}(\delta U_{\mathrm{t}}+\delta U_{\mathrm{p}})=0.$$

$$\delta {\rho }_{\mathrm{t}}(\mathbf{x},t)={\mu }_{\mathrm{t}}{\delta }_{\mathrm{D}}(\mathbf{x}-{\mathbf{x}}_{\mathrm{t}}(t)),$$

$$a_{\mathrm{t}}^{\alpha }(t)=-\int \mathrm{d}\mathbf{x}\delta {\rho }_{\mathrm{t}}(\mathbf{x},t){\psi }^{(\alpha )\ast }(\mathbf{x})=-{\mu }_{\mathrm{t}}{\psi }^{(\alpha )\ast }({\mathbf{x}}_{\mathrm{t}}(t)).$$

$$\delta {\rho }_{\mathrm{p}}(\mathbf{x},t)=\sum _{\mathrm{b}}\int \mathrm{d}\mathbf{v}\delta F^{\mathrm{b}}(\mathbf{x},\mathbf{v},t).$$

$$a_{\mathrm{p}}^{\alpha }(t)=-\int \mathrm{d}\mathbf{x}\delta {\rho }_{\mathrm{p}}(\mathbf{x},t){\psi }^{(\alpha )\ast }(\mathbf{x})=-\sum _{\mathrm{b}}\int \mathrm{d}\mathbf{x}\mathrm{d}\mathbf{v}\delta F^{\mathrm{b}}(\mathbf{x},\mathbf{v},t){\psi }^{(\alpha )\ast }(\mathbf{x}).$$

$${\psi }^{(\alpha )}(\mathbf{x})=\sum _{\mathbf{m}}{\psi }_{\mathbf{m}}^{(\alpha )}(\mathbf{J}){\mathrm{e}}^{\mathrm{i}\mathbf{m}\cdot \mathbf{\theta }}\text{with}{\psi }_{\mathbf{m}}^{(\alpha )}(\mathbf{J})=\frac{1}{(2\pi {)}^d}\int \mathrm{d}\mathbf{\theta }{\psi }^{(\alpha )}(\mathbf{x}(\mathbf{\theta },\mathbf{J})){\mathrm{e}}^{-\mathrm{i}\mathbf{m}\cdot \mathbf{\theta }},$$

$$a_{\mathrm{t}}^{\alpha }(t)=-{\mu }_{\mathrm{t}}\sum _{\mathbf{m}}{\psi }_{\mathbf{m}}^{(\alpha )\ast }({\mathbf{J}}_{\mathrm{t}}){\mathrm{e}}^{-\mathrm{i}\mathbf{m}\cdot ({\mathbf{\theta }}_{\mathrm{t}}^0+{\mathbf{\Omega }}_{\mathrm{t}}t)}.$$

$$\hat{f}(\omega )={\int }_{-\infty }^{+\infty }\mathrm{d}tf(t){\mathrm{e}}^{\mathrm{i}\omega t};f(t)=\frac{1}{2\pi }{\int }_{-\infty }^{+\infty }\mathrm{d}\omega \hat{f}(\omega ){\mathrm{e}}^{-\mathrm{i}\omega t}.$$

$${\hat{a}}_{\mathrm{t}}^{\alpha }(\omega )=-2\pi {\mu }_{\mathrm{t}}\sum _{\mathbf{m}}{\psi }_{\mathbf{m}}^{(\alpha )\ast }({\mathbf{J}}_{\mathrm{t}}){\mathrm{e}}^{-\mathrm{i}\mathbf{m}\cdot {\mathbf{\theta }}_{\mathrm{t}}^0}{\delta }_{\mathrm{D}}(\omega -\mathbf{m}\cdot {\mathbf{\Omega }}_{\mathrm{t}}).$$

$$a_{\mathrm{p}}^{\alpha }(t)=-\sum _{\mathrm{b}}\sum _{\mathbf{m},{\mathbf{m}}^{\prime }}\int \mathrm{d}\mathbf{\theta }\mathrm{d}\mathbf{J}\delta F_{\mathbf{m}}^{\mathrm{b}}(\mathbf{J},t){\psi }_{{\mathbf{m}}^{\prime }}^{(\alpha )\ast }(\mathbf{J}){\mathrm{e}}^{\mathrm{i}(\mathbf{m}-{\mathbf{m}}^{\prime })\cdot \mathbf{\theta }}.$$

$${\hat{a}}_{\mathrm{p}}^{\alpha }(\omega )=-(2\pi {)}^d\sum _{\mathrm{b}}\sum _{\mathbf{m}}\int \mathrm{d}\mathbf{J}\delta {\hat{F}}_{\mathbf{m}}^{\mathrm{b}}(\mathbf{J},\omega ){\psi }_{\mathbf{m}}^{(\alpha )\ast }(\mathbf{J}).$$

$$\delta U_{\mathrm{t}}(\mathbf{x},t)=a_{\mathrm{t}}^{\alpha }(t){\psi }^{(\alpha )}(\mathbf{x});\delta U_{\mathrm{p}}(\mathbf{x},t)=a_{\mathrm{p}}^{\alpha }(t){\psi }^{(\alpha )}(\mathbf{x}).$$

$$\delta\widehat{F}^{\mathrm{a}}_{\bm{m}}(\bm{J},\omega)=-\frac{\bm{m}\!\cdot\!\partial_{\bm{J}}F^{\mathrm{a}}}{\omega\!-\!\bm{m}\!\cdot\!\bm{\Omega}}\bigg{[}\delta\widehat{U}^{\mathrm{t}}_{\bm{m}}(\bm{J},\omega)\!+\!\delta\widehat{U}^{\mathrm{p}}_{\bm{m}}(\bm{J},\omega)\bigg{]}\,.$$

$$\varepsilon_{\alpha\beta}(\omega)\,\widehat{a}^{\beta}_{\mathrm{p}}(\omega)=\big{[}\delta_{\alpha\beta}\!-\!\varepsilon_{\alpha\beta}(\omega)\big{]}\,\widehat{a}^{\beta}_{\mathrm{t}}(\omega)=\widehat{\mathbf{M}}_{\alpha\beta}(\omega)\,\widehat{a}^{\beta}_{\mathrm{t}}(\omega)\,,$$

$$\varepsilon_{\alpha\beta}(\omega)=\delta_{\alpha\beta}-\sum_{\bm{m}}(2\pi)^{d}\!\!\int\!\!\mathrm{d}\bm{J}\,\frac{\bm{m}\!\cdot\!\partial_{\bm{J}}\big{[}\sum_{\mathrm{b}}F^{\mathrm{b}}(\bm{J})\big{]}}{\omega\!-\!\bm{m}\!\cdot\!\bm{\Omega}}\,\psi^{(\alpha)*}_{\bm{m}}(\bm{J})\,\psi^{(\beta)}_{\bm{m}}(\bm{J})\,.$$

$${\hat{a}}_{\mathrm{d}}^{\alpha }(\omega )={\hat{a}}_{\mathrm{t}}^{\alpha }(\omega )+{\hat{a}}_{\mathrm{p}}^{\alpha }(\omega )={\epsilon }_{\alpha \beta }^{-1}(\omega ){\hat{a}}_{\mathrm{t}}^{\beta }(\omega ),$$

$$\delta {\hat{U}}_{\mathrm{d}}(\mathbf{x},\omega )={\psi }^{(\alpha )}(\mathbf{x}){\epsilon }_{\alpha \beta }^{-1}(\omega ){\hat{a}}_{\mathrm{t}}^{\beta }(\omega )=-2\pi {\mu }_{\mathrm{t}}\sum _{\mathbf{m}}{\psi }^{(\alpha )}(\mathbf{x}){\epsilon }_{\alpha \beta }^{-1}(\omega ){\psi }_{\mathbf{m}}^{(\beta )\ast }({\mathbf{J}}_{\mathrm{t}}){\mathrm{e}}^{-\mathrm{i}\mathbf{m}\cdot {\mathbf{\theta }}_{\mathrm{t}}^0}{\delta }_D(\omega -\mathbf{m}\cdot {\mathbf{\Omega }}_{\mathrm{t}}).$$

$$\delta U_{\mathrm{d}}(\mathbf{x},t)=-{\mu }_{\mathrm{t}}\sum _{\mathbf{m}}{\psi }^{(\alpha )}(\mathbf{x}){\epsilon }_{\alpha \beta }^{-1}(\mathbf{m}\cdot {\mathbf{\Omega }}_{\mathrm{t}}){\psi }_{\mathbf{m}}^{(\beta )\ast }({\mathbf{J}}_{\mathrm{t}}){\mathrm{e}}^{-\mathrm{i}\mathbf{m}\cdot ({\mathbf{\theta }}_{\mathrm{t}}^0+{\mathbf{\Omega }}_{\mathrm{t}}t)}.$$

$$\delta U_{\mathrm{bare}}(\mathbf{x},t)=-{\mu }_{\mathrm{t}}\sum _{\mathbf{m}}{\psi }^{(\alpha )}(\mathbf{x}){\psi }_{\mathbf{m}}^{(\alpha )\ast }({\mathbf{J}}_{\mathrm{t}}){\mathrm{e}}^{-\mathrm{i}\mathbf{m}\cdot ({\mathbf{\theta }}_{\mathrm{t}}^0+{\mathbf{\Omega }}_{\mathrm{t}}t)}.$$

$${\psi }^{(\alpha )}({\mathbf{x}}_1)=\sum _{{\mathbf{m}}_1}{\psi }_{{\mathbf{m}}_1}^{(\alpha )}({\mathbf{J}}_1){\mathrm{e}}^{\mathrm{i}{\mathbf{m}}_1\cdot {\mathbf{\theta }}_1},$$

$$\delta U_{2|1}^{\mathrm{d}}(1,t)=-\sum _{{\mathbf{m}}_1,{\mathbf{m}}_2}{\mu }_2{\mathrm{e}}^{\mathrm{i}({\mathbf{m}}_1\cdot {\mathbf{\theta }}_1-{\mathbf{m}}_2\cdot {\mathbf{\theta }}_2)}{\psi }_{{\mathbf{m}}_1}^{(\alpha )}({\mathbf{J}}_1){\epsilon }_{\alpha \beta }^{-1}({\mathbf{m}}_2\cdot {\mathbf{\Omega }}_2){\psi }_{{\mathbf{m}}_2}^{(\beta )\ast }({\mathbf{J}}_2).$$

$$\frac{1}{{\mathcal{D}}_{{\mathbf{m}}_1,{\mathbf{m}}_2}({\mathbf{J}}_1,{\mathbf{J}}_2,\omega )}={\psi }_{{\mathbf{m}}_1}^{(\alpha )}({\mathbf{J}}_1){\epsilon }_{\alpha \beta }^{-1}(\omega ){\psi }_{{\mathbf{m}}_2}^{(\beta )\ast }({\mathbf{J}}_2),$$

Full text

††pagerange: Dressed diffusion and friction coefficients in inhomogeneous multicomponent self-gravitating systems–F.2

1 Introduction and motivation

The kinetic theory of stellar systems was initiated by Chandrasekhar (1942). He first described the motion of a star in a stellar system semi-heuristically by using an analogy with Brownian motion (Chandrasekhar, 1943a). He argued that the force acting on a star has two components: a mean field component due to the smooth distribution of the system and a fluctuating component arising from discreteness effects. Discreteness effects (also called finite${-N}$ effects, granularities, graininess…) account for gravitational encounters. For a spatially homogeneous system, the mean field force vanishes so that only gravitational encounters (“collisions”) can produce an evolution. The fluctuating force acting on a star gives rise to diffusive motion in velocity space. However, a purely diffusive motion would lead to a divergence of the kinetic energy of the star and would not establish a statistical equilibrium state at late time. As a result, Chandrasekhar realised that something was “missing” in his description and that the diffusive process must be accompanied by a dissipative process. He concluded that the star must also experience a dynamical friction. Using a phenomenological Langevin equation incorporating a friction force proportional and opposite to the velocity of the star and a random force modeled as a white noise, he derived a Fokker-Planck equation describing the evolution of the velocity distribution function ${f(\bm{v},t)}$ of the star. Then, requiring that the Maxwell-Boltzmann distribution must be a stationary state of this Fokker-Planck equation, he showed that the coefficients of diffusion and friction are related to each other by an Einstein relation. This is the manifestation of the fluctuation-dissipation theorem. In parallel, using a more rigorous kinetic theory, he directly computed the coefficients of diffusion and friction in the approximation of close binary encounters and explicitly checked that, for a Maxwellian distribution, the Einstein relation is indeed satisfied. He then used his kinetic theory to estimate the rate of escape of stars from clusters (Chandrasekhar, 1943b), and found observational evidence for the operation of dynamical friction.

In the approach of Chandrasekhar, further developed by Cohen et al. (1950); Rosenbluth et al. (1957); King (1960); Binney & Tremaine (2008), the diffusion and friction coefficients are obtained independently, from two different calculations, and are then injected into the Fokker-Planck equation. When the distribution function is allowed to change self-consistently, one obtains an integrodifferential equation describing the evolution of the distribution function of the stellar system as a whole. This integrodifferential equation is formally equivalent to the Landau equation that was introduced earlier in the different context of plasma physics (Landau, 1936).111The Landau equation was sometimes misunderstood. For example, the paper of Cohen et al. (1950) comments: “A similar but incomplete approach [to the work of Chandrasekhar] was made somewhat earlier by Landau. In this reference, the important terms representing dynamical friction, which should appear in the diffusion equation, are set equal to zero as a result of certain approximations.” This claim is misleading since the Landau equation includes both terms of diffusion and friction.

Indeed, a few years before Chandrasekhar’s seminal papers on stellar dynamics, Landau derived an integrodifferential kinetic equation for homogeneous Coulombian plasmas interacting via inverse-square forces. He obtained this equation as a weak deflection approximation of the Boltzmann equation. In this equation, the terms of diffusion and friction arise simultaneously from a unique formalism. The Landau equation has a form similar to the Fokker-Planck equation except that the diffusion tensor is placed between the two velocity derivatives, ${\partial_{t}f\!=\!\partial_{v_{i}}(D_{ij}\partial_{v_{j}}f)\!-\!\partial_{v_{i}}(F^{\rm pol}_{i}\!f)}$, while in the standard Fokker-Planck equation it is placed after the two velocity derivatives, ${\partial_{t}f\!=\!\partial_{v_{i}}\partial_{v_{j}}(D_{ij}f)\!-\!\partial_{v_{i}}(F^{\rm fric}_{i}\!f)}$. As a result, the friction force that appears in the Landau equation differs from the friction force that appears in the Fokker-Planck equation of Chandrasekhar. We shall call it the “friction by polarisation” $\bm{F}_{\rm pol}$, to distinguish it from the “true friction” $\bm{F}_{\rm fric}$. Of course, one can immediately transform the Landau equation into the standard Fokker-Planck equation and find the relation ${F^{\rm fric}_{i}\!=\!F^{\rm pol}_{i}\!+\!\partial_{v_{j}}D_{ij}}$ between the two friction forces. One can then check (see, e.g. Chavanis, 2013) that the Landau equation is fully equivalent to the Fokker-Planck equation of Chandrasekhar even if the equations do not appear in the same form.222A virtue of the Landau equation resides in its symmetric structure from which we can immediately derive the conservation laws (mass, energy, impulse, angular momentum) and the ${H-}$theorem for the Boltzmann entropy.

The force of dynamical friction was calculated by Marochnik (1968); Kalnajs (1971); Kandrup (1983); Bekenstein & Zamir (1990); Chavanis (2008) from a linear response theory based on the Liouville or on the Klimontovich equation. In these calculations, it results from a polarisation process. A test star perturbs the distribution of the field stars and the retroaction of the field stars on the test star leads to a friction force that decelerates the test star. The force of dynamical friction calculated by these authors differs from that calculated by Chandrasekhar by a factor of two (or by a factor ${(m\!+\!m_{\rm f})/m}$ if the test star and the field stars have different masses). This is because, as noted in Chavanis (2013), they actually calculated the “friction by polarisation” (the one that arises in the Landau equation), not the “true friction” (the one that arises in the Fokker-Planck equation of Chandrasekhar).

The kinetic theory of Chandrasekhar is based on two simplifying assumptions. He assumed that the system is spatially homogeneous and neglected collective effects (the fact that a star is surrounded by a cloud of other stars that tends to enhance the gravitational attraction). Various generalisations and improvements of the theory of Chandrasekhar (1943a) were subsequently proposed. Some authors tried to deal with spatial inhomogeneity. Kandrup (1983); Bekenstein & Maoz (1992); Maoz (1993); Nelson & Tremaine (1999); Chavanis (2008) reconsidered the interaction of a test particle with a background stochastic force in the context of the fluctuation-dissipation theorem, and showed that the friction force depends on the global structure of the system. This was also investigated in Del Popolo & Gambera (1999); Del Popolo (2003) which extended Chandrasekhar & von Neumann (1943)’s analysis to the case where the background particles are inhomogeneously distributed with a density profile decaying as $\rho\sim r^{-p}$. Kinetic theories for spatially inhomogeneous systems were developed in Severne & Haggerty (1976); Parisot & Severne (1979); Kandrup (1981); Chavanis (2008, 2013) using position-velocity variables. Collective effects were first taken into account in plasma physics, where the system is spatially homogeneous because of electroneutrality. Balescu (1960) and Lenard (1960) derived a generalisation of the Landau equation by accounting properly for collective effects (the fact that a charge is surrounded by a cloud of opposite charges that tends to shield the electrostatic interaction). In the Balescu-Lenard equation, the bare potential of interaction is replaced by a dressed potential of interaction that takes into account the dressing of the particles by their polarisation cloud. As a result, the Debye length appears naturally without any ad hoc assumptions and regularises the logarithmic divergence at large scales that occurs in the Landau equation when collective effects are neglected. Balescu (1960) started from the diagram technique introduced by Prigogine & Balescu (1959), and Lenard (1960) started from the Bogoliubov (1946) equations (now known as the BBGKY hierarchy). Independently, Hubbard (1961) used a Fokker-Planck approach and directly calculated the coefficients of diffusion and friction by taking collective effects into account.333If one substitutes these expressions into the Fokker-Planck equation and performs minor transformations (a substitution that Hubbard did not explicitly make), one obtains the Balescu-Lenard equation. Inversely, from the Balescu-Lenard equation, one can recover the expressions of the diffusion and friction coefficients obtained by Hubbard. These considerations show that the self-consistent Fokker-Planck equation of Hubbard is fully equivalent to the Balescu-Lenard equation (see, e.g. Chavanis, 2012b).

In order to take collective effects into account in stellar systems, one must simultaneously account for their spatial inhomogeneity otherwise the kinetic equation presents a strong (algebraic) divergence at large-scales related to the Jeans instability (Weinberg, 1993; Chavanis, 2013). Kinetic theories for spatially inhomogeneous systems with position-velocity variables taking collective effects into account were developed by Miller (1966); Thorne (1968); Gilbert (1968, 1970); Lerche (1971). Their kinetic equations are very complicated but they managed to show that collective effects are equivalent to increasing the effective mass of the stars, hence reducing the relaxation time. Spatial inhomogeneity can conveniently be dealt with by using angle-action variables for integrable systems (Goldstein, 1950; Binney & Tremaine, 2008). On the other hand, collective effects can be dealt with by introducing a biorthogonal basis of potentials and densities and using Kalnajs (1976)’s matrix method. The friction force and the energy exchange rate for spatially inhomogeneous stellar systems with angle-action variables were calculated by Lynden-Bell & Kalnajs (1972); Tremaine & Weinberg (1984); Palmer & Papaloizou (1985); Weinberg (1986, 1989); Rauch & Tremaine (1996). The corresponding inhomogeneous Balescu-Lenard equation has been derived by Heyvaerts (2010) from the BBGKY hierarchy and by Chavanis (2012c) from the Klimontovich approach.444The inhomogeneous Balescu-Lenard equation with angle-action variables was previously derived by Luciani & Pellat (1987), where the response matrix is introduced at a formal level without an explicit representation. Chavanis (2012c) developed a Fokker-Planck approach and directly calculated the diffusion and friction coefficients with angle-action variables by taking collective effects into account, thereby generalising the results of Hubbard (1961) to inhomogeneous systems. He also developed a test particle approach relying on a (thermal) bath approximation in which the integrodifferential Balescu-Lenard equation is transformed into a differential Fokker-Planck equation. Neglecting collective effects, one recovers the results obtained from the inhomogeneous Landau equation (Chavanis, 2007, 2013). Making a local approximation, one recovers the original results of Landau (1936) and Chandrasekhar (1942, 1943a).

The inhomogeneous Balescu-Lenard equation was recently implemented for the first time in astrophysics. In Fouvry et al. (2015b); Fouvry et al. (2015c, 2017b), it was applied to razor-thin and thickened stellar discs, and proved useful to probe complex secular regimes of diffusion. These works showed in particular how collective effects cause cool discs to have two-body relaxation times much shorter than naively expected. In addition, they showed how this relaxation introduces small-scale structures in the disc, which secularly destabilise it at the collisionless level. The inhomogeneous Balescu-Lenard equation was also recently applied in the context of the ${1D}$ inhomogeneous Hamiltonian Mean Field (HMF) model (Benetti & Marcos, 2017). Finally, it was specialised to dynamically degenerate systems, such as quasi-Keplerian systems (galactic centers, protostellar discs) in Fouvry et al. (2017a) and Sridhar & Touma (2017) (without collective effects). Following these recent successes, the inhomogeneous Balescu-Lenard equation now appears as a powerful and predictive framework. It could also be used as a valuable check of the accuracy of ${N-}$body integrators on secular timescales.

In this paper, we present a derivation of the coefficients of diffusion and dynamical friction for self-gravitating stellar systems, while taking into account both inhomogeneity and collective effects. We also present the set of stochastic Langevin equations, dual to the Fokker-Planck equation, that can be a powerful alternative to numerically solving the inhomogeneous Balescu-Lenard equation through stochastic ${N-}$body techniques. Compared to the derivation of Chavanis (2012c), the present formalism presents several advantages: it shows more clearly the dressing of particles by their polarisation cloud, as we calculate from the start an explicit expression of the dressed potential of a moving particle. Collective effects are obtained by solving the linearised Klimontovich equation. As a result, the approximation where collective effects are neglected is straightforwardly recovered. In addition, this calculation also takes into account the possibility for particles to have different masses, a generalisation that is particularly important for astrophysical applications, where mass segregation is deemed to play a significant role on secular timescales. In the single component case, the results of Chavanis (2012c) are recovered.

The paper is organised as follows. Section 2 presents the dressed potential of a test particle. Section 3 computes the diffusion coefficient in action space. Section 4 focuses on the coefficient of dynamical friction. Section 5 relates the corresponding Fokker-Planck equation to the Balescu-Lenard equation. Section 6 provides the equivalent stochastic Langevin equations and stresses their practical interest. Section 7 discusses the different stages in the evolution of stellar systems in the light of kinetic theory, and finally, section 8 wraps up. Appendix A computes the drift vector. Appendix B gives a direct derivation of the friction force by polarisation. Appendix C derives the main properties of the multicomponent Balescu-Lenard equation (steady states, energy conservation, ${H-}$theorem). Appendix D develops the corresponding test particle approach in a bath formalism and considers the thermal bath and the sinking satellite problems. Appendices E and F provide mathematical results needed in the calculations.

2 Dressed potential of a moving particle

When considering self-gravitating systems, two main difficulties have to be overcome. First, self-gravitating systems are spatially inhomogeneous, which makes the trajectories of individual particles intricate. Assuming the mean system to be integrable, one can deal with inhomogeneity thanks to the use of angle-action coordinates. The second difficulty arises from the system’s self-gravity, i.e. its ability to amplify perturbations. Dealing with these collective effects requires to study the dressing of fluctuations on dynamical times. Let us therefore first compute the gravitational polarisation induced by a moving particle in a self-gravitating system.

2.1 Notations

We consider a test particle moving in a self-gravitating inhomogeneous system, composed of various components denoted with “$\mathrm{a}$”, “$\mathrm{b}$”, etc. We note as ${F^{\mathrm{a}}_{\rm tot}(\bm{x},\bm{v})}$ the distribution function (DF) of particles of component “$\mathrm{a}$” with individual mass $\mu_{\mathrm{a}}$. The DFs are normalised such that ${\!\int\!\mathrm{d}\bm{x}\mathrm{d}\bm{v}\,F^{\mathrm{a}}_{\rm tot}\!=\!M_{\rm tot}^{\mathrm{a}}}$, where ${M_{\rm tot}^{\mathrm{a}}\!=\!N_{\mathrm{a}}\mu_{\mathrm{a}}}$ is the total active mass of the component “$\mathrm{a}$” composed of $N_{\mathrm{a}}$ particles. The test particle is denoted by a subscript “$\mathrm{t}$”, which should not be mixed up with the time $t$. The system’s total gravitational potential is written as ${U(\bm{x},t)}$. The DFs and the potential can be decomposed as

[TABLE]

where $F^{\mathrm{a}}$ only varies on secular times and stands for the unperturbed part of the DF, while ${\delta F^{\mathrm{a}}}$ depends on time and is associated with the polarisation cloud surrounding the test particle. Here, $U_{0}$ stands for the system’s mean potential and carries on secular timescales, ${\delta U_{\mathrm{t}}}$ is the potential created by the test particle, and finally ${\delta U_{\mathrm{p}}}$ is the polarisation cloud induced by the test particle. We assume that all perturbations are small, so that

[TABLE]

The dynamics of a particle of component “$\mathrm{a}$” is given by the total Hamiltonian $H_{\mathrm{a}}$ reading

[TABLE]

For simplicity, the phase coordinates are denoted by ${\Gamma\!=\!(\bm{x},\bm{v})}$. We assume the system to be quasi-stationary and integrable. As a consequence, one can introduce angle-action coordinates ${(\bm{\theta},\bm{J})}$ (Binney & Tremaine, 2008). Within these coordinates, the system’s mean potential becomes ${U_{0}\!=\!U_{0}(\bm{J})}$ while, following Jeans’ theorem, the mean quasi-stationary DFs $F^{\mathrm{a}}$ depend also only on the actions $\bm{J}$. These coordinates are canonical and conserve infinitesimal volumes, so that

[TABLE]

The gradients w.r.t. to a vectorial variable such as $\bm{x}$ are denoted by $\partial_{\bm{x}}$, while the derivative w.r.t. the time $t$ will sometimes be written as $\partial_{t}$ for brevity. We also rely on the matrix method (Kalnajs, 1976) and introduce a biorthogonal basis of potentials and densities ${(\psi^{(\alpha)},\rho^{(\alpha)})}$ satisfying

[TABLE]

where ${u(\bm{x})}$ corresponds to the interaction potential, i.e. ${u(|\bm{x}|)\!=\!-G/|\bm{x}|}$ in the gravitational case. Any perturbation in the system can then be represented using these basis elements. For example, a density field and the corresponding potential can be written as ${\delta\rho(\bm{x},t)\!=\!a^{\alpha}(t)\,\rho^{(\alpha)}(\bm{x})}$ and ${\delta U(\bm{x},t)\!=\!a^{\alpha}(t)\,\psi^{(\alpha)}(\bm{x})}$, where the sum over $\alpha$ is implied.

2.2 Dressing of a test particle

The Klimontovich equation (Klimontovich, 1967) reads symbolically for the one particle DF, ${F_{\rm tot}^{\mathrm{a}}\!=\!\sum_{i}\mu_{\mathrm{a}}\delta_{\mathrm{D}}(\bm{x}\!-\!\bm{x}_{i}(t))\delta_{\mathrm{D}}(\bm{v}\!-\!\bm{v}_{i}(t))}$:

[TABLE]

where individual motions are given by Hamilton’s equations reading ${\mu_{\mathrm{a}}\dot{\bm{x}}_{\mathrm{a}}\!=\!\partial H_{\mathrm{a}}/\partial\bm{v}_{\mathrm{a}}}$ and ${\mu_{\mathrm{a}}\dot{\bm{v}}_{\mathrm{a}}\!=\!-\partial H_{\mathrm{a}}/\partial\bm{x}_{\mathrm{a}}}$. When the motion of the stars is described by the angle-action variables ${(\bm{J},\bm{\theta})}$, associated with the mean potential $U_{0}$, the components of the Hamiltonian $H_{\mathrm{a}}$ from equation (3) read

[TABLE]

Following equation (6), the Klimontovich equation for the total DF ${F^{\mathrm{a}}_{\rm tot}\!=\!F^{\mathrm{a}}\!+\!\delta F^{\mathrm{a}}}$ becomes

[TABLE]

Taking into account that ${\partial_{\bm{\theta}}H_{0}\!=\!0}$, ${\partial_{\bm{\theta}}F^{\mathrm{a}}\!=\!0}$, and ${\partial_{\bm{J}}H_{0}\!=\!\mu_{\mathrm{a}}\bm{\Omega}}$, and retaining only linear terms (quasilinear approximation) (here ${\partial_{t}F^{\mathrm{a}}}$ is a second order term that can be neglected), equation (8) becomes

[TABLE]

Equation (9) is the linearised Klimontovich equation and describes the amplification of perturbations on dynamical timescales.

2.3 Dressed potential of a test particle

In this section, let us compute the dressed potential generated by a given test particle. The density ${\delta\rho_{\mathrm{t}}}$ generated by the test particle is straightforwardly expressed as

[TABLE]

where $\mu_{\mathrm{t}}$ stands for the mass of the test particle, and ${\bm{x}_{\mathrm{t}}(t)}$ for its position at time $t$. Here $\delta_{\mathrm{D}}$ is the Dirac function. One may then expand this density on the basis of densities so that ${\delta\rho_{\mathrm{t}}(\bm{x},t)\!=\!a_{\mathrm{t}}^{\alpha}(t)\,\rho^{(\alpha)}(\bm{x})}$. One has

[TABLE]

The polarisation potential ${\delta U_{\mathrm{p}}}$ is self-consistently generated by the contributions of the perturbation DFs, ${\delta F^{\mathrm{b}}}$, so that the associated density perturbation ${\delta\rho_{\mathrm{p}}}$ reads

[TABLE]

The associated coefficients in the potential-density expansion ${\delta\rho_{\mathrm{p}}(\bm{x},t)\!=\!a_{\mathrm{p}}^{\alpha}(t)\,\rho^{(\alpha)}(\bm{x})}$ read

[TABLE]

Relying on the ${2\pi-}$periodicity of the angles $\bm{\theta}$, one can perform harmonic (Fourier) expansions of these expressions. Let us first decompose the basis elements as

[TABLE]

where we introduced $d$ as the dimension of the physical space. In the angle-action coordinates, the unperturbed motion of the test particle is given by ${\bm{\theta}_{\mathrm{t}}\!=\!\bm{\theta}_{\mathrm{t}}^{0}\!+\!\bm{\Omega}_{\mathrm{t}}t}$ and ${\bm{J}_{\mathrm{t}}\!=\!\text{cst.}}$, so that equation (11) becomes

[TABLE]

Let us now introduce the temporal Fourier transform with the convention

[TABLE]

Starting from equation (15), the temporal Fourier transform of ${a^{\alpha}_{\mathrm{t}}(t)}$ reads

[TABLE]

The coefficients produced by the polarisation cloud of the test particle, $a^{\alpha}_{\mathrm{p}}$, can similarly be expressed via the harmonic Fourier transforms of the DFs, ${\delta F^{\mathrm{b}}_{\bm{m}}(\bm{J},t)}$, so that equation (13) gives

[TABLE]

The integral over ${\mathrm{d}\bm{\theta}}$ yields ${(2\pi)^{d}\delta_{\bm{m}}^{\bm{m}^{\prime}}}$. The Fourier transform of equation (18) w.r.t. time then reads

[TABLE]

Finally, the potentials ${\delta U_{\mathrm{t}}}$ and ${\delta U_{\mathrm{p}}}$ can straightforwardly be expressed in terms of the $a_{\mathrm{t}}^{\alpha}$ and $a_{\mathrm{p}}^{\alpha}$ coefficients as

[TABLE]

The next step involves computing the coefficients $\widehat{a}^{\alpha}_{\mathrm{p}}$ as a function of the parameters of the test particle. This is possible by solving equation (9), which describes the self-gravitating amplification of perturbations on dynamical timescales. When Fourier transformed w.r.t. the angles and time, it yields

[TABLE]

Substituting equation (20) into equation (21), ${\delta\widehat{F}^{\mathrm{a}}_{\bm{m}}(\bm{J},\omega)}$ can be expressed as a function of the $\widehat{a}_{\mathrm{p}}^{\alpha}$ coefficients. Equation (19) then becomes a self-consistent equation for these coefficients. After a few transformations, it can be recast as

[TABLE]

where the sum over $\beta$ is implied. In equation (22), we introduced the usual Kronecker symbol $\delta_{\alpha\beta}$, and the system’s susceptibility ${\varepsilon_{\alpha\beta}}$ as

[TABLE]

Finally, in equation (22), we also introduced the system’s response matrix ${\widehat{\mathbf{M}}(\omega)}$, which controls the strength of the self-gravitating amplification in the system. The total dressed potential, ${\delta U_{\mathrm{d}}}$, defined as ${\delta U_{\mathrm{d}}\!=\!\delta U_{\mathrm{p}}\!+\!\delta U_{\mathrm{t}}}$ immediately follows. Its coefficients ${a_{\mathrm{d}}\!=\!a_{\mathrm{t}}\!+\!a_{\mathrm{p}}}$ read

[TABLE]

where we assumed that the system is linearly stable so that ${\varepsilon_{\alpha\beta}(\omega)}$ can be inverted. In conclusion, thanks to the coefficients ${\widehat{a}^{\beta}_{\mathrm{t}}(\omega)}$ from equation (17), the total dressed potential perturbation ${\delta U_{\mathrm{d}}}$ reads

[TABLE]

Written as a function of time, equation (25) becomes

[TABLE]

Note that the bare potential perturbation ${\delta U_{\rm bare}\!=\!\delta U_{\mathrm{t}}}$ has the same functional form with $\varepsilon$ being replaced by the identity, so that

[TABLE]

The dressed potential, ${\delta U_{2|1}^{\mathrm{d}}}$, created by particle 2 and felt by particle 1 with action $\bm{J}_{1}$ and angle $\bm{\theta}_{1}$, is given by equation (26) when taking for ${\psi^{(\alpha)}(\bm{x}_{1})}$ the sum

[TABLE]

and replacing the test particle “$\mathrm{t}$” by the particle 2. After some reordering, equation (26) gives

[TABLE]

Let us finally introduce the system’s dressed susceptibility coefficients ${1/\mathcal{D}_{\bm{m}_{1},\bm{m}_{2}}(\bm{J}_{1},\bm{J}_{2},\omega)}$ as

[TABLE]

so that one can write

[TABLE]

With the shorthand notations

[TABLE]

equation (31) becomes

[TABLE]

The total potential created by the sum over all discrete particles in the system, each of them having its own angle-action variables, is then simply the sum over all the individual contributions given by equation (33), so that

[TABLE]

Under the assumptions made in the present derivation, the discrete noise is small which means that the fluctuating part in $U^{\rm di}_{1}$ is small (compared to the typical kinetic energy of a given particle). This implies that ${\sum_{\bm{m}_{1}}\sum_{\bm{m}_{2}\neq 0}\Lambda_{\bm{m}_{1},\bm{m}_{2}}(\bm{J}_{1},\bm{J}_{2},\omega_{2})}$ should be treated like a small contribution. The smallness in the noise level should therefore be measured as a function of the number of $\Lambda$ in an expression. For example, $U^{\rm di}_{1}$ in equation (34), is expressed to first order in the noise level.

3 Diffusion coefficients in action space

The evolution of the action of particle 1 under the influence of the force applied by all other discrete particles is captured by the Hamiltonian ${H_{1}\!=\!\mu_{1}U_{1}^{\rm di}}$. The associated Hamilton’s equation for the action reads

[TABLE]

Given equation (34), it becomes

[TABLE]

which can be integrated for a time ${\Delta t}$ to give

[TABLE]

Equation (37) requires the full knowledge of the particles’ motion in order to account for the explicit time dependencies in the integral term.

The frequency spectrum of fluctuations generated by the particles 2 is associated will all harmonics ${\bm{m}_{2}\!\cdot\!\bm{\Omega}_{2}}$. Provided the system is not dynamically degenerate, its zero frequency component corresponds to ${\bm{m}_{2}\!=\!\bm{0}}$. Its fluctuating part corresponds to non-zero values of $\bm{m}_{2}$ in equation (37). Summation over all but null vectors will be represented with a dash. The action diffusion tensor is an average which will be written as {\big{<}\displaystyle\Delta\bm{J}_{1}\!\otimes\!\Delta\bm{J}_{1}/\Delta t\big{>}}. This average is carried on the dynamical variables of the particles 2 and on the angles of particle $1$. It enters the Fokker-Planck equation for the DF in action space for particle 1 and component “$\mathrm{a}$”. As always in Fokker-Planck theory, this diffusion coefficient must be evaluated at second order in the noise level (Risken, 1996). Since equation (37) is by construction first order, and, since the diffusion tensor is quadratic in ${\Delta\bm{J}}$, one should compute other involved expression, e.g., ${\mathrm{e}^{z_{12}(t)}\!=\!\mathrm{e}^{\mathrm{i}(\bm{m}_{1}\cdot\bm{\theta}_{1}(t)-\bm{m}_{2}\cdot\bm{\theta}_{2}(t))}}$ to zeroth order. For this term, this corresponds to the uniform angular motion, for which one can write

[TABLE]

Similarly, the susceptibility coefficient ${\Lambda(...;t)}$ can be evaluated at ${t\!=\!0}$. At the level of this approximation, one therefore gets

[TABLE]

While omitting that slowly varying variables should be evaluated at ${t\!=\!0}$, it follows from equation (39) that

[TABLE]

One may now average this expression over the initial angles of the particles $1$, $2$, and $3$, and on the action distribution of particles $2$ and $3$. Here, one should pay attention to the fact that particle $1$ acts as our test star, while particles $2$ and $3$ both run over the field stars, i.e. over all stars except particle $1$. Let us first perform an average over the initial angles of the test and field stars. Keeping only the dependencies w.r.t. the initial angles, equation (40) requires to consider a term generically of the form

[TABLE]

where it is important to note that the sum on particles $2$ and $3$ is restricted to all the field stars for particle $1$, i.e. particles $2$ and $3$ are always different from particle $1$. Because only non-zero values of $\bm{m}_{2}$ and $\bm{m}_{3}$ contribute to the fluctuations, equation (41) immediately imposes for particle $2$ and $3$ to be the same, so that the sum $\sum_{3}$ can straightforwardly be executed. Averaging over $\bm{\theta}_{1}^{0}$, $\bm{\theta}_{2}^{0}$, and over the action distribution of particle $2$, then amounts to performing in equation (40) the replacement

[TABLE]

where the presence of the prefactor ${1/\mu_{\mathrm{b}}}$ is associated with the chosen normalisation of the DFs, ${\!\int\!\mathrm{d}\bm{x}\mathrm{d}\bm{v}F^{\mathrm{b}}\!=\!M_{\rm tot}^{\mathrm{b}}\!=\!N_{\mathrm{b}}\mu_{\mathrm{b}}}$. In equation (40), the average over $\bm{\theta}_{1}^{0}$ yields $\delta_{\bm{m}_{1}}^{-\bm{m}_{1}^{\prime}}$, while the integration over the initial angle $\bm{\theta}_{2}^{0}$ gives ${(2\pi)^{d}\delta_{\bm{m}_{2}}^{-\bm{m}_{3}}}$. Equation (40) becomes

[TABLE]

Now symmetries imply that (see Appendix F)

[TABLE]

Finally, the time limit can be carried using the relation (see Appendix E)

[TABLE]

It follows that

[TABLE]

Equation (46) is the final expression of the diffusion tensor acting on the action vector of the test particle.

4 Coefficients of dynamical friction

4.1 Mean potential and shot noise

In section 2, the dressed potential created by all discrete particles was computed. This potential can be split into a quasi-stationnary part, $U_{\rm st}$, which is essentially the mean ensemble potential, $U_{0}$, and a fluctuation, ${\delta\widetilde{U}}$, coming from the discrete particles. The real potential, including the noise contributions, which is felt by particle 1 is given by

[TABLE]

The potential $U_{0}$ is the average, computed from the one point DF of the zero frequency component of the discrete potential, i.e. $U_{\rm st}$. This averaging may differ subtly from this zero frequency part, in as much as it might display coarse grained features induced by textures in the action distribution which are not taken into account by the DF. Here, we will assume that the two potentials can be identified. The power spectrum of fluctuations caused by particle 2 is given by all the harmonics ${\bm{m}_{2}\!\cdot\!\bm{\Omega}_{2}}$. Its zero frequency contribution is given by the ${\bm{m}_{2}\!=\!\bm{0}}$ component, provided the system’s potential is not dynamically degenerate. The fluctuating part induced by discrete particles corresponds then to the non-zero values of $\bm{m}_{2}$ in equation (34). The components corresponding to ${\bm{m}_{1}\!=\!\bm{0}}$ do not impact the variation of the actions, as can be seen in equation (36). In contrast, they could impact the angles, which do play a role in what follows. The potential fluctuations corresponding to the discrete nature of the particles are accounted for in equation (34) by the ${\bm{m}_{2}\!\neq\!\bm{0}}$ and any $\bm{m}_{1}$ contribution in the sum. As previously, summations over all but null vectors will be represented by a dash. When identifying the mean potential with the zero frequency potential, the Hamiltonian for the motion of particle 1 in the presence of shot noise from the other particles then reads

[TABLE]

where we relied on equation (34). Without any ensemble average at this stage, the dynamical evolution of particle $1$ under the influence of the Hamiltonian (48) is given by the following differential equations

[TABLE]

4.2 Mean friction at second order

Let us now compute the mean drag {\big{<}\Delta\bm{J}_{1}\big{>}} applied onto particle 1 during ${\Delta t}$. This change in the vector action must be computed to second order in the level of the noise. Since equation (49) is clearly first order only, this implies that we cannot rely on zeroth order approximation for the ${\Lambda_{\bm{m}_{1},\bm{m}_{2}}(\bm{J}_{1},\bm{J}_{2},\bm{m}_{2}\!\cdot\!\bm{\Omega}_{2};t)}$ and ${\mathrm{e}^{\mathrm{i}z_{12}(t)}}$ factors.555Should we do so, {\big{<}\Delta\bm{J}_{1}\big{>}} would vanish identically. They must then be computed at the next order following, e.g., Ecker (2013).

The calculation at second order of the change in action {\big{<}\Delta\bm{J}_{1}\big{>}} is a somewhat technical calculation, that we present in detail in Appendix A. In the same Appendix, we also detail how one may average this drift vector over the initial angles of the involved particles. One finally obtains in equation (93) that the averaged drift vector is given by

[TABLE]

Equation (50) is the final expression of the dynamical friction coefficient acting on the action vector of the test particle.

5 From Fokker-Planck to Balescu-Lenard

Let us now show how the Fokker-Planck equation based on the friction and diffusions coefficients obtained previously is in fact fully consistent with the Balescu-Lenard equation (see also Appendix D for a discussion of the test particle approach). The Fokker-Planck equation for the system’s DF can be obtained from the Master equation of a Markov process by using the Kramers-Moyal expansion for the transition probability (Risken, 1996). If the expansion stops after the second term, one gets the Fokker-Planck equation (also called the forward Kolmogorov equation), reading

[TABLE]

It involves only the drift and diffusion coefficients

[TABLE]

where one should note that these coefficients depend on “${\mathrm{a}}$”, the considered component. In the present context, $\bm{F}_{\rm fric}$ represents the friction force and $\bm{D}$ is the diffusion matrix in action space. In general, for complex systems, it is not possible to determine the Kramers-Moyal coefficients $\bm{D}^{(n)}$ from first principles. However, in the present case, this could be achieved, in the two previous sections, by considering an expansion of the equations of Hamiltonian dynamics in powers of ${1/N}$ in the limit ${N\!\rightarrow\!+\infty}$.666It can be shown that the quasilinear approximation amounts to neglecting terms of order ${1/N^{2}}$ or smaller (Chavanis, 2012c). At order ${1/N}$, the diffusion and drift coefficients were obtained in equations (46) and (50), and read

[TABLE]

where we note that the sums ${\sum_{\bm{m}_{1},\bm{m}_{2}}^{\prime}}$ are restricted to non-zero values of $\bm{m}_{1}$ and $\bm{m}_{2}$. It can be shown that the higher order Kramers-Moyal coefficients are negligible at order ${1/N}$. This fully justifies the Fokker-Planck equation (51).

Using the notations from equation (52), the Fokker-Planck equation (51) can be rewritten as

[TABLE]

or, equivalently, as

[TABLE]

In order to make the connection with the Balescu-Lenard equation, let us rewrite equation (56) under a form in which the diffusion coefficient is “sandwiched” between the two action derivatives, i.e.

[TABLE]

where we defined

[TABLE]

Here, $\bm{F}_{\rm pol}$ represents the friction force by polarisation (Chavanis, 2012c). It differs from the true friction force $\bm{F}_{\rm fric}$ by a term involving the derivatives of the diffusion tensor $\bm{D}$. Integrating equation (54) by parts and comparing the resulting expression with equations (53) and (58), we finally get

[TABLE]

The friction force by polarisation can also be obtained from a direct calculation based on a linear response theory (see Appendix B). Substituting equations (53) and (59) into equation (57), we immediately obtain

[TABLE]

This is the inhomogeneous Balescu-Lenard equation (Heyvaerts, 2010; Chavanis, 2012c). Hence we have demonstrated that the Balescu-Lenard equation is equivalent to the Fokker-Planck equation. These are the kinetic equations describing the secular evolution of dressed particles in inhomogeneous systems via resonant binary interactions.

It is straightforward to specialise the previous expressions to the case where collective effects are not accounted for. This amounts to replacing the dressed potential perturbation ${\delta U_{\mathrm{d}}}$ from equation (26) by the bare potential perturbation ${\delta U_{\rm bare}}$ from equation (27), while all the following calculations remain the same. The dressed susceptibility coefficients ${1/\mathcal{D}_{\bm{m}_{1},\bm{m}_{2}}}$ from equation (30) then become the bare susceptibility coefficients ${1/\mathcal{D}_{\bm{m}_{1},\bm{m}_{2}}^{\rm bare}}$ reading

[TABLE]

where the bare susceptibility coefficients ${A_{\bm{m}_{1},\bm{m}_{1}}(\bm{J}_{1},\bm{J}_{2})}$ (Lynden-Bell, 1994; Pichon, 1994; Chavanis, 2013) are given by the Fourier transform in angles of the interaction potential $u$, so that

[TABLE]

The detailed calculations leading to the third equality of equation (61) are given in Appendix B of Fouvry et al. (2015b). Because of these strong similarities, the bare analogs of the drift and diffusion coefficients from equations (53), (54) and (59) are immediately given by

[TABLE]

Similarly, the inhomogeneous Balescu-Lenard equation (60), when neglecting collective effects, becomes the inhomogeneous Landau equation (Chavanis, 2013) reading

[TABLE]

One can easily check that, in the single species case, all the results presented in this section agree with those obtained in Chavanis (2012c) via a different method.

6 From Fokker-Planck to Langevin

The Kramers-Moyal coefficients appearing in the Fokker-Planck equation (51) may be derived from stochastic Langevin equations (Risken, 1996). In the present context, such a Langevin equation describes the evolution of the action ${\bm{{J}}(t)}$ of a given (test) star. Let us consider a general Langevin equation of the form

[TABLE]

where $\bm{h}(\bm{{J}},t)$ is a vector, $\bm{g}(\bm{{J}},t)$ is a tensor, and ${\bm{\Gamma}(t)}$ is a Gaussian white noise (Langevin force) whose statistics satisfy

[TABLE]

where $\bm{I}$ is the identity matrix. When the tensor $\bm{g}(\bm{{J}},t)$ explicitly depends on the action $\bm{{J}}$ of the particle, we say that the noise is multiplicative. Using the Stratonovich picture (ordinarily used by physicists), the drift and diffusion coefficients are given by

[TABLE]

The other Kramers-Moyal coefficients are zero. The last term in the expression of $\bm{D}^{(1)}$ is the noise-induced drift or spurious drift. The drift and diffusion coefficients determine the Fokker-Planck equation (55) which describes the evolution of the probability density. The drift and diffusion coefficients $\bm{D}^{(1)}$ and $\bm{D}^{(2)}$ are uniquely determined by the functions $\bm{h}$ and $\bm{g}$ of the Langevin equations as given by equation (69).

Let us now consider the inverse problem, i.e., the determination of the Langevin equations from the Fokker-Planck equation. As discussed in Risken (1996), in the multidimensional case, the functions $\bm{h}$ and $\bm{g}$ are not uniquely determined by the drift and diffusion coefficients $\bm{D}^{(1)}$ and $\bm{D}^{(2)}$. One particular solution obtained by diagonalising the positive definite matrix $\bm{D}^{(2)}$ is given by

[TABLE]

The general solution can then be obtained by multiplying the matrix \big{(}\sqrt{{D}}^{(2)}\big{)}_{ij} with arbitrary orthogonal matrices. However, the expression of equation (70) is sufficient for our purposes.

It is not easy to numerically solve the Balescu-Lenard equation (60). However, since we have established that the Balescu-Lenard equation is equivalent to the Fokker-Planck equation, it may be more convenient to solve numerically the stochastic Langevin equations for each individual star (characterised by its action $\bm{{J}}$) and make an ensemble average to reconstruct the system’s DF. In the theory of Brownian motion, this is the so-called molecular dynamics method. The main idea is to simulate the Langevin force on a computer, integrate the equations of motion with the simulated Langevin force and then take the average for a large number of realisations (Risken, 1996). Of course, in the present context, the diffusion and friction coefficients $\bm{D}^{(1)}$ and $\bm{D}^{(2)}$, and therefore the Langevin coefficients $\bm{h}$ and $\bm{g}$, which describe the self-induced noisy environment, must be updated self-consistently as the system’s DF ${F(\bm{J},t)}$ changes on secular timescales. Even if we are led back to a discrete ${N-}$particles system (recall that we started from a Hamiltonian system of $N$ stars), the gain of the stochastic approach is to allow for a time discretisation of the particles’ trajectories with a timestep ${\Delta t}$ that is orders of magnitude larger than the timestep required to solve the Hamiltonian dynamics, since the complicated effects of collisions are encapsulated in the stochastic force and in the drift.

In addition, this method may be more flexible for generalisation. For example, for non integrable systems, equation (67) could be extended to account for chaotic stochasticity as

[TABLE]

where the stochastic Langevin force, ${\bm{\Gamma}_{\mathrm{c}}(t)}$, and its action varying amplitude, $\bm{g}_{\mathrm{c}}$, are set to match the orbital diffusion induced by the chaotic sea within action space.

7 The different stages in the evolution of stellar systems

We are now in a position to describe accurately the different stages that occur in the evolution of a stellar system. Fundamentally, a stellar system is a Hamiltonian system of $N$ stars in gravitational interaction. If we are interested in the evolution of the DF, ${F\!=\!F(\bm{x},\bm{v},t)}$, one can identify different dynamical regimes, each of them characterised by a different kinetic equation:

(i) For sufficiently “short” times (that can be astronomical in practice!), the evolution of the DF of a stellar system is governed by the Vlasov-Poisson equations (Jeans, 1915; Vlasov, 1938). The Vlasov equation is a mean field equation which describes the “collisionless” evolution of the system. Mathematically speaking, it is valid in the limit ${N\!\rightarrow\!+\infty}$ with ${\mu\!\sim\!1/N}$.777For self-gravitating systems, there are mathematical difficulties to rigorously justify the Vlasov equation because of the $r^{-1}$ divergence of the gravitational potential as ${r\!\rightarrow\!0}$. A stellar system described by the Vlasov equation that is initially in an unsteady state, or in a dynamically unstable steady state, generically undergoes a process of violent relaxation  (Lynden-Bell, 1967) and reaches a quasi-stationary (virialised) state on a coarse-grained scale. This process takes place in a few dynamical times $t_{\mathrm{d}}$. Violent relaxation is a complex process associated with large potential fluctuations, phase mixing and nonlinear Landau damping (Mouhot & Villani, 2011). The quasi-stationary state resulting from violent relaxation is difficult to predict in general. However, it must be (close to) a stable steady state of the Vlasov-Poisson equations, as found in numerical simulations. According to Jeans’ theorem, the DF of a stellar system trapped in a quasi-stationary state is generically a function of the actions only, ${F\!=\!F(\bm{J})}$.

(ii) On a secular timescale ${\!\sim\!Nt_{\mathrm{d}}}$, gravitational encounters between stars (departures from the mean field dynamics, granularities, graininess, finite${-N}$ effects,…) come into play and must be taken into account in the dynamics. Because of gravitational encounters (via resonances), the system’s DF slowly changes by evolving through a succession of quasi-stationary states, ${F\!=\!F(\bm{J},t)}$. This self-induced “collisional” evolution of the DF is described by the inhomogeneous Balescu-Lenard or Fokker-Planck equation (Heyvaerts, 2010; Chavanis, 2012c), which is a refinement of the homogeneous Chandrasekhar (Chandrasekhar, 1942) and Landau (Landau, 1936) equations, taking into account spatial inhomogeneity and collective effects. We note that gravitational encounters between stars need not be local (Lynden-Bell & Kalnajs, 1972) but can be distant, e.g. capturing the mechanism of resonant relaxation (Rauch & Tremaine, 1996). An alternative description of the system’s dynamics can be written in terms of ${N-}$body stochastic Langevin equations associated with the inhomogeneous Fokker-Planck equation (see section 6).

(iii) It may happen that, during the collisional evolution, the system’s DF, ${F\!=\!F(\bm{J},t)}$, becomes dynamically (Vlasov) unstable. In that case, one has to come back to the Vlasov-Poisson equations to describe its evolution. This drives a dynamical phase transition from this unstable state to a new stable state. This has been found in Sellwood (2012) and explicity demonstrated in Fouvry et al. (2015c) in the case of stellar discs, showing the transition between a disc-phase (axisymmetric) and a bar-phase (non-axisymmetric). A similar dynamical phase transition was evidenced previously for a toy model of particles with long-range interactions (Campa et al., 2008).

(iv) It can be shown that the (inhomogeneous) Balescu-Lenard equation conserves mass and energy and satisfies a ${H-}$theorem for the Boltzmann entropy (Chavanis, 2007; Heyvaerts, 2010), see Appendix C. As a result, one expects that the DF relaxes for ${t\!\rightarrow\!+\infty}$ towards the Boltzmann DF which maximises the entropy at fixed mass and energy. However, for self-gravitating systems, in most cases, the Boltzmann entropy has no maximum (e.g. Padmanabhan, 1990), so that there exists no statistical equilibrium state in a strict sense. For example, the late time evolution of globular clusters proceeds through stellar evaporation (Spitzer, 1940). According to the virial theorem, the central density increases as the system expands. When the system becomes sufficiently centrally condensed, an instability develops and leads to core collapse. This instability, called the gravothermal catastrophe (Lynden-Bell & Wood, 1968), arises from the negative specific heat of the inner part of the cluster. Core collapse can be stopped by the formation of a binary star that can release an enormous amount of energy able to reverse the collapse and drive a re-expansion of the whole cluster until the next collapse takes place. This can lead to a series of gravothermal oscillations (Bettwieser & Sugimoto, 1984).

(v) Finally, even in the collisionless regime ${N\!\rightarrow\!+\infty}$, the DF may evolve under the effect of external perturbations, again passing through a succession of quasi-stationary states, ${F\!=\!F(\bm{J},t)}$. The kinetic equation that governes this dynamics is the secular collisionless diffusion equation introduced in Binney & Lacey (1988); Weinberg (2001); Pichon & Aubert (2006); Fouvry et al. (2015a) for inhomogeneous systems and in Nardini et al. (2012); Chavanis (2012b) for homogeneous systems. We note that, contrary to the Balescu-Lenard equation which has no free parameter, this equation needs an input which is the power spectrum of the external potential fluctuations. The effect of the external environment may or may not outrun that of the self-induced evolution: this is the classical conundrum of “nature” and “nurture” driven secular evolutions.

8 Conclusion

The derivation of the coefficients of diffusion and dynamical friction in a stable, inhomogeneous, multicomponent, self-gravitating system was presented. The method followed here is based on the detailed study of the dynamics of a test particle, when perturbed by the dressed potential perturbations induced by a discrete bath of background particles. It was shown in particular how the averaged coefficients of diffusion and dynamical friction are fully consistent with those involved in the associated inhomogeneous Balescu-Lenard equation. As a result, the Balescu-Lenard equation can be interpreted as a Fokker-Planck equation in which the diffusion and friction coefficients evolve self-consistently (i.e. they depend on the DF itself).

The present derivation has several advantages. First of all, it clarifies the physical content of the Balescu-Lenard equation by showing its equivalence with the traditional Fokker-Planck equation that was introduced initially in the seminal work of Chandrasekhar (1943a) and that has been adopted by most astrophysicists. This approach confirms that the force acting on a star can be decomposed into a smooth component due to the mean field of the whole system and fluctuations due to finite${-N}$ effects (encounters). In turn, the fluctuations have a completely random part that can be described by a multiplicative Gaussian white noise and a systematic part corresponding to the effect of dynamical friction. This is in complete agreement with the physical picture given by Chandrasekhar (1943a). However, considerable progress has been made in the calculation of the coefficients of diffusion and dynamical friction with respect to early approaches that focused on spatially homogeneous stellar systems (making a local approximation) and neglected collective effects (the dressing of a star by its gravitational wake). We are now in a position to account for spatial inhomogeneity and collective effects accurately. This results in self-consistent expressions of the coefficients of diffusion and dynamical friction (given by equations (53) and (54)) at order ${1/N}$, that encompass previous results obtained in the literature.888Actually, there can remain a logarithmic divergence at small scales (in particular for ${3D}$ spherical systems) due to strong collisions. This divergence can be solved by developing a treatment “à la Boltzmann” or by taking into account the bending of trajectories like in the work of Chandrasekhar (1943a). Alternatively, one may introduce in the diverging expression a cut-off at the Landau length which corresponds to the impact parameter producing a deflection at $90^{\circ}$ of the particle’s trajectory. There is no such divergence in the case of stellar discs. We refer the reader to Chavanis (2013) for a thorough and detailed discussion of the links between self-consistent kinetic equations such as the Balescu-Lenard and Landau equations, and other approaches, such as the two-body encounters theory introduced in Chandrasekhar (1943a). The present formalism also allows for a self-consistent description of a spectrum of masses, with a proper accounting of the induced secular mass segregation, which should be of interest to various astrophysical contexts, from galactic centers to protostellar discs. Another advantage of the present derivation is practical. Instead of numerically solving the Balescu-Lenard equation or the Fokker-Planck equation, it may be more convenient to solve a system of $N$ Langevin equations describing the stochastic trajectories of stars on an intermediate timescale with the coefficients of diffusion and dynamical friction obtained from the Fokker-Planck approach. This procedure may be useful in stellar dynamical simulations since it allows one to use larger timesteps compared to the ones used in the original ${N-}$body Hamiltonian equations, as the encounters between stars have been taken into account in the kinetic parametrisation. Yet, computing the diffusion flux of such kinetic equations remains a challenge. In order to deal with the system’s inhomogeneity, one has to construct a set of angle-action coordinates ${(\bm{\theta},\bm{J})}$. In order to characterise the self-gravitating amplification, one may rely on the matrix method to construct a biorthogonal basis of potential and density elements ${(\psi^{(p)},\rho^{(p)})}$ and estimate the system’s global response matrix ${\widehat{\mathbf{M}}(\omega)}$. Finally, the secular evolution being driven by resonant encounters, one has to solve the non-local resonance condition, ${\delta_{\mathrm{D}}(\bm{m}_{1}\!\cdot\!\bm{\Omega}_{1}\!-\!\bm{m}_{2}\!\cdot\!\bm{\Omega}_{2})}$, present in the Balescu-Lenard equation (60). See, e.g., Fouvry et al. (2015c) for an illustration of how these various difficulties may be solved in the context of razor-thin axisymmetric stellar discs.

Although the formalism and discussions are presented in the context of self-gravitating systems, this approach is actually valid for arbitrary systems with long-range interactions in any dimension of space. One just has to introduce a proper biorthogonal basis of potentials and densities associated with the specific interaction potential, as defined by equation (5), and introduce the set of angle-action variables associated with the unperturbed Hamiltonian $H_{0}$. In addition, there are strong analogies between two-dimensional point vortices and stellar systems (see, e.g. Chavanis, 2002). In the same manner that a test star in a star cluster has a diffusion motion in velocity space due to the fluctuations of the gravitational force and experiences a dynamical friction due to a polarisation process, a point vortex evolving in a sea of field vortices has a diffusion motion in position space due to the fluctuations of the velocity field and experiences a systematic drift (Chavanis, 1998) due to a polarisation process. The evolution of the probability density of its position is governed by a Fokker-Planck equation that can be written in the form of a Balescu-Lenard equation (Chavanis, 2012a) in complete parallel with the Fokker-Planck and Balescu-Lenard equations of stellar systems. In the thermal bath approach, the friction and drift coefficients are related to the diffusion coefficients by a form of Einstein relation expressing the fluctuation-dissipation theorem.

Acknowledgements

This paper was initially drafted by Jean Heyvaerts in 2010, independently of Chavanis (2012c), and completed posthumously by his co-authors to honour his memory. JBF thanks Scott Tremaine for insightful comments. Support for Program number HST-HF2-51374 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. This research is part of ANR grant Spin(e) (ANR-13-BS05-0005, http://cosmicorigin.org).

Appendix A Computing the drift vector

In this Appendix, we compute the mean drag {\big{<}\Delta\bm{J}_{1}\big{>}} acting on particle $1$ during the time ${\Delta t}$. Using the shorthand notation introduced in equation (32), equation (49) can be rewritten as

[TABLE]

where in the last equation, the gradient w.r.t. $\bm{J}_{2}$ also includes the dependency of ${\omega_{2}\!=\!\bm{m}_{2}\!\cdot\!\bm{\Omega}_{2}(\bm{J}_{2})}$. The solution to the previous system must be sought to second order in the noise, which is given by the number of $\Lambda$ factors. The change in ${\Delta\bm{J}_{1}}$ during ${\Delta t}$ is formally given through integration of equation (72) between [math] and ${\Delta t}$, so that

[TABLE]

where one must note the time dependence of ${\Lambda(\,\cdot\,;t)}$ which has to be accounted for at this order in the noise.

Progress can be made towards solution accurate to second order by computing ${\Lambda_{\bm{m}_{1},\bm{m}_{2}}(1,2,\omega_{2};t)}$, after substituting the expressions for $\dot{\bm{J}}_{1}$ and $\dot{\bm{J}}_{2}$ given by equation (72). Then equation (74) gives

[TABLE]

Via time integration, we may then obtain an expression for ${\Lambda_{\bm{m}_{1},\bm{m}_{2}}(1,2,\omega_{2};t)}$ which is explicitly second order (recalling that $\Lambda$ is already first order). At this order, all the involved elements can be evaluated to zeroth order, so that in the r.h.s. of equation (76), all the occurences of $\Lambda$ can be evaluated for ${t\!=\!0}$. In the upcoming calculations, this is no more explicitly written to simplify the notations. The time integration of equation (76) gives

[TABLE]

where we insist on the fact that the occurences of ${\Lambda}$ in the r.h.s. of equation (77) are evaluated at zeroth order, i.e. for ${t\!=\!0}$. Let us now obtain ${z_{12}(t_{1})}$, which enters equation (75), thanks to the time integration of equation (73). One gets

[TABLE]

where ${\Lambda}$ is evaluated at zeroth-order, i.e. for ${t\!=\!0}$. The frequencies ${\bm{\Omega}_{1}(t_{2})}$ and ${\bm{\Omega}_{2}(t_{2})}$ follow by time integration of equation (49) and read

[TABLE]

Then, the relative angular velocity $g_{12}$ at time $t_{2}$ reads

[TABLE]

Equation (78) finally becomes

[TABLE]

The lines 2 and 3 of equation (81) correspond to the first order correction relative to the zeroth expression on the first line. Hence the complex exponential ${\mathrm{e}^{\mathrm{i}z_{12}(t_{1})}}$ can be expanded as

[TABLE]

where we relied on the usual development ${\mathrm{e}^{\mathrm{i}\varepsilon}\!\simeq\!1\!+\!\mathrm{i}\varepsilon}$ at first order in $\varepsilon$. In equation (82), the expressions of ${\dot{\bm{J}}_{1}(t_{3})}$ and ${\dot{\bm{J}}_{2}(t_{3})}$ are given by Hamilton’s equation (72). The expressions (77) and (82) are then replaced in equation (75). It yields the expression of ${\Delta\bm{J}_{1}/\Delta t}$ which now reaches the required second order level, so that

[TABLE]

Let us first insist on the fact that equation (83) is explicitly second order in the noise. Indeed, the only linear term in equation (83), arising from ${\mathrm{e}^{\mathrm{i}z_{12}(0)}\Lambda_{\bm{m}_{1},\bm{m}_{2}}(1,2,\omega_{2})}$ will vanish when averaged over the initial angle $\bm{\theta}_{2}^{0}$. We also recall that at the present level of approximation, the complex exponentials $\mathrm{e}^{\mathrm{i}z(t)}$ should be evaluated to zeroth order, that corresponding to the uniform angular motion at fixed frequency. For example, in equation (83), one should read ${\mathrm{e}^{\mathrm{i}z_{1^{\prime}3^{\prime}}(t_{2})}}$ as

[TABLE]

Similarly, at the order considered here, the susceptibility coefficients ${\Lambda(...;t)}$ should be evaluated at ${t\!=\!0}$.

One may then follow the same method as the one presented in equation (40) to obtain the averaged action diffusion tensor. Indeed, one can average equation (83) over the initial angles of the particles $1$, $2$, $3$, and $4$, as well as on the action distribution of particles $2$, $3$, and $4$. Here, one should pay attention to the fact that particle $1$ acts as our test star, while particles $2$ and $3$ both run over the field stars associated with particle $1$, i.e. over all stars except particle $1$. The situation is slightly different for particle $4$ which runs over the field stars associated with particle $2$, i.e. over all stars except particle $2$. Let us first perform an average over the initial angles of all stars. Considering only second order terms (i.e. involving two factors $\Lambda$) and keeping only the dependencies w.r.t. the initial angles, equation (83) requires to study two different generic terms

[TABLE]

where it is important to note that the sums on particles $2$ and $3$ are restricted to all stars except particle $1$, while the sum on particle $4$ is restricted to all stars except particle $2$. Because only non-zero values of $\bm{m}_{2}$, ${\bm{m}_{3}^{\prime}}$, and ${\bm{m}_{4}^{\prime}}$ contribute to the fluctuations, equation (85) therefore immediately imposes for particles $3$ and $2$ to be the same and for particles $4$ and $1$ to be the same. As a consequence, the sums $\sum_{3}$ and $\sum_{4}$ can be straightforwarldy executed. We may then average equation (83) over $\bm{\theta}_{1}^{0}$, $\bm{\theta}_{2}^{0}$, and over the action distribution of particle $2$, following the same substitution as in equation (42). All in all, equation (83), when averaged and restricted to second order terms, becomes

[TABLE]

where one should pay attention to the different mass prefactors ${\mu_{\mathrm{a}}\!=\!\mu_{1}}$ and $\mu_{\mathrm{b}}$, whose role is essential to induce mass segregation in multicomponent systems. In equation (86), we also performed the change of notations ${\bm{m}_{3}^{\prime}\!\to\!\bm{m}_{2}^{\prime}}$ and ${\bm{m}_{4}^{\prime}\!\to\!\bm{m}_{1}^{\prime}}$. The averaging process over the initial angles $\bm{\theta}_{1}^{0}$ and $\bm{\theta}_{2}^{0}$ yields

[TABLE]

where we recall that ${g_{12}\!\equiv\!\bm{m}_{1}\!\cdot\!\bm{\Omega}_{1}\!-\!\bm{m}_{2}\!\cdot\!\bm{\Omega}_{2}}$, and is evaluated at zeroth order. Equation (86) then becomes

[TABLE]

The various double and triple time integrals of complex exponentials occuring in equation (88) can be replaced by distributions in the limit ${\Delta t\!\to\!+\infty}$. They are computed in Appendix E, and one has

[TABLE]

These relations allow us to rewrite equation (88) as

[TABLE]

One can straightforwardly show that for any function ${G(g_{12})}$, one has

[TABLE]

Applying these formulae to equation (90), one can collect quite a few terms in this equation, to get

[TABLE]

The lack of symmetry of this equation might look troublesome. In fact at the second order level, one can identify ${\Lambda_{\bm{m}_{2},\bm{m}_{1}}(2,1,\omega_{1})}$ with the complex conjugate of ${\Lambda_{\bm{m}_{2},\bm{m}_{1}}(2,1,\omega_{1})}$, as demonstrated in Appendix F. Similarly, ${\Lambda_{-\bm{m}_{1},-\bm{m}_{2}}(2,1,-\omega_{2})}$ is the complex conjugate of ${\Lambda_{\bm{m}_{1},\bm{m}_{2}}(1,2,\omega_{2})}$ (see Appendix F). Equation (92) can therefore be rewritten as

[TABLE]

Appendix B Direct calculation of the friction force by polarisation

The friction force $\bm{F}_{\rm pol}$ that appears in the Fokker-Planck equation (57) is called the “friction by polarisation” (Chavanis, 2013). This is just one component of the total friction force $\bm{F}_{\rm fric}$ that appears in the Fokker-Planck equation (56). Physically, $\bm{F}_{\rm pol}$ is the force resulting directly from the retroaction of the field stars to the perturbation caused by the test star, like in a polarisation process. Some particular, or formal, expressions of the friction force by polarisation have been derived in Marochnik (1968); Kalnajs (1971); Kandrup (1983); Bekenstein & Zamir (1990); Chavanis (2008) from a linear response theory based on the Liouville equation or on the Klimontovich equation. Let us derive here its general expression in angle-action variables taking spatial inhomogeneity and collective effects into account and check that its expression agrees with equation (59). This calculation is inspired by section 3.3 of Chavanis (2012c). The perturbation induced by the test particle is determined by the coefficients $\widehat{a}^{\alpha}_{\mathrm{p}}(\omega)$. Following equation (24), they can be rewritten as

[TABLE]

The temporal Fourier transform of the potential perturbation ${\delta U_{\mathrm{p}}}$ associated with the gravitational wake reads

[TABLE]

where the coefficients ${\widehat{a}^{\beta}_{\mathrm{t}}(\omega)}$ of the test particle potential are given by equation (17). Therefore

[TABLE]

Written as a function of time, equation (96) becomes

[TABLE]

Replacing ${\psi^{(\alpha)}(\bm{x})}$ by the sum

[TABLE]

we obtain

[TABLE]

This formula completely specifies the polarisation cloud. The corresponding force in action, ${\delta\bm{F}_{\mathrm{p}}(\bm{x},t)\!=\!-\partial_{\bm{\theta}}\delta U_{\mathrm{p}}(\bm{x},t)}$, is

[TABLE]

One may then evaluate this force at the location of the test particle, and average it over the orbit of the test particle. One obtains the friction force by polarisation, ${\bm{F}_{\rm pol}\!=\!\!\int\!\mathrm{d}\bm{\theta}_{\mathrm{t}}/(2\pi)^{d}\delta F_{\mathrm{p}}(\bm{x}(\bm{\theta}_{\mathrm{t}},\bm{J}_{\mathrm{t}}))}$, reading

[TABLE]

Introducing the system’s dressed and bare susceptibility coefficients ${1/\mathcal{D}_{\bm{m}_{1},\bm{m}_{2}}(\bm{J}_{1},\bm{J}_{2},\omega)}$ and ${1/\mathcal{D}_{\bm{m}_{1},\bm{m}_{2}}^{\rm bare}(\bm{J}_{1},\bm{J}_{2})}\!=\!-{A_{\bm{m}_{1},\bm{m}_{1}}(\bm{J}_{1},\bm{J}_{2})}$ from equations (30) and (61), one gets

[TABLE]

Finally, using identity (54) of Chavanis (2012c), i.e.

[TABLE]

we obtain

[TABLE]

which coincides with equation (59). If we neglect collective effects, we recover equation (65). The fact that $\bm{F}_{\rm pol}$ is just one component of the true friction force $\bm{F}_{\rm fric}$ is clear from equation (58) (see also section 3 of Chavanis (2012c)).

Note that this calculation remains valid if the test particle is of different nature from the field particles. In particular, it could be a satellite of mass $\mu_{\mathrm{t}}$ moving in a collisionless fluid of stars with mass ${\mu_{\mathrm{b}}\!\sim\!1/N_{\mathrm{b}}}$ governed by the Vlasov equation when ${N_{\mathrm{b}}\!\rightarrow\!+\infty}$. In the limit ${\mu_{\mathrm{t}}\!\gg\!\mu_{\mathrm{b}}}$, the friction by polarisation is the only force that acts on the test particle, and the satellite sinks at the center of the system (see Appendix D.2).

Appendix C Properties of the multicomponent Balescu-Lenard equation

In this Appendix, we derive the main properties of the inhomogeneous Balescu-Lenard equation (60) following the works of Chavanis (2007); Heyvaerts (2010); Chavanis (2012c) and references therein. This extends their results to the important class of astrophysical systems containing different components.

C.1 Conservative form

The Balescu-Lenard equation (60) can be written in the conservative form

[TABLE]

where

[TABLE]

is the flux of particles of component “$\mathrm{a}$”. Under that form, it is immediately clear that the Balescu-Lenard equation conserves the total mass of each component.

C.2 Equilibrium state: Boltzmann distribution

It is straightforward to check that the Boltzmann distribution

[TABLE]

where ${\epsilon(\bm{J})}$ is the energy of a star by unit of mass, ${\beta\!=\!1/T}$ is the inverse “temperature” and $A_{\mathrm{a}}$ is a normalisation constant, is a steady state of the Balescu-Lenard equation. Indeed, using

[TABLE]

we get

[TABLE]

Substituting this relation into the Balescu-Lenard flux from equation (106), we find

[TABLE]

The integrand involves the term ${\delta_{\mathrm{D}}(\bm{m}_{1}\!\cdot\!\bm{\Omega}_{1}\!-\!\bm{m}_{2}\!\cdot\!\bm{\Omega}_{2})(\bm{m}_{1}\!\cdot\!\bm{\Omega}_{1}\!-\!\bm{m}_{2}\!\cdot\!\bm{\Omega}_{2})}$, which is obviously equal to zero, so that one has ${\partial_{t}F^{\mathrm{a}}\!=\!0}$ for each component. Note that the temperature in the Boltzmann distribution of equation (107) is the same for all the components. This corresponds to an equipartition of energy, and usually implies that heavy particles sink at the center of the system while light particles wander around.

The Boltzmann distribution can be obtained by maximising the Boltzmann entropy $S_{\mathrm{B}}$ defined by equation (120) while conserving the total energy $E$ and the total mass $M_{\mathrm{a}}$ of each species of particles. Writing the variational principle as

[TABLE]

where $\beta$ (inverse temperature) and $\alpha_{\mathrm{a}}$ (chemical potentials) are the Lagrange multipliers associated with energy and mass conservation, we obtain

[TABLE]

Since this condition must be satisfied for arbitrary variations ${\delta F^{\mathrm{a}}}$ the term in braces must vanish, leading to the Boltzmann distribution (107). The Boltzmann distribution is therefore a critical point of entropy at fixed mass and energy. However, let us stress once again that a statistical equilibrium state does not always exist for self-gravitating systems (notably for ${3D}$ spherical systems). The Boltzmann entropy may not have a (global or local) maximum. Even worse, the Boltzmann distribution from equation (107) may not be normalisable (i.e. the Boltzmann entropy may not have any critical point).

C.3 Energy conservation

The total energy of the system is

[TABLE]

Taking its time derivative and using equation (105), we get

[TABLE]

Integrating by parts, assuming that boundary terms do not contribute, and using equation (108), we obtain

[TABLE]

Substituting equation (106) into equation (115), we get

[TABLE]

Interchanging the dummy variables ${(\mathrm{a},\mathrm{b})}$, ${(\bm{m}_{1},\bm{m}_{2})}$ and ${(\bm{J}_{1},\bm{J}_{2})}$, and using the property ${\mathcal{D}_{\bm{m}_{2},\bm{m}_{1}}(\bm{J}_{2},\bm{J}_{1},\omega)\!=\!\mathcal{D}_{\bm{m}_{1},\bm{m}_{2}}(\bm{J}_{1},\bm{J}_{2},\omega)^{*}}$ (see Appendix F.2), we obtain

[TABLE]

Taking the half-sum of these equations, we get

[TABLE]

As in equation (110), the integrand involves the term ${\delta_{\mathrm{D}}(\bm{m}_{1}\!\cdot\!\bm{\Omega}_{1}\!-\!\bm{m_{2}}\!\cdot\!\bm{\Omega}_{2})(\bm{m}_{1}\!\cdot\!\bm{\Omega}_{1}\!-\!\bm{m}_{2}\!\cdot\!\bm{\Omega}_{2})}$, which is identically zero. One therefore has ${\dot{E}\!=\!0}$, the total energy of the system is conserved. Note, however, that the energy of each species is not individually conserved. Using equations (106) and (115), one has

[TABLE]

C.4 ${H-}$theorem

The multicomponent Boltzmann entropy is

[TABLE]

Taking its time derivative, and using equation (105), we get

[TABLE]

Integrating by parts and assuming that boundary terms do not contribute, we obtain

[TABLE]

Substituting equation (106) into equation (122), we get

[TABLE]

Interchanging the dummy variables ${(\mathrm{a},\mathrm{b})}$, ${(\bm{m}_{1},\bm{m}_{2})}$ and ${(\bm{J}_{1},\bm{J}_{2})}$, and using the property $\mathcal{D}_{\bm{m}_{2},\bm{m}_{1}}(\bm{J}_{2},\bm{J}_{1},\omega)=\mathcal{D}_{\bm{m}_{1},\bm{m}_{2}}(\bm{J}_{1},\bm{J}_{2},\omega)^{*}$ (see Appendix F.2), we obtain

[TABLE]

Taking the half-sum of these equations, we get

[TABLE]

from which we obtain the ${H-}$theorem

[TABLE]

It establishes that the Boltzmann entropy is monotonically increasing. Note however that the Boltzmann entropy of each species does not individually satisfy a ${H-}$theorem. Using equations (106) and (122), one has

[TABLE]

For neutral plasmas, that are spatially homogeneous, the Boltzmann entropy is bounded from above. In that case, one can show from the ${H-}$theorem that the homogeneous Balescu-Lenard relaxes, for ${t\!\rightarrow\!+\infty}$, towards the Boltzmann distribution (107). This is the maximum entropy state at fixed mass and energy. We note that the Balescu-Lenard equation singles out the Boltzmann distribution among all possible steady states of the Vlasov equation. For self-gravitating systems, that are spatially inhomogeneous, the Boltzmann entropy is typically not bounded from above (this is notably the case for ${3D}$ spherical or ${2D}$ flat systems). In that case, the inhomogeneous Balescu-Lenard equation does not relax towards an equilibrium state. It can describe stellar evaporation and core collapse (gravothermal catastrophe) as discussed in section 7. It cannot, however, account for the formation of binary stars and for gravothermal oscillations since the formation of binaries results from three-body collisions that are neglected in the Balescu-Lenard treatment. For other systems with long-range interactions for which the Boltzmann entropy is bounded from above, the Boltzmann distribution is always a steady state of the inhomogeneous Balescu-Lenard equation (see section C.2). However, the inhomogeneous Balescu-Lenard equation does not necessarily relax towards this distribution. In very specific situations, the system may remain blocked in another state if the resonance condition cannot be fulfilled. This “kinetic blocking” is illustrated in the case of ${2D}$ point vortices in Chavanis & Lemou (2007) when the profile of angular velocity is monotonic. In that case, no resonance is possible whatever the distribution function and the system remains frozen in a distribution that is different from the Boltzmann distribution. Only higher order correlations (three-body, four-body…) may unblock the system. Such correlations have their own higher order (e.g. ${1/N^{2}}$) kinetic equations.

C.5 Initial flux in non-thermalised mix

Let us consider an initial condition in which all the particles have a Boltzmann distribution of the form

[TABLE]

However, we do not assume equipartition of energy: the temperature ${T_{\mathrm{a}}\!=\!1/\beta_{\mathrm{a}}}$ of each species may be different. This may correspond to the situation in a stellar disc where old and young populations of stars co-exist. The collisions, described by the Balescu-Lenard equation, change the energy of each component so that their temperature is the same at equilibrium (equipartition of energy). The initial value of the flux defined by equation (106) is

[TABLE]

It can be written as

[TABLE]

where $\bm{D}^{\mathrm{a}\mathrm{b}}$ is the diffusion tensor of component “$\mathrm{a}$” caused by the collisions with the particles of component “$\mathrm{b}$”. The total diffusion tensor defined by equation (53) can be written as ${\bm{D}\!=\!\sum_{\mathrm{b}}\bm{D}^{\mathrm{a}\mathrm{b}}}$. The initial rate of change of energy, obtained from equations (115) and (130) is given by

[TABLE]

It can be written more explicitly as

[TABLE]

The rate of change of energy of species “$\mathrm{a}$” scales like the mass weighted sum of inverse temperature differences. We can similarly compute the initial rate of change of entropy from equations (122) and (130). We find that

[TABLE]

Appendix D Test particle approach

Different interpretations can be given to the Balescu-Lenard equation. The Balescu-Lenard equation (60) is an integro-differential equation that describes the evolution of an ensemble of particles (e.g. stars) in interaction. In this interpretation, all the particles are treated on the same footing and their distribution function ${F^{\mathrm{a}}(\bm{J}_{1},t)}$ evolves self-consistently according to equation (60). In a second interpretation, one can select a particular test particle with mass $\mu_{\mathrm{t}}$ and study the evolution of its probability density ${P(\bm{J}_{1},t)}$ in a cloud of field particles with a static distribution function ${F^{\mathrm{b}}(\bm{J}_{2})}$. In this interpretation,999As discussed in Appendix F of Chavanis (2013), this interpretation is valid either for a single component system or for a multicomponent system. The test particle may represent just one particle or an ensemble of non-interacting particles of the same species. In all cases, one has to assume that the collisions between the test particle(s) and the field particles do not alter the distribution of the field particles. the Balescu-Lenard equation (60) is transformed into a differential equation

[TABLE]

usually referred to as the Fokker-Planck equation. The diffusion and friction coefficients are given by

[TABLE]

The Fokker-Planck equation (134) can be rewritten as

[TABLE]

The rates of change of energy ${E_{\mathrm{t}}(t)\!=\!\!\int\!\mathrm{d}\bm{J}_{1}\,P(\bm{J}_{1},t)\mu_{\mathrm{t}}\epsilon(\bm{J}_{1})}$ and Boltzmann entropy ${S_{\mathrm{t}}(t)\!=\!-\!\int\!\mathrm{d}\bm{J}_{1}P(\bm{J}_{1},t)\ln P(\bm{J}_{1},t)}$ of the test particle are

[TABLE]

[TABLE]

For the self-consistency of this interpretation, the field stars must be at statistical equilibrium with the Boltzmann distribution (see section D.1), in a “blocked state” (see section C.4), or have a very long relaxation time (see section D.2) so that their distribution function ${F^{\mathrm{b}}(\bm{J}_{2})}$ does not change under the effect of collisions among themselves on the timescale over which the test particle evolves. A possible astrophysical situation of relevance is a hot halo embedding a cold galactic disc. At zeroth order, the halo can be taken to remain unresponsive to its own fluctuations and to that of its disc.

D.1 Thermal bath: Einstein relation and fluctuation-dissipation theorem

In this subsection, let us consider the relaxation of a test star described by a probability density ${P(\bm{J}_{1},t)}$ in a cloud of field stars at statistical equilibrium with the Boltzmann distribution given by equation (107). Using equation (109), we find that the friction by polarisation (137) experienced by the test star becomes

[TABLE]

Using the Dirac delta function to replace ${\bm{m}_{2}\!\cdot\!\bm{\Omega}_{2}}$ by ${\bm{m}_{1}\!\cdot\!\bm{\Omega}_{1}}$ in the last parenthesis, we find that

[TABLE]

where $\bm{D}$ is the diffusion tensor given by equation (135), in which the field stars have the Boltzmann distribution (107). We may then define the friction tensor, $\bm{\xi}$, as

[TABLE]

Equation (143) is the appropriate form of the Einstein relation for this problem. It relates the friction coefficient $\bm{\xi}$ to the diffusion coefficient $\bm{D}$, the inverse temperature ${\beta\!=\!1/T}$, and the mass of the test particle $\mu_{\mathrm{t}}$. This is a manifestation of the fluctuation-dissipation theorem. As previously emphasised (Chavanis, 2012c), the Einstein relation is valid for the friction by polarisation, not for the true friction. The Fokker-Planck equation (134), when written similarly to equation (138), takes the form of a generalised Kramers equation

[TABLE]

At equilibrium, the probability density of the test particle relaxes towards the Boltzmann distribution

[TABLE]

which is the steady state of equation (144). These results generalise those of Chandrasekhar (1943a).

For a thermal bath, the rates of energy and entropy of the test particle are

[TABLE]

[TABLE]

We have the general relation ${\dot{S}_{\rm pol}\!=\!-\beta\dot{E}_{\rm diff}}$ between the rate of entropy due to the polarisation and the rate of energy due to the diffusion. We note that the Fokker-Planck equation (144) does not conserve the energy and does not satisfy a ${H-}$theorem for the Boltzmann entropy contrary to the Balescu-Lenard equation (60). This is because it describes the evolution of a test particle in a thermal bath with a fixed temperature $T$. This corresponds to a canonical description while the Balescu-Lenard equation corresponds to a microcanonical one. We can obtain a form of ${H-}$theorem for the Fokker-Planck equation (144) by introducing the free energy ${F_{\mathrm{t}}\!=\!E_{\mathrm{t}}\!-\!TS_{\mathrm{t}}}$ which is the Legendre transform of the entropy $S_{\mathrm{t}}$ w.r.t. the energy $E_{\mathrm{t}}$, with conjugate parameter $T$. Taking the time derivative of $F_{\mathrm{t}}$, substituting the Fokker-Planck equation (144), and integrating by parts, we obtain the (canonical) ${H-}$theorem

[TABLE]

It establishes that the Boltzmann free energy is monotonically decreasing. At equilibrium, ${\dot{F}_{\mathrm{t}}\!=\!0}$, leading to the Boltzmann distribution from equation (145). The Boltzmann distribution minimises the Boltzmann free energy $F_{\mathrm{t}}$ while accounting for the normalisation condition ${\!\int\!\mathrm{d}\bm{J}_{1}P(\bm{J}_{1})\!=\!1}$. The first variations write ${\delta F_{\mathrm{t}}\!-\!\alpha\,\delta\!\int\!\mathrm{d}\bm{J}_{1}P(\bm{J}_{1})\!=\!0}$, leading to the equilibrium state from equation (145).

Let us now consider an initial condition in which the test particle has a Boltzmann distribution of the form

[TABLE]

where ${T_{\mathrm{t}}\!=\!1/\beta_{\mathrm{t}}}$ may be different from the temperature ${T\!=\!1/\beta}$ of the bath. The initial rate of change of the test particle energy, $E_{\mathrm{t}}$, is given by

[TABLE]

It can be written more explicitly as

[TABLE]

The initial rate of change of entropy is

[TABLE]

The initial rate of change of free energy ${F_{\mathrm{t}}\!=\!E_{\mathrm{t}}\!-\!TS_{\mathrm{t}}}$ can be written as ${\dot{F}_{\mathrm{t}}\!=\!(1\!-\!T/T_{\mathrm{t}})\dot{E}_{\mathrm{t}}}$.

D.2 Sinking satellite

Let us now assume that the test particle has a mass $\mu_{\mathrm{t}}$ much larger than the mass $\mu_{\mathrm{b}}$ of the field particles. More precisely, we assume that ${\mu_{\mathrm{t}}\!\sim\!1}$ while ${\mu_{\mathrm{b}}\!\sim\!1/N_{\mathrm{b}}}$ with ${N_{\mathrm{b}}\!\gg\!1}$. For ${N_{\mathrm{b}}\!\rightarrow\!+\infty}$, ${\mu_{\mathrm{b}}\!\rightarrow\!0}$ and the field particles form a collisionless fluid of stars. Since their distribution function ${F^{\mathrm{b}}(\bm{J}_{2})}$ does not evolve under the effect of collisions (the relaxation time of the field particles scales like ${N_{\mathrm{b}}t_{\mathrm{d}}\!\rightarrow\!+\infty}$), it can have an arbitrary shape provided that it is Vlasov stable. In this fluid limit for which ${N_{\mathrm{b}}\!\rightarrow\!+\infty}$, the diffusion coefficient from equation (135) vanishes and the friction force from equation (136) reduces to the component proportional to $\mu_{\mathrm{t}}$, which corresponds to the friction by polarisation given by equation (137). Finally, the Fokker-Planck equation (134) reduces to

[TABLE]

where ${F\!=\!\sum_{\mathrm{b}}F^{\mathrm{b}}}$ is the total distribution function of the field particles. Equation (153) can be rewritten as

[TABLE]

This is just the continuity equation corresponding to the deterministic (not stochastic) equation of motion

[TABLE]

We note that the mass $\mu_{\mathrm{b}}$ of the field particles does not appear explicitly in this expression since they form a collisionless “fluid of stars” entirely determined by its total distribution function $F$. As a result, the friction force given by equation (155) can be obtained from the calculation of Appendix B, which is solely based on the Vlasov equation (when ${N_{\mathrm{b}}\!\rightarrow\!+\infty}$) without reference to a discrete Hamiltonian system. Equation (155) can be viewed as the correct generalisation of the famous Chandrasekhar’s formula of dynamical friction when the mass of the test particle is much larger than the mass of the field particles (Chavanis, 2013). One can easily compute the average energy lost by the test particle. We have

[TABLE]

Substituting equation (155) into equation (156), we get

[TABLE]

When ${F\!=\!F(\epsilon)}$, we obtain

[TABLE]

For a thermal bath, considering the limit ${\mu_{\mathrm{b}}\!\sim\!1/N_{\mathrm{b}}\!\rightarrow\!0}$ and ${\beta\!\rightarrow\!+\infty}$ such that ${\beta\mu_{\mathrm{b}}\!\sim\!1}$, we may finally write

[TABLE]

This expression is consistent with equation (151) when ${\beta\!\rightarrow\!+\infty}$. It shows the somehow paradoxical resurgence of the diffusion coefficient from equation (135) in the friction term from equation (155) due to the fluctuation-dissipation theorem from equation (142), although diffusion caused by finite${-N}$ effects is neglected (at leading order) in the present approach.

The Balescu-Lenard equation (60), the Fokker-Planck equation (134) and the generalised Kramers equation (144) describe a competition between diffusion and friction. As a result, the probability density of the test particle relaxes towards a statistical equilibrium state in which the two effects balance each other establishing the Boltzmann distribution from equations (107) or (145).101010This is the usual picture of Brownian theory. As we previously indicated, the situation is more complicated for self-gravitating systems since a statistical equilibrium state does not always exist. In the situation described by the deterministic equation (155), the test particle just feels a friction force and sinks towards the center of the system. In astrophysics, this is traditionally referred to as the “sinking satellite” problem. For small but non-zero values of $\mu_{\mathrm{b}}$, this result is consistent with the multicomponent Balescu-Lenard equation (60) and with the multicomponent Boltzmann distribution (107) that account for a segregation of mass. Heavy particles have the tendency to sink at the center of the system while light particles move around. Various descriptions of this collisionless dynamical friction have already been proposed in astrophysics (Kalnajs, 1971; Tremaine & Weinberg, 1984; Weinberg, 1986, 1989). For a thorough illustration of this process, we refer in particular to Weinberg (1989), which provides a detailed numerical analysis of the sinking of satellites in spherical galaxies, when accounting for or neglecting collective effects.

Note that the diffusion coefficient (135) is proportional to the mass ${\mu_{\mathrm{b}}\!\sim\!1/N_{\mathrm{b}}}$ of the field particles because it is due to the fluctuations of the force produced by these particles (see section 3). When ${N_{\mathrm{b}}\!\rightarrow\!+\infty}$, the fluctuations vanish (the field particles form a collisionless fluid) and the diffusion coefficient tends to zero. On the other hand, the friction force (136) involves a term proportional to ${\mu_{\mathrm{b}}\!\sim\!1/N_{\mathrm{b}}}$ related to the derivative of the diffusion coefficient ${\partial\bm{D}/\partial\bm{J}}$ and a term proportional to ${\mu_{\mathrm{t}}\!\sim\!1}$, which corresponds to the friction by polarisation $\bm{F}_{\rm pol}$ from equation (137) (see the decomposition from equation (58)). When ${N_{\mathrm{b}}\!\rightarrow\!+\infty}$, the first term tends to zero while the second term remains finite. The term in ${\partial\bm{D}/\partial\bm{J}}$ is proportional to the mass $\mu_{\mathrm{b}}$ of the field particles because it is due to the fluctuations of the force produced by these particles. The force by polarisation from equation (137) is proportional to the mass $\mu_{\mathrm{t}}$ of the test particle because it is due to the retroaction of the fluid of stars to the perturbation caused by the test particle (see Appendix B). Therefore, this force remains finite even when the field particles have no fluctuation (${N_{\mathrm{b}}\!\rightarrow\!+\infty}$), i.e. even if the test particle is evolving within a collisionless fluid. This is the collisionless resonant dynamical friction, which accounts for both the system’s inhomogeneity and collective effects. We refer to Chavanis (2013) for a detailed discussion of the links between this formalism and other approaches, in particular the two-body encounters theory pioneered by Chandrasekhar (1943a).

D.3 An ilustration of dynamical friction in a Mestel disc

The calculation of the force of dynamical friction acting on a star is a problem of considerable interest, initiated by the seminal work of Chandrasekhar (1943a). Chandrasekhar & von Neumann (1943) attempted to derive this force from a purely stochastic formalism in the case where the system is infinite and homogeneous and the stars have a Maxwellian velocity distribution. Their calculations were extended by Del Popolo & Gambera (1999); Del Popolo (2003) to an inhomogeneous medium with a density decaying as ${\rho\!\sim\!r^{-p}}$, assuming again that the velocity distribution of the stars is Maxwellian. In order to treat in a self-consistent manner more general situations of spatial inhomogeneity, and take into account collective effects, the formalism developed in the present paper is necessary.111111We emphasise that this formalism recovers in the appropriate limit the two-body encounters theory of Chandrasekhar (1943a) for a ${3D}$ homogeneous medium (this is discussed in detail in Chavanis (2013)). Conversely, the situation presented in e.g.  Del Popolo & Gambera (1999); Del Popolo (2003) cannot be treated with this formalism because the corresponding distribution (albeit physically interesting) is not a quasi-stationary state of the Vlasov equation and, as such, does not possess suitable angle-action coordinates. Explicit applications were carried for razor-thin (Fouvry et al., 2015b; Fouvry et al., 2015c), thickened (Fouvry et al., 2017b) stellar discs, and Keplerian discs (Fouvry et al., 2017a). Let us illustrate this kinetic theory with the calculation of the friction force acting on a sinking satellite in a collisionless fluid of stars.

Computing the resonant collisionless dynamical friction acting on a massive perturber requires the construction of the angle-action coordinates ${(\bm{\theta},\bm{J})}$, the specification of the biorthogonal basis ${(\psi^{(p)},\rho^{(p)})}$, the computation of the system’s response matrix ${\widehat{\mathbf{M}}(\omega)}$, and the resolution of the non-local resonance condition ${\delta_{\mathrm{D}}(\bm{m}_{1}\!\cdot\!\bm{\Omega}_{1}\!-\!\bm{m}_{2}\!\cdot\!\bm{\Omega}_{2})}$. Fouvry et al. (2015c) presented this calculation for razor-thin axisymmetric stellar discs and recovered the spontaneous self-consistent formation on secular timescales of a narrow resonant ridge of orbits in action space, as first observed in the numerical simulations of Sellwood (2012).

Relying on the calculations performed in Fouvry et al. (2015c), let us estimate the friction force by polarisation that a massive perturber undergoes when embedded in a collisionless disc (i.e. in the limit of an infinite number of bath particles in the disc). Let us specifically consider an infinitely thin Mestel disc for which the circular speed is a constant $V_{0}$ independent of the radius. The stationary background potential $\psi_{\rm M}$ and its associated surface density $\Sigma_{\rm M}$ are given by

[TABLE]

where $R_{\rm max}$ is a scale parameter of the disc. Following Binney & Tremaine (2008), a self-consistent DF for this system is given by

[TABLE]

where the exponent $q$ is given by $q={V_{0}^{2}}/{\sigma_{r}^{2}}\!-\!1\,,$ with $\sigma_{r}$ being the constant radial velocities spread within the disc. In equation (161), $C_{\rm M}$ is a normalisation constant. In order to ensure linear stability, the DF from equation (161) is additionally tapered in the inner and outer regions of the disc. We refer to Fouvry et al. (2015c) for further details on the physical system considered. Razor-thin axisymmetric discs are explicitly integrable, so that one naturally introduces the action coordinates ${\bm{J}\!=\!(J_{\phi},J_{r})}$. Here, $J_{\phi}$ is the azimuthal action of the particle, its angular momentum, and encodes the typical distance of the particle to the center. The second action is $J_{r}$, the radial action, which captures the amplitude of the particle’s radial libration. The larger $J_{r}$, the wider the radial oscillations and the hotter the orbit. Because the radial action satisfies ${J_{r}\!\geq\!0}$, exactly circular orbits (i.e. orbits with ${J_{r}\!=\!0}$) remain circular during the secular evolution and can therefore only diffuse along the ${J_{\phi}-}$direction. For simplicity, the calculation will be restricted to massive perturbers on circular orbits in the collisionless limit, assuming that the razor-thin disc is made of an infinite number of particles. See also Weinberg (1989) for an example of calculation of the dressed dynamical friction acting on circular orbits in the case of ${3D}$ spherical systems. Following equation (155), the evolution of a massive perturber on a circular orbit of angular momentum $J_{\phi}^{\mathrm{t}}$ is given by the one-dimensional differential equation

[TABLE]

where $F_{\rm pol}^{\phi}$ stands for the component along the ${J_{\phi}-}$direction of the ${2D}$ friction force vector $\bm{F}_{\rm pol}$. Through equation (58), such a friction force is straightforwardly estimated once the disc’s self-consistent first- and second-order diffusion coefficients are computed. These coefficients were computed in Fouvry et al. (2015c) for a razor-thin tapered Mestel disc. In the collisionless limit (i.e. ${\mu_{\star}\!\to\!0}$, with $\mu_{\star}$ the individual mass of the particles forming the disc), the friction force ${F_{\rm pol}^{\phi}(J_{\phi}^{\mathrm{t}})}$ acting on a massive perturber of mass $\mu_{\mathrm{t}}$ on the circular orbit $J_{\phi}^{\mathrm{t}}$ follows. This is illustrated in figure 1, which performs the calculations respectively within the Balescu-Lenard and Landau frameworks, i.e. with or without collective effects.

Figure 1 exhibits a complex behaviour for the friction force by polarisation. First note the importance of collective effects, which significantly hasten the dynamical friction. Around ${J_{\phi}^{\mathrm{t}}\!\simeq\!2.8}$ the Balescu-Lenard formalism predicts ${F_{\rm pol}^{\phi}\!>\!0}$, i.e. a diffusion towards higher angular momentum. This abnormal diffusion is also found in the lower panel of figure 8 in Sellwood (2012), where the individual diffusion of stars also follows this trend. This can be interpreted as a consequence of the corotation resonance, whose inner Lindblad resonance is at the central ridge (${J_{\phi}^{\mathrm{t}}\!\simeq\!1.2}$). When collective effects are neglected, this region of positive friction vanishes. One should also keep in mind that the geometrical constraints associated with the razor-thin geometry (compared to ${3D}$ spherical systems) impact the expected properties of the friction force, leading in some regimes to its complete cancellation (Kalnajs, 1971). Understanding the properties of the resonant collisionless dynamical friction in razor-thin discs in more details could be the subject of future works.

Appendix E Some distribution formulae

Let us compute here the few distributions encountered in the main text.

E.1 Calculation of ${\eta(x)}$

Let us first see how the distribution

[TABLE]

operates on a given function ${f(x)}$. Writing ${u\!=\!x\Delta t}$, one has

[TABLE]

The integral on $u$ converges and can be computed noting that

[TABLE]

Hence

[TABLE]

E.2 Calculation of ${\kappa(x)}$

We now compute

[TABLE]

Let us introduce ${\tau\!=\!t_{1}\!-\!t_{2}}$ and reverse integration between $t_{1}$ and $\tau$ while keeping track of the triangular shape of the ${(t_{1},\tau)}$ integration domain. One gets

[TABLE]

We recall the identity

[TABLE]

Hence

[TABLE]

E.3 Calculation of ${\gamma(x)}$

Let us finally compute

[TABLE]

Reversing integration between $t_{2}$ and $t_{3}$ yields

[TABLE]

Appendix F Symmetries

F.1 Relation between $\Lambda_{-\bm{m}_{1},-\bm{m}_{2}}(\bm{J}_{1},\bm{J}_{2},-\omega)$ and $\Lambda_{\bm{m}_{1},\bm{m}_{2}}(\bm{J}_{1},\bm{J}_{2},\omega)$

Let us first show that

[TABLE]

Recall that

[TABLE]

Following equations (25) and (A25) from Heyvaerts (2010), one can write

[TABLE]

where ${\widehat{\alpha}}$ is an element of the basis, which is in general different from $\alpha$ (see Heyvaerts (2010)). This immediately gives

[TABLE]

F.2 Relation between $\Lambda_{\bm{m}_{2},\bm{m}_{1}}(\bm{J}_{2},\bm{J}_{1},\omega)$ and $\Lambda_{\bm{m}_{1},\bm{m}_{2}}(\bm{J}_{1},\bm{J}_{2},\omega)$

Let us now demonstrate that to second order in the noise level, one can assume that ${\Lambda_{\bm{m}_{2},\bm{m}_{1}}(\bm{J}_{2},\bm{J}_{1},\omega)\!=\!\Lambda_{\bm{m}_{1},\bm{m}_{2}}^{*}(\bm{J}_{1},\bm{J}_{2},\omega)}$. First note that in equation (92) the two frequencies $\omega_{1}$ and $\omega_{2}$ coincide because of the Dirac Delta factor ${\delta_{\mathrm{D}}(\omega_{1}\!-\!\omega_{2})}$. Now we have

[TABLE]

where, in the last line, we introduced $\varepsilon^{-1\dagger}$ as the hermitian conjugate of $\varepsilon^{-1}$. The anti-hermitic part of $\varepsilon^{-1}$, which appears in the square bracket, was computed by Heyvaerts (2010) and can be written as

[TABLE]

Now we have (see equation (25) in Heyvaerts (2010))

[TABLE]

Injecting this relation into equation (177) yields

[TABLE]

This relation shows that the difference between ${\Lambda_{\bm{m}_{2},\bm{m}_{1}}(\bm{J}_{2},\bm{J}_{1},\omega)}$ and ${\Lambda_{\bm{m}_{1},\bm{m}_{2}}^{*}(\bm{J}_{1},\bm{J}_{2},\omega)}$ is a term involving two $\Lambda$ factors, which implies that it is of higher order with respect to the noise, as discussed in the end of section 2. This difference is proportional to the anti-hermitic part of $\varepsilon$, and therefore corresponds to a particular form of the fluctuation-dissipation theorem. Since equation (92) aims for a second order expression in the noise, the difference {\delta_{\mathrm{D}}(\omega_{1}\!-\!\omega_{2})\big{[}\Lambda_{\bm{m}_{2},\bm{m}_{1}}(\bm{J}_{2},\bm{J}_{1},\omega_{1})\!-\!\Lambda_{\bm{m}_{1},\bm{m}_{2}}^{*}(\bm{J}_{1},\bm{J}_{2},\omega_{2})\big{]}}, because it would introduce a third order correction, can be neglected here.