Chaotic Dynamics in a Quantum Fermi–Pasta–Ulam Problem

We investigate the emergence of chaotic dynamics in a quantum Fermi—Pasta—Ulam problem for anharmonic vibrations in atomic chains applying semi-quantitative analysis of resonant interactions complemented by exact diagonalization numerical studies. The crossover energy separating chaotic high energy phase and localized (integrable) low energy phase is estimated. It decreases inversely proportionally to the number of atoms until approaching the quantum regime where this dependence saturates. The chaotic behavior appears at lower energies in systems with free or fixed ends boundary conditions compared to periodic systems. The applications of the theory to realistic molecules are discussed.


Introduction
Understanding vibrational energy flow in molecules is one of the challenges in modern science and technology [1,2]. Vibrational energy flows control energetics of chemical reactions, determine heat balance in modern nano-devices [1,[3][4][5][6] and can be manipulated similarly to electrons and photons and used to carry and process quantum information [7][8][9]. Intramolecular energy relaxation and transport are dramatically sensitive to the molecule's ability to attain the thermal equilibrium [4,10].
After seminal work by Stewart and McDonalds [11], it has been realized that the internal vibrational relaxation can be absent or proceed very slowly in small enough molecules and/or at low temperature. Based on these observations, the concept of localization of low energy anharmonic vibrational states of poly-atomic molecules within the manifold of harmonic product states of almost independent normal modes was put forward by Logan and Wolynes [12]. In earlier [13][14][15] and later [16][17][18] work, similar ideas have been developed for particle and spin systems. Theory was further extended combining random matrix theory methods [19][20][21] and Bose Statistics Triangle Rule approach [22][23][24] and this extension was reasonably consistent with the experimental observations [11].
This development is qualitatively consistent with the investigations of the classical counterpart problem of anharmonic vibrational dynamics. Its simplest realization in atomic chains probed as a modeling system for irreversible dynamics was considered in the celebrated work by Fermi, Pasta and Ulam [25] (FPU), where the quasi-periodic behavior has been discovered for the evolution of the initial excitation instead of irreversible energy equipartition. Despite over sixty years of investigations of the FPU problem, its complete understanding remains a challenge [26][27][28][29].
Both quantum and classical non-linear vibrational dynamics can be characterized by a critical energy separating low energy integrable (localized) and high energy chaotic (delocalized) behaviors. In the chaotic regime, each part of the system can be thermalized due to its interaction with the rest

Model
The FPU model of anharmonic atomic chain with different common boundary conditions including periodic, fixed ends and free ends (see Figure 1) can be described by the Hamiltonians defined as 24 , periodic, , fixed ends, 4 24 , free ends. (1) Below, we set mass, harmonic force constant and Planck constant to unity,h = M = k = 1. Force constants A and B describe relative strengths of third-and fourth-order anharmonic interactions. The fixed ends problem has been studied in the classical FPU paper [25].
Anharmonic interactions should be weak for the system energy E of interest to justify the applicability of the series expansion for non-linear terms. Assuming approximate energy equipartition, one can estimate (x i − x i+1 ) 2 ∼ E/N, which leads to anharmonic interaction estimates V 3 ∼ AE 3/2 / √ N and V 4 ∼ BE 2 /N for the third-and fourth-order anharmonic interactions, respectively. Comparing harmonic and anharmonic interactions, we end with the restrictions for energy density in the form However, we impose a stronger constraint on the energy requiring the stability with respect to the dissociation. Consequently, the total energy should be less than the dissociation energy, E d . This energy can be estimated for the single bond assuming that anharmonic energy becomes comparable to the harmonic one which suggests E d ∼ 1 A 2 ∼ 1 B . For instance, for the Morse potential [41] often used to model atomic interactions, one has We assume that the system energy is always smaller than the dissociation energy so the molecule is stable with respect to large coordinate displacements. Equation (4) can be satisfied for a quantum system only if it is satisfied at least for the minimum energy that can be estimated as a quantization energy E ∼ 1. Consequently, the anharmonic interactions should be weak, which requires These requirements are well satisfied in real molecules because of the small amplitudes of vibrations of heavy atoms. For instance, using the Morse potential for C-C bond, one can express the dimensionless parameters A and B as The conditions of the weakness of anharmonicity are better satisfied since the dissociation energy contains the additional large numerical factor (see Equation (3)).
These frequencies are identical to normal mode quantization energies (remember that we set the Planck constanth to unity). The mode with zero quasi-momentum p = 0 can be excluded from the consideration because it corresponds to identical displacements of all atoms that cannot modify a system energy. Therefore, only N − 1 normal modes are significant and the harmonic part of the Hamiltonian can be re-expressed in the diagonal form with respect to these normal modes as where operators b † p and b p describe creation or annihilation of one quantum of vibration of normal mode (phonon) p. The harmonic problem is obviously integrable since the system breaks into N independent oscillators (phonons) and each phonon population number operator ν p = b † p b p represents a local integral of motion [42] in the momentum representation. Each many-body eigenstate |S > of the harmonic system can then be represented by an arbitrary sequence of integer population numbers S = {ν p }.
Anharmonic interaction mixes up these states because it breaks down the conservation of individual phonon population numbers. In the periodic system with only fourth-order anharmonicity (A = 0 in Equation (1)), this interaction can be expressed as The factor ∆(p 1 + p 2 + p 3 + p 4 ) is equal to unity if the sum of all four momenta is equal to zero or integer fraction of N (due to Unklamp processes); otherwise it is equal to zero giving rise to a quasi-momentum conservation.
Because of the above conservation law, the basis states of the system with given normal mode population numbers can be split into N subsystems with the total quasi-momenta Q = 0, 1, −1, ..., N/2 for even N or Q = 0, 1, −1, ..., −(N − 1)/2 for odd N determined with the accuracy to the addition of integer number of Ns. Each subsystem should be studied separately since the states from different subsystems do not interact with each other. In addition, the states with Q = 0 and Q = N/2 for even N possess a mirror reflection symmetry with respect to replacement all states S = ν p with the states S − = ν −p . Then, the states with Q = 0, N/2 can be split into two subgroups symmetric or antisymmetric with respect to the mirror reflection symmetry. Consequently, all many-body states can be split into N + 1 subgroups for odd N and N + 2 subgroups for even N that should be considered separately.
Similarly, one can consider the interacting normal modes for free and fixed boundary atoms (see Equation (1)). One can similarly introduce normal modes for this problem and their anharmonic coupling. In these two cases, one cannot introduce quasi-momenta because of the lack of translational symmetry. However, there is a mirror reflection symmetry with respect to the middle of the chain. Then, all states can be separated into two subgroups of symmetric and anti-symmetric states with respect to that symmetry. The states belonging to the different subgroups can be considered separately.

Localization-Delocalization Transition: Qualitative Analytical Consideration
There are over sixty years of history of the investigation of chaos in the classical FPU problem and this problem remains a challenge [26,27]. The situation with the quantum mechanical problem is even more complicated [43,44]. Below, we summarize the established results for the system of a few atoms [34] and attempt to extend them to atomic chains having many atoms N 1 first for the β FPU problem containing only fourth-order anharmonic interactions and then extend it to the mixed α + β problem.

Localization-Chaos Transition in the Small
Localization-chaos transition in systems with the small number of atoms N ∼ 1 can be well understood following Ref. [34]. In this paper, the critical energy E c separating localized and chaotic states is estimated using a dimensionality arguments since the only value having the dimension of energy can be constructed using the parameters in Equation (1) in the form The numerical studies of both classical and quantum mechanical problems in Ref. [34] have confirmed these expectations provided that the system is semiclassical, i.e. the system energy E c is much larger than the quantization energyhω ∼ 1. This requires B 1. Transitions in quantum and classical systems occur under almost identical conditions since the system is semiclassical because the maximum quantization energy max(ω p ) is of order of unity (Equation (7)). Consequently, it is much less than the energy per the mode E c /N expressing the thermal energy k B T (see Equation (10), remember that B ≤ 1 and N ∼ 1).
Similarly, one can estimate the critical energy for the α FPU problem with the third-order anharmonic interaction as Equations (10) and (11) differ from the expectations of the analysis exploiting resonances for many-body transitions that has been successfully applied to problems of interacting spins [15,[45][46][47] or electrons [16,48]. According to this criterion, chaos emerges in the presence of approximately one resonance per the many-body state under the condition that the diagonal interaction of resonant modes is larger or comparable to their resonant coupling [15,[47][48][49] that is needed to avoid destructive interference between consecutive resonant transitions. Such interaction is present naturally for the fourth-order anharmonicity in Equation (9) (for instance, the terms with p 1 = p 2 and p 3 = p 4 are diagonal in the phonon product state representation). There is no such interaction in the case of the third-order anharmonic interactions (α FPU problem), which changes the definition of delocalization transition as it is discussed in Section 3.3.
However, in the β FPU problem under consideration, the matrix elements M of the four phonon interactions in Equation (9) grow proportionally to the squared population numbers M ∼ Bν 2 p ∼ BE 2 for the system energy E exceeding the quantization energy. The typical energy change in a four phonon process is of order of their quantization energy that is of order of unity. Consequently, the amount of resonances approaches unity at E ∼ B −1/2 in contrast to Equation (10).
This conflict can be resolved at the qualitative level modifying the definition of the resonances in accordance with Ref. [50] where a single-particle localization problem has been considered for harmonically coupled vibrations. For example, two unit mass oscillators with frequencies ω a and ω b coupled by the interaction k ab u a u b are in resonance under the condition |ω 2 a − ω 2 b | < k ab , while in terms of matrix elements the resonance takes place at |ω 2 a − ω 2 b | < k ab (ν a + 1)(ν b + 1). To define the resonance correctly one can consider the energy change not for a single resonant transition but for the whole set of possible transitions involving these four phonons, which will increase the typical energy change due to the transition (hω ∼ 1) by the factor of a typical phonon population number ν p . Then, the resonance criterion can be written as Bν 2 p ∼ ν p . Setting δE ∼ ν p ≈ E c we end up with Equation (10). The problem of interest with large number of atoms needs a special consideration given in the next section.

Classical Regime
Consider the localization-chaos transition in the case of a large number of atoms assuming that the system is semiclassical (population of each vibrational mode exceeds unity). One can still consider resonant interactions similarly to the previous section. In the periodic system of N atoms, one can find N 3 possible four phonon processes for a typical state (the fourth phonon mode is fixed by the quasi-momentum conservation law in Equation (9)). Consequently, the minimum energy difference between two modes coupled by the fourth-order anharmonic interaction is given by δE ∼ ν p /N 3 ∼ E/N 4 (factor ν p is added similarly to Section 3.1). The interaction matrix element scales as M ∼ BE 2 /N 3 . Here, the factor 1/N comes from the definition of anharmonic interactions in Equation (9) and the factor ν 2 2 while the thermal energy k B T for N classical oscillators is given by E/N [51]. Setting δE ∼ M to ensure the presence of resonant interactions we estimate the localization threshold as Similar dependence can be obtained for the atomic chains with fixed or free ends boundary conditions where there is no quasi-momentum conservation. In those systems, one has N 4 possible four phonon transitions and 1/N 2 (instead of 1/N, Equation (9)) scaling of anharmonic interaction matrix element. Then, extra factors N −1 are canceled out on both sides of the criterion of resonance leading to Equation (12).
It is noticeable that the estimated behavior of localization threshold in Equation (12) is qualitatively consistent with the earlier estimates [39,40,52,53] obtained using the stability analysis of the classical dynamics of a non-linear FPU chain in the form Since this equation agrees with numerical studies in Section 4 for free and fixed ends boundary conditions, we used it for quantitative estimates.
In our qualitative analysis of resonant interactions, we considered only typical phonons with energy close to unity, while the low frequency phonons were ignored. Based on the present understanding of localization-chaos transition, it is hard to expect that they can suppress the chaotic dynamics because the typical phonons form the ergodic spot normally capable to equilibrate the rest of the system [48,54]. It is hard to expect that they can give additional support to the chaotic dynamics since they are coupled weakly to the rest of the system compared to typical phonons.
On the other hand, one can imagine marginal states with the only low frequency phonons being excited. These states can possibly show anomalously strong localization behavior as predicted for the classical systems in Refs. [38,55]. There are other suggestions for classical systems [27][28][29]56] that the crossover in Equation (13) does not describe the transition to a truly integrable (localized) behavior but separates strongly ergodic and weakly ergodic regimes at high and low energies, respectively. Since the numerical simulations in Section 4 show the pure localization transition, we did not see any evidence for such behavior in a quantum regime. It cannot be excluded that at larger number of atoms some additional channels for chaotic behavior can emerge.
The criterion in Equation (12) is valid until the system remains semiclassical, meaning that the phonon population numbers exceed unity. This requires the thermal energy E c (N)/N to exceed the quantization energy, which is of order of unity. Thus, the classical regime takes place at sufficiently small number of atoms The crossover energy E c expresses the minimum threshold energy in the classical regime of vibrations. At larger N, the system should be treated quantum mechanically as considered in Section 3.2.2.

Quantum Mechanical Regime
We begin the consideration with the analysis of the problem in terms of resonant interactions. Imagine that the system energy is spread between phonons of energy such that √ E/N < ≤ 1. In our case of small energy E < N, the thermal energy is given by k B T ≈ √ E/N and the lower limit for the energy qualitatively represents the typical thermodynamic equilibrium.
For an arbitrary energy , the total number of phonons is given by n ∼ E/ . This number is smaller than the number of quantum states with energy of order of that is given by N . The modeling system is non-degenerate so typical populations of vibrational states do not exceed unity and one can describe the emergence of chaos requiring a single resonant interaction per a many-body quantum state (cf. [47,48]). The typical anharmonic interaction strength for periodic boundary conditions scales as M ∼ B 2 /N. The energy difference to the adjacent state coupled to the given state and having the same number of phonons can be estimated as /N c where N c is the number of anharmonically coupled states with the same number of phonons. This number can be estimated considering the number of possible anharmonic transitions including n 2 possible double annihilations of phonons and N creations (the fourth phonon is fixed by the quasi-momentum conservation law and we consider only processes conserving the number of phonons). The resonant interactions exist under the condition B 2 /N < /(n 2 N). Then, the critical energy E c = n can be estimated as The generalization to the non-periodic boundary conditions can be done similarly to that in the previous section.
This answer is universal and insensitive to the number of atoms. It predicts the saturation of the dependence of critical energy on the number of atoms in the quantum regime. Based on the numerical results in Section 4, we assume that the saturation takes place at E c = N. Then, combining Equation (13) with Equation (15), one can write the summary of the predicted behaviors as In Section 4, it is verified for the minimum division of energy E into phonons with energies of order of 1. The more accurate numerical analysis of the problem is postponed for the future.
The consideration ignores correlations between phonon energies and momenta, that can take place due to quasi-momentum conservation in a periodic system [44] or some trace of its conservation in the system with fixed end boundary conditions. These processes are fully suppressed for free ends boundary conditions where the above consideration is most applicable. It is less applicable for the periodic system where these correlations can be significant. For very small system energies comparable to the maximum quantization energy 1, the periodic system becomes integrable [44] so the consideration fails. We still believe that our consideration is valid even for a periodic system where the hot spot [54] can be formed by several excited phonons with nearly maximum energy. These phonons, indeed, form chaotic state (see Section 4) and can equilibrate other parts of the system. The accurate numerical verification should resolve the raised questions.

α FPU Problem
Here, we consider the effect of the third-order anharmonic interaction on the state of the system. Let us begin the consideration with the classical regime, E > N. For the small number of atoms, the dimension based arguments lead to the estimate A c ∼ 1/ √ E (Equation (11)). For a large number of atoms N in the periodic system, one can find N 2 possible three phonon processes so the minimum energy shift can be estimated as δE ∼ ν p /N 2 ∼ E/N 3 . Remember that the third phonon state is fixed by the quasi-momentum conservation law in Equation (9). The interaction matrix element scales as M ∼ AE 3/2 /N 2 . Consequently, there are resonant interactions in the case of sufficiently large energy This estimate is consistent with Ref. [53]; however, we do not think it describes the localization breakdown correctly because of the lack of the diagonal interaction. In this case, resonant transitions are independent of each other [47,48], which prevents the system from delocalization similar to the XY model, where there is no diagonal interaction [57]. Following Ref. [57], one can consider the induced resonant interaction in higher orders anharmonicity following the Schrieffer and Wolff method [58]. In the first non-vanishing order, the fourth-order anharmonic interaction will be generated. This generated interaction is similar to the one in the β FPU problem with the effective interaction constant B * ∼ A 2 if expressed in the momentum space. However, the induced diagonal interaction is much less than the third-order resonant interaction because A 1, thus it cannot enhance delocalization due to three phonon transitions. Following Refs. [17,18], one can suggest the weaker delocalization criterion of one resonance per each normal mode. This leads to the criterion E c3 ∼ 1/A 2 that is insensitive to the number of atoms.
However, the more efficient delocalization should take place due to the induced fourth-order interaction characterized by the interaction strength B * ∼ A 2 . In that case, one can expect the chaotic behavior following the estimate of Equation (10) that reads Similar to Section 3.2.2, this criterion is valid in the classical regime realized at N < 1/A while in the opposite regime this dependence saturates at Following Ref. [28], one can expect that this prediction should be valid to the same extent as Equation (10). Indeed, if one considers the combined α + β problem containing both third-and fourth-order anharmonic interactions then the chaotic state formation is dramatically suppressed at B = 4A 2 /9 because under these conditions the non-linear interaction would be identical to power series expansion of the integrable Toda model [59]. Consequently, in this regime, one can expect that the third-order anharmonic interaction characterized by the constant A should produce similar delocalization effect to the fourth-order problem, characterized by the interaction constant B ∼ A 2 in a full accord with the estimate of Equation (18).
Following the recipes of Ref. [28], one can extend the above consideration to the general α + β problem, which can be reduced to the β FPU problem with the interaction constant B * defined as Consequently, one can predict the localization threshold energy for the α + β problem in the form of generalized Equation (16) This result is not applicable if the denominator in Equation (21) is very close to zero. In the case of nearly zero denominator, the problem can be effectively described by the sixth-order anharmonic interaction with the interaction constant C = B 2 [28]. In this regime, the similar analysis of resonant interactions can be applied leading to the threshold energy behaviors E c ∼ 1/(B √ N) in the classical regime and E c ∼ 1/B 2/3 in the quantum regime, where N > B −2/3 .
Equation (21) is the main result of the present work. In the next section, some numerical justification is given based on the diagonalization of Hamiltonians in Equation (1) within the reduced basis of many-body states.

Numerical Analysis of the Transition Localization-Chaos
The numerical analysis is limited to the β FPU problem to avoid overcomplexity. Below, the numerical studies attempting to justify the analytical predictions of Section 3 are reported. In Section 4.1, we define the numerical criterion of the chaotic behavior. Since the basis of many-body states is infinitely large, one should restrict the phase space. In Section 4.2, we introduce the method of basis restriction considering the states with the fixed number of phonons. In Section 4.3, we investigate the dependence of the localization threshold on the system energy (number of phonons) and number of atoms.

Level Statistics
The chaotic and integrable (or localized) phase of quantum systems can be identified using the statistics of energy levels. It is expected that in the chaotic phase all states substantially overlap with each other which leads to their energy level repulsion and, consequently, Wigner-Dyson level statistics [60,61] suggesting zero probability density for nearest eigenstates energy difference approaching zero. In the localized phase the overlaps of a majority of states are negligible so their energies are independent, which results in the Poisson statistics for energy level differences. In numerical studies exploiting exact diagonalization of the system Hamiltonian, the energy level statistics can be probed directly and used to identify the state of the system.
Other methods including the analysis of correlation functions [46,62], entanglement entropy [63] or local integrals of motion [42] can also be used to study the delocalization with respect to the specific basis. However, the results depend on the choice of the basis. For instance, the basis of single particle states can be defined in the coordinate or momentum representations and localization in the coordinate space suggests delocalization in the momentum space and vice versa. Eigenstates of the FPU problem at very low energies [44] are delocalized in the basis of product states composed by independent phonon states, while the problem remains integrable [43]. The level statistics based definition is basis independent and therefore it seems to be the most objective criterion to distinguish localized and chaotic phases.
The level statistics have been characterized using the averaged ratio of successive gaps, < r >, defined as [61] < r >= min(δ n , δ n+1 ) max(δ n , δ n+1 ) , where δ n = E n+1 − E n is the energy difference of adjacent energy levels of the system, Equation (1), obtained by means of exact diagonalization of the Hamiltonian. According to Ref. [61], in the chaotic regime characterized by Wigner-Dyson statistics, one has < r >≈ 0.5307 for the Gaussian Orthogonal Ensemble of interacting states, while in the case of localization where the Poisson statistics is expected one has < r >≈ 0.3863.
If a system has integrals of motion, which takes place for our system of interest (see Section 2) the states with different values of the related integrals do not repel each other that would lead to the inevitable deviation from the Wigner-Dyson statistics even in the delocalized regime. To avoid this problem, the states should be split into subgroups with a certain value of all integrals of motion [64]. For the periodic system and even number of atoms N, one can introduce N + 2 such subgroups characterized by the quasi-momentum including N − 2 subgroups with quasi-momenta different from 0 or N/2 (Q = −N/2 + 1, ..., −1, 1, ..., N/2 − 1) and four subgroups with Q = 0, N/2 either symmetric or anti-symmetric with respect to the reflection transformation. For an odd N, one has N + 1 subgroups including N − 1 subgroups with quasi-momenta (Q = −(N − 1)/2, −(N − 1)/2 + 1, ..., −1, 1, ..., (N − 1)/2) and two subgroups with Q = 0 either symmetric or anti-symmetric with respect to the reflection transformation. For other boundary conditions, one has only two subgroups either symmetric or anti-symmetric with respect to the center of the chain. Since the results for symmetric and anti-symmetric states are quite similar, the level statistics reported below are related to the symmetric states.
In contrast to spin or particle systems [46,61], we cannot directly apply exact diagonalization method to Hamiltonians in Equation (1) because they have infinite basis of states since the vibration population numbers can take infinite number of values. Therefore, the basis states should be restricted to the finite number of states as described in the following section.

Basic Approximation
For any boundary conditions and specific subgroup of states, the Hamiltonian in Equation (1) cannot be exactly diagonalized since the total number of possible basis states is infinite. To avoid this complexity, the off-diagonal anharmonic interaction is restricted to the terms conserving the total number of excited quanta, n t , similar to Ref. [34]. This means that only terms having two b † and two b operators are taken into consideration. Similar terms are left for other boundary conditions. This approximation should be valid at least qualitatively if the annahrmonic interaction is weak. Consequently, the anharmonic interaction energy Bn 2 t /N should be less than the harmonic interaction energy n t that yields The modified Hamiltonian has a finite basis set for each specific number of atoms N and number of phonons n t so it can be studied using the full diagonalization of the problem. The representative level statistics for the chain of N atoms with free ends boundary conditions, total number of phonons n t = 14 and the strength of anharmonic interaction B = 0.2 is shown in Figure 2. In contrast to the problems of particles or spins placed in a random potential [61] where the ratio parameter r can be averaged over many disorder realizations, here we have the only one realization of the system. In this realization the ratio r itself represents quasi-random number ranging between 0 and 1 in a chaotic manner as shown by the dashed dark blue line in Figure 2. However, averaging the data over 972 adjacent states (5% of the total number of states) leads to the smooth curve clearly approaching chaotic limit of 0.53 near the middle of the spectrum (similar value < r > ∼ 0.53 at the upper edge of the spectrum E ≈ 30 is probably a random coincidence, since it differs for other parameters N, n t and B at spectrum edges). The average ratio < r > for the given set of parameters N, n t , B has been determined taking the arithmetic average of this minimum ratio over the middle half of the system eigenstates, as shown in Figure 2. This procedure describes how the data were collected to analyze the transition between the localized and chaotic regimes as a function of the number of atoms and phonons and the strength of anharmonic interactions. In other calculations, the same averaging of the ratio parameter r was performed.  (22)) including r as it is and average ratio < r > over 972 adjacent eigenstates that clearly tends to the chaotic behavior < r > ≈ 0.53 at the maximum density of states. The average level statistics < r > = 0.5087 for the chain containing 8 atoms, states with the number of phonons n t = 14 and anharmonicity strength B = 0.2 averaged for mean 1/2 of eigenstates located between two vertical lines.
The typical harmonic energy can be estimated for the given number n t of phonons using their sinusoidal dispersion law, Equation (7), as E h = 2n t < |sin(x)| >= 4n t /π. In the case of Figure 2, this energy can be estimated as 17.83. This energy is smaller than the typical average energy by around 10% because of the anharmonic correction to the energy, which is still small.
The chosen representative states having maximum density at fixed number of phonons do not perfectly represent the true thermodynamic states of the system at the given energy. In the classical regime E > N (see Equation (14)), the thermodynamic average number of phonons scales as E ln(N) due to the contribution of low frequency phonons. We believe that this difference is not crucial since the logarithmic factor is related to low frequency phonons, which have substantially reduced anharmonic interaction strength and therefore can be ignored in the consideration of resonant interactions as discussed in Section 3. The other reason is that the investigated states coexists with the "thermodynamic" states at the same energy. If the states under consideration are chaotic, the other states at the same energy should be usually chaotic as well [54].
In the quantum regime, E < N, the representation of the typical configuration by E phonons with typical energy of order of unity is much less relevant than for the classical regime since the typical phonon energy is given by the thermal energy √ E/N that is much less than 1. However, since the delocalization criterion, Equation (15), is universal and does not depend on the number of phonons we also believe that the theory should be applicable to the whole system at least quantitatively.
Thus, the numerical results reported below are preliminary and need improvement that is postponed for the future.
The validity of the approach has been checked for the classical regime extending the basis to all states with the number of phonons less or equal to n t . The results for this extension are consistent with those for the phonon number just equal to n t . However, the calculations are much faster in the latter case and they permit us to obtain more conclusive results. The approach that seems to be more "natural" restricting the basis to the states with energies less than a certain maximum energy E max works much worse and requires E max ∼ 2E to give a reasonable estimate for the level statistics at energy E, which substantially limits our abilities to obtain conclusive results. This could be the consequence of broken connections due to the exclusion of significant states naturally present in the theory conserving the number of phonons.

Dependence of Localization Transition on the Boundary Conditions and the Numbers of Phonons and Atoms
Since our results below for the level statistics (< r >, see Equation (22)) are expressed as a function of anharmonic interaction strength B (remember that only β FPU problem is considered numerically), it is convenient to re-express the criterion in Equation (13) in terms of the critical strength B c dependence on the number of atoms N and phonons n t . Using Equation (14) for the classical and quantum regimes, we get.

Effect of Boundary Conditions
To examine the effect of boundary conditions, we consider some representative data obtained for the level statistics parameter < r > vs. the strength of anharmonic interaction B following the technique described in Section 4.2 (see Figure 2) for the chain of N = 10 atoms with all possible boundary conditions and for quasi-momenta Q = 0, 1 and 2 in the case of periodic conditions. These dependencies are shown in Figure 3. According to Figure 3, it is clear that at large anharmonicity B > 0.4 the system is chaotic, while it is integrable at small anharmonicity B < 0.15. The Poisson statistics limit < r > ≈ 0.38 is not reached even at small B < 0.4 possibly because the system under consideration is not random in contrast to electronic systems [60] and some correlations between energies are significant in integrable states.
For N = 10 and n c = 10 the criterion of Equation (24) predicts B c ≈ 0.15. This number is in an excellent agreement with the numerical results for free ends or fixed ends boundary conditions as indicated by the vertical line in Figure 3. Assuming that B c ≈ 0.15 characterizes the transition in the case of fixed or free ends boundary conditions, one can estimate the threshold anharmonicity in the periodic system using similar criterion (see dashed vertical line in Figure 3) as B c,per ∼ 0.28. This estimates the approximate difference between two critical anharmonicities as B c,per /B c ≈ 1.87 and we use this ratio for quantitative estimates of transition parameters below. Probably, theoretical analysis of Refs. [39,52] is less applicable to the periodic system because of additional integrals of motion there. The emergence of chaotic phase in periodic system at larger anharmonicity compared to other boundary conditions can be the consequence of the smaller effective phase space in the former case due to additional integrals of motion lacking in the latter case. On the other hand, the results for periodic conditions are almost insensitive to the quasi-momentum, and the results for fixed and free ends boundary conditions are also quite similar to each other. Therefore, in our predictions for localization threshold, we do not distinguish between the free and fixed end boundary conditions as well as between different quasi-momenta for periodic boundary conditions. However, it is necessary to distinguish between periodic boundary conditions and others. We suggest the simplest form of difference redefining Equation (24) reasonably valid for fixed and free ends regimes by multiplying its left hand side by the factor 0.28/0.15 in agreement with the numerical results in Figure 3. Then, for the periodic system, Equation (24) should be modified as Consequently, one should modify the critical energy behavior predicted by Equation (21) as

Dependence of Localization Threshold on Numbers of Atoms and Phonons
Consider the dependence of the threshold anharmonicity on the energy expressed through the number of phonons. Most of the data are presented for periodic chains because the large number of integrals of motion there reduces the total number of states permitting us to investigate larger numbers of atoms and phonons compared to other boundary conditions. For a demonstration of the method, we consider periodic chain for N = 10 atoms with possible numbers of phonons n t = 9, 10, 11 and 12. The subgroup of symmetric states with quasi-momentum Q = 0 is considered.
The choice of possible parameters n t is limited because of the poor data averaging for small size of phase space less than 5000 states and exponential increase of the number of states with increasing n t . Indeed, for n t = 9 the basis contains 4420 states that is insufficient for good averaging of level statistics, as shown in Figure 4a, while for n t = 12 the basis contains 26,720 states that is close to the maximum matrix size where exact diagonalization can still be performed using standard MATLAB algorithms.
To determine the algebraic dependence of localization threshold on the number of phonons, we use the data rescaling procedure similarly to the earlier work in spin systems [46,57,64]. This procedure attempts to attain the maximum match between different data rescaling the x axis. As shown in Figure 4b, the reasonable match can be attained rescaling the data for different n t with respect to those for maximum N = 12 by the n t -dependent parameter η shown in Figure 4b. The scaling of parameter η(n t ) is related to that of a critical anharmonicity B c (n t ) as B c (n t ) = B c (n t,max )η(n t ). (27) In our case, n t,max = 12. Consequently, we end up with the dependence of the critical anharmonicity B c on the number of phonons, n t , as shown in Figure 5.   The observed dependence is in between two predictions of Equation (24) that is not surprising because the calculations are made for n t ∼ N = 10 near the crossover between classical and quantum regimes. This justifies our definition of that crossover in Equations (21) and (26). Similar behavior takes place for the same number of atoms and quasi-momentum Q = 1, as shown in Figure 6.
Consider limiting quantum (n t < N) and classical (n t > N) regimes. Representative results for quantum regime are shown in Figure 7 in the case of N = 13 atoms and a number of phonons, n t , ranging from 6 to 10. The observations clearly agree with Equation (24) for n t < N. The opposite, classical limit, n t > N, is represented by the periodic chain of N = 6 atoms considered for the number of phonons, n t , ranging from 18 to 36. The obtained dependence shown in Figure 8 is very close to the inverse proportionality, Equation (24), valid in this limit. The growing deviation at small n t is probably caused by quantum effects, significant for n t ∼ N.
Similarly, one can consider the dependence of the threshold anharmonicity, B c , on the number of atoms N at fixed number of phonons n t . This dependence is expected to be an inverse proportionality in the classical regime of a large number of phonons, n t > N, while no dependence is expected in the opposite, quantum limit of a small number, see Equation (24). These expectations are consistent with the results given in Figure 9 in the quantum (n t = 6, Figure 9a) and classical (n t = 12, Figure 9b) limits.
The results for other boundary conditions are also consistent with theoretical predictions, as illustrated in Figure 10, both in classical (N = 6) and quantum (N = 11) regimes.
Thus, the numerical investigation of localization-chaos transition supports the theory predictions, Equations (21) and (26). A nearly perfect match of rescaled dependencies of average ratio < r > on anharmonic interactions (see Figures 6a-8a) for different numbers of atoms or phonons suggests that the width of the transition is proportional to the transition energy. However, this conclusion is very preliminary because of the relatively narrow domain of numerical studies and it needs further verification.

Discussion
Here, we reformulate the results in terms of standard notations in Table 1, and attempt to apply them to organic molecules. We predicted the threshold energies for emergence of chaotic dynamics for combined α + β FPU problem as a function of anharmonic interaction strengths and system sizes as given by Equations (21) and (16).
For practical application of these results, it is convenient to re-express them in terms of the dimensional force constant k, atomic mass M and Planck constanth. This requires changing the anharmonic interaction constants as B → B/k 2 and A → A 2 /k 3 in classical estimates and modifying the critical energy as N → Nh √ k/M. The results are presented in Table 1 in the standard notations. One can attempt to apply these results to organic molecules using parameters for C-C bond extracted from the Morse potential [41] that can be defined in terms of bond dissociation energy E d = 5.78 × 10 −19 J and inverse interaction radius α = 3.45 × 10 10 m −1 as Consequently, the expressions for the threshold energy in classical and quantum regimes for either free or fixed ends boundary conditions can be written as E c,cl = 4π 2 E d /(3N) and where M is the atomic mass. The transition between two regimes takes place at the number of atoms N c ≈ 2π The chaos can take place in the stable molecular state at energy less than the dissociation energy that is true only for sufficiently long molecules containing N ≥ 14 atoms. This is the result for atomic interactions determined by the Morse potential.
Considering the specific parameters for C-C bond, one can estimate the minimum crossover energy to the chaotic state as E c,q ≈ 0.7E d and transition to the quantum regime is expected at N > N c ≈ 20 atoms.
It is interesting to find how long the chain of carbons should be to attain the chaotic state at room temperature. Since in the Morse potential model room temperature is much smaller than the characteristic quantization energyh √ k/M ∼ 1300 cm −1 , we should use the quantum expression for the energy E ≈ Nπ 2 (k B T) 2 /(h √ k/M). Setting E ∼ E c,q and k B T ∼ 4 × 10 −21 J at room temperature, we get the estimate This number is very large. The energy of molecules studied in the experiments [11] is of order of 3000 cm −1 , which is much lower than the minimum energy E c,q needed to reach the chaotic state; however, some of them show a fast internal relaxation. Therefore, a Morse potential based model of the FPU atomic chain seems to be not quite relevant there. Perhaps this is because real molecules (e.g., alkane chains) are not perfectly linear but have a zigzag shape, making them much softer. In addition, transverse and optical modes have been ignored, while their effect can be significant [12,65]. Accurate studies of molecules thus require more accurate definitions of their parameters.

Materials and Methods
Analytical estimates use the analysis of resonant interactions. These methods can be qualitatively justified by the similarity of the problem to the exactly solvable localization problem on the Bethe lattice [12,47,48]. We ignore logarithmic factors appearing in these considerations being concentrated on the power law dependencies.
The numerical study exploits the exact diagonalization of Hamiltonian matrices using the standard MatLab software facilities [66].

Conclusions
Here, we briefly summarize the results of the present work. The semi-quantitative theory is developed to determine the critical energy separating localized (integrable) and chaotic behaviors in the quantum FPU chain of atoms with different boundary conditions. The criterion of delocalization has been suggested considering resonant interactions for combined α + β FPU problem. It is predicted that the critical energy decreases with the number of atoms inversely proportionally to this number until the effective thermal energy exceeds the normal mode quantization energy in agreement with previous analysis of the classical β FPU problem (see Equations (21) and (26) and Table 1). At larger numbers of atoms, the critical energy does not depend on this number.
The qualitative behavior predicted by Equation (21) is obtained using resonant language similarly to in Ref. [47], where the matching Bethe lattice problem has been used to justify the results. Similar to Bethe lattice problem, one can expect the appearance of the additional logarithmic factor in Equation (21). However, in our specific case, it is of order of 1 since the argument of logarith is determined by the ratio of diagonal and off-diagonal interactions [47], which have same order of magnitude in the problem under consideration. The quantitative expression in Equation (21) gives a reasonable estimate by order of magnitude but does not pretend to be the accurate expression.
The attempt of numerical verification of the results has been made in the oversimplified model with conserving number of phonons. This model shows that the chaos emerges at smaller energies for free and fixed ends boundary conditions compared to the system with periodic boundary conditions because of the smaller phase space in the latter case. The behaviors obtained are consistent with theory predictions but more realistic models need to be studied for accurate theory verification.
The application of the theory to atomic chains of carbon atoms described by the Morse potential predicts the occurrence of chaotic behavior for very long chains and high system energy that does not agree with experimental observations. Most probably this is because the model describes perfectly linear chains, while realistic (e.g., alkane) chains have more complicated structure and should be modeled with modified parameters. Acknowledgments: Authors acknowledges Sergei Flach, Ivan Khaymovich and Igor Rubtsov for stimulating discussion.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.