Statistical Mechanics of an Integrable System

Abstract

We provide here an explicit example of Khinchin’s idea that the validity of equilibrium statistical mechanics in high dimensional systems does not depend on the details of the dynamics, as it is basically a matter of choosing the “proper” observables. This point of view is supported by extensive numerical simulation of the one-dimensional Toda chain, an integrable non-linear Hamiltonian system where all Lyapunov exponents are zero by definition. We study the relaxation to equilibrium starting from very atypical initial conditions and focusing on energy equipartion among Fourier modes, as done in the original Fermi-Pasta-Ulam-Tsingou numerical experiment. We consider other indicators of thermalization as well, e.g. Boltzmann-like probability distributions of energy and the behaviour of the specific heat as a function of temperature, which is compared to analytical predictions. We find evidence that in the general case, i.e., not in the perturbative regime where Toda and Fourier modes are close to each other, there is a fast reaching of thermal equilibrium in terms of a single temperature. We also find that equilibrium fluctuations, in particular the behaviour of specific heat as function of temperature, are in agreement with analytic predictions drawn from the ordinary Gibbs ensemble. The result has no conflict with the established validity of the Generalized Gibbs Ensemble for the Toda model, which is on the contrary characterized by an extensive number of different temperatures. Our results suggest thus that even an integrable Hamiltonian system reaches thermalization on the constant energy hypersurface, provided that the considered observables do not strongly depend on one or few of the conserved quantities. This suggests that dynamical chaos is irrelevant for thermalization in the large-N limit, where any macroscopic observable reads of as a collective variable with respect to the coordinate which diagonalize the Hamiltonian. The possibility for our results to be relevant for the problem of thermalization in generic quantum systems, i.e., non-integrable ones, is commented at the end.

Appendix

1.1 Integrals of Motion

Just for the sake of self-consistency we briefly remind here the conservation laws of the Toda system. In the case of a periodic lattice with N particles where $q_0=q_N$ the first four integrals of motion read as:

$$\begin{aligned} \mathscr {I}_1&= \sum _{i=1}^N p_i \nonumber \\ \mathscr {I}_2&= \sum _{i=1}^N ~\frac{p_i^2}{2} + \mathscr {U}_i \nonumber \\ \mathscr {I}_3&= \sum _{i=1}^N ~\frac{p_i^3}{3} + (p_i+p_{i+1})\mathscr {U}_i \nonumber \\ \mathscr {I}_4&= \sum _{i=1}^N ~\frac{p_i^4}{4} + (p_i^2+p_{i+1}^2+p_ip_{i+1})\mathscr {U}_i + \frac{\mathscr {U}_i^2}{2}+\mathscr {U}_i\mathscr {U}_{i+1}, \end{aligned}$$

(23)

where $\mathscr {U}_i = e^{-(q_{i+1}-q_i)}$. Note that $\mathscr {I}_2$ is nothing but the energy (apart for an irrelevant additive constant). The expressions in Eq. (23), see, e.g., [39], can be extended to the case of fixed boundary conditions by conventionally defining a periodic chain with $2N+2$ particles where the coordinates from $q_0$ to $q_{N+1}$ are identical to the original one, while those such as $q_{N+1+i}$ with $i=1,\ldots ,N$ are defined in order to have an odd function with respect to the center of the lattice:

$$\begin{aligned} q_{N+1+i}&= -q_{N+1-i} \nonumber \\ p_{N+1+i}&= -p_{N+1-i} \end{aligned}$$

(24)

1.2 Specific Heats

The first step to compute the specific heat is the calculation of the partition function. For consistency with the numerical simulations, we compute it with fixed boundary conditions. That is, we have N variables $q_i$ with $i=1,\ldots ,N$ to integrate over, while

$$\begin{aligned} q_0 = q_{N+1} = 0 \end{aligned}$$

(25)

By identifying with the same variable the two fixed boundaries, i.e., $q_0=q_{N+1}$, we can thus change variables to

$$\begin{aligned} r_i = q_{i+1} - q_i, \end{aligned}$$

(26)

with, in particular

$$\begin{aligned} r_{N}&= q_{N+1} - q_N, \nonumber \\ r_{N+1}&= q_1 - q_{N+1}. \end{aligned}$$

(27)

We can thus safely integrate over the $N+1$ variables $r_i$ under the global constraint $\sum _{i=1}^{N+1}r_i = 0$:

$$\begin{aligned} \mathscr {Z}_N^{(\text {P})}(\beta ) = e^{\beta N} \int _{-\infty }^\infty \prod _{i=1}^{N+1} dr_i~e^{-\beta \sum _{i=1}^{N+1} \exp (-r_i)}~\delta \left( \sum _{i=1}^{N+1} r_i \right) . \end{aligned}$$

(28)

In order to unfold the global constraint it is useful to exploit the inverse Laplace rapresentation of the partition function:

$$\begin{aligned} \mathscr {Z}_N^{(\text {P})}(\beta )&= e^{\beta N}\int _{s_0-i\infty }^{s_0+i\infty } \frac{ds}{2\pi i} \int _{-\infty }^\infty \prod _{i=1}^{N+1} dr_i~ \exp \left( -\beta \sum _{i=1}^{N+1} \exp (-r_i) + s \sum _{i=1}^{N+1} r_i \right) \nonumber \\&= e^{\beta N}\int _{s_0-i\infty }^{s_0+i\infty } \frac{ds}{2\pi i} \left[ \int _{-\infty }^\infty dr~\exp \left( -\beta e^{-r} + s r\right) \right] ^{N+1} \nonumber \\&= e^{\beta N}\int _{s_0-i\infty }^{s_0+i\infty } \frac{ds}{2\pi i}~\exp \left\{ N \log z_\beta (s) \right\} \end{aligned}$$

(29)

where we have introduced the partition function per degree of freedom:

$$\begin{aligned} z_\beta (s) = \int _{-\infty }^\infty dr~\exp \left( -\beta e^{-r} + s r\right) . \end{aligned}$$

(30)

It is convenient to regard $\beta$ as a parameter in Eq. (30) and consider $z_\beta (s)$ as a function of s in the complex s plane. Clearly the integral expression in Eq. (30) is well defined only for $\text {Re}(s)< 0$, so that it is convenient to change variable in the integral of Eq. (29) from s to $\mu = -s$:

$$\begin{aligned} \mathscr {Z}_N^{(\text {P})}(\beta ) = e^{\beta N}\int _{\mu _0-i\infty }^{\mu _0+i\infty } \frac{d\mu }{2\pi i}~ \exp \left\{ N \log \hat{z}_\beta (\mu ) \right\} , \end{aligned}$$

(31)

where now the contour in the complex plane passes through $\text {Re}(\mu _0)>0$ and we have the partition function per degree of freedom defined as:

$$\begin{aligned} \hat{z}_\beta (\mu ) = \int _{-\infty }^\infty dr~\exp \left( -\beta e^{-r} - \mu r\right) . \end{aligned}$$

(32)

The partition function of the whole system $Z_N(\beta )$ can be computed from Eq. (31) by means of a saddle-point approximation in the large-N limit. But first we need to compute explicitly $\hat{z}_\beta (\mu )$ in Eq. (32). By going through two changes of variables, i.e.,

$$\begin{aligned} x\rightarrow u&= e^{-x} \nonumber \\ u \rightarrow t&= \beta u, \end{aligned}$$

(33)

we obtain

$$\begin{aligned} \hat{z}_\beta (\mu )&= \int _0^\infty du~u^{\mu -1}e^{-\beta u} = \frac{1}{\beta ^\mu } \int _0^\infty dt~t^{\mu -1}~e^{-t} \nonumber \\ \hat{z}_\beta (\mu )&= \frac{\Gamma (\mu )}{\beta ^\mu }. \end{aligned}$$

(34)

In order to compute the integral in Eq. (31) by means of the saddle-point approximation we just need to solve the saddle-point equation

$$\begin{aligned}&\frac{\partial }{\partial \mu }\left( \log \Gamma (\mu ) - \mu \log (\beta ) \right) = 0 \nonumber \\&\Longrightarrow ~\frac{\Gamma '(\mu )}{\Gamma (\mu )} = \log (\beta ) \end{aligned}$$

(35)

At this stage, without doing any approximation and indicating as $\hat{\mu }(\beta )$ the function implicitly determined by the condition in Eq. (35), one has that the partition function of the system reads:

$$\begin{aligned} \mathscr {Z}_N^{(\text {P})}(\beta ) = e^{\beta N} \frac{\Gamma (\hat{\mu }(\beta ))}{\beta ^{\hat{\mu }(\beta )}}. \end{aligned}$$

(36)

For instance the specific heat $C_V$, which can be obtained from Eq. (36) using the standard formula in Eq. (18) of the main text, reads as

$$\begin{aligned} C_V = - \frac{1}{\psi '(\hat{\mu }(\beta ))} - \left[ \psi (\hat{\mu }(\beta ))- \log (\beta ) \right] \left( \frac{\psi ^{''}(\hat{\mu }(\beta ))}{[\psi '(\hat{\mu }(\beta ))]^3} + \frac{1}{\psi '(\hat{\mu }(\beta ))} \right) + \hat{\mu }(\beta ), \end{aligned}$$

(37)

being $\psi (x)$ the so-called digamma function $\psi (x) = \Gamma '(x)/\Gamma (x)$. The expression in Eq. (37) is not very insightful, but is necessary to plot the full analytic behavior of the function. Similarly one can write down the expression for the temperature, but we leave this exercise to the willing reader.

On the contrary, what can be obtained explicitly, are the asymptotic behaviors for both the average potential energy $\langle u \rangle$ and the specific heat $C_V$ in the limit of small and large temperatures T. This is obtained by the knowledge of the two asymptotic behaviors of the digamma function $\psi (\mu )$:

$$\begin{aligned} \mu \rightarrow \infty ~~~&\Longrightarrow ~~~\psi (\mu ) \sim \log (\mu ) \nonumber \\ \mu \rightarrow 0~~~&\Longrightarrow ~~~\psi (\mu ) \sim - \frac{1}{\mu }. \end{aligned}$$

(38)

We thus have two different approximations for our saddle-point equations in the low-temperature, $\beta \rightarrow \infty$, or the high-temperature, $\beta \rightarrow 0$, limit:

$$\begin{aligned} \text {Low T:}~~~\Longrightarrow ~~~ \log (\mu )&= \log (\beta )~~~\Longrightarrow ~~~\mu =\beta \nonumber \\ \text {High T:}~~~\Longrightarrow ~~~ -\frac{1}{\mu }&= \log (\beta )~~~\Longrightarrow ~~~\mu =-\frac{1}{\log (\beta )} \nonumber \\ \end{aligned}$$

(39)

Let us first consider the low-temperature case, where Toda is well-approximated by a system of weakly coupled harmonic oscillators, for which the specific heat is constant. In this case we get:

$$\begin{aligned} Z_N(\beta ) \approx \exp \left( (N+1) \left[ \beta + \log \Gamma (\beta ) -\beta \log (\beta ) \right] \right) . \end{aligned}$$

(40)

By using the Stirling’s approximation

$$\begin{aligned} \log \Gamma (\beta ) = \beta \log (\beta ) - \beta -\frac{1}{2} \log (\beta ), \end{aligned}$$

(41)

we thus get

$$\begin{aligned} Z_N(\beta ) \approx \exp \left[ - \frac{N}{2}\log (\beta ) \right] . \end{aligned}$$

(42)

From Eq. (42) and from the definition of $C_V$ in Eq. (18) in the main text we get

$$\begin{aligned} C_V = \beta ^2\frac{\partial ^2}{\partial \beta ^2}\left( \frac{N}{2}\log (\beta ) \right) = \frac{1}{2}. \end{aligned}$$

(43)

In the light of the high-temperature expansion, we get

$$\begin{aligned} \beta \rightarrow 0 \Longrightarrow \mathscr {Z}_N(\beta ) \approx \left\{ N\left[ \beta +\log \Gamma \left( -\frac{1}{\log (\beta )}\right) \right] \right\} . \end{aligned}$$

(44)

By exploiting the identity of the Euler gamma function $z\Gamma (z)=\Gamma (z+1)$ one has that

$$\begin{aligned} z \rightarrow 0 ~~~\Longrightarrow ~~~\Gamma (z) \sim \frac{1}{z}, \end{aligned}$$

(45)

from which we have

$$\begin{aligned}&\log \mathscr {Z}_N(\beta ) \approx \left\{ N\left[ \beta +\log \left( -\log (\beta ) \right) \right] \right\} \nonumber \\&\Longrightarrow ~~~C_V \sim \frac{1}{\log (T)} \end{aligned}$$

(46)

In a similar manner one obtains the behavior of the average energy $\langle u \rangle$, which is reported in the main text.