S De Bièvre

P E Parris
email: parris@mst.edu

Equilibration, generalized equipartition, and diffusion in dynamical Lorentz gases

We demonstrate approach to thermal equilibrium in the fully Hamiltonian evolution of a dynamical Lorentz gas, by which we mean an ensemble of particles moving through a d-dimensional array of fixed soft scatterers that each possess an internal harmonic or anharmonic degree of freedom to which moving particles locally couple. We analytically predict, and numerically confirm, that the momentum distribution of the moving particles approaches a Maxwell-Boltzmann distribution at a certain temperature T , provided that they are initially fast and the scatterers are in a sufficiently energetic but otherwise arbitrary stationary state of their free dynamics-they need not be in a state of thermal equilibrium. The temperature T to which the particles equilibrate obeys a generalized equipartition relation, in which the associated thermal energy k B T is equal to an appropriately defined average of the scatterers' kinetic energy. In the equilibrated state, particle motion is diffusive.

Introduction

While it is a well established fact that large systems converge to thermal equilibrium starting from a more or less arbitrary initial state, few or no general dynamical mechanisms responsible for this approach to equilibrium have been identified. An exception is the case where a small system S with Hamiltonian H S and only a few degrees of freedom is weakly coupled to a large "reservoir" R with Hamiltonian H R and many degrees of freedom. Under these circumstances, the small system will converge from an arbitrary initial state ρ 0 to a state ρ S ∼ exp(-H S /k B T ) of thermal equilibrium provided that the reservoir is itself initially in a thermal state ρ R ∼ exp(-H R /k B T ). This phenomena of return to equilibrium has been proven rigourously in a number of (relatively simple) systems [START_REF] Jakšić | Ergodic properties of classical dissipative systems[END_REF][START_REF] Bach | Return to equilibrium[END_REF][START_REF] Dereziński | Return to equilibrium for pauli-fierz systems[END_REF]. In such systems return to equilibrium occurs because, roughly speaking, the smaller subsystem acts as a small perturbation of the larger one. Strong stability properties of the thermal states ρ ∼ exp(-H/k B T ) then force the coupled system to a joint equilibrium state characterized by the initial temperature of the reservoir. Studies of this kind however, do not address the question of how the past evolution of the reservoir dynamically led to its degrees of freedom being in thermal equilibrium in the first place.

In this paper we demonstrate approach to equilibrium within a class of fully Hamiltonian models that we refer to as dynamical Lorentz gases. In these models, each member of an ensemble of particles moves through an array of localized independent scatterers, each possessing an internal harmonic or anharmonic degree of freedom to which the particle locally couples [see (1.1)]. The degrees of freedom of the scatterers are initially drawn from a stationary but not necessarily thermal state of their uncoupled dynamics. Our main result, which we predict analytically and confirm numerically, is that asymptotically in time, as a result of repeated independent scattering events between a particle and the degrees of freedom of the medium, the particle momentum distribution ρ t (p) will, starting from an arbitrary initial distribution ρ 0 (p) of sufficiently large average mean speed, inevitably converge to a Maxwell-Boltzmann distribution characterized by a well-defined, nonzero temperature T . This constitutes a true approach (rather than return) to equilibrium since this emergence of the Maxwell-Boltzmann distribution occurs even when the scattering medium, which serves as a reservoir for the particle, is not itself in a state of thermal equilibrium to begin with. Our result therefore helps to explain, from a purely dynamical point of view, the robustness and ubiquity of the Boltzmann factor that characterizes states of thermal equilibrium.

In our analysis, the effective temperature T that dynamically emerges from the interaction between the moving particle and the scatterers leads to a generalized equipartition relation [see (4.1)], in which the thermal energy k B T /2 per degree of freedom of the particle is equal to a suitably-defined coupling-weighted average of the kinetic energy of the scatterers. In two particular cases it reduces to the standard (i.e., un-weighted) equipartition result: (i) when the coupling of the particle to a scatterer is linear in the displacement from mechanical equilibrium of the latter's coordinate (independent of the particular stationary distribution of the scatterers), and (ii) when the scatterers are themselves initially in thermal equilibrium at a well-defined temperature T (independent of the specific coupling between the scatterers and the particle). As an example of a situation in which the standard equipartition result is not obeyed, we predict analytically, using the generalized equipartition relation, and confirm numerically that when the particle is quadratically coupled to "harmonic scatterers" that are each initially in their own micro-canonical state of fixed energy, the final thermal energy of the particle is exactly one-half the value expected from standard equipartition.

We show finally that in these models the asymptotic mean-squared displacement of an initially localized ensemble of moving particles grows linearly in time, with a diffusion constant D that depends on the effective temperature T through the relation [see (5.8)]

D ∼ T ν , with 1 2 ≤ ν < 5 2 ,
where the power ν depends on the nonlinearity of the coupling and on the anharmonicity of the scatterers. A more strongly nonlinear coupling leads to a lower value of ν and slower diffusion. Stronger scatterer anharmonicity, on the other hand, leads to a higher value of ν and faster diffusion. Our demonstration of these results relies on an approximate random walk description of the motion of the particle between and during successive collisions, which are treated as independent scattering events. It depends in particular upon a careful analysis of the interaction between the moving particles and the local scatterers, and notably of the energy exchange that occurs between them. Based on the Hamiltonian nature of the dynamics we show that, for sufficiently fast particles, although random fluctuations in the energy change of the particle (per scattering event) have a magnitude that is of order p -1 , the average energy change is negative at high energies, and of order p -2 [see (3.9), (3.7) and (3.8)]. Thus, while random fluctuations tend to induce a stochastic growth in the particle's speed to ever higher values, a weaker (but systematic) average energy loss acts as a source of dynamical friction [START_REF] Chandrasekhar | Dynamical friction. I. General considerations: the coefficient of dynamical friction[END_REF] that tends to reduce it. The relative strengths of these two processes, as encoded in the precise power laws above, are of the exact form required to dynamically drive the kinetic energy distribution of the particle asymptotically in time to a Boltzmann distribution, with a temperature T determined by the ratio of the magnitude of the fluctuations to that of the dynamical friction [see (3.24)]. We stress that, in our analysis, this fluctuation-dissipation-like relation emerges naturally from the microscopic Hamiltonian dynamics, and acts as a defining property of the temperature of the limiting Boltzmann distribution, rather than as an a posteriori property of thermal equilibrium, as it does in most treatments (see for example [START_REF] Kubo | Statistical physics[END_REF]).

The results presented here constitute a generalization in various ways of our previous work [START_REF] Silvius | Adiabaticnonadiabatic transition in the diffusive hamiltonian dynamics of a classical holstein polaron[END_REF] in collaboration with A. Silvius. There, diffusive behavior was established for a particle moving through a one-dimensional lattice of harmonic scatterers to which the particle was coupled linearly with a very particular form factor. Approach to equilibrium, although observed numerically, was not addressed in that paper. Our work here benefits from the insights gained since in [START_REF] Aguer | Classical motion in force fields with short range correlations[END_REF], where the related but very different problem of the motion of fast particles in a random time-dependent potential was studied. The potentials considered in that paper, while also generated by an array of isolated time-evolving scatterers, are non-reactive (or inert), in that they do not respond to the particle's passage. The total energy of such a system is not conserved, and the random scattering events experienced by the particle cause its kinetic energy to grow slowly but indefinitely, in the phenomenon of stochastic acceleration. As discussed above, however, when the scatterer degrees of freedom are treated dynamically, as they are in the current paper, the Hamiltonian interaction between the scatterers and the moving particle provides a way for the particle to dissipate excess energy to the medium; this mechanism completely suppresses stochastic acceleration and allows the particle to approach thermal equilibrium. In fact, in the Hamiltonian models considered here, the frictional component of the force cannot be independently made small compared to the strength of its fluctuating part. As a result, in such systems stochastic acceleration cannot be observed, even on intermediate time-scales, before equilibration sets in.

We now describe more precisely the models that form the focus of our analysis. After a dimensional re-scaling of the dynamical variables, masses, coupling constants, and time into appropriate dimensionless forms, the particlescatterer system in our model is assumed to obey the following equations of motion

ÿ(t) = -α N η(Q N (t))∇σ(y(t) -x N ) M QN (t) + U (Q N (t)) = -αη (Q N (t))σ(y(t) -x N ).      (1.1)
In these equations 1 , y(t) ∈ R d , d > 1 represents the particle position at time t, and Q N ∈ R the displacement of the internal degree of freedom associated with the scatterer centered at the fixed point x N ∈ R d . The potential U governing the uncoupled dynamics of the scatterer (i.e., of its internal degree of freedom) as well as the coupling function η are assumed to be smooth functions of their argument. More precisely, we will always assume that U and η exhibit polynomial growth of the type U (Q) ∼ |Q| r , |η(Q)| ∼ |Q| r , for some integers 0 < r ≤ r, so that in particular U is confining. The simplest case, in which η is a linear function of its argument will be referred to as "linear coupling". The locations x N , N ∈ Z d of the scattering centers can be chosen either randomly (with uniform density) or lying on a regular lattice. The form factor σ(•) appearing in the interaction terms is assumed to be a rotationally invariant smooth function, bounded by one in absolute value and supported in a ball of radius 1/2, so that the particle interacts with the scatterer at x N only when y(t) -x N ≤ 1/2. We suppose min N =M x N -x M > 1, so that the interaction regions associated with different scatterers do not overlap, and further assume that the system has a finite horizon, so that the distance over which a particle can freely travel without encountering a scatterer, is less than some fixed distance L * > 0, uniformly in time and space and independently of the direction in which it moves. Finally, α, α ≥ 0 are coupling constants that we assume to be small. A typical numerically computed trajectory y(t) is presented in Fig. 1, for a system of the above type, described more fully in Sec. 4, in which the scatterers are harmonic oscillators, the coupling is linear, and the form factor σ is equal to unity inside the circular interaction regions indicated, and vanishes everywhere else.

In the inert models studied in [START_REF] Aguer | Classical motion in force fields with short range correlations[END_REF], α = 0 so that the scatterers do not respond to the presence of the particle [see (1.1)]. In this paper we are interested in the case where they do, and take α = α. In this case the equations (1.1) are, with p = ẏ and P N = M QN , equivalent to the energy conserving equations of motion generated by the Hamiltonian

H = p 2 2 + N H scat (Q N , P N ) + α N η(Q N )σ(y -x N ),
being dimensionless. They can be obtained from a dimensional model

m d 2 ŷ dτ 2 = -α N η( QN (τ )/ )∇σ((ŷ(t) -xN )/L) M d 2 QN dτ 2 + kU ( QN (τ )/ ) = -αη ( QN (τ )/ )σ((ŷ(t) -xN )/L),
where and L are lengths, m and M masses, k, α, α energies, by introducing

ω 2 = k m 2 , M = M L 2 m 2 , α = α mω 2 2 , α = α mω 2 2
and performing the change of variables

t = ωτ, y(t) = ŷ(τ )/L, Q(t) = Q(τ )/ , xN = xN /L,
which leaves a dimensionless model governed by three independent dimensionless parameters, M, α, and α. The left panel shows an extended part of a typical trajectory, starting from the origin, computed for a particle moving through a hexagonal lattice of harmonic scatterers, with unit lattice spacing and linear coupling, as described more fully in Sec. 4. The right panel depicts a close-up view of the region indicated by an arrow in the left panel. In the close-up, the points where the particle enters and exits each circular interaction region are indicated with black dots.

where

H scat (Q, P ) = P 2 2M + U (Q). (1.2)
Obviously, it is not strictly necessary in our analysis to introduce the separate coupling parameter α. However, as we will see, retention of this extra parameter in the analysis helps to clearly bring out the source of the dynamical friction exerted on the particle by the scatterers, and to aid in its identification with the back-reaction of the scatterer. It should also be obvious that, in this class of models, evolution in time of the internal degree of freedom Q N (t) of a scatterer does not lead to a corresponding change in its location x N , nor in the size or the shape of the interaction region experienced by a moving particle that encounters it. It does, however, lead to a time-dependent change in the magnitude (and possibly the sign) of the interaction energy between the particle and the scatterer during the time that the particle traverses the interaction region associated with it. When α = 0 and αη = +∞, the scatterers form isolated and impenetrable hard spherical obstacles; in this limit the model reduces to a standard Lorentz gas. In our generalization of such a gas, αη is assumed to be finite, so the scatterers are soft; a particle with high enough energy can then enter the interaction region associated with each one. Moreover, as noted above, with α = α, the internal degree of freedom of each scatterer responds dynamically to the presence of the moving particle in an energy-conserving manner, providing a mechanism for the moving particle to dissipate excess energy to the scatterer degrees of freedom. It is these features, taken together, that motivate our reference to this class of models as dynamical Lorentz gases.

The rest of the paper is organized as follows. In Sec. 2 we introduce a coupled random walk in the particle's momentum and position which provides a good approximation to the full dynamics, and which forms the foundation upon which the rest of our analysis is based. In Sec. 3 we show how this leads to an uncoupled random walk in the particle's energy, which we use to obtain a Fokker-Planck equation for the particle's speed distribution function. This allows us to establish our main result regarding the approach of the momentum distribution to equilibrium. In Sec. 4 we further explore the implications of our analytical result for the final particle temperature, and use the result of that analysis to make predictions for specific dynamical models. We then present numerical results, based on computations of the sort appearing in Fig. 1, that confirm the basic predictions of our analysis. In Sec. 5 we consider the particle's motion in position space, show it to be diffusive, and calculate the dependence of the diffusion constant on the effective temperature T to which the particle equilibrates. We also discuss how the stochastic acceleration of the particle, present when α = 0, is suppressed in a reactive medium, with α = α. Additional numerical results supporting our calculation of the diffusion constant are also presented in that section. Section 6 contains a summary, and a discussion of the connection of our results to those appearing in early work of Chandresekhar [START_REF] Chandrasekhar | Dynamical friction. I. General considerations: the coefficient of dynamical friction[END_REF][START_REF] Chandrasekhar | Dynamical friction. II. The rate of escape of stars from clusters and the evidence for the operation of dynamical friction[END_REF][START_REF] Chandrasekhar | Dynamical friction. III. A more exact theory of the rate of escape of stars from clusters[END_REF]. The Appendix contains details of the perturbative expansions that underly our analysis of Sec. 3.

Particle in a field of scatters: a random walk description

We first recall the random walk description developed in [START_REF] Aguer | Classical motion in force fields with short range correlations[END_REF] to describe a particle moving through an array of random but non-reacting ( α = 0) scatterers, and generalize it to the present situation in which the internal degree of freedom of the scatterer responds dynamically to the presence of the moving particle (α = α). To that end, we consider a typical trajectory generated by the equations of motion (1.1) for a particle that successively encounters scattering centers y n := x Nn at a sequence of instants t n , with incoming momenta p n and impact parameters b n . These quantities are related by (see Fig. 2) Note that we view the impact parameter b n as a vector perpendicular to the incoming direction. We choose t 0 = 0 and x N 0 = 0, with the particle starting near a scatterer located at the origin. The initial displacements and momenta of the scatterers are assumed to be identically and independently distributed according to a stationary state

y - n := y (t n ) = x Nn - 1 2 e n + b n , e n = p n p n b n • e n = 0, b n ≤ 1/2.
ρ(Q, P ) = ρ(H scat (Q, P )) = ρ(Q, -P ) (2.1)
of the scatterer Hamiltonian (1.2). At each scattering event the momentum of the particle undergoes a change ∆p n = p n+1 -p n . Since no forces act on the particle in the region between scattering centers, the momentum p n+1 with which it leaves the n-th scatterer is also the one with which it impinges on the next. The momentum change ∆p n depends on the impact parameter b n and on the arrival time t n at the nth scattering center, through the displacement Q n := Q Nn (t n ) and momentum P n := P Nn (t n ) of that scatterer at the time the particle encounters it. Specifically, we can write

∆p n = R(p n , b n , Q n , P n ), (2.2) where R(p, b, Q, P ) = -α +∞ 0 dτ η(Q(τ ))∇σ(q(τ )) (2.3)
in which q(τ ), Q(τ ) (i.e., without subscripts) denote the unique solution of the single scatterer problem

q(t) = -αη(Q(t))∇σ(q(t)), M Q(t) = -U (Q(t)) -ασ(q(t))η (Q(t)), (2.4)
with initial conditions

q(0) = b - 1 2 p/ p , p(0) = p, Q(0) = Q and P (0) = P. (2.5)
After leaving the support of the nth scatterer, the particle travels a distance n to the next, which it reaches after a time ∆t n = t n+1 -t n = n / p n+1 with impact parameter b n+1 , that both depend on the geometry of the array of scatterers, and on the dynamics of the nth scattering event through the precise point y + n = y n (t + n ) from which the particle leaves the n-th scatterer, and through the outgoing direction e n+1 . The particle finds the next scatterer in the state (Q n+1 , P n+1 ). The process then repeats itself.

Based upon this description of the dynamics, and ignoring the role of recollisions, we argue, as in [START_REF] Aguer | Classical motion in force fields with short range correlations[END_REF], that the motion of an ensemble of particles through an array of scatterers of this type is well approximated by a coupled discrete-time random walk in momentum and position space. Each step of the walk corresponds to one scattering event, where the variables (Q n , P n ) that characterize the scatterer, and the variables n , b n that characterize the approach of the particle onto the scatterer, are drawn from the distributions that govern them in the actual system of interest. In this manner, starting from a given initial condition (q 0 , p 0 ), the velocity p n+1 , the location y n+1 , and the time t n+1 of the particle immediately before scattering event n + 1 are iteratively determined through the relations

p n+1 = p n + R (p n , κ n ) t n+1 = t n + * p n+1 y n+1 = y n + * e n+1    (2.6)
where κ n = (b n , Q n , P n ) and where R (p, κ) is defined by (2.3) and (2.4)-(2.5). In this process, the parameters (Q n , P n ) are assumed to be independently and identically chosen from the assumed stationary distribution ρ (H scat (Q, P )), while the b n are independently chosen at each step uniformly from the d -1 dimensional ball of radius 1/2 perpendicular to p n . Note that the impact parameter b n is the only one of these parameters that is not already independent. We are implicitly assuming therefore that any correlations in the impact parameter die off quickly with successive collisions, and that those portions of a long trajectory that occur at the same speed sample all possible impact parameters in this way. Note also that, without the loss of any essential physics, we have in (2.6) replaced the random variable n at each time step with the average distance * < L * between scattering events. With these assumptions, (2.6) defines a Markovian random walk determined completely by p 0 and the random process (κ n ) n∈N . Note that, for each n, when p n is defined by (2.6), κ n is independent of p n , since the latter only depends on values of κ k with k < n. In what follows we write • for averages over all realizations of the random process κ n . In addition, for a function f depending on p and κ = (b,

Q, P ), b • p = 0, b ≤ 1/2
we denote the average over κ by

f (p) = b ≤1/2 db C d dQ dP ρ(Q, P )f (p, b, Q, P ), (2.7)
where C d is the volume of the ball of radius 1/2 in R d-1 .

Equilibration

In this section, we show that, for suitable initial conditions described below, the distribution of particle speeds for an ensemble of particles evolving according to the random walk (2.6) asymptotically approaches a Maxwellian distribution, and derive expressions for the temperature T to which it equilibrates. To this end, we note that the first equation defining the coupled random walk in (2.6) is independent of the remaining two. This allows us to temporarily ignore the spatial motion of the particle, and the time, and independently focus on the random walk that occurs in the momentum p n and/or the kinetic energy E n = p n 2 /2 of the particle as a function of the step (i.e., collision) number n. Thus, for example, the energy exchanged

∆E (p, κ) = 1 2 (p + R (p, κ)) 2 -p 2
during a single collision between a particle of fixed incoming momentum p and a scatterer, can be written in terms of the solutions to (2.4) as

∆E = α t + 0 Λ(s)σ(q(s))ds where Λ(s) = η (Q(s)) Q(s), (3.1)
and where t + , the instant of time at which the particle emerges from the interaction region ( q ≤ 1/2) at the end of the scattering event, can be replaced in the integral by +∞.

In the following we will suppose that the scatterer distribution ρ has a finite mean energy

E * = H scat (Q, P )ρ(H scat (Q, P ))dQdP
and that the probability of finding a scatterer with an energy much higher than E * is negligibly small. Both the normalized Liouville measure on the energy surface H scat (Q, P ) = E * , associated with the microcanonical distribution, and the Boltzmann distribution Z -1 β exp(-H scat /k B T ), for T large enough, satisfy this condition. In addition, we will assume that the particles are almost always both energetic and fast. By energetic, we mean that the particle has a kinetic energy p n 2 /2 well above the typical interaction potential αη(Q Nn )σ(y(t n ) -x Nn ) that it encounters in any scattering event.

Since η(Q) ∼ Q r , and U (Q) ∼ |Q| r , with r ≤ r, the typical size of η(Q) is of order E r /r * at most, so this first condition can be expressed as

p n 2 >> αE r /r * . (3.2)
In addition, we need the particles to be fast, which means they cross the interaction region in a time p -1 that is short compared to the typical time over which the scatterer evolves. Since the period τ E * of H scat at energy E * behaves as

τ E * ∼ E 1 r -1 2 *
, this second requirement is equivalent to the relation

p n 2 >> E 1-2 r * . (3.3)
The conditions (3.2) and (3.3) can of course be imposed on the initial distribution of particle momenta, given the distribution ρ of scatterer energies. Provided that, in the distribution to which the particle eventually equilibrates, almost all of the particles continue to satisfy these conditions, the following analysis will provide an accurate description of the equilibration process. Since, as we will show, the particle distribution equilibrates to a final temperature T for which p 2 is of the order k B T ∼ P 2 ∼ E * [see (3.24)], both of these conditions do indeed continue to be satisfied during the particle's approach to equilibrium.

Since we will focus on fast particles, we expand the function R (p, κ) in inverse powers of p . This leads to a corresponding expansion

∆E (p, κ) = K =0 β ( ) (κ) p + O p -K-1 , (3.4)
for the energy transferred in a single collision, in which the β ( ) (κ) are explicit scalar functions of the incoming collision parameters. We first remark that, since t + is of order p -1 , it is clear that β (0) = 0. A straightforward computation of β (1) and β (2) , worked out in the Appendix, leads to the following results:

β (1) (κ) =αη (Q) P M L 0 ( b ) β (2) (κ) = α M 2 P 2 η (Q) -M η (Q)U (Q) L 1 ( b ) - α α 2M (η (Q)) 2 L 2 0 ( b ),              (3.
β (1) = 0, Σ 2 1 := β (1) 2 = α 2 M 2 (η P ) 2 L 2 0 , (3.7)
since ρ(Q, P ) = ρ(Q, -P ). Also, since we assumed that ρ is stationary for the free dynamics of the scatterer generated by H scat [see (2.1)], one readily checks that, for any choice of η,

[P 2 η (Q) -M η (Q)U (Q)] = 0. Hence β (2) = - αα 2M (η L 0 ) 2 . (3.8)
Assuming p n is large in the sense of (3.2)-(3.3), and dropping for the moment all higher order terms, we conclude from (3.4) that the energy change undergone by the particle during the nth scattering event is approximately given by:

∆E n := ∆E (p n , κ n ) = β (1) (κ n ) p n + β (2) (κ n ) p n 2 . (3.9)
One recognizes here a dominant fluctuating term (in p n -1 ), independent of the value of α, which is of zero average in view of (3.7), and a smaller subdominant term (in p n -2 ) which fluctuates about a negative non-zero average when α = α. Whereas the random fluctuations of the first term have a tendency to increase the particle's speed without bound, the second term, while weaker, is systematic and has a tendency to reduce it. It is the competition between these two effects when α = α that eventually leads the particle to equilibrate, as we show below. If, on the other hand, α = 0, then, as shown in [START_REF] Aguer | Classical motion in force fields with short range correlations[END_REF], the particle undergoes a stochastic acceleration, with it's speed increasing as p n ∼ n 1/3 . Note, therefore, that for this class of systems, we have been able to explicitly separate the force acting on the particle as the result of its passage through the medium into a frictional and a random part, a notoriously difficult problem of statistical mechanics in general (see [START_REF] Kubo | Statistical physics[END_REF], p.37). It is furthermore clear from the above discussion that the frictional part is due to a back-reaction effect: it results from the change induced in the particle's motion by the change in the medium's motion, which is itself brought about by the passage of the particle. This effect is entirely absent when α = 0.

We now prove that the two competing effects described above balance out so as to drive the particle's momentum distribution precisely to a Maxwell-Boltzmann distribution. It is convenient for analyzing the asymptotic behavior to focus on the random walk associated with a new scaled dynamical variable

ξ n = p n 3 3Σ 1 . (3.10)
For that purpose, we first note, using (3.4), that

p n+1 2 p n 2 = 1 + 4 =1 2β ( ) n p n +2 + O p n -7
, where β ( )

n := β ( ) (κ n ).
From this last expression we obtain the relations

p n+1 p n = 1 + 3 =1 β ( ) n p n +2 + O p n -6 p n+1 -p n = 3 =1 β ( ) n p n +1 + β (4) n -1 2 β (1) n 2 p n 5 + O p n -6 ,
and consequently

∆ p n 3 = p n 2 ∆ p n 1 + p n+1 p n + p n+1 2 p n 2 = 3 =1 3β ( ) n p n -1 + 3 β (4) n + 1 2 β (1) n 2 p n 3 + O p n -4 . (3.11)
In terms of ξ n defined in (3.10) and the parameters

n = β (1) n Σ 1 , µ n = -(3Σ 4 1 ) -1/3 β (2) n , ν n = 3 -2/3 Σ -5/3 1 β (3) n , γ n = 1 3 β (4) n Σ 2 1 + 1 2 2 n ,
equation (3.11) takes the form

ξ n+1 = ξ n + n -µ n ξ -1/3 n + ν n ξ -2/3 n + γ n ξ -1 n + O(ξ -4/3 n ).
(3.12)

From (3.7)-(3.8) one finds n = 0, 2 n = 1, µ n := µ = αα 2(3Σ 4 1 ) 1/3 M (η L 0 ) 2 ≥ 0. (3.13)
The rather involved computation of the higher order coefficients β ν n = 0, and that

γ n := γ = 1 3   β (4) n Σ 2 1 + 1 2   = 1 6 (d -2) + δβ (4) a + δβ (4) b 3Σ 2 1 , (3.14)
where δβ are defined in (A.20) and (A.14) of the Appendix. Using these results (3.12) can be written

ξ n+1 = ξ n + n -µξ -1/3 n + γξ -1 n + O 0 (ξ -1/3 n ) + O(ξ -4/3 n ), (3.15)
where the notation O 0 ξ n -1/3 means the term is O ξ n -1/3 and of zero average. Again focusing on particles that over the course of their evolution spend an overwhelming amount of time at high speeds, we drop the asymptotically small error terms in this last equation to obtain a one-dimensional random walk

ξ n+1 = ξ n + n - µ ξ 1/3 n + γ ξ n . (3.16)
Comparing this to the random walk for p n 2 in (3.9), one recognizes again in the second and third term on the right hand side the effect of the dominant fluctuating part of the force and of its subdominant frictional part. The fact that, in (3.16), the dominant fluctuating term is independent of the random variable ξ n itself constitutes an advantage over (3.9) that will simplify the following analysis. It arises from the fact that ξ n is proportional to the third power of the particle's speed. Again, we stress the all-important strict positivity of µ [see (3.13)] that occurs when α = α and that, as we will see, is needed to assure that the third term on the right hand side of (3.16) acts as a source of dynamical friction capable of balancing the diffusive growth of ξ n generated by the independent and identically distributed terms n of zero mean. The crucial role of the much smaller last term in (3.16) will become evident below. We note that essentially the same equation was obtained in [START_REF] Aguer | Classical motion in force fields with short range correlations[END_REF] for the non-reactive case for which α = 0, and so µ = 0. It was then proven for that case that ξ n ∼ √ n, leading to unbounded growth of the energy of the particle. This stochastic acceleration is completely suppressed when µ > 0, as we will see.

The value of γ in (3.14) is obtained as the result of an involved computation of the coefficient β (4) in the Appendix. In particular it is shown there that, up to corrections that are small provided α is small and that the bath is sufficiently energetic (i.e., that E * is sufficiently high), γ depends on the dimensionality d of the system through the simple relation

γ = 1 6 (d -2), (3.17)
which is independent of the model parameters, the interaction η, and the potential U . We now study the asymptotic behavior of the random walk (3.16) executed by the variable ξ n for the case in which µ > 0. For fast particles relevant to the present analysis, and small coupling α, the dynamical quantity ξ n ∝ p n 3 /α 2 will be much greater than unity, and the changes ∆ξ n that occur during any given collision will be small compared to the value of ξ n itself. Thus, over a large number of collisions during which ξ n does not appreciably change, we may take n to be a quasi-continuous variable, in terms of which the random walk (3.16) takes the form of a stochastic differential equation ξ

(n) = (n) -v (ξ(n)) , (3.18)
where (n) is the random process satisfying

(n) = 0, (n) (n ) = δ(n -n ),
and where

v(ξ) = 3 2 µ ξ 2/3 -γ ln ξ. (3.19)
The corresponding forward Kolmogorov equation for the momentum density ρ(ξ, n) of the particle then reads

∂ n ρ(ξ, n) = ∂ ξ v (ξ)ρ(ξ, n) + 1 2 ∂ ξ ρ(ξ, n) . (3.20)
We point out that (3.18) can alternatively be interpreted as a Langevin-type equation describing the rate of change of the position ξ of an overdamped particle moving in a confining potential v and subject to a fluctuating random force (n). In that interpretation the variable n in (3.18) plays the role of time, and (3.20) is the Fokker-Planck equation for the probability distribution ρ(ξ, n) associated with such a Langevin process. It is then clear that the solutions of (3.20) converge for large n (i.e., after many collisions) to the limiting distribution

ρeq (ξ) = ρeq (ξ 0 ) exp -2 ξ ξ 0 dξ v (ξ ) = N exp [-2v(ξ)] = N ξ 2γ exp(-3µ ξ 2/3 ), (3.21)
where N is a normalization constant. Physically, ρeq (ξ)dξ gives, in equilibrium, the fraction of collisions for which the final particle speed p is such that the quantity ξ = 1 3 p n 3 /Σ 1 has a value lying between ξ and ξ + dξ. Changing variable from ξ back to the speed p , one has, in general,

ρ( p , n) = p 2 Σ -1 1 ρ(ξ, n), which gives the distribution ρeq ( p ) = p 2 Σ 1 ρeq (ξ) (3.22)
associated with the fraction of collisions, in equilibrium, in which the final particle speed is p . The average time that a particle emerging from a collision with speed p remains at that speed after the collision is just the fraction ρ eq ( p )dp of time that a particle spends in equilibrium with a speed lying between p and p + d p , is found from the relation

ρ eq ( p ) ∝ * p ρeq ( p ).
Using (3.17), we thus find that the distribution ρ eq ( p , t) of particle speeds p asymptotically approaches a Maxwellian

ρ eq ( p ) = N 0 p 1+6γ exp(- p 2 2k B T ) = N 0 p d-1 exp(- p 2 2k B T ), (3.23)
with an effective temperature T given by the expression

k B T = - (3Σ 1 ) 2/3 6µ = 1 2 (β (1) ) 2 -β (2) = 2 1 2M η 2 P 2 η 2 , (3.24)
the meaning of which we shall further analyze in Sec. 4. Note the crucial role played by the precise value of γ = 1 6 (d -2) appearing in the last term of (3.16); it produces the correct "density of states" prefactor in (3.23), without which the particle speed distribution P eq ( p ) would not strictly be Maxwellian. To finally show that the distribution of (vector) momenta p approaches the corresponding Maxwellian distribution

ρ eq (p) = 1 (2πk B T ) d 2 exp(-p 2 /2k B T ), (3.25)
it is sufficient to prove that moving particles, as a result of repeated collisions, asymptotically lose any memory of their initial direction, so that the corresponding momentum distribution becomes asymptotically isotropic. Such a demonstration is given in Sec. 5, which contains an analysis of the rate at which random collision-induced deflections turn the particle, and in which that information is then further used to calculate the diffusion constant that characterizes the growth of the particle's mean-squared displacement in equilibrium. Before turning to an analysis of the spatial motion, however, we consider in Sec. 4 some specific consequences of our analysis of the particle's approach to equilibrium, and illustrate some of the predictions of our analysis with appropriate numerical calculations.

Generalized equipartition and numerical simulations

We begin by noting that the expression for the effective temperature T in (3.24) gives rise to a generalized equipartition relation: for i = 1 . . . d

1 2 p 2 i T = 1 2 k B T = P 2 2M η := 1 2M η 2 P 2 η 2 . (4.1)
Here • T stands for an average with respect to the Maxwell distribution (3.25) and • η denotes a weighted average of the kinetic energy of the scatterers in the medium, calculated using the weighting function [η (Q)] 2 ρ (Q, P ) , rather than simply using the stationary scatterer distribution ρ (Q, P ) alone.

Clearly, in the case in which the coupling is linear, meaning that η(Q) = Q, (4.1) reduces to the standard equipartition result

1 2 p 2 i T = 1 2 k B T = P 2 2M , (4.2)
according to which all degrees of freedom of the system have the same mean kinetic energy. Note that this result for linear coupling is completely independent of the stationary scatterer distribution ρ : the medium does not itself need to be in thermal equilibrium for (4.2) to obtain. If the coupling is non-linear, so that η = 1, the situation is more complicated, as we have indicated. Nevertheless, for the special case in which the scatterers are themselves in thermal equilibrium at some temperature T * , the corresponding canonical distribution Z -1 scat exp(-H scat (Q, P )/k B T * ) factors into separate distributions for Q and P , and hence, in this case, so does the average

η 2 P 2 = η 2 P 2 .
Thus, in this way the dependence on η disappears, and one again finds that the particle momentum distribution converges to a Maxwellian with a temperature T = T * given by (4.2), independent of the precise form of the coupling. As a result, standard equipartition holds, meaning that asymptotically in time all degrees of freedom of the system have the same average kinetic energy, given by 1 2 k B T . As an interesting example in which the effective temperature T that is reached does not lead to the standard equipartition relation, consider a situation in which the scatterers are oscillators with U (Q) = 1 2 Q 2 , each drawn from a uniform ensemble with phase space density ρ(Q, P ) ∼ δ(H scat (Q, P )-E * ) of fixed energy E * , and quadratically coupled to the particle with η(Q) = 1 2 Q 2 . A simple computation shows that for this case, for all i = 1 . . . d, In other words, the mean kinetic energy per degree of freedom of the particle equals one half the mean kinetic energy of the scatterer: standard equipartition does therefore not hold in this case.

1 2 p 2 i T = 1 2 k B T = 1 2M Q 2 P 2 Q 2 = E * 4 = 1 2 P 2 2M . ( 4
To demonstrate the general features that emerge from the preceding analysis, we have performed numerical calculations for particle-scatterer systems evolving according to the actual equations of motion (1.1), for the two cases of linear coupling η(Q) = Q and quadratic coupling η(Q) = 1 2 Q 2 discussed above. All numerical calculations were performed with particles moving through a two-dimensional hexagonal array of harmonic scatterers (U = 1 2 Q 2 ), centered at the points x N = N 1 u + N 2 v, where u = (a * , 0), v = (a * /2, √ 3a * /2), and N = (N 1 , N 2 ) ∈ Z 2 . For these calculations the form factor was chosen to be σ(q) = χ (1 -2 q ) where χ(x) is the usual step function, equal to unity for x > 0 and to zero for x ≤ 0, and the lattice parameter a * was set equal to 10/9, ensuring a finite horizon. Note that with this choice of a form factor, as can be seen in the close-up in Fig. 1, the particle undergoes (easily computed) impulsive changes in its momentum at the instants it encounters the edge of an interaction region, but at all other times follows straight line trajectories at fixed speed.

In the numerical results presented in Figs. 3456, for a given set of parameters, each particle in an ensemble of 10 4 trajectories was given the same initial speed, at a random point on the boundary of the scatterer at the origin, with an initial velocity drawn with equal probability from all physically A similar sequence is depicted in Fig. 4 for the same system, but with the initial state of each oscillator drawn from a uniform distribution on the energy surface H scat (Q, P ) = E * = 2. The last panel of that figure confirms the prediction that for linear coupling the final distribution is the same as when the oscillators are thermally distributed. Additional numerical results (not shown here) confirm that this limiting distribution is independent of the strength of the coupling parameter α as predicted by (3.24). Figures 5 and6 show corresponding results for the case of quadratic coupling [η(Q) = Q 2 /2], again with α = 0.1, M = 1, and E * = 2. With quadratic coupling, the preceding analysis predicts different final temperatures according to whether the oscillators are distributed thermally with average energy E * or uniformly on their individual energy surfaces H scat (Q, P ) = E * = 2. The numerical results confirm this prediction. The limiting temperature of the particles that is obtained when the scatterers are uniformly distributed, given in (4.3), is one-half the value that emerges when the scatterers are drawn from a canonical distribution having the same average kinetic energy. The standard equipartition relation does not hold, in this case, since the average kinetic energy associated with the particle's degrees of freedom equals one half that associated with the oscillator.

Diffusion and suppression of stochastic acceleration

We now study the motion in position space of an ensemble of particles in thermal equilibrium at high temperature, each of which undergoes the random walk (2.6) under the same hypotheses on the model as in the previous section, in particular that α = α. Our analysis closely follows the approach of [START_REF] Aguer | Classical motion in force fields with short range correlations[END_REF], where the case α = 0 was treated. At the end of the section we compare the behavior of the mean-squared displacement in these two very different situations. We begin by noting that, since y n+1 = y n + * e n and e n = pn pn , an understanding of the particle's trajectory requires more than just a knowledge of how its speed changes in time; one also needs to understand how the random changes induced in the particle's momentum by the scatterers cause it to turn. To explore this question, we use perturbation theory in (2.4), to expand the function

R (p, κ) = K k=1 α (k) (e, κ) p k + O p -K-1 , e = p p .
(5.1) in inverse powers of p . In particular, one finds that

α (1) (e, κ) = -αη(Q) +∞ -∞ dλ ∇σ(b + λe), so that e • α (1) = 0. (5.2)
Starting from the first equation of (2.6) and (5.1), a simple computation [START_REF] Aguer | Classical motion in force fields with short range correlations[END_REF] yields the result

e n+1 = e n + δ n , (5.3)
where

δ n = α (1) n p n 2 + O( p n -3 ), α (1) n := α (1) (e n , κ n ).
Thus, as suggested in Sec. 3, an initially fast particle of momentum p 0 undergoes small random and independent deflections that cause the unit vectors e n to diffuse isotropically over the unit sphere. To determine the number of collisions required for the diffusing unit vectors to cover the sphere, we compute the value m = m * (p 0 ) for which the correlation function

e m -e 0 2 = m-1 k=0 m-1 k =0 δ k • δ k
becomes of order 1. Under the assumptions of the random walk (2.6), the off-diagonal terms of this last expression clearly have zero average so that e m -e 0 2 = m α

(1) 0 2 p 0 4 , provided that, over n ≤ m collisions, the particle's speed does not appreciably change. A more careful computation, carried out in [START_REF] Aguer | Classical motion in force fields with short range correlations[END_REF] verifies that this is indeed the leading order contribution to the correlation function. Hence, we write

m * (p 0 ) = p 0 4 α (1) 2 .
(5.4)

To justify the assumption that the particle's speed does not change appreciably during the first m * (p 0 ) collisions, recall that, as we have seen, the energy distribution of the particle approaches a Boltzmann distribution with a temperature T given by (4.1); since E * is assumed to be large, the same is true for k B T and hence the bulk of the particles in the limiting equilibrium distribution, which have speeds comparable to the thermal speed v T = √ dk B T , are both fast and energetic, in the sense of (3.2) and (3.3). The reasoning above then applies if we can argue that during m T := m * (v T ) collisions, such particles do not appreciably change their speed. To see this, note that from (3.4), and from that fact that β (1) = 0 = β (3) , we can infer that

∆E(p) = β (2) p 2 + O p -4 = - α α(η P ) 2 L 2 0 2M p 2 + O p -4 . (5.5)
Hence, over a large number n of collisions, an initially fast particle with momentum p will, to leading order, decelerate on average at a rate

d p dn = d p dE ∆E(p) = β (2) p 3 .
On average, therefore, it will take on the order of

n * (p) = p 4 4|β (2) |
(5.6) collisions for the particle to slow down to speeds comparable to the critical speed

v * = max{ αE r /r * , E 1-2 r * } obtained from (3.2)-(3.
3). For particles in equilibrium, which have speeds of the order v T , we find, setting n T := n * (v T ), that

n T m T = α (1) 2 (β (1) ) 2 ∼ η 2 (η ) 2 ∼ Q 2 ∼ T 2/r .
Consequently, at high temperatures the average time it takes a typical particle to slow down is much longer than the one it needs to turn, thus justifying the assumptions underlying the derivation of (5.4). We now conclude. In equilibrium, such particles, which will provide the dominant contribution to the growth of the mean-squared displacement of the ensemble, travel for roughly m T collisions before changing direction significantly. After this many collisions the direction of motion of such a particle will be uncorrelated with its initial direction. The process then repeats itself. We therefore expect the large-scale motion of the particle to be well approximated by the following random walk:

y (k+1)m T = y km T + * m T e km T .
In other words, over large enough length scales, the particle motion is essentially the same as if it traveled on a straight line path for m T collisions, before turning in a random direction and traveling again along a straight line over the same distance. Thus, in this last equation, y km T denotes the position of the particle after a certain number k of these larger excursions. With this picture we can then write

y 2 km T = km 2 T 2 * .
Extrapolating to all values of n, we thus find that y 2 n = nm T 2 * , which implies y 2 t = m T * v T t, from which we obtain, using (5.4), the diffusion constant

D = m T * v T = * v 5 T α (1) 2 .
(5.7)

The diffusive growth in time of the numerically computed mean-squared displacement of an ensemble of particles moving according to the fully Hamiltonian dynamics (1.1), for the model described in Sec. 4, is illustrated in Fig. 7 for the case of linear and quadratic coupling. In the numerical work presented in that figure, straight lines indicate linear fits to the numerical data, which are represented by open symbols, and we have again taken α = 0.1 and M = 1. The different curves appearing in that figure correspond to temperatures such that k B T = 1, 2, 5, 10, 15, 20 and 25. We conclude by providing numerical results that show that this oversimplified argument, which neglects the variation in the speed of the particle

D ∼ T ν , 1 2 ≤ ν = 5 2 - 2r r ≤ 5 2 - 2 r .
(5.8)

As mentioned in the introduction, this shows that diffusion is enhanced by reducing the non-linearity in the coupling (i.e., by lowering r ) and by more strongly confining the scatterer degree of freedom (increasing r). In particular, for quadratic potentials (r = 2), our analysis predicts that D ∼ T 3/2 if the coupling is linear, and that D ∼ √ T if the coupling is quadratic; those predictions are, in fact, as seen in the numerical results presented for these two cases in Fig. 8, which presents the slopes of the linear fits in Fig. 7 as a function of the thermal energy at which they were evaluated.

We conclude this section with a remark on the relation between the inert Lorentz gas (α = 0), treated in [START_REF] Aguer | Classical motion in force fields with short range correlations[END_REF], and the reactive one treated here ( α = α). Note that n * (p) in (5.6) behaves as α-1 . This is as expected, since it was shown in [START_REF] Aguer | Classical motion in force fields with short range correlations[END_REF] that, when α = 0, not only does the particle's speed not decrease, it actually increases without bound as p n ∼ n 1/6 as a function of the collision number n. But, as shown in [START_REF] Aguer | Classical motion in force fields with short range correlations[END_REF], for initially fast particles of momentum p, this stochastic acceleration begins to manifest itself only after n s (p) ∼ p 6 collisions. Thus, the time scale n * (p) ∼ p 4 over which dynamical friction manifests itself by slowing down fast particles, and which characterizes the equilibration process, is always much shorter than is required by the relatively slow process of stochastic acceleration. Thus, in the reactive models considered here, stochastic acceleration is completely suppressed and cannot manifest itself before the system equilibrates. This is a result of the fluctuation-dissipation theorem, which implies that the frictional component of the force cannot be independently made small compared to the strength of its fluctuating part.

Discussion

In this paper we have identified a dynamical mechanism that induces approach to equilibrium in a classical Hamiltonian system. We have demonstrated numerically and analytically that as classical particles scatter repeatedly off the internal degrees of freedom of an array of independent oscillators with initial conditions drawn from a fixed invariant distribution, their momentum distribution approaches a Maxwellian. The effective temperature that characterizes this distribution is obtained from a generalized equipartition relation that involves details of the nature of the coupling between the particle and the scatterers, as well as of the initial invariant distribution of the scatterers. Since our analytical results are derived in a high temperature and weak coupling regime, they strictly speaking only establish the precise behavior of the equilibrium distribution at large values of the momentum.

The numerical data presented in Figs. 3456, on the other hand, clearly indicate that the momentum distribution converges to a Maxwell distribution over the full range of particle speeds, including very low ones. Identification of the dynamical mechanism responsible for this agreement with the Maxwell distribution down to the very lowest speeds is an interesting open problem.

We have also, in our analysis, implicitly focused on the behavior of an ensemble of particles each of which has an entire medium of scatterers to itself, as it were. We have, in particular, not considered effects that might arise in a single system possessing a finite density of moving particles that would, as they equilibrate by exchanging energy with the scatterers, cause the scatterer distribution itself to evolve (presumably) towards a limiting distribution. Indeed, one anticipates that the state of an ensemble of scatterers would, upon being subjected to repeated collisions by a sequence of independent itinerant particles, similarly evolve towards a thermal distribution. The demonstration of such a result would provide an obvious next step towards a more complete dynamical demonstration of approach to equilibrium for the dynamical Lorentz gas, treated as a whole.

The term dynamical friction used in the present work was adopted from a series of papers [START_REF] Chandrasekhar | Dynamical friction. I. General considerations: the coefficient of dynamical friction[END_REF][START_REF] Chandrasekhar | Dynamical friction. II. The rate of escape of stars from clusters and the evidence for the operation of dynamical friction[END_REF][START_REF] Chandrasekhar | Dynamical friction. III. A more exact theory of the rate of escape of stars from clusters[END_REF] by Chandrasekhar, in which he studies the escape rate of a tracer star from a galactic cluster, the latter being treated as a gas of stars. Chandrasekhar first argues, on general grounds, that if there were no dynamical friction, defined as a deceleration of the star along its direction of motion, galactic clusters could not exist for as long as they do; stars would escape from them too frequently as a result of stochastically induced fluctuations in their velocities to values exceeding the escape velocity. Assuming the velocity distribution of the stars to be Maxwellian, he then shows that a friction term is indeed produced as a result of successive two-body scattering events between the tracer star and the other stars of the galactic cluster. Finally, he argues that it has a value compatible with the hypothesis of a Maxwellian velocity distribution.

Our analysis and concerns in the present work, although carried out on very different models, are clearly very similar in spirit to those of Chandrasekhar's. The main differences between our analysis and that of Chandresekhar can be identified as follows. Chandrasekhar computes the average momentum change parallel to the incoming velocity p of the star, per collision, finding ∆p ∼ -C/ p 3 , where C > 0 is a constant that depends on the various parameters of the system. The average momentum change along the direction of motion is therefore negative and of order p -3 , and it is this mechanism that Chandrasekhar identifies with dynamical friction.

In the notation of the current work, this corresponds to a calculation of the quantity e • R(p), and indeed for the current models we obtain a similar result: e

• R(p) = e • α (3) p -3 + . . . with e • α (3) = β (2) - 1 2 α (1) • α (1) = - αα 2M (η ) 2 L 2 0 - 1 2 α 2 η 2 ∇L 0 2 < 0. (6.1)
We argue however that strict negativity of e • α (3) does not, by itself, allow one to conclude that the system can maintain, let alone approach, a Maxwellian momentum distribution. Indeed, as we saw, for that purpose one needs to establish a stronger condition: it is the coefficient β (2) , that in our analysis clearly governs the effect of dynamical friction on the particle's speed (or energy), which must be strictly negative.

To emphasize this point we note that, in our models, even when the Lorentz gas is inert so that α = 0, the value of e • α (3) is strictly negative, as it is in Chandresekhar's analysis. But in the inert Lorentz gas there is no dynamical friction, because β (2) = 0 (see [START_REF] Aguer | Classical motion in force fields with short range correlations[END_REF]); as a result, particles in the ensemble do not equilibrate, they undergo an unbounded stochastic acceleration of the sort that Chandresekhar sought to avoid. In fact, in our models, and probably in Chandresekhar's as well, a negative contribution to the average change in the component of momentum along the direction of motion arises, at least in part, because particles turn (as we have discussed in Sec. 5), an effect that has nothing to do with friction or the dissipation of excess energy. and similarly for other functions of s. Here, the notation q (k) (s) designates any function of s such that q(s) -q (k) (s) = O(s k+1 ); it is, in particular, not unique, an observation that will be helpful in the computations below. Also, since we always have s ≤ t + , which is of order p -1 , it follows that any term of O(s k+1 ) is automatically of O( p k+1 ). Then, with q = b -1 2 e, we have

q (3) (s) = q + ps -αη(Q)∇K 0 (s p , b ) p -2 -α P M η (Q)∇K 1 (s p , b ) p -3 ,
where, for b • e = 0, we define, [compare (3.6)],

L k (µ, b ) = µ 0 dλ σ(b + (λ - 1 2 )e) λ k , K k (µ, b ) = µ 0 dλ L k (λ, b ).        (A.1)
Hence ∆E defined in (3.1) can be written as follows:

∆E = ∆E a + ∆E b + ∆E c + O( p -5 ) (A.2) where, with Λ(s) = η (Q(s)) Q(s), ∆E a = α t + 0 ds Λ(s)σ(q + ps) ∆E b = - α 2 η(Q) p 2 t + 0 ds Λ(s)∇K 0 (s p , b ) • ∇σ(q + ps) ∆E c = - α 2 P η (Q) M p 3 t + 0 ds Λ(s)∇K 1 (s p , b ) • ∇σ(q + ps).
We need to expand each of these terms to order p -4 included in order to determine β ( ) , = 0 . . . 4. In fact, as we have noted in Sec. 3, β (0) = 0 because t + is of order p -1 . Furthermore, since 0 ≤ s ≤ t + and t + is of the order p -1 , to obtain the contributions β

( ) a , β ( ) b , β ( )
c to β ( ) from ∆E a , ∆E b , and ∆E c requires that we expand Λ(s), respectively, to order 3, 1, and 0.

We first determine β (1)) = β

(1)

a , β (2) = β (2)
a , to which ∆E b and ∆E c do not contribute. To expand ∆E a , we compute, from the equations of motion for Q

(t), Q(3) (s) = P M + S(s), where
S(s) = - 1 M s 0 dτ U (Q (2) (τ )) + αη (Q (2) (τ ))σ(q (2) (τ )) , (A.3) potential U . As a result, Q(t)
is periodic in time, and so, therefore, is the interaction potential αη(Q(t))σ(q). The corresponding results of [START_REF] Aguer | Classical motion in force fields with short range correlations[END_REF], which were derived under conditions which include this case, can then be applied to this situation. In particular, the coefficients β (i) (α = 0), for i = 0 . . . 4 were computed in [START_REF] Aguer | Classical motion in force fields with short range correlations[END_REF] and found to satisfy the equality2

β (4) (α = 0) = 1 2 (d -3)(β (1) ) 2 (α = 0). (A.11)
Using (3.5) this yields

β (4) (α = 0) = 1 2 (d -3)(β (1) ) 2 . (A.12) Hence β (4) = β (4) ( α = 0) + δβ (4) = 1 2 (d -3)(β (1) ) 2 + δβ (4) a + δβ (4) b (A.13)
where we have used the result δβ (4) c = 0 which follows from (A.7). Moreover, δβ

b can be read off from (A.8), i.e., δβ

b = - α 2 αη(Q)η (Q) 2 M +∞ 0 dλ L 0 (λ, b )∇K 0 (λ, b ) • ∇σ(q + λe). (A.14) (4)
It is therefore sufficient to determine δβ

a , which turns out to be less difficult than to determine β (4) a directly. Indeed, and this is our second observation, all expansion coefficients that we compute, and in particular β (4) a , are polynomial in α and α. Thus, to compute δβ (4) a it suffices to compute only those terms that are of non-vanishing order in α. The required terms turn out to be of first or second order at the most. Taken together, these two observations allow us to compute far fewer terms, and result in a considerable simplification that is ultimately brought about by the fact that (A.12) is neither obvious nor easily established. To prove it by direct expansion in the present context would, in fact, necessitate combining terms from β To see this, it suffices to examine the various error terms and use the fact that |Q| ν ∼ E ν/r * and P 2 ∼ E * , for all ν ∈ N. In the case of linear coupling η(Q) = Q and an even confining potential U (-Q) = U (Q), the computation considerably simplifies since one then readily sees that δβ << E * ∼ P 2 so that indeed the last term is negligible for large E * . More explicitly still, for the model studied numerically in this paper, where d = 2, U (Q) = 1 2 Q 2 and σ(q) = χ(1 -2 q ), one finds after an explicit computation of the last term above that for this particular case differing by p 3π 32k B T from its Maxwell-Boltzmann value. Even at the relatively small value of the thermal energy (k B T = 2) used in the numerical results of Fig. 3, a deviation in the tails of the distribution of this magnitude in the power law multiplying the Boltzmann factor would be difficult to observe without considerably better statistics.

  Figure1: The left panel shows an extended part of a typical trajectory, starting from the origin, computed for a particle moving through a hexagonal lattice of harmonic scatterers, with unit lattice spacing and linear coupling, as described more fully in Sec. 4. The right panel depicts a close-up view of the region indicated by an arrow in the left panel. In the close-up, the points where the particle enters and exits each circular interaction region are indicated with black dots.

Figure 2 :

 2 Figure 2: A particle at time t n impinging with momentum p n and impact parameter b n on the nth scatterer, centered at the point y n = x Nn .

  5) where for k ∈ N, and b, e ∈ R d , with e • e = 1, b • e = 0, L k ( b ) := 1 0 dλ λ k σ(b + (λ -that the right hand side in (3.6) is a function of b only, as a result of the rotational invariance of σ. It follows that [see (2.1) and (2.7)]

n

  is performed in the Appendix and leads to the result that [see (A.10) and (A.13)]

Figure 3 :

 3 Figure3: Approach to equilibrium of the particle speed distribution for a particle moving through a lattice of thermally equilibrated harmonic oscillators with average energy E * = 2, with linear coupling. The solid curve in the last panel is a Maxwellian velocity distribution with the predicted final thermal energy k B T = 2, equal to that of the oscillators.

Figure 4 :

 4 Figure4: Approach to equilibrium of the particle speed distribution for a particle moving through a lattice of harmonic oscillators each drawn from their own microcanonical distribution with average energy E * = 2, and linear coupling. The solid curve in the last panel is a Maxwellian distribution with the predicted final thermal energy k B T = 2.

Figure 5 :

 5 Figure5: Approach to equilibrium of the particle speed distribution for a particle moving through a lattice of thermally distributed harmonic oscillators with average energy E * = 2, with quadratic coupling. The solid curve in the last panel is a Maxwellian velocity distribution with a thermal energy k B T = 2, equal to that of the oscillators.

Figure 6 :

 6 Figure6: Approach to equilibrium of the particle speed distribution for a particle moving through a lattice of harmonic oscillators each drawn from their own microcanonical distribution with average energy E * = 2, and quadratic coupling. The solid curve in the last panel is a Maxwellian velocity distribution with a thermal energy k B T = 1 predicted by analysis, one-half the value found in Figs.3-5.

Figure 7 :

 7 Figure7: Numerically computed mean squared displacements for an ensemble of particles moving through a lattice of thermally equilibrated harmonic oscillators, for linear and quadratic couplings as indicated, at thermal energies k B T = 1, 2, 5, 10, 15, 20, 25. Steeper slopes in each panel correspond to higher temperatures.

Figure 8 :

 8 Figure 8: Temperature dependence of the diffusion constant D obtained from the data in Fig. 7 for the case of linear coupling (r = 1) and quadratic coupling (r = 2).

c=dλ K 1

 1 (α = 0). A different, more efficient approach is used in[START_REF] Aguer | Classical motion in force fields with short range correlations[END_REF].To compute δβ (4) a , we start by rewriting ∆E a in the form ∆E a = αA( p , κ, α, α)P M + αB( p , κ, α, α) + O( p -5 B 1 ( p , κ, α, α) + B 2 ( p , κ, α, α) + O( p -5 ), (A.17) with B 1 ( p , κ, α, + λe) (A.18) B 2 ( p , κ, α, α) = α α M p 4 (η(η ) 2 )(Q) λ , b ) • ∇σ(q + λ e). (A.19)Then, using (A.15), we find that 3 δβ (4) a = δ∆E (4) a = αδA (4) (κ, α, α)P M + αδB (4) (κ, α, α). (A.20)We first computeδA( p , κ, α, α) = A( p , κ, α, α) -A( p , κ, α, (λ, b )σ(q + λe) -α M 2 p 4 P η (Q) 2 1 0 dλ K 1 (λ, b )σ(q + λe) + O( p -5 ). (A.21)One then easily finds thatδA (4) = -α M 2 P η (Q)η (Q) 1 0 dλ λK 0 (λ, b )σ(q + λe) (λ, b )σ(q + λe) . (A.22)Next, we compute the contributions from δB in (A.15). The term coming from B 2 can be read off immediately from (A.19): λ , b ) • ∇σ(q + λ e).(A.23)For δB 1 , we writeδB 1 ( p , κ, α, α) = 1 p 1 0 dλ I(λ)σ(q + eλ),whereI(λ) = δ η (Q (3) ) × S 1 ( 3) (λ) + Ĵ(3) (λ) σ(q + λe). (A.24)It follows from (A.4) thatJ(λ) = -α M p 2 η Q (3) λ p , α = 0 × η (Q)K 0 (λ, b ) + P M p η (Q)K 1 (λ, b ) S 1 λ p +O( p -5 ) = -α M p 2 η (Q)η (Q)K 0 (λ, b )S 1 (Q) + αη (Q)L 0 (λ, b ) + O( p -2 ),this finally yieldsJ (3) (λ) = α M 2 η (Q)η (Q)K 0 (λ, b ) λU (Q) + αη (Q)L 0 (λ, b ) . (A.25)To compute Ĵ(3) (λ), we note that, since S 1 (λ/ p ) is of order p -1 , we haveĴ(λ) = η Q (2) Q (2) (τ )) -U (Q (2) (τ, α = 0) +αη (Q (2) (τ ))σ(q + pτ ) = α M p 2 η (Q) λ p 0 dτ U (Q (2) (τ, α = 0))K 0 (τ p , b ) Q (2) (τ ))σ(q + pτ ) + O( p -4 ) = α M p 3 (η U )(Q) 3 (U η )(Q)L 2 (λ, b ) α2 M 2 p 3 (η η )(Q) λ 0 dλ K 0 (λ , b )σ(q + λ e) -α M 3 p 3 P 2 η (Q)L 2 (λ, b ) + O( p -4 ) (U η )(Q) -P 2 η (Q) (λ) = η (Q)δS

  (U η )(Q) -P 2 η (Q) λ 2 δS (1) 1 . (A.26)One can now compute β (4) in (A.13) by using (A.14), (A.20), (A.22), (A.24), (A.23), (A.25) and (A.26). Obviously, a general condensed expression for β (4) and hence of γ, defined in (3.14), is not readily obtainable. However, one can easily check that, under the hypotheses of this paper, namely α = α small, E * large and|η(Q)| ∼ |Q| r , U (Q) ∼ |Q| r , with 0 < r ≤ r, γ = 1 6 (d -2) + (α, E -1 * ).

  δA (4) = δB

  3) β (1) 2 + αδB

  Since U ∼ |Q| r-2 , it follows that U ∼ E 1-2 r *

All variables and constants appearing in this equation of motion can be thought of as

* / p . Thus, the equilibrium probability density ρ eq ( p ), which governs

Here and in what follows, whenever a function f depends on α, α, and possibly on other variables such as p and κ, we shall write f ( α = 0) for the values of the function f on the hyperplane α = 0, and δf := f -f ( α = 0).

Here and in what follows, we write f (k) for the coefficient of p -k in the expansion of a function f in powers of p -1

Acknowledgments This work was supported by Ministry of Higher Education and Research, Nord-Pas de Calais Regional Council and FEDER through the Contrat de Projets Etat Region (CPER) 2007-2013. P. E. P. thanks the Research Centre INRIA-Futurs for its hospitality and support for his stay at SIMPAF, where part of this work was performed.

A The high momentum expansion of the energy transfer

In this Appendix we determine the expansion coefficients β ( ) in (3.4), for = 0 . . . 4, as well as their means and variances, which are crucial ingredients of the main results derived in the bulk of the paper. We obtain these from (3.1) and a perturbative calculation of the solutions to (2.4) to sufficiently high order. To this end, it is convenient to introduce a notation in which, for any k ≥ 0 q(s) = q (k) (s) + O(s k+1 ), and hence

From this, one finds, with L 0 as defined in (A.1),

Now using the result

and

one easily obtains (3.5).

We now turn to

the calculation of which is more involved. Fortunately, in order to establish (3.16) only the averages of these quantities, in the sense of (2.7), are needed; these are much easier to compute. As a first step in this calculation we note that, since to lowest order

Next, using (A.5), we easily determine the leading terms of

29

where

Now, to compute β

(3) a we need to identify the s 2 contribution to Λ(s). For that purpose, we first write (A.3) as

where

Because S 1 (s) is of order s and S 2 (s) is of order s 3 , it follows that

Identifying the s 2 term in this expression is now easily done using the relevant terms of (A.4); one observes that they are linear or cubic in P and hence of zero average in any stationary distribution ρ. In conclusion, since β

a . To simplify the computation, we make the following two observations. First, when α = 0 in the equations of motion (1.1), the particle moves in a time-dependent potential αη(Q(t))σ(q) where Q(t) solves the one-dimensional equations of motion of a particle of mass M moving under the influence of the confining