{'fraction_non_alphanumeric': 0.05366354207110766, 'fraction_numerical': 0.013838815242109934, 'mean_word_length': 4.206323996534796, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 3, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Using Hilbert transform and classical chains to simulate quantum walks

We propose a simulation strategy which uses a classical device of linearly coupled chain of springs to simulate quantum dynamics, in particular the quantum walks. Through this strategy, we obtain the quantum wave function from classical evolution. Specially, this goal is achieved with the classical momenta of the particles on the chain and their Hilbert transform, from which we construct the many-body momentum and Hilbert transformed momentum pair correlation functions yielding the real and imaginary parts of the wave function, respectively. With such wave function, we show that the classical chain's energy and heat spreading densities can be related to the wave function's modulus square. This relation indicates a concept of"phonon random walks", and thus it provides a new perspective to understand ballistic heat transport. The results here may give a definite answer to Feynman's idea of using a classical device to simulate quantum physics.


I. INTRODUCTION
In his pioneering work entitled "Simulating Physics with Computers" [1], R. Feynman posed two important questions. First: What kind of computers are we going to use to simulate physics? This led to the fundamental concept of a quantum computer [2]. The second was: Can a quantum system be simulated by a classical computer? The answer to this in the words of Feynman is: "No! Since this is called the hidden-variable problem: it is impossible to represent the results of quantum mechanics with a classical universal device". Feynman also conditions this statement (see precise details in [1]) and writes that such rather far reaching conclusion is valid provided that there is no "hocus-pocus". The aim of the present work is to theoretically build a classical device that can be used to simulate quantum dynamics. Our device is a system of springs initially prepared at thermal equilibrium. The main difficulty recognized by Feynman is, as he wrote: "... we cannot simulate ψ in the normal way", where ψ is the wave function. Indeed, quantum mechanics is built on interference and a complex valued field. Of course, the interference is a wave property and for that reason we can imagine that vibrations are useful. But how can we get the complex valued wave function? Well, as we show below, the hocus-pocus is based on the Hilbert transform [3]. The Hilbert transform takes a real function f (t) and extends it to the complex plane in such a way that it satisfies the Cauchy-Riemann equations [4]. It has been used extensively extending real signals to the complex plane. Our goal below is to identify what kind of classical object which together with its Hilbert transform gives ψ. In that sense the main focus of the paper is Feynman's second problem on how to go from a classical device to quantum mechanics. However, our work is also potentially related to his first vision. Namely, the * Electronic address: phyxiongdx@fzu.edu.cn † Electronic address: thiel@posteo.de ‡ Electronic address: Eli.Barkai@biu.ac.il field of phononics [5] suggests to use phonons as a source of computation, thus our work shows precisely how this strategy can be used to simulate quantum physics [6][7][8][9][10][11].
In the early days of quantum theory, there have been several relevant attempts to interpret quantum reality within a classical framework [12][13][14][15][16][17]. However, this is not the goal of the present paper which is focused on building a device. These attempts, Bohm's causal interpretation [12,13] and Nelson's stochastic approach [14], yield a recipe for quantum dynamics based on classical concepts. These while sharing some general themes with our construction, as we will show below, are never the less very different from what we propose (see discussions in the conclusion part). In fact, from the start, our goal is not to claim that quantum mechanics is not needed, or to replace its concepts, but rather we wish to theoretically build a machine made of springs that mimics certain aspects of quantum theory with the hope that this will be both of academic interest but also that it will advance the field of phonon computation [5]. The main motivation of this work was to construct the so-called tight-binding quantum walk using a classical device (see below), the former concept is not only relatively modern but also claimed to be useful for quantum search algorithms [18][19][20].
In what follows, we first outline the relevant theory of quantum walks and provide the corresponding wave function. Sec. III then demonstrates how to construct this wave function via the Hilbert transform. The simulation strategy to verify our theory and the scope and limitation of this strategy are also discussed. In Sec. IV, the physical meaning of the wave function in classical mechanics is analyzed, based on which, we show in Sec. V that the wave function's modulus square can intriguingly represent the classical energy and heat spreading densities. Finally, Sec. VI draws our conclusion, followed by several appendices of additional details.

II. THEORY
At first we briefly review a particular type of quantum dynamics that our classical device will be able to simulate. We consider a particle on a one-dimensional (1D) lattice described by a Hamiltonian operator H. Lattice sites m are integers and the system is infinite and translation invariant, stretching from −∞ to ∞. This implies , see Refs. [6][7][8][9][10][11]. α j are coupling constants with dimensions of an energy. The fact that the system is translation invariant is crucial in all of our analysis, which is based on Fourier transform (we could treat also the system on, say, a ring). This is also one of the main requirements from the classical device which we will soon construct. Then for a particle initially located on the origin, the amplitude ψ m (t) at time t is where q is the wave number, e imq is a free wave, and ω q is the frequency of a free wave, determined by He imq = ω q e imq . We use units whereh = 1. ω q is also called the dispersion relation. Note that i dψm(t) dt = ω q ψ m (t), which is easily verified from Eq. (1). The tight-binding Hamiltonian with nearest-neighbour (NN) jumps serves as an example. Here α 1 = 1 and all other α j = 0. Accordingly i dψm(t) dt = ψ m+1 + ψ m−1 . Inserting the eigenstate solution e imq , one gets the dispersion relation ω q = 2 cos(q), which in turn bears ψ m (t) = 1 2π π −π e i[mq−2 cos(q)t] dq = i −m J m (2t). Its modulus square is the density ρ(m, t) = |ψ m (t) is the Bessel function of the first kind [6,11], for a schematic presentation see Fig. 1(a). This density has three well-known properties: ballistic scaling tρ(m, t) ρ(m/t, t), U-like shape and oscillations (interference) (see Appendix A). The system is also called the tight-binding quantum walk [6]; the particle's amplitude is initially localized and may tunnel to adjacent lattice sites. Experimentally such quantum walks can be observed e.g. in waveguide lattices [21].

III. CONSTRUCTION OF THE WAVE FUNCTION
We next explain how to use classical particles connected with springs to construct this wave function. In theory we provide a classical analog-computer whose output is ψ m (t). This also gives the exact definition and method of measurement of ψ m (t). Before doing that, let us first divide the wave function into real and imaginary parts. By doing so, we use ω q = ω −q to obtain: and So that the density is Our device is a chain of classical particles whose labels are m, arranged on a ring of size N → ∞ with periodic boundary conditions. The particles are interacting via linear springs which give the phonon dispersion relation ω q (see examples in the Hamiltonians below). Initially, the system is in contact with a heat bath of temperature T , so the initial condition is drawn from a Boltzmann-Gibbs distribution. After preparation, we solve the Newtonian equation of motion. The particles evolve classically until some time t. We will show that the momenta of the particles and their Hilbert transform (see definition below) can give the wave function. For finite T , the classical particles' momenta are initially random, we will employ the correlation functions of these quantities to reproduce the quantum dynamics. Eventually the temperature will not play any role, in the sense that the wave function can be obtained at all temperatures. This implies that in principle our device can operate at room temperatures. Like the quantum system, the main ingredient of our classical device is that the interactions are translation invariant and linear, in that sense all particles are identical. Essentially the label of classical particles corresponds to the lattice site in the quantum problem, hence our method is focused on the quantum discrete space systems (the classical system is not space discretized). We will consider a chain of infinite size, but all of our results can be extended also for the finite size rings (see the discussion in Sec. III D and F). We treat systems with NN coupling and the Hamiltonian There are N = L+1 particles; all of them have unit mass. p m is the momentum of the mth particle; r m is its displacement from equilibrium position; ∆r m = r m+1 − r m denotes the NN stretch. Applying periodic boundary conditions, we recover a ring. More generally, the interaction potential is of the form L m=0 n 1 2 A n (r m+n −r m ) 2 , where A n is the spring constant between n-next nearest neighbors (see examples below).

B. The real part of the wave function
The two-body momentum correlation function is defined as [22] ρ p (m, t) = with χ * denoting the conjugate of χ and · the spatiotemporal average. This correlation function represents the momentum correlation function of any two particles whose separation is m with a time lag t since the system is translation invariant [of course, classical momenta are real and p * = p in (6)]. Following the method proposed by Montroll and Mazur [22], the Hamiltonian (5) is equivalent to H = 1 2 L k=0 |P k | 2 + ω 2 k |R k | 2 when applying the following normal transformation Here ω k is the kth normal mode's frequency. For example, if one considers the harmonic chain, i.e., Hamiltonian (5) with V (ξ) = ξ 2 /2, ω k = 2| sin(kπ/N )|. P k and R k are the normal coordinates; the matrix C has the form and satisfies L m=0 C m,k C * m,l = δ k,l (10) with δ representing the Kronecker symbol. So, for each normal mode, we obtain the evolution equation which determines R k and P k : and Substitute Eqs. (7-13) into (6) and use the following equipartition conditions: P k (0)P * l (0) = k B T δ k,l ; R k (0)R * l (0) = k B T δ k,l / ω 2 k , and R k (0)P * l (0) = 0 (k B is the Boltzmann constant), one finally obtains [22] Taking N (or L) → ∞ and accordingly 2πk/N → q, we have ω k → ω q . In particular we obtain the classical dispersion relation of the harmonic chain as a specific example: For any classical chain's dispersion relation we have: which is just Re[ψ m (t)] from (2). Note that the classical dispersion of a system of springs is controlled by the interactions of a bead with its neighbors. If these interactions are controllable, we may obtain rather general forms of the classical dispersion, as we demonstrate in examples below.
C. Imaginary part of ψm(t): the π/2-shifted momentum correlation function Define the following cross-correlation function where p m (t) is what we call the π/2-shifted momentum.
Since p m (t) is a linear combination of normal modes, p m (t) can be obtained by shifting the underlying normal modes and performing the transformation p m (t) = L k=0 C m,k P k (t) [see Eq. (7)]. Now use (18) and repeat the steps from ρ p (m, t) to ρ p (m, t). One gets Finally take N → ∞, to obtain which is Im[ψ m (t)] in (3).

D. Hilbert transform
The π/2-shifted momentum is mathematically equivalent to the negative Hilbert transform [3] of the momentum. In fact, the Hilbert transform's effect is a π/2 phase shift of each frequency components. For an arbitrary p m (t), one defines its Hilbert transform H [p m (t)] as (see Ref. [3]) where p.v. indicates the Cauchy principal value.
form a complex conjugate pair that defines the so-called analytic signal The last term of (23) arises because H [P k (t)] = − P k (t) = ω k R k (t) in view of (12) and (18). We now conjecture that the normalized correlation function ψ m (t) = which is just the Schrödinger equation in Fourier space that we mentioned at the beginning. dR k (t) dt and dP k (t) dt were obtained from (12) and (13). Turning back to lattice space, we get the wave function's real and imaginary parts. One might think that the relation between ψ m and Z m is just a curious mathematical coincidence. However, the quantum and the classical systems share wave-like properties. For this reason, we searched for this intriguing exact correspondence between the classical and quantum worlds. The appearance of the Schrödinger equation confirms our idea that classical systems' correlation functions can be used to obtain quantum dynamics. It is worthwhile to note that the theory works also for finite number of particles. Translation invariance of both, the classical and corresponding quantum system, ensures the applicability of Fourier analysis. Thus, the limit of N → ∞ is not a general request.

E. Construction of the wave function by simulations
In Fig. 1 we compared the wave function obtained from a quantum walk with the one obtained from a harmonic chain. To do so, we first computed the wave function's modulus square |ψ m (t)| 2 by plugging the dispersion relation (15) into Eqs. (2-4), see Fig. 1(a). Then we simulated a harmonic chain with N = 4001 particles, initially connected to a heat bath with T = 0.5 (see simulation detail in Appendix B). The same values are used below. The results are also verified for other temperatures. Employing the correlation functions, ρ p (m, t) = pi+m(t)pi(0) The result is depicted in Fig. 1(b). In practice we apply the Hilbert transform on each particle's momentum to obtain ρ p (m, t); the details are also described in Appendix B.

F. Scope and limitation
We focus on simple systems with translational invariance, both for the quantum dynamics and the classical device. The classical particles are linearly coupled. Their label m corresponds to a lattice site of the quantum system. The linear chain describes a single quantum particle on a lattice. It can be a finite ring of elements, or stretch to infinity. Furthermore, the initial state of the classical system is a thermal one. This implies uncorrelated momenta. At t = 0, the momentum correlation function is a Kronecker delta, consequently the quantum particle is initially localized. These are certainly constraints on the generality of our device and we are still far from a universal classical machine capable of simulating all aspects of quantum mechanics. Still the Hilbert transform technique is encouraging and hopefully further research will unravel more general devices. For that reason we proceed to show that the wave function can, in the classical context, be used to describe the kinetic energy, stretch, and total energy and heat correlation functions.

IV. THE PHYSICAL MEANING OF THE WAVE FUNCTION IN CLASSICAL MECHANICS
So far, we have shown how to construct the wave function from a linear chain's correlation functions. Now we would like to point out the physical meaning of the wave function in the classical domain. It is a physically significant observable describing various equilibrium correlation functions beyond the momentum and its Hilbert transformed correlation functions (see Appendix C). First, one can rigorously prove that the square of the wave function's real part {Re[ψ m (t)]} 2 is the kinetic energy correlation function (see Appendix C1). Second, the non-normalized stretch correlation function C ∆r (m, t), defined by ∆r which gives Furthermore, we will demonstrate that the potential energy and the total energy correlation functions can also be related to the wave function (see Appendix C4). Based on all of these facts, below we will provide evidences that the density ρ(m, t) = |ψ m (t)| 2 describes both the normalized (total) energy and heat correlation functions, in the long time limit [23]. Thus, we see, that the proposed wave function contains rich physical information on the classical chain.

V. ENERGY AND HEAT SPREADING DENSITIES
A. U shaped density for harmonic chain As demonstrated in Fig. 1, given a phonon dispersion relation, one can simulate a quantum walk by observing The rescaled densities for the harmonic chain (t = 600), obtained from plugging Eq. (15) into Eqs. (2)(3)(4), are compared with the prediction from the Arcsine distribution and simulations.
the wave function's modulus square [see Eq. (4)]. Here we compare the prediction of Eqs. (2)(3)(4) to the simulations of energy and heat spreading densities for a harmonic chain in the long-time limit. The energy and heat spreading densities are usually obtained from the following correlation functions, i.e., ρ E (m, t) = ∆Ei+m(t)∆Ei (0) ∆Ei (0)∆Ei (0) and and Q i (t) are the energy and heat densities at location i and time t. ∆E i (t) and ∆Q i (t) are their fluctuations, respectively (see Refs. [24][25][26][27][28][29] and Appendix B for detailed definitions).
Intriguingly, we find that both densities nicely match and converge to an U-shape (see Fig. 2). This U-shaped solution is Lévy's well-known Arcsine law that also describes the occupation times of an 1D Brownian particle in a half space (see [30][31][32]).

B. Quasi particle velocity approach
To understand this U-shape, we here propose a method we call "quasi particle velocity approach" (QPVA), which perfectly gives rise to the Arcsine law in the harmonic chain.
Our main idea is that in a linear chain, phonons can be understood as quasi particles. They start on the origin and will spread out ballistically with m = vt. In this case, given the probability distribution function (PDF) of the particle's velocity, h(v), the density is given by To obtain the velocity PDF, one first takes the Fourier which is just the characteristic function. h(v) is obtained by inverse Fourier transform of this characteristic function. Thus, we need to find e iµv . In our study, we consider highly localized initial conditions. This corresponds to an uniform distribution of wave vectors, q. A free wave's velocity v q can be identified with the group  velocity (dω q /dq). Hence in the harmonic chain with dispersion relation ω q = 2 |sin (q/2)| we have then hence Finally, substituting (32) into (28), we finally get the rescaled Arcsine distribution which has the predicted U-shape and the correct ballistic scaling in the long time limit.

Models
Dispersion relation II ω ± q = k1 + k2 ± k 2 1 + k 2 2 + 2k1k2 cos(2q) III ωq = 4 sin 2 (q/2) + 1 IV ωq = 2 sin 2 (q/2) + γ sin 2 (q) where ω − q (ω + q ) denotes the frequency of acoustic (optical) phonons; k1 and k2 are the strength of the adjacent couplings; γ represents the comparative strength of the NNN to the NN couplings. We use k1 = 1/3, k2 = 2/3 and set γ = 0.25 here. Figure 3(b,g,l) show the prediction for ρ E (m, t) and ρ Q (m, t) from Eqs. (2-4) and Fig. 3(c,h,m,d,i,n) depicts the simulated values. All features of the simulation are nicely captured by our prediction. This clearly demonstrates the generality of our approach. We note that extending our theory to Model III (IV) is straightforward, one just needs to integrate Eqs. (2)(3)(4) with the ω q given in Table I. However, for Model II, one should consider contributions from both acous- The solution then is ρ(m, t) = |ψ − m (t) + ψ + m (t)| 2 (see detailed analysis in Appendix D). The excellent agreement between simulation and theory indicates that our theory also works well for systems with two branches of phonons. This unusual fact may stimulate the conception of new phononics devices [37], since one may be able to identify independently the contributions of acoustic and optical phonons.
From Fig. 3 one sees that the ballistic transport exhibits non-universal features (unlike the super-diffusion investigated for example in [38]). Motivated by this, we use the proposed QPVA to find the long time asymptotic behaviour of the packet, and want to understand how ex-  In the inset of (d) the y-axis is logarithmic.
actly it depends on ω q . Similar to the harmonic chain, one may view the problem as a packet of quasi particles initially localized in space. Each of the quasi particles travels with speed v q = dωq dq . All values of −π < q < π are uniformly distributed. We find For the harmonic chain with the dispersion relation (15), this yields the mentioned Arcsine law. For other ω q , Fig. 3 shows excellent agreement between simulations and our formula. We note that this picture only yields the long time limit of ρ(m, t). It does not feature the fine oscillations so typical of the ballistic dynamics (These oscillations are presented in Fig. 1, and they are important in the short time limit).

D. Discussion on non-linear integrable Toda system
As a last example, we consider the celebrated Toda chain with Hamiltonian (5) and V (ξ) = ae −bξ /b + aξ + c, which also bears ballistic transport [39]. The aim here is to discuss whether our ideas could apply to a non-linear integrable system. The Toda chain's dispersion relation has the explicit form (see Ref. [39]). Inserting this relation into Eqs. (2)(3)(4) gives predictions for the packet shape, which can be compared with simulations. Here Sn is the Jacobian elliptic function with modulus κ. κ is some constant determined by a, b and T that describes the non-linearity, 0 ≤ κ < 1.
is the complete elliptic integral of the first (second) kind. In contrast to the linear systems, ω q is no longer a single line, it rather depends on κ, and hence on T . For κ close to zero, we get harmonic behaviour. As κ → 1, v q (the group velocity defined by dωq dq ) vanishes for large q and increases for small q [see Fig. 4(a)]. We show below how these features affect the density's shape.
Predictions for the densities are shown in Fig. 4(b)-(d) for several values of κ. For small κ, the Toda chain behaves like a linear system: an U-shape can be observed [ Fig. 4(b)]. As κ increases, the density's central parts become humped together and the two side peaks emerge at |m/t| > 1 [see Fig. 4

(c)-(d) and inset].
Numerical results for ρ Q (m, t) are given for comparison. The simulations are performed with V (ξ) = e −ξ + ξ − 1, i.e. a = b = −c = 1 (see details in Appendix B). Since we lack information about the T -dependence of the dispersion relation, the comparison is just indirect and we are unable to provide predictions for a given T . From Fig. 4(f)-(h) we can indeed verify the two trends of ρ(m, t): for small T , the central parts of ρ Q (m, t) are very similar to the U-shape and slight side peaks at |m/t| > 1 can already be identified [ Fig. 4(f)]. As T increases, more and more front parts emerge and also a hump in the central part appears [ Fig. 4(g)-(h)]. Evidently the trends in simulations and predictions coincide, although both do not match precisely.
We also examined the front parts located at |m/t| > 1 and find that they can be fitted quite well with a Gaussian distribution [ Fig. 4(e)]. As T increases, so do the mean and variance of the Gaussian. This is in good agreement with the velocity fluctuations conjecture that was suggested by the Lévy walks approach for predicting the non-linear non-integrable Fermi-Pasta-Ulam-β chain's density [38].

VI. CONCLUSION
In summary, we have demonstrated how to use classical chains of springs to simulate quantum dynamics (in particular, the quantum walk). To do this, we have suggested to use the momenta's and their Hilbert transform's pair correlation function to construct a quantum like wave function. Such a strategy successfully solves the challenges addressed by Feynman on this topic. Therefore it provokes the general idea of making a classical machine to reproduce quantum dynamics. We leave to future work if the Hilbert transform technique can be used to model other aspects of quantum mechanics, e.g., spin, magnetic field, and many-body systems. We will show in a future publication that our device can work also for non-translation invariant systems and can thus model transport in e.g. disordered systems. In that case the Fourier analysis used all along this text does not work, and other tools are needed to solve the problem.
Previous methods [12][13][14] used the solution of the Schrödinger equation ψ = Re iS and then constructed a classical ensemble that yields back ψ. Their paths are generated using a deterministic lawẋ = ∂ x S [12,13] or the Langevin equation dx = (ν∂ x R/R + ∂ x S)dt + ν 1/2 dη [14], where x in their case is the particle's displacement, ν is proportional toh and η is Gaussian white noise. Also we use ensemble of particles to construct the wave packet, and we combine both stochastic (our initial conditions are drawn from the Boltzmann-Gibbs distribution) and deterministic (solution of Newtonian equations) tools. But here the resemblance ends. In our approach we used the Hilbert transform which is the most natural way to extend a signal to the complex plane and then we construct the classical correlation functions which yield the real and imaginary parts of the wave function. Furthermore, in our examples, the transformation is used on the individual trajectory level and only for momentum, so our approach is based on extending the classical trajectory to the complex plane. Therefore, the wave function in our case has a classical interpretation in terms of classical observables without invoking fictitious forces which have no classical analogue.
With above strategy, we have found that this proposed wave function's modulus square corresponds to the classical energy and heat spreading densities. Such densities have been found to exhibit non-universal shapes (dependent on the phonon dispersion relation), showing quantitative agreements with the simulation results of ballistic heat transport in many integrable systems. We have also proposed the quasi-particle velocity approach to understand the long time asymptotic behaviour of these ballistic wave packets. This similarity to quantum walks together with the picture of the quasi particle velocity approach provides a new perspective to understand ballistic heat transport. An extension of these ideas for general nonlinear, nonintegrable systems will provide a "hydrodynamic foundation" of these models [31,43]. As mentioned, the density of the nearest-neighbor (NN) tight-binding quantum system has three key properties: ballistic scaling, U-shape and oscillations. All of these are nicely demonstrated in Fig. 5, where the rescaled density tρ(m, t) is plotted against m/t. The oscillations are better visible in the inset.

Appendix B: Simulation detail
We mainly focus on the following four correlation functions. (i) The momentum correlation function and (ii) the cross-correlation function between momentum and its Hilbert transform In relation to heat transport, we use (iii) the correlation function for energy fluctuations [24][25][26][27][28][29] (iv) and the correlation function for heat energy fluctuations Here m = j − i; · represents the spatio-temporal average; ∆χ ≡ χ − χ is the corresponding quantity's fluctuations. For (i)-(iii) the labels i and j correspond to the labels of particles. The energy E i is defined by the sum of kinetic energy p 2 i /2 and potential energy V (which may also depend on the position of other particles).
The heat density Q i is defined for a finite volume (bin). Its expression can be derived from basic thermodynamics. For details we refer to textbooks [24,25] and to other publications [26][27][28][29]. The indices for the heat fluctuation density correspond to bin labels rather than particle labels. In each bin, we calculate the number of particles in the bin M i , the energy in the bin E i and the pressure F i exerted on the bin. Finally the heat is obtained from Since the system is one dimensional the pressure is equal to the force and can be calculated from the gradient of the potential.
For the Toda chain, its general potential is V (ξ) = a b e −bξ + aξ + c, from which one can derive a modulus dependent dispersion relation [39]. For simulations at finite temperature, we employ the simple form V (ξ) = e −ξ + ξ − 1, i.e., we simply set a = b = −c = 1. The modulus κ shown in the Toda chain's dispersion relation [39] is some constant determined by the parameters a, b and T describing the non-linearity.
For all the simulations, we set both the equilibrium distance between the particles as well as the lattice constant to unity. So the number of particles N is equal to the system size. All systems except the Toda chain have symmetric potentials. Therefore the average pressure F calculated from simulations is always zero. For the Toda chain with an asymmetric potential, the average pressure depends on temperature. Our simulations indicate F ≈ 0.05, F ≈ 0.48 and F ≈ 1.34 for T = 0.05, T = 0.5 and T = 1.5, respectively.
We consider a chain of size N = 4001 with periodic boundary conditions. To allow heat and energy fluctuations to actually spread out, we compute the correlation function up to a lag time of t = 600. The heat correlation function was calculated from discretised chain with 2000 bins.
We use the stochastic Langevin heat baths [40,41] to thermalize the system and to prepare a canonical equilibrium state with a given temperature. We employ the Runge-Kutta algorithm of 7-th to 8-th order with a time step of 0.05 to evolve the system. Each canonical equilibrium system is prepared by evolving the system for a long enough time (> 10 7 time units) from properly assigned random initial states. Then the heat bath is switched off and we switch to the microcanonical ensemble with fixed energy. Finally the system is evolved in isolation and we obtain the correlation information. We used ensembles of circa 8 × 10 9 data points to compute the correlation functions. We also note that to get the final results of the correlation functions, one should consider a correction term as suggested in [27]. For more details on the implementation and techniques one can refer to [42].
To simulate the correlation functions of ρ p (m, t) and ρ p (m, t), we first record the time series of p m (t). Then p m (t) is obtained by using the numerical Hilbert transform [3]. The pairs are then used to calculate the correlation functions. As an example, Fig. 10 plots the results 6: pm(t) and pm(t) versus t for a harmonic chain, here we use the particle with the label of m = 1000 for example. of p m (t) and p m (t) for the harmonic chain; Fig. 7 depicts the correlation functions of ρ p (m, t) and ρ p (m, t) from simulations and compares them with the wave function's real and imaginary parts. The kinetic energy of the mth particle is E k m = p 2 m /2, whence its correlation function ρ E k (m, t) is defined by (C1) Since the system is initially in contact with heat baths, the initial momenta and the initial positions are jointly normal distributed. The equipartition conditions stated in the main text give the variance/covariance of these Gaussian random variables: P k (0)P * k (0) = δ k,k k B T , R k (0)R * k (0) = δ k,k k B T / ω 2 k , and P k (0)R * k (0) = 0. In particular Wick's theorem applies and the fourth moments can be replaced with the second moments, e.g. P 4 k (0) = 3 P 2 k (0) 2 = 3(k B T ) 2 . This allows us to compute the denominator, which is the normalization condition. First note that it can be simplified to Hence we have (by also using translational invariance in time and space): Now we need to find the numerator. We use Eqs. (7)(8)(9)(10), and employ Einstein's notation: the summation is carried out over all k's and l's: Using Eq. (9), we can substitute C 0,l C * 0,l with 1/N . Furthermore, since R k (0) and P l (0) are uncorrelated, a lot of the occurring terms actually vanish. We are left with two big sums: We apply Wick's theorem on the remaining expectations: P k (0)P k (0)P l (0)P l (0) = (k B T ) 2 (δ l,l δ k,k + δ k,l δ k ,l + δ k,l δ k ,l ). The other expectation is easier, because R k (0) and P l (0) are not correlated: R k (0)R k (0)P l (0)P l (0) = δ k,k δ l,l (k B T ) 2 / ω 2 k . We obtain: The sum over l and l cancels some of the 1/N factors. The terms with δ k,k lead to cos 2 ( ω k t) + sin 2 ( ω k t) = 1, which bears the first summand of the result. We obtain where we used C m,k = 1 √ N exp 2πi mk N . Finally, taking N = L + 1 → ∞, we find that (C6) is actually (C7) the square of the wave function's real part {Re[ψ m (t)]} 2 . To verify our proof, we compared the real part with the kinetic energy's correlation function in Fig. 8 for the harmonic chain.
An alternative way to prove the equivalence between {Re[ψ m (t)]} 2 and ρ E k (m, t) relies only on the Gaussian property of the momenta. Considering two Gaussian random variables X, Y with zero mean ( X = 0; Y = 0) and variance ( X 2 = σ 2 X ; Y 2 = σ 2 Y ). We write their joint distribution as which then helps us to calculate |p m (t)| 2 |p 0 (0)| 2 : Obviously, |p m (t)| 2 = |p 0 (0)| 2 = k B T , and = pm(t)p0(0) k B T = ρ p (m, t) is just the momentum correlation. So, we obtain

Stretch correlation function
Define the NN stretch as s m = ∆r m = r m+1 − r m , the stretch correlation function is Repeating the steps of the last section, we get: The non-normalized numerator is Now it is interesting to find that C s (m, t) is related to the wave function via (C15) Take the harmonic chain as an example. Its dispersion relation is ω q = 2 sin q 2 = 2 − 2cos(q). Inserting this expression one gets ρ s (m, t) = 1 2π π −π cos(qm) cos(ω q t)dq which is verified by simulations in Fig. 9(a) and (c).

Stretch-momentum cross-correlation function
The stretch-momentum cross-correlation function is defined as Following the similar steps of the above, one finds which is related to the wave function by d dt Inserting the harmonic chain's dispersion relation ω q = 2 sin q 2 = 2 − 2cos(q), we obtain Similarly, it is easy to find that the momentum-stretch cross-correlation function is In Fig. 9(b) and (d), we compare the prediction of C sp (m, t) with simulation. As expected, they agree well with each other.

Potential energy and total energy correlation functions
Following reference [23], one can define the potential energy for harmonic chain as Under this definition, the potential energy correlation function is [s 0 (0)] 4 − [s 0 (0)] 2 2 = ρ s 2 (m, t). (C23) Since r m (t) and r 0 (0) are the Gaussian random variables, so are s m (t) and s 0 (0). Given ρ s (m, t), ρ s 2 (m, t) can be derived by using the alternative argument of appendix C1. From above, we already know |s 0 (0)| 2 = which is related to the wave function. Next, we deal with the total energy correlation function. The total energy is defined by Under this definition, its non-normalized correlation function is The terms of k B T and k B T /2 are because the equipartition conditions tell us E 0 (0) = E m (t) = k B T , E k 0 (0) = E k m (t) = k B T /2, and E p 0 (0) = E p m (t) = k B T /2. Now C E (m, t) can be divided into the following four terms: with [C sp (m, t)] 2 . Therefore, the non-normalized total energy correlation function is [23]  which then is related to the wave function. Figure 10 presents both the predictions and simulations for ρ Ep (m, t) and C E (m, t). As can be seen, they agree well with each other.  , respectively, then set ρ T (m, t) = ρ − (m, t) + ρ + (m, t). In Fig. 13 we compare the result of ρ T (m, t) with ρ(m, t) and simulations. Fortunately, we find that in this particular model, the shape of ρ T (m, t) [ Fig. 13 . We therefore argue that one would be able to separate the contributions of acoustic and optical phonons. This may stimulate possible applications for the design of phononics devices [37].