Quantization of Integrable and Chaotic Three-Particle Fermi–Pasta–Ulam–Tsingou Models

We study the transition from integrability to chaos for the three-particle Fermi–Pasta–Ulam–Tsingou (FPUT) model. We can show that both the quartic β-FPUT model (α=0) and the cubic one (β=0) are integrable by introducing an appropriate Fourier representation to express the nonlinear terms of the Hamiltonian. For generic values of α and β, the model is non-integrable and displays a mixed phase space with both chaotic and regular trajectories. In the classical case, chaos is diagnosed by the investigation of Poincaré sections. In the quantum case, the level spacing statistics in the energy basis belongs to the Gaussian orthogonal ensemble in the chaotic regime, and crosses over to Poissonian behavior in the quasi-integrable low-energy limit. In the chaotic part of the spectrum, two generic observables obey the eigenstate thermalization hypothesis.


Introduction
The Fermi-Pasta-Ulam-Tsingou (FPUT) model [1] has widely been investigated with the aim of understanding the approach to thermodynamic equilibrium in an isolated system [2,3]. Since the beginning of these studies, the relation with integrability in Hamiltonian mechanics and the Kolmogorov-Arnold-Moser theorem [4][5][6] has been stressed [7]. More recently, the "closeness" of the model to the integrable Toda lattice has carefully been examined [8].
Here, we take a different approach to integrability in the FPUT model, with the aim of discussing its quantization. We study the N-particle model for periodic boundary conditions and we introduce a specific Fourier transform, which allows us to conveniently write the nonlinear interaction terms among the Fourier mode variables [9]. We then consider the three-particle case N = 3 and prove that it is integrable when either the cubic nonlinearity vanishes or the quartic nonlinearity does. Interestingly, the model with only quartic nonlinearity, dubbed β-FPUT model, has an additional integral of the motion besides energy and momentum. This integral, previously discovered in Ref. [10] using an analogy with celestial mechanics, can be associated to a cylindrical symmetry, which is only manifested in Fourier modes. Already for N = 4, this integral is no longer present, and what remains for higher values of N, even in the thermodynamic limit, is the presence of invariant subspaces in the Fourier mode space [9], whose stability has carefully been examined by Chechin and coworkers (see Ref. [11] and Refs. therein). It should be here mentioned that a study of the Einstein-Brillouin-Keller (EBK) quantization of the integrable N = 3 Toda lattice was undertaken by Isola, Kantz and Livi [12]. The complete (α + β)-FPUT model, with both non-zero cubic and quartic nonlinearity is non-integrable and displays strong chaotic motion in an intermediate energy range, as we show numerically by displaying Poincaré sections.
The quantization of classical chaotic systems is still an extremely active field of research, which has its roots in the discovery of the phenomenon of dynamical localization by Casati, Chirikov, Ford and Izrailev [13] and in the pioneering investigation of level statistics for Sinai billiards by Berry [14]. This latter topic was further investigated by Berry and Tabor [15], who characterized specifically the Poisson statistics in the integrable limit. These studies led to the famous universality hypothesis by Bohigas, Giannoni and Schmit [16], who first suggested that the spectra of quantum chaotic systems are in agreement with the predictions of random matrix theory (an early account of these developments can be found in Tabor's book [17]). The analysis of spectral statistics for smooth Hamiltonians continued with the work of Seligman, Verbaarschot and Zirnbauer [18], who studied the transition from regular to chaotic phase-space and fitted the levels with Poisson or Gaussian orthogonal ensemble (GOE) statistics, respectively, for an even parity potential.
In this paper, we quantize the three-particle FPUT model in the formalism of first quantization, and discuss its transition from integrability to chaos using methods from random matrix theory. Moreover, we check the eigenstate thermalization hypothesis (ETH) [19,20] for two relevant observables, kinetic energy and the additional integral of the motion. This complements works about the bosonic FPUT model [21][22][23], and in particular a recent study [24] investigating the localization transition in the thermodynamic limit.
The paper is organized as follows. In Section 2, we introduce the appropriate Fourier mode variables and prove the integrability of both the α-and β-FPUT N = 3 model. In Section 3, we quantize the three-particle FPUT model and study its level statistics. Section 4 is devoted to conclusions.

The FPUT Model and the Integrable Three-Particle Case
The Hamiltonian of the (α + β)-Fermi-Pasta-Ulam-Tsingou coupled oscillator model is where q i is the relative displacements of the i-th oscillator from its equilibrium position and p i the conjugate momentum. We choose periodic boundary conditions q N+1 = q 1 , p N+1 = p 1 . The unperturbed Hamiltonian can be written in quadratic form where t denotes the transposed of a column vector and A is the NxN symmetric matrix whose elements are Due to the periodic boundary conditions, matrix A is degenerate. Its eigenvalues are are the harmonic frequencies of the linear oscillator chain. The normalized eigenvectors are u 0 n = 1/ √ N, for k = 0, u k n = √ 2/N sin(2πkn/N + γ), for k = 1, . . . , [(N − 1)/2], u N/2 n = (−1) n / √ N, for k = N/2 (N even), u N−k n = √ 2/N sin(2πkn/N + γ + π/2), u N−k n = √ 2/N sin(2πkn/N + γ + π/2), for k = N/2 + 1, . . . , N − 1. Here, • denotes the integer part, the index n runs from 1 to N and the free parameter γ can be fixed arbitrarily, due to the degeneracy of A. The NxN matrix S, which maps Fourier modes Q 0 , . . . , Q N−1 into q 1 , . . . , q N , is orthogonal, and its columns are given by the eigenvectors of A. The momenta in the two dual spaces, P 0 , . . . , P N−1 and p 1 , . . . , p N , are related by the same linear transformation S. Although the choice of γ does not affect the quadratic part of the Hamiltonian in Fourier modes, it does have consequences on the expression of the nonlinear interaction terms among the Fourier modes. A particularly convenient choice is γ = π/4, for which all columns of the transformation matrix S take the same expression 1/ √ N[sin(2πkn/N) + cos(2πkn/N)], with n = 1, . . . , N as the row index and k = 0, . . . , N − 1 as the column index, respectively. Using this transformation, the nonlinear cubic term of the Hamiltonian is transformed into where U γ n,k = cos The sum over n in Equation (5) can be explicitly performed using the properties of the U s (see Appendix A). The result is where and ∆ r = (−1) m if r = mN and 0 otherwise, with m an integer. The B k 1 ,k 2 ,k 3 coupling coefficients incorporate momentum conservation. The same procedure can be used to derive the quartic term in the Hamiltonian, obtaining where again reflecting momentum conservation. This method is fully general and can be applied to derive the Hamiltonian in Fourier modes Q for all values of N.
The first nontrivial value of N is N = 3, since the N = 2 case, being a two-body problem, is totally integrable. For N = 3, we obtain the following potential: where, for convenience, the Hamiltonian is rewritten in slightly different variables where Q 1 is swapped with Q 2 and, then, the sign of Q 1 is changed, leaving the Hamiltonian invariant. Therefore, we can represent the dynamics in Fourier space as the motion of a particle in the two-dimensional potential (11). Let us first discuss the potential (11) with α = 0 and β = 0. The kinetic term of the Hamiltonian is and hence, the center of mass motion is free, (P 0 = const.). The potential is invariant under rotations in the (Q 1 , Q 2 ) plane; therefore, the pseudo-angular momentum is conserved. Going back to the original coordinates with an inverse Fourier transform, it can be shown that which was known to be a constant of motion, see [10]. In Fourier coordinates, this unexpected constant of motion acquires an explicit physical meaning, by associating it to a rotational symmetry. Thinking of the three coordinates (q 1 , q 2 , q 3 ) as positions of a particle in three dimensions, this symmetry corresponds to a cylindrical rotational symmetry of the potential around the axis (u 0 1 , u 0 2 , u 0 3 ) = (1, 1, 1)/ √ 3, the first eigenvector of A. This symmetry is related to the component of the angular momentum, L 0 , along this axis. The two other constants of motion for this three-dimensional phase space are energy and the momentum along the axis P 0 = (p 1 + p 2 + p 3 )/ √ 3. The case α = 0, β = 0 is also integrable. This can be shown by introducing the variables s = Q 1 + Q 2 and d = Q 1 − Q 2 . In terms of these variables, the potential is written as V(s, d) = (3/4)(s 2 + d 2 ) + (α/2)s 3 and is therefore separable along the directions s and d. Along d, the motion is harmonic, while along s, the motion is the one of a cubic oscillator. The motion in the full space is therefore a composition of free, harmonic and cubic anharmonic motions along the three directions.
The intermediate case where both α = 0 and β = 0 is non-integrable. A typical potential surface is displayed in Figure 1: it shows two valleys of different depths. If the motion is within each valley, at low energy, it is quasi-regular and very weakly chaotic, close to integrable. On the contrary, when the motion runs across the two valleys, at higher energies, it is strongly chaotic, displaying intersections on the Poincaré plane that cover an area, see Figures 2 and 3 with two different section planes (Q 1 , P 1 ) and (Q 1 , Q 2 ) (the numerical integration algorithm is described in Appendix B). More precisely, the potential minima have coordinates and the corresponding values of the potential energy are Two valleys exist for weak to moderate quartic interactions (β < 3α 2 ). The boundary value of the potential separating motions within the valleys and above them is given by Notice that in the case of strong quartic interactions (β > 3α 2 ), only one minimum occurs, and the motion is again quasi-integrable.

Figure 2.
Poincaré maps (Q 1 (t * ), P 1 (t * )) for the particle in the potential (11) at fixed times t * defined by P 2 (t * ) = 0. The particle is initially at rest, and in the crater of the upper potential valley, i.e.,  Summarizing, this analysis shows that the FPUT three-particle model is integrable in the two limits α → 0, β = 0 and α = 0, β → 0. Moreover, the model displays a strongly chaotic behavior for intermediate values and α and β and for energies that are close to a saddle in the two-dimensional phase-space of a fictitious particle, which represents the motion in Fourier space. These results call for an analysis of the quantum behavior of such a model, which is pursued in the following section.

The Quantum Three-Particle FPUT Model
In the remainder of the paper, we extend the study to a quantum mechanical version of the three-particle FPUT model. To this end, we promote classical coordinates p j , q j to quantum operatorsp j ,q j , which satisfy canonical commutation relations [q i ,p j ] = iδ i,j (h = 1). The quantum Hamiltonian of the three-particle FPUT model is then given by Analogously to the classical part, we can reduce the degrees of freedom by defining normal modesQ 0 ,Q 1 andQ 2 . The modeQ 0 , corresponding to center of mass motion, decouples from the other two, and we are thus left with the problem of obtaining the eigenvalues and eigenstates for the relevant normal modesQ 1 andQ 2 of the potential

The Integrable Cases
As for the classical FPUT model, if α or β are equal to zero, the model is integrable. In particular, if α = 0, as already discussed in Section 2, the model has a manifest cylindrical symmetry with respect to the Q 0 axis. Therefore, it is convenient to pass to cylindrical coordinates (r, φ, z) defined as In terms of these new variables, the time-independent Schrödinger equation reads where we introduced the wave function ψ = ψ(r, φ, z). Moreover, the Hamiltonian commutes with the operators corresponding to the two integrals of motionP 0 and L 0 =Q 1P2 −Q 2P1 . As a consequence, the energy levels can be labeled by three quantum numbers E = E n,l 0 ,p 0 corresponding to the three commuting observablesĤ,P 0 ,L 0 , and we can search for the solutions of the Schrödinger equation having the separable form The normalized eigenfunctions of theP 0 andL 0 operators factorize, and f n,l 0 (r) is the radial wave function which remains to be found. L is the size of the system along the axial direction z. We notice that square modulus of the total wave function must be normalized to unity by requiring that drdφdzr 2 ψ n,p 0 ,l 0 (r, φ, z) 2 = drr 2 f n,l 0 (r) 2 = 1 .
It is then convenient to introduce the new radial function u n,l 0 (r) = √ r f n,l 0 (r). Finally, inserting this ansatz in Equation (19), we obtain the radial equation − 1 2 d 2 u n,l 0 (r) dr 2 + V eff (r)u n,l 0 (r) = E n,l 0 u n,l 0 (r) , where E n,l 0 = E n,l 0 ,p 0 − p 2 0 /2, and we introduced the effective potential which, in addition to the interaction potential V(r) = (3/2)r 2 + (9/8)βr 4 , also contains a centrifugal barrier term. Summarizing, we mapped the problem into a single particle moving in a one-dimensional potential. We are not aware of the exact analytical solutions of Equation (22) with given boundary conditions, but the equation could be of course solved numerically to determine the radial spectrum quantitatively. This concludes the proof of the integrability of the quantum version of the three-particle β-FPUT model. The quantum integrability of the α-FPUT model is even simpler because the corresponding Schrödinger equation is separable in two independent one-dimensional Schrödinger equations, describing a harmonic and an anharmonic oscillator, respectively.

The Non-Integrable (α+β)-FPUT Three-Particle Model
The rotational symmetry of potential (11) for α = 0, leading to integrability, is lost when the general case α = 0 is considered. In the classical case, we hence found chaotic behavior, as shown by the Poincarè sections. Quantum mechanically, we should expect a situation similar to the one studied in Ref. [18], where the Poisson level spacing statistics was found in the integrable limit, while these statistics were seen to belong to the GOE universality class in the case of the mixed regular-chaotic phase space. In the quantum FPUT model, we should therefore find the GOE statistics at intermediate energies and Poisson statistics at both low and high energy.
We adopt a numerical diagonalization procedure of Hamiltonian Equation (16) [25]. We discretize the Hamiltonian on a two-dimensional square lattice with grid size [−L/2, L/2] 2 and N e grid points assuming open boundary conditions. The grid is carefully chosen such that the minima of the potential are well resolved, and the boundaries of the grid are completely dominated by the quartic term. Rescaled numerical spectra are presented in Figure 4, with ascending order of energy levels. Eigenvalues stabilize as N e is increased, and we checked that the uncertainty in the numerical value of each eigenvalue is small compared to the mean level spacing ∆. The central object in our analysis is the distribution function which is obtained from consecutive level spacings Here, ∆ denotes the mean level spacing. The spectral statistics we present are in line with the Gaussian orthogonal ensemble (GOE) of random matrices with distribution function [17,26,27] This poses a clear signature of a chaotic quantum system with real, normally distributed entries of the Hamiltonian. In contrast, integrability will be associated with an absence of level repulsion and a Poissonian level distribution function Another measure for chaos is given by the level-statistics ratio Figure 5a shows the level distribution function P(s) for an energy window E ∈ [1.9, 2.2] with n = 300 levels. With apparent level repulsion in place, level spacings are distributed according to the GOE prediction, indicating chaotic behavior. This is contrasted with a distribution obtained in an energy window E ∈ [7, 9] (n = 300 levels) in Figure 5b, with the Poissonian level statistic. The two discussed cases confirm the physical picture presented above, though with noticeable fluctuations due to the number of considered levels. The level spacing statistics in different regions of the spectrum is strongly affected by the shape of the potential. We will consider representative parameter values of α and β for which the potential landscape has two non-equivalent minima separated by a barrier of height V b , see Figure 6. At very low energies, the wave function localizes in one of the potential minima, which gives rise to harmonic levels E n ω(n + 1 2 ). At very high energies, the quartic term dominates over the cubic, effectively restoring rotational symmetry, and therefore integrability. From a semi-classical point of view, we expect most of the chaos to appear near the separatrix for levels at energies E V b .
The transition from Poissonian to GOE statistics as a function of energy is traced via the level ratio r, see Figure 7. For this, we sweep the energy window with n eigenvalues centered around an energy E. At high energies (E V b ), we find a level ratio consistent with the Poisson distribution, r ≈ 0.38. At energies near to the barrier height but sufficiently high so that it is not strongly affected by the presence of the potential minima (V b ≤ E ≤ 20V b ), the level ratio peaks at the GOE value r ≈ 0.52, signaling a transition to a fully quantum chaotic behavior of the system. The transition between these two extremal cases occurs on the same energy scale on which rotational symmetry is broken, and hence appears visible in Figure 7. For low energies, the r-value decreases again but not fully to the Poissonian value since the minima support harmonic oscillator eigenfunctions. This behavior should be compared with the one of the integrable model (α = 0) shown in Figure 7b. In this case, the level ratio is never compatible with chaotic statistics independently from the energy interval considered.
Real-space wave functions support the picture of chaos in distinct regions of the spectra, see Figure 6a-c. At low energies (E < V b ), back to the picture of quasi-integrability, wave functions in the potential minimum attain the shape of two-dimensional orbitals which can be labeled by radial and angular quantum numbers, Figure 6a. In contrast, the chaotic behavior for states with intermediate energies (V b < E < 20V b ) is reflected in irregularly shaped wavefunctions with random direction of the nodal curves, Figure 6b. At high energy (E V b ), close to the other integrable limit, wavefunctions have a structure similar to the one observed at low energy, Figure 6c. In fact, the nodal curves strongly repel each other, and do not cross, up to our numerical resolution, see Figure 8b. This was predicted by Berry [14] and first verified for a free particle in a stadium [28].

Eigenstate Thermalization Hypothesis
According to the eigenstate thermalization hypothesis (ETH) of Deutsch [19] and Srednicki [20] (see Ref. [29] for a review), if a quantum system is chaotic, the average of a generic observable taken on the eigenstates of the energy is a smooth function of the energy eigenvalue itself. In connection with ETH, and focusing on the relation between quantum and classical systems (e.g., for the FPUT model), the adiabatic gauge potential is introduced to characterize suppressed heating by a slow drive [30][31][32]. We checked the ETH property, restricting to the stationary solutions of the Schrödinger equation, without any attempt to analyze the time evolution. Indeed, it is also guessed that the relaxation time to the average is controlled by the size of the system, i.e., the number of particles. Therefore, we cannot expect a fast relaxation in such a low-dimensional phase space as the one of the three-particle FPUT model. The averages on the energy eigenstates of kinetic energy and angular momentum are shown in Figures 9 and 10 as a function of energy in three different regions of the spectrum. Both at low and at high energy, the kinetic energy and angular momentum are irregularly scattered in a wide region of values, and change abruptly as the energy is varied. It is only in the central chaotic region of the spectrum that they appear to fluctuate stochastically around an average. We believe that this result is a particularly transparent numerical verification of the ETH hypothesis.

Conclusions
In this paper, we investigated the classical and quantum Fermi-Pasta-Ulam-Tsingou model with three particles. Classically, the model is integrable in the two limits where the cubic or the quartic term of the Hamiltonian vanishes. In the intermediate region, the model shows a mixed phase space with both regular and chaotic trajectories, in agreement with the picture suggested by the Kolmogorov-Arnold-Moser theorem. We computed the energy spectrum and showed that, as expected from the Bohigas-Giannoni-Schmit universality hypothesis, the level spacing satisfies the Gaussian orthogonal ensemble statistics in the chaotic part of the spectrum. In the two integrable limits, the statistics become Poissonian. We also showed how two generic observables obey the eigenstate thermalization hypothesis of Deutsch and Srednicki in the chaotic part of the spectrum.

Institutional Review Board Statement: Not applicable.
Data Availability Statement: Data available upon reasonable request.