Relaxation Processes in Long-Range Lattices

TL;DR

This paper investigates how long-range interactions affect the relaxation times in lattice systems, revealing polynomial scaling with system size and the importance of non-Markovian effects in the relaxation process.

Contribution

It introduces a kinetic equation for long-range lattices that explains the observed timescales and interaction thresholds, highlighting the role of non-Markovian effects.

Findings

Relaxation times scale polynomially with system size.

A kinetic equation explains the interaction range threshold.

Non-Markovian effects significantly influence relaxation dynamics.

Abstract

The relaxation to equilibrium of lattice systems with long-range interactions is investigated. The timescales involved depend polynomially on the system size, potentially leading to diverging equilibration times. A kinetic equation for long-range lattices is proposed, which explain these timescales as well as a threshold in the interaction range reported in [Phys. Rev. Lett. 110, 170603 (2013)]. Non-Markovian effects are shown to play an important role in the relaxation of systems of up to thousands of particles.

Equations

$$H(\mathbf{p},\theta )=\sum _{k=1}^N\frac{p_k^2}{2}+\frac{1}{2\tilde{N}}\sum _{\begin{array}{c}j,k=1\\ (j\ne k)\end{array}}^N\frac{v({\theta }_j-{\theta }_k)}{r_{jk}^{\alpha }},$$

$$\frac{\partial f_N}{\partial t}$$

$$f_s(i_1,\dots ,i_s)\equiv \int \prod _{j\ne i_1,\dots ,i_s}\mathrm{d}j{f}_N(1,\dots ,N).$$

$$\frac{\partial }{\partial t}{f}_s(J_1)$$

$$f_2(j,k)$$

$$f_3(j,k,l)$$

$$\frac{\partial }{\partial t}{f}_1(j)=-p_j\frac{\partial }{\partial {\theta }_j}{f}_1(j)+\frac{1}{\tilde{N}}\sum _{\begin{array}{c}k=1\\ (k\ne j)\end{array}}^N\int \mathrm{d}k\frac{1}{r_{jk}^{\alpha }}\frac{\partial v_{jk}}{\partial {\theta }_j}\frac{\partial }{\partial p_j}\left[f_1(j)f_1(k)+g_2(j,k)\right].$$

$$\frac{\partial }{\partial t}{f}_2(j,k)$$

$$\left(\frac{\partial }{\partial t}+p_j\frac{\partial }{\partial {\theta }_j}+p_k\frac{\partial }{\partial {\theta }_k}\right)g_2(j,k)=\frac{1}{\tilde{N}}\frac{1}{r_{jk}^{\alpha }}\frac{\partial v_{jk}}{\partial {\theta }_j}{\partial }_{jk}\left[f_1(j)f_1(k)+g_2(j,k)\right]$$

$$g_2(j,k)=\lambda g_2^{(1)}(j,k)+{\lambda }^2{g}_2^{(2)}(j,k)+\mathcal{O}\left({\lambda }^3\right),$$

$$g_3(j,k,l)={\lambda }^2{g}_3^{(2)}(j,k,l)+\mathcal{O}\left({\lambda }^3\right).$$

$$\frac{\partial }{\partial t}{f}_1(j)$$

$$\frac{\partial }{\partial t}$$

$$\frac{\partial }{\partial t}{f}_1(j)=\frac{1}{\tilde{N}}\sum _{\begin{array}{c}k=1\\ (k\ne j)\end{array}}^N\frac{1}{r_{jk}^{\alpha }}\frac{\partial }{\partial p_j}\int \mathrm{d}k\frac{\partial v_{jk}}{\partial {\theta }_j}{g}_2(j,k).$$

$$\left(\frac{\partial }{\partial t}+p_j\frac{\partial }{\partial {\theta }_j}+p_k\frac{\partial }{\partial {\theta }_k}\right)g_2^{(1)}(j,k)=\frac{1}{\tilde{N}}\frac{1}{r_{jk}^{\alpha }}\frac{\partial v_{jk}}{\partial {\theta }_j}{\partial }_{jk}{f}_1(j)f_1(k).$$

$$g_2^{(1)}(j,k;t)=e^{(-p_j\partial /\partial {\theta }_j-p_k\partial /\partial {\theta }_k)t}{g}_2^{(1)}(j,k;0)$$

$$f_1(j;t-\tau )f_1(k;t-\tau )=e^{(p_j\partial /\partial {\theta }_j+p_k\partial /\partial {\theta }_k)\tau }{f}_1(j;t)f_1(k;t)+\mathcal{O}\left(\lambda \right),$$

$$e^{(-p_j\partial /\partial {\theta }_j-p_k\partial /\partial {\theta }_k)\tau }$$

$$g_2^{(1)}(j,k;t)=\frac{1}{\tilde{N}}{\int }_0^t\mathrm{d}\tau \frac{1}{r_{jk}^{\alpha }}{v}^{\prime }({\theta }_{jk}-p_{jk}\tau ){\partial }_{jk}{f}_1(j;t-\tau )f_1(k;t-\tau ).$$

$$g_2^{(1)}(j,k)=g_2^{(1)}(k,j).$$

$$\frac{\partial }{\partial t}{f}_1(p_j;t)=\frac{\lambda }{{\tilde{N}}^2}\sum _{\begin{array}{c}k=1\\ (k\ne j)\end{array}}^N\frac{1}{r_{jk}^{2\alpha }}{\int }_0^t\mathrm{d}\tau \int \mathrm{d}{p}_k{\partial }_{jk}⟨{\mathcal{F}}_{jk}(0){\mathcal{F}}_{jk}(\tau )⟩{\partial }_{jk}{f}_1(p_j;t-\tau )f_1(p_k;t-\tau ),$$

$$⟨{\mathcal{F}}_{jk}(0){\mathcal{F}}_{jk}(\tau )⟩={\int }_{-\pi }^{\pi }\mathrm{d}{\theta }_k\frac{\partial }{\partial {\theta }_j}v({\theta }_{jk})\frac{\partial }{\partial {\theta }_j}v({\theta }_{jk}-p_{jk}\tau ).$$

$$v(\theta )=\sum _{n=-\infty }^{\infty }\tilde{v}(n)e^{-in\theta },$$

$$\tilde{v}(n)=\frac{1}{2\pi }{\int }_{-\pi }^{\pi }\mathrm{d}\theta v(\theta )e^{in\theta },$$

$${\tilde{g}}_2^{(l)}(n,m;p_i,p_j;t)={\int }_{-\pi }^{\pi }\mathrm{d}{\theta }_i{\int }_{-\pi }^{\pi }\mathrm{d}{\theta }_j{g}_2^{(l)}(i,j;t)e^{in{\theta }_i}{e}^{im{\theta }_j}.$$

$${\tilde{g}}_2^{(1)}(n,m;p_j,p_k;t)$$

$$\frac{\partial }{\partial t}{f}_1(p_j;t)$$

$${\int }_0^{\infty }\mathrm{d}\tau e^{ia\tau }=\pi {\delta }_{+}(a)=\pi \delta (a)+i\mathcal{P}\left(\frac{1}{a}\right),$$

$$\frac{\partial }{\partial t}{f}_1(p;t)=\frac{2{\pi }^2\lambda }{\tilde{N}}\sum _{\begin{array}{c}l=1\\ (l\ne j)\end{array}}^N\frac{1}{|x_j-x_l{|}^{2\alpha }}\sum _{n=-\infty }^{\infty }{\left|n\tilde{v}(n)\right|}^2$$

$$\times \frac{\partial }{\partial p}{\int }_{-\infty }^{\infty }\mathrm{d}{p}^{\prime }\delta (n(p-p^{\prime }))\left(\frac{\partial }{\partial p}-\frac{\partial }{\partial p^{\prime }}\right)f_1(p;t)f_1(p^{\prime };t).$$

Full text

I Introduction

Systems with long-range interactions are present at various scales in nature, from astrophysics to atomic scales Campa, Dauxois, and Ruffo (2009). Apart from the peculiarities in their statistical equilibrium, such as negative specific heat Lynden-Bell and Lynden-Bell (1977) or ensemble inequivalence Barré, Mukamel, and Ruffo (2001), an important feature of these systems are the evolution time scales that arise. For example, governed by the gravitational force, galactic dynamics is characterized by its inefficiency at redistributing kinetic energy Binney and Tremaine (2008), resulting in relaxation times which increase with the system size. At the microscopic scale, long-range effects rise in the presence of light-mediated (dipole-dipole) interactions between cold atoms Chalony et al. (2013); Barré, Marcos, and Wilkowski (2014); Aspelmeyer, Kippenberg, and Marquardt (2014), and size-dependent equilibration times were also reported Guerin, Araújo, and Kaiser (2016).

Long-range interactions may manifest between either the internal or external degrees of freedom of the particles, which leads to different scaling laws for the relaxation times Gabrielli, Joyce, and Marcos (2010). In this context, a new platform for the study of long-range lattices, i. e. for particles at fixed positions but with dynamical internal variables, has emerged with the realization of ion chains with tunable interactions Britton et al. (2012); Senko et al. (2015). There, the coupling between the particles decay as an adjustable power-law of their mutual distance, which makes ion traps a versatile platform to study the crossover from short- to long-range, including in the quantum realm Richerme et al. (2014); Jurcevic et al. (2014).

From a theoretical point of view, explicit scaling laws were derived for the a few specific interactions: Astrophysical systems present relaxation times that scale as $N/\log N$, with $N$ the number of interacting bodies; For the so-called Hamiltonian mean-field model Antoni and Ruffo (1995), with infinite-range interactions, a scaling with $N$ for non-homogeneous and $N^{2}$ for homogeneous states Yamaguchi (2003); Bouchet and Dauxois (2005); da Rocha Filho et al. (2014) was derived. The existence of size-dependent time scales was reported in the more general context of long-range lattices with $1/r^{\alpha}$ interactions ($\alpha$ smaller than $d$ the dimension of the system for the long-range case), both in the classical and in the quantum context Kastner (2011); Bachelard and Kastner (2013). In particular, a puzzling threshold in the scaling laws was reported at $\alpha=d/2$, which does not correspond to the transition from long-range to short-range systems.

In this work we derive explicitly the scaling of the relaxation time of one-dimensional ($d=1$) long-range classical lattices. The threshold at $\alpha=d/2$ originates in two-particle correlation terms, and it manifests both in the non-Markovian (small size) and Markovian (large size) regimes. Indeed, the derivation of a kinetic equation for long-range lattices allows to identify the specific scaling of the contribution of non-Markovian terms, which are all the more important for fast-decaying interactions terms ($\alpha\to d$).

II Kinetic equation for long-range lattices

II.1 Long-range lattices

Let us consider a one-dimensional lattice of $N$ particles at coordinates $x_{k}=k\xi$, $k=1,..N$ and $\xi$ the lattice step. Each particle has an internal degree of freedom represented by an angular coordinate $\theta_{k}$ and its canonically conjugate momentum $p_{k}$, with a Hamiltonian of the form:

[TABLE]

where $v(\theta_{j}-\theta_{k})/r_{jk}^{\alpha}$ is the pair interaction potential for particles $j$ and $k$, separated by the distance $r_{jk}=\|x_{j}-x_{k}\|$, with $v(0)=0$. The factor $1/\tilde{N}\equiv 1/N^{1-\alpha}$ in front of the potential energy term is introduced (e. g. by a renormalization of time) in order to obtain an extensive total energy and to properly define the passage to a Vlasov equation description. Such interaction was investigated as a generalization of the Hamiltonian Mean-Field model Anteneodo and Tsallis (1998), and describes more generally lattices of fixed particles interacting through their internal degrees of freedom, as for example in light scattering problems Lehmberg (1970); Máximo et al. (2015); Hill, Sanders, and Deng (2017).

The evolution of the system and its relaxation are studied introducing the $N$-particle distribution function $f_{N}(1,2,\ldots,N;t)$, defined such that $f_{N}\,d1\cdots dN$ is the probability for particle $k$ ($k=1,2,\ldots N$) to be in the volume element $dk\equiv d\theta_{k}dp_{k}$ that contains the point $(\theta_{k},p_{k})$ in phase-space, at time $t$. We note that since all particles are localized on the lattice points, the distribution $f_{N}$ is not required to be invariant with respect to particle permutations, which must be taken into account in the determination of the generalized form of the BBGKY hierarchy.

The Liouville equations for the lattice system then writes as

[TABLE]

with $v_{kl}\equiv v(\theta_{k}-\theta_{l})$.

II.2 Generalized BBGKY hierarchy

Since each particle is distinguished by its position on the lattice, in order to define the $s$-particle reduced distribution function we have to specify which particle positions and momenta are integrated out. Thus, the $s$-particle reduced distribution function depending on particles of indices $i_{1},\ldots,i_{s}$ is defined by:

[TABLE]

In Eq. (3) and from now on we keep the time dependence implicit, except where necessary. For the one-particle distribution function ($s=1$), the specification on which particle phase-space coordinates it depends is not required, provided boundary effects are negligible. Note that these reduced distribution functions are not invariant by permutation of their arguments in a lattice, in general.

Let us now consider two disjoint sets of different particle indices $J_{1}$ and $J_{2}$, such that $J_{1}$ has $s$ indices and $J_{2}$ has $N-s$ indices. By integrating Eq. (2) over the position and momentum variables of particles in $J_{2}$ we obtain:

[TABLE]

where we introduced the cross derivative $\partial_{jk}\equiv\partial/\partial p_{j}-\partial/\partial p_{k}$. By eliminating vanishing surface terms in Eq. (4) and using Eq. (3) we obtain the generalized form of the BBGKY hierarchy:

[TABLE]

where $J_{1}^{(l)}\equiv J_{1}\cup\{l\}$. The next step towards the kinetic equation is the introduction of the $s$-particle correlation functions, as discussed below.

II.3 Irreducible cluster representation and prototypical kinetic equation

The $s$-particle reduced distribution function can be decomposed into a purely uncorrelated part and contributions from $s$-particle correlation functions $g_{s}(j_{1},\ldots,j_{s})$ as Balescu (1997):

[TABLE]

and so on, where $P(j,k,l)$ stands for all different permutation of particles $j,k,l$. All permutations must be considered due to the absence of permutation invariance. Since we are considering correlation among different particles, in the sum in Eq. (6) only terms with $k\neq l$ must be considered for $g_{2}(k,l)$. In order to determine the order of magnitude of the functions $g_{s}$, we note that the correlation between two particles require the interaction between these particles, the correlation function $g_{2}$ is of order $\tilde{N}^{-1}$, the order of the interaction. Similarly, $g_{3}$ requires at least two pair interactions, and is therefore of order $\tilde{N}^{-2}$. More generally, the correlation function $g_{s}$ is of order $\tilde{N}^{-s+1}$.

The kinetic equation is a closed equation for the one-particle reduced function, and therefore as a preliminary step we replace Eq. (6) into Eq. (5) for $s=1$ ($J_{1}={j}$) to obtain:

[TABLE]

Proceeding similarly for $s=2$ (with $J_{1}=\{j,k\}$) leads to:

[TABLE]

Hence, in order to have a closed equation for $f_{1}$ the correlation function $g_{2}$ must be expressed in terms of $f_{1}$. Replacing the irreducible cluster expansion in Eq. (6) into Eq. (8) yields:

[TABLE]

where we used $\int_{-\pi}^{\pi}\mathrm{d}\theta_{l}(\partial v_{kl}/\partial\theta_{l})=0$.

By considering only terms up to order $1/\tilde{N}$, one obtains a closed form equation for $g_{2}$. Its solution, inserted in Eq. (7), and after some simplifying assumptions (see below) leads to the generalization of the Balescu-Lenard equation Balescu (1997). It is known that, in the mean-field case ($\alpha=0$), the collision term in the kinetic equation vanishes at order $1/\tilde{N}$ for a homogeneous state Sano (2012). Consequently, one must consider terms of order $1/\tilde{N}^{2}$ by including three-particle correlations to describe the relaxation process. Nevertheless, solving the resulting equation for $g_{3}$ is a daunting task. Here we rather follow the weak-coupling approach of Ref. da Rocha Filho et al. (2014): it consists in expanding the hierarchy in orders of a weak coupling, considering an inter-particle potential of order $\lambda\ll 1$, but at the same time retaining only dominant terms in $1/\tilde{N}$. In the case of long-range systems, where each particle is coupled to all others, it refers to disordered regimes, where the macroscopic degrees of freedom that describe the coupling (such as the magnetization) are very weak on average. For a detailed discussion of the weak coupling approach, the reader is referred to Ref. da Rocha Filho et al. (2014); Fouvry, Bar-Or, and Chavanis (2019). This leads to a generalization of the Landau equation for lattice systems. The correlations functions are then expanded as:

[TABLE]

Before determining a kinetic equation at leading order in $\lambda$ and $1/\tilde{N}$, let us first consider the mean-field limit for the lattice system.

II.4 Vlasov equation

The mean field description is obtained in the limit $\tilde{N}\rightarrow\infty$ and is equivalent to neglecting two-particle correlations $g_{2}={\cal O}(\tilde{N}^{-1})$ in Eq. (7). It results in the generalized form of the Vlasov equation Bachelard et al. (2011):

[TABLE]

after eliminating a vanishing surface term. The continuous limit on the lattice is obtained by performing $d\rightarrow 0$ at constant lattice length ${\cal N}=Nd$:

[TABLE]

In Eq. (II.4) collisional (granularity) effects are neglected. For states homogeneous in $\theta$, one has $\partial f_{1}(p)/\partial t=0$, i. e., the distribution function has no dynamics. An evolution is recovered by introducing the corrections due to collisions, which is the purpose of the next sections. The less trivial case of non–homogeneous states can be addressed using action-angle variables Chavanis (2007), yet it will not be considered here.

II.5 Weak coupling limit

For the homogeneous states under consideration, the prototypical kinetic equation in Eq. (7) assumes the form:

[TABLE]

We now proceed to obtained a closed form for the two-particle correlation function $g_{2}$ in the weak-coupling limit by performing the expansion in power series on the parameter $\lambda$ given in Eq. (10). This represents a generalization of the Landau equation for lattice systems.

By replacing Eq. (10) into Eq. (9) and retaining only dominant terms in $\lambda$ we obtain:

[TABLE]

As discussed in Ref. da Rocha Filho et al. (2014), this approximation is justified by the fact that the total effective force on a given particle is weak in a non-magnetized state. Equation (14) can be solved in the form of a convolution as:

[TABLE]

The first term on the right-hand side of this equation is a transient contribution from the correlations at the initial time, and can be discarded after a short transient Balescu (1997). Since for a homogeneous state, the angle variables $\theta_{j}$ evolve in a free (ballistic) motion, up to corrections of order $\lambda$, we obtain

[TABLE]

where we have used the free time propagator $\exp(-t\>\partial/\partial\theta)$. Using the identities

[TABLE]

where $\theta_{jk}\equiv\theta_{j}-\theta_{k}$ and $p_{jk}\equiv p_{j}-p_{k}$, we rewrite equation (II.5) as

[TABLE]

In particular, the latter equation implies that

[TABLE]

Plugging Eq. (18) into Eq. (13), and using that fact that the mean-field force cancels for a homogeneous state, we obtain:

[TABLE]

where the force auto-correlation function has been introduced:

[TABLE]

A kinetic equation in the form of Eq. (20) is clearly non-Markovian, and one usually goes a step further: if the force auto-correlation decays to zero faster than any significant change in the one-particle distribution function, then the time integral in Eq. (20) can be extended to infinity and we can set $f_{1}(p;t-\tau)\rightarrow f_{1}(p;t)$. This results in a Markovian dynamical evolution Balescu (1997); Lourenço and da Rocha Filho (2015). The validity of this approximation is discussed in details in Sec. III.

Let us first show that, assuming the Markovian regime is reached, the first-order contribution to the kinetic equation is dominant. To this end, expressions (20) and (21) are cast in the Fourier space, using the following series for the potential $v(\theta)$:

[TABLE]

where the Fourier coefficients are given by:

[TABLE]

and the correlation function by:

[TABLE]

From Eq. (18) we obtain:

[TABLE]

where $\delta_{n,m}$ refers to the Kronecker delta. As $g^{(1)}_{2}(j,k)$ is the correlation function between two different particles, it is only defined for $j\neq k$. Equation (13) then converts into:

[TABLE]

where we have used the property $v(-n)=v^{*}(n)=v(n)$ and the Cauchy integral:

[TABLE]

with ${\cal P}(x)$ the Cauchy principal part of $x$ (an odd function). The imaginary part of the right-hand side of Eq. (26) vanishes and setting $p_{j}=p$ and $p_{l}=p^{\prime}$, we obtain:

[TABLE]

The right-hand side of Eq. (28) vanishes due to the Dirac delta function, just as for the Balescu-Lenard equation for the mean-field ($\alpha=0$) case da Rocha Filho et al. (2014). Thus, under the hypothesis of a Markovian dynamics, one needs to go one order further in $\lambda$ to determine the kinetic equation.

II.6 Higher-order contributions

Neglecting the two last terms in the right-hand side of Eq. (9), i. e., retaining only terms up to order $1/\tilde{N}$ and neglecting collective effects, leads to an integral equation for $g_{2}$. Its solution, inserted into Eq. (13), will provide a generalization of the Balescu-Lenard equation for lattice systems Balescu (1997); Lenard (1960). Once more, it results in a vanishing collisional integral. Indeed, the only contributions to $g_{2}$ with a non-vanishing contribution to the kinetic equation originate in the three-particle correlation function $g_{3}$ of Eq. (13), which are of order $1/\tilde{N}^{2}$.

Using algebra package MAPLE map , the leading order contribution for $g_{3}$ are determined to be

[TABLE]

where we have introduced

[TABLE]

Inserting the three-particle correlation $g_{3}$ from Eq. 29 into Eq. (9), one can obtain a closed form for the kinetic equation, to the second order in the coupling. This approach was used in Ref. da Rocha Filho et al. (2014) to obtain a kinetic equation for the Hamiltonian mean-field model. Since we are here interested in the scaling of the relaxation time, rather than its exact expression, we now proceed with discussing the Markov approximation.

III Relaxation times and the Markov approximation

III.1 The Markovian hypothesis

The Markov approximation consists in assuming that the force auto-correlation function $C_{p}$ in Eq. (21) vanishes over a time such that the one-particle distribution function does not change significantly. If it does not apply, auto-correlation terms in the force distribution makes that the right-hand term of Eq. (20) contributes substantially, and may even dominate over the higher-order contributions discussed in Sec. II.6.

The (non-)Markovian nature of the dynamics strongly depends on the model, especially on the system size Balescu (1997), but also on the initial conditions. Let us here consider the cosine potential from the $\alpha$XY chain Anteneodo and Tsallis (1998); Bachelard and Kastner (2013):

[TABLE]

for which we study numerically the evolution of the force auto-correlation. The simulations are realized by integration of the equation of motion, where the initial state is a random realization of a state homogeneous in angles and bounded in momentum:

[TABLE]

The constant $p_{0}$ allows to choose an energy such that the system remains in a non-magnetized phase at all times Bachelard et al. (2011).

The Markovian nature of the dynamics is tested by comparing the dynamics of the average auto-correlation function of the force $C_{F}$ with the macroscopic evolution of the distribution function, and is assessed from molecular dynamics simulations. This average auto-correlation function $C_{F}$ is defined by:

[TABLE]

so it is normalized at time zero. Here we implemented the numerical solutions of the Hamiltonian equations of motion using a fourth-order symplectic integrator Yoshida (1990) in a parallel implementation as described in Ref. da Rocha Filho (2014). Because we are considering homogeneous states, we monitor the first moments of the momentum distribution $M_{k}\equiv\langle p^{k}\rangle$. The odd moments fluctuate around zero, and the second moment $M_{2}$ does not evolve substantially due to energy conservation. We thus focus on $M_{4}$ (other even momentum $M_{6}$, $M_{8}$, $\ldots$ support the same conclusions but higher flutuations with higher order), which typically present a slow dynamics during the relaxation process Bachelard and Kastner (2013). Figure 1 shows the force auto-correlation as a function of time for an increasing system size $N$ (ranging from $256$ to $16384$) and several values of $\alpha$. The decay of the force auto-correlation does not appear to depend significantly on the system size, and little on the interaction range $\alpha$. We thus consider that it cancels around $t_{FC}\approx 30$, and consider that the Markovian approximation is valid if the fourth momentum $M_{4}$ does not vary significantly in this time interval, i. e. that $M(t_{FC})-M_{4}(0)$ is small.

Differently from the force auto-correlation, the momentum $M_{4}$ presents an evolution that depends strongly on both the system size and the interaction range. Figure 2 shows the variation of $M_{4}$ from its initial value as a function of time. Larger $N$ and smaller $\alpha$ values are associated to a much slower evolution of the momentum distribution. Thus, the Markov approximation is reached for increasing system sizes, and the larger the range of the interaction (that is, the smaller $\alpha$), the smaller the system size required to reach it.

The (non-)Markovian nature of the dynamics has been characterized in details for the mean-field case ($\alpha=0$) in Ref. Lourenço and da Rocha Filho (2015). We here simply evaluate the typical relative change of the fourth momentum during the time over which the force auto-correlation cancels, by characterizing the average growth of the momentum (through a linear fit of $M_{4}(t)$) as given by the normalized quantity:

[TABLE]

where $\left\langle\cdots\right\rangle$ stands for a time average between $t=0$ and $t_{FC}$. Its evolution is shown in Fig. 3. Thus, close to the threshold to short-range interactions ($\alpha=1$), the dynamics requires huge system sizes to reach the Markovian regime whereas close to the infinite-range case ($\alpha=0$), the Markov approximation is already valid for modest system sizes. We note that such results are highly dependent on the model, its parameters, dimensionality, etc.

III.2 Non-Markovian regime

Hence, if over the momentum distribution changes significantly over the time scale during which the force auto-correlation is non-zero, the dynamics must be considered non-Markovian: The lower order term in the right-hand side of Eq. (20) does not vanish and contributes significantly to the single-particle distribution evolution. The timescales of the evolution of the single-particle distribution $f_{1}$ with time originates in the sum $\sum_{k}r^{-2\alpha}_{jk}$, whose scaling with the system size changes with the interaction range. Assuming periodic boundary conditions, for which the distance is given by $r_{jk}=\xi\min(|j-k|,N-|j-k|)$, in the large-$N$ limit one obtains

[TABLE]

Inserting the $1/\tilde{N}^{2}\sim N^{2\alpha-2}$ in the above expression leads to the following scaling for the relaxation time, as due to the non-Markovian contribution:

[TABLE]

The threshold at $\alpha=d/2$ reported in Ref. Bachelard and Kastner (2013) is thus already present in the non-Markovian contribution, from two-particle collisions term (lower order term). This threshold is also clearly visible in the variation of the momentum distribution presented in Fig. 3.

III.3 Markovian regime

Prediction (36) fits only partially the numerical findings for the one-dimensional classical chain reported in Ref. Bachelard and Kastner (2013), where a scaling of the order of $\tau\sim N^{1.5}$ was observed for $\alpha<0.5$, before it decays roughly as $\tau\sim N^{2.5-2\alpha}$ for $1/2<\alpha<1$. In other words, the relaxation times appears to scale with a factor of order $\sqrt{N}$ larger than as predicted by the non-Markovian contribution.

The origin of this discrepancy can be found in the competition between non-Markovian and Markovian contributions. Indeed, as discussed extensively for the mean-field case in Ref. Lourenço and da Rocha Filho (2015), this competition leads to the observation of scaling laws that are intermediate between the two regimes: $\tau\sim N$ for the non-Markovian terms and $\tau\sim N^{2}$ for the Markovian terms, for $\alpha=0$. The same effect happens for long-range lattices with any $\alpha$, for the system sizes achievable by numerical simulations.

The inspection of Eqs. (9) and (29) reveals that similar $\sum_{k}r^{2\alpha}_{jk}$ terms appear for the Markovian contribution to the dynamics, although the variety of terms prevents a thorough analysis of the contribution of $g_{3}$. Nevertheless, the presence of these sums is an indicator that the threshold at $\alpha=d/2$ will be preserved at larger system sizes, as supported by the variation of the momentum $M_{4}$ presented in Fig. 3. There, simulations realized with $N=262144$ (a factor $64$ larger than the relaxation study presented in Ref. Kastner (2011)) still exhibit the threshold, although such a large system size prevents from realizing a full equilibration study, instead restricting us to the initial stage of the dynamics. Hence, the non-Markovian and Markovian terms present different scaling laws, both exhibiting the $\alpha=d/2$ threshold, and systems of thousands of particles being in-between these two regimes.

We note that while a rigorous evaluation of the Markovian nature of the dynamics would require a dedicated, and arduous, study of each reduced distribution and auto-correlation function involved in Eq.(9). This task is beyond the scope of this manuscript, and the good agreement between the prediction of Eq.(36) and the numerical results suggest that these terms do not contain scalings that change significantly the (non-)Markovian nature of the dynamics.

IV Conclusion

In this work, we have derived a kinetic equation for long-range lattices with power-law interactions. This allowed us to identify contributions that change of scaling with the system size at $\alpha=d/2$, both in the non-Markovian and in the Markovian terms. Lattices of hundreds to thousands of particles are subject to a competition between these two kinds of terms, at least in the context of the homogeneous states of the $\alpha$XY chain considered here.

Establishing a kinetic equations for these classical systems is an important step toward the understanding of their relaxation to equilibrium. The presence of a threshold may lead to a further distinction between the classification of the interaction range, just as it was done between short- and long-range Campa, Dauxois, and Ruffo (2009), or between lattices and moving bodies Gabrielli, Joyce, and Marcos (2010).

The tools presented here could in principle be generalized to address quantum systems, by using the quantum analog of the BBGKY hierarchy, e. g. from the Wigner function representation of the quantum equation of motion Liboff (2003), and study thermalization processes in quantum systems. In this context, the emergence of flexible experimental platforms such as cold atom or trapped ion setups, where the engineering of the Hamiltonians and the reduction of decoherence channels make constant progresses, is a promising tool to investigate the equilibration processes of both classical and quantum systems.

Acknowledgements.

We thank fruitful discussions with M. Kastner. R. B. benefited from Grants from São Paulo Research Foundation (FAPESP) (Grant Nos. 2014/01491-0 and 2018/01447-2). R. B. and T. M. R. F benefited from Grants from the National Council for Scientific and Technological Development (CNPq) Grant No.  305842/2017-0.