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-13 for position and momentum correlations and as t-23 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-14 for position and momentum correlators and with rate t-12 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.
