Violent relaxation in the Hamiltonian Mean Field model: I. Cold collapse and effective dissipation

TL;DR

This paper studies violent relaxation in the Hamiltonian Mean Field model, proposing an approximation scheme that captures damping of collective oscillations and effective dissipation, aligning well with numerical simulations for cold initial conditions.

Contribution

It introduces a novel approximation method based on moments of the distribution function, inspired by the Caldeira-Leggett model, to describe collisionless relaxation in long-range interacting systems.

Findings

Derived a dissipative equation for virial oscillations

Predicted damping and effective potential renormalization

Achieved good agreement with numerical simulations

Abstract

In $N$-body systems with long-range interactions mean-field effects dominate over binary interactions (collisions), so that relaxation to thermal equilibrium occurs on time scales that grow with $N$, diverging in the $N\to\infty$ limit. However, a faster and non-collisional relaxation process, referred to as violent relaxation, sets in when starting from generic initial conditions: collective oscillations (referred to as virial oscillations) develop and damp out on timescales not depending on the system's size. After the damping of such oscillations the system is found in a quasi-stationary state that survives virtually forever when the system is very large. During violent relaxation the distribution function obeys the collisionless Boltzmann (or Vlasov) equation, that, being invariant under time reversal, does not "naturally" describe a relaxation process. Indeed, the dynamics is moved to smaller and smaller scales in phase space as time goes on, so that observables that do not depend on small-scale details appear as relaxed after a short time. We propose an approximation scheme to describe collisionless relaxation, based on the introduction of moments of the distribution function, and apply it to the Hamiltonian Mean Field (HMF) model. To the leading order, virial oscillations are equivalent to the motion of a particle in a one-dimensional potential. Inserting higher-order contributions in an effective way, inspired by the Caldeira-Leggett model of quantum dissipation, we derive a dissipative equation describing the damping of the oscillations, including a renormalization of the effective potential and yielding predictions for collective properties of the system after the damping in very good agreement with numerical simulations. Here we restrict ourselves to "cold" initial conditions; generic initial conditions will be considered in a forthcoming paper.

Figures

Equations

$$\mathcal{H}=\frac{1}{2}\sum _{i=1}^N{p}_i^2-\frac{1}{2N}\sum _{i,j=1}^N\cos \left({\vartheta }_i-{\vartheta }_j\right),$$

$$U\left({\vartheta }_1,\dots ,{\vartheta }_N\right)=\frac{1}{2N}\sum _{i,j=1}^{N}v\left({\vartheta }_i-{\vartheta }_j\right)=-\frac{1}{2N}\sum _{i,j=1}^N\cos \left({\vartheta }_i-{\vartheta }_j\right).$$

$$m_x$$

$$m_y$$

$$U=-\frac{1}{2N}\sum _{i,j=1}^N\cos \left({\vartheta }_i-{\vartheta }_j\right)=-\frac{1}{2}{m}^2,$$

$${\dot{\vartheta }}_i$$

$${\dot{p}}_i$$

$$m_x[f]$$

$$m_y[f]$$

$$U[f](\vartheta )=-{\int }_{-\infty }^{+\infty }d{p}^{\prime }{\int }_{-\pi }^{\pi }d{\vartheta }^{\prime }f({\vartheta }^{\prime },p^{\prime },t)\cos \left(\vartheta -{\vartheta }^{\prime }\right)=-\cos \vartheta ⟨\cos {\vartheta }^{\prime }⟩-\sin \vartheta ⟨\sin {\vartheta }^{\prime }⟩=-m_x[f]\cos \vartheta -m_y[f]\sin \vartheta ,$$

$$\frac{\partial f}{\partial t}+p\frac{\partial f}{\partial \vartheta }-\left(m_x[f]\sin \vartheta -m_y[f]\cos \vartheta \right)\frac{\partial f}{\partial p}=0.$$

$$U[f](\vartheta )=A_0[f]+\sum _{n=1}^{\infty }{A}_n[f]\cos \left(n\vartheta \right)+\sum _{n=1}^{\infty }{B}_n[f]\sin \left(n\vartheta \right),$$

$$U[f](\vartheta )=A_1[f]\cos \vartheta +B_1[f]\sin \vartheta .$$

$$A_1[f]$$

$$B_1[f]$$

$$U[f](\vartheta )={\int }_{-\infty }^{+\infty }d{p}^{\prime }{\int }_{-\pi }^{\pi }d{\vartheta }^{\prime }f({\vartheta }^{\prime },p^{\prime })v(\vartheta -{\vartheta }^{\prime }),$$

$$U[f](\vartheta )=-m_x[f]\cos \vartheta -m_y[f]\sin \vartheta ,$$

$$f(\vartheta ,p,t)=f(-\vartheta ,-p,t)$$

$$\frac{\partial f}{\partial t}+p\frac{\partial f}{\partial \vartheta }-m[f]\sin \vartheta \frac{\partial f}{\partial p}=0.$$

$$\frac{\partial f}{\partial t}+p\frac{\partial f}{\partial \vartheta }-m[f]\left[\sum _{j=0}^J\frac{(-1{)}^j}{(2j+1)!}{\vartheta }^{2j+1}\right]\frac{\partial f}{\partial p}=0.$$

$$\frac{\partial f}{\partial t}+p\frac{\partial f}{\partial \vartheta }-m[f]\left[\sum _{j=0}^J\frac{(-1{)}^j}{(2j+1)!}{\vartheta }^{2j+1}\right]\frac{\partial f}{\partial p}+pf\left[\delta (\vartheta +\pi )-\delta (\vartheta -\pi )\right]=0.$$

$$I_{k,n}(t)=⟨{\vartheta }^k{p}^n⟩={\int }_{-\infty }^{+\infty }dp{\int }_{-\pi }^{\pi }d\vartheta f(\vartheta ,p,t){\vartheta }^k{p}^n.$$

$${\dot{I}}_{k,n}=k{I}_{k-1,n+1}-nm\sum _{j=0}^J\frac{(-1{)}^j}{(2j+1)!}{I}_{k+2j+1,n-1},$$

$$m=\sum _{j=0}^{J+1}\frac{(-1{)}^j}{(2j)!}{I}_{2j,0}.$$

$${\dot{I}}_{0,2}$$

$${\dot{I}}_{1,1}$$

$${\dot{I}}_{0,2}$$

$$m=1-\frac{1}{2}{I}_{2,0}.$$

$$\dot{x}$$

$$\dot{y}$$

$$\dot{z}$$

$$\dot{z}=-\left(1-\frac{1}{2}x\right)\dot{x}=-\frac{d}{dt}{\left(1-\frac{1}{2}x\right)}^2,$$

$$\epsilon =\frac{z}{2}-\frac{1}{2}{\left(1-\frac{1}{2}x\right)}^2=\frac{⟨p^2⟩}{2}-\frac{m^2}{2},$$

$$\dot{x}$$

Full text

Here we propose an approximation scheme to describe the collisionless relaxation process, based on the introduction of suitable moments of the distribution function, and apply it to a simple toy model, the Hamiltonian Mean Field (HMF) model. To the leading order, virial oscillations are equivalent to the motion of a particle in a one-dimensional potential. Inserting higher-order contributions in an effective way, inspired by the Caldeira-Leggett model of quantum dissipation, we derive a dissipative equation describing the damping of the oscillations, including a renormalization of the effective potential and yielding predictions for collective properties of the system after the damping in very good agreement with numerical simulations. Here we restrict ourselves to “cold” initial conditions, i.e., where the velocities of all the particles are set to zero: generic initial conditions will be considered in a forthcoming paper.

Long-range interactions; Vlasov equation; Hamiltonian Mean Field model; violent relaxation; cold collapse

I Introduction

Systems with long-range interactions, i.e., interactions that decay with the distance $r$ between the interacting bodies slower than $r^{-d}$, where $d$ is the dimension of space, exhibit peculiar features. A striking peculiarity of these systems is that relaxation to thermal equilibrium occurs on time scales $\tau_{\text{rel}}$ that grow with the number $N$ of constituents of the system, diverging in the $N\to\infty$ limit Campa et al. (2014, 2009). Thus, in practice, a many-body long-range-interacting system never reaches equilibrium and is typically found in a non-thermal state. The paradigmatic example of a system with long-range interactions that can be found in nature is a self-gravitating system where111For self-gravitating systems a proper thermal equilibrium state in the usual sense does not exist, at least in three dimensions, because gravity is non-confining so that a Maxwellian velocity distribution would lead to the evaporation of the system, unless the latter has infinite mass. The relaxation time here is the time scale over which binary encounters induce a loss of memory of the initial conditions and the Boltzmann entropy grows. Gravity in spaces with dimension less than three (where a thermal equilibrium state is well defined) also has relaxation times diverging with $N$. $\tau_{\text{rel}}\propto N/\ln N$: for a typical elliptical galaxy the relaxation time is of the order of $10^{17}$ years, i.e., several orders of magnitude larger than the age of the Universe Binney and Tremaine (2008). Other examples of long-range interactions are given by dipolar forces in three-dimensional condensed matter systems Campa et al. (2014), unscreened electromagnetic interactions in plasmas Nicholson (1983), and effective interactions mediated by the electromagnetic field in systems of cold atoms in an optical cavity Schütz and Morigi (2014); Gupta and Casetti (2016). Such a striking dynamical behaviour is due to the fact that in long-range-interacting systems mean-field interactions dominate over binary encounters (collisions), so that the time evolution of the one-particle distribution function $f(\mathbf{q},\mathbf{p},t)$ is given by the collisionless Boltzmann—also referred to as Vlasov—equation for times $t<\tau_{\text{rel}}$. Being such an equation invariant under time reversal, the Boltzmann entropy is a constant of motion and no relaxation towards thermal equilibrium can occur as long as $f$ obeys the Vlasov equation. We might therefore expect that, in the typical situation, the distribution function $f$ will not stop its evolution until the finite-$N$ collisional effects become important, bringing the system towards the thermal equilibrium state corresponding to the maximum of the Boltzmann entropy.

However, the typical evolution is not the naively expected one: when starting from generic initial conditions, the system undergoes collective oscillations and apparently settles down after much shorter time scales, not depending on its size, remaining trapped in a non-thermal state virtually forever, when $N$ is sufficiently large. This first part of the dynamical evolution has been called “violent relaxation” by Donald Lynden-Bell who first proposed a statistical-mechanical approach to this phenomenon Lynden-Bell (1967). The prototypical example is the collapse of an isolated cloud of self-gravitating particles in three dimensions, where starting from the pioneering numerical works by Hénon Hénon (1964) and van Albada van Albada (1982) (see Sylos Labini (2012) and references therein for more recent results) violent relaxation has been observed in the form of strong collective oscillations that damp out on a short time scale, of the order of a few dynamical times222The dynamical time $\tau_{D}$ of a self-gravitating system is defined as $\tau_{D}\approx\sqrt{\pi R^{3}/(GM)}$, where $R$ is the size if the system, or $\tau_{D}\approx\sqrt{1/(G\varrho)}$ where $\varrho$ is the average mass density, and is an estimate of the crossing time, i.e., the time needed for a typical particle to cross the entire system., leaving the system in a non-thermal but virialized stationary state. The collective oscillations during violent relaxation of a self-gravitating system are usually referred to as “virial oscillations”, since the virial ratio $2K/U$, where $K$ is the total kinetic energy and $U$ is the total gravitational potential energy, oscillates around, and eventually sets to, the stationary value $2K/U=-1$. In the last decades it has become clear that violent relaxation is not a peculiarity of three-dimensional self-gravitating systems: for instance, such a phenomenon has been observed in two-dimensional self-gravitating systems, when they collapse from a dilute state as well as when an initially thermal state is strongly perturbed Di Cintio et al. (2018), in mean-field toy models, again when an initially thermal state is strongly perturbed Teles et al. (2015), and in semiclassical models of cold atoms in an optical cavity, after a quench from a thermal equilibrium state Gupta and Casetti (2016). It is now clear that virial oscillations and violent relaxation are universal properties of long-range-interacting systems, occurring with essentially the same features as in the gravitational collapse case Campa et al. (2014). After the damping of the virial oscillations the system is found in a non-thermal state that is usually referred to as a quasi-stationary state (QSS), to indicate that such a state is not a true equilibrium of the system, being not stable under the action of collisions; in the $N\to\infty$ limit, however, quasi-stationary states become true stationary states of the Vlasov equation. Given that long-range-interacting systems spend virtually an infinite amount of time in quasi-stationary states, their prediction and characterization becomes maybe the most relevant task concerning this class of physical systems. However, we are still not able to predict the actual state in which the system is going to settle down after violent relaxation as, for example, the final state reached after a gravitational collapse, in the self-gravitating case. Indeed, to make a physically relevant example, elliptical galaxies are thought to be in a quasi-stationary state—also on purely observational grounds, the fact that most elliptical galaxies do share the same kind of luminosity profiles is an indication that some non-collisional relaxation process has occurred—but we are not able to predict such a state from first principles, given the initial conditions, and it is not even clear actually how pronounced the dependence on the initial conditions is Binney and Tremaine (2008); Bertin (2000).

Theoretical approaches to violent relaxation started with Lynden-Bell’s theory Lynden-Bell (1967), motivated by the attempt to predict the density profile of elliptical galaxies by means of a statistical-mechanical-like approach; however, it appears that Lynden-Bell’s theory only works well when virial oscillations are somehow suppressed Levin et al. (2014); Sylos Labini (2012). Despite many advances since then (see e.g. Campa et al. (2014); Levin et al. (2014) for reviews and Leoncini et al. (2009); Campa and Chavanis (2010); Chavanis and Campa (2010); Assllani et al. (2012); Benetti et al. (2014) for applications to particular systems), a satisfactory theory of violent relaxation is still lacking for any long-range-interacting systems, even the simpler ones, not only for galaxies. The phenomenology of virial oscillations and violent relaxation suggests that the mechanism behind such a collisionless relaxation is similar to Landau damping in plasmas, where the dynamics is reversible but the collective oscillations lose their energy towards the motions of the particles via resonant interaction (see e.g. the discussion in Kandrup (1998)). However, this kind of Landau damping would be fully nonlinear and would occur in non-homogeneous states; at variance with the linear and homogeneous case that admits an exact solution in terms of uncoupled normal modes, the nonlinear and non-homogeneous case presents technical challenges that have not been overcome yet, although progress has been made in particular cases Barré et al. (2010, 2011). It is worth mentioning that a phase-space interpretation of collisionless relaxation is possible (see e.g. Tremaine et al. (1986), or Mohout and Villani (2011) for a more mathematically-oriented discussion), and would be useful in the following. As time progresses, the Vlasov evolution splits the distribution function $f$ in finer and finer filaments that fold on themselves. As a consequence the dynamics relocates on smaller and smaller scales in phase space: then, rigorously speaking, the evolution in phase space does not admit any asymptotic state. However, any coarse-grained distribution function $\tilde{f}$ will not be affected by this fine dynamics and will eventually stop its evolution, so that the coarse graining actually introduces an effective time arrow in the system. Despite the dynamics would never really stop, any physically meaningful observable related to not too small scales will settle to a pseudo-equilibrium value.

In the present paper we aim at giving a contribution to the understanding of this problem. In order to reduce the complexity, we shall restrict ourselves to a very simplified toy model, the Hamiltonian Mean Field (HMF) model (see Sec. II), that, despite being one of the simplest models of a long-range-interacting systems, does exhibit all the above mentioned phenomenology. Moreover, we shall restrict ourselves to “cold” initial conditions, i.e., where the velocities of all the particles are set to zero. We shall introduce a hierarchy of moments of the distribution function, in order to study the dynamics of the lowest-order moments by taking into account the contribution of the higher-order ones in an effective way, in close analogy to the Caldeira-Leggett model of quantum dissipation Caldeira and Leggett (1981, 1983): as we shall see, this will allow us to write an effective equation of motion for the low-order moments where the contribution of the other moments results in an explicit dissipation as well as a renormalization of the low-order dynamics. In a forthcoming paper Giachetti et al. (2019) we shall generalize the approach introduced here to generic initial conditions.

The paper is organized as follows: in Sec. II we introduce the HMF model, in Sec. III we introduce our theoretical approach, apply it to the cold collapse of the HMF model and compare its predictions with numerical simulations, and in Sec. IV we draw our conclusions and discuss open points and future developments. Two appendices are devoted to technical and/or side aspects.

II Hamiltonian Mean Field model

The Hamiltonian Mean Field (HMF) model is a toy model that has become paradigmatic for the study of equilibrium as well as nonequilibrium properties of systems with long-range interactions. According to Chavanis and Campa Chavanis and Campa (2010), it was introduced by Messer and Spohn Messer and Spohn (1982)—who called it the “cosine model”—after a suggestion by Battle Battle (1977); the model was then made popular by Antoni and Ruffo, who also introduced the name and the acronym HMF, in a seminal paper Antoni and Ruffo (1995). Many works have been devoted to the study of the HMF model since then: a review of the main results can be found in Campa et al. (2014, 2009). The HMF model is a classical system with $N$ degrees of freedom, defined by the Hamiltonian333Note that our Hamiltonian (1) differs—by an additive constant—from the one used in many works on the HMF model where the interaction energy is chosen such as to be always larger than zero.

[TABLE]

where $\vartheta_{i}\in[-\pi,\pi]$ are angular coordinates, $p_{i}$ are their conjugated momenta and we have chosen natural units such that the coupling constant and the inertia of each degree of freedom are equal to unity. We note that we chose the coupling constant to be positive, so that the interaction is attractive (but a repulsive version of the HMF has been considered too, see e.g. Antoni and Ruffo (1995); Campa et al. (2014)) and the Kac scaling has been used, such that the Hamiltonian (1) is extensive, although clearly non additive.

The HMF model admits (at least) two different interpretations. It can be seen either as a fully connected system of planar ($XY$) classical spins, each with a finite inertia and parametrized by one of the angles $\vartheta_{i}$, interacting via a classical, ferromagnetic Heisenberg exchange interaction444The acronym HMF can also stand for Heisenberg Mean Field., or as a system of classical particles moving on a circle, whose center coincides with the origin of the reference frame and where the angles $\vartheta_{i}$ are the angular polar coordinates of each particle, interacting via the potential

[TABLE]

Although the particle interpretation is the most natural in our context, the terminology commonly used for the HMF model stems from the magnetic interpretation. Introducing the magnetization vector $\mathbf{m}=\left(m_{x},m_{y}\right)$ such that

[TABLE]

the interaction energy (2) becomes

[TABLE]

where $m=\left|\mathbf{m}\right|=\sqrt{m_{x}^{2}+m_{y}^{2}}$. The vector $\mathbf{m}$ is indeed the (instantaneous) magnetization of the system (1); in the particle interpretation, it is a vector pointing to the direction of the maximum particle density on the circle and whose modulus is related to the degree of clustering of the particles: $m=0$ corresponds to a uniform density and $m=1$ to a completely collapsed (delta-like) distribution of particles555The square modulus $m^{2}$ of the magnetization plays an analogous rôle to the reciprocal of the gravitational radius in a self-gravitating system.. The Hamilton equations of motion derived from (1) can be written as

[TABLE]

where it is apparent that the coupling between the degrees of freedom comes only from the magnetization, whence the “mean field” in the name of the system.

In the limit $N\to\infty$ we can introduce the single-particle distribution function $f(\vartheta,p,t)$ and replace the sums $\frac{1}{N}\sum_{i=1}^{N}$ with the phase space averages weighted by $f$, i.e., the integrals $\int_{-\infty}^{+\infty}dp\int_{-\pi}^{\pi}d\vartheta\,f$. The magnetization components become functionals of $f$,

[TABLE]

the mean-field potential is given by

[TABLE]

so that the distribution function evolves in time according to the Vlasov equation

[TABLE]

It is worth noticing that the mean-field interaction of the HMF model can be seen as the lowest-order expansion of a generic mean-field interaction for particles on a circle. Indeed, on a circle any interaction must be periodic, so that we can expand it in a Fourier series,

[TABLE]

and neglecting the constant term $A_{0}$ by redefining the zero level of the interaction and keeping only the lowest order in the expansion, we get

[TABLE]

Requiring rotational invariance on the circle, the only possible preferred direction is that specified by the mean position of the particles, implying

[TABLE]

Being the mean-field potential defined as

[TABLE]

it is linear in $f$, so that the quantity $\mu$ in Eqs. (11) cannot depend of $f$ and is therefore a scalar constant that can be absorbed in the choice of units and set equal to $-1$, since it has to be negative if we want the potential to be attractive. We finally get

[TABLE]

i.e., the mean-field potential of the HMF model, that can thus be regarded as a set of particles interacting with the lowest Fourier modes of a generic attractive collective interaction on a circle. Indeed, it has been shown that the HMF interaction stems from a softened gravitational interaction between particles constrained on a circle, in the limit of an infinitely large softening length Tatekawa et al. (2005).

From now on we shall only consider initial conditions that are symmetric around $\vartheta=0$ and with a total momentum equal to zero, i.e., such that $f(\vartheta,p,0)=f(-\vartheta,-p,0)$. Being such invariance conserved by the equations of motion, in their discrete version (5) as well as in the continuum Vlasov limit (8), the distribution function will be such that

[TABLE]

at any time $t$; this implies not only a vanishing total momentum but also $m_{y}\equiv 0$, so that the magnetization will be always along the $\vartheta=0$ axis and $m\equiv m_{x}$. The Vlasov equation then becomes

[TABLE]

The periodic boundary conditions imply that $f(\pi,p)=f(-\pi,p)$, so that the single-particle phase space has the topology of a cylinder; this, together with the condition (14), means that we can also see it as the sheet $\vartheta\in[0,\pi]$ in the $(\vartheta,p)$ Euclidean plane, with two perfectly elastic walls in $\vartheta=0$ and $\vartheta=\pi$.

As a final remark for this Section, since we are mainly concerned with the HMF model as a toy model of a self-gravitating system, it is worth mentioning that the HMF model has been considered in an astrophysical context in Pichon (1994) and, more recently, in Fouvry and Bar-Or (2018).

III Cold collapse in the HMF model

Let us now consider a cold initial distribution function, i.e., such that all the velocities are equal to zero, peaked (even very weakly) around $\vartheta=0$ (so that $m(0)=m_{0}>0$, with $m_{0}$ possibly very small). As one can easily guess, the attractive nature of the potential will immediately lead to a collapse that, in turn, results in the growth of $m$ and then in the loss of potential energy that is turned into kinetic energy. After reaching the minimum of the collective potential, particles will start to turn kinetic energy back into potential energy as the system bounces back to a less clustered state. We thus expect a behaviour that reminds of a Jeans-like instability, where an initial collapse is followed by collective oscillations, analogous to the virial oscillations observed in a self-gravitating system. Indeed, we may expect the dynamical evolution of an HMF model after a cold initial condition to be very similar to a gravitational collapse, and numerical simulations show that this is true. Moreover, simulations show that, as in the gravitational case, after a relatively short time violent relaxation sets in, the virial oscillations damp out and the system appears to relax to a quasi-stationary state, where collective quantities as $m$ are nearly constant. An example of the behaviour of the magnetization $m(t)$ of an HMF model with $N=2\times 10^{6}$ particles in a cold collapse starting with very small magnetization ($m_{0}\approx 10^{-2}$) is shown in Fig. 1. The simulation whose results are shown in Fig. 1 as well as all the other simulations of the HMF model whose results are shown in the following have been conducted by integrating Eqs. (5) by means of a third-order bilateral symplectic algorithm Casetti (1995).

III.1 Moments of inertia of the distribution function

Let us now introduce our theoretical approach to virial oscillations and collisionless relaxation in a cold collapse of the HMF model. First of all, let us observe that all the states of the system are expected to become rather tightly collapsed, even after starting with small magnetization, as shown in Fig. 1. We can thus replace the $\sin\vartheta$ term in the Vlasov equation (15) with its Taylor expansion up to a finite order $2J+1$, so that the time evolution of $f$ is given by

[TABLE]

For finite and not very large values of $J$ this approximation will work better when the state is tightly collapsed, as noted before: in spite of that, for nearly homogeneous states $m\approx 0$ and the depth of the potential well $-m\cos\vartheta$ is small, so that the actual form of the potential is quite irrelevant (as long as the minimum remains in $\vartheta=0$).

The truncated potential from which the polynomial force given by the finite expansion in (16) arises is somewhat analogous to the Villain approximation to the XY model Villain (1975), although the latter is usually truncated at the harmonic order, while we consider a generic finite order. However, such a force is discontinuous in $\vartheta=\pm\pi$ on the circle; to avoid this complication, that would make the distribution function nondifferentiable, we “open” the ring and insert elastic walls in $\vartheta=\pm\pi$, replacing Eq. (16) with

[TABLE]

We now define the generalized moments of inertia of the distribution function $f$ as follows

[TABLE]

Higher order moments probe the finer grain of the distribution, so that we expect the lower order moments, blind to the fine structure, to settle before the others. To obtain an equation of motion for the moments $I_{k,n}$ we take the time derivative of the right-hand-side of Eq.  (18); by exchanging integration and differentiation, and using Eq. (17) to express $\partial_{t}f$, after integrations by parts and assuming that $f$ decays sufficiently fast for large $p$’s we get rid of the boundary terms and obtain a system of equations that contains only the moments,

[TABLE]

where $m$ can be expressed in terms of the moments as the following series up to order $2J+2$,

[TABLE]

The quantity $m$ introduced in Eq. (20) is the truncated magnetization, and the truncation is made at an order consistent to that of the finite polynomial expansion of the $\sin\vartheta$ term in the Vlasov equation (16). Hence, Eqs. (19) and (20) yield a closed system of equations at any given order $J$: $J=0$ involves only moments with $k+n\leq 2$, $J=1$ moments with $k+n\leq 4$, and so on. We note that $I_{0,0}\equiv 1$ and that, due to the fact that $f$ is even in $\vartheta$ and $p$, $I_{k,n}\not=0$ only if $k+n$ is even. Equations (19) and (20) clearly show the hierarchy in the interactions: lower order moments strongly interact with each other while are only weakly affected by the higher order ones, which in turn are strongly forced by the lower order moments. In particular, the dynamics of the low-order moments will describe the virial oscillations, and we expect damping to be a result of the back-action of the higher-order inertia moments, that interfere incoherently with the virial oscillations as depicted in the scheme drawn in Fig. 2.

In the following we shall elaborate on this idea.

III.2 Leading order description of virial oscillations

The leading order ($J=0$, $k+n\leq 2$) in our approximation scheme involves a self-consistent harmonic approximation for the mean-field potential. The only moments to be considered are $I_{2,0}=\langle\vartheta^{2}\rangle$, related to the width of the spatial distribution of the particles, $I_{0,2}=\langle p^{2}\rangle$, proportional to the average kinetic energy, and $I_{1,1}=\langle\vartheta p\rangle$, i.e., the covariance between space and velocity variables. Their time evolution is given by Eqs. (19) and (20) that reduce, to the leading order, to

[TABLE]

and

[TABLE]

By setting, to ease the notation, $x=I_{0,2}$, $y=I_{1,1}$ and $z=I_{0,2}$, Eqs. (21) and (22) can be rewritten as

[TABLE]

From the above equations we get

[TABLE]

so that Eqs. (23) admit the integral of motion

[TABLE]

that is, the energy per particle: as it should be, $\varepsilon\geq-1/2$. Inserting Eq. (25) into Eq. (23b) we get

[TABLE]

that is, Hamilton equations of motion of a particle of mass $1/2$ moving in an effective potential $V_{\text{eff}}(x)$; indeed Eqs. (26) can be cast as

[TABLE]

where

[TABLE]

Clearly also the energy of this motion,

[TABLE]

is a constant. We note that

[TABLE]

so that $\Lambda=y^{2}-xz$; since the definition of the $I_{k,n}$ implies $y^{2}\leq xz$, one has

[TABLE]

The case of completely cold initial conditions, i.e., of vanishing initial kinetic energy, translates into the condition $z_{0}=z(0)=0$. Let us then considering the initial conditions on $y$ and $x$. As far as $x_{0}=x(0)$ is concerned, we have

[TABLE]

where $m_{0}$ is the initial magnetization. Since $y^{2}\leq xz$, Eq. (32) implies $y_{0}=y(0)=0$, while Eq.  (31) implies $\Lambda=0$, so that a cold collapse corresponds to a motion with zero total energy and vanishing initial velocity (thus, with zero initial potential energy) in the effective potential (28). The total energy is completely specified by the initial magnetization according to

[TABLE]

so that $\varepsilon\in\left[-\frac{1}{2},0\right]$ and $V_{\text{eff}}$ can be expressed as

[TABLE]

and is plotted as a function of $x$ for some values of $m_{0}$ in Fig. 3. The particle starts at $x_{0}$ with zero velocity and, since $V_{\text{eff}}(x_{0})=0$, oscillates in the allowed region given by $V_{\text{eff}}\leq 0$, that is

[TABLE]

In the oscillations the particle always reaches $x=0$, that is $m=1$. This is coherent with the expectation, and the numerical evidence, that in a cold collapse the system performs large-amplitude oscillations and always comes close to a maximally collapsed configuration. Initially larger magnetizations (corresponding to smaller $x_{0}$) give rise to smaller-amplitude oscillations, being the latter always bounded between [math] and $x_{0}$, again a feature that is observed in simulations.

When $m_{0}=1$ the potential has a minimum in $x_{0}=0$ and no evolution is possible, since the allowed region for $x$ shrinks to the point $x_{0}$ itself: this makes sense because the whole system is in the bottom of the $N$-body potential well. Also the initial position $x_{0}=2$, corresponding to $m_{0}=0$, is a stationary state of the potential, but it is a maximum, so that any $m_{0}$ slightly larger than zero would lead to a collapse. As expected, this is an unstable equilibrium, that corresponds to the instability of a cold uniform distribution. Hence, the leading-order description we have considered so far captures many properties of the initially cold virial oscillations. However, the oscillations predicted by Eqs. (26) would never damp: damping is not included in this description666This feature is shared by other effective descriptions of virial oscillations like the so-called “envelope equations” derived by Levin et al. Levin et al. (2014) to model the virial oscillations after water-bag initial conditions.. This notwithstanding, one may assume that the damping, due to the back-action of the higher-order moments neglected so far, can be modeled by just adding the simplest dissipative term to Eqs. (26), that is, by adding a term proportional to $-\dot{x}$ to Eq. (27). Although, as we shall see in Sec. III.3, such an effective dissipation should be accompanied by a renormalization of the effective potential, this crude approximation indeed yields reasonable results for the asymptotic magnetization after the damping of the oscillations, that would simply coincide with the value $\overline{m}$ of $m$ corresponding to the minimum of $V_{\text{eff}}$. The first derivative of $V_{\text{eff}}(x)$ has a root in

[TABLE]

so that, using $x=2\left(1-m\right)$, we get

[TABLE]

This result implies that a cold collapse will always lead to a state with magnetization larger than a finite value that is predicted as $\overline{m}_{\text{min}}=\frac{2}{3}$; the latter lower limit corresponds to the initial condition $m_{0}=0$. Moreover, $\overline{m}$ increases with $m_{0}$, eventually reaching the already discussed fixed point $\overline{m}=m_{0}=1$. All these properties are observed in the numerical simulations. The estimate for the asymptotic magnetization given by Eq. (37), that contains no adjustable parameters, is compared with a more refined prediction to be derived in Sec. III.3 and with the results obtained in numerical simulations in Fig. 5. It is apparent that Eq. (37) systematically overestimates all the numerical results, including $\overline{m}_{\text{min}}$, but the trend is the correct one.

As a further consistency check of what we have derived so far we note that the leading-order prediction (37) for the magnetization after the damping of the oscillations is not only consistent with the virial theorem, as it should be for any stationary state, but it could have been directly derived from the virial theorem itself, exploiting the fact that the leading order corresponds to a harmonic approximation for the interparticle potential (see appendix A for details).

III.3 Effective dissipation

Let us now discuss how to effectively tackle the contribution of all the higher-order moments we have neglected so far and to show that such a contribution indeed yields an effective damping as well as a renormalization of the effective potential (34).

We start by observing that, were the magnetization a given constant $m=\mu$ instead of a self-consistent function of time, and in the harmonic approximation that amounts to set $J=0$ in Eqs. (19), the equations of motion for the moments would become a set of independent linear systems of differential equations; each system would be composed of $2L+1$ coupled equations, where $2L=k+n$ is a non-negative even number, given by

[TABLE]

The real part of the eigenvalues of the matrix777Being the interaction among the $I$’s in Eqs. (38) only between nearest neighbours, these matrices are tridiagonal and thus diagonalizable. associated to each system vanishes. In fact Eqs. (38) stem from the Vlasov equation for the HMF model (15) in the harmonic approximation (i.e., with $\sin\vartheta\approx\vartheta$) and with $m=\mu$. In this case the evolution in phase space reduces to a rigid rotation, i.e., to a periodic motion. Also the dynamics of the $I$’s must then be periodic and this rules out the possibility of a non-imaginary eigenvalue. Moreover, being the equations real, the conjugate of any eigenvalue is still an eigenvalue and since their number is odd there must be a zero mode. This means that the corresponding left eigenvector does not evolve and we have a first integral, linear in the $I_{k,n}$’s. Summarizing, we have a $2L$-dimensional effective linear dynamics whose eigenvalues are couples of imaginary conjugate numbers: we can thus consider any of the systems of $2L+1$ equations with $k+n=L$ in Eqs. (38) as the Hamilton equations of $2L$ decoupled harmonic oscillators. In addition, the periodicity implies that all the frequencies of the oscillators are integer multiples of a common fundamental one. When we take into account the dependence of $m$ on the moments $I_{k,0}$ the equations become nonlinear and this simple picture is no longer valid. Nonetheless, being the equations invariant under time reversal, it is reasonable to assume the $L>1$ moments to act as coupled mechanical systems exhibiting an oscillatory behaviour in the proximity of equilibrium. However, since the periodicity is lost, no assumption can be made on the values of the actual frequencies.

In close analogy with the works by Caldeira and Leggett Caldeira and Leggett (1981, 1983) aiming at modeling a quantum dissipative environment in order to describe quantum brownian motion, we assume that the action of the high-order moments can be considered as that of a bath of independent harmonic oscillators coupled to the low-order ones. In this classical counterpart of the Caldeira-Leggett mechanism, dissipation emerges as the result of the loss of coherence of the set of oscillators, whose dynamics is individually reversible, when the number of oscillators diverges. The effective dissipation mechanism is thus seen as a “classical decoherence”. It is worth noting that this feature is not exclusive of the Caldeira-Leggett model: at the classical level also other approaches exist where an effective dissipation may emerge from the coupling with an infinite assembly of oscillators (see e.g. Zwanzig (2001)).

Let us now go into the details. We retain the same notation as before, i.e., $x=I_{0,2}$, $y=I_{1,1}$ and $z=I_{0,2}$, and we set

[TABLE]

where $\xi$ sums up all the contributions of the higher-order moments and is implicitly assumed small. The equations of motion (23) for $x$, $y$ and $z$ would then become

[TABLE]

but being $\xi$ an (unknown) function of time we cannot obtain an integral of motion by manipulating these equations as above. Still, we impose the conservation of energy

[TABLE]

so that we can express $z$ as a function of $\varepsilon$, $x$ and $\xi$ and eliminate it from the equations (40), that become

[TABLE]

Introducing

[TABLE]

such that

[TABLE]

and neglecting $\mathcal{O}(\xi^{2})$ terms, Eqs. (42) are equivalent to

[TABLE]

where $V_{\text{eff}}(x)$ is still given by Eq. (28). We have now to describe the dynamics of $\xi$. We assume that $\xi$ is a sum of many oscillating terms,

[TABLE]

where the $q_{k}$’s are the coordinates of a set of $M$ harmonic oscillators, with $M$ large; $\zeta_{0}$ and $c_{1},\ldots,c_{M}$ are real constants. The first term in Eq. (46) implies that the coupling between low- and high-order modes may modulate the average values of the oscillations. As we shall see in the following, we can think of such a term as a counter-term that tunes the renormalization of the effective potential, allowing to make a simplifying assumption on the spectrum of the oscillators. Equation (45) then becomes

[TABLE]

To derive the equations of motion for the oscillators, we assume the above equation derives from a Lagrangian $L$ given by

[TABLE]

where

[TABLE]

describes the $x$ sector, i.e., the dynamics of the lowest-order modes,

[TABLE]

describes the bath of harmonic oscillators modeling the higher-order moments,

[TABLE]

describes the interaction between the bath and the lowest-order moments888Such an interaction is linear in the $q_{k}$’s as in the classical Caldeira-Leggett model, but at variance with the latter model the coupling with the $x$ degree of freedom is nonlinear and realized via the function $g(x)$., and

[TABLE]

corresponds to the counter-term. From the Lagrangian (48) not only Eq. (47) follows, but also the equations of motion of the bath degrees of freedom, as

[TABLE]

with $k=1,\ldots,M$. From Eq. (53) we see that the low-order modes act on the oscillators of the bath by displacing their equilibrium position. Let us now solve the dynamics of the oscillators of the bath: this can be done thanks to the linearity in $q_{k}$ of Eqs. (53). As shown in Appendix B, we get, for $k=1,\ldots,M$,

[TABLE]

where the first term on the r.h.s. is the homogeneous solution

[TABLE]

while the other two terms come from the interaction. Using Eq. (54) we can write

[TABLE]

We now want to take the limit of a large number of oscillators, $M\to\infty$. In order to replace the sums with integrals over the frequencies we introduce the spectral density $J(\omega)$ defined as

[TABLE]

so that Eq. (56) becomes

[TABLE]

where

[TABLE]

To go further we need some assumptions on the spectral density. As in the standard Caldeira-Leggett model we assume a white spectrum up to a cutoff $\Omega$, that is, an Ohmic spectral function,

[TABLE]

and substituting the latter into Eq. (58) we get

[TABLE]

The second term in the r.h.s. of the above equation, that is, the displacement of the oscillators due to the coupling, is proportional to $\Omega$, so that it would diverge in the $\Omega\to\infty$ limit. However, the quantity given by Eq. (61) is not an observable quantity, at variance with $\xi$; inserting Eq. (61) into Eq. (46) we get

[TABLE]

so that we can reabsorb the divergence into the bare coupling constant $\zeta_{0}$, and by taking the limit $\Omega\to\infty$ we get

[TABLE]

where $\zeta$ is the renormalized (and thus unknown) coupling

[TABLE]

and where we have used the integral representation of the Dirac $\delta$ distribution. By interchanging the derivative and the integral in Eq. (63) we obtain

[TABLE]

so that, after substituting Eq.  (65) into (45), the equation of motion for $x$ finally reads

[TABLE]

Being $\eta\left[g^{\prime}(x)\right]^{2}\geq 0$, the third term in the r.h.s. of Eq. (66) is indeed a damping term, although depending also on $x$. Given that $\xi_{0}(t)$ is an infinite sum of oscillating terms with zero mean and a white spectrum, each term being given by Eq. (55), it is equivalent to a white noise with zero mean. Hence Eq. (66) is a Langevin equation with multiplicative noise, at variance with the classical Caldeira-Leggett model where the noise is purely additive. The r.m.s. amplitude of $\xi^{0}(t)$ depends on the initial conditions of the oscillators; were the latter in thermal equilibrium at a given temperature $T$, as in the classical Caldeira-Leggett model, one would have $\langle\xi^{0}(t)\xi^{0}(t^{\prime})\rangle\propto\eta T\delta(t-t^{\prime})$, in agreement with the fluctuation-dissipation theorem Caldeira and Leggett (1981, 1983). However, in our case there is no thermal bath, and the oscillators only model the contribution of the higher-order moments, that is assumed to be small for our theory to be consistent. We thus make the further assumption that the initial conditions of the oscillators are such that the amplitude of the noise (or the effective temperature of the bath) is very small, so that we can neglect the multiplicative noise in the following.

The multiplicative nature of the noise and the dependence on $x$ of the damping term are not the only differences between Eq.  (66) and the Langevin equation of the classical Caldeira-Leggett model: there is also the term involving the coupling $\zeta$, that can be seen as a “dissipative correction” that renormalizes the effective potential $V_{\text{eff}}$. After neglecting the noise term and recalling that $g(x)=x(2-x)$ we can indeed write Eq. (66) as

[TABLE]

where

[TABLE]

is the renormalized effective potential. As observed above, the value of the parameter $\zeta$, as well as that of the friction coefficient $\eta$, can not be predicted by the theory. Yet, they both have to be small for the theory to be consistent. Moreover, Eq. (68) implies that the dissipative correction to the effective potential vanishes when $x=0$ and $x=2$, regardless of the value of $\zeta$, and this is physically sound, because these values of $x$ correspond to $m=1$ and $m=0$, respectively, and are fixed points of the dynamics of the system, that are thus left unchanged by the dissipative correction. A comparison between the bare and renormalized effective potentials is shown in Fig. 4.

Equation (67) implies that the asymptotic values of $x$, that determine those of the magnetization $m$, are the minima $\overline{x}^{(R)}$ of $V^{(R)}_{\text{eff}}(x)$. To calculate them we exploit the fact that the dissipative correction is small, writing $\overline{x}^{(R)}=\overline{x}+\delta\overline{x}$ where $\delta\overline{x}=\mathcal{O}(\zeta)$ and expanding $V^{\prime(R)}_{\text{eff}}(x)$ around $\overline{x}$ up to first order in $\zeta$, obtaining

[TABLE]

Equation (69) implies that the dissipative correction $\delta\overline{m}$ to the leading-order estimate of the magnetization $\overline{m}=1-\overline{x}/2$ is

[TABLE]

so that if $\zeta>0$ we have a negative correction, that is the correct one since, as already observed, $\overline{m}$ overestimates the numerical results. Using the expression of $\overline{x}$ given by Eq. (36), the the leading-order prediction of the asymptotic magnetization $\overline{m}^{(R)}$, including the dissipative correction, as a function of the initial magnetization $m_{0}$ is thus

[TABLE]

where $\overline{m}$ is given by Eq. (37) and

[TABLE]

with

[TABLE]

The leading-order prediction of the asymptotic magnetization $\overline{m}$, given by Eqs. (71) and (72) is plotted in Fig. 5 against the initial magnetization $m_{0}$ and compared with the leading-order prediction without the dissipative correction given by Eq. (37) and the results of numerical simulations of cold collapse of an HMF model with $N=2\times 10^{5}$ particles, for two different classes of initial conditions. Choosing $\zeta=0.11$ the agreement between the theoretical prediction and the simulated data is very good, although the leading-order prediction depends only on $m_{0}$ so that it cannot resolve the fine differences in the asymptotic magnetization obtained with different initial conditions corresponding to the same $m_{0}$.

Going beyond the prediction of the asymptotic value of the magnetization, we do not expect the solutions of Eq. (67) to be able to accurately reconstruct the actual time evolution $m(t)$ during the cold collapse; although Eq. (67) contains the damping, it still has only the information on the lowest-order moments. Indeed, when $m_{0}$ is small the agreement between theoretical and numerical outcomes of $m(t)$ is only qualitative, but when $m_{0}$ is large enough, so that the system is almost always in a tightly collapsed state, the Fourier power spectra of the solutions of Eq. (67) are not very different from the ones of the $m(t)$ obtained in numerical simulations: an example is shown in Fig. 6. This is a further confirmation of the fact that our effective description captures a non-negligible part of the actual dynamics already at the leading order.

It is worth mentioning that the relation between dissipation in Vlasov dynamics and the Caldeira-Leggett model had been previously discussed by Hagstrom and Morrison Hagstrom and Morrison (2011) but only in the special context of linear Landau damping in a homogeneous background. More precisely, they found an explicit mapping between the Vlasov-Poisson system and the Hamilton equations of the Caldeira-Leggett model, by expressing both systems in terms of normal modes. Although less rigorous than Hagstrom and Morrison’s, our results suggest that the analogy between collisionless dissipation and the Caldeira-Leggett mechanism can be pushed considerably forward to include the fully nonlinear case of violent relaxation towards non-homogeneous states.

III.4 Next-to-leading order description of cold collapse

The time evolution of the inertia moments of the distribution function can be studied, in principle, up to any desired order. Going to the next-to-leading order in the dynamics of the $I_{k,n}$ means considering $J=1$ and $k+n\leq 4$. This involves eight moments, including the three already appearing at the leading order, and the equations of motion of these moments, neglecting all the higher-order ones, read as

[TABLE]

where

[TABLE]

Equations (74) are invariant under time reversal, so that they cannot exhibit any damping and relaxation towards a fixed point. However, as in the case of the leading order, we expect the higher-order terms here neglected to provide an effective dissipation that we can model, in the simplest way, by adding dissipative terms to Eqs. (74) as follows:

[TABLE]

where $m$ is still given by Eq. (75) and the real positive constants $\gamma_{(L)}$, with $L=2$ and $L=4$, respectively, are the effective friction coefficients that are expected to be different for moments of different order: in particular, since the dissipation time scale of the higher-order moments should be longer than that of the lower-order moments, $\gamma_{(2)}>\gamma_{(4)}$. In Eqs. (76) we could have considered also off-diagonal dissipative terms, but since this is only an effective description and we do not expect the precise form of the dissipative terms to dramatically affect the main features of the dynamics, we chose the simplest possible form. We note that even with the inclusion of the friction terms the energy per particle $\varepsilon$ is conserved, because

[TABLE]

In principle we expect that a proper effective treatment of the higher-order moments would give rise not only to dissipative terms but also to a renormalization of the dynamics, like the dissipative correction to the effective potential arising at the leading order: however, since an analytical approach like that carried out in Sec. III.3 appears extremely difficult if not unfeasible at the next-to-leading order, we neglected this aspect and numerically solved Eqs. (76) for different initial conditions, estimating the asymptotic values of $m$ from the numerical solution. Therefore we cannot expect the prediction of Eqs. (76) to accurately match numerical data: yet, if our approach is consistent, such a prediction should be closer to the numerical data than the bare leading-order prediction given by Eq. (37). Moreover, while the asymptotic magnetization predicted at the leading order depends only on the initial magnetization $m_{0}$, here a dependence on finer details of the initial conditions may show up.

The space of the initial conditions of the system of equations (76) is eight-dimensional. However, considering only cold initial conditions implies that all the moments containing $p$ initially vanish, that is, the only moments that can be nonzero in the initial state are $I_{2,0}$ and $I_{4,0}$. The initial condition can thus be parametrized by two independent quantities. Observing that $\left\langle\left(\vartheta^{2}-\langle\vartheta\rangle^{2}\right)^{2}\right\rangle\geq 0$ implies $I_{4,0}\geq I^{2}_{2,0}$ we can define a parameter $b\geq 1$ such as

[TABLE]

and parametrize the initial condition by means of $m_{0}$ and $b$. Since

[TABLE]

given $m_{0}$ and $b$ the value of $I_{2,0}(0)$ is

[TABLE]

According to Eq. (80), if $b<3/2$ all the values of $m_{0}$ are allowed; if $b>3/2$ some values of $m_{0}$ are forbidden. For a uniform distribution $b=1.8$, and for the initial conditions of the simulations considered in Fig. 5 (those obtained with a Gaussian overdensity superimposed on a uniform distribution) $b\in(1.8,9)$ as long as $m_{0}\lesssim 0.8$ and becomes much larger for larger initial magnetizations.

Numerical results confirm that the next-to-leading order predictions of the asymptotic magnetization do not match the numerical data, yet systematically improve the bare leading-order $\overline{m}$ given by Eq.  (37). Moreover, a dependence on $b$ shows up, although rather weak: this is consistent with the fact that in a cold collapse the dominant effect is that of the lowest-order moments, so that the outcome is mostly determined by $m_{0}$. in Fig. 7 we plot the values of $\overline{m}_{\text{NLO}}$ obtained by means of the numerical solution of Eqs. (76) with $\gamma_{(2)}=0.1$ and $\gamma_{(4)}=0.05$ for some values of $b$, i.e., $b\in[1,4]$.

In order to appreciate the quantitative improvement with respect to the bare leading-order result, in Fig. 8 we show the difference $\overline{m}_{\text{NLO}}-\overline{m}$ between the next-to-leading order and the bare leading-order prediction for the appropriate values of $b$ relative to the initial conditions of a subset of the simulations results already shown in Fig. 5, together with the difference between the outcome of the simulations and the bare leading order value. It is apparent that for $m_{0}\gtrsim 0.4$ the next-to-leading order prediction substantially improves the base leading-order prediction, giving more than half of the needed correction to match the simulation data; for smaller initial magnetizations the improvement is less substantial.

We note that for small values of $b$ and small values of $m_{0}$ we obtain a uniform asymptotic state, i.e., $m=0$ (data not shown in Figs. 7 and 8). Given that when $\varepsilon<0$ the magnetization can not vanish, this is a shortcoming of our theoretical approach999Or better of the numerical implementation we are considering now at the next-to-leding order. that, as discussed above, is better suited for describing collapsed states, so that it may not work well for very small initial magnetizations. However, this concerns only a very small subset of the space of initial conditions.

Finally, we expect that going to the next-to-leading order should improve also the description of the actual time evolution $m(t)$, because finer details on the dynamics are now taken into account. Indeed, Fig. 9 shows that for the same value of $m_{0}$ considered in Fig. 6 the Fourier power spectrum of $m(t)$ obtained by means of a numerical simulation of the HMF model is better reproduced by the solution of Eqs. (76) than by the leading-order solution.

IV Concluding remarks

We have shown that an effective theory of collisionless (violent) relaxation after a cold collapse is possible, at least in the case of a very simple toy model with long-range interactions, the HMF model. Exploiting the fact that during collisionless relaxation the dynamics relocates of finer and finer scales, a description of the evolution of the collective properties of the system, blind to the fine structure, can be done by looking at the lowest-order moments of the distribution function. To the leading order we have shown that a description of the virial oscillations (without damping) naturally emerges as equivalent to the motion of a fictive particle in a one-dimensional potential. Then, inserting the contribution of the higher-order moments in an effective way, we have explicitly derived a dissipative equation describing virial oscillations and their damping, including a renormalization of the effective potential and yielding predictions for the value of the magnetization $m$ after the damping in very good agreement with numerical simulations. The derivation of the effective dissipation is inspired by Caldeira and Leggett’s treatment of open quantum systems Caldeira and Leggett (1981, 1983). As already mentioned in Sec. III.3, Hagstrom and Morrison Hagstrom and Morrison (2011) already discussed a relation between collisionless relaxation and the Caldeira-Leggett model: in the linear and homogeneous case they found an explicit mapping between the Vlasov-Poisson system and the Hamilton equations of the Caldeira-Leggett model. Although finding an analogous transformation valid in the fully nonlinear and non-homogeneous case is surely a very difficult task, it would put our approach on a much more rigorous basis, so that it is a line that is worth pursuing. Work is in progress along this direction.

In this paper we have considered only completely cold initial conditions, to which our theoretical approach seems best suited. However, with some small modifications this approach proves useful and effective also beyond the cold collapse: the extension to generic initial conditions for the HMF model will be presented in a forthcoming paper Giachetti et al. (2019). Work is also in progress on going beyond the HMF model, extending the approach to more complicated systems.

Apart from the possibility of making predictions on the quasi-stationary state reached after a cold collapse in the HMF model, the approach we have presented here can be considered a first, preliminary step towards a theoretical approach to collisionless relaxation that fully exploits the hierarchy of scales that is produced during the Vlasov evolution in order to derive an effective dynamical evolution, valid on sufficiently large scales, able to elucidate the relation between initial conditions and quasi-stationary states. We are admittedly still far from a comprehensive theory of violent relaxation predicting the shape of galaxies, yet the approach we have described here might give some hints towards this goal. Finally, it is worth recalling that the quasi-stationary state reached after the violent relaxation might not be the whole story, as it has been recently shown by means of numerical simulations of self-gravitating particles that structures like spiral arms and rings very similar to those observed in real galaxies may appear as long-living transients during the violent relaxation itself, when the initial conditions are not completely symmetric and the initial angular momentum does not vanish Benhaiem et al. (2017). This is a further indication that a deeper theoretical understanding of violent relaxation in long-range-interacting system is needed and could prove fruitful.

Appendix A Virial theorem, thermal equilibrium and the leading-order prediction of the magnetization

As mentioned in Sec. III.2, the bare leading-order prediction $\overline{m}$ of the magnetization in the quasi-stationary state after a cold collapse given by Eq. (37) could have been found using the virial theorem. This follows from the fact that at the leading order the mean-field potential $U(\vartheta)$ is harmonic,

[TABLE]

applying the virial theorem to the quasi-stationary state we thus have

[TABLE]

and using $\overline{m}=1-\langle\vartheta^{2}\rangle/2$ and the conservation of energy $\varepsilon=\langle p^{2}\rangle/2-\overline{m}^{2}/2$ we can write

[TABLE]

In a cold collapse $\varepsilon=-m_{0}^{2}/2$, and inserting the latter into Eq. (83) we have

[TABLE]

whose only positive solution is

[TABLE]

that is, Eq. (37).

It is interesting to compare the latter result and the actual values of $\overline{m}$ found in numerical simulations with the magnetization in thermal equilibrium. As shown in Antoni and Ruffo (1995); Campa et al. (2014) the magnetization $m_{\text{eq}}$ in thermal equilibrium is implicitly given as a function of the total energy $\varepsilon$ by the solution of the equation

[TABLE]

where $I_{k}(x)$ is the modified Bessel function of order $k$ and $\beta$ is defined by

[TABLE]

The magnetization in thermal equilibrium after a cold collapse as a function of the initial magnetization, $m_{\text{eq}}(m_{0})$, is thus obtained imposing $\varepsilon=-m_{0}^{2}/2$ in Eq. (87). The thermal equilibrium prediction is plotted in Fig. 10 together with the leading-order prediction $\overline{m}$ given by Eq. (37) and compared to the simulation data. It is apparent that the numerical results obtained starting from initial conditions where a Gaussian overdensity is superimposed on a uniform background yield magnetizations in the quasi-stationary state after the cold collapse that are definitely different from the thermal predictions, although the difference is not very large, while cold waterbag initial conditions lead to quasi-stationary states where the magnetization might be consistent with the thermal equilibrium value, at least for sufficiently high or sufficiently low initial magnetizations. This does not mean, however, that such quasi-stationary states are close to thermal: it only means that the differences between the quasi-stationary state and the thermal state show up when we look at finer scales of the distribution function and are not apparent at the larger scale captured by $m$.

Appendix B Solution of the equations of motion of the auxiliary oscillators

Here we describe how to solve the equations of motion of the auxiliary oscillators introduced in Sec. III.3 to model the contribution of the higher-order moments to the dynamics of the low-order ones, i.e., Eqs. (53), that we rewrite here for convenience,

[TABLE]

where $k=1,\ldots,M$.

The solution is a direct adaptation of that used to solve the analogous equations for the classical Caldeira-Leggett model, where $g(x)=x$, to the case of a generic function $g(x)$. We define the Laplace transform $\tilde{f}(s)$ of a generic function of time $f(t)$ as follows

[TABLE]

so that

[TABLE]

where $s\in\mathbb{R}$ is larger than the real part of any pole of $\tilde{f}(s)$. If we Laplace-transform Eqs. (88) we get

[TABLE]

where $\mathcal{L}[g(x)]$ is the Laplace transform of the composite function $g(x(t))$. Solving for $\tilde{q}_{k}(s)$ we obtain

[TABLE]

so that, coming back to the time domain,

[TABLE]

Recalling that, for $t>0$,

[TABLE]

we can write

[TABLE]

where

[TABLE]

The second term on the r.h.s. of Eq. (95) can be written as

[TABLE]

and since, using the convolution theorem, we can write

[TABLE]

where $\Theta(x)$ is the Heaviside step function, substituting the above result into Eq. (97) and then back into Eq. (95) we finally obtain

[TABLE]