Correlation Functions for a Chain of Short Range Oscillators

We consider a system of harmonic oscillators with short range interactions and we study their correlation functions when the initial data is sampled with respect to the Gibbs measure. Such correlation functions display rapid oscillations that travel through the chain. We show that the correlation functions always have two fastest peaks which move in opposite directions and decay at rate $t^{-\frac{1}{3}}$ for position and momentum correlations and as $t^{-\frac{2}{3}}$ for energy correlations. The shape of these peaks is asymptotically described by the Airy function. Furthermore, the correlation functions have some non generic peaks with lower decay rates. In particular, there are peaks which decay at rate $t^{-\frac{1}{4}}$ for position and momentum correlators and with rate $t^{-\frac{1}{2}}$ for energy correlators. The shape of these peaks is described by the Pearcey integral. Crucial for our analysis is an appropriate generalisation of spacings, i.e. differences of the positions of neighbouring particles, that are used as spatial variables in the case of nearest neighbour interactions. Using the theory of circulant matrices we are able to introduce a quantity that retains both localisation and analytic viability. This also allows us to define and analyse some additional quantities used for nearest neighbour chains. Finally, we study numerically the evolution of the correlation functions after adding nonlinear perturbations to our model. Within the time range of our numerical simulations the asymptotic description of the linear case seems to persist for small nonlinear perturbations while stronger nonlinearities change shape and decay rates of the peaks significantly.


Introduction
In this manuscript we consider a system of N " 2M`1 particles interacting with a short range harmonic potential with Hamiltonian of the form H " where 1 ď m ! N , κ 1 ą 0, κ m ą 0, and κ s ě 0 for 1 ă s ă m. In order to make sense of (1.1) we need to introduce boundary conditions. Throughout this paper we consider periodic boundary conditions. By that we mean that the indices j are taken from Z{N Z and therefore q N`j " q j , p N`j " p j holds for all j. The Hamiltonian (1.1) can be rewritten in the form Hpp, qq :" 1 2 xp, py`1 2 xq, Aqy, (1.2) where p " pp 0 , . . . , p N´1 q, q " pq 0 , . . . , q N´1 q, x ., .y denotes the standard scalar product in R N and where A P MatpN, Rq is a positive semidefinite symmetric circulant matrix generated by the vector a " pa 0 , . . . , a N´1 q namely A kj " a pj´kqmod N or where a 0 " 2 m ÿ s"1 κ s , a s " a N´s "´κ s , for s " 1, . . . , m and a s " 0 otherwise. (1.4) Due to the condition κ 1 ą 0 we have xq, Aqy " 0 iff all spacings q j`1´qj vanish. Therefore the kernel of A is one-dimensional with the constant vector p1, . . . , 1q providing a basis. This also implies that the lattice at rest has zero spacings everywhere. Observe, however, that one may introduce an arbitrary spacing ∆ for the lattice at rest by the canonical transformation Q j " q j`j ∆, P j " p j which does not change the dynamics. The periodicity condition for the positions Q j then reads Q N`j " Q j`L with L " N ∆ (see e.g. [17,Sec. 2]). The harmonic oscillator with only nearest neighbour interactions is recovered by choosing a 0 " 2κ 1 , a 1 " a N´1 "´κ 1 , and the remaining coefficients are set to zero. The equations of motion for the Hamiltonian H take the form d 2 dt 2 q j " m ÿ s"1 κ s pq j`s´2 q j`qj´s q, j P Z{N Z.
The integration is obtained by studying the dynamics in Fourier space (see e.g. [9]). In this paper we study correlations between momentum, position and local versions of energy. Following the standard procedure in the case of nearest neighbour interactions we replace the vector of position q by a new variable r so that the Hamiltonian takes the form H " 1 2 xp, py`1 2 xr, ry.
Such a change of variables may be achieved by any linear transformation r " T q, (1.5) with an NˆN matrix T that satisfies A " T T, (1.6) where T denotes the transpose of T . In the case of nearest neighbour interactions one may choose r j " ? κ 1 pq j`1´qj q corresponding to a circulant matrix T generated by the vector τ " ? κ 1 p´1, 1, 0, . . . , 0q. We that satisfies (1.6). The crucial point here is that T is not the standard (symmetric) square root of the positive semidefinite matrix A but a localized version generated by some vector τ with zero entries everywhere, except possibly in the first m`1 components. Hence the j-th component of the generalized elongation r defined through (1.5) depends only on the components q s with s " j, j`1, . . . , j`m. It is worth noting that 1 " p1, . . . , 1q satisfies T 1 " 0 since x1, A1y " 0. This implies m ÿ s"0 τ s " 0 , r j " m ÿ s"1 τ s pq j`s´qj q and N´1 ÿ j"0 r j " p1, . . . , 1qT q " 0.
The local energy e j takes the form e j " 1 2 p 2 j`1 2 r 2 j . The goal of this manuscript is to study the behaviour of the correlation functions for the momentum p j , the generalized elongation r j and the local energy e j when N Ñ 8 and t Ñ 8. Due to the spatial translation invariance of the Hamiltonian Hpp, qq " Hpp, q`λ1q, λ P R, that corresponds to the conservation of total momentum, we reduce the Hamiltonian system by one degree of freedom to obtain a normalizable Gibbs measure. This leads to the reduced phase space We endow M with the Gibbs measure at temperature β´1, namely: dµ " Z N pβq´1δ˜N´1 ÿ k"0 p k¸δ˜N´1 ÿ k"0 q k¸e´β Hpp,qq dpdq (1.9) where Z N pβq is the norming constant and δpxq is the delta function. For convenience we introduce the vector upj, tq " pr j ptq, p j ptq, e j ptqq.
For the harmonic oscillator with nearest neighbor interactions such limits have been calculated in [10].
In an interesting series of papers, (see e.g. [18], and also the collection [8]) several researchers have considered the evolution of space-time correlation functions, for "anharmonic chains", which are nonlinear nearest-neighbor Hamiltonian systems of oscillators. The authors consider the deterministic evolution from random initial data sampled from a Gibbs ensemble, with a large number of particles and study the correlation functions S N αα 1 . In addition to intensive computational simulations [6], [13], Spohn and collaborators also propose and study a nonlinear stochastic conservation law model [17], [18]. Using deep physical intuition, it has been proposed that the long-time behaviour of space-time correlation functions of the deterministic Hamiltonian evolution from random initial data is equivalent to the behaviour of correlation functions of an analogous nonlinear stochastic system of PDEs. Studying this stochastic model, Spohn eventually arrives at an asymptotic description of the "sound peaks" of the correlation functions in normal modes coordinates which are related to S αα 1 by orthogonal transformation: using the notation of [Formula (3.1)] [17]. Here f KPZ is a universal function that first emerges in the Kardar-Parisi-Zhang equation and it is related to the Tracy-Widom distribution, [20], (for a review see [1] and also [4]). A common element to the above cited papers is the observation that such formulae should hold for non-integrable dynamics, while the correlation functions of integrable lattices of oscillators will exhibit ballistic scaling, which means the correlation functions decay as 1 t for t large. For example, in [6] the authors present the results of simulations of the Toda lattice in 3 different asymptotic regimes (the harmonic oscillator limit, the hard-particle limit, and the full nonlinear system). They present plots of the quantity tSpx, tq as a function of the scaled spatial variable x{t (here Spx, tq represents any of the correlation functions). The numerical results support the ballistic scaling conjecture in some of the asymptotic scaling regimes. Further analysis in [19] gives a derivation of the ballistic scaling for the Toda lattice. The decay of equilibrium correlation functions show similar features as anomalous heat transport in one-dimensional systems [3], [7] [2] which leads to conjecture that the two phenomena are related [8].
In [12] the authors also pursue a different connection to random matrices, and in particular to the Tracy-Widom distribution. Over the last 15 years, there has emerged a story originating in the proof that for the totally asymmetric exclusion process on a 1-D lattice (TASEP), the fluctuations of the height function are governed (in a suitable limit) by the Tracy-Widom distribution. Separately, a partial differential equations model for these fluctuations emerged, which takes the form of a stochastic Burgers equation: where ζ is a stationary spatio-temporal white noise process. (The mean behaviour of TASEP is actually   described by the simpler Euler equation  Bu  Bt  "´λu  Bu  Bx  .). From these origins there have now emerged proofs, for a small collection of initial conditions, that the fluctuations of the solution to (1.12) are indeed connected to the Tracy-Widom distribution (see [1] and the references contained therein). In [12], the authors considered continuum limits of anharmonic lattices with random initial data, in which there are underlying conservation laws describing the mean behaviour that are the analogue of the Euler equation associated to (1.12). By analogy with the connection between TASEP and (1.12), they proposed that the time-integrated currents are the analogue of the height function, and should exhibit fluctuations about their mean described by the Tracy-Widom distribution, again based on the use of the nonlinear stochastic pde system as a model for the deterministic evolution from random initial data. As one example, they consider the quantity where upx, tq arises as a sort of continuum limit of a particle system obeying a discrete analogue of a system of conservation laws taking the form B t upx, tq`B x jpx, tq " 0, in which jpx, tq is a local current density for upx, tq. The authors suggest a dual interpretation of Φpx, tq as the height function from a KPZ equation, and thus arrive at the proposal that where a 0 and Γ are model-dependent parameters, and ξ T W is a random amplitude with Tracy-Widom distribution.
Our main result is the analogue of the relations (1.11) for the harmonic oscillator with short range interactions and (1.14) for the harmonic oscillator. For stating our result, we first calculate the dispersion relation |ωpkq| for the harmonic oscillator with short range interaction in the limit N Ñ 8 obtaining see (2.21). The points k " 0, 1 contribute to the fastest moving peaks of the correlation functions that have a velocity˘v 0 where v 0 " a ř m s"1 s 2 κ s " f 1 p0q{p2πq. If f 2 pkq ă 0 for all 0 ă k ď 1{2 then as t Ñ 8 the following holds uniformly in j P Z (cf. Theorem 2.6 and Figure 1): where Aipwq " 1 π ş 8 0 cospy 3 {3`wyqdy, w P R, is the Airy function, and λ 0 : . The above formula is the linear analogue of the Tracy-Widom distribution in (1.11).
Furthermore we can tune the spring intensities κ s , s " 1, . . . , m in (1.15) so that we can find an pm´1qparameter family of potentials such that for j "˘v˚t, with 0 ď v˚ă v 0 , one has In this case the local behaviour of the correlation functions is described by the Pearcey integral (see Theorem 2.7 and Figures 2, 3 below). For example a potential with such behaviour is given by a spring interaction of the form κ s " 1 s 2 for s " 1, . . . , m and m even (see Example 2.8 below). In Section 3.3 we study numerically small nonlinear perturbations of the harmonic oscillator with short range interactions and our results suggest that the behaviour of the fastest peak has a transition from the Airy asymptotic (1.16) to the Tracy-Widom asymptotic (1.11), depending on the strength of the nonlinearity. Namely the asymptotic behaviour in (1.11) that has been conjectured for nearest neighbour interactions seems to persist also for sufficiently strong nonlinear perturbations of the harmonic oscillator with short range interactions. Remarkably, our numerical simulations indicate that the non generic decay in time of other peaks in the correlation functions persists under small nonlinear perturbations with the same power law t´1 {4 as in the linear case, see e.g. Figures 4 and 6.
So as not to overlook a large body of related work, we observe that the quantities we consider here are somewhat different than those considered in the study of thermal transport, though there is of course overlap. (We refer to the Lecture Notes [7] for an overview of this research area and also the seminal paper [15].) As mentioned above, we study the dynamical evolution of space-time correlation functions and the statistical description of random height functions, where the only randomness comes from the initial data. By comparison, in the consideration of heat conduction and transport in low dimensions, anharmonic chains are often connected at their ends to heat reservoirs of different temperatures, and randomness is present primarily in the dynamical laws, not only in fluctuations of initial data.
This manuscript is organized as follows. In Section 2 we study the harmonic oscillator with short range interactions and we introduce the necessary notation and the change of coordinates q Ñ r that enables us to study correlation functions. We then study the time decay of the correlation functions via steepest descent analysis and we show that the two fastest peaks travelling in opposite directions originate from the points k " 0 and k " 1 in the spectrum. Such peaks have a decay described by the Airy scaling. We then show the existence of potentials such that the correlation functions have a slower time decaying with respect to "Airy peaks". In Section 3 we show that the harmonic oscillator with short range interactions has a complete set of local integrals of motion in involution and the correlation functions of such integrals have the same structure as the energy-energy correlation function. Finally, we show that the evolution equations for the generalized position, momentum can be written in the form of conservation laws which have a potential function. For the case of the harmonic oscillator with nearest neighbour interaction, we show that this function is a Gaussian random variable and determine the leading order behaviour of its variance as t Ñ 8. This may be viewed as the analogue of formula (1.14) for the linear case. Technicalities and a description of our numerics are deferred to the Appendix.

The harmonic oscillator with short range interactions
As it was explained in the introduction we rewrite the Hamiltonian for the harmonic oscillator with short range interactions . The matrix T will also play a role in constructing a complete set of integrals that have a local density in the sense that we just described for the energy. In order to state our result we have to introduce some notation. First of all, a matrix A of the form (1.3) with a P R N is called a circulant matrix generated by the vector a.
x k "0, otherwise, while the vectorx P R N is called half-m-physical generated by y P R m`1 if y 0 "´ř m s"1 y s and x k "y k , for 0 ď k ď m x k "0, for m ă k ď N´1.
Following the proof of a classic lemma by Fejér and Riesz, see e.g. [16, pg. 117 f], one can show that a circulant symmetric matrix A of the form (1.2) generated by a m-physical vector a always has a circulant localized square root T that is generated by a half-m physical vector τ .
Proposition 2.2. Fix m P N. Let the circulant matrix A be generated by an m-physical vector a, then there exist a circulant matrix T generated by an half-m-physical vector τ such that: A " T T . (2.1) Moreover, we can choose τ such that ř m s"1 sτ s ą 0. Then one has For example, if we consider m " 1, and a 0 " 2κ 1 and a 1 " a N´1 "´κ 1 . The matrix T is generated by the vector τ " pτ 0 , τ 1 q with τ 0 "´?κ 1 and τ 1 " ? κ 1 . When m " 2 and a 0 " 2κ 1`2 κ 2 , a 1 " a N´1 "´κ 1 , so that the quantities r j are defined as Next we integrate the equation of motions. The Hamiltonian Hpp, qq represents clearly an integrable system that can be integrated passing through Fourier transform. Let F be the discrete Fourier transform with entries F j,k :" 1 Thanks to the above properties, the transformation defined by is canonical. Furthermore s p p j " p p N´j and s p q j " p q N´j , for j " 1, . . . , N´1, while p p 0 and p q 0 are real variables. The matrices T and A are circulant matrices and so they are reduced to diagonal form by F: Let ω j denote the eigenvalues of the matrix T ordered so that FT F´1 " diagpω j q. Then |ω j | 2 are the (non negative) eigenvalues of the matrix A and whereã is the m-physical vector generated by a andτ is the half m-physical vector generated by τ according to Definition 2.1. It follows that The Hamiltonian H, can be written as the sum of N´1 oscillators There are no terms involving p p 0 , p q 0 since the conditions defining M (1.8) imply that p p 0 " 0 and p q 0 " 0. The Hamilton equations are (2.7) Thus the general solution reads: and p q 0 ptq " 0 and p p 0 ptq " 0. Inverting the Fourier transform, we recover the variables q " F´1p q, p " F p p and r " F´1p r where Correlation Decay. We now study the decay of correlation functions for Hamiltonian systems of the form (1.2). We recall the definition (1.9) of the Gibbs measure at temperature β´1 on the reduced phase space M, namely: where Z N pβq is the norming constant of the probability measure. For a function f " f pp, qq we define its average as We first compute all correlation functions (1.10), then we will evaluate the limit N Ñ 8. We first observe that (1.9) in the variables pp p, p qq :" p s Fp, Fqq becomes where dp p j dp q j " d p p j d p p j d p q j d p q j and we recall that p p j " p p N´j , p q j " p q N´j , p r j " ω j p q j , for j " 1, . . . , N´1. From the evolution of p p j and p q j in (2.8) and (2.9), we arrive at the relations (2.14) Now we are ready to compute explicitly the correlation functions in the physical variables. We show the computation for the case S N 11 pj, tq, and we leave to the reader the details for the other cases: S N 11 pj, tq " xr j ptqr 0 p0qy " cos p|ω l |tq cosˆ2π lj N˙" S N 22 pj, tq .
(2. 15) In the same way we have that: The dispersion relation given by (2.4) takes the form where we substitute for the a s their values as in (1.4). We are interested in obtaining the continuum limit of the above correlation functions. We first define ωpkq to provide continuum limits of ω and |ω | 2 , namely where the variable {N has been approximated with k P r0, 1s. One may use equation (A.1) to check the consistency of the two equations of (2.21). To this end observe that ωpkq " Qpe´2 πik q, Ę ωpkq " Qpe 2πik q, and |ωpkq| 2 " pe 2πik q. Lemma 2.3. Let ωpkq be defined as in (2.21), set f pkq :" |ωpkq|, and denote θpkq :" argpωpkqq for 0 ď k ď 1, where the ambiguity in the definition of θ is settled by requiring θ to be continuous with θp0q P p´π, πs. Then, for all k P r0, 1s we have Furthemore, the functions f and θ´π 2 are C 8 on r0, 1s and they both possess odd C 8 -extensions at k " 0 which implies in particular θp0q " π 2 . Proof. The symmetries follow directly from the definition of ω in (2.21). From (2.21) we also learn that |ωpkq| 2 ě 2κ 1 p1´cosp2πkqq ą 0 for k P p0, 1q. Thus the smoothness of f and θ only needs to be investigated for k P t0, 1u. By symmetry we only need to study the case k " 0. The smoothness of the function θ may be obtained from the expansion near k " 0 together with ř m s"1 sτ s ą 0 (see Proposition 2.2). Since cotpθp0qq " 0 and ωpkq ą 0 for small positive values of k we conclude that θp0q " π 2 from the requirement θp0q P p´π, πs. This also implies the existence of a smooth odd extension of θ´π 2 at k " 0 because cotpθpkqq has such an extension. For the function f the claims follow from the representation near k " 0 where sincpxq " sinpxq x denotes the smooth and even sinus cardinalis function. Lemma 2.4. In the limit N Ñ 8 the correlation functions have the following expansion and θpkq " arg ωpkq with ωpkq as in (2.21).
Proof. For any periodic C 8 -function g on the real line with period 1, gpkq " ř nPZĝ n e 2πikn , one has It follows from Lemma 2.3 that the integrands in (2.25)-(2.27) can be extended to 1-periodic smooth functions because we have for small positive values of k that cos pf p´kqtq cos p´2πkjq " cos pf pkqtq cos p´2πkjq " cos pf p1´kqtq cos p2πp1´kqjq , sin pf p´kqtq cos p´2πkj˘θp´kqq "´sin pf pkqtq cos p´2πkj˘pπ´θpkqqq " sin pf p1´kqtq cos p2πp1´kqj˘θp1´kqq .
Observing in addition that the summands corresponding to " 0 are missing in (2.15)-(2.17) the first claim is proved. Together with (2.19) this also implies the second claim.
Next we analyse the leading order behaviour (as t Ñ 8) of the limiting correlation functions S αα 1 pj, tq using the method of steepest descent. In order to explain the phenomena that may occur we start by discussing S 11 . Denote ξ :" j t and φ˘pk, ξq :" f pkq˘2πξk . (2.29) With these definitions and using the symmetry (2.23) we may write The leading order behaviour (t Ñ 8) of such an integral is determined by the stationary phase points k 0 P r0, 1s, i.e. by the solutions of the equation B Bk φ´pk 0 , ξq " 0 which depend on the value of ξ. Such stationary phase points do not need to exist. In fact, as we see in Lemma 2.5 b) below, the range of f 1 is given by some interval r´2πv 0 , 2πv 0 s so that there are no stationary phase points for |ξ| ą v 0 . As in the proof of Lemma 2.4 one can argue that the integrand e itφ´pk,j{tq can be extended to a periodic smooth function of k on the real line with period 1. It then follows from integration by parts that S 11 pj, tq decays rapidly in time. More precisely, for every fixed δ ą 0 we have S 11 pj, tq " O`t´8˘as t Ñ 8, uniformly for |j| ě pv 0`δ qt. (2.31) This justifies the name of sound speed for the quantity v 0 .
In the case |ξ| ď v 0 there always exists at least one stationary phase point k 0 " k 0 pξq P r0, 1s. Each stationary phase point may provide an additive contribution to the leading order behaviour of ş 1 0 e itφ´pk, j{tq dk for j near ξt. However, the order of the contribution depends on the multiplicity of the stationary phase point. For example, let k 0 be a stationary phase point of φ´p¨, ξq, i.e. B Bk φ´pk 0 , ξq " 0. Denote by the smallest integer bigger than 1 for which B Bk φ´pk 0 , ξq ‰ 0. Then k 0 contributes a term of order t 1{ to the t-asymptotics of ş 1 0 e itφ´pk, j{tq dk for j in a suitable neighbourhood of ξt. Before treating the general situation let us recall the case of nearest neighbour interactions. There we have f pkq " f 1 pkq " a 2κ 1 p1´cosp2πkqq " 2 ? κ 1 sinpπkq , k P r0, 1s . The range of f 1 1 equals r´2πv 0 , 2πv 0 s with v 0 " ? κ 1 . For every |ξ| ď v 0 there exists exactly one stationary phase point k 0 pξq P r0, 1s of φ´p¨, ξq that is determined by the relation cospπk 0 pξqq " ξ{v 0 . A straight forward calculation gives Moreover, we have k 0 pv 0 q " 0 and k 0 p´v 0 q " 1 and therefore B 3 Bk 3 φ´pk 0 p˘v 0 q,˘v 0 q "¯2π 3 v 0 ‰ 0. This implies that in addition to (2.31) we have S 11 pj, tq " Opt´1 {2 q, except for j near˘v 0 t where S 11 pj, tq " Opt´1 {3 q. In order to determine the behaviour near the least decaying peaks that travel at speeds˘v 0 we expand f 1 near the stationary phase points. Let us first consider ξ " v 0 with k 0 " 0. Introducing Substituting y " 2πλ 0 t 1{3 k leads for k close to 0 to the asymptotic expression Using the well-known representation Aipwq " 1 π ş 8 0 cospy 3 {3`wyqdy, w P R, of the Airy function and performing a similar analysis around the stationary phase point k 0 "´1 for ξ "´v 0 one obtains an asymptotic formula for the region not covered by (2.31) {2¯, t Ñ 8, uniformly for |j| ă pv 0`δ qt (2.32) for δ ą 0 (see e.g. [14]). Observe that due to the decay of Aipwq for w Ñ˘8, the Airy term is dominant roughly in the regions described by v 0 t´optq ă |j| ă v 0 t`oppln tq 2{3 q.
From the arguments just presented it is not difficult to see that the derivation of (2.32) only uses the following properties of f " f 1 :  2.5 b). Then the decay rate of S 11 pj, tq for j near vt is at most of order t´1 {3 . The decay is even slower (at least of order t´1 {4 ) if f 3 pk 0 q " 0 holds in addition. We show in Theorem 2.7 that this may happen for κ in some submanifold of R m of codimension 1 (see also Examples 2.8 and 2.9). Nevertheless, if κ 2 , . . ., κ m are sufficiently small in comparison to κ 1 then condition (2.33) is always satisfied as we show in Theorem 2.6 c).
Proof. Statement a) follows directly from the last formula in the proof of Lemma 2.3 and from the expansion sinc 2 pxq " 1´x This representation also settles statement c). As we know already f 1 p0q " 2πv 0 "´f 1 p1q we may establish statement b) by verifying that |f 1 pkq| ă 2πv 0 holds for all k P p0, 1q. To this end we write f " p ř m s"1 h 2 s q 1{2 with h s pkq " 2 ? κ s sinpπskq. Using the Cauchy-Schwarz inequality we obtain for 0 ă k ă 1 that where the last inequality follows from | cospπkq| ă 1 and κ 1 ą 0.
Proof. The rapid decay claimed in statement a) can be argued in the same way as (2.31) for S 11 " S 22 . Due to relations (2.18) and (2.28) one only needs to consider S 12 and S 21 . Indeed, using Lemma 2.3 one may show that the imaginary parts of the integrands used in the representation of S 12 and S 21 in (2.39) below have smooth extensions to all k P R that are 1-periodic. This is all that is needed because | B Bk φ˘pk, j{tq| ą 2πδ by Lemma 2.5 b) uniformly for k P r0, 1{2s and |j| ą pv 0`δ qt.
Statement c) follows from the continuous dependence of the derivatives f 2 and f 3 on the parameters pκ 2 , . . . , κ m q (see Lemma 2.5 c) and from simple facts for the case of nearest neighbour interactions f 1 pkq " 2 ? κ 1 sinpπkq discussed above. Indeed, from f 2 p0q " 0 and from f 3 1 p0q ă 0 it follows that there exists such an ε ą 0 such that f 3 pkq ă 0 and hence also f 2 pkq ă 0 for k in some region p0, δq uniformly in pκ 2 , . . . , κ m q P r0, εq m´1 . As f 2 1 pkq ă´2π 2 ? κ 1 sinpπδq for all k P rδ, 1{2s we may prove the claim in this region by reducing the value of ε if necessary. Theorem 2.6 provides the leading order asymptotics of the limiting correlations S αα 1 pj, tq for t Ñ 8 in the simple situation that the second derivative of the dispersion relation is strictly negative on the open interval p0, 1q (cf. condition (2.33)). Moreover, statement c) shows that there is a set of positive measure in parameter space κ P R m where this happens. For general values of κ, however, different phenomena may appear. In particular, there might exist stationary phase points of higher order leading to slower time-decay of the correlations (see discussion before the statement of Lemma 2.5). By a naive count of variables and equations one might expect that decay rates t´1 {p3`pq occur on submanifolds of parameter space of dimension m´p. Theorem 2.7 shows that this is indeed the case for p " 1. Moreover, we present in this situation a formula for the leading order contribution of the corresponding stationary phase points to the asymptotics of S αα 1 pj, tq. Despite being non-generic in parameter space it is interesting to note that decay rates t´1 {4 can be observed numerically (see Figures 2 and 3). There is also a second issue that may arise if condition (2.33) fails. Namely, for v P p´v 0 , v 0 q there can be several values of k P p0, 1 2 s satisfying f 1 pkq˘2πv " 0 so that the contributions from all these stationary points need to be added to describe the leading order behaviour for j near vt. a) For m ě 3 there is an pm´1q-parameter family of potentials for which there exist k˚" k˚pκq P p0, 1 2 q with f 2 pk˚q " 0, f 3 pk˚q " 0, f pivq pk˚q ‰ 0, and 0 ă v˚:" f 1 pk˚q 2π ă v 0 , (2.40) Figure 1. Correlation functions S αα 1 for the harmonic oscillator with nearest neighbour interaction with κ 1 " 1 (top left) and the harmonic potential with κ s " 1 s 2 for s " 1, 2 in Example 2.8 (center left) and the potential of Example 2.9 in the bottom left. In the second column the Airy scaling (2.36) of the corresponding fastest moving peaks. The Airy asymptotic is perfectly matching the fastest peak and capturing several oscillations.
with v 0 as in (2.35). Set λ˚:" 1 2π p|f pivq pk˚q|{4!q 1 4 ą 0. Then for j Ñ 8 and t Ñ 8 in such a way that v˚t´j λ˚t 1 4 is bounded, the contribution of the stationary phase point k˚to the correlation functions is given by: and P˘has to be chosen according to the sign of f pivq pk˚q. If j Ñ´8 with bounded pv˚t`jq{pλ˚t 1{4 q the contributions of the stationary point k˚can be obtained from the ones presented in (2.41)-(2.43) by replacing φ´by φ`, θ by´θ, and j in the argument of P˘by´j. b) When k˚" 1 2 one has f 1 p1{2q " 0 and f 3 p1{2q " 0 by the symmetry (2.23). For each m ě 2 there is an pm´1q-parameter family of potentials so that f 2 p1{2q " 0 and f pivq p1{2q ‰ 0 holds in addition. In this case the contribution of the stationary phase point k˚" 1{2 to the correlation functions in the asymptotic regime t Ñ 8 with bounded j{t  itpf pkq`2πk j t q`eitpf pkq´2πk j t q¯d k. (2.46) In order to compute the contribution of the stationary phase point k˚to the large t asymptotics of the integral in (2.46) we expand Introducing the change of variables In the situation k˚" 1{2 of statement b) one uses in addition that tφ˘p1{2, j{tq " tf p1{2q˘jπ, ωp1{2q " ř m s"1 τ s p1´cospπsqq "´2 . In this way and with the help of (2.28) all relations of (2.45) can be deduced. We now show the existence of a codimension 1 manifold in parameter space that exhibits such higher order stationary phase points in the situation of b) where k˚" 1{2. As we have f 3 p1{2q " 0 by symmetry (2.23) we only need to solve The solution of the above equation is It is clear from the above relation that for m even, choosing κ 1 sufficiently big one has κ m ą 0 while for m odd, it is sufficient to choose κ s`1 ą s 2 ps`1q 2 κ s ą 0, s odd and 1 ď s ď m´2. Note that in the situation of (2.48) f pivq p 1 2 q ‰ 0 holds iff ř m s"1 κ s s 4 p´1q s`1 ‰ 0. This condition simply removes an pm´2q-dimensional plane from our manifold (2.48) which defines a hyperplane in the positive cone of the m-dimensional parameter space. Therefore we have found an pm´1q-parameter family of potentials such that the correlation functions decay as in (2.45).
We observe that the first term in the above expression is always positive, while the second term is positive if we require that κ s ą ps`2q 2 s 2 κ s`2 ą 0 for s " 6, 10, 14, . . . , m´2. The remaining two conditions f 1 p1{4q ą 0 and f pivq p1{4q ‰ 0 are easy to satisfy: The sign of f 1 p1{4q agrees with the sign of ř s odd sκ s p´1q s´1 2 and can be made positive by choosing κ 1 sufficiently large. In the situation where (2.49) and (2.50) hold the fourth derivative f pivq p1{4q does not vanish iff ř s even s 4 κ s p´1q s 2 ‰ 0. This can be achieved by adjusting, for example, the value of κ 2 . We have now shown that there exists κ˚P R m such that the first four derivatives of f have all desired properties at k " 1{4. In order to obtain the pm´1q-dimensional solution manifold in parameter space, we apply the Implicit Function Theorem to F pk, κq :" pf 2 pk, κq, f 3 pk, κqq. By a straight forward computation on sees that We can therefore solve F pk, κq " 0 near p1{4, κ˚q by choosing pk, κ 4 q as functions of the remaining parameters κ j with j ‰ 4. Example 2.8. m even. Choosing κ s " 1 s 2 for s " 1, . . . , m one has that conditions (2.47) are satisfied and f pivq`1 2˘ă 0. For κ s " 1 s α , s " 1, . . . , m´1, 2 ă α ă 3, and κ m given by (2.48), there is α " αpmq such that κ m ă κ m´1 . m odd. Choosing κ s " 1 s , for s " 1, . . . m´1, one has from (2.48) κ m " m´1 2m 2 ă κ m´1 and f pivq p 1 2 q ą 0. In all these examples the correlation functions S αα 1 pj, tq, α, α 1 " 1, 2 decrease as t´1 4 near j " 0. Example 2.9. We consider the case m " 3 and we want to get a potential that satisfies (2.40) with v˚ą 0. We chose as a critical point of f pkq the point k˚" 1 3 thus obtaining the equations The speed of the peak is v˚" ? 2κ 1 4 and f pivq p 1 3 q "´6 8 ? 6 6 π 4 ? κ 1 . The correlation functions S αα 1 pj, tq, α, α 1 " 1, 2 decrease as t´1 4 and S 33 pj, tq decreases like t´1 2 as t Ñ 8 and j " v˚t, see Figure 3. Note that one may obtain a 2-parameter family of solutions of (2.40) by picking, for example, the particular solution related to κ 1 " 1 and by showing that the system of equations pf 2 , f 3 qpk, κq " 0 can be solved near (1/3, 1, 1/8, 7/72) by choosing k and κ 3 as functions of κ 1 and κ 2 using the Implicit Function Theorem in the same way as at the end of the proof of Theorem 2.7. 3. Complete set of integrals with local densities, currents and potentials, and some numerics for nonlinear versions 3.1. Circulant hierarchy of integrals. In this section we construct a complete set of conserved quantities that have local densities. The harmonic oscillator with short range interaction is clearly an integrable system. A set of integrals of motion is given by the harmonic oscillators in each of the Fourier variables: However, when written in the physical variables p and q, the quantities depend on all components of the physical variables. We now construct integrals of motion each having a density that involves only a limited number of components of the physical variables and this number only depends on the range m of interaction.
For this purpose we denote by te k u N´1 k"0 the canonical basis in R N . Theorem 3.1. Let us consider the Hamiltonian

1)
with the symmetric circulant matrix A as in (1.2), (1.3). Define the matrices tG k u M k"1 to be the symmetric circulant matrix generated by the vector 1 2 pe k`eN´k q and tS k u M k"1 to be the antisymmetric circulant matrix generated by the vector 1 2 pe k´eN´k q. Then the family of Hamiltonians defined as together with H 0 :" H forms a complete family pH j q 0ďjďN´1 of integrals of motion that, moreover, is in involution.
Proof. Observe first that the Hamiltonian H 0 " H is included in the description of formula (3.2) as G 0 equals the identity matrix. Using the symmetries G k " G k , 0 ď k ď pN´1q{2, the Poisson bracket tF, Gu " x∇ q F, ∇ p Gy´x∇ q G, ∇ p F y may be evaluated in the form All these expressions vanish. To see this, it suffices to observe that multiplication is commutative for circulant matrices and, for the bottom line, that S is skew symmetric: S "´S . We present explicit formulas for the limits S k,n in Appendix C from which one can deduce that they have the same scaling behaviour as the energy-energy correlation function S 33 when t Ñ 8.

Currents and potentials.
In this subsection we write the evolution with respect to time of r j , p j and e j in the form of a (discrete) conservation law by introducing the currents. Each conservation law has a potential function that is a Gaussian random variable. In the final part of this subsection we determine the leading order behaviour of the variance of this Gaussian random variable as t Ñ 8 in the case of nearest neighbour interactions. For introducing the currents we recall that r " T q with T as in (1.7). Then one has τ pr j´rj´ q, p j`N " p j , j " 0, . . . , N´1. We define a potential function for the above conservation law Φpj, tq :" Then it is straightforward to verify that Φ t pj, tq " J pj, tq and Φpj, tq´Φpj´1, tq " upj, tq. The quantities Φ 1 pj, tq and Φ 2 pj, tq can be expressed as sums of independent centered Gaussian random variables and are therefore also Gaussian random variables with zero mean and variance xpΦ 1 pj, tqq 2 y and xpΦ 2 pj, tqq 2 y, where all the averages are taken with respect to the distribution (1.9), see also (2.10). We calculate the variance for the case of the harmonic oscillator with nearest neighbour interactions. In this particular case (3.11) After some lengthy calculations one obtains: Evaluating the r.h.s. of the above expressions in the limit t Ñ 8 we arrive to the following theorem.
Theorem 3.2. In the limit N Ñ 8 and t Ñ 8 the quantities Φ 1 pj, tq and Φ 2 pj, tq defined in (3.11) are Gaussian random variables that have the following large t behaviour: (3.14) The leading order behaviour of the variances σ 2 1 and σ 2 2 agrees. In the physically interesting region |j| t ď ? κ 1 it is given by The proof of the above theorem relies on steepest descent analysis of the oscillatory integrals in (3.13). But because the integrand is actually quite large ( " Ct 2 ) near k " 0, we consider the following Cauchy-type integral instead, which gives the leading order asymptotic behaviour of the integrals appearing in (3.13), since |ωpkq|´2p1´cos p|ωpkq|tqq cos p2πpj`1qkqdk´F 0 p0q Ñ 0 as t, j Ñ 8 . For |j| t ă p1´ q ? κ 1 , ą 0, the analysis of F 0 pzq is quite straightforward -a standard stationary phase calculation combined with a contour deformation to permit the evaluation at z " 0. For t and j growing to 8 such that |j| t « ? κ 1 , the analysis is more complicated because the point of stationary phase is encroaching upon the origin, where the integrand itself is actually large as t Ñ 8. For this case, one must construct a local parametrix, following quite closely the analysis presented in [5], and we omit the details of this analysis. In order to analyse Φ 1 observe that the difference of the integrals in relations (3.12) and (3.13) is given by ş 1 0 |ωpkq|´2 r1´cos p2πpj`1qkqs dk which can also be treated by a stationary phase calculation combined with a contour deformation.
We numerically compute and study the correlatios functions for these systems sampling the initial conditions according to the Gibbs measures of just their harmonic part at temperature β´1 " 1.
In the weakly nonlinear case, the fastest peaks of the correlation functions scale numerically according to the Airy parametrices (cf. Theorem 2.6) as can be deduced from the top pictures in Figures 4, 5 while for stronger nonlinearity the fastest peaks seem to scale like t  Figure 4, still scales in time like t´1 4 . Indeed performing a regression analysis of the log-log plot one can see a scaling like t´0 .267 that is slightly faster then t´1 4 (see Figure 6).  The existence of the τ j 's is a consequence of the Fejér-Riesz lemma. For the convenience of the reader we present a proof following the presentation in [16, pg. 117 f]. Denote by P the polynomial of degree 2m given by P pzq :" z m pzq. Observe that for all x P R we have By the positivity of κ 1 equality holds in the inequality above iff cospxq " 1. This implies that P has no zeros on the unit circle |z| " 1 except for z " 1. We denote by η k , 1 ď k ď r ă , the zeros of P that lie within the unit disc |η k | ă 1 and by ξ k , 1 ď k ď r ą , the zeros of P with |ξ k | ą 1, recorded repeatedly according to their multiplicities, so that P pzq "´κ m pz´1q r0 Using the uniqueness of such a factorization for any polynomial together with the relation z 2m P pz´1q " P pzq one obtains that r ă " r ą and that the zeros can be listed in such a way that η k " ξ´1 k for all 1 ď k ď r ă . Moreover, we learn that r 0 is even with 1 ď :" r 0 {2 " m´r ă . Now it follows from formula (A.2) that lpzq " z´mP pzq " c pz´1´1q pz´1q  The speed ξ 0 of the fastest peak is determined numerically. The decay rate of the slower moving peaks that are scaling like t´1 {4 in the linear case (see Figure 1), is not very clear due to their highly oscillatory behaviour.
Choosing d P C with d 2 " c we see that Qpzq :" dpz´1q ś ră k"1 pz´ξ k q satisfies (A.1). Next we show that the coefficients of the polynomial Q are real. To this end observe that P has real coefficients and therefore all non-real zeros of P come in complex conjugate pairs with equal multiplicities. Therefore the polynomial d´1Qpzq " ř m j"0 s j z j has only real coefficients s j . Relation (A.1) implies a 0 " d 2 ř m j"0 s 2 j . Consequently, d 2 is the quotient of two positive numbers and d must be real. Thus we have τ j " ds j P R for all 0 ď j ď m. We complete the proof by arguing that ř m s"0 τ s " 0 and p ř m s"1 sτ s q 2 " ř m s"1 s 2 κ s hold true. This can be deduced from (A.1) via Qp1q 2 " p1q " 0 and´2Q 1 p1q 2 " 2 p1q "´ř m s"1 2s 2 κ s .

B. Pearcey integral
The general Pearcey integral is defined as P pb, aq :" ż 8 8 e ipt 4`b t 2`a tq dt, 0 ď arg b ď π, a P R. (B.1) This integral decribes cusp singularities in physical phenomena, like the semiclassical limit of the linear Schrödinger equation. The integral (B.1), after a rotation of the integration path through an angle of π{8 that removes the rapidly oscillatory term e it 4 , can be written in the formP pb, aq " 2e iπ{8 P pbe´i π{4 , ae iπ{8 q, with P pb, aq :"  Figure 4 for S 11 pj, tq and S 21 pj, tq and several values of times. The peak is highly oscillatory and the oscillations are interpolated by the red line that suggests a scaling of the correlation function S 11 pj, tq and S 21 pj, tq near j " 0 compatible with t´1 4 .
We are interested in the case b " 0. The corresponding integral is absolutely convergent for all complex values of a and represents the analytic continuation of the Pearcey integral. For the Pearcey integral P`paq defined in (2.44) we obtain P`paq " 2e iπ{8 P p0, ae iπ{8 q .
Note that the integral P´paq, also defined in (2.44), can be related to the function P by rotating the integration path by an angle of´π{8. This gives P´paq " 2e´i π{8 P p0, ae´i π{8 q. From this we learn that P´pāq " P`paq. On the reals we therefore have P´paq " P`paq , a P R .

D. Numerical Computation
The numerical computations have been implemented with Python software, all codes are available on GitHub [11]. Fig. 1-3 are the result of the numerical evaluation via the standard routine numpy.trapz of the integrals in (2.25)-(2.28) for various values of j and t and then we just added the Airy function (1.16) and the Pearcey integral (2.45).
To obtain Fig. 4 we proceed in the following way. First we sampled a random initial data according to the Gibbs measure defined by the corresponding harmonic part of (3.18), namely the Hamiltonian of Example 2.8 with m " 2. We let these data evolve according to the Hamilton equations of (3.18) and compute the values of the correlations function. Then we repeated this procedure 4ˆ10 6 times and we averaged the values of the correlations functions. On the left panel we plot the correlation functions, instead on the right one we focus on the extreme peak and we guess a proper scaling depending on the size of the perturbation. Fig. 5 is made in a similar way, where now the nonlinear potential has the same harmonic part as Example 2.9.
In Fig. 6 we focus our attention on the central peak of the chain with potential as is Fig. 4. We follow the same procedure as before and plot in logarithmic scale the average scaling of the highest peak in the center of the chain. We decide to plot the average height of this peak since it is highly oscillatory and it is difficult to precisely track the oscillations.