Equilibration and Thermalization of Classical Systems

It is demonstrated that the canonical distribution for a subsystem of a closed system follows directly from the solution of the time-reversible Newtonian equation of motion in which the total energy is strictly conserved. It is shown that this conclusion holds for both integrable or nonintegrable systems even though the whole system may contain as little as a few thousand particles. In other words, we demonstrate that the canonical distribution holds for subsystems of experimentally relevant sizes and observation times.


Introduction
Boltzmann's postulate of equal a priori probability is the cornerstone of statistical mechanics. This postulate eliminates the difficulties of deriving the statistical properties of a many-body system from its dynamical evolution by introducing a probabilistic description in terms of the (micro)canonical ensemble. In classical mechanics, the equivalence of these two fundamentally different levels of description remains elusive [1,2].
Considerable progress has been made through the development of ergodic theory, but more than a hundred years after its conception, the difficulty of proving that dynamical systems are mixing constitutes a considerable obstacle for explicit demonstrations that the canonical distribution follows directly from (nondissipative) Newtonian dynamics of systems containing a finite number of particles [1,2].
In connection with the mixing property, the discovery of 'deterministic chaos' has made a huge impact on old discussions about the relations between classical and statistical mechanics [3][4][5][6]. The irreversibility in the macroworld is often related with unstable motion in generic mechanical systems characterized by positive Lyapunov exponents and a nonzero Kolmogorov entropy. Integrable systems are considered from this point of view as exceptional (not to say pathological) and irrelevant for the problem of justification of statistical physics. Indeed, there is no tendency to equilibration in a system of noninteracting entities (particles of an ideal gas, normal modes in harmonic oscillator systems, solitons in Korteweg-de Vries (KdV) systems, etc). Moreover, the famous Fermi-Pasta-Ulam paradox [7][8][9][10][11] and its interpretation in terms of closeness of their model to the completely integrable KdV model (see [4]) has emphasized (probably, overemphasized) this point.
An alternative route to thermalization is to couple the dynamical system to a fictitious thermostat and to analyze the relaxation of the system + thermostat to a stationary state [12]. This is standard practice in molecular dynamics simulations [13]. Such thermostats come in different flavors, may have different dynamics and bring the system into the thermal equilibrium state or into another stationary state. The thermostated methods can be analyzed to obtain the relaxation to equilibrium via a fluctuation relation involving dissipation [12,14], basically an extension of the Boltzmann H-theorem. Although these methods are powerful, due to the introduction of dissipation this approach has no bearing on the original problem, namely to show that reversible Newtonian dynamics alone leads to canonical distributions.
In the foregoing discussion, the focus was on the dynamic properties of the whole closed system, whereas in reality, the system of interest is always a subsystem of a larger system. Therefore, in this paper, we study the equilibration and thermalization of a subsystem S, containing N s particles, of an isolated (closed) system S + E with N = N S + N E particles. All particles not existing in S are considered to be in the environment E, containing N E particles. The system S + E evolves under time-reversible Newtonian dynamics.
For the dynamical properties, integrability is an important issue. The integrability of the isolated Hamiltonian system implies that its Liouville operator L S+E is diagonal in the representation of angle-action variables [5]. If we choose a subsystem S of the isolated integrable system which is also integrable, then its Liouville operator L S is diagonal in the angle-action variables of the subsystem. However, these variables can be different from those of the isolated Hamiltonian system. L S+E and L S do not commute and cannot be diagonalized simultaneously. In this situation, the subsystem S may be expected to equilibrate but surprisingly this very natural issue has not been clarified so far and thus requires additional studies.
From the probabilistic formalism of statistical mechanics, see for example [1,2], it is quite natural to expect that under certain weak assumptions (e.g. a large number of degrees of freedom, weak interaction between S and E and the total energy fixed) the configurations of the subsystem S are distributed according to the canonical ensemble for S [1,2]. However, from the dynamical point of view, the relevant question is whether a single time-dependent trajectory in the phase space of S will produce configurations that are in concert with the canonical distribution of statistical mechanics. To the best of our knowledge this question has not been answered as yet.
In this paper, we present clear and unambiguous evidence that the canonical distribution for S follows directly from the solution of the time-reversible Newtonian equation of motion of the entire system. Evidently, the time for S to reach thermal equilibrium may depend strongly on the initial configuration. For some initial configurations, equilibration of S may be extremely slow or even not occur at all. However, such initial configurations are exceptional and may, from a physical point of view, be regarded as pathological. In this paper, we only consider 'typical' initial configurations for S and E which are abundant and show that the energy of the subsystem and of the environment equilibrate on a relatively short, microscopic time scale, even if the whole system contains only of the order of a few thousands of degrees of freedom. The key feature of our demonstration is that we follow the time evolution of a closed (isolated) system S + E with a fixed energy and consider only a subsystem S. We do not perform ensemble averaging, nor do we invoke arguments based on (non)ergodicity [1,2,[15][16][17] or on the thermodynamic limit [1,2]. Nevertheless, we observe that S is governed by the canonical distribution.
Our demonstration is not a mathematical proof but is based on the exact numerical solution of what is perhaps the simplest of all interacting many-body systems: a one-dimensional harmonic oscillator model of a solid, see figure 1(a). By solving the Newtonian equation of motion of the whole, integrable system S + E and analyzing only a single trajectory of the subsystem in phase space, we show that the number of times that the subsystem S is observed to possess a certain energy is distributed according to canonical ensemble theory, implying for example that within the subsystem equipartition holds [18][19][20][21]. Repeating the analysis for a one-dimensional model of classical magnetic moments, see figure 1(b), which is known to exhibit chaotic behavior [22], it is found that this conclusion does not depend on whether or not the motion is chaotic. For both models, the distributions extracted from the single-trajectory Newtonian dynamics are found to be in excellent agreement with the corresponding results of microcanonical Monte Carlo simulations.
In a landmark paper, Bogoliubov [23] showed that, under two assumptions of which one is insignificant [24], a single oscillator coupled linearly to a harmonic environment described by the canonical distribution relaxes to the canonical distribution. Our results suggest that this conclusion generalizes to subsystems containing more than one oscillator and to (large) environments that initially are not in thermal equilibrium. Moreover, it is not necessary to consider ensembles. One particular trajectory of the whole system suffices to yield the canonical distribution. Livi et al [25] numerically studied the chaotic behavior of two different nonlinear models and computed the time-averaged, ensemble-averaged potential energy and specific heat of subsystems with up to 32 particles out of 1024 particles. With some exceptions, good agreement with the statistical mechanical results was found. In contrast to Livi et al we do not invoke ensemble averaging but extract all information from a single trajectory.
The paper is structured as follows. Section 2 introduces the harmonic lattice model ( figure 1(a)), the computational techniques used to calculate the time-dependent trajectory in configuration phase space and the microcanonical Monte Carlo method which we employ to provide a final check on the distributions extracted from the time-dependent solution. Section 3 follows the same structure as section 2, discussing a nonintegrable model of classical magnetic moments ( figure 1(b)). The simulation procedure is specified in section 4. In section 5, we recall some well-known results from statistical mechanics and explain how we test the hypothesis that the distributions extracted from a single, time-dependent trajectory agree with the canonical distribution. Our simulation results are presented in section 6 and in the appendix. Finally, section 7 gives some estimates for the time scales of equilibration and contains our conclusion.

Harmonic oscillators
We consider what is presumably the simplest model of interacting particles, namely the harmonic model for vibrations of atoms in a solid. The coordinate of the jth atom is written as r j = R j + x j , where R j denotes the equilibrium position of the atom and x j represents the displacement of the atom from its equilibrium position which is assumed to be much smaller than the interatomic distance. For the present purpose it follows from our numerical work (results not shown) that the dimensionality of the underlying lattice and boundary conditions are not important. Therefore, we confine the discussion to atoms arranged on a ring, see figure 1(a).
The Hamiltonian of the system reads where x j and p j are the displacement and momentum of the jth particle, m 2 is the spring constant and N is the total number of particles. For convenience, m and are set to 1 and all quantities such as momenta and displacements are taken to be dimensionless. Because of the periodic boundary conditions p j = p j+N and x j = x j+N , equation (1) can be brought in diagonal form by a Fourier transformation. Changing variables for k = 0, . . . , N − 1, the Hamiltonian reads where ω k = 2|sin π k/N |, P k = P N −k and X k = X N −k . For k = 0, . . . , N − 1, the variables change in time according to where P k (0) and X k (0) are the initial values of the momentum and displacement in Fourier space, respectively. Procedures to choose the initial values are discussed in section 2.3. The motion of each set of coordinates {X k , P k } is described by a single sinusoidal oscillation and is decoupled from the motion of all other sets. Equation (4) directly shows that the model equation (1) is integrable.

Newtonian time evolution
Instead of numerically solving the equations of motion for the real-space coordinates for a given initial configuration {x 0 (0), . . . , x N −1 (0), p 0 (0), . . . , p N −1 (0)}, we perform two fast Fourier transformations to compute {X 0 (0), . . . , X N −1 (0)} and {P 0 (0), . . . , P N −1 (0)} and then use equation (4) to find {X 0 (t), . . . , X N −1 (t), P 0 (t), . . . , P N −1 (t)}. Application of the inverse fast Fourier transformation then yields {x 0 (t), . . . , The number of arithmetic operations is of the order 4N log 2 N ; hence the procedure is very efficient. Most importantly, it does not suffer from systematic or cumulative numerical errors. The accuracy of the solution is close to machine precision (14 digits), independent of the time t. Therefore, the energy of the whole system is conserved for all times, up to machine precision.

Microcanonical Monte Carlo simulation
In order to compare with the Newtonian time evolution of section 2.1, we construct a microcanonical Monte Carlo procedure for the oscillator model. We want to generate statistically independent sets of variables {x 0 , . . . , x N −1 , p 0 , . . . , p N −1 } such that the total energy equation (1) is fixed to a chosen value E 0 . This is most conveniently done by introducing new variables Z k = ω k X k and generating Gaussian random variables Note that we exclude the k = 0 variables because a nonzero P 0 merely causes an insignificant uniform translation of the particles on the whole ring and a nonzero X 0 does not contribute to the energy because ω k=0 = 0. The next step is to introduce new variablesP such that This procedure generates uniformly distributed points of the hypersphere with radius E 0 [26]. Applying the inverse fast Fourier transformation to the data sets {0,Z 1 /ω 1 , . . . ,Z N −1 /ω N −1 } and {0,P 1 , . . . ,P N −1 } yields the desired real-space variables.

Initial configuration
It is trivial to construct initial configurations for which the subsystem will not show relaxation to some stationary state at all. For instance, if only one of the normal modes of the oscillator system is excited (meaning that the energy of all other modes is zero), the subsystem energy does not change with time. As mentioned in the introduction, we disregard such initial configurations simply because they are exceptional and physically irrelevant from the viewpoint of statistical mechanics. In mathematical terms, the exceptional initial conditions usually form a set of measure zero of all initial conditions.
We have used various procedures to construct the 'typical' initial configurations of the subsystem and the environment: is the specified initial energy (not the kinetic energy!) per particle of the environment. Initially, all the energy in the environment is kinetic, the potential energy being zero. Similarly, but independently, we generate the momenta { p 0 , . . . , p N S −1 } (setting x 0 = · · · = x N S −1 = 0) of the particles in the subsystem. Then, we remove the zero-energy translational mode of the whole system. This procedure is flexible in that it allows us to prepare the environment and subsystem at different energies. A simple calculation shows that the Kullback-Leibler distance [27] between a probability distribution generated by this procedure and a canonical distribution (with the corresponding average energy) is infinite, suggesting that the initial configurations have little in common with those drawn from the thermal equilibrium distribution. To demonstrate equilibration, the initial energies of the environment E E /N E and subsystem E S /N S need to be different. 2. We employ the same procedure as in procedure 1, but instead of drawing from a Gaussian distribution we generate the momenta by drawing from a uniform distribution over the interval [−a, a] where a = √ 6E E /N E and √ 6E S /N S for the environment and subsystem, respectively. The procedure is used to check that the precise form of the distribution is not important. 3. We use the same procedure as in procedure 2 with a = √ 3E E /N E but generate both the momenta and coordinates. This method can generate configurations that show equipartition without being drawn from the canonical distribution and provides yet another check. 4. We generate random values for the normal-mode variables {X 0 (0), . . . , X N −1 (0), P 0 (0), . . . , P N −1 (0)} using the canonical distribution with temperature T and set P 0 (0) = 0 to remove the zero-energy translational mode. This procedure does not allow independent control of the subsystem energy and is only useful for studying the dynamics in equilibrium. We use it for comparison with the data obtained by methods 1 and 2.

Magnetic moments
As a nonlinear model that exhibits chaotic motion, we consider a system of classical magnetic moments interacting with nearest neighbors only, see figure 1(b). The Hamiltonian is given by where S i = (sin θ i cos φ i , sin θ i sin φ i , cos θ i ) is a three-dimensional unit vector representing the angular momentum of the classical magnetic moment at lattice site i, J is the interaction strength and N is the total number of spins. This model may be viewed as the classical limit of the one-dimensional spin-1/2 Heisenberg model [28]. Our motivation to consider the nonlinear model equation (7) instead of say model equation (1) extended by adding anharmonic terms is that, as discussed below, the time evolution of the former can be computed numerically with strict conservation of the energy of the entire system, whereas for the latter, we do not know how to construct time integrators that preserve the energy of the entire system. Strict conservation of the energy of the whole system is an essential requirement for any time integrator that aims to keep the whole system in the microcanonical state throughout the time evolution.
The anharmonic character of model equation (7) can be made explicit by introducing the canonically conjugate variables P i ≡ cos θ i and Q i ≡ φ i , yielding For numerical purposes, it is much more convenient to work with the variables S i than with P i and Q i . In the canonical ensemble, for a chain of N + 1 spins, the internal energy is given by [28] where β = 1/k B T and is the Langevin function. The specific heat reads [28] In our numerical work, we adopt periodic boundary conditions S i = S i+N and then, in a strict sense, equations (9) and (11) should be replaced by more complicated expressions in terms of spherical Bessel functions [29] but the finite-size corrections vanish exponentially with N and are, for our purposes, negligible. In the following, we use units such that k B = 1 and J = 1.

Newtonian time evolution
In general, the equations of motion d dt of the unit vectors {S i } are nonlinear. Only in the low-temperature limit, the Hamiltonian equation (7) can be approximated by a harmonic oscillator model of the form equation (1). and (a, b, c) form a right-handed set of orthogonal unit vectors [30,31]. More generally, equation (12) has simple analytical solutions for N = 2 and 3 [22]. The motion of N = 4 magnetic moments arranged on a ring is regular [22]. For N > 4, the system exhibits chaotic motion [22], except for special initial conditions such as the spin-wave and soliton solutions [32].
We integrate the nonlinear equations of motion, equation (12), using a second-or fourthorder Suzuki-Trotter product-formula method with a time step τ = 0.02 which conserves (i) the volume of the phase space, (ii) the length of each magnetic moment and (iii) the total energy [33]. Due to the chaotic character of the motion of the magnetic moments, their trajectories are unstable with respect to rounding and time-integration errors. Nevertheless, the numerical method used guarantees that the motion of the magnetic moments strictly conserves the energy, as required for a reversible Newtonian dynamic or a microcanonical ensemble simulation.

Microcanonical Monte Carlo simulation
As an alternative to the numerical integration of the classical equations of motion (section 3.1), we also perform microcanonical Monte Carlo simulations. Given a configuration of magnetic moments with energy E 0 , we can sample the hyperspace of configurations as follows. We randomly select a particular magnetic moment, say j, and with probability 1/2 choose to perform a random rotation of S j about S j−1 + S j+1 or to carry out a reflection of the vector S j about the vector S j−1 + S j+1 in the plane formed by these two vectors. This process is repeated many times. Obviously, neither the rotations nor the reflections alter the local energy E j ; hence the total energy is strictly conserved. Nevertheless, this microcanonical procedure is ergodic.

Initial configuration
The initial configurations of the magnetic moments in the environment and subsystem are generated using techniques similar to those used for the oscillator system: 1. We set all θ i = π/2 and employ the Metropolis Monte Carlo method [34] for generating a configuration of the φ i of the environment according to the distribution where A is a normalization constant. Similarly but independently, we generate a configuration of φ's of the subsystem using the distribution where A is another normalization constant. This procedure is flexible in that it allows us to prepare the environment and subsystem at different energies (controlled by b E and b S ). The Kullback-Leibler distance [27] between a probability distribution generated by this procedure and a canonical distribution (with the corresponding average energy) with Hamiltonian equation (7) is infinite, suggesting that the initial configurations have little in common with those drawn from the thermal equilibrium distribution. 2. The initial configurations for the environment and subsystem are drawn from canonical distributions with Hamiltonian equation (7) at different temperatures. This method is used to test the dependence of the results on the initial configuration only.

Simulation procedure
The whole system is divided into two parts, a subsystem S with N S particles and an environment with for the oscillator and magnetic system, respectively. In our numerical work, the initial values of the coordinates or magnetic moments are generated according to the procedures specified earlier.
We are primarily interested in the motion of the particles contained in the subsystem defined by the real-space indices j ∈ {0, . . . , N S − 1}, 1 N S < N . The whole system, consisting of the subsystem and the environment, evolves in time according to the Newtonian equations of motion, thereby conserving the energy of the whole system. During the time evolution, we compute the subsystem energy E S and count how many times this energy falls within a certain energy range, yielding the histogram P(E s ) of subsystem energies. The timedependent data are sampled with a time interval t = 1 (for the definition of the units used, see below), the bin size E s = 0.001N S and the number of samples is 10 6 -10 8 . The precise choices for t and E s are not important, and the number of samples affects only the quality of the histogram. The average and variance of E S are calculated as The histograms P(E S ) are constructed in the same way for both models. Guided by the expectation that in the thermodynamic limit the time-dependent fluctuations of the subsystem energy are related to the specific heat as defined by statistical mechanics, we define where T S denotes the subsystem temperature which, up to this point, is yet to be determined. Although statistical mechanics teaches that in the thermodynamic limit, for infinitely long times and weak coupling, T S is expected to agree with the temperature T as defined in statistical mechanics [1,2], there is, to our knowledge, no proof that this equivalence holds for experimentally relevant system sizes and observation times. In section 5, we show that the parameter T S can be extracted from the distribution P(E S ) and that for 20 N S , T S ≈ T to a very good approximation.

Statistical mechanics
The hypothesis that Newtonian dynamics causes S to visit points in phase space with frequencies that match with those of the canonical probability distribution is tested as follows.
According to statistical mechanics, the distribution of energy in the canonical ensemble is given by [2] p where T , E, g(E), S(E) = ln g(E) and Z are the temperature (in units of k B = 1), the total energy, the density of states, the (microcanonical) entropy and the partition function, respectively. The function p(E) has a maximum at some energy E * , the most probable energy at the temperature We remark that even though the right-hand side of equation (19) denotes a microcanonical temperature, our work uses a unique canonical temperature T . The proper choice of E = E * in equation (19) removes any ambiguity. In the vicinity of E * , we have [2] p where a n = 1 n!
and A is a normalization constant. In the thermodynamic limit (N E → ∞ before N S → ∞), all but the quadratic term in the exponential can be neglected and the distribution is Gaussian [2]. For a system of N S harmonic oscillators, we have S(E) = (N S − 1)ln E + S 0 where S 0 is a constant [2], yielding a n = (−1) n−1 nT n (N S − 1) n−1 , 2 n.
For the one-dimensional model of magnetic moments, we do not have a closed-form expression for g(E) but we can use the equivalence of the microcanonical and canonical ensemble in the thermodynamic limit to derive closed-form expressions for the coefficients a n . The exact result, see equation (9) for a chain of N S spins, can at least, in principle, be inverted to yield T = T (E). Then, the calculation of the coefficients a n is straightforward and we find that where c = 1 − 1/[T sinh(1/T )] 2 and L = L(1/T ). As n increases, the number of terms in the expression of a n rapidly increases and, therefore, we only give the expressions for n 4. In our analysis, we used equation (20) including all terms n 8. For both models and for fixed subsystem size N S , the expansion coefficients a n are functions of only one parameter, namely the temperature T . This parameter, the position of the maximum E * and the number of terms in the expansion equation (20) completely determine the approximate form of p(E), as given by statistical mechanics. Fitting this approximate form to the distribution P(E s ) obtained by solving the Newtonian equations of motion yields the values of the fitting parameters E * and T S .
For the two models considered in this paper, the coefficients a n are simple functions of T and N S . Therefore, using T and E * as adjustable parameters a fit of equation (20) to the histogram P(E S ) obtained from the dynamical evolution of the subsystem yields an estimate of the temperature T S of the subsystem.  Poincaré-type map of the differences X = (x 0 − x 1 )/2 and P = ( p 0 − p 1 )/2 √ 2 for a ring of N = 512 oscillators, an integrable system, (left) and of the differences X = S x 1 − S x 2 and P = S z 1 − S z 2 for a ring of N = 512 magnetic moments, a chaotic system [22] (right). Each graph contains 5000 points.
If the estimate T S is indeed the temperature of the subsystem, the second central moment of P(E S ) should be related to the specific heat of the subsystem. To check this, we define where E S and E 2 S are the first and second moments of P(E S ), respectively.

Poincaré maps
The Poincaré maps for a system of four oscillators (figure 2 (left)) and magnetic moments (figure 2 (middle and right)) illustrate the difference in the dynamical behavior between the oscillator and magnetic moment model. As expected, the oscillators trace out a regular pattern in the (x, p) plane. Also, the trajectory of a magnetic moment of a ring of four sites is regular but, in contrast, if we open the ring to make a chain, the behavior changes from regular to chaotic, in agreement with theory [22].
However, the Poincaré maps of the difference of nearest-neighbor coordinates are very similar (see figure 3). In the normal-mode representation, each of the oscillators traces out an ellipse in the (X k , P k ) plane but in the original coordinates, this trace is much more complex, as  Clearly, it is difficult to detect some regularity in this map. Moreover, this map is very similar to figure 3 (right), the Poincaré-type map of a one-dimensional classical model of magnetic moments. Therefore, it is this similarity, not the presence of chaos, that seems to be responsible for the equilibration and thermalization of S (to be demonstrated below) in these two different systems.

Equilibration
We study the dynamic evolution of the subsystem S when it is brought in contact with the environment E. Initially, using one of the procedures described in sections 2.3 and 3.3, S and E are prepared such that they have a different energy. As the whole system S + E evolves in time, strictly conserving the total energy, we monitor the energy of the subsystem as a function of time.
Some representative results are shown in figure 4. From figure 4, it is clear that for both models, the energy of the subsystem E S rapidly approaches the energy of the whole system. These simulations also demonstrate that equipartition is established, a necessary but not a sufficient condition for the subsystem to be in thermal equilibrium.

Thermalization
Establishing that the subsystem is in thermal equilibrium requires a demonstration that the (normalized) frequency distribution P(E S /N S ) of subsystem energies, obtained by monitoring the energy E S of a single trajectory in the phase space of the subsystem at regular time intervals, and the statistical mechanical distribution p(E) are the same, meaning that the Kullback-Leibler distance [27] between these two distributions is very small. As shown in figure 5, the simulation data for P(E S /N S ) (red lines) and fitted p(E S /N S ) (black lines) are in excellent agreement, for both models alike. In the thermodynamic limit (N E → ∞ before N S → ∞), all but the quadratic term in the exponential can be neglected and the distribution is Gaussian [2]. Therefore, for large N S , T S and E * can be obtained by fitting a Gaussian to P(E S /N S ) but from figure 5, it is clear that for small subsystem sizes, P(E S /N S ) deviates significantly from a Gaussian. However, taking into account the higher-order terms in the expansion equation (20), the agreement between simulation data and the prediction of statistical mechanics is excellent. Repeating the simulations with different initial conditions (including different initial energies for the subsystem or the environment) strongly suggests that this agreement is generic (see the appendix).
As a conclusive test, we perform microcanonical Monte Carlo simulations for both models and obtain the distributions P mc (E S ) of the subsystems (see the appendix). The microcanonical Monte Carlo simulation generates statistically independent configurations of the whole system strictly according to the microcanonical distribution but samples the phase space in a completely different manner than does Newtonian dynamics. In all cases, we find that the difference between P mc (E S ) and P(E S ) is very small. For instance, for the system of oscillators with N S = 20, N = 65 536 and 10 8 samples, we find that the Kullback-Leibler distance D(P mc (E S ); P(E S )) ≈ 4 × 10 −2 , indicating that the probability that the two distributions are different is very small.

Discussion
Having established that the interaction of the environment with the subsystem causes both the environment and the subsystem to equilibrate and also drives the latter to its thermal equilibrium, it becomes possible to derive from the Newtonian dynamics alone estimates for the equilibration time. To this end we express the equilibration time estimated from the simulations in physical units. Typical frequencies of vibration in a solid are of the order of 10 11 Hz. Using this number to set the scale of the frequency in our model, we find that equilibration takes of the order of 10 −9 s. Similarly, for the system of magnetic moments, a realistic value of J/k B is of the order of 10 K, yielding an equilibration time of the order of 10 −8 s. Classical spin systems with Hamiltonians that encode frustration and/or disorder of regular or random kind are however, expected to exhibit larger, possibly much larger, equilibration time scales. The dynamical properties for subsystems of such systems under Newtonian evolution are beyond the scope of this paper.
Even though the subsystems and the environments which we have simulated are very small in the thermodynamic sense, the subsystem and environment equilibrate on a nanosecond time scale. Therefore, for an isolated nanoparticle of even a few thousands of atoms, an experimental probe that concentrates on only a few of those atoms should yield data that follow the canonical distribution. is derived on the basis of purely statistical considerations. Hence, the subsystem is in thermal equilibrium. Figure A.1 (right) shows the distribution P(E S /N S ) for the oscillator model but instead of analyzing a single trajectory in the phase space, microcanonical Monte Carlo simulation was used to sample the phase space. These data have also been obtained by taking 4 × 10 6 samples. Visual comparison of figure A.1 (left) and (right) suggests that these two very different methods of sampling points in phase space yield very similar results.
In figure A.2, we replot the distributions for the oscillator system ( figure A.1 (left)) as a function of rescaled and displaced energies and also plot the data for the variances as a function of the subsystem size N S . Both plots confirm that the data are in excellent agreement with the predictions of statistical mechanics, this in spite of the fact that the number of particles in the subsystem is only of the order of 100.
For reference, in table A.1 we collect the relevant numerical data extracted from the distributions shown in figure A.1 and also present results for the Kullback-Leibler distance between the data obtained by Newtonian dynamics, a Gaussian distribution and the results of the microcanonical Monte Carlo sampling. From table A.1, it is clear that sampling from a single trajectory in phase space yields results that are in excellent agreement with those of the canonical ensemble, even though the subsystem is rather small and the whole system belongs to the class of integrable models.   figure A.1 (left). The total energy per oscillator is E 0 /N = 0.998. D 1 is the Kullback-Leibler distance between the histogram P(E S ) and the theoretical prediction equation (20). D 2 is the Kullback-Leibler distance between the histogram P(E S ) and the Gaussian distribution with average E S and variance σ (E S ). D 3 is the Kullback-Leibler distance between the histogram P(E S ) obtained by solving the Newtonian equations of motion (see figure A.1 (left)) and the histogram obtained by the microcanonical simulation (see figure A.1 (right)). The results for the subsystem temperature, energy and specific heat extracted from the time-dependent data are in excellent agreement with the corresponding results of the canonical ensemble. In figure A.3 we present the results for the relaxation of only one oscillator brought out of equilibrium in an environment of many equilibrated oscillators. Figure A.3 (left) shows results for a system with 65 536 oscillators, which is prepared at equilibrium except that the momentum of oscillator j = 32 768 is increased by a factor of 100 such that the kinetic energy of this oscillator is about a factor 10 000 larger than the average kinetic energy at equilibrium. The main figure shows that the histograms for the momenta, computed according to the same procedure as the one used to obtain the histogram P(E S ), of two well-separated oscillators j = 32 768 (red line) and j = 1 (blue squares) are almost identical. The small differences can be further reduced by enlarging the environment and increasing the number of samples. The inset shows the initial time evolution of p 1 (blue line) and p 32 768 (red line) and demonstrates that p 32 768 equilibrates as a function of time. From figure A.3 it is clear that equilibration is re-established after bringing one oscillator out of equilibrium. Note that this is only the case if the environment is large enough (N E → ∞ before p j → ∞), as can be seen by comparing figure A.3 (left) and (right) for a system with 65 536 and 16 384 oscillators, respectively. A similar finding was reported by Bogoliubov [23], who showed that under the condition N → ∞, a single oscillator coupled linearly to a harmonic environment (containing N oscillators) described by the canonical distribution relaxes to the canonical distribution. In figure A.4, we present some representative results for the nonintegrable classical magnetic moment model with initial configurations of the subsystem and environment constructed using procedure 2 described in section 2.3. The solution of the Newtonian equations of motion is much more time consuming than in the case of the oscillator model and, therefore, we have limited the number of samples to 10 6 . Qualitatively, the histograms produced by the magnetic moment model show the same features as those of the oscillator model. However, due to the fact that the environment is much smaller (about 2000 instead of about 60 000 particles) and the statistics are less good, we may expect that there are minute differences in the tails of the distributions obtained by Newtonian dynamics and microcanonical Monte Carlo sampling. Although hard to see in the graphs, this expectation is clearly confirmed by the Kullback-Leibler distance D 3 , given in table A.2, which increases with the subsystem size.