Coarse-grained collisionless dynamics with long-range interactions

We derive an effective evolution equation for a coarse-grained distribution function in the case of long-range-interacting systems by performing a coarse graining that preserves the symplectic structure of the collisionless Boltzmann, or Vlasov, equation obeyed by the fine-grained distribution function. We first derive a general form of such an equation based on symmetry considerations and very general assumptions only. We then restrict ourselves to the case of one-dimensional systems, performing the coarse-graining and obtaining an explicit version of the equation. Finally, we use such an equation to predict the dependence of the damping times on the coarse-graining scale and check them against numerical results for the Hamiltonian Mean Field (HMF) model.

Introduction. Long-range interactions, whose potential energy decays with the distance r between the interacting bodies slower than r −d , where d is the dimension of space [1], are relevant to astrophysics and plasma physics [2,3], since gravitational and unscreened Coulomb forces are long-ranged, as well as to condensed matter, given that dipolar interactions in three dimensions or effective interactions mediated by the electromagnetic field in systems of cold atoms in a cavity [4,5] are long-ranged, and occur also in two-dimensional fluids [6]. Systems with longrange interactions exhibit peculiar features both at equilibrium and out of equilibrium [1,7,8]. The dynamics of such systems is dominated by collective, mean-field effects, rather than by binary collisions; as a consequence, the relaxation time towards equilibrium τ R diverges with the number of particles N , as already noted by Chandrasekhar for self-gravitating systems [2,9]. In the N → ∞ limit, or for times t < τ R for a large but finite long-range-interacting system, the dynamics obeys the non-collisonal Boltzmann equation [2], also referred to as the Vlasov equation [1]. Let us consider for simplicity a system of particles with unit masses and introduce the single-particle Hamiltonian where q = (q 1 , . . . , q d ), p = (p 1 , . . . , p d ), U is the self-consistent interaction potential U (q, t) = dp dq f (q , p , t)V (|q − q |) , f (q, p, t) is the single-particle distribution function defined on the 2d-dimensional phase space, V (r) is the pair potential energy between two particles at distance r and dp dq = d i=1 dp i dq i . This allows us to make explicit the symplectic structure underlying the Vlasov equation, that can be written as where {·, ·} is the Poisson bracket on the single-particle phase space. The Hamiltonian character of the Vlasov equation (3) has many important consequences [10][11][12]. Here we just stress that Eq. (3) is time-reversal-invariant and its dynamics is constrained by an infinite number of conservation laws: the Casimirs, i.e., the functionals are conserved for any choice of the function C(f ). Remarkably, the Boltzmann entropy is nothing but a particular Casimir, corresponding to C(f ) = −f ln f , so that it is a constant of motion and no H theorem holds. All these properties seem to suggest that no relaxational dynamics is possibile: any time dependence of f should survive forever in the N → ∞ limit, and at least up to t ≈ τ R when collisional effects set in for a large but finite system.
Numerical results depict a totally different scenario: starting from a generic initial condition, a given observable exhibits oscillations that damp out on a rather fast time scale not depending on N (at variance with τ R ) until it attains a nearly constant value. The paradigmatic example is gravitational collapse [13][14][15], where the relevant observable is either the gravitational radius or the virial ratio, so that the damped oscillations are termed "virial oscillations", and this purely noncollisional relaxation is referred to as "violent relaxation" [16]. Violent relaxation is a universal phenomenon, occurring in any long-range-interacting system, and the state reached by the system after this process is referred to as a quasi-stationary state (QSS), since in a finite system it will eventually relax to equilibrium for t > τ R [1]. Despite many advances in the last decades [1,7,8], a theory able to predict these states given a generic initial condition is still missing. It is widely believed that the mechanism at the basis of violent relaxation is similar to Landau damping [12,17,18]. Basically this means that the Vlasov dynamics does never actually stop: rather it trickles down towards smaller and smaller scales until it no longer affects the behavior of any reasonable, coarse-grained, observable. This agrees with the results found in [19], where it was proven that, given a coarse-grained distribution functionf , defined by averaging out all the details below some finite scale Λ in phase space, and any convex function C(x), the corresponding Casimir C[f ] will decrease in time. Despite this, a convincing quantitative picture of this process is still missing: our aim is then to contribute to the filling of this gap by providing an effective evolution equation for the coarse-grained distribution functionf . First, we propose a general form of this equation, up to coefficients, based on symmetry considerations only. Then we show how to actually derive its complete expression in the case of d = 1 bounded (or periodic) systems, by explicitly performing the coarse graining, and use such an equation to predict the dependence of damping times on the coarse-graining scale, checking the results against numerical simulations of a simple model. Symplectic coarse graining. We aim at writing an effective evolution equation for the coarse-grained distribution functionf . Many properties of such an equation can be derived from symmetry considerations and very general assumptions, that define what we will refer to as symplectic coarse graining. First of all, we observe that iff is normalized to unity then a coarse-grained single-particle Hamiltonian H[f ] is defined too as in Eq. (1), withf in place of f ; to ease the notation, from now on we simply write H in place of H[f ]. We then assume that the coarse graining procedure does not depend on the particular choice of the canonical coordinates, i.e., it preserves the symplectic structure. This means that the dynamical evolution off can be expressed in terms of Poisson brackets. Moreover, we assume that such Poisson brackets contain functions of H andf alone and are linear inf ; physically, this means that the particles interact only via the mean field Hamiltonian, as in the Vlasov equation (3). As a consequence of these assumptions, the evolution equation forf can be written as where L H (f ) is an operator depending on H and acting linearly onf , whose most general form is a linear combination of nested Poisson brackets wheref appears only once, that is, of terms of the form λ 1 (H), λ 2 (H), · · · λ k (H),f · · · , where the λ k 's are generic functions of H. By repeatedly using the identity {λ k (H), ·} = λ k (H) {H, ·} and by denoting with {H, ·} nf the n nested Poisson brackets H, H, · · · H,f · · · n times we may thus write where the µ n 's are generic functions that absorb the coefficients of the linear combination. In Eq. (6) we have highlighted the first term of the sum, assuming µ 1 (H) ≡ 1, as is reasonable sincef → f and Eq. (6) must reduce to Eq. (3) when 1 Λ → 0. We note that both the normalization off and the total energy E[f ] are conserved by Eq. (6), as required by a physically sound evolution and as shown in the Supplemental material. Equation (6) is thus the most general outcome of symplectic coarse graining. The non-Vlasov-like contributions to this equation (i.e., the terms of the sum on the r.h.s.) can be split into two classes: the even and the odd ones (that is, containing an even/odd number of Poisson brackets). The odd terms do not break time-reversal invariance, so that they renormalize the time-reversible Vlasov evolution, while the even terms break the time-reversal invariance, so that they may account for the dissipative character of the coarse-grained evolution. However, we expect the physical evolutions to be a subset of all the possible ones described by Eq. (6), i.e., that not all the µ n 's are admissible. For instance, as already mentioned, all the convex Casimirs defined by the coarse-grained distribution functionf must be decreasing functions of time. It is not easy to impose such a constraint on Eq. (6), but it is interesting to note that the lowest-order truncation of the latter equation, obtained by setting µ n = 0 ∀ n > 2, with the additional constraint µ 2 (x) > 0 ∀ x, does satisfy this constraint (see Supplemental material), and actually describes a Vlasov-like evolution with added diffusive effects, hence admitting a relaxational behavior. Indeed, {H, ·} is proportional to the directional derivative along the Hamiltonian flow generated by H, so that {H, {H, ·}} is a sort of anisotropic Laplacian and the second term on the r.h.s. of Eq. (7) describes a diffusion taking place along the Hamiltonian flow, whose strength depends on µ 2 , that in turn will depend on the coarse-graining scale. This will become apparent in the one-dimensional case that we are going to tackle in the following.
One-dimensional systems. Let us now show how to perform a symplectic coarse graining of the Vlasov equation (3) and obtain an explicit evolution equation for the coarse-grainedf in the case of systems living in one dimension, bounded in space or with periodic boundary conditions on the linear coordinate. In this case the single-particle Hamiltonian has one degree of freedom. Therefore, at a given time, the Hamiltonian is integrable and a canonical transformation (p, q) → (J, ϑ) exists, where (J, ϑ) are action-angle variables. Being the single-particle Hamiltonian time-dependent in general, such a transformation leads to a Hamiltonian independent of the angle ϑ only at a given time t. The instantaneous flow will be such as to keep J constant, and the angle will linearly evolve in time, where ω(J) = dH/dJ. Let us consider a (small) interval of actions ∆J = J 2 − J 1 and define the average frequency ω over ∆J, and consequently δω = ω − ω, and a distribution function coarse-grained along J as where J is such that ω = ω(J). To get a truly coarse-grained distribution function one should average also over an interval of angles ∆ϑ, but it is better to consider such an average as carried over a time interval ∆t, that is, to assume we are blind to changes of the coordinates of the particles occurring on time scales smaller than ∆t. This means that we neglect the time dependence of the Hamiltonian on a time scale ∆t, and we can use action-angle coordinates for times between t and t + ∆t, defining a non-constant coarse-graining scale on ϑ, namely ∆ϑ = ω∆t.
Our coarse-grained distribution functionf is then the functionf given by Eq. (9), further averaged over an interval of angles of width ∆ϑ centered in ϑ. As a consequence, we are not able to distinguish any change off on scales smaller than the coarse-graining volume ∆Γ = ∆ϑ∆J. For times between t and t + ∆t we have approximated the flow in phase space with a stationary one, so that its evolution operator should be written as U ∆t = e ∆t{H,·} . The latter is not constant over ∆Γ, but within such volume we can consider J and ϑ as uniformly distributed random variables, so that (up to very unlikely initial conditions) we can replace the evolution operator U ∆t with the coarse-grained onẽ The evolution dictated by Eq. (10) satisfies the constraint on the evolution of convex Casimirs and can be translated into a differential equation for the coarse-grained distribution functionf (t) (see Supplemental material) as where we have introduced the diffusion coefficients D n (J) = (∆t) n−1 n! κ n (δω) (12) in terms of the cumulants κ n of the probability distribution of δω. Eqs. (11) and (12) are written as such only at the time t chosen to define the action-angle coordinates. However, this problem can be overcome by rewriting the equation in a coordinate-independent way, by noting that −ω(J)∂ ϑ = {H, ·}, so that Eq. (11) is nothing but Eq. (6) with the coefficients µ n (H) explicitly given 2 as µ n = D n (J(H))[ω(J(H))] −n . Note that Eq. (12) implies that all the diffusion coefficients vanish if ω does not depend on J, so that the coarse-grainedf obeys the Vlasov equation as the fine-grained f does. This is coherent with our picture, because no randomness is present if ω does not depend on J and all the particles coherently drift in ϑ with the same frequency. Indeed, in the harmonic case where H is linear in J and ω is constant the motion can be described in terms of normal coordinates without any damping. As already discussed, the even terms are those responsible for the breaking of the time-reversal symmetry; we can estimate their order of magnitude as while the odd coefficients are even more suppressed with n, being the distribution of δω even at the leading order. Hence, as long as ∆J and ∆t are not too large, only the very first terms of the sum in Eq. (11) will give a non-negligible contribution. To the leading order in ∆J and ∆t the evolution equation forf becomes a Fokker-Planck equation (see Supplemental material) where the leading-order diffusion coefficient is and (14) can be cast in the covariant form that is a special case of Eq. (7).
Scaling of damping times. The choice of the coarse graining scale is crucial to obtain an effective description of the observable we are interested in. A small coarse graining scale would be of no utility if we are interested in an observable whose characteristic scale is large, because the resulting effective equation will describe its evolution essentially in the same way as the original Vlasov equation. According to the lowest-order truncation of the effective evolution equation, i.e., the Fokker-Planck equation (14), the characteristic damping time of the coarse-grained distribution functionf is Let us now ask how τ depends on the coarse graining scale. In order for our coarse graining to be independent of the choice of the coordinates in phase space we have to define its scale in terms of phase space volumes, that are invariant under canonical transformations. Let ∆Γ be the coarse graining volume in action-angle variables; we have ∆Γ = ∆J∆ϑ, so that we may assume that 3 ∆J ∝ √ ∆Γ and ∆ϑ ∝ √ ∆Γ. Equation (15) contains ∆t instead of ∆ϑ, because of the way we have constructed our coarse-grained equation; however, ∆t ∝ ∆ϑ ∝ √ ∆Γ, so that Eqs. (15) and (17) imply τ ∝ (∆Γ) −3/2 . If S is the smallest volume (i.e., surface, since our phase space is two-dimensional) we want to probe, neglecting the dynamics occurring on scales smaller that S, we have to choose ∆Γ ≈ S, so that the damping time at scale S obeys the scaling relation The time given by (18) is the time after which the dynamics of the fine-grained f has moved to scales smaller than S in phase space, so that our coarse-grained description is no longer able to detect it. One way to probe different scales is to look at the Fourier modes of the distribution function in phase space. Defining the wave vector k = (k p , k q ), where k q and k p are its components along any couple of canonical coordinates (p, q), the k-th mode of the distribution function is and probes a strip in phase space of width proportional to k −1 , where k = k 2 p + k 2 q . Therefore, to describe the evolution of f k we have to choose S ∝ k −1 , so that we expect to see the damping of f k on a time scale To check the scaling law (20) we have considered one of the most studied models with long-range interactions, the Hamiltonian Mean Field (HMF) model [20]. We have performed numerical simulations of the dynamics of an HMF model with N = 2 × 10 5 particles, computing the evolution of the f k starting from nonstationary configurations (see Supplemental material for further details on the simulations and for results obtained with different initial conditions). In Fig. 1 we plot the envelope of the real part of f k minus its averaged value in the quasi-stationary state reached after the violent relaxation f QSS k as a function of the rescaled time t/k 3/2 . It is apparent that rescaled damping time is around 20; moreover, we obtain a very good data collapse, meaning that not only the damping times but actually the entire dynamics, for not too short times, follows the scaling law (20).
Conclusions. We have presented a derivation of an effective evolution equation for a coarse-grained distribution function in the case of long-range-interacting systems, whose dynamics obeys the collisionless Boltzmann (or Vlasov) equation in the limit of a large number of particles. A general form of the equation has been given based on symmetry considerations, in particular requiring the conservation of the symplectic structure, and an explicit derivation of the complete equation has been presented for the bounded (or periodic) one-dimensional case. The effective equation has a zero-order Vlasov term, and two kind of additional terms: those that break the time-reversal symmetry, thus describing damping, and those that renormalize the Vlasov dynamics. We note that the same kind of structure of an effective equation for the moments of the distribution function had been found, at the leading order and based on heuristic considerations, in [21]. The lowest-order non-time-reversal-invariant term is a diffusion along the Hamiltonian flow: the idea of a diffusion along the J = constant lines in phase space had already been exploited for the HMF model in [22], again based on heuristic considerations, but only in the stationary case. Our results provide a solid quantitative picture of the mechanism underlying violent relaxation in long-range-interacting systems and shed light on the rôle of the coarse-graining scale: numerical results for the HMF model are in very good agreement with our predictions. As recently discussed in [23], a relaxation occurring on time scales smaller than the collisional one might also be induced by a finite number of particles N ; however such mechanism seems not to be relevant to violent relaxation, that occurs in finite systems as well as in the Vlasov N → ∞ limit, and whose timescale does not depend on N . Although we have explicitly derived the effective equation only in the one-dimensional case, we expect the extension of our procedure to systems of higher dimensionality that are integrable at a given time, like e.g. the self-gravitating case with imposed spherical symmetry, to be possible. Moreover, one may think of applying a truncated form of the general equation, as for instance the lowest-order given by Eq. (7), supplemented by some ansatz for the unknown function µ 2 , to describe violent relaxation in generic long-range-interacting systems. Acknowledgment. This work is part of MIUR-PRIN2017 project Coarse-grained description for non-equilibrium systems and transport phenomena (CO-NEST) n. 201798CZL whose partial financial support is acknowledged.

I. INTRODUCTION
We present here proofs of some results put forward in the main paper and details on how the numerical result presented in Fig. 1 of the main paper has been obtained as well as other results on numerical simulations of the Hamiltonian Mean Field (HMF) model.

A. Poisson brackets and phase space integrals
First of all let us recall some properties of Poisson brackets. We consider a Hamiltonian system with d degrees of freedom and Hamiltonian H(q 1 , . . . , q d , p 1 , . . . , p d ). Being by definition that is a divergence in phase space, we have for any f ,g decaying sufficiently fast for large values of coordinates and momenta. In Eq. (22) the integral is extended to the whole 2d-dimensional phase space and we have used the shorthand notation dΓ = dp dq = d i=1 dp i dq i that we shall continue to use from now on. From Eq. (22), considering three functions f , g and h again decaying sufficiently fast at infinity, and applying the Leibnitz rule we get the integration by parts formula

II. GENERAL FORM OF THE EFFECTIVE EVOLUTION EQUATION
The general evolution equation for the coarse-grained distribution function-Eq. (6) in the main paper-can be written as where {H, ·} nf = H, H, · · · H,f · · · n times (26) and the µ n 's are generic functions of H.
where the M n 's are such as M n (H) = µ n (H). Integrating over the whole phase space and using Eq. (22) we have and being ∂ tf dΓ = d dt f dΓ Eq. (29) implies the conservation of the norm.

Conservation of the energy
Working in one dimension to ease the notation, the energy of the system in a state defined by the coarse-grained distribution can be written as The time derivative of the energy is where the M n (H) are such that M n (H) = Hµ n (H); integrating the above equation over the whole phase space and using Eq. (33) we get dE/dt = 0.

B. Time evolution of convex Casimirs
The fine-grained Vlasov evolution has infinite conserved quantities (Casimirs) obtained by integrating a generic function C(f ) of the distribution function f over the whole phase space. Replacing f by a coarse-grained one the Casimirs ar no longer cosntant of motion. However, among all the Casimirs defined using any coarse-graining distribution functionf , that is those corresponding to a convex C (that will be referred to as "convex Casimirs" from now on) must be non-increasing functions of time [19]. In the case of one-dimensional systems, to be considered below, we are able to prove that our version of the coarse-grained dynamics does agree with such a constraint (see §III A below). We did not succeed in proving such a result for the most general form (25) of the coarse-grained evolution equation, but we can show that convex Casimirs do not increase with time if we restrict ourselves to the lowest-order truncation of Eq. (25), that is and using Eq. (27), with µ 2 (H) = M 2 (H), we get The first integral in the r.h.s. of the above equation vanishes due to Eq. (22), while integrating by parts the second term using Eq. (24) we get It is interesting to note that Eq. (39) tells us that C[f ] does not reach its minimum: its evolution eventually stops whenf approaches a stationary solution, that is, such that H,f = 0. The latter is a necessary feature of a consistent evolution, because it would not be possible, in general, to reach a state where all the (infinite) convex Casimirs are simultaneously minimized (see also the discussion in Ref. [19]).

III. ONE-DIMENSIONAL SYSTEMS
In the case of one-dimensional systems, where the effective Hamiltonian is integrable at a given time, it is possible to perform the symplectic coarse graining and obtain an explicit evolution equation for the coarse-grained distribution functionf , that can be written asf where, as in Eq. (10) of the main paper,Ũ where the average over ∆Γ means a double average, first over the action variable on the interval ∆J and then over the angle variable on the interval ∆ϑ that can be expressed in terms of the time scale ∆t as ∆ϑ = ω∆t.

A. Evolution of convex Casimirs
Let us now show that the evolution defined by Eqs.
where we dropped the average over the angle variable since the integrand only depends on J . Considering now a sufficiently small ∆t such as to expand the exponential in the above equation up to first order and applying the resulting evolution operator to the coarse-grained distribution function evaluated at time t,f t , we get To first order accuracy we can replace the derivative in the above equation with the ratio off (ϑ, J)−f (ϑ−ω(J )∆t, J) and ω(J )∆t, so that Eq. (43) becomesf On the other hand, for any convex function C(x) and for any random variable x, so that or, explicitly writing the average over ∆J once again, To obtain a condition on the Casimir functional at time t + ∆t we have to integrate the above relation in ϑ and J all over the phase space, obtaining butf is a periodic function of ϑ, so that that no longer depends on J ; the average over ∆J is thus trivial and Eq. (48) becomes that is, what we wanted to prove.

B. Differential equation forf
Let us now show that the finite-time evolution defined by Eqs. (40) and (41) can be cast into a differential equation forf , that is, Eq. (11) in the main paper. We begin once again by writing the evolution operator in action-angle variables,Ũ ∆t = e ∆t{H,·} ∆Γ = e −ω(J )∆t∂ ϑ ∆J . (51) Then, since operators at the exponent evaluated at different points commute, we can apply the usual cumulant expansion and findŨ where κ n is the n-th cumulant of the probability distribution of the frequencies ω. The time evolution becomes so that Let δω = ω − ω, where ω = ω(J) = κ 1 (ω) is the average of the distribution of the frequencies ω. Then κ 1 (δω) = 0 and κ n (δω) = κ n (ω) for any n > 1, so that we can extract the first term from the sum in Eq. (54) and obtain Eq. (11) of the main paper, that is with the diffusion coefficients D n defined as D n (J) = (∆t) n−1 n! κ n (δω) .
that is, Eq. (14) of the main paper. The latter equation can be cast in covariant form we note that so that Eq. (58) becomes that is, Eq. (16) in the main paper.

D. Langevin equation
It is interesting to note that the leading-order effective evolution equation in the one-dimensional case (60) can be cast in the form of a Fokker-Planck equation and interpreted, in turn, in the corresponding Langevin formalism. The nested Poisson bracket in Eq. (60) can be written as where i = 1, 2, x 1 = q, x 2 = p, v i = ε ij ∂ xi ω with ε ij the totally antisymmetric Levi-Civita symbol and we have used the Einstein summation convention over repeated indices, so that Eq. (60) becomes in which we recognize the general form of a Fokker-Planck equation with a non-isotropic and non-uniform diffusion coefficient. Such an equation is in turn equivalent, in the Langevin formalism, to the Stratonovich differential equation [24]ẋ ξ(t) being a white noise with correlation function Exploiting the definition of v j we can writeq as expected, this couple of equations can be derived from the stochastic Hamiltoniañ where once again we are using the Stratonovich formalism.

IV. NUMERICAL SIMULATIONS
We report here in the following some details on the numerical simulations used to obtain the results shown in Fig.  1 of the main paper as well as some results obtained with different initial conditions, that confirm the scaling of the damping times with the coarse-graining scale found for the results shown in the main paper.

A. HMF model
The Hamiltonian Mean Field (HMF) model [20] has been the workhorse of the studies on long-range-interacting systems in the last decades. Its Hamiltonian is where q i ∈ [−π, π] ∀ i = 1 . . . N and p 1 , . . . , p N are the conjugated momenta. It can be seen either as a system of N globally coupled XY spins or as a systems of N globally coupled particles with unit mass moving on a ring. In the following we shall set J = 1, thus considering only attractive (ferromagnetic, in the spin language) interactions.

B. Details of the simulations
We numerically integrated the equations of motion derived from the Hamiltonian (67) using a third order symplectic algorithm [25]. We always generated initial conditions symmetric in position and momentum space, i.e., such as being such a symmetry conserved by the exact dynamical evolution, following [26] we actually computed the dynamics of N/2 particles only and imposed the symmetry by generating at each timestep the other N/2 positions and momenta such as for i = 1, . . . , N/2. The observables we have evaluated as a function of time, after starting from non-stationary initial conditions, are the Fourier modes of the distribution function, defined as where k = (k q , k p ). A typical time evolution of such an observable is shown in Fig. 2. Here initial conditions are such that the q's are drawn from a uniform distribution in [−1, 1] while the p's are distributed according to a Gaussian with zero mean and standard deviation σ = 0.1. The latter are the same initial conditions yielding the data shown in Fig. 1 of the main paper, where the envelope of f k was plotted as a function of the rescaled time t/k 3/2 after a subtraction of its asymptotic value. The inset of Fig. 2 here shows which values have been selected as envelope ones.

C. Other results
Let us now show some numerical results of the same kind as those shown In Fig. 1 of the main paper, that is, the envelopes of f k plotted as a function of the rescaled time t/k 3/2 after a subtraction of their asymptotic values, but obtained starting from different initial conditions. We always considered N = 2 × 10 5 . The results shown in Fig. 3 were obtained starting from "waterbag" initial conditions, while those reported in Fig. 4 refer to "cold" initial conditions, i.e., where the kinetic energy is zero, as considered in [21]. In both cases the data collapse shows that the scaling law τ ∝ k 3/2 for the damping times predicted in the main paper agrees well with the numerical results, not only for the damping times themselves but also for the whole dynamical evolution.