Prethermalization in Fermi-Pasta-Ulam-Tsingou chains

Since (1) has $N$ degrees of freedom its harmonic limit has a set of $N$ eigenmodes, denoted $\{P_{j},Q_{j}\}_{j}$, with which we can rewrite the position and momentum in (1) at any time $t$ as

each mode having eigenfrequency and eigenenergy

In these canonical coordinates, the Hamiltonian (1) takes the form [19]

where [28]

Note that $B_{i,j,k}$ is invariant with respect to permutations of its indices.

Eq. (1) is much easier to deal with, numerically speaking, than (II.1). However, we are interested in several properties that need to be computed in the eigenmode basis, such as the eigenenergies themselves and the [normalized] spectral entropy

where the Shannon entropy and normalized eigenenergy are given, respectively, by

Spectral entropy has been used several times in order to quantify equipartition in nonlinear systems [23, 4, 29, 30]. The logic for employing it is simple: If ideal equipartition is reached, $\rho_{j}(t)$ will be stationary and the same for every mode $j$, freezing $\eta(t)$ at some value $\langle\eta\rangle$. If one assumes a Gibbs distribution at equilibrium, $\langle\eta\rangle$ can be calculated analytically as a function of the Shannon entropy of the initial state, being given by

where $\gamma{\approx}0.5772$ is the Euler constant [31]. Since $\eta(0){=}1$ for any initial state, one expects $\eta(t)$ to decrease and follow a complicated and possibly highly oscillatory evolution until eventually hitting its stationary value. Once this happens the dynamics cannot, at least on average, depart significantly from thermal oscillations around $\langle\eta\rangle$. It is then natural to consider that equipartition is reached once $\eta(t)$ hits its equilibrium value for the first time [30]. Such first-passage time furnishes our definition of equipartition time, $\tau_{\mathrm{eq}}$, and is appealing because it does not require any ad hoc assumptions.

Equipartition time as defined in Sec. II.2 is, evidently, attached to the observables being monitored, namely the mode energies. Different observables will generally reach equipartition at different times (if at all), such that observable equipartition does not actually measure properties intrinsic to the system at hand. Nevertheless, since mode energies are conserved in the unperturbed harmonic chain, one knows well how these observables should behave upon approaching an integrable limit, i.e., equipartition should slow down and eventually stop for sufficiently small energy densities.

There is, however, an observable-independent time scale that is completely intrinsic to the system and should not depend on the initial state, given by the Lyapunov time (the inverse of the maximal Lyapunov exponent). Since, like the maximal Lyapunov exponent (MLE), all exponents in the Lyapunov spectrum (LS) are invariant with respect to homeomorphisms [27], the LS provides a valuable set of time scales related to the average “pull” along each 1-dimensional submanifold that forms the hyperbolic skeleton of the system’s dynamics in phase space. The information in the LS can also be condensed in the form of a scalar, the Kolmogorov-Sinai entropy, which by Pesin’s theorem can be obtained as the sum of all positive exponents in the LS [32].

In order to compute the LS we employ the well-known prescription of [33, 34], which amounts to computing the spectrum of the symmetric matrix

In the above $\mathcal{M}$ represents the system’s stability matrix, which is obtained by solving Hamilton’s equations in tangent space:

where $\mathrm{Hess}$ is the Hessian with respect to $(\boldsymbol{p},\boldsymbol{x})$. Since the stability matrix has to be computed along a trajectory, the above equation is solved in parallel with Hamilton’s equations for an initial phase-space point $(\boldsymbol{p}_{0},\boldsymbol{x}_{0})$ until a final time that is long enough to result in a converged LS. During this evolution, one must also perform QR-diagonalizations of $\mathcal{M}$ in order to avoid numerical denegenacies, as is well-known and described in, e.g., [34].

As can be seen in (9), the stability matrix will generally depend on the trajectory along which it is calculated, i.e., on the initial phase-space point $(\boldsymbol{p}_{0},\boldsymbol{x}_{0})$. However, for a sufficiently large number of degrees of freedom and energy density, phase space will consist of a single metrically intransitive set that is equally accessible to all trajectories in the system [35, 36]. Equivalently, one can formulate this previous property as stating that, for sufficiently large $N$ and $h{=}E/N$, the system lies outside the KAM regime in which phase space can be split into a chaotic web and pockets of near-integrability. Thus, since all trajectories are “free to roam” around the entire accessible phase space, the LS is expected to be independent of the trajectory along which it is computed. This independence is a stronger evidence that the system at hand is ergodic than computing observables, and can be ascertained directly from the verification that the KS entropy is independent of the trajectory along which it is computed.

Initializing the system in a root $Q_{r}$ means setting $x_{j}{=}0$ unless $j{=}r$. The first excitation cycle, therefore, starts with a lone mode $Q_{r}$, which by (III) excites modes $i=2r$ and $i=2(N-r+1)$. If the first condition is met, then the second delta will be one iff $r=(N+1)/2$, with $N+1$ even. Thus, if the chain has an odd number of sites and the root chosen is $Q_{(N+1)/2}$, then all coefficients in the coupling tensor vanish and the energy remains forever localized in the root, i.e. this mode is a q-breather.

Now, assuming that the second Kronecker delta in (III) vanishes, i.e. that $i\neq 2(N-r+1)$, then the only mode that is excited starting from $Q_{r}$ is $Q_{2r}$. This constitutes the first excitation cycle, namely the first run over $j$ and $k$ indices in (11). Note that $Q_{r}$ excites $Q_{2r}$ only if $2(r-1)<N$, since $i$ cannot be negative or zero. This justifies splitting excitation sequences into two types, namely low- and high-mode excitations, depending on whether the root $r$ fulfills $2(r-1)<N$ or not. The only difference is that for low-mode excitations the initial sequence will be ascending, while it is descending in the high-mode case. Although these scenarios are essentially identical, we will assume low-mode excitations here in order for sequences to start in the ascending phase.

The second cycle for low-mode excitations starts with two excited modes, $Q_{r}$ and $Q_{2r}$. Evidently, $Q_{r}$ will re-excite $Q_{2r}$, which on the other hand will hit $Q_{4r}$ due to (III). According to (III) mode $Q_{3r=(2r+r)}$ will also be excited. Thus, the second excitation cycle results in the sequence $\{(Q_{r}),(Q_{2r}),(Q_{3r},Q_{4r})\}$. For the third cycle, $Q_{6r}$ and $Q_{8r}$ will be present due to (III), together with all terms that can be formed from $j+k$ sums as dictated by (III). In this ascending phase, terms created from $|j-k|$ are not important, since any $j$ and $k$ will be multiples of $r$ that were already excited in previous cycles. Thus, by induction, the first ascending sequence starting from a root $Q_{r}$, with $2(r-1)<N$, is simply

until a reflection takes place in (III), i.e. until $j+k=2(N+1)$ for some $x_{j}$ and $Q_{k}$. The reflection, however, will only excite new modes if $2(N+1)-(j+k)$ is not a multiple of $r$, since otherwise the terms would have been already excited during the ascending sequence. This will be explored in the following.

At the end of the first ascending sequence every term is a multiple of the root, $r$. Thus, at the reflection, substituting $j+k=a\,r$ in the last delta of (III) singles out $r=2(N+1)/a$, allowing us to state two important facts:

1. 1.

   If $r$ divides $2(N+1)$ the first ascending sequence will be reflected into itself (for low modes) or be equivalent to a low mode excitation sequence that reflects into itself (for high modes). This shows, in particular, that $Q_{2}$ is nonthermal for all $N$, since $2(N+1)$ is a always divisible by 2. It also shows that $Q_{1}$ is a trivial thermal root, since it excites all modes in the first ascending sequence for any $N$;

2. 2.

   The reflection takes place when a multiple of $r$ is larger than $N$ for the first time, i.e. when $x_{j}=Q_{(a+1)r}$. Evidently this mode lies outside mode space and the excited mode is actually $Q_{2(N+1)-(a+1)\,r}$ after the reflection. A consequence is that even roots can never excite odd modes, since for even modes $2(N+1)-(a+1)\,r$ is always even. This shows that $Q_{2r}$ are nonthermal roots for all $r$ and all $N$.

Once a reflection occurs, the selection rule for $|j-k|$ generates a descending sequence starting from $Q_{2(N+1)-(a+1)\,r}$, namely,

The descent continues in multiples of $r$ until the next mode would exit the chain, i.e. until $2(N+1)-(a+b)\,r<0$.

Once the first descending sequence is over, another reflection takes place, now going from left to right, mediated by $\delta_{i,|j-k|}$ in (III). Substituting the last mode number from the descending sequence in this delta, we see that the excited mode after this reflection is $Q_{2(N+1)-(a+b)\,r+r}$. The second ascending sequence will start from this mode number, once again ascending in steps of $r$ until $2(N+1)-(a+b-c)r>N$. Then, we have another reflection from left to right that excites mode $2(N+1)-[2(N+1)-(a+b-c)r-r]=(a+b-c+1)r$, etc. Instead of writing full abstract sequences, we now provide a few examples showing how easy such sequences are to write in practice.

Let us start with a chain with $N=11$ masses. Since $2(N+1)=24$, the only three nontrivial thermal roots in this chain are $Q_{5}$, $Q_{7}$ and $Q_{11}$. The root $Q_{3}$ is evidently nonthermal since it divides 24, but let’s take a look at its excitation sequence. The first ascending sequence is $\{(Q_{3}),(Q_{6}),(Q_{9})\}$, and we stop here because the next term falls outside mode space. The reflection term comes from $i=2(N+1)-(6+9)=9$, which is already excited. This will be now propagated backwards to $Q_{6}$ and $Q_{3}$ accoding to (III), retracing the ascending sequence and never exciting any new mode. Now, for the leftmost reflection, the excited mode is $Q_{|3-3|}$, which does not exist. Thus, the bush of $Q_{3}$ is given by $\{(Q_{3}),(Q_{6}),(Q_{9})\}$. This is verified numerically in Fig. 1(a).

Let us now consider what happens to $Q_{3}$ for $N=12$ instead. The first ascending sequence is now updated to $\{(Q_{3}),(Q_{6}),(Q_{9},Q_{12})\}$. The reflection gives $i=2(N+1)-(3+12)=11$, which was not previously excited and is not a multiple of $r=3$. This will then propagate backwards to $11-3=8$, $11-6=5$, etc. The full excitation sequence is then given by

where the arrows indicate whether the sequence is ascending ($\nearrow$) or descending ($\swarrow$). This excitation sequence is seen to match numerics in Fig. 1(b).

As another example, for $Q_{3}$ and $N=13$ we have

In both (III.4) and (III.4) the root $Q_{3}$ resulted in three excitation sequences: two ascending and one descending. If $Q_{5}$ is chosen as the root, then the full excitation sequence for $N=13$ is

In general, it is easy to see that since every step in the sequence is performed in multiples of $r$, one needs $r$ ascents and descents (or $r-1$ reflections) to cover mode space for a thermal root $Q_{r}$, no matter what $N$ is chosen. Since large $r$ requires more cycles to cover mode space, it means less modes are excited at the end of each descending or ascending sequence. One can once again use the case of $N=13$ as an example: Here, at the end of the first ascending sequence, the root $Q_{1}$ has already excited all modes in the chain, which is a unique property of the fundamental root; $Q_{3}$, however, has only excited 5 of the 13 available modes, and $Q_{5}$ only 3 out of 13. In Sec. IV we will show that this effect strongly impacts any quantity computed from evolving roots.

Sec. III described how energy can remain localized within a small subset of modes before spreading throughout mode space: The higher excited the root, the smaller the initial subset of excited modes, and also the slower the energy transfer. During this transfer time the dynamics is mostly localized within the excited modes, which essentially “wait” while slowly transferring energy to the ones in the next cycles, forming the metastable plateaus well known in FPUT literature [17, 38]. The longer it takes for the energy to reach the end of the sequence, the longer a metastable plateau survives.

The well-known fact that the metastable plateaus associated to roots become more visible at lower energy densities is unsurprising, since in the limit of zero energy density such modes are normal modes of the associated harmonic lattice. However, the behavior of plateaus for excited roots was not previously addressed. To this end, Fig. 2 shows a comparison of evolving mode energies and spectral entropies for the prototypical $\alpha$-FPUT chain (1) computed starting from the first 5 odd roots in a chain with $N=63$ masses, together with a random initial state. Fig. 2(a) displays the standard mode energy dynamics found in literature: An excitation of $Q_{1}$ monotonically spreads to all other modes and eventually results in energy equipartition, with all $E_{k}$’s converging to the mean energy $E/N=h$ (black dashed line). Metastable plateaus are clearly visible in the prethermal regime, in which a measurement would indicate that the system is not ergodic and cannot be described by equilibrium statistical mechanics. However, such a statement is evidently incorrect, as waiting a bit longer would result in a system that does display statistical behavior.

The approach from an out-of-equilibrium initial state towards a thermal one is most clearly seen in tracking the spectral entropy as a function of time, shown in the insets: The moment it touches the red line, which represents its equilibrium value (8), is where we consider equipartition to have been achieved. The associated equipartition times are displayed in the main plots as an orange triangle and slowly increase as we move from $Q_{1}$ up to $Q_{9}$, with the shortest time corresponding to the random initial condition. Evidently, this latter case is less interesting since the equipartition time depends strongly on the initial sampling, but for roots the structures formed by mode energies as they evolve in time show that, indeed, higher-excited roots isolate an increasingly smaller subset of modes in mode space. Fig. 2(a) also shows a dense subset of modes homogeneously concentrated around the mean, while Fig. 2(e) clearly shows that two modes, namely $Q_{9}$ and $Q_{18}$, remain relatively isolated from the others for much longer times. To see this more clearly, in Fig. 3 we take vertical slices of Fig. 2 at three different fixed times for $Q_{1}$, $Q_{3}$ and $Q_{9}$. At $t{=}\tau_{\mathrm{eq}}$ low frequency modes that are multiples of 3 and 9 can still be seen to have above-average energies for roots $Q_{3}$ and $Q_{9}$. At $t{=}10\tau_{\mathrm{eq}}$ this effect has already dissipated for $Q_{3}$, but the energies of modes $Q_{9}$ and $Q_{18}$ remain slightly above the mean energy when evolving $Q_{9}$, as visible in Fig. 3(c).

Evidently, the longer it takes for a subset of modes to reach equipartition, the longer the prethermal regime. This is clearly seen in Fig. 2: The relatively obfuscated metastability associated to the $Q_{1}$-root seen in panel (a) is already maximized in (b), where $Q_{3}$ is evolved, and becomes more and more blatant as higher-excited modes are used as roots. Evidently, the random initial condition in Fig. 2(e) does not display metastability, since here all modes are excited at once and the slow, selection-rule mediated transferring of energy from excited modes to previously unexcited ones does not take place. Since the mean energy is a constant of motion, the prethermal regime’s duration can also be visualized by plotting the standard deviation of $E_{k}(t)$, $\sigma[E_{k}(t)]$, as a function of time, which will be approximately constant while the energy is still strongly localized in a subset of normal modes. This is indeed seen to be true in Fig. 4, where it is also clear that the metastability’s lifetime increases monotonically as a function of root excitation for this energy density.

Spectral entropy (6) is a function of the initial state, so it is no surprise that the equipartition times computed from it in Fig. 2 depend on it. These times have no direct relation to the duration of prethermalization, except for the fact that they are necessarily longer than the length of metastable plateaus. Since random initial conditions do not, at least generally, undergo prethermal regimes, it is expected that the equipartition times of such typical initial conditions will be shorter than the ones of roots for any energy density chosen. This expectation is confirmed in Fig. 5, where we compare equipartition times for 280 typical conditions with the ones obtained from odd roots as a function of the energy density. Clearly, the former are bounded from above by the latter, which slowly increase as a function of excitation number.

If equipartition times reach the order of the propagation time used to compute them, namely $\mathcal{O}(10^{8})$, convergence can no longer be achieved for these parameters, which is clear in Fig. 5 by the presence of missing data for some values of energy density. However, the “disappearance” of points does not increase fully monotonically with mode excitation number, with some data displaying what could be described as “erratic behavior”. Since dependence on initial conditions is tantamount to ergodicity breaking, investigating quantities that are proved to be invariant in any ergodic system will shed light into the source of the instabilities seen in Fig. 5. The most important of such quantities is the MLE, $\lambda_{1}$, which indicates the presence of chaotic behavior and is not only invariant with respect to diffeomorphisms but also constant for ergodic dynamical systems [27]. From this exponent one can extract the Lyapunov time, $\tau_{1}=1/\lambda_{1}$, which measures the time a vector tangent to a trajectory takes to align with the maximally chaotic, one-dimensional submanifold in the tangent bundle [39, 40]. Since the Lyapunov time is computed from the divergence of two initially close trajectories, it takes a finite time to converge and it is often useful to observe its evolution as a function of time, as done in Fig. 6.

Panels (a) and (b) of Fig. 6 display the evolution of $\tau_{1}(t)$ extracted from propagating odd roots and random initial conditions for an energy density of $h{\approx}0.085$. Not only is convergence achieved for all initial states employed, but the time-dependent portrait of both typical and atypical initial conditions is identical. Upon decreasing the energy density, however, Fig. 6(c) shows that atypical initial conditions start behaving very differently not only with respect to typical ones, but also among themselves. While the Lyapunov times obtained from $Q_{1}$ and $Q_{3}$ are similar to the ones of the random initial conditions in Fig. 6(d), the data obtained from evolving $Q_{5}$ is seen to start forming a “belly” before convergence. Thus, it shows that $Q_{5}$ lies in a region of phase space that takes longer to align with the maximally chaotic direction, similarly to the trapping times described in [23]. The difference is that here we observe longer trapping times for a fixed energy density by changing only the initial conditions.

By evolving even higher excited roots such as $Q_{7}$, Fig. 6(c) shows that the “belly” increases when compared to the one of $Q_{5}$, such that $Q_{7}$ is trapped for a time one order of magnitude longer. Extrapolating to $Q_{9}$ this transient trapping regime, which is a manifestation of prethermalization in an invariant quantity and therefore must disappear for sufficiently long times, ends up being longer than the propagation time and does not converge in our computations. The reason for the longer transient times as we go up from $Q_{1}$ to $Q_{9}$ is that the higher the excitation number, the smaller the subspace of mode space in which dynamics is approximately restricted to during the prethermal regime, as discussed in Sec. III. Indeed, Fig. 3 shows that a large fraction of the energy is still concentrated in multiples of the initial excitation number at equipartition time, which as seen by comparing Figs. 6 and 7 is always longer than Lyapunov time when dealing with atypical initial conditions. Thus, as long as the root is a low lying mode (that is, $Q_{r}$ with $2(r-1)<N$), the trapping time will increase with the excitation number. The breather at $Q_{32}$ marks the boundary between low-mode and high-mode excitations, with $Q_{31}$ being the slowest mode to thermalize for $N{=}63$–so slow that there is no point in attempting to obtain data from it, since this is already computationally hard for $Q_{9}$. It is interesting to note that components of the trajectory that “scan” phase space in order to obtain convergence in $\tau_{1}$ are the low-energy ripples that are not multiples of $r$ for a root $Q_{r}$. Evidently, these ripples might not have enough energy to go very far, requiring longer transient times. The energy is spread among all modes most efficiently when evolving $Q_{1}$, which explains why this state behaves essentially as a random state despite its extremely long equipartition time when compared to the latter.

Fig. 7, which shows the final (converged) Lyapunov times as a function of energy density as obtained from typical and atypical initial conditions, is clear evidence that the system we are dealing with is numerically ergodic: The times are the same, no matter what initial conditions are chosen (although the larger error bars and spread for $h{\approx}10^{-2}$ show that $t_{f}{=}10^{9}$ is barely long enough to reach numerical ergodicity for the smallest values of $h$ used here). This does not mean, as discussed before, that all data for roots converges, with the main reason for failure being the presence of extremely long prethermal regimes as seen in Fig. 6.

Given the imprint prethermalization has in a dynamical invariant such as the Lyapunov time, in Fig. 8 we display the time evolution of one of the most meaningful, yet computationally expensive, invariant quantities in dynamical systems, namely the Kolmogorov-Sinai entropy, $\kappa_{\mathrm{KS}}$. This quantity, just like the Lyapunov time and the associated maximal exponent, displays a clear dip when computed by propagating roots, while for a typical initial condition it approaches its converged value from above. Since $\kappa_{\mathrm{KS}}$ is the sum of all positive Lyapunov exponents, what Fig. 8 shows is that the prethermal regime is present in all phase-space directions, not only the maximal one associated to the MLE. Upon decreasing the energy density we also see a tendency of a larger dip for higher-excited modes, in accordance with our expectation that the prethermal regime is more pronounced for higher excitation numbers. We also note that convergence of Kolmogorov-Sinai entropies to the same final value is likely the strongeast display of numerical ergodicity achievable in simulations.