1 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

$$\begin{aligned} H=\sum _{j=0}^{N-1}\frac{p_j^2}{2}+\sum _{s=1}^m\frac{\kappa _s}{2}\sum _{j=0}^{N-1}(q_j-q_{j+s})^2\,, \end{aligned}$$
(1.1)

where \(1\le m\ll N\), \(\kappa _1>0\), \(\kappa _m>0\), and \(\kappa _s\ge 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 \({\mathbb Z}/ N{\mathbb Z}\) and therefore

$$\begin{aligned} q_{N+j}=q_j,\quad p_{N+j}=p_j \end{aligned}$$

holds for all j. The Hamiltonian (1.1) can be rewritten in the form

$$\begin{aligned} H(\mathbf{p},\mathbf{q}) := \frac{1}{2} \langle \mathbf{p}, \mathbf{p}\rangle +\frac{1}{2}\langle \mathbf{q},A \mathbf{q}\rangle , \end{aligned}$$
(1.2)

where \(\mathbf{p}=(p_0,\ldots ,p_{N-1}) \), \(\mathbf{q}=(q_0,\ldots ,q_{N-1}),\) \(\langle \,.\,,\,.\rangle \) denotes the standard scalar product in \({\mathbb R}^N\) and where \(A\in \text{ Mat }(N,{\mathbb R})\) is a positive semidefinite symmetric circulant matrix generated by the vector \(\mathbf{a}=(a_0,\ldots ,a_{N-1})\) namely \(A_{kj}=a_{(j-k)\text{ mod } N}\) or

$$\begin{aligned} A = {\begin{bmatrix} a_{0}&{} a_{1}&{}\ldots &{} a_{N-2}&{} a_{N-1}\\ a_{N-1}&{} a_{0}&{} a_{1}&{}&{} a_{N-2}\\ \vdots &{} a_{N-1}&{} a_{0}&{}\ddots &{}\vdots \\ a_{2}&{}&{}\ddots &{}\ddots &{} a_{1} \\ a_{1}&{} a_{2}&{}\ldots &{} a_{N-1}&{} a_{0}\\ \end{bmatrix}}\, , \end{aligned}$$
(1.3)

where

$$\begin{aligned} \begin{aligned}&a_0=2\sum _{s=1}^m\kappa _s,\quad a_s=a_{N-s}=-\kappa _s,\quad \text{ for } s=1,\ldots ,m \text { and }a_s=0 \text { otherwise.} \end{aligned} \end{aligned}$$
(1.4)

Due to the condition \(\kappa _1 >0\) we have \(\langle \mathbf{q},A \mathbf{q}\rangle =0\) iff all spacings \(q_{j+1}-q_j\) vanish. Therefore the kernel of A is one-dimensional with the constant vector \((1,\ldots , 1)^\intercal \) providing a basis. This also implies that the lattice at rest has zero spacings everywhere. Observe, however, that one may introduce an arbitrary spacing \(\Delta \) for the lattice at rest by the canonical transformation \(Q_j = q_j + j\Delta \), \(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 \Delta \) (see e.g.[17, Sect. 2]).

The harmonic oscillator with only nearest neighbour interactions is recovered by choosing

$$\begin{aligned} a_0=2\kappa _1,\quad a_1=a_{N-1}=-\kappa _1, \end{aligned}$$

and the remaining coefficients are set to zero.

The equations of motion for the Hamiltonian H take the form

$$\begin{aligned} \dfrac{d^2}{dt^2}q_j=\sum _{s=1}^m \kappa _s (q_{j+s}-2q_j+q_{j-s}),\quad j\in {\mathbb Z}/ N{\mathbb Z}. \end{aligned}$$

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 \(\mathbf{q}\) by a new variable \(\mathbf{r}\) so that the Hamiltonian takes the form

$$\begin{aligned} H=\frac{1}{2}\langle \mathbf{p},\mathbf{p}\rangle +\frac{1}{2}\langle \mathbf{r},\mathbf{r}\rangle . \end{aligned}$$

Such a change of variables may be achieved by any linear transformation

$$\begin{aligned} \mathbf{r}=T\mathbf{q}, \end{aligned}$$
(1.5)

with an \(N\times N\) matrix T that satisfies

$$\begin{aligned} A = T^\intercal T, \end{aligned}$$
(1.6)

where \( T^\intercal \) denotes the transpose of T. In the case of nearest neighbour interactions one may choose \(r_j = \sqrt{\kappa _1} (q_{j+1}-q_j)\) corresponding to a circulant matrix T generated by the vector \(\varvec{\tau }=\sqrt{\kappa _1}(-1, 1,0,\ldots ,0)\). We show in Proposition 2.2 below that short range interactions given by matrices A of the form (1.3), (1.4) also admit such a localized square root. More precisely, there exists a circulant \(N\times N\) matrix T of the form

$$\begin{aligned} T ={\begin{bmatrix} \tau _0&{} \tau _{1}&{}\ldots &{} \tau _{m}&{}0 &{}\ldots &{}0 \\ 0&{} \tau _{0}&{} \tau _{1}&{}\ldots &{}\tau _m&{}0&{} \\ &{}\ddots &{} \ddots &{} \ddots &{}\ddots &{}&{}\\ \tau _m&{}0&{}\ddots &{} \ddots &{} \ddots &{}\ddots &{}\\ &{}\ddots &{} \ddots &{} \ddots &{}\ddots &{} \ddots &{}\ddots \\ \tau _2&{}\ldots &{}\tau _m&{}0&{}\ldots &{}\tau _0&{}\tau _1\\ \tau _1&{}\tau _2&{}\ldots &{}\tau _m &{}0&{}0&{}\tau _{0}\\ \end{bmatrix}}\, . \end{aligned}$$
(1.7)

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 \(\varvec{\tau }\) with zero entries everywhere, except possibly in the first \(m+1\) components. Hence the j-th component of the generalized elongation \(\mathbf{r}\) defined through (1.5) depends only on the components \(q_s\) with \(s=j,j+1,\ldots ,j+m\). It is worth noting that \({\varvec{1}}=(1,\ldots , 1)^\intercal \) satisfies \(T {\varvec{1}}=0\) since \(\langle {\varvec{1}}, A{\varvec{1}}\rangle =0\). This implies

$$\begin{aligned} \sum _{s=0}^m\tau _s=0\,, \quad r_j= \sum _{s=1}^{m} \tau _s (q_{j+s}-q_j)\, \quad \text {and} \quad \sum _{j=0}^{N-1}r_j= (1,\ldots , 1) T \mathbf{q}= 0. \end{aligned}$$

The local energy \(e_j\) takes the form

$$\begin{aligned} e_j=\dfrac{1}{2}p_j^2+\dfrac{1}{2}r_j^2\,. \end{aligned}$$

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\rightarrow \infty \) and \(t\rightarrow \infty \). Due to the spatial translation invariance of the Hamiltonian \(H(\mathbf{p},\mathbf{q}) = H(\mathbf{p},\mathbf{q}+\lambda {\varvec{1}})\), \(\lambda \in {\mathbb 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

$$\begin{aligned} {\mathcal M}:= \left\{ (\mathbf{p},\mathbf{q})\in {\mathbb R}^{N}\times {\mathbb R}^{N} \, : \, \sum _{k=0}^{N-1} p_k =0\, ;\, \sum _{k=0}^{N-1} q_k = 0 \right\} . \end{aligned}$$
(1.8)

We endow \({\mathcal M}\) with the Gibbs measure at temperature \(\beta ^{-1}\), namely:

$$\begin{aligned} \mathrm{d}\mu = Z_N(\beta )^{-1} \delta \left( \sum _{k=0}^{N-1} p_k\right) \delta \left( \sum _{k=0}^{N-1} q_k \right) e^{-\beta H(\mathbf{p},\mathbf{q})}\mathrm{d}\mathbf{p}\mathrm{d}\mathbf{q}\end{aligned}$$
(1.9)

where \(Z_N(\beta )\) is the norming constant and \(\delta (x)\) is the delta function.

For convenience we introduce the vector

$$\begin{aligned} {\varvec{u}}(j,t)=(r_j(t),p_j(t),e_j(t)). \end{aligned}$$

We consider the correlation functions

$$\begin{aligned} \begin{aligned}&S^N_{\alpha \alpha '}(j,t) = \left\langle u_{\alpha }(j,t)u_{\alpha '}(0,0) \right\rangle - \left\langle u_{\alpha }(j,t)\right\rangle \left\langle u_{\alpha '}(0,0) \right\rangle ,\;\;\alpha ,\alpha '=1,2,3, \end{aligned} \end{aligned}$$
(1.10)

where the symbol \(\left\langle \,.\,\right\rangle \) refers to averages with respect to \(d\mu \) . We calculate the limits

$$\begin{aligned} \lim _{N\rightarrow \infty }S^N_{\alpha \alpha '}(j,t)=S_{\alpha \alpha '}(j,t)\,. \end{aligned}$$

For the harmonic oscillator with nearest-neighbour interactions such limits have been calculated in [10].

In an interesting series of papers, (see e.g. [18], and also the collection [7]) several researchers have considered the evolution of space-time correlation functions, for “anharmonic chains”, which are nonlinear nearest-neighbour 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_{\alpha \alpha '}^{N}\).

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_{\alpha \alpha '}\) by orthogonal transformation:

$$\begin{aligned} \tilde{S}_{\alpha \alpha } \cong \left( \lambda _{s} t \right) ^{-2/3} f_{\text{ KPZ }} \left( ( \lambda _{s} t )^{-2/3} ( x - \alpha c t ) \right) \ , \end{aligned}$$
(1.11)

using the notation of [17, Formula (3.1)]. Here \( f_{\text{ 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 \(\frac{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 tS(xt) as a function of the scaled spatial variable x/t (here S(xt) 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 [2, 3, 7] which leads to conjecture that the two phenomena are related [8].

There is one aspect regarding all the correlations we consider that we wish to point out. For linear force laws they are oscillating rapidly in the region between the sound peaks, except for those regions that display slower decay rates (see Figs. 1 and 3). In agreement with [6] we observe numerically that these oscillations persist in the weakly nonlinear regime but disappear for fully nonlinear systems (see Figs. 4, 5).

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:

$$\begin{aligned} \frac{\partial u}{\partial t} = \nu \frac{\partial ^{2} u}{\partial x^{2} } - \lambda u \frac{\partial u}{\partial x} + \frac{\partial \zeta }{\partial x} \ , \end{aligned}$$
(1.12)

where \(\zeta \) is a stationary spatio-temporal white noise process. (The mean behaviour of TASEP is actually described by the simpler Euler equation \( \displaystyle \frac{\partial u}{\partial t} = - \lambda u \frac{\partial u}{\partial x}\).). 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

$$\begin{aligned} \Phi (x,t)=\int _{0}^{t} \mathfrak {j}(x,t') dt' - \int _{0}^{x} \mathfrak {u}(x',0)dx' \ , \end{aligned}$$
(1.13)

where \(\mathfrak {u}(x,t)\) arises as a sort of continuum limit of a particle system obeying a discrete analogue of a system of conservation laws taking the form \(\partial _{t} \mathfrak {u}(x,t) + \partial _{x} \mathfrak {j}(x,t) = 0\), in which \(\mathfrak {j}(x,t)\) is a local current density for \(\mathfrak {u}(x,t)\). The authors suggest a dual interpretation of \(\Phi (x,t)\) as the height function from a KPZ equation, and thus arrive at the proposal that

$$\begin{aligned} \Phi (x,t) \simeq a_{0} t + \left( \Gamma t \right) ^{1/3} \xi _{TW} \ , \end{aligned}$$
(1.14)

where \(a_{0}\) and \(\Gamma \) are model-dependent parameters, and \(\xi _{TW}\) is a random amplitude with Tracy-Widom distribution.

Fig. 1
figure 1

Correlation functions \(S_{\alpha \alpha '}\) for the harmonic oscillator with nearest-neighbour interaction with \(\kappa _1=1\) (top left) and the harmonic potential with \(\kappa _s=\frac{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

Fig. 2
figure 2

Correlation function \(S_{33}(j,t)\) for the potential \(\kappa _s=1/s^2\) for \(m=2\) in Example 2.8 for several values of time on the left. On the right one sees that the Pearcey scaling provided in (2.45) matches perfectly for the central peak of \(S_{33}(j,t)\)

Fig. 3
figure 3

Potential of Example 2.9. The top left figure displays \(S_{33}(j,t)\) for several values of t. The scaling of \(S_{33}\) according to the Airy function in Theorem 2.6 for the fastest moving peak and the scaling of the slower moving peak according to the Pearcey integral are shown top right and bottom left, respectively. The corresponding critical points of the derivative of the dispersion function can be seen in the bottom right figure

Fig. 4
figure 4

Correlation function \(S^{(N)}_{11}(j,t)\) for several values of times and for the perturbed Hamiltonian (3.18) with \(\kappa _s\) as in Example 2.8, \(\chi =0.01\) and \(\gamma =0.001\) in the top figure and \(\chi =0.1\) and \(\gamma =0.01\) in the lower figure. On the right top figure, the scaling of the fastest peak according to Airy parametrix (see Theorem 2.6 and Fig. 1) and according to \(t^{-2/3}\) in the lower figure. The speed \(\xi _0\) of the fastest peak is determined numerically. One can see that the central peak has a low decay in the top left figure, while in the left bottom figure it is destroyed by the relatively stronger nonlinearity

Fig. 5
figure 5

Correlation function \(S^{(N)}_{22}(j,t)\) for several values of times and for the perturbed Hamiltonian (3.18) with \(\kappa _s\) as in Example 2.9, \(\chi =0.01\) and \(\gamma =0.001\) in the top figure and \(\chi =0.1\) and \(\gamma =0.01\) in the lower figure. The right top figure shows the scaling of the fastest peak compatible with the Airy parametrix and according to \(t^{-2/3}\) in the lower figure. The speed \(\xi _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 Fig. 1), is not very clear due to their highly oscillatory behaviour

Fig. 6
figure 6

Logarithmic plot of the central peak of the example in Fig. 4 for \(S_{11}(j,t)\) and \(S_{21}(j,t)\) 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}(j,t)\) and \(S_{21}(j,t)\) near \(j\sim 0\) compatible with \(t^{-\frac{1}{4}}\)

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 \(|\omega (k)|\) for the harmonic oscillator with short range interaction in the limit \(N\rightarrow \infty \) obtaining

$$\begin{aligned} f(k)=|\omega (k)|=\sqrt{2\sum _{s=1}^m\kappa _s\left( 1-\cos (2\pi k s)\right) }\,, \end{aligned}$$
(1.15)

see (2.21). The points \(k= 0,1\) contribute to the fastest moving peaks of the correlation functions that have a velocity \(\pm v_0\) where \(v_0= \sqrt{\sum _{s=1}^ms^2\kappa _s}=f'(0)/(2\pi )\). If \(f''(k)<0\) for all \(0<k\le 1/2\) then as \(t \rightarrow \infty \) the following holds uniformly in \(j \in {\mathbb Z}\) (cf. Theorem 2.6 and Figs. 1 and 3):

$$\begin{aligned} \begin{aligned}&S_{\alpha \alpha '}(j,t)= \frac{1}{2 \beta \lambda _0 t^{1/3}}\\&\left[ (-1)^{\alpha +\alpha '}\text{ Ai }\left( \dfrac{j-v_0t}{ \lambda _0 t^{1/3}}\right) + \text{ Ai }\left( -\dfrac{j+v_0t}{ \lambda _0 t^{1/3}}\right) \right] +{\mathcal O}\left( t^{-1/2}\right) ,\quad \alpha ,\alpha '=1,2\\ S_{33}(j,t)&= \frac{1}{2 \beta ^2 \lambda _0^2 t^{2/3}} \left[ \text{ Ai}^2\left( \dfrac{j-v_0t}{ \lambda _0 t^{1/3}}\right) \right. \\&\quad \left. + \text{ Ai}^2\left( -\dfrac{j+v_0t}{ \lambda _0 t^{1/3}}\right) \right] +{\mathcal O}\left( t^{-5/6}\right) \,, \end{aligned} \end{aligned}$$
(1.16)

where \(\text{ Ai }(w) = \frac{1}{\pi } \int _0^\infty \cos (y^3/3+wy) d y\), \(w \in {\mathbb R}\), is the Airy function, and \(\lambda _0 := \frac{1}{2}\Big (\frac{1}{v_0} \sum _{s=1}^m s^4\kappa _s \Big )^{1/3}.\) The above formula is the linear analogue of the Tracy-Widom distribution in (1.11).

Furthermore we can tune the spring intensities \(\kappa _s\), \(s=1,\ldots ,m\) in (1.15) so that we can find an \((m-1)\)-parameter family of potentials such that for \(j\sim \pm v^*t\), with \(0\le v^*<v_0\), one has

$$\begin{aligned} S_{\alpha \alpha '}(j,t)= {\mathcal {O}}\left( \frac{1}{t^{\frac{1}{4}}}\right) ,\;\;\alpha ,\alpha '=1,2\,,\quad S_{33}(j,t)= {\mathcal {O}}\left( \frac{1}{t^{\frac{1}{2}}}\right) ,\quad \text{ as } t\rightarrow \infty \,. \end{aligned}$$

In this case the local behaviour of the correlation functions is described by the Pearcey integral (see Theorem 2.7 and Figs. 2, 3). For example a potential with such behaviour is given by a spring interaction of the form \(\kappa _s=\dfrac{1}{s^2}\) for \(s=1,\ldots , m\) and m even (see Example 2.8 below).

In Sect. 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. Figs. 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 Sect. 2 we study the harmonic oscillator with short range interactions and we introduce the necessary notation and the change of coordinates \(\mathbf{q}\rightarrow \mathbf{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 Sect. 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\rightarrow \infty \). 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.

2 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

$$\begin{aligned} H(\mathbf{p},\mathbf{q})=\sum _{j=0}^{N-1}\frac{p_j^2}{2}+\sum _{s=1}^m\frac{\kappa _s}{2}\sum _{j=0}^{N-1}(q_j-q_{j+s})^2 = \sum _{j=0}^{N-1}\left( \frac{p_j^2}{2} +\frac{1}{2} \Big (\sum _{s=1}^m \tau _s (q_{j+s}-q_j) \Big )^2\right) \end{aligned}$$

so that we may define a Hamiltonian density

$$\begin{aligned} e_j=\frac{p_j^2}{2}+\frac{1}{2} \Big (\sum _{s=1}^m \tau _s (q_{j+s}-q_j) \Big )^2, \end{aligned}$$

which is local in the variables \((\mathbf{p},\mathbf{q})\) for fixed m. Namely, if we let \(N\rightarrow \infty \), the quantity \(e_j\) involves a finite number of physical variables \((\mathbf{p},\mathbf{q})\). Recall that the coefficients \(\tau _s\) are the entries of the circulant localized square root T of the matrix A by which we mean a solution of the equation (1.6) of the form (1.7). 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 \(\mathbf{a}\in {\mathbb R}^N\) is called a circulant matrix generated by the vector \(\mathbf{a}\).

Definition 2.1

(m-physical vector and half-m-physical vector) Fix \(m \in {\mathbb N}\). For any odd \( N > 2m\), a vector \(\tilde{\mathbf{x}}\in {\mathbb R}^{N}\) is said to be m-physical generated by \(\mathbf{x}=(x_0,x_1,\ldots , x_m)\in {\mathbb R}^{m +1}\) if \(x_0=-2\sum _{s=1}^mx_s\) and

$$\begin{aligned} \tilde{x}_0 =&x_0 \, ,\\ \tilde{x}_1=&\tilde{x}_{N-1} = x_1<0 ,\;\;\tilde{x}_m= \tilde{x}_{N-m} = x_m<0,\\ \tilde{x}_k=&\tilde{x}_{N-k} = x_k \le 0, \text{ for } 1< k < m,\\ \tilde{x}_{k} =&0, \text{ otherwise, } \end{aligned}$$

while the vector \(\tilde{\mathbf{x}}\in {\mathbb R}^{N}\) is called half-m-physical generated by \(\mathbf{y}\in {\mathbb R}^{m +1}\) if \(y_0=-\sum _{s=1}^my_s\) and

$$\begin{aligned} \tilde{x}_k=&y_k, \text{ for } 0 \le k \le m \\ \tilde{x}_k=&0, \text{ for } m< k \le N-1. \end{aligned}$$

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 \(\mathbf{a}\) always has a circulant localized square root T that is generated by a half-m physical vector \({\varvec{\tau }}\).

Proposition 2.2

Fix \(m \in {\mathbb N}\). Let the circulant matrix A be generated by an m-physical vector \(\mathbf{a}\), then there exist a circulant matrix T generated by an half-m-physical vector \({\varvec{\tau }}\) such that:

$$\begin{aligned} A = T^\intercal T\,. \end{aligned}$$
(2.1)

Moreover, we can choose \({\varvec{\tau }}\) such that \(\sum _{s=1}^m s \tau _s > 0\). Then one has \(\sum _{s=1}^m s \tau _s=\sqrt{\sum _{s=1}^ms^2\kappa _s}\).

The proof of the proposition is contained in Appendix A:.

For example, if we consider \(m=1\), and \(a_0=2\kappa _1\) and \(a_1=a_{N-1}=-\kappa _1\). The matrix T is generated by the vector \({\varvec{\tau }}=(\tau _0,\tau _1)\) with \(\tau _0=-\sqrt{\kappa _1}\) and \(\tau _1= \sqrt{\kappa _1}\). When \(m=2\) and \(a_0=2\kappa _1+2\kappa _2\), \(a_1=a_{N-1}=-\kappa _1\), \(a_2=a_{N-2}=-\kappa _2\). The matrix T is generated by the vector \({\varvec{\tau }}=(\tau _0,\tau _1,\tau _2)\) with

$$\begin{aligned}&\tau _0=- \frac{\sqrt{\kappa _1}}{2}-\frac{1}{2}\sqrt{\kappa _1+4\kappa _2}, \;\;\tau _1=\sqrt{\kappa _1},\\&\tau _2=-\frac{\sqrt{\kappa _1}}{2}+\frac{1}{2}\sqrt{\kappa _1+4\kappa _2}, \end{aligned}$$

so that the quantities \(r_j\) are defined as

$$\begin{aligned} r_j=\tau _1(q_{j+1}-q_j)+\tau _2 (q_{j+2}-q_j)\,,\quad j\in {\mathbb Z}/ N{\mathbb Z}\,. \end{aligned}$$

Next we integrate the equation of motions. The Hamiltonian \(H(\mathbf{p},\mathbf{q})\) represents clearly an integrable system that can be integrated passing through Fourier transform. Let \({\mathcal F}\) be the discrete Fourier transform with entries \({\mathcal F}_{j,k}:= \frac{1}{\sqrt{N}} e^{- 2\mathrm{i}\pi j k/N}\) with \(j,k=0,\ldots , N-1\). It is immediate to verify that

$$\begin{aligned} {\mathcal F}^{-1} = \bar{{\mathcal F}} \qquad {\mathcal F}^\intercal = {\mathcal F}. \end{aligned}$$
(2.2)

Thanks to the above properties, the transformation defined by

$$\begin{aligned} ({\widehat{p}}, {\widehat{q}}) = (\bar{\mathcal F}p, {\mathcal F}q) \end{aligned}$$
(2.3)

is canonical. Furthermore \(\bar{{\widehat{p}}}_j = {\widehat{p}}_{N-j}\) and \(\bar{{\widehat{q}}}_j = {\widehat{q}}_{N-j}\), for \(j=1,\ldots ,N-1\), while \({\widehat{p}}_0\) and \({\widehat{q}}_0\) are real variables. The matrices T and A are circulant matrices and so they are reduced to diagonal form by \({\mathcal F}\):

$$\begin{aligned} {\mathcal F}A {\mathcal F}^{-1}= {\mathcal F}T^\intercal T{\mathcal F}^{-1}=\overline{({\mathcal F}T{\mathcal F}^{-1})}^\intercal ({\mathcal F}T{\mathcal F}^{-1})\;. \end{aligned}$$

Let \(\omega _j\) denote the eigenvalues of the matrix T ordered so that \({\mathcal F}T{\mathcal F}^{-1}=\) diag\((\omega _j)\). Then \(|\omega _j|^2\) are the (non negative) eigenvalues of the matrix A and

$$\begin{aligned} |\omega _j|^2=\sqrt{N}( \overline{ {\mathcal F}} \tilde{\mathbf{a}})_j,\quad \omega _j=\sqrt{N}( \overline{{\mathcal F}}\tilde{ \varvec{\tau }})_j, \quad j=0,\ldots , N-1, \end{aligned}$$
(2.4)

where \(\tilde{\mathbf{a}}\) is the m-physical vector generated by \(\mathbf{a}\) and \( \tilde{\varvec{\tau }}\) is the half m-physical vector generated by \( \varvec{\tau }\) according to Definition 2.1. It follows that

$$\begin{aligned} \omega _0=0,\quad \omega _j=\overline{\omega }_{N-j}\,,\quad j=1,\ldots ,N-1, \end{aligned}$$
(2.5)

which implies \(|\omega _{j}|^2=|\omega _{N-j}|^2\), \(j=1,\ldots ,N-1\). The Hamiltonian H, can be written as the sum of \(N-1\) oscillators

$$\begin{aligned} H({\widehat{\mathbf{p}}}, {\widehat{\mathbf{q}}}) = \frac{1}{2}\left( \sum _{j=1}^{N-1} |{\widehat{p}}_j |^2 + |\omega _j|^2 |{\widehat{q}}_j|^2 \right) \, = \sum _{j=1}^{\frac{N-1}{2}} |{\widehat{p}}_j |^2 + |\omega _j|^2 |{\widehat{q}}_j|^2 \, . \end{aligned}$$
(2.6)

There are no terms involving \({\widehat{p}}_0, {\widehat{q}}_0\) since the conditions defining \({\mathcal M}\) (1.8) imply that \({\widehat{p}}_0=0\) and \({\widehat{q}}_0=0\). The Hamilton equations are

$$\begin{aligned} {\left\{ \begin{array}{ll} \dfrac{d}{dt}{\widehat{q}}_{j} =\overline{{\widehat{p}}_{j}}\vspace{.2cm}\\ \dfrac{d}{dt}\overline{{\widehat{p}}}_{j} =- |\omega _j|^2{\widehat{q}}_{j}\, .\\ \end{array}\right. } \end{aligned}$$
(2.7)

Thus the general solution reads:

$$\begin{aligned} \begin{aligned}&{\widehat{q}}_{j}(t) = {\widehat{q}}_{j}(0 )\cos (|\omega _{j}|t) + \frac{\overline{{\widehat{p}}_{j}(0)}}{|\omega _j|}\sin (|\omega _j| t)\, , \\&\overline{{\widehat{p}}_{j}(t)} = \overline{{\widehat{p}}_{j}(0)}\cos (|\omega _{j}|t) - |\omega _j|{\widehat{q}}_{j}(0)\sin (|\omega _j| t) \,,\quad j=1,\ldots , N-1, \end{aligned} \end{aligned}$$
(2.8)

and \({\widehat{q}}_{0}(t)=0\) and \({\widehat{p}}_{0}(t)=0\). Inverting the Fourier transform, we recover the variables \(\mathbf{q}={\mathcal F}^{-1}{\widehat{\mathbf{q}}}\), \(\mathbf{p}={\mathcal F}{\widehat{\mathbf{p}}}\) and \(\mathbf{r}={\mathcal F}^{-1}{\widehat{\mathbf{r}}}\) where

$$\begin{aligned} {\widehat{r}}_j =\omega _j{\widehat{q}}_j, \, j = 0, \ldots , N-1\,. \end{aligned}$$
(2.9)

2.1 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 \(\beta ^{-1}\) on the reduced phase space \({\mathcal M}\), namely:

$$\begin{aligned} \mathrm{d}\mu = Z_N(\beta )^{-1} \delta \left( \sum _{k=0}^{N-1} p_k\right) \delta \left( \sum _{k=0}^{N-1} q_k \right) e^{-\beta H(\mathbf{p},\mathbf{q})}\mathrm{d}\mathbf{p}\mathrm{d}\mathbf{q}\end{aligned}$$

where \(Z_N(\beta )\) is the norming constant of the probability measure. For a function \(f=f(\mathbf{p},\mathbf{q})\) we define its average as

$$\begin{aligned} \left\langle f \right\rangle := \int _{{\mathbb R}^{2N}} f(\mathbf{p},\mathbf{q}) \, \, \mathrm{d}\mu . \end{aligned}$$

We first compute all correlation functions (1.10), then we will evaluate the limit \(N\rightarrow \infty \). We first observe that (1.9) in the variables \(({\widehat{\mathbf{p}}}, {\widehat{\mathbf{q}}}) := (\bar{\mathcal F}\mathbf{p}, {\mathcal F}\mathbf{q})\) becomes

$$\begin{aligned} \mathrm{d}\mu = Z_N(\beta )^{-1} \prod _{j=1}^{\frac{N-1}{2}}e^{-\beta (|{\widehat{p}}_j|^2 +|\omega _j|^2|{\widehat{q}}_j|^2 )}\mathrm{d}{\widehat{p}}_j \mathrm{d}{\widehat{q}}_j \end{aligned}$$
(2.10)

where \(\mathrm{d}{\widehat{p}}_j \mathrm{d}{\widehat{q}}_j =\mathrm{d}\mathfrak {R}{\widehat{p}}_j \mathrm{d}\mathfrak {I}{\widehat{p}}_j \mathrm{d}\mathfrak {R}{\widehat{q}}_j \mathrm{d}\mathfrak {I}{\widehat{q}}_j\) and we recall that \({\widehat{p}}_j = \overline{{\widehat{p}}}_{N-j}\), \({\widehat{q}}_j = \overline{{\widehat{q}}}_{N-j}\), \({\widehat{r}}_j =\omega _j{\widehat{q}}_j,\) for \(j=1,\ldots , N-1\).

From the evolution of \({\widehat{p}}_j\) and \({\widehat{q}}_j\) in (2.8) and (2.9), we arrive at the relations

$$\begin{aligned} \left\langle {\widehat{p}}_j(t)\overline{{\widehat{p}}_k(0)} \right\rangle&= \left\langle \overline{{\widehat{p}}_k(0)}\left( {\widehat{p}}_j (0)\cos (|\omega _j|t) - |\omega _j|\overline{{\widehat{q}}_{j}(0)}\sin \left( |\omega _j| t\right) \right) \right\rangle = \delta _{j,k} \frac{1}{\beta }\cos (|\omega _j|t), \end{aligned}$$
(2.11)
$$\begin{aligned} \left\langle {\widehat{p}}_j(t){\widehat{r}}_k(0) \right\rangle&= \left\langle \omega _k {\widehat{q}}_k(0)\left( {\widehat{p}}_j (0)\cos (|\omega _j|t)- |\omega _j |\overline{ \hat{q}_{j}(0)}\sin \left( |\omega _j| t\right) \right) \right\rangle \nonumber \\&= -\delta _{j,k} \frac{\omega _j}{|\omega _j|\beta }\sin \left( |\omega _j| t\right) , \end{aligned}$$
(2.12)
$$\begin{aligned} \left\langle {\widehat{r}}_j(t){\widehat{p}}_k(0)\right\rangle&= \left\langle \omega _j{\widehat{p}}_k(0)\left( {\widehat{q}}_j (0)\cos (|\omega _j|t) + \frac{\overline{{\widehat{p}}_{j}(0)}}{|\omega _j|}\sin \left( |\omega _j| t\right) \right) \right\rangle = \delta _{j,k}\frac{\omega _j}{|\omega _j|\beta }\sin (|\omega _j|t) \end{aligned}$$
(2.13)
$$\begin{aligned} \left\langle {\widehat{r}}_j(t)\overline{{\widehat{r}}_k(0)}\right\rangle&= \left\langle \overline{\omega }_k\omega _j\overline{{\widehat{q}}_k(0)}\left( {\widehat{q}}_j (0)\cos (|\omega _j|t) + \frac{\overline{{\widehat{p}}_{j}(0)}}{|\omega _j|}\sin \left( |\omega _j| t\right) \right) \right\rangle = \delta _{j,k} \frac{1}{\beta }\cos (|\omega _j|t). \end{aligned}$$
(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}(j,t) \), and we leave to the reader the details for the other cases:

$$\begin{aligned} \begin{aligned} S^{N}_{11}(j,t)&= \left\langle r_{j}(t)r_0(0) \right\rangle = \frac{1}{N}\left\langle \sum _{k,l=1}^{N-1}{\widehat{r}}_k(t){\widehat{r}}_l(0)e^{2\pi \imath \frac{jk}{N}}\right\rangle \\&= \frac{1}{N\beta } \sum _{l=1}^{N-1}\cos \left( |\omega _l|t\right) \cos \left( 2\pi \frac{lj}{N}\right) =S^{N}_{22}(j,t)\,. \end{aligned} \end{aligned}$$
(2.15)

In the same way we have that:

$$\begin{aligned}&S^{N}_{12}(j,t)= \frac{1}{N\beta }\sum _{l=1}^{N-1}\sin (|\omega _l|t)\cos \left( 2\pi \frac{lj}{N} + \arg (\omega _l) \right) \end{aligned}$$
(2.16)
$$\begin{aligned}&S^{N}_{21}(j,t)= - \frac{1}{N\beta }\sum _{l=1}^{N-1}\sin (|\omega _l|t)\cos \left( 2\pi \frac{lj}{N}- \arg (\omega _l) \right) \end{aligned}$$
(2.17)
$$\begin{aligned}&S^{N}_{31}(j,t)= S^{N}_{32}(j,t) = S^{N}_{13}(j,t) = S^{N}_{23}(j,t) = 0 \end{aligned}$$
(2.18)
$$\begin{aligned}&S^{N}_{33}(j,t) =\frac{1}{2}((S^N_{11})^2 + (S^N_{22})^2 + (S^N_{12})^2+(S^N_{21})^2)+ \frac{3(N-1)}{2N^2\beta ^2}\,. \end{aligned}$$
(2.19)

The dispersion relation given by (2.4) takes the form

$$\begin{aligned} \begin{aligned}&\omega _\ell =-\sum _{s=1}^m\tau _s\left( 1-\cos \left( 2\pi \frac{s\ell }{N}\right) \right) +\mathrm{i}\sum _{s=1}^m\tau _s\sin \left( 2\pi \frac{s\ell }{N}\right) \\&|\omega _\ell |^2= \sum _{s=0}^{N-1}a_se^{-2\pi \mathrm{i}\frac{s\ell }{N}}= 2\sum _{s=1}^m\kappa _s\left( 1-\cos \left( 2\pi \frac{s\ell }{N}\right) \right) \,, \end{aligned} \end{aligned}$$
(2.20)

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 \(\omega (k)\) to provide continuum limits of \(\omega _\ell \) and \(|\omega _\ell |^2\), namely

$$\begin{aligned} \begin{aligned} \omega (k)&:=-\sum _{s=1}^m\tau _s\left( 1-\cos \left( 2\pi sk\right) \right) +\mathrm{i}\sum _{s=1}^m\tau _s\sin \left( 2\pi sk\right) \\ |\omega (k)|^2&= 2\sum _{s=1}^m\kappa _s\left( 1-\cos (2\pi k s)\right) \,, \end{aligned} \end{aligned}$$
(2.21)

where the variable \(\ell /N\) has been approximated with \(k\in [0,1]\). One may use equation (A.1) to check the consistency of the two equations of (2.21). To this end observe that \(\omega (k)=Q(e^{-2\pi ik})\), \(\overline{\omega (k)}=Q(e^{2\pi ik})\), and \(|\omega (k)|^2=\ell (e^{2\pi ik})\).

Lemma 2.3

Let \(\omega (k)\) be defined as in (2.21), set \(f(k):=|\omega (k)|\), and denote \(\theta (k):=\) arg\((\omega (k))\) for \(0\le k \le 1\), where the ambiguity in the definition of \(\theta \) is settled by requiring \(\theta \) to be continuous with \(\theta (0) \in (-\pi , \pi ]\). Then, for all \(k\in [0,1]\) we have

$$\begin{aligned}&\omega (1-k)=\overline{\omega (k)}, \end{aligned}$$
(2.22)
$$\begin{aligned}&f(1-k)=f(k), \end{aligned}$$
(2.23)
$$\begin{aligned}&\theta (1-k)\equiv - \theta (k)\quad (\text{ mod } \;\;2\pi ). \end{aligned}$$
(2.24)

Furthemore, the functions f and \(\theta -\frac{\pi }{2}\) are \(C^\infty \) on [0, 1] and they both possess odd \(C^\infty \)-extensions at \(k=0\) which implies in particular \(\;\theta (0)=\frac{\pi }{2}\,\).

Proof

The symmetries follow directly from the definition of \(\omega \) in (2.21). From (2.21) we also learn that \(|\omega (k)|^2 \ge 2 \kappa _1 (1-\cos (2\pi k)) > 0\) for \(k\in (0, 1)\). Thus the smoothness of f and \(\theta \) only needs to be investigated for \(k\in \{0,1\}\). By symmetry we only need to study the case \(k=0\). The smoothness of the function \(\theta \) may be obtained from the expansion near \(k=0\)

$$\begin{aligned} \cot (\theta (k)) = -k\pi \frac{\sum _{s=1}^m s^2\tau _s}{\sum _{s=1}^m s \tau _s} + {\mathcal O}(k^3) \end{aligned}$$

together with \(\sum _{s=1}^m s \tau _s > 0\) (see Proposition 2.2). Since \(\cot (\theta (0))=0\) and \(\mathfrak {I}\omega (k) > 0\) for small positive values of k we conclude that \(\theta (0)=\frac{\pi }{2}\) from the requirement \(\theta (0) \in (-\pi , \pi ]\). This also implies the existence of a smooth odd extension of \(\theta -\frac{\pi }{2}\) at \(k=0\) because \(\cot (\theta (k))\) has such an extension. For the function f the claims follow from the representation

$$\begin{aligned} f(k)=2\pi k \left( \sum _{s=1}^m s^2\kappa _s \text{ sinc}^2(\pi s k)\right) ^{1/2} \end{aligned}$$

near \(k=0\) where sinc\((x)=\frac{\sin (x)}{x}\) denotes the smooth and even sinus cardinalis function.

\(\square \)

Lemma 2.4

In the limit \(N\rightarrow \infty \) the correlation functions have the following expansion

$$\begin{aligned}&S^{N}_{\alpha \alpha '}(j,t)+\frac{\delta _{\alpha \alpha '}}{N\beta }= S_{\alpha \alpha '}(j,t)+ {\mathcal {O}} \left( N^{-\infty }\right) ,\quad \alpha , \alpha '=1,2,\\&S^{N}_{33}(j,t)= S_{33}(j,t)+ {\mathcal {O}} \left( N^{-1}\right) , \end{aligned}$$

where \(\delta _{\alpha \alpha '}\) denotes the Kronecker delta,

$$\begin{aligned} S_{11}(j,t)&=S_{22}(j,t) =\frac{1}{\beta }\int _{0}^1\cos \left( |\omega (k)|t\right) \cos \left( 2\pi kj\right) \mathrm{d}k \end{aligned}$$
(2.25)
$$\begin{aligned} S_{12}(j,t)&=\frac{1}{\beta }\int _{0}^1\sin \left( |\omega (k)|t\right) \cos \left( 2\pi kj +\theta (k)\right) \mathrm{d}k, \end{aligned}$$
(2.26)
$$\begin{aligned} S_{21}(j,t)&=-\frac{1}{\beta }\int _{0}^1\sin \left( |\omega (k)|t\right) \cos \left( 2\pi kj -\theta (k)\right) \mathrm{d}k, \end{aligned}$$
(2.27)
$$\begin{aligned} S_{33}(j,t)&=\frac{1}{2}(S_{11}^2 + S_{22}^2 + S_{12}^2+S_{21}^2), \end{aligned}$$
(2.28)

and \(\theta (k)=\arg \omega (k)\) with \(\omega (k)\) as in (2.21).

Proof

For any periodic \(C^\infty \)-function g on the real line with period 1, \(g(k)=\sum _{n\in {\mathbb Z}}\hat{g}_n e^{2\pi ikn}\), one has

$$\begin{aligned} \frac{1}{N} \sum _{\ell =0}^{N-1} g\left( \frac{\ell }{N}\right) = \sum _{m\in {\mathbb Z}} \hat{g}_{mN} = \int _0^1 g(k) \mathrm{d}k + {\mathcal {O}} \left( N^{-\infty }\right) \,. \end{aligned}$$

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

$$\begin{aligned} \cos \left( f(-k)t\right) \cos \left( -2\pi kj\right)= & {} \cos \left( f(k)t\right) \cos \left( -2\pi kj\right) \;\;\\= & {} \;\; \cos \left( f(1-k)t\right) \cos \left( 2\pi (1- k)j\right) \,,\\ \sin \left( f(-k)t\right) \cos \left( -2\pi kj \pm \theta (-k)\right)= & {} -\sin \left( f(k)t\right) \cos \left( -2\pi kj \pm (\pi -\theta (k))\right) \\= & {} \sin \left( f(1-k)t\right) \cos \left( 2\pi (1-k)j \pm \theta (1-k)\right) \,. \end{aligned}$$

Observing in addition that the summands corresponding to \(\ell =0\) are missing in (2.15)-(2.17) the first claim is proved. Together with (2.19) this also implies the second claim. \(\square \)

Next we analyse the leading order behaviour (as \(t \rightarrow \infty \)) of the limiting correlation functions \(S_{\alpha \alpha '}(j,t)\) using the method of steepest descent. In order to explain the phenomena that may occur we start by discussing \(S_{11}\). Denote

$$\begin{aligned} \xi :=\dfrac{j}{t} \quad \text{ and } \quad \phi _\pm (k, \xi ) := f(k) \pm 2 \pi \xi k\,. \end{aligned}$$
(2.29)

With these definitions and using the symmetry (2.23) we may write

$$\begin{aligned} S_{11}(j,t)= \frac{1}{2\beta }\mathfrak {R}\int _{0}^{1}\left( e^{it(f(k)+2\pi \xi k)} +e^{it(f(k)-2\pi \xi k)}\right) d k = \frac{1}{\beta }\mathfrak {R}\int _{0}^{1} e^{it\phi _- (k, \xi )} d k\,.\nonumber \\ \end{aligned}$$
(2.30)

The leading order behaviour (\(t \rightarrow \infty \)) of such an integral is determined by the stationary phase points \(k_0 \in [0, 1]\), i.e. by the solutions of the equation \(\frac{\partial }{\partial k} \phi _-(k_0, \xi ) = 0\) which depend on the value of \(\xi \).

Such stationary phase points do not need to exist. In fact, as we see in Lemma 2.5(b) below, the range of \(f'\) is given by some interval \([-2\pi v_0, 2 \pi v_0]\) so that there are no stationary phase points for \(| \xi | > v_0\). As in the proof of Lemma 2.4 one can argue that the integrand \(\mathfrak {R}e^{it\phi _- (k, j/t)}\) 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}(j,t)\) decays rapidly in time. More precisely, for every fixed \(\delta > 0\) we have

$$\begin{aligned} S_{11}(j,t)= {\mathcal O}\left( t^{-\infty }\right) \quad \text{ as } t \rightarrow \infty , \text{ uniformly } \text{ for } |j|\ge (v_0 + \delta )t. \end{aligned}$$
(2.31)

This justifies the name of sound speed for the quantity \(v_0\).

In the case \(| \xi | \le v_0\) there always exists at least one stationary phase point \(k_0=k_0(\xi ) \in [0, 1]\). Each stationary phase point may provide an additive contribution to the leading order behaviour of \(\int _{0}^{1} e^{it\phi _- (k, \, j/t)} d k\) for j near \(\xi 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 \(\phi _- (\cdot , \xi )\), i.e. \(\frac{\partial }{\partial k} \phi _- (k_0, \xi ) = 0\). Denote by \(\ell \) the smallest integer bigger than 1 for which \(\frac{\partial ^\ell }{\partial k^\ell } \phi _-(k_0, \xi ) \ne 0\). Then \(k_0\) contributes a term of order \(t^{1/\ell }\) to the t-asymptotics of \(\int _{0}^{1} e^{it\phi _- (k, \, j/t)} d k\) for j in a suitable neighbourhood of \(\xi t\).

Before treating the general situation let us recall the case of nearest-neighbour interactions. There we have

$$\begin{aligned} f(k) = f_{1}(k) =\sqrt{2 \kappa _1 (1- \cos (2\pi k))} = 2 \sqrt{\kappa _1} \sin (\pi k)\,, \quad k \in [0, 1]\,. \end{aligned}$$

The range of \(f_{1}'\) equals \([-2\pi v_0, 2 \pi v_0]\) with \(v_0= \sqrt{\kappa _1}\). For every \(|\xi | \le v_0\) there exists exactly one stationary phase point \(k_0(\xi )\in [0, 1]\) of \(\phi _- (\cdot , \xi )\) that is determined by the relation \(\cos (\pi k_0(\xi )) = \xi /v_0\). A straight forward calculation gives

$$\begin{aligned} \frac{\partial ^2}{\partial k^2} \phi _-(k_0(\xi ), \xi ) = f_{1}''(k_0(\xi )) = -2 \pi ^2 \sqrt{v_0^2-\xi ^2} = 0 \; \Leftrightarrow \; \xi =\pm v_0\,. \end{aligned}$$

Moreover, we have \(k_0(v_0)=0\) and \(k_0(-v_0)=1\) and therefore \(\frac{\partial ^3}{\partial k^3} \phi _-(k_0(\pm v_0), \pm v_0) = \mp 2 \pi ^3 v_0 \ne 0\). This implies that in addition to (2.31) we have \(S_{11}(j,t)= {\mathcal O}(t^{-1/2})\), except for j near \(\pm v_0 t\) where \(S_{11}(j,t)= {\mathcal O}(t^{-1/3})\). In order to determine the behaviour near the least decaying peaks that travel at speeds \(\pm v_0\) we expand \(f_{1}\) near the stationary phase points. Let us first consider \(\xi =v_0\) with \(k_0=0\). Introducing \(\lambda _0=\frac{1}{2\pi }|f_{1}'''(0)/2|^{1/3}=\frac{1}{2}v_0^{1/3}\) we obtain

$$\begin{aligned} f_{1}(k)=2\pi v_0 k - \frac{1}{3} (2 \pi \lambda _0 k)^3 + {\mathcal O}( k^5 )\,, \quad \text{ as } k \rightarrow 0. \end{aligned}$$

Substituting \(y = 2 \pi \lambda _0 t^{1/3} k\) leads for k close to 0 to the asymptotic expression

$$\begin{aligned} t \phi _- (k, \, j/t)= \frac{v_0 t-j}{\lambda _0 t^{1/3}} y - \frac{1}{3} y^3 + {\mathcal O}( t^{-2/3} )\,, \quad \text{ as } t \rightarrow \infty . \end{aligned}$$

Using the well-known representation \(\text{ Ai }(w) = \frac{1}{\pi } \int _0^\infty \cos (y^3/3+wy) d y\), \(w \in {\mathbb R}\), of the Airy function and performing a similar analysis around the stationary phase point \(k_0=-1\) for \(\xi =-v_0\) one obtains an asymptotic formula for the region not covered by (2.31)

$$\begin{aligned} S_{11}(j,t)= & {} \frac{1}{2 \beta \lambda _0 t^{1/3}} \left[ \text{ Ai }\left( \dfrac{j-v_0t}{ \lambda _0 t^{1/3}}\right) + \text{ Ai }\left( -\dfrac{j+v_0t}{ \lambda _0 t^{1/3}}\right) \right] \nonumber \\&+{\mathcal O}\!\!\left( t^{-1/2}\right) ,t \rightarrow \infty , \text{ uniformly } \text{ for } |j|<(v_0 + \delta )t \end{aligned}$$
(2.32)

for \(\delta >0\) (see e.g. [14]). Observe that due to the decay of Ai(w) for \(w \rightarrow \pm \infty \), the Airy term is dominant roughly in the regions described by \(v_0 t - o (t)<|j|<v_0 t + o ((\ln t)^{2/3})\).

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}\):

$$\begin{aligned} f''(k)< 0 \quad \text{ for } \text{ all } \quad 0 < k \le \frac{1}{2}\,, \end{aligned}$$
(2.33)

together with

$$\begin{aligned} \quad f''(0) = 0\,,\quad f'''(0)< 0\,, \quad \text{ and } \quad f(1-k)=f(k)\quad \text{ for } \text{ all } \quad 0 \le k < \frac{1}{2}\,. \end{aligned}$$
(2.34)

Conditions (2.33) and (2.34) imply that statements (2.31) and (2.32) hold with \(v_0=\frac{f'(0)}{2 \pi } > 0\) and \(\lambda _0=\frac{1}{2\pi }|f'''(0)/2|^{1/3}\).

It follows from equation (2.23) and from statement (a) of Lemma 2.5 below that the conditions of (2.34) are always satisfied in our model. Condition (2.33), however, might fail. Indeed, it is not hard to see that there exist open regions in the \({\varvec{\kappa }}\)-space \({\mathbb R}_+^m\) where there always exist stationary phase points \(k_0 \in (0, 1) \) of higher multiplicity, i.e. with \(f''(k_0)=0\). In this situation the value of \(v:= \frac{f'(k_0)}{2 \pi }\) lies in the open interval \((-v_0, v_0)\) (cf. Lemma 2.5b). Then the decay rate of \(S_{11}(j,t)\) 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'''(k_0)=0\) holds in addition. We show in Theorem 2.7 that this may happen for \({\varvec{\kappa }}\) in some submanifold of \({\mathbb R}_+^m\) of codimension 1 (see also Examples 2.8 and 2.9). Nevertheless, if \(\kappa _2\), \(\ldots \), \(\kappa _m\) are sufficiently small in comparison to \(\kappa _1\) then condition (2.33) is always satisfied as we show in Theorem 2.6(c).

Before stating our main results of this section, Theorems 2.6 and 2.7, we first summarize some more properties of the function f.

Lemma 2.5

Given \((\kappa _1,\ldots ,\kappa _m)\) with \(\kappa _1>0\), \(\kappa _m>0\), and \(\kappa _j\ge 0\) for \(1<j<m\). Denote \(f(k)=|\omega (k)|\) for \(0\le k \le 1\) as introduced in Lemma 2.3 and define \(v_0:=(\sum _{s=1}^m s^2\kappa _s)^{\frac{1}{2}}\). Then the following holds:

  1. (a)

    \(\;f(0)=f''(0)=0\,\), \(\;f'(0)=2\pi v_0\,\), and \(\;f'''(0)=-\frac{2 \pi ^3}{v_0}\sum _{s=1}^m s^4\kappa _s\,\).

  2. (b)

    \(f'([0, 1]) = [-2\pi v_0,2\pi v_0]\). \(f'\) attains its maximum only at \(k=0\) and its minimum only at \(k=1\).

  3. (c)

    Fix \(\kappa _1 >0\). Then the map f can be extended as a \(C^\infty \)-function of the variables \((k,\kappa _2,\ldots ,\kappa _m)\) on the set \([0,1]\times [0, \infty )^{m-1}\).

Proof

Statement (a) follows directly from the last formula in the proof of Lemma 2.3 and from the expansion \(\hbox {sinc}^2(x) = 1 -\frac{x^2}{3} +{\mathcal O}(x^4)\) for small values of x:

$$\begin{aligned} f(k)=2\pi k \left( \sum _{s=1}^m s^2\kappa _s \text{ sinc}^2(\pi s k)\right) ^{1/2}=\; 2 \pi v_0 k-\frac{\pi ^3}{3v_0}\left( \sum _{s=1}^m s^4\kappa _s\right) k^3 + {\mathcal O}(k^5)\,. \end{aligned}$$

This representation also settles statement (c). As we know already \(f'(0)=2\pi v_0=-f'(1)\) we may establish statement (b) by verifying that \(|f'(k)|<2\pi v_0\) holds for all \(k\in (0, 1)\). To this end we write \(f= (\sum _{s=1}^m h_s^2)^{1/2}\) with \(h_s(k) = 2\sqrt{\kappa _s}\sin (\pi s k)\). Using the Cauchy-Schwarz inequality we obtain for \(0<k<1\) that

$$\begin{aligned} |f'(k)|= & {} \frac{|\sum _{s=1}^m h_s(k) h_s'(k)|}{\left( \sum _{s=1}^m h^2_s(k)\right) ^{1/2}} \le \left( \sum _{s=1}^m (h'_s)^2(k)\right) ^{1/2} \\&= 2\pi \left( \sum _{s=1}^m s^2\kappa _s \cos ^2(\pi s k)\right) ^{1/2} < \; 2\pi v_0 \,, \end{aligned}$$

where the last inequality follows from \(|\cos (\pi k)|<1\) and \(\kappa _1 > 0\). \(\square \)

We are now ready to state our first main result in this section.

Theorem 2.6

Let \(m \in {\mathbb N}\), fix \(\delta >0\), denote \(f(k)=|\omega (k)|\) as introduced in Lemma 2.3, and set

$$\begin{aligned} v_0 := \sqrt{\sum _{s=1}^m s^2\kappa _s}, \quad \lambda _0 := \frac{1}{2}\left( \frac{1}{v_0} \sum _{s=1}^m s^4\kappa _s \right) ^{1/3}. \end{aligned}$$
(2.35)
  1. (a)

    For all \(\alpha \), \(\alpha '=1,2,3\) we have rapid decay as \(t \rightarrow \infty \), uniformly for \(|j|>(v_0 + \delta )t\), i.e.

    $$\begin{aligned} S_{\alpha \alpha '}(j,t)= {\mathcal O}\left( t^{-\infty }\right) . \end{aligned}$$
  2. (b)

    If \(f''(k)<0\) for all \(0<k\le 1/2\) then as \(t \rightarrow \infty \) the following holds uniformly for \(|j|<(v_0 + \delta )t\):

    $$\begin{aligned} S_{11}(j,t)= & {} \frac{1}{2 \beta \lambda _0 t^{1/3}} \left[ \text{ Ai }\left( \dfrac{j-v_0t}{ \lambda _0 t^{1/3}}\right) + \text{ Ai }\left( -\dfrac{j+v_0t}{ \lambda _0 t^{1/3}}\right) \right] +{\mathcal O}\left( t^{-1/2}\right) \;= \; S_{22}(j,t)\,, \end{aligned}$$
    (2.36)
    $$\begin{aligned} S_{12}(j,t)= & {} \frac{1}{2\lambda _0t^{1/3}\beta }\left( \text{ Ai }\left( -\frac{j+v_0t}{\lambda _0 t^{1/3 } }\right) -\text{ Ai }\left( \frac{j-v_0t}{\lambda _0 t^{1/3 } }\right) \right) +{\mathcal O}(t^{-\frac{1}{2}})= \; S_{21}(j,t)\,, \end{aligned}$$
    (2.37)
    $$\begin{aligned} S_{33}(j,t)= & {} \frac{1}{2 \beta ^2 \lambda _0^2 t^{2/3}} \left[ \text{ Ai}^2\left( \dfrac{j-v_0t}{ \lambda _0 t^{1/3}}\right) + \text{ Ai}^2\left( -\dfrac{j+v_0t}{ \lambda _0 t^{1/3}}\right) \right] +{\mathcal O}\left( t^{-5/6}\right) \,. \end{aligned}$$
    (2.38)
  3. (c)

    For every \(\kappa _1>0\) there exists \(\varepsilon = \varepsilon (\kappa _1) > 0\) such that for all \((\kappa _2, \ldots , \kappa _m) \in [0, \varepsilon )^{m-1}\) we have \(f''(k) < 0\) for all \(0<k\le 1/2\).

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\in {\mathbb R}\) that are 1-periodic. This is all that is needed because \(|\frac{\partial }{\partial k} \phi _\pm (k,\,j/t)|>2\pi \delta \) by Lemma 2.5b) uniformly for \(k\in [0,1/2]\) and \(|j|>(v_0 + \delta )t\).

We have already argued above that conditions (2.33), (2.34) suffice to derive the first claim of statement (b) with \(v_0=\frac{f'(0)}{2 \pi } > 0\) and \(\lambda _0=\frac{1}{2\pi }|f'''(0)/2|^{1/3}\). The expressions for \(f'(0)\) and \(f'''(0)\) stated in Lemma 2.5 (a) justify the definitions of (2.35).

Using the symmetry relations (2.23) and (2.24) we derive a representation for \(S_{12}\) and \(S_{21}\) that is suitable for a steepest descent analysis

$$\begin{aligned} \begin{aligned} S_{12}(j,t) =&\frac{1}{\beta }\int _{0}^{1/2}\Big (\sin (f(k)t - 2\pi kj-\theta (k)) + \sin (f(k)t + 2\pi kj+\theta (k))\Big ) \mathrm{d}k \\ =&\frac{1}{\beta }\mathfrak {I}\int _{0}^{1/2} \Big (e^{it\phi _- (k,\,j/t)}e^{-i\theta (k)} + e^{it\phi _+ (k,\,j/t)} e^{i\theta (k)} \Big )d k \\ S_{21}(j,k)=&- \frac{1}{\beta }\mathfrak {I}\int _{0}^{1/2} \Big (e^{it\phi _- (k,\,j/t)}e^{i\theta (k)} + e^{it\phi _+ (k,\,j/t)} e^{-i\theta (k)} \Big )d k \\ \end{aligned} \end{aligned}$$
(2.39)

where \(\phi _\pm (k,\xi ) =f(k)\pm 2\pi \xi k\) as in (2.29) above. Expanding for k close to zero one obtains \(\phi _\pm (k,j/t)=2\pi v_0 k-\frac{1}{3}(2\pi )^3\lambda _0^3k^3\pm 2\pi k\frac{j}{t} +{\mathcal O}(k^5)\). Substituting \(y = 2 \pi \lambda _0 t^{1/3} k\) leads to the asymptotic expression

$$\begin{aligned} t\phi _\pm (k,j/t)=\frac{v_0t\pm j}{\lambda _0t^{\frac{1}{3}}}y-\frac{1}{3}y^3 + {\mathcal O}( t^{-\frac{2}{3}})\, \quad \text{ as } t\rightarrow \infty . \end{aligned}$$

Keeping in mind that \(\theta (0)=\frac{\pi }{2}\) we obtain

$$\begin{aligned} S_{12}(j,t)= & {} \frac{1}{2\lambda _0t^{1/3}\beta }\left( \text{ Ai }\left( -\frac{j+v_0t}{\lambda _0 t^{1/3 } }\right) -\text{ Ai }\left( \frac{j-v_0t}{\lambda _0 t^{1/3 } }\right) \right) +{\mathcal O}(t^{-\frac{1}{2}})=S_{21}(j,t)\,. \end{aligned}$$

Regarding the expansion for \(t\rightarrow \infty \) of \(S_{33}(j,t)\) it follows immediately from the expression (2.28) and the expansions of \(S_{\alpha \alpha '}(j,t)\) with \(\alpha ,\alpha '=1,2\).

Statement (c) follows from the continuous dependence of the derivatives \(f''\) and \(f'''\) on the parameters \((\kappa _2, \ldots , \kappa _m)\) (see Lemma 2.5 c) and from simple facts for the case of nearest-neighbour interactions \(f_{1}(k)=2\sqrt{\kappa _1}\sin (\pi k)\) discussed above. Indeed, from \(f''(0)=0\) and from \(f'''_{1}(0) < 0\) it follows that there exists such an \(\varepsilon >0\) such that \(f'''(k) < 0\) and hence also \(f''(k) < 0\) for k in some region \((0, \delta )\) uniformly in \((\kappa _2, \ldots , \kappa _m) \in [0, \varepsilon )^{m-1}\). As \(f''_{1}(k) < -2\pi ^2\sqrt{\kappa _1}\sin (\pi \delta )\) for all \(k \in [\delta , 1/2]\) we may prove the claim in this region by reducing the value of \(\varepsilon \) if necessary. \(\square \)

Theorem 2.6 provides the leading order asymptotics of the limiting correlations \(S_{\alpha \alpha '}(j,t)\) for \(t\rightarrow \infty \) in the simple situation that the second derivative of the dispersion relation is strictly negative on the open interval (0, 1) (cf. condition (2.33)). Moreover, statement (c) shows that there is a set of positive measure in parameter space \({\varvec{\kappa }}\in {\mathbb R}_+^m\) where this happens. For general values of \({\varvec{\kappa }}\), 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/(3+p)}\) 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_{\alpha \alpha '}(j,t)\). Despite being non-generic in parameter space it is interesting to note that decay rates \(t^{-1/4}\) can be observed numerically (see Figs. 2 and 3). There is also a second issue that may arise if condition (2.33) fails. Namely, for \(v \in (-v_0, v_0)\) there can be several values of \(k \in (0,\frac{1}{2}]\) satisfying \(f'(k)\pm 2\pi 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.

Theorem 2.7

Recall from (2.21) the formula for the dispersion relation

$$\begin{aligned} f(k)=|\omega (k)| =\sqrt{2\sum _{s=1}^m\kappa _s\left( 1-\cos (2\pi k s)\right) }\,. \end{aligned}$$

(a) For \(m\ge 3\) there is an \((m-1)\)-parameter family of potentials for which there exists \(k^*=k^*({\varvec{\kappa }})\in (0,\frac{1}{2})\) with

$$\begin{aligned} f''(k^*)=0,\;\;f'''(k^*)=0,\;\;f^{(iv)}(k^*)\ne 0,\;\; \text{ and } \;\;0< v^*:=\frac{f'(k^*)}{2\pi }<v_0, \end{aligned}$$
(2.40)

with \(v_0\) as in (2.35). Set \(\lambda ^{*}:=\dfrac{1}{2\pi }(|f^{(iv)}(k^*)|/4!)^{\frac{1}{4}}>0\). Then for \(j\rightarrow \infty \) and \(t\rightarrow \infty \) in such a way that

$$\begin{aligned} \frac{v^*t - j}{\lambda ^{*}t^{\frac{1}{4}}} \end{aligned}$$

is bounded, the contribution of the stationary phase point \(k^*\) to the correlation functions is given by:

$$\begin{aligned} S_{11}(j,t),S_{22}(j,t)\;:&\quad \dfrac{1}{2\beta \pi \lambda ^{*}t^{\frac{1}{4}}} \mathfrak {R}\left( e^{it\phi _-(k^*,j/t)}{{\mathcal {P}}}_\pm \left( \frac{v^*t-j}{\lambda ^{*}t^{\frac{1}{4}}}\right) \right) +{\mathcal O}(t^{-\frac{1}{2}})\,, \end{aligned}$$
(2.41)
$$\begin{aligned} S_{12}(j,t) \;:&\quad \dfrac{1}{2\beta \pi \lambda ^{*}t^{\frac{1}{4}}} \mathfrak {I}\left( e^{it\phi _-(k^*,j/t)-i\theta (k^*)}{{\mathcal {P}}}_\pm \left( \frac{v^*t-j}{\lambda ^{*}t^{\frac{1}{4}}}\right) \right) + {\mathcal O}(t^{-\frac{1}{2}})\,, \end{aligned}$$
(2.42)
$$\begin{aligned} S_{21}(j,t) \;:&\quad -\dfrac{1}{2\beta \pi \lambda ^{*}t^{\frac{1}{4}}} \mathfrak {I}\left( e^{it\phi _-(k^*,j/t)+i\theta (k^*)}{\mathcal P}_\pm \left( \frac{v^*t-j}{\lambda ^{*}t^{\frac{1}{4}}}\right) \right) +{\mathcal O}(t^{-\frac{1}{2}})\,, \end{aligned}$$
(2.43)

where \(\phi _\pm (k,\xi ) =f(k)\pm 2\pi \xi k\), \(\theta (k)=\arg \omega (k)\) as defined in Lemma 2.3, \({\mathcal P}_\pm (a)\) denote the Pearcey integrals, cf. Appendix B:,

$$\begin{aligned} {{\mathcal {P}}}_\pm (a)=\int _{-\infty }^{\infty }e^{i(\pm y^4+ay)}dy,\;\; \;\;a\in {\mathbb {R}}, \end{aligned}$$
(2.44)

and \( {{\mathcal {P}}}_\pm \) has to be chosen according to the sign of \(f^{(iv)}(k^*)\). If \(j\rightarrow -\infty \) with bounded \((v^*t+j)/(\lambda ^{*} t^{1/4})\) the contributions of the stationary point \(k^*\) can be obtained from the ones presented in (2.41)-(2.43) by replacing \(\phi _-\) by \(\phi _+\), \(\theta \) by \(-\theta \), and j in the argument of \({{\mathcal {P}}}_\pm \) by \(-j\).

(b) When \(k^*=\frac{1}{2}\) one has \(f'(1/2)=0\) and \(f'''(1/2)=0\) by the symmetry (2.23). For each \(m \ge 2\) there is an \((m-1)\)-parameter family of potentials so that \(f''(1/2)=0\) and \(f^{(iv)}(1/2)\ne 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\rightarrow \infty \) with bounded \(j/t^{\frac{1}{4}}\) is given by (\(\lambda ^{*}\) defined as in statement (a) with \(k^*=\frac{1}{2}\))

$$\begin{aligned} \begin{aligned} S_{11}(j,t), S_{22}(j,t) \;:&\quad \frac{(-1)^j}{2 \beta \pi \lambda ^{*}t^{\frac{1}{4}}}\mathfrak {R}\left( e^{i t f(\frac{1}{2})}{{\mathcal {P}}}_\pm \left( \frac{j}{\lambda ^{*}t^{\frac{1}{4}}}\right) \right) +{\mathcal O}(t^{-\frac{1}{2}})\\ S_{12}(j,t), S_{21}(j,t) \;:&\quad -\text{ sgn }\,(\!\sum \limits _{s \,odd}\tau _s)\frac{(-1)^j}{2 \beta \pi \lambda ^{*}t^{\frac{1}{4}}}\mathfrak {I}\left( e^{i t f(\frac{1}{2})}{{\mathcal {P}}}_\pm \left( \frac{j}{\lambda ^{*}t^{\frac{1}{4}}}\right) \right) +{\mathcal O}(t^{-\frac{1}{2}})\\ S_{33}(j,t) \;:&\quad \frac{1}{4 \beta ^2\pi ^2(\lambda ^{*})^2t^{\frac{1}{2}}}\left| {\mathcal P}_\pm \left( \frac{j}{\lambda ^{*}t^{\frac{1}{4}}}\right) \right| ^2+{\mathcal O}(t^{-\frac{3}{4}})\,. \end{aligned}\nonumber \\ \end{aligned}$$
(2.45)

Proof

We begin by proving formula (2.41) for the momentum or position correlations \(S_{22}(j,t)=S_{11}(j,t)\) under the assumption that we have found a \(k^{*} \in (0,1/2)\) for which all the relations of (2.40) are satisfied.

From (2.30) and Lemma 2.3 we obtain

$$\begin{aligned}&S_{11}(j,t)=S_{22}(j,t) =\frac{1}{\beta }\mathfrak {R}\int _{0}^{\frac{1}{2}}\left( e^{it(f(k)+2\pi k\frac{j}{t})} +e^{it(f(k)-2\pi k\frac{j}{t})}\right) d k. \end{aligned}$$
(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

$$\begin{aligned} f(k)=f(k^*)+2\pi v^*(k-k^*)+f^{(iv)}(k^*)(k-k^*)^4/4!+O((k-k^*)^5)\,. \end{aligned}$$

Introducing the change of variables

$$\begin{aligned} y=2\pi \lambda ^{*}(k-k^*)t^{\frac{1}{4}},\quad \lambda ^{*}=\dfrac{1}{2\pi }(|f^{(iv)}(k^*)|/4!)^{\frac{1}{4}} \end{aligned}$$

one obtains

$$\begin{aligned} tf(k)-2\pi j k=tf(k^*)-2\pi jk^*+y\frac{v^*t-j}{\lambda ^{*}t^{\frac{1}{4}}}\pm y^4+{\mathcal O}(t^{-\frac{1}{4}}) \end{aligned}$$

where the ± sign is determined by the sign of \(f^{(iv)}(k^*)\). Then using the Pearcey integral (2.44), the expansion (2.41) can be derived in a straightforward way from (2.46). In a similar way the expansions (2.42) and (2.43) are obtained by applying the above analysis to the expression (2.39).

In the situation \(k^{*}=1/2\) of statement (b) one uses in addition that \(t\phi _{\pm }(1/2, j/t)=tf(1/2) \pm j\pi \), \(\omega (1/2)=-\sum _{s=1}^m\tau _s (1-\cos (\pi s))=-2\sum _{s \,odd}\tau _s\), see (2.21), and consequently \(e^{\pm i\theta (1/2)} = -\) sgn\((\sum _{s \,odd}\tau _s)\). The leading order contribution of the stationary phase point \(k^{*}=1/2\) to the integral representation of, say, \(S_{12}\) in (2.39) is then given by

$$\begin{aligned} -\text{ sgn }\,\left( \!\sum \limits _{s \,odd}\tau _s\right) \frac{(-1)^j}{2 \beta \pi \lambda ^{*}t^{\frac{1}{4}}}\mathfrak {I}\left( e^{i t f(\frac{1}{2})}\left( \int _{-\infty }^{0}e^{i(\pm y^4-wy)}dy + \int _{-\infty }^{0}e^{i(\pm y^4+wy)}dy\right) \right) \end{aligned}$$

with \(w=\frac{j}{\lambda ^{*}t^{\frac{1}{4}}}\). 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'''(1/2)=0\) by symmetry (2.23) we only need to solve

$$\begin{aligned} f''\left( \frac{1}{2}\right) =0\quad \text{ which } \text{ is } \text{ equivalent } \text{ to } \quad \sum _{s=1}^ms^2(-1)^{s+1}\kappa _s=0\,. \end{aligned}$$
(2.47)

The solution of the above equation is

$$\begin{aligned} \kappa _m=\frac{(-1)^m}{m^2}\sum \limits _{s=1}^{m-1}s^2(-1)^{s+1}\kappa _s. \end{aligned}$$
(2.48)

It is clear from the above relation that for m even, choosing \(\kappa _1\) sufficiently big one has \(\kappa _m>0\) while for m odd, it is sufficient to choose \(\kappa _{s+1}>\frac{s^2}{(s+1)^2}\kappa _s > 0\), s odd and \(1\le s \le m-2\). Note that in the situation of (2.48) \(f^{(iv)}(\frac{1}{2})\ne 0\) holds iff \(\sum _{s=1}^m\kappa _s s^4(-1)^{s+1}\ne 0\). This condition simply removes an \((m-2)\)-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 \((m-1)\)-parameter family of potentials such that the correlation functions decay as in (2.45).

Finally, we show for \(m \ge 4\) our claim about the solution set of (2.40). The case \(m=3\) is treated in Example 2.9. Our strategy is to first show that there exists a \({\varvec{\kappa }}^{*}\) that satisfies \(f''(1/4, {\varvec{\kappa }}^{*}) = 0\), \(f'''(1/4, {\varvec{\kappa }}^{*}) = 0\), \(f'(1/4, {\varvec{\kappa }}^{*}) > 0\), and \(f^{(iv)}(1/4, {\varvec{\kappa }}^{*}) \ne 0\). We then invoke the Implicit Function Theorem to show the existence of the \((m-1)\)-dimensional solution manifold, where the stationary phase point \(k^{*}\sim 1/4\) may and will depend on the parameters. The conditions \(f''(\frac{1}{4}, {\varvec{\kappa }})=0\) and \(f'''(\frac{1}{4}, {\varvec{\kappa }})=0\) imply

$$\begin{aligned}&f'''\left( \frac{1}{4}\right) =0\rightarrow \sum \limits _{s \,odd}(-1)^{\frac{s-1}{2}}s^3\kappa _s=0, \end{aligned}$$
(2.49)
$$\begin{aligned}&f''\left( \frac{1}{4}\right) =0\rightarrow \left( 2\sum \limits _{s \,odd}\kappa _s + 2\sum _{s \,even}\kappa _s(1-(-1)^{\frac{s}{2}})\right) \sum _{s \,even}s^2\kappa _s(-1)^{\frac{s}{2}} \nonumber \\&\quad -\left( \sum _{s \,odd}s\kappa _s(-1)^{\frac{s-1}{2}}\right) ^2=0\,. \end{aligned}$$
(2.50)

One needs to treat the case m odd and even separately. Here we consider only the case m even. The odd case can be treated in a similar way. Equation (2.49) gives

$$\begin{aligned} \kappa _{m-1}=\dfrac{(-1)^{\frac{m}{2}}}{(m-1)^3}\sum \limits _{s \,odd,s=1}^{m-3}(-1)^{\frac{s-1}{2}}s^3\kappa _s. \end{aligned}$$

If \(m=2\ell \) with \(\ell \) even, a positive solution \(\kappa _{m-1}\) exists, provided that \(\kappa _1\) is sufficiently big. If \(m=2\ell \) with \(\ell \) odd then one needs to require \(0<\kappa _{s}<\frac{ (s+2)^3}{s^3}\kappa _{s+2}\) for \(s=1,5,9,\ldots ,m-5\).

The equation (2.50) is a linear equation in \(\kappa _4\) and we solve it for \(\kappa _4\) obtaining

$$\begin{aligned} \kappa _4= & {} \frac{1}{32}\dfrac{\left( \sum \limits _{s \,odd,s=1}^{m-3}\kappa _s(-1)^{\frac{s-1}{2}}s(1-\frac{s^2}{(m-1)^2})\right) ^2}{\sum \limits _{s \,odd,s=1}^{m-3}\kappa _s(1+\frac{s^3(-1)^{\frac{m+s-1}{2}}}{(m-1)^3}) + \sum \limits _{s \,even,s=2}^{m}\kappa _s(1-(-1)^{\frac{s}{2}})} \\&\quad +\frac{1}{16}\sum _{s \,even, s\ne 4, s=2}^{m}s^2\kappa _s(-1)^{\frac{s-2}{2}}\,. \end{aligned}$$

We observe that the first term in the above expression is always positive, while the second term is positive if we require that \(\kappa _{s}>\frac{ (s+2)^2}{s^2}\kappa _{s+2}>0\) for \(s=6,10, 14,\ldots ,m-2\). The remaining two conditions \(f'(1/4) > 0\) and \(f^{(iv)}(1/4)\ne 0\) are easy to satisfy: the sign of \(f'(1/4)\) agrees with the sign of \(\sum _{s \,odd}s\kappa _s(-1)^{\frac{s-1}{2}}\) and can be made positive by choosing \(\kappa _1\) sufficiently large. In the situation where (2.49) and (2.50) hold the fourth derivative \(f^{(iv)}(1/4)\) does not vanish iff \(\sum _{s \,even}s^4\kappa _s(-1)^{\frac{s}{2}} \ne 0\). This can be achieved by adjusting, for example, the value of \(\kappa _2\). We have now shown that there exists \({\varvec{\kappa }}^{*} \in {\mathbb R}_+^m\) such that the first four derivatives of f have all desired properties at \(k=1/4\). In order to obtain the \((m-1)\)-dimensional solution manifold in parameter space, we apply the Implicit Function Theorem to \(F(k, {\varvec{\kappa }}) := (f''(k, {\varvec{\kappa }}), f'''(k, {\varvec{\kappa }}))\). By a straight forward computation on sees that

$$\begin{aligned} \det \left[ \frac{\partial F}{\partial (k, \kappa _4)} (1/4, {\varvec{\kappa }}^{*}) \right] = - f^{(iv)}(1/4, {\varvec{\kappa }}^{*}) \, \frac{\partial f''}{\partial \kappa _4}(1/4, {\varvec{\kappa }}^{*}) \, \ne \, 0\,. \end{aligned}$$

We can therefore solve \(F(k, {\varvec{\kappa }})=0\) near \((1/4, {\varvec{\kappa }}^{*})\) by choosing \((k, \kappa _4)\) as functions of the remaining parameters \(\kappa _j\) with \(j\ne 4\). \(\square \)

Example 2.8

m even. Choosing \(\kappa _s=\frac{1}{s^2}\) for \(s=1,\ldots , m\) one has that conditions (2.47) are satisfied and \(f^{(iv)}\left( \frac{1}{2}\right) <0\).

For \(\kappa _s=\frac{1}{s^{\alpha }}\), \(s=1,\ldots , m-1\), \(2<\alpha <3\), and \(\kappa _m\) given by (2.48), there is \(\alpha =\alpha (m)\) such that \(\kappa _{m}<\kappa _{m-1}\).

m odd. Choosing \(\kappa _s=\frac{1}{s}\), for \(s=1,\ldots m-1\), one has from (2.48) \(\kappa _m=\frac{ m-1}{2m^2}<\kappa _{m-1}\) and \(f^{(iv)}(\frac{1}{2})>0\).

In all these examples the correlation functions \(S_{\alpha \alpha '}(j,t)\), \(\alpha ,\alpha '=1,2\) decrease as \(t^{-\frac{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(k) the point \(k^*=\frac{1}{3}\) thus obtaining the equations

$$\begin{aligned} \kappa _2=\frac{1}{8}\kappa _1,\quad \kappa _3=\frac{7}{72}\kappa _1\,. \end{aligned}$$

The speed of the peak is \(v^*=\dfrac{\sqrt{2\kappa _1}}{4}\) and \(f^{(iv)}(\frac{1}{3})=-\frac{68\sqrt{6}}{6}\pi ^4\sqrt{\kappa _1}\).

The correlation functions \(S_{\alpha \alpha '}(j,t)\), \(\alpha ,\alpha '=1,2\) decrease as \(t^{-\frac{1}{4}}\) and \(S_{33}(j,t)\) decreases like \(t^{-\frac{1}{2}}\) as \(t\rightarrow \infty \) and \(j\sim v^*t\), see Fig. 3. Note that one may obtain a 2-parameter family of solutions of (2.40) by picking, for example, the particular solution related to \(\kappa _1 = 1\) and by showing that the system of equations \((f'',f''')(k,{\varvec{\kappa }})=0\) can be solved near (1/3, 1, 1/8, 7/72) by choosing k and \(\kappa _3\) as functions of \(\kappa _1\) and \(\kappa _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: \({\widehat{H}}_j=\frac{1}{2}\left( |{\widehat{p}}_j| ^2 + |\omega _j|^2 |{\widehat{q}}_j|^2 \right) \), \(j=0,\ldots \frac{N-1}{2}\). However, when written in the physical variables \(\mathbf{p}\) and \( \mathbf{q}\), the quantities

$$\begin{aligned} {\widehat{H}}_j=\frac{1}{2}\sum _{k,l=0}^{N-1} {\mathcal F}_{j,k} \overline{{\mathcal F}_{j,l} }( p_k p_l+|\omega _j|^2 q_k q_l) \end{aligned}$$

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 \(\{\mathbf{e}_k\}_{k=0}^{N-1}\) the canonical basis in \({\mathbb R}^N\).

Theorem 3.1

Let us consider the Hamiltonian

$$\begin{aligned} H(\mathbf{p},\mathbf{q}) = \frac{1}{2} \mathbf{p}^\intercal \mathbf{p}+ \frac{1}{2}\mathbf{q}^\intercal A \mathbf{q}\,, \end{aligned}$$
(3.1)

with the symmetric circulant matrix A as in (1.2), (1.3). Define the matrices \(\{G_k\}_{k=1}^{M}\) to be the symmetric circulant matrix generated by the vector \(\frac{1}{2}(\mathbf{e}_k + \mathbf{e}_{N - k})\) and \(\{S_k\}_{k=1}^{M}\) to be the antisymmetric circulant matrix generated by the vector \(\frac{1}{2}(\mathbf{e}_k - \mathbf{e}_{N - k})\). Then the family of Hamiltonians defined as

$$\begin{aligned} H_k(\mathbf{p},\mathbf{q})&= \frac{1}{2} \mathbf{p}^\intercal G_k \mathbf{p}+ \frac{1}{2}\mathbf{q}^\intercal T^\intercal G_kT \mathbf{q}=\frac{1}{2}\sum _{j=0}^{N-1}[p_jp_{j+k}+r_jr_{j+k}]\, ,\end{aligned}$$
(3.2)
$$\begin{aligned} H _{k+ \frac{N-1}{2}}(\mathbf{p},\mathbf{q})&=\mathbf{p}^\intercal T^\intercal S_k T \mathbf{q}=\frac{1}{2}\sum _{j=0}^{N-1}\left[ \left( \sum _{\ell =0}^m\tau _\ell p_{j+\ell }\right) (r_{j+k}-r_{j-k})\right] \,,\, \quad \end{aligned}$$
(3.3)

for \(k=1,\ldots ,\frac{N-1}{2}\) together with \(H_0:=H\) forms a complete family \((H_j)_{0\le j \le 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^\intercal =G_k\), \(0\le k\le (N-1)/2\), the Poisson bracket \(\{F, G\} = \langle \nabla _{\mathbf{q}} F, \nabla _{\mathbf{p}} G\rangle - \langle \nabla _{\mathbf{q}} G, \nabla _{\mathbf{p}} F\rangle \) may be evaluated in the form

$$\begin{aligned} \begin{array}{lll} \{H_k, H_\ell \} &{}= \mathbf{q}^\intercal \big ( T^\intercal G_{k} T G_{\ell } - T^\intercal G_{\ell } T G_{k} \big ) \mathbf{p}\,,\qquad &{}\text{ for } 0 \le k, \ell \le \frac{N-1}{2}, \\ \{H_k, H_\ell \} &{}= \mathbf{p}^\intercal \big ( T^\intercal S_k T T^\intercal S_\ell T - T^\intercal S_\ell T T^\intercal S_k T \big ) \mathbf{q}\,,\qquad &{}\text{ for } \frac{N+1}{2} \le k, \ell \le N-1, \\ \{H_k, H_\ell \} &{}= \mathbf{q}^\intercal T^\intercal G_{k} T T^\intercal S_\ell T \mathbf{q}- \mathbf{p}^\intercal T^\intercal S_\ell T G_k \mathbf{p}\,,\qquad &{}\text{ for } 0 \le k \le \frac{N-1}{2}, \frac{N+1}{2} \le \ell \le N-1. \end{array} \end{aligned}$$

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_\ell \) is skew symmetric: \(S^\intercal _\ell = -S_\ell \). \(\square \)

Now we introduce the local densities corresponding to the just defined integrals of motion

$$\begin{aligned} e_j^{(k)}= {\left\{ \begin{array}{ll} &{} \frac{1}{2}\left( p_jp_{j+k} + r_jr_{j+k}\right) \, , \text{ for } k=1, \ldots ,\frac{N-1}{2} \\ &{}\left( \sum _{l=0}^m \tau _l p_{j+l}\right) \left( r_{j+k} - r_{j-k}\right) ,\;\; \text{ for } k=\frac{N+1}{2}, \ldots ,N\,. \end{array}\right. } \end{aligned}$$

together with their correlation functions

$$\begin{aligned} S_{(k+3,n+3)}^{(N)}(j,t) := \left\langle e^{(k)}_j(t)e^{(n)}_0(0) \right\rangle - \left\langle e^{(k)}_j(t)\right\rangle \left\langle e^{(n)}_0(0)\right\rangle \,. \end{aligned}$$
(3.4)

and limits

$$\begin{aligned} S_{k,n}(j,t)=\lim _{N\rightarrow \infty }S^{(N)}_{k,n}(j,t). \end{aligned}$$
(3.5)

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\rightarrow \infty \).

3.2 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\rightarrow \infty \) in the case of nearest neighbour interactions.

For introducing the currents we recall that \(\mathbf{r}=T\mathbf{q}\) with T as in (1.7). Then one has

$$\begin{aligned} \begin{aligned}&\dot{r}_j=\sum _{\ell }T_{j\,\ell } p_\ell =\sum _{\ell =1}^{{m}}\tau _\ell ( p_{j+\ell }-p_j),\quad r_{j+N}=r_j\\&\dot{p}_j=-\sum _{\ell }T_{\ell \,j} r_\ell =\sum _{\ell =1}^{{m}}\tau _\ell (r_j- r_{j-\ell }),\quad p_{j+N}=p_j,\;\;j=0,\ldots ,N-1. \end{aligned} \end{aligned}$$
(3.6)

To write the above equation in the form of a discrete conservation law we introduce the local currents

$$\begin{aligned}&{\mathcal J}_j^{(r)} :=\sum _{s=0}^{m-1}p_{j+1+s}\sum _{\ell =s+1}^m\tau _\ell \, ,\quad {\mathcal J}_j^{(p)} := \sum _{s=1}^{m}r_{j+1-s}\sum _{\ell =s}^m\tau _\ell . \end{aligned}$$
(3.7)

Then the equations of motion (3.6) can be written in the form

$$\begin{aligned}&\dot{r}_j={\mathcal J}_{j}^{(r)}-{\mathcal J}_{j-1}^{(r)} \end{aligned}$$
(3.8)
$$\begin{aligned}&\dot{p}_j={\mathcal J}_{j}^{(p)}-{\mathcal J}_{j-1}^{(p)},\quad j=0,\ldots ,N-1. \end{aligned}$$
(3.9)

From the above equations it is clear that the momentum \(p_j\) and the generalized elongation \(r_j\) are locally conserved. The evolution of the energy \(e_j:=\dfrac{1}{2}p_j^2+\dfrac{1}{2}r_j^2\) at position j takes the form

$$\begin{aligned} \dot{e}_j={\mathcal J}_j^{(e)}-{\mathcal J}_{j-1}^{(e)}\,, \quad {\mathcal J}_j^{(e)}=\sum _{s=1}^{m}\tau _s\sum _{\ell =0}^{s-1}r_{j+1-s+\ell }p_{j+1+\ell }. \end{aligned}$$
(3.10)

We remark that all the currents \({\mathcal J}_{j}^{(r)}\), \({\mathcal J}_{j}^{(p)}\) and \({\mathcal J}_{j}^{(e)}\) are local quantities in the variables \(\mathbf{q}\) and \(\mathbf{p}\). Observe furthermore that the total current \(\sum _j {\mathcal J}_{j}^{(p)}\) for the momentum is a multiple of the total stretch and that the total current \(\sum _j {\mathcal J}_{j}^{(r)}\) for the stretch is a multiple of the total momentum. Therefore these two total currents are conserved quantities. However, the total current for the energy \(\sum _j {\mathcal J}_{j}^{(e)}\) is not a conserved quantity. Of course, one could easily remedy this by simply subtracting the mean total current from each \({\mathcal J}_{j}^{(e)}\) at the cost of losing the local nature of the current. It is an interesting question whether one can find more integrals of motion (besides momentum and stretch) with corresponding currents that are both local quantities in the variables \(\mathbf{q}\) and \(\mathbf{p}\) so that the currents are conserved. To the best of our knowledge this is an open question even in the simplest case of linear nearest-neighbour interactions.

We recall the notation of the introduction

$$\begin{aligned} {\varvec{u}}(j,t)=(r_j(t),p_j(t),e_j(t)), \end{aligned}$$

and we introduce the vector of currents \( {\varvec{J}}(j,t)=({\mathcal J}_{j}^{(r)}(t), {\mathcal J}_{j}^{(p)}(t), {\mathcal J}_{j}^{(e)}(t))\,.\) The equations of motion take the compact form

$$\begin{aligned} \dfrac{d}{dt}{\varvec{u}}(j,t)={\varvec{J}}(j,t)-{\varvec{J}}(j-1,t). \end{aligned}$$

We define a potential function for the above conservation law

$$\begin{aligned} \varvec{\Phi }(j,t):=\int _0^t{\varvec{J}}(j,t')dt'+\sum _{\ell =0}^j\mathbf{u}(\ell ,0). \end{aligned}$$

Then it is straightforward to verify that \(\varvec{\Phi }_t(j,t)={\varvec{J}}(j,t)\) and \(\varvec{\Phi }(j,t)-\varvec{\Phi }(j-1,t)=\mathbf{u}(j,t)\). The quantities \(\Phi _1(j,t)\) and \(\Phi _2(j,t)\) can be expressed as sums of independent centered Gaussian random variables and are therefore also Gaussian random variables with zero mean and variance \(\langle (\Phi _1(j,t))^2\rangle \) and \(\langle (\Phi _2(j,t))^2\rangle \), 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

$$\begin{aligned} \begin{aligned} \Phi _1(j, t)&=\sqrt{\kappa _{1}}\int _{0}^{t} p_{j+1}(t')dt' +\sum _{\ell =0}^{j} r_{\ell }(0) =\sqrt{\kappa _{1}}(q_{j+1}(t)-q_0(0)) \,\\ \Phi _{2}(j,t)&=\sqrt{\kappa _{1}}\int _{0}^{t}r_{j}(t')dt' + \sum _{\ell =0}^{j} p_{\ell }(0) \ . \\ \end{aligned} \end{aligned}$$
(3.11)

After some lengthy calculations one obtains:

$$\begin{aligned}&\lim _{N \rightarrow \infty } \langle \left( \Phi _1(j,t) \right) ^{2}\rangle = \frac{2\kappa _{1}}{\beta } \int _{0}^{1} |\omega (k)|^{-2} \left[ 1 - \cos {\left( | \omega (k)|t \right) } \cos {\left( 2 \pi (j+1) k \right) } \right] dk \end{aligned}$$
(3.12)
$$\begin{aligned}&\lim _{N \rightarrow \infty }\langle \left( \Phi _{2}(j,t) \right) ^{2} \rangle = \frac{2\kappa _{1}}{\beta } \int _{0}^{1} |\omega (k)|^{-2} ( 1 - \cos {(|\omega (k)|t)}) \cos {\left( 2\pi (j+1)k\right) }dk + \frac{j+1}{\beta } \,. \end{aligned}$$
(3.13)

Evaluating the r.h.s. of the above expressions in the limit \(t\rightarrow \infty \) we arrive to the following theorem.

Theorem 3.2

In the limit \(N\rightarrow \infty \) and \(t\rightarrow \infty \) the quantities \(\Phi _1(j, t)\) and \(\Phi _2(j, t)\) defined in (3.11) are Gaussian random variables that have the following large t behaviour:

$$\begin{aligned} \lim _{N \rightarrow \infty }\Phi _1(j, t) ={{\mathcal {N}}}(0,\sigma _1^2) \qquad \text{ and } \qquad \lim _{N \rightarrow \infty }\Phi _2(j, t) ={\mathcal {N}}(0,\sigma _2^2)\,. \end{aligned}$$
(3.14)

The leading order behaviour of the variances \(\sigma _1^2\) and \(\sigma _2^2\) agrees. In the physically interesting region \(\frac{|j|}{t}\le \sqrt{\kappa _1}\) it is given by

$$\begin{aligned} \sigma _1^2=\frac{ t \sqrt{\kappa _{1}}}{\beta } + {\mathcal O}\big (t^{\frac{1}{3}} \big ) = \sigma _2^2 \,. \end{aligned}$$
(3.15)

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 ( \( \sim C t^{2}\)) near \(k=0\), we consider the following Cauchy-type integral instead,

$$\begin{aligned} F_{0}(z) =\frac{1}{2 \pi ^{2} \beta } \int _{-1/2}^{1/2} \frac{ 1 - \cos {(|\omega (k)|t)}}{(k-z)^{2}} \cos {\left( 2\pi (j+1)k\right) }dk \ , \end{aligned}$$
(3.16)

which gives the leading order asymptotic behaviour of the integrals appearing in (3.13), since

$$\begin{aligned}&\frac{ 2 \kappa _{1}}{\beta } \int _{-1/2}^{1/2} |\omega (k)|^{-2} ( 1 - \cos {(|\omega (k)|t)}) \cos {\left( 2\pi (j+1)k\right) }dk - F_{0}(0) \rightarrow \ 0 \ \ \text{ as } t, j \rightarrow \infty \ .\nonumber \\ \end{aligned}$$
(3.17)

For \(\frac{|j|}{t} < (1 - \epsilon ) \sqrt{\kappa _{1}}\), \(\epsilon > 0\), the analysis of \(F_{0}(z)\) 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 \(\infty \) such that \(\frac{|j|}{t} \approx \sqrt{\kappa _{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 \rightarrow \infty \). 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 \(\Phi _1\) observe that the difference of the integrals in relations (3.12) and (3.13) is given by \(\int _{0}^{1} |\omega (k)|^{-2} \left[ 1 - \cos {\left( 2 \pi (j+1) k \right) } \right] dk\) which can also be treated by a stationary phase calculation combined with a contour deformation.

3.3 Nonlinear Regime

In this section we consider a nonlinear perturbation of the harmonic oscillators with short range interactions of the form

$$\begin{aligned} H(\mathbf{p},\mathbf{q})= & {} \sum _{j=0}^{N-1}\frac{p_j^2}{2}+\sum _{s=1}^m\kappa _s\left( \frac{1}{2}\sum _{j=0}^{N-1}(q_j-q_{j+s})^2 \right. \nonumber \\&\left. + \frac{\chi }{3}\sum _{j=0}^{N-1}(q_j-q_{j+s})^3 + \frac{\gamma }{4}\sum _{j=0}^{N-1}(q_j-q_{j+s})^4\right) \,. \end{aligned}$$
(3.18)

We consider Examples 2.8 and 2.9 with different strengths of nonlinearity namely

$$\begin{aligned}&{ m=2, \, \kappa _1 = 1, \kappa _2 = \frac{1}{4},}\;\; {\left\{ \begin{array}{ll} &{}{\chi =0.01 \hbox { and } \gamma =0.001}\\ &{}{\chi =0.1 \hbox { and } \gamma =0.01 } \end{array}\right. } \\&{ m=3, \, \kappa _1 = 1, \kappa _2 = \frac{1}{8}, \kappa _2 = \frac{7}{72}, }\;\; {\left\{ \begin{array}{ll} &{}{\chi =0.01 \hbox { and } \gamma =0.001}\\ &{}{\chi =0.1 \hbox { and } \gamma =0.01 } \end{array}\right. }\, . \end{aligned}$$

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 \(\beta ^{-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 Figs. 4, 5 while for stronger nonlinearity the fastest peaks seem to scale like \(t^{\frac{2}{3}}\) in equation (1.11), see bottom figures in Figs. 4, 5. The non generic peaks that are present in the linear cases and scale like \(t^{1/4}\) have a fast decay in the case of strong nonlinearity. However for weak nonlinearities, the central peak in the top left Fig. 4, still scales in time like \(t^{-\frac{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^{-\frac{1}{4}}\) (see Fig. 6).