  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.
