Kinetic theory of ${1D}$ homogeneous long-range interacting systems sourced by ${1/N^{2}}$ effects

TL;DR

This paper derives a new kinetic equation for 1D homogeneous long-range systems that relax very slowly under 1/N^2 effects, explaining their long-term evolution and relaxation to equilibrium.

Contribution

It introduces a closed, explicit kinetic equation for 1D homogeneous long-range systems that accounts for 1/N^2 effects and matches N-body simulations.

Findings

Kinetic blocking prevents relaxation under 1/N effects in 1D homogeneous systems.

The new kinetic equation satisfies an H-Theorem ensuring relaxation to equilibrium.

Quantitative agreement between the kinetic equation and N-body simulations.

Abstract

The long-term dynamics of long-range interacting $N$-body systems can generically be described by the Balescu-Lenard kinetic equation. However, for ${1D}$ homogeneous systems, this collision operator exactly vanishes by symmetry. These systems undergo a kinetic blocking, and cannot relax as a whole under ${1/N}$ resonant effects. As a result, these systems can only relax under ${1/N^{2}}$ effects, and their relaxation is drastically slowed down. In the context of the homogeneous Hamiltonian Mean Field model, we present a new, closed and explicit kinetic equation describing self-consistently the very long-term evolution of such systems, in the limit where collective effects can be neglected, i.e. for dynamically hot initial conditions. We show in particular how that kinetic equation satisfies an $H$-Theorem that guarantees the unavoidable relaxation to the Boltzmann equilibrium distribution. Finally, we illustrate how that kinetic equation quantitatively matches with the measurements from direct $N$-body simulations.

Equations

$$H=\frac{1}{2}\sum _{i=1}^N{v}_i^2+\mu \sum _{i<j}^{N}U({\theta }_i-{\theta }_j),$$

$$U({\theta }_i-{\theta }_j)=-U_0\cos ({\theta }_i-{\theta }_j),$$

$$\frac{\partial F(v_1)}{\partial t}=$$

$$\times$$

$$\frac{\partial F(v_1)}{\partial t}$$

$$\displaystyle\times\!\!\int\!\!\mathrm{d}v_{3}\bigg{\{}$$

$$+$$

$${\mathbf{k}}_1=(2,-1,-1);{\mathbf{k}}_2=(1,-2,1).$$

$$\frac{\partial F(v_{1})}{\partial t}=\frac{\partial}{\partial v_{1}}\bigg{[}\mathcal{F}(v_{1})\bigg{]},$$

$$\frac{\partial F(v_{1})}{\partial t}=\frac{\pi^{3}}{2}U_{0}^{4}\mu^{2}\frac{\partial}{\partial v_{1}}\bigg{[}\mathcal{P}\!\!\int\!\!\frac{\mathrm{d}v_{2}}{(v_{1}-v_{2})^{4}}\,K(v_{1},v_{2})\bigg{]},$$

$$K(v_{1},v_{1}+\delta v)=K_{3}(v_{1})\,(\delta v)^{3}+\mathcal{O}\big{(}(\delta v)^{4}\big{)},$$

$$K_3=F^{(4)}FF-F^{(3)}{F}^{\prime }F+3F^{\prime \prime }{F}^{\prime }{F}^{\prime }-3F{F}^{\prime \prime }{F}^{\prime \prime }.$$

$$\mathcal{P}\int \frac{\mathrm{d}{v}_2}{(v_1-v_2{)}^4}K(v_1,v_2)$$

$$\displaystyle\,\sim\mathcal{P}\!\!\int\!\!\mathrm{d}\delta v\,\bigg{\{}\frac{K_{3}}{\delta v}+\mathcal{O}(1)\bigg{\}},$$

$$F_{\mathrm{B}}(v)=A{\mathrm{e}}^{-\beta (v-v_0{)}^2},$$

$$\frac{\partial F_{\mathrm{B}}(v_1)}{\partial t}$$

$$=0.$$

$$M(t)$$

$$P(t)$$

$$E(t)$$

$$\frac{\mathrm{d}M}{\mathrm{d}t}=\!\!\int\!\!\mathrm{d}v_{1}\,\frac{\partial}{\partial v_{1}}\bigg{[}\mathcal{F}(v_{1})\bigg{]}=0.$$

$$\frac{\mathrm{d}P}{\mathrm{d}t}$$

$$\frac{\mathrm{d}E}{\mathrm{d}t}=-\int \mathrm{d}{v}_1{v}_1\mathcal{F}(v_1).$$

$$\frac{\mathrm{d}E}{\mathrm{d}t}$$

$$\displaystyle\times\bigg{\{}$$

$$+$$

$$=0,$$

$$S(t)=-\int \mathrm{d}{v}_{1}s(F(v_1,t)),$$

$$\frac{\mathrm{d}S}{\mathrm{d}t}=\int \mathrm{d}{v}_1\frac{1}{F(v_1)}\frac{\partial F(v_1)}{\partial v_1}\mathcal{F}(v_1).$$

$$\frac{\mathrm{d}S}{\mathrm{d}t}=\frac{{\pi }^3}{8}{U}_0^4{\mu }^2$$

$$\displaystyle\times\bigg{\{}$$

$$+$$

$$\frac{\mathrm{d}S}{\mathrm{d}t}\ge 0.$$

$$\forall v,v^{\prime}\,:\;G\bigg{(}\frac{v+v^{\prime}}{2}\bigg{)}=\frac{G(v)+G(v^{\prime})}{2}.$$

$$G(v)=-2\beta (v-v_0)\text{ with }\beta >0,$$

$${\sigma }^2=\frac{1}{M_{\mathrm{tot}}}\int \mathrm{d}\theta \mathrm{d}v{v}^{2}F(v).$$

Full text

I Introduction

The evolution of long-range interacting systems (Campa et al., 2014) is generically composed of two stages. First, the system undergoes a rapid (collisionless) violent relaxation (Lynden-Bell, 1967), which, owing to strong potential fluctuations, allows for the system to reach a quasistationary state (i.e. a steady state of the Vlasov equation). Following these drastic orbital rearrangements, the system has become dynamically frozen for the mean-field dynamics. It is then only through finite-$N$ effects, i.e. Poisson shot noise due to the finite number of particles, that the system can keep evolving. The system undergoes then a slow (collisional) relaxation that drives it to thermodynamical equilibrium.

The efficiency of the long-term dynamical relaxation of a system therefore depends on $N$, the system’s number of particles. This dependence is generically accounted for by adding collisional operators to the collisionless Vlasov equation, i.e. by computing the appropriate kinetic equation. As an example, in the limit where ${1/N}$ effects dominate the dynamics, such an evolution can generically be described by the Balescu-Lenard equation Balescu (1960); Lenard (1960), whose generalisation to inhomogeneous systems is only recent (Heyvaerts, 2010; Chavanis, 2012a). In that context, two body (resonant) encounters lead to a collision operator of order ${1/N}$, so that the relaxation occurs on a timescale scaling like ${Nt_{\mathrm{d}}}$, with $t_{\mathrm{d}}$ the system’s dynamical time222One important exception is ${3D}$ self-gravitating systems, where the relaxation time is of order ${Nt_{\mathrm{d}}/\log(N)}$, owing to the Coulomb logarithm Binney and Tremaine (2008)..

Yet, for ${1D}$ homogeneous systems, such kinetic equations exactly vanish, i.e. two-body encounters are unable to lead to an overall relaxation of the system, as highlighted in the context of ${1D}$ plasmas (Eldridge and Feix, 1963; Dawson, 1964; Rouet and Feix, 1991), the ${1D}$ Hamiltonian Mean Field (HMF) model (Zanette and Montemurro, 2003; Yamaguchi et al., 2004; Bouchet and Dauxois, 2005; Campa et al., 2007, 2008; Rocha Filho et al., 2014), the dynamics of long-range coupled particles on the unit sphere Gupta and Mukamel (2011); Barré and Gupta (2014); Lourenço and Rocha Filho (2015); Fouvry et al. (2019a), or even the dynamics of point vortices (Chavanis and Lemou, 2007). All these systems are said to being undergoing a kinetic blocking. For such systems, it is only three-body effects and higher-order correlations that can drive the system’s relaxation to thermodynamical equilibrium (Bouchet et al., 2010; Chavanis, 2010, 2012b; Sano, 2012), making the relaxation time much longer than ${Nt_{\mathrm{d}}}$.

As it must originate from perturbations to the system’s dynamics of increasing order in ${1/N}$, it is natural to expect that the timescale for the collisional relaxation of a ${1D}$ homogeneous system would scale like ${N^{2}t_{\mathrm{d}}}$. Indeed, this ${N^{2}t_{\mathrm{d}}}$ scaling of the dynamics was observed for ${1D}$ plasmas (Eldridge and Feix, 1963; Dawson, 1964; Rouet and Feix, 1991), or for long-range coupled particles on the sphere (Lourenço and Rocha Filho, 2015; Fouvry et al., 2019a). In the case of the HMF model, differents scalings proportional to ${N^{1.7}t_{\mathrm{d}}}$ (Zanette and Montemurro, 2003; Yamaguchi et al., 2004), or even ${\mathrm{e}^{N}t_{\mathrm{d}}}$ (Campa et al., 2007) were reported. But, in Rocha Filho et al. (2014), these results were convincingly interpreted as being side effects associated with a too small value of $N$, and a scaling in $N^{2}$ was recovered through a careful analysis of simulations with larger values of $N$.

In the present paper, building upon Rocha Filho et al. (2014), we set out to study such a very long-term dynamics of the HMF model in the homogeneous limit. In the limit where collective effects are neglected (i.e. the neglect of the ability of the mean system to amplify perturbations), we present a new, closed kinetic equation describing the collisional relaxation of that system on ${N^{2}t_{\mathrm{d}}}$ timescales, as driven by three-body correlations. We explore the generic properties of this new collision operator, and quantitatively compare its predictions to direct $N$-body simulations of that system.

The paper is organised as follows. In Section II, we briefly present the considered HMF model, and the kinetic equation describing its relaxation at order ${1/N^{2}}$, as given by Eq. (4). The detailed procedure followed to obtain that equation is presented in Appendices A–D, while the effective analytical calculations were performed using a computer algebra system (Fouvry et al., 2019b). In Section III, we explore some of the fundamental properties of that kinetic equation, in particular its well-posedness, its conservation properties, and its $H$–Theorem that guarantees the relaxation to the homogeneous Boltzmann equilibrium, provided that it is linearly stable. Finally, in Section IV, we quantitatively illustrate how this kinetic equation matches with direct measurements from numerical simulations, for hot enough initial distributions. Finally, we conclude in Section V.

II The kinetic equation

We are interested in the long-term dynamics of the HMF model (Antoni and Ruffo, 1995). It is composed of $N$ particles of individual mass ${\mu=M_{\mathrm{tot}}/N}$, with $M_{\mathrm{tot}}$ the system’s total mass. The canonical phase space coordinates are ${(\theta,v)}$, and the total Hamiltonian reads

[TABLE]

where the pairwise interaction potential is given by

[TABLE]

with ${U_{0}>0}$ the amplitude of the pairwise coupling.

In the homogeneous limit, the instantaneous statistical state of the system can be described by the velocity distribution function (DF), ${F=F(v,t)}$, which, following Eq. (36), is taken to be normalised as ${\!\int\!\mathrm{d}\theta\mathrm{d}vF(v,t)=M_{\mathrm{tot}}}$. Describing the long-term relaxation of such a system, amounts then to describing the long-term evolution of that DF, as driven by a closed kinetic equation.

When limiting oneself only to ${1/N}$ effects, the dynamics of that DF is generically given by the Balescu-Lenard equation, which in the present context reads (Benetti and Marcos, 2017)

[TABLE]

where, to shorten the notations, we do not write explicitly the time-dependence of the DFs. We also introduced the dielectric coefficient, ${\varepsilon_{k}(\omega)}$, explicitly spelled out in Eq. (65). As already mentioned in Introduction, owing to the Dirac delta, ${\delta_{\mathrm{D}}(v_{1}-v_{2})}$, such a collision operator exactly vanishes by symmetry. The homogeneous HMF system undergoes a kinetic blocking, and its overall relaxation is immune to ${1/N}$ two-body effects.

As a result, it is only by being driven by weaker three-body correlations, associated with ${1/N^{2}}$ effects, that the present system can relax to its thermodynamical equilibrium. This is the dynamics of interest in this paper. Following an approach similar to Rocha Filho et al. (2014), we present in Appendices A–D our approach to derive such a closed kinetic equation accounting for ${1/N^{2}}$ effects.

The main steps of this derivation are as follows. (i) As detailed in Appendix A, we first derive the usual BBGKY equations, by obtaining the coupled evolution equations for the system’s $1$-, $2$-, and $3$-body distribution functions. (ii) As shown in Appendix B, using the cluster expansion, these coupled evolution equations are written as coupled evolution equations for the system’s $1$-body DF, ${F(v)}$, and the $2$-, and $3$-body correlation functions. The main gain of this rewriting is that these equations are now sorted by increasing order in ${1/N}$ corrections. (iii) As presented in Appendix C, these equations are then truncated at order ${1/N^{2}}$. At this stage, a key simplification comes from our neglect of the contributions from collective effects, i.e. the system’s ability to amplify perturbations, as is justified for dynamically hot initial distributions. (iv) Finally, in Appendix D, we show how this set of (well-posed) coupled partial differential equations can be solved, allowing for an explicit and closed expression for the collision operator. While not intrinsically challenging, such calculations are made cumbersome because of the large number of terms involved. These calculations were therefore carried out using Mathematica and are spelled out in detail in Fouvry et al. (2019b).

All in all, the final kinetic equation derived in that fashion reads

[TABLE]

where, to shorten the notations, we introduced the velocity vector ${\mathbf{v}=(v_{1},v_{2},v_{3})}$, as well as the two resonance vectors

[TABLE]

In that equation, we also introduced $\mathcal{P}$ as the Cauchy principal value, which acts on the integral ${\!\int\!\mathrm{d}v_{2}}$. We postpone to Section III.1 the justification of the well-posedness of such a principal value. Finally, we note that Eq. (4) is tightly related to the ${1/N^{2}}$ kinetic equation already put forward in Eq. (23) of Rocha Filho et al. (2014). The differences are some corrections in the overall prefactor, and the sign of the second resonant term. Compared to Rocha Filho et al. (2014), in Eq. (4), we also performed additional rewritings and manipulations, that offer a simpler collision operator, involving only one principal value, and only up to first-order gradients in the system’s DF, as detailed at the end of Appendix D.

As usual, it is possible to rewrite Eq. (4) under the form of a continuity equation, reading

[TABLE]

where the instantaneous flux in velocity space333With such a convention, the flux is opposite to the direction effectively followed by individual particles during their diffusion., ${\mathcal{F}(v_{1})}$, follows directly from Eq.(4).

As expected, Eq. (4) is proportional to ${\mu^{2}\simeq 1/N^{2}}$, i.e. it describes a collisional relaxation on ${N^{2}t_{\mathrm{d}}}$ timescales. We also note that the collision operator in the r.h.s. involves the system’s DF three times, highlighting the fact that this kinetic equation describes a dynamics sourced by three-body correlations. These correlations are matched through two different resonance conditions on the velocities, namely ${\delta_{\mathrm{D}}(\mathbf{k}_{1/2}\!\cdot\!\mathbf{v})}$. We also note that the two resonance terms are opposite one to another, provided one makes the change ${(v_{1}\leftrightarrow v_{2})}$ in the last integrand. This will prove important to ensure some of the equation’s conservation properties, as detailed in Section III.3. Equation (4) is the main result of this section, as this closed kinetic equation is the appropriate kinetic equation to describe the long-term evolution of a dynamically hot ${1D}$ homogeneous systems, sourced by ${1/N^{2}}$ effects, and tailored here to the particular case of the HMF model. We finally note that Eq. (4) only holds as long as the system’s mean DF remains linearly stable, see the end of Section IV for a more detailed discussion.

III Properties

In this section, we now explore some of the key properties of the kinetic Eq. (4).

III.1 Well-posedness

Given the presence of a high-order resonance denominator in Eq. (4), it is not obvious that the kinetic equation is well-defined, i.e. that there are no divergences when ${v_{2}\to v_{1}}$. As a result, let us study the behaviour of the integrand in the limit ${v_{2}\to v_{1}}$. In order to shorten the notations, we temporarily rewrite Eq. (4) as

[TABLE]

where the function ${K(v_{1},v_{2})}$ immediately reads from Eq. (4). Assuming that ${F(v)}$ is a smooth function, it is straightforward to perform a limited development of ${K(v_{1},v_{1}+\delta v)}$ for ${\delta v\to 0}$. One gets

[TABLE]

where the first non-zero coefficient, ${K_{3}(v)}$, reads

[TABLE]

As a consequence, in the vicinity of ${v_{2}\to v_{1}}$, the integral from Eq. (7) takes the form

[TABLE]

which is a meaningful and well-posed integral in terms of a principal value.

III.2 Boltzmann distribution

The thermodynamical equilibrium states resulting from the collisional relaxation of a homogeneous $N$-body system are expected to be (shifted) homogeneous Boltzmann distributions of the form

[TABLE]

where $\beta$ is the inverse temperature, and $A$ a normalisation constant.

Owing to the explicit form of the collision operator from Eq. (4), it is straightforward to check that such DFs are indeed equilibrium solutions of the kinetic equation. Indeed, noting that the resonance vectors, $\mathbf{k}_{1}$ and $\mathbf{k}_{2}$, from Eq. (5) are of zero sum, one has

[TABLE]

This highlights that the diffusion flux for homogeneous Boltzmann distributions exactly vanishes, i.e. these DFs are equilibrium solutions of the ${1/N^{2}}$ kinetic Eq. (4). In Section III.4, owing to an $H$–Theorem, we will strengthen this result by proving that the homogeneous Boltzmann DFs from Eq. (11) are the only equilibrium solutions of the present kinetic equation.

III.3 Conservation laws

The kinetic equation (4) satisfies various conservation properties, in particular the conservation of the total mass ${M(t)}$, the total momentum, ${P(t)}$, and the total energy, ${E(t)}$, as we will now briefly justify. Ignoring numerical prefactors, they are respectively defined as

[TABLE]

where the total energy only contains the kinetic energy, because we assumed that the system remains homogeneous on average.

Following the rewriting from Eq. (6), the conservation of the total mass follows from the absence of any boundary contributions, so that

[TABLE]

A similar calculation can be pursued for the total momentum, and one gets

[TABLE]

using the symmetrisation ${(v_{1}\leftrightarrow v_{2})}$ in Eq. (4).

Finally, regarding the conservation of energy, following an integration by parts of Eq. (13), one writes

[TABLE]

Using the definition of the flux from Eq. (6), this expression will then involve an integral of the form ${\!\int\!\mathrm{d}v_{1}\mathrm{d}v_{2}\mathrm{d}v_{3}}$, which allows us to use symmetrisations w.r.t. the integration variables. First, we symmetrise all the terms w.r.t. the permutation ${(v_{1}\leftrightarrow v_{2})}$. Then, for the subsequent expression, we perform two additional symmetrisations, namely (i) ${(v_{2}\leftrightarrow v_{3})}$ for the terms involving the resonance condition ${\delta_{\mathrm{D}}(\mathbf{k}_{1}\cdot\mathbf{v})}$, and (ii) ${(v_{1}\leftrightarrow v_{3})}$ for the terms involving the resonance condition ${\delta_{\mathrm{D}}(\mathbf{k}_{2}\cdot\mathbf{v})}$. By doing so, the resonant denominator from Eq. (4) remains a sole function ${(v_{1}\!-\!v_{2})}$, and this avoids the creation of any other type of resonance conditions. All these calculations are straightforward and carried out in detail in Fouvry et al. (2019b). Forgetting prefactors, one gets

[TABLE]

which exactly vanishes owing to the resonance conditions.

III.4 $H$–Theorem

We define the system’s instantaneous entropy as

[TABLE]

with ${s(F)=F\ln(F)}$ the entropy functional. Starting from the rewriting of Eq. (6), it is straightforward to show that the system’s entropy evolves according to

[TABLE]

Following Eq. (6), this expression involves an integral of the form ${\!\int\!\mathrm{d}v_{1}\mathrm{d}v_{2}\mathrm{d}v_{3}}$, allowing for symmetrisations w.r.t. the integration variables. We perform the exact same symmetrisations as the one performed in Eq. (17) to check for energy conservation. All these calculations are straightforward, and carried out in detail in Fouvry et al. (2019b). One gets

[TABLE]

where we used the shortening notation ${F_{1}=F(v_{1})}$, and ${F_{1}^{\prime}=\partial F/\partial v_{1}}$. Given that all the terms involved in this integral are positive, the kinetic Eq. (4) satisfies an $H$–Theorem, i.e. one has

[TABLE]

The expression of the entropy increase from Eq. (20) allows us then to tackle the question of determining which DFs are equilibrium states for the diffusion, i.e. which DFs satisfy ${\mathrm{d}S/\mathrm{d}t=0}$. Provided that one uses the symmetrisation ${(v_{1}\leftrightarrow v_{2})}$, the constraints associated with the two resonance conditions from Eq. (20) are identical, so that we only need to consider one. Recalling the expression of the resonance vector $\mathbf{k}_{1}$ from Eq. (5), and introducing the function ${G(v)=F^{\prime}(v)/F(v)}$, we note that a DF is stationary if it satisfies

[TABLE]

Because this constraint has to be satisfied for all $v$ and $v^{\prime}$, we can conclude that ${v\mapsto G(v)}$ has to be a line, i.e. one has

[TABLE]

where the constraint ${\beta>0}$ stems from the fact that ${\!\int\!\mathrm{d}vF(v)=M_{\mathrm{tot}}<+\infty}$, i.e. the DF has to be normalised, and cannot get infinitely large for ${v\to+\infty}$. Equation (23) immediately translates to the differential equation ${F^{\prime}(v)/F(v)=-2\beta(v-v_{0})}$, which naturally integrates to the (shifted) homogeneous Boltzmann DF introduced in Eq. (11).

As a conclusion, the only equilibrium DFs of the kinetic Eq. (4) are the (shifted) homogeneous Boltzmann distributions. This is an important result of this section. Indeed, while any (stable) DF ${F(v)}$ was an equilibrium distribution for the ${1/N}$-dynamics of an homogeneous long-range interacting system, only homogeneous Boltzmann DFs are equilibrium distributions for these systems’ ${1/N^{2}}$-dynamics.

III.5 Dimensionless rewriting

In order to have a better grasp at the scalings of Eq. (4), let us finally rewrite it under a dimensionless form.

Following the conservation of total energy obtained in Eq. (17), we introduce the system’s (conserved) velocity dispersion as

[TABLE]

This typical velocity entices us then to define a dimensionless velocity as ${u=v/\sigma}$, and a dimensionless time as ${\overline{t}=t/t_{\mathrm{d}}}$ (time) with ${t_{\mathrm{d}}=1/\sigma}$ the dynamical time. Similarly, we define the system’s dimensionless probability distribution function (PDF) as

[TABLE]

that satisfies the normalisation condition ${\!\int\!\mathrm{d}u\,\overline{F}(u)=1}$. Finally, in order to assess the “dynamical temperature” of the system and the strength of the associated collective effects, we introduce the dimensionless stability parameter

[TABLE]

following a notation similar to (Toomre, 1964). The larger $Q$, the more stable the system, and the weaker the collective effects. In Appendix E, we motivate the definition of $Q$, and directly relate it to the system’s dielectric function.

Using these conventions, one can rewrite Eq. (4) as

[TABLE]

where, similarly to Eq. (4), we introduced the velocity vector ${\mathbf{u}=(u_{1},u_{2},u_{3})}$, and the resonance vectors, ${(\mathbf{k}_{1},\mathbf{k}_{2})}$, are given by Eq. (5). Similarly to Eq. (6), we can rewrite Eq. (27) under the form of a continuity equation, reading

[TABLE]

where the dimensionless instantaneous flux, ${\overline{\mathcal{F}}(u_{1})}$, follows from Eq. (27).

Equation (27) is the appropriate dimensionless writing to understand the expected relaxation time of a given system. Indeed, assuming that the collision operator within brackets is of order unity, we find therefore that the relaxation time scales like

[TABLE]

It is interesting to note that one recovers that dynamically colder systems, i.e. systems with smaller values of $Q$, relax faster than hotter systems. However, because collective effects were neglected in the derivation of Eq. (27), one has to place oneself in the regime ${Q\gg 1}$ for the present kinetic equation to apply. In that dynamically hot regime, collective effects are indeed unimportant, but, because of the factor $Q^{4}$ in Eq. (29), relaxation will only occur on very long timescales.

IV Applications

In order to investigate the validity of the kinetic Eq. (4), we now set out to explore numerically the long-term relaxation of such systems, and compare it with the kinetic prediction.

For clarity, all the details of our numerical implementation are given in Appendix F. The main difficulty with such a numerical exploration comes from our neglect of collective effects in the derivation of the kinetic Eq. (4). As defined in Eq. (26), this asks therefore for the consideration of initial conditions with ${Q\gg 1}$, for which, following Eq. (29), the relaxation can only occur on very late timescales, making the $N$-body simulations more challenging. The larger $Q$, the weaker the collective effects (e.g., as can be seen in Fig. 4), and therefore the better should be the match between the kinetic prediction and the $N$-body measurements.

As a first illustration, we present in Fig. 1, an example of a system’s relaxation towards equilibrium, for an initial condition following the non-Gaussian PDF from Eq. (70).

As expected, even if any homogeneous DF, ${F=F(v)}$, is submitted to a kinetic blocking, and undergoes no relaxation through the ${1/N}$ Balescu-Lenard Eq. (3), it can still relax as a result of higher-order correlation effects, e.g., as captured by the kinetic Eq. (4), whose detailed predictions we may now compute.

In Fig. 2, we illustrate the initial dimensionless flux, ${\overline{\mathcal{F}}(u,t\!=\!0)}$, as defined in Eq. (28), using on the one hand direct measurements from $N$-body simulations (following the method presented in Appendix F), and on the other hand computing the prediction from the kinetic Eq. (27).

As highlighted, the larger $Q$, the hotter the system, i.e. the weaker the collective effects, and therefore the better the matching between the $N$-body measurements and the kinetic prediction. For systems with smaller velocity dispersions, Eq. (27) does not apply anymore, and asks to be generalised in order to account for the contribution of collective effects to hasten or slow down the system’s relaxation. Finally, there are (at least) two possible origins for the slight mismatch still observed in Fig. 2 between the measured fluxes and the predicted one: (i) remaining contributions associated with collective effects, that are expected to slowly fade out as one increases $Q$; (ii) some non-vanishing kinetic contributions from the term in ${G_{2}^{(1)}\!\times\!G_{2}^{(1)}}$ that was neglected in Appendix C, when deriving the system’s truncated BBGKY evolution equations. Such generalisations are beyond the scope of this paper.

As can be noted from the overall prefactor in Eq. (27), one expects the timescale for the system’s relaxation to scale like ${N^{2}t_{\mathrm{d}}}$, w.r.t. $N$ the number of particles. This is investigated in Fig. 3, where we illustrate the dependence of the system’s relaxation efficiency, ${\overline{\mathcal{E}}=\!\int\!\mathrm{d}u\,|\overline{\mathcal{F}}(u)|}$, as a function of the number of particles.

In that figure, we observe that the dimensionless relaxation efficiency is independent of $N$, so that the scaling of the relaxation time from Eq. (29) is indeed recovered.

Before concluding, let us finally briefly describe the system’s possible dynamics depending on the value of the dynamical temperature $Q$. As derived in Eq. (67), we have shown that for ${Q\leq Q_{\mathrm{c}}}$, the homogeneous system is initially linearly unstable. This then leads to the following possible behaviours.

(i) For ${Q\leq Q_{\mathrm{c}}}$, the homogeneous system is initially linearly unstable, so that it rapidly becomes inhomogeneous. For such an inhomogeneous configuration, the kinetic blocking from Eq. (3) does not hold anymore. Provided that the system remains inhomogeneous, its dynamics is driven by the inhomogeneous Balescu-Lenard equation, and drives a relaxation, whose associated relaxation time scales linearly with $N$ (Benetti and Marcos, 2017). As noted in Campa et al. (2008), an interesting case is given by an initial configuration satisfying ${1<Q<Q_{\mathrm{c}}}$. While such an initial condition is initially unstable, its final homogeneous Boltzmann equilibrium is stable (since ${1<Q}$). As a consequence, at some point, the system has to evolve from a inhomogeneous distribution to a homogeneous one, which accordingly delays the relaxation.

(ii) For ${Q\gtrsim Q_{\mathrm{c}}}$, because it is linearly stable, the system remains initially homogeneous, and evolves according to a homogeneous kinetic equation in ${1/N^{2}}$. However, because it is close to the stability threshold, collective effects, i.e. the system’s ability to amplify perturbations, have to be taken into account. This asks for a kinetic equation more general than Eq. (27), where we neglected collective effects. As noted in Campa et al. (2008), an interesting case is given by an initial configuration such that ${Q_{\mathrm{c}}<Q<1}$. Such a system is initially stable so that it will first remain homogeneous and undergo a slow ${1/N^{2}}$ relaxation. Yet, since ${Q<1}$, the associated homogeneous Boltzmann equilibrium is unstable. As a consequence during its (homogeneous) relaxation, the system will unavoidably become unstable at some point. This will drive a dynamical phase transition rapidly making it inhomogeneous. The final stages of the relaxation are then the ones of an inhomogeneous relaxation, that scales in ${1/N}$. Owing to this dynamical phase transition, one expects therefore the system’s overall relaxation time to have an intermediate scaling between $N$ and $N^{2}$.

(iii) Finally, for ${Q\gg Q_{\mathrm{c}}}$, one recovers the case considered in the present paper. Because it is dynamically so hot, the system will not undergo any instability. It will therefore remain homogeneous throughout its relaxation. On the same grounds, collective effects are unimportant and can be neglected. The system’s dynamics is therefore described by Eq. (27), and leads to a relaxation time scaling like $N^{2}$, as illustrated in Fig. 3.

V Conclusion

In the present paper, we focused our attention on the description of the very long-term dynamics of the HMF model in the homogeneous limit, one particular example of a ${1D}$ long-range interacting system. As highlighted in Eq. (3), such systems are generically submitted to a kinetic blocking that prevents their relaxation as a whole under ${1/N}$ resonant effects. As such, their evolution is drastically slowed down, and is only made possible by the cumulative contributions of higher-order ${1/N^{2}}$ effects.

Placing ourselves within the dynamically hot limit, for which collective effects can be neglected, and following an approach similar to Rocha Filho et al. (2014), we showed how one could explicitly solve the BBGKY hierarchy of equations, truncated at order ${1/N^{2}}$. This led us to Eq. (4), a closed, explicit, and self-consistent kinetic equation describing the long-term relaxation of the system’s homogeneous DF as driven by ${1/N^{2}}$ effects, as long as the mean system remains linearly stable.

We put forward the main properties of that new kinetic equation, in particular the fact that it satisfies an $H$–Theorem, that guarantees the unavoidable relaxation of the system towards the homogeneous Boltzmann thermodynamical equilibrium (provided that it is linearly stable). This result highlights the fundamental importance of Boltzmann’s $H$–Theorem that keeps being satisfied in the present ${1/N^{2}}$ context where the relaxation is sourced by the product of three DFs. This therefore further extends the validity of Boltzmann’s $H$–Theorem beyond the traditional kinetic equations such as the collisional Boltzmann, Landau and Balescu-Lenard equations, which are sourced only by product of two DFs.

In Eq. (29), we subsequently detailed how the present formalism predicts a relaxation time scaling like ${t_{\mathrm{r}}\simeq Q^{4}N^{2}t_{\mathrm{d}}}$, with $t_{\mathrm{d}}$ the dynamical time, and $Q$ the system’s stability parameter, as defined in Eq. (26). In particular, this implied a relaxation time scaling like $N^{2}$ w.r.t. the number of particles, a scaling already thoroughly checked in Rocha Filho et al. (2014).

Finally, in Section IV, we presented explicit comparisons of this new kinetic equation with numerical measurements from direct $N$-body simulations. We illustrated in Fig. 1 how at this ${1/N^{2}}$ order, the system does not suffer anymore from a kinetic blocking, and can indeed relax to the homogeneous Boltzmann equilibrium, provided that it is linearly stable. We quantitatively showed in Fig. 2 how the numerically measured diffusion fluxes converge to the kinetic prediction, as the system is made hotter so that collective effects become more and more negligible. We also illustrated in Fig. 3 how the $N^{2}$ scaling of the relaxation time is also recovered numerically.

The kinetic equation presented in Eq. (4) is only a first step towards the detailed characterisation of the (very) long-term dynamics of long-range interacting systems. In the present context, calculations were made more tractable through the following assumptions: (i) the HMF model contains only one harmonic, ${k=\pm 1}$, in its pairwise interaction, reducing drastically the allowed resonances; (ii) we neglected contributions associated with collective effects, which prevented us from having to solve and include the linear response theory of the system; (iii) finally, in the evolution equation for the three-body correlation function, ${\partial G_{3}/\partial t}$, we neglected the contributions from the source term in ${G_{2}^{(1)}\!\times\!G_{2}^{(1)}}$ in Eq. (45), that would have led to an additional collision term proportional to $F^{4}$, instead to $F^{3}$ for the dominant term included in Eq. (4). Further work should try to alleviate these shortcomings, by allowing for more complex resonances, accounting for collective effects to describe dynamically colder systems that are linearly more responsive, and finally by accounting for possible contributions from higher order terms in the system’s DF.

More generally, one should investigate the structure of the collision operators for even higher-order kinetic equations, e.g., at order ${1/N^{3}}$. For example, as recovered in the classical ${1/N}$ Landau and Balescu-Lenard equations, and as recovered here for the ${1/N^{2}}$ kinetic equation, Boltzmann distributions are always found, a posteriori, to be equilibrium states of the collision operator. It would be of interest to investigate whether or not such a property generically holds for higher order expansions, and, if so, understand why.

Acknowledgements.

JBF acknowledges support from Program number HST-HF2-51374 which 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. BB is supported by membership from Martin A. and Helen Chooljian at the Institute for Advanced Study.

Appendix A The BBGKY hierarchy

In this Appendix, we briefly repeat the derivation of the BBGKY hierarchy, to describe the dynamics of a long-range coupled $N$-body system. Notations and normalisations are inspired from the ones considered in Balescu (1997).

We assume that the system is composed of $N$ identical particles of individual mass ${\mu=M_{\mathrm{tot}}/N}$, with $M_{\mathrm{tot}}$ the system’s total mass. We introduce the system’s $N$-body PDF, ${P_{N}(\mathbf{w}_{1},...,\mathbf{w}_{N},t)}$, with ${\mathbf{w}=(\theta,v)}$ the phase space coordinates, normalised so that ${\!\int\!\mathrm{d}\mathbf{w}_{1}...\mathrm{d}\mathbf{w}_{N}P_{N}=1}$. The dynamics of $P_{N}$ is governed by Liouville’s equation

[TABLE]

where we introduced the full $N$-body Hamiltonian

[TABLE]

with ${U(\theta_{i}-\theta_{j})}$ the considered pairwise interaction. In Eq. (30), we also introduced the Poisson bracket over $N$ particles as

[TABLE]

We can subsequently define the system’s reduced PDFs as

[TABLE]

Integrating Eq. (30) w.r.t. all particles but the $n$ first, we obtain the BBGKY hierarchy of equations, namely

[TABLE]

where we used the symmetry of $P_{N}$ w.r.t. exchanges of particles. Similarly to Eq. (31), the $n$-body Hamiltonian $H_{n}$ (resp. ${[\,\cdot\,,\,\cdot\,]_{n}}$ the Poisson bracket over $n$ particles) naturally follows from Eq. (31) (resp. Eq. (32)), provided that one replaces $N$ by $n$. In Eq. (34), we also introduced ${\delta H_{n+1}}$ as the specific interaction energy of the ${(n+1)^{\mathrm{th}}}$ particle with the $n$ first particles. It reads

[TABLE]

As usual, we note that the BBGKY hierarchy from Eq. (34) is not closed, as the evolution equation for ${\partial P_{n}/\partial t}$ involves the higher-order PDF, $P_{n+1}$.

In order to simplify the combinatorial prefactors appearing in Eq. (34), we finally introduce the reduced DFs, $F_{n}$, as

[TABLE]

With such a choice, these DFs scale as ${F_{n}\sim 1}$, w.r.t. $N$ the total number of particles. We can then rewrite the BBGKY hierarchy Eq. (34) under the simple form

[TABLE]

The three first equations of the BBGKY hierarchy, i.e. the evolution equations for $F_{1}$, $F_{2}$, and $F_{3}$, will be the starting point of the derivation of the kinetic equation presented in Eq. (4).

Appendix B The cluster expansion

In Appendix A, we briefly rederived the BBGKY hierarchy of evolution equations for the system’s reduced DFs. In order to be able to perform perturbative developments w.r.t. $N$ the total number of particles, we now introduce the cluster representation of the DFs, following an approach similar to the one presented in Balescu (1997).

We introduce the system’s $2$-body correlation function, ${G_{2}(\mathbf{w}_{1},\mathbf{w}_{2})}$, as

[TABLE]

where we used the shortened notation ${F_{1}(1)=F_{1}(\mathbf{w}_{1})}$. This correlation function characterises how much the statistics of the distribution of $2$ particles differs from being separable. Similarly, we introduce the system’s $3$-body correlation function, ${G_{3}(\mathbf{w}_{1},\mathbf{w}_{2},\mathbf{w}_{3})}$, as

[TABLE]

Finally, we introduce the $4$-body correlation function, ${G_{4}(\mathbf{w}_{1},\mathbf{w}_{2},\mathbf{w}_{3},\mathbf{w}_{4})}$, as

[TABLE]

The best way to check for the sanity of the previous definitions is to compute the normalisation of the correlation functions, and their scaling w.r.t. $N$. Integrating Eqs. (38), (39) and (40) w.r.t. their phase space coordinates, one obtains

[TABLE]

where we used the shortening notation ${\mathrm{d}1=\mathrm{d}\mathbf{w}_{1}}$. Owing to these scalings, one can therefore use the correlation functions to perform perturbative expansions w.r.t. the small parameter ${1/N}$.

The next step of the calculation is now to inject the previous decompositions into the three first equations of the BBGKY hierarchy, as given by Eq. (37), in order to obtain the evolution equations for ${\partial F_{1}/\partial t}$, ${\partial G_{2}/\partial t}$ and ${\partial G_{3}/\partial t}$. Such equations can be cumbersome to derive, and were obtained using computer algebra in Fouvry et al. (2019b). Writing the system’s $1$-body DF as ${F=F_{1}}$, its time evolution is given by

[TABLE]

The second equation of the hierarchy, for ${\partial G_{2}/\partial t}$, reads

[TABLE]

where we introduced the symmetrising notation

[TABLE]

Finally, the third equation of the hierarchy, for ${\partial G_{3}/\partial t}$, reads

[TABLE]

Here, similarly to Eq. (44), we introduced the symmetrising notation

[TABLE]

where it is assumed that the function $G$ is symmetric w.r.t. its two last indices, i.e. one should have ${G(1,2,3)=G(1,3,2)}$, as is the case for the term in Eq. (45).

Equations (42), (43) and (45) are the starting blocks to obtain a self-consistent set of coupled evolution equations describing the system’s entire dynamics up to order ${1/N^{2}}$.

Appendix C Truncating the BBKGY hierarchy

In this Appendix, we detail how one may truncate Eqs. (42), (43), and (45) to lay the groundwork to derive the closed kinetic equation presented in Eq. (4).

The first step of these simplifications is to perform a truncation at order ${1/N^{2}}$ of these three evolution equations. In Eq. (42), we also note that the collision term for ${\partial F/\partial t}$ only involves ${G_{2}(1,2)}$, whose norm scales like ${1/N}$, as given by Eq. (41). As a consequence, if one aims at deriving a kinetic equation at order ${1N^{2}}$, it is essential to account for the corrections of order ${1/N^{2}}$ that can arise in $G_{2}$. To perform the truncation at order ${1/N^{2}}$, we therefore introduce explicitly the small parameter ${\varepsilon=1/N}$. Following the definition ${\mu=M_{\mathrm{tot}}/N}$, and the scalings from Eq. (41), we perform the replacements

[TABLE]

Using this rewriting, we then keep in the evolution equations only terms up to order $\varepsilon^{2}$. Moreover, owing to the split of $G_{2}$ in two components, we can split the associated evolution Eq. (43) in two components, namely at order ${1/N}$ (resp. ${1/N^{2}}$) that will govern the dynamics of ${\partial G_{2}^{(1)}/\partial t}$ (resp. ${\partial G_{2}^{(2)}/\partial t}$).

A subsequent simplification arises from the homogeneous assumption, i.e. the assumption that system’s mean DF remains a function of $v$ only. This allows us to get rid of the phase mixing term, ${v_{1}\partial F/\partial\theta_{1}}$, in Eq. (42), and also get rid of all the mean-field potential components, i.e. terms involving ${\!\int\!\mathrm{d}2F(2)\,U^{\prime}(\theta_{1}\!-\!\theta_{2})=0}$.

In order to ease the derivations of the kinetic equation, we also assume that the mean system is sufficiently dynamically hot for collective effects to be negligible444In the context of ${1/N}$ dynamics, such an assumption gets the Balescu-Lenard equation to reduce to the Landau kinetic equation (Heyvaerts, 2010; Chavanis, 2012a).. Such an assumption amounts to neglecting the backreaction of a correlation function on the perturbating potential in which that same correlation function is evolving. As a result, we perform the following simplifications

[TABLE]

Finally, we perform three last approximations: (i) in the evolution equation for ${\partial F/\partial t}$, we may neglect the contribution from ${G_{2}^{(1)}}$ that is responsible for the usual ${1/N}$ Landau equation, which identically vanishes for ${1D}$ homogeneous systems, (as highlighted in Eq. (3)); (ii) in the evolution equation for ${\partial G_{2}^{(2)}/\partial t}$, we may safely neglect the contributions from the source term proportional to ${\mu\partial G_{2}^{(1)}(1,2)/\partial v_{1}}$, as one can check that it does not contribute to the kinetic equation (see Fouvry et al. (2019b)); (iii) in the evolution equation for ${\partial G_{3}/\partial t}$, we neglect the contributions from the term proportional to ${G_{2}^{(1)}\!\times\!G_{2}^{(1)}}$, as it will lead to a collision operator proportional to $F^{4}$, while the other source terms will lead to a collision operator proportional to $F^{3}$, that dominates for sufficiently dynamically hot systems.

Following these various truncations and simplifications, we now have at our disposal a set of four coupled evolution equations, that jointly describe the long-term dynamics of the considered system at order ${1/N^{2}}$. The dynamics of ${\partial F(1)/\partial t}$ reads

[TABLE]

The dynamics of ${\partial G_{2}^{(1)}(1,2)/\partial t}$ is given by

[TABLE]

while the dynamics of the second-order correction, ${\partial G_{2}^{(2)}(1,2)/\partial t}$ takes the form

[TABLE]

Finally, the dynamics of the $3$-body correlation function, ${\partial G_{3}(1,2,3)/\partial t}$, reads

[TABLE]

All together, Eqs. (49), (50), (51), and (52) form the starting point to derive the kinetic Eq. (4), as we describe in Appendix D.

Appendix D Deriving the kinetic equation

In this Appendix, we detail the protocol followed to obtain the ${1/N^{2}}$ kinetic equation put forward in Eq. (4), following an approach similar to Rocha Filho et al. (2014). Here, we only present the overall approach and the key steps, while the detailed effective (and cumbsersome) computations were performed using symbolic calculations in Mathematica, as detailed in Fouvry et al. (2019b). From the technical point of view, the main difficulty is to deal, without mistake, with the large number of terms that appear in the successive resolutions of the evolution equations, hence the need for a numerical implementation of this calculation.

Luckily, the four truncated evolution equations, Eqs. (49), (50), (51), and (52) form a closed and well-posed hierarchy of coupled partial differential equations. In particular, owing to the absence of any collective effects, that would require for the explicit characterisation of the system’s linear response, the evolution equations can easily be solved in sequence. The first step is to solve for the time evolution of ${G_{2}^{(1)}(1,2)(t)}$, as governed by Eq. (50). This explicit solution may then be used as a (time-dependent) source term in Eq. (52) to obtain the time evolution of ${G_{3}(t)}$. This function can then be used as a (time-dependent) source term in Eq. (51) to derive the time evolution of ${G_{2}^{(2)}(t)}$.

In each of these three steps, we rely on two main assumptions: (i) Bogoliubov’s ansatz, so that we may take ${F(1,t)=\text{cst.}}$ when solving for the time-evolution of a correlation function; (ii) we neglect the transients associated with any specific initial conditions in the system’s correlations, i.e. we solve these differential equations with the initial conditions ${G_{2}^{(1)}(t\!=\!0)\!=\!G_{3}(t\!=\!0)\!=\!G_{2}^{(2)}(t\!=\!0)\!=\!0}$. Moreover, in order to easily deal with phase mixing terms of the form ${v_{1}\partial G_{2}^{(1)}/\partial\theta_{1}}$, we perform Fourier developments of all the correlation functions w.r.t. their $\theta$-dependence. Similarly, the interaction potential is also expanded in its Fourier harmonics. As imposed by Eq. (2), in the present case of the HMF model, the interaction potential takes the simple form ${U(\theta)\!=\!\sum_{k=\pm 1}\!\tfrac{-U_{0}}{2}\mathrm{e}^{\mathrm{i}k\theta}}$, so that only the harmonics ${k\!=\!\pm 1}$ can support the interaction, which offers a drastic reduction in the total number of resonant terms that can contribute to the system’s dynamics.

Following these three successive resolutions, we now have at our disposal an explicit solution for the time dependence of ${G_{2}^{(2)}(t)}$. Owing to Bogoliubov’s ansatz, we may then consider the limit ${t\to+\infty}$ of that expression, in order to obtain the asymptotic behaviour of $G_{2}^{(2)}$, and inject it in Eq. (49) to obtain the closed ${1/N^{2}}$ collision operator driving the long-term evolution of ${\partial F/\partial t}$. At this stage, a typical time integral appearing in the expression of ${G_{2}^{(2)}(t)}$ takes the form

[TABLE]

where the frequencies $\omega_{1}$, $\omega_{2}$, and $\omega_{3}$ are some linear functions of the velocities $v_{1}$, $v_{2}$, and $v_{3}$, i.e. the resonances involved in the dynamics, while some additional gradients w.r.t. the velocities can get intertwined with the time integrals. Now, our goal is to estimate the asymptotic limit ${t\to+\infty}$ of that expression in order to estimate the collision operator driving the dynamics of ${\partial F/\partial t}$. To do so, we use the asymptotic formula

[TABLE]

with ${\delta_{\mathrm{D}}(\omega)}$ the Dirac delta, and ${\mathcal{P}(1/\omega)}$ the Cauchy principal value (see the expression for ${\delta_{+}(x)}$ in Eq. (6.40) of Balescu (1997)). For nested integrals as in Eq. (53), we apply consecutively the formula from Eq. (54)555One could be concerned by the nested bounds from the three successive integrals of Eq. (53). Even if one has ${\!\!\int_{0}^{t}\!\!\mathrm{d}t_{1}\!\!\int_{0}^{t_{1}}\!\!\mathrm{d}t_{2}\!\!\int_{0}^{t_{2}}\!\!\mathrm{d}t_{3}\!=\!\tfrac{1}{6}\!\!\int_{0}^{t}\!\!\mathrm{d}t_{1}\!\!\int_{0}^{t}\!\!\mathrm{d}t_{2}\!\!\int_{0}^{t}\!\!\mathrm{d}t_{3}}$, when applying successively the formula from Eq. (54), the ${\tfrac{1}{6}}$ volume prefactor does not have to be accounted for.. Doing so, one still prevents for now the evalutation of the gradients w.r.t. the velocities, so that such gradients would only act on the Dirac deltas and the Cauchy principal values.

Once all the time integrals have been replaced by their asymptotic behaviours, the derived kinetic equation takes the form

[TABLE]

where we introduced the shortening notations ${\mathcal{P}=\mathcal{P}(\tfrac{1}{v_{1}-v_{2}})}$, ${\mathcal{P}^{\prime}=\mathcal{P}^{\prime}(\tfrac{1}{v_{1}-v_{2}})}$, and ${\mathcal{P}^{\prime\prime}=\mathcal{P}^{\prime\prime}(\tfrac{1}{v_{1}-v_{2}})}$. When making the substitution ${(v_{1}\leftrightarrow v_{2})}$ in Eq. (55), it is important to note that ${(\mathcal{P},\mathcal{P}^{\prime},\mathcal{P}^{\prime\prime})\to(-\mathcal{P},\mathcal{P}^{\prime},-\mathcal{P}^{\prime\prime})}$. Finally, in Eq. (55), we also introduced the differential operators

[TABLE]

and

[TABLE]

We note that Eq. (55) is almost identical to the ${1/N^{2}}$ kinetic equation already put forward in Eq. (23) of Rocha Filho et al. (2014) for the same physical system. The differences are some corrections in the overall prefactor, and the overall sign of the ${(v_{1}\leftrightarrow v_{2})}$ term.

Luckily, the result from Eq. (57) can be significantly simplified, by using integration by parts, as well as the parity symmetries of the Dirac deltas, the principal values and their derivatives, leading to the final result from Eq. (4). The detailed steps for these calculations can be found in Fouvry et al. (2019b). We briefly present them below for completeness.

The key step is to perform an integration by parts in Eq. (57) w.r.t. the integration variable $v_{2}$, using the formula

[TABLE]

At this stage, the derivatives of the Dirac deltas that appear are subsequently integrated using an integration by parts w.r.t. the integration variable $v_{3}$, so that

[TABLE]

and similarly for ${\delta_{\mathrm{D}}^{\prime}(2v_{2}-v_{1}-v_{3})}$. Proceeding that way allows us not to create any higher order derivatives of the Cauchy principal values. The kinetic equation then becomes simpler, as it reads

[TABLE]

In that equation, we introduced the differential operators

[TABLE]

and

[TABLE]

with the resonance vectors $\mathbf{k}_{1}$ and $\mathbf{k}_{2}$ already defined in Eq. (5).

At this stage, we finally note that the term in ${\mathcal{P}\mathcal{P}^{\prime}}$ in Eq. (60) will not contribute to the dynamics. Indeed, from Eq. (61), we note that $M_{\mathrm{I}}$ is invariant under the change ${(v_{2}\leftrightarrow v_{3})}$. This symmetry can be leveraged to get rid of this term, as follows. Owing to the presence of the double integral ${\!\int\!\mathrm{d}v_{2}\mathrm{d}v_{3}}$, one can perform the symmetrisation ${(v_{2}\leftrightarrow v_{3})}$ for that term. This leaves the Dirac delta, ${\delta_{\mathrm{D}}(\mathbf{k}_{1}\cdot\mathbf{v})}$, invariant. From that same resonance condition, we note that the arguments of the Cauchy principal values are transformed as ${\tfrac{1}{v_{1}-v_{2}}\!\to\!\tfrac{1}{v_{1}-v_{3}}\!=\!-\tfrac{1}{v_{1}-v_{2}}}$. Given the parities of $\mathcal{P}$ and $\mathcal{P}^{\prime}$, we can therefore conclude that the term in ${\mathcal{P}\mathcal{P}^{\prime}}$ in Eq. (60) is antisymmetric under the change ${(v_{2}\leftrightarrow v_{3})}$ so that the overall contribution of this term vanishes. The last step of the calculation is finally to perform the replacement {\!\int\!\mathrm{d}v_{2}\,\big{(}\mathcal{P}^{\prime}\big{)}^{2}\!\to\!\mathcal{P}\!\int\!\mathrm{d}v_{2}/(v_{1}\!-\!v_{2})^{4}}. All in all, one finally obtains the closed kinetic equation spelled out in Eq. (4).

Appendix E Linear response theory

As highlighted in Appendix C, in order to obtain the kinetic Eq. (4), we had to neglect the contributions from collective effects in the BBGKY evolution equations. As a consequence, in order to test that kinetic equation, it is essential to place ourselves in regimes where collective effects are indeed unimportant. Luckily the strength of the self-gravitating amplification can be directly estimated by solving the linear response theory of the system. This is what we briefly reproduce in that Appendix.

A systematic approach to perform this calculation is to proceed by analogy starting from the generic result regarding the linear stability analysis of (inhomogeneous) long-range interacting systems. Following Eq. (5.94) of Binney and Tremaine (2008) (similarly Eq. (G3) of Fouvry and Bar-Or (2018) and references therein), a system’s stability is governed by the response matrix

[TABLE]

where the angle-action coordinates are ${(\theta,J)=(\theta,v)}$, the orbital frequencies are ${\Omega(J)=v}$. A system is then linearly unstable if there exists a complex frequency ${\omega=\omega_{0}\!+\!\mathrm{i}\eta}$, with ${\eta>0}$, such that ${\widehat{\mathbf{M}}(\omega)}$ has an eigenvalue equal to $1$. The natural basis elements follow from the pairwise interaction from Eq. (2), that can be written under the separable form

[TABLE]

and it is straightforward to check that ${(\psi^{(+)},\psi^{(-)})}$ indeed form a biorthogonal basis, as defined, e.g., in Eq (G1) of Fouvry and Bar-Or (2018). Similarly their Fourier transform w.r.t. the angle can easily be computed. It is independent of the action $v$, and reads ${\psi^{(p)}_{k}=\delta_{k}^{p}\sqrt{U_{0}/2}}$. Owing to this Kronecker delta, the ${2\times 2}$ response matrix from Eq. (63) is then diagonal. We may finally introduce the susceptibility matrix (or dielectric function), ${\bm{\varepsilon}=\mathbf{I}-\widehat{\mathbf{M}}}$, that is also diagonal with the coefficients

[TABLE]

Let us emphasise that the result from Eq. (65) is identical to the result presented in Eq. (9) of Chavanis and Delfini (2009).

Relying on the same de-dimensionalisation as in Eq. (27), one can rewrite the susceptibility coefficient from Eq. (65) under the form

[TABLE]

with ${\overline{\omega}=\omega t_{\mathrm{d}}}$ a dimensionless frequency. We also recall that the dimensionless PDF, ${\overline{F}(u)}$, and the stability parameter, $Q$, were given in Eqs. (25) and (26).

Luckily, it is possible to further precise our characterisation of the system’s stability for the range of DFs that will be considered in this paper. Following the generic shape of the test PDF from Eq. (70), let us therefore assume that the system’s DF, ${F(v)}$, is single-humped, i.e. it possesses a single maximum. We will also assume that the DF’s is even, so that this maximum is reached in ${v=0}$. Owing to this parity, in Eq. (66), we note that ${\varepsilon_{+}(\overline{\omega})=\varepsilon_{-}(\overline{\omega})}$, so that we may limit ourselves to only studying $\varepsilon_{+}(\overline{\omega})$.

Because the DF is single-humped in ${v=0}$, following Nyquist’s criterion (see, e.g., Section 2.6 in Chavanis and Delfini (2009)), such a DF is linearly stable if, and only if, ${\varepsilon_{+}(0)>0}$. As a consequence, following Eq. (66), such a DF is linearly stable if, and only if, one has

[TABLE]

Conveniently, it is straightforward to compute these stability thresholds for the generic test DF from Eq. (70). This DF is parametrised by the power index $\alpha$, such that ${\alpha=2}$ corresponds to the Gaussian distribution. One finds

[TABLE]

In summary, a homogeneous Gaussian PDF (i.e. a homogeneous Boltzmann distribution) is therefore linearly stable if, and only if, it satisfies ${Q>1}$.

For more generic PDFs, the susceptibility coefficient from Eq. (66) cannot always be computed analytically. One has to resort to numerical evaluations, e.g., following the method presented in Appendix D of Fouvry et al. (2019a). We illustrate this method in Fig. 4 by representing the Nyquist contours associated with the test PDF from Eq. (70).

As expected, we recover in Fig. 4 that the larger is $Q$, the weaker is the system’s susceptibility, i.e. dynamically hot systems are less efficient at amplifying perturbations through collective effects. Figure 4 also illustrates that the test case considered in Fig. 2 is sufficiently hot for the system to be linearly stable. In addition, collective effects should prove sufficiently negligible for the kinetic Eq. (4) to be in its applicability regime.

Appendix F Numerical applications

In this Appendix, for completeness, we briefly present our numerical approach to perform $N$-body simulations of the HMF model. Following Eq. (1), the specific Hamiltonian of a test particle embedded in that system reads

[TABLE]

where we introduced the system’s instantaneous magnetisations as ${M_{x}(t)\!=\!U_{0}\mu\sum_{i=1}^{N}\cos(\theta_{i}(t))}$, as well as ${M_{y}(t)\!=\!U_{0}\mu\sum_{i=1}^{N}\sin(\theta_{i}(t))}$. Two important remarks should be made w.r.t. Eq. (69). First, because the pairwise interaction does not diverge at zero angular separation, the test Hamiltonian from Eq. (69) and the associated evolution equations can also be used to obtain correct evolution equations for each of the system’s particles, treating the magnetisations as external, i.e. not taking any derivatives of it. Second, as the magnetisations involve a sum over the $N$ particles, they should be interpreted as global and shared quantities, that need to be computed only once for each timestep. This allows for the computational complexity of integrating for one timestep to scale like ${\mathcal{O}(N)}$, rather than ${\mathcal{O}(N^{2})}$ as in the usual $N$-body problem.

To compute the velocity fields at a given time, we proceed as follows: (i) we compute and store ${(\cos,\sin)}$ for all particles; (ii) we reduce these quantities to compute the instantaneous magnetisations ${(M_{x},M_{y})}$; (iii) we compute the velocity fields ${\mathrm{d}\theta_{i}/\mathrm{d}t}$ and ${\mathrm{d}v_{i}/\mathrm{d}t}$. Owing to the fact that the Hamiltonian from Eq. (69) is separable, particles are then advanced using a fourth-order symplectic integrator (see Eq. (4.6) in (Yoshida, 1990)). The numerical simulations presented in Section IV were all performed using an integration timestep equal to ${\delta t=1/(2\sigma)}$ that guaranteed a relative error in the total energy of the order of $10^{-7}$.

Following (Nardon and Pianca, 2009), the initial distribution of the system is taken to be a generalised Gaussian distribution. For a given index $\alpha$ and velocity dispersion $\sigma$, its PDF reads

[TABLE]

This PDF satisfies the normalisation condition ${\!\int\!\mathrm{d}u\,P(u)=1}$, is of zero mean and variance $\sigma^{2}$. For ${\alpha=2}$, this corresponds to the Gaussian distribution, while larger values of $\alpha$ are associated with less peaked distributions. Luckily, the PDF from Eq. (70) can also easily be sampled (see Eq. (9) in Nardon and Pianca (2009)).

To measure fluxes, as in Figs. 2 and 3, we proceeded as follows. For each setup, we perform $N_{\mathrm{real}}$ realisations, only changing the initial conditions. The dimensionless velocity space, ${u\in[-3,3]}$, is truncated in $50$ equal size bins. For each realisation, each velocity bin location, and each timestep, we compute the proportion of particles left to that location, subsequently averaged over all the available realisations. For each velocity bin, the associated time series are then fitted with a linear time-dependence, whose slope is the local diffusion flux, ${\overline{\mathcal{F}}(u,t\!=\!0)}$. To estimate the associated measurement errors, we follow the exact same approach for $N_{\mathrm{boot}}$ measurements, except that the sample of $N_{\mathrm{real}}$ realisations over which the ensemble average is performed allows for repetitions, i.e. for the same realisations to be used more than once. The measurement is then given by the median value, while the errors are given by the ${16\%}$ and ${84\%}$ confidence levels. In Figs. 1 and 2, we used the values ${N=10^{3}}$, ${N_{\mathrm{real}}\!=\!N_{\mathrm{boot}}\!=\!10^{3}}$, and simulated the systems up to ${t=5\!\times\!10^{8}t_{\mathrm{d}}}$. With our implementation, running one such realisation asked for about ${18\mathrm{h}}$ of computation on a single core. In Fig. 3, we used the values ${0.5\!\times\!10^{3}\!\leq\!N\!\leq\!1.7\!\times\!10^{3}}$, ${N_{\mathrm{real}}\!=\!N_{\mathrm{boot}}\!=\!224}$, and simulated the systems up to ${t=5\!\times\!10^{8}t_{\mathrm{d}}}$.