Integrable hydrodynamics of Toda chain: case of small systems

Passing from a microscopic discrete lattice system with many degrees of freedom to a mesoscopic continuum system described by a few coarse-grained equations is challenging. The common folklore is to take the thermodynamic limit so that the physics of the discrete lattice describes the continuum results. The analytical procedure to do so relies on defining a small length scale (typically the lattice spacing) to coarse grain the microscopic evolution equations. Moving from the microscopic scale to the mesoscopic scale then requires careful approximations. In this work, we numerically test the coarsening in a Toda chain, which is an interacting integrable system, i.e., a system with a macroscopic number of conserved charges. Specifically, we study the spreading of fluctuations by computing the spatio-temporal thermal correlations with three different methods: (a) using microscopic molecular dynamics simulation with a large number of particles; (b) solving the generalized hydrodynamics equation; (c) solving the linear Euler scale equations for each conserved quantities. Surprisingly, the results for the small systems (c) match the thermodynamic results in (a) and (b) for macroscopic systems. This reiterates the importance and validity of integrable hydrodynamics in describing experiments in the laboratory, where we typically have microscopic systems.

The hydrodynamic theory is a cornerstone for understanding coarse-grained phenomena in many-body systems.The local evolution of slowly changing fields of conservation laws manifests the large-scale space-time dynamics of the system.This old theory has found new interest to predict transport in low-dimensional systems relevant to the current technologies.In nonintegrable systems, with few relevant conservation laws, the non-linear fluctuating hydrodynamic theory helps to understand diffusive and non-diffusive transport [1].In integrable systems with many conservation laws, transport is historically pictured as an underlying quasiparticle undergoing deterministic scattering.Recent experiments have made novel classical and quantum dynamics in these systems [2].However, such a formalism can only sometimes allow for the analytical tractability of physically relevant quantities measured in the lab or computer experiments.Integrable models such as harmonic chain and hard point particle gas (or models reducible to them) have been studied within the standard kinetic theory formalism to connect this quasiparticle picture and to understand coarse-grained transport properties [3].The assumed simplicity of a many-body integrable system only sometimes implies the analytical tractability of physically relevant quantities.The complexity arises for interacting systems such as the XXZ chain or the Toda chain [4].The solvability in these systems are addressed using Lax pairs in classical [27] and using the Bethe ansatz [5] in quantum setup but characterizing ,dynamical properties remain challenging using these formalisms.It is typically expected that, because of the lack of chaos in an integrable system, the transport is ballistic in nature.This conjecture is verified in the context of energy diffusion in integrable systems, for example, in Toda chain, Lieb-Linegar gas, sinh-gordan model etc. [6].However, exceptions can be found for example, in the integrable XXZ model which is integrable but the spin transport is non-ballistic [7].
Early attempts to describe large-scale dynamics included merging soliton theory with kinetic theory [9,10].The general picture of quasi-particles has recently been incorporated in the development of integrable hydrodynamics that address the full spatio-temporal structure of fluctuations and beyond [8].There is a serious effort to extend hydrodynamic theory to the quantum domain with growing interest in both experimental and theoretical perspectives [11,12].Curiously, experiments which contain few atoms are well described by the hydrodynamic theory initially developed in the thermodynamic limit.This simple question has no obvious answer and has deep connections to coarse-graining and taking the hydrodynamic limit [13,14].For example, the 1D Bose gas in the experiment [2] composed of 40 − 250 atoms may not seem enough for the thermodynamic limit, but the experimental results can be explained using integrable hydrodynamics.
In this work, we theoretically explore the question: can a small system be described by hydrodynamic theory?Our analysis is based on the classical integrable hydrodynamics of the classical Toda chain as developed in [15][16][17][18].The framework can be generalized to other integrable systems, classical or quantum [19].In the quantum domain, numerical simulations with time-dependent density matrix renormalization group of a thermal expansion in the XXZ model lead to a similar conclusion: hydrodynamic predictions with smooth initial conditions are accurate, even for small system sizes [10].
Here, we construct the theoretical framework to study spatio-temporal correlation in the Toda chain at equilibrium with three different but related approaches.The first uses a microscopic molecular dynamics (MD) simulation similar to [20], the second is by solving the generalized hydrodynamics equation for the Toda chain [15,16] and the third by solving the linear Euler scale equations considering all the conserved quantities.The three methods have different FIG.1: Schematic representation of increasing non-locality of conserved quantities of interacting integrable models with local interaction.The orange lines show the support of two sites for a local energy, while the grey line represents a larger support for higher order conserved quantity with support of four sites.levels of difficulty / accessibility using finite computational resources.The MD simulations can be performed for a relatively large chain (say N = 2048 particles) for ∼ 3000 computational hours, while solving the integrodifferential GHD equation in the thermodynamic limit is almost instantaneous, but requires a large time to setup the framework.Solving the Euler equation is limited to a few atoms only and is of medium complexity.En route to testing the accuracy of the three methods, we find that a small system with N = 8 particles has a spatio-temporal correlation at equilibrium almost indistinguishable from that obtained using the GHD equation in the thermodynamic limit.

I. EULER SCALE HYDRODYNAMICS
We discuss the classical many-body interacting integrable system in a general setting, formulate the conserved quantities, and give the prescription to construct the Euler scale hydrodynamics using their respective currents.
Integrable ingredients: Consider a system of N particles with position q i and momentum p i described by the Hamiltonian H({p i , q i }).The system is integrable if the equation of motion generated by the Poisson flow Ȧ = {A, H} can be written in terms of the Lax equation, where the dots refer to the time derivative.The common form of the Lax equation with the Lax pairs denoted by the N × N matrices L ± ≡ L ± ({p i , q i }) and M ≡ M ({p i , q i }), is Using the cyclic property of the trace, the global conserved quantities of the system are I k = Tr(L + L − ) k .Note that these conserved quantities are not unique and are defined to a constant prefactor.These are chosen such that the equilibrium correlation matrices for energy have an order of magnitude governed by the temperature T of the system.M is typically a skew-symmetric square matrix, i.e. a matrix whose transpose equals its negative.We are interested in determining the local variations of these globally conserved quantities that describe the current flowing in and out of a mesoscopic unit cell in the hydrodynamic regime.
Conserved quantities and currents: The local evolution of the conserved quantity (I k = i I k i ) can be written in the continuity form by observing that they follow the same evolution under the Lax flow.Associated with the evolution, locally conserved quantities I k i , is given as, where we define the local net currents flowing in and out of the site i as (M I k ) ii = j k i and (I k M ) ii = j k i+1 respectively.Typically for a H with only nearest-neighbour interaction, M is a two-sided off-diagonal matrix, resulting in the current flowing in and out of neighboring sites.However, for a long-range interacting system, the current flows in from all bonds.Another interesting observation is that the I k i are increasingly non-local in real space in a short-ranged system, i.e.I k i depend on the product of lattice variables in between i ± k, pictorially shown in Fig. 1.This will be important to understand our results later.The eigenvalues of the matrix L + L − is are also conserved quantities, but, however, they are not local in real space.
However, the eigenvalues of the L + L − matrix only capture some of the conserved quantities for translation-invariant Hamiltonian systems.There is an additional conserved quantity called the stretch variable, which is related to the inverse of the system density and plays a major role in transport.For the nearest-neighbour interaction, if the translation operator is defined as {T, H} = 0, then I 0 = i (T − 1)q i and the total momentum I 1 = i p i are the additional conserved quantity for systems with periodic boundary conditions.It is important to exercise some caution: Although these look like all the conserved quantities of the system, there can be others that go unnoticed, i.e. any function of the conserved quantity is also conserved.See, for example, [21] for the role in the non-linear Schrodinger equation also in the context of Toda lattice, which was investigated in [22].
General Gibbs ensemble (GGE): The integrable system at equilibrium is described formally by GGE, where µ k are the dual thermodynamic variables corresponding to each conserved quantity.For example, a system at Gibbs equilibrium is described by dual variable temperature (µ 2 = β, corresponding to Hamiltonian H ) and pressure (µ 0 = P , corresponding to stretch or the inverse density).In the rest of the text, any average .is understood to be taken with respect to the GGE characterized by β, P .
Euler scale (fluctuating) hydrodynamics at equilibrium: The coarse-grained large-scale space-time behaviour of the system describing the macroscopic evolution is governed by local fluctuations of the globally conserved quantities.We are interested in the system prepared in a GGE characterized by canonically conjugate thermodynamic variables such as temperature, pressure, etc.In equilibrium, each conserved quantity has a well-defined average that does not change with time, i.e. ∂ t I k i = 0.This contributes to the trivial shift of the conserved quantity of the system and is typically subtracted.The local fluctuation of the total conserved quantity The linearized hydrodynamic Euler equations for the evolution of u k (x, t) in the continuum limit are where we have switched from i → x in going to the continuum and j k (x, t) = j k (x, t) − j k x is the local current fluctuation associated with the kth conserved quantity.Following the standard prescription of hydrodynamics to the first order in small fluctuations, we express the current in terms of the local fluctuations of the conserved fields i.e. j k (x, t) = A kl (x, t, x , t )u l (x , t )dx dt , where the kernel A kl (x, t, x , t ) = δj k (x,t) δu l (x ,t ) encapsulates all information of the evolution.We have neglected the higher-order expansions which are not important in the ballistic scaling limit addressed in the current manuscript but become important for super-diffusive transport [1,23].The kernel is generally hard to compute and much simplification is obtained in the averaged current at equilibrium, where the time dependence can be dropped due to time translation invariance at equilibrium.The dependence on positions is dropped using spatial invariance and local interactions.
We switch to vector notation, with u the row vector consisting of the local variation of conserved quantities and A a square matrix with dimensions of the number of conserved quantities.The average current over the GGE at equilibrium is Here A = δj δu can be interpreted as the projection of the current to the conserved quantity and ζ a Gaussian The diffusion matrix D and the noise matrix B are added phenomenologically and cannot be easily derived from a microscopical picture.It may seem surprising that these terms are not zero here, as the chaos required for the origin of diffusion and noise is absent in an integrable system.Typically, the origin of these terms is argued as follows: The system is as a fluid cell where each cell is a 'mesoscopic' region, i.e. a region of finite extent, which is large compared to the distance between particles, but small compared to the macroscopic spatial variation scales L. The dynamics inside the fluid cell can be split into two parts: the first being the net current flowing through the fluid cell.Each of these currents are projected on the basis of conserved quantities as discussed before.The second part comes from the dynamics in the rapid variation in the fluid cell due to interactions with neighbour cells.This cannot be taken into account through the projection to the conserved quantities and is modelled as dissipation and noise.This gives the Euler scale hydrodynamic equation The above equation maintains conservation laws as in the microscopic system, keeping the total fluctuation across the system constant in time.Although this looks formal, a more rigorous derivation using large deviation theory can be obtained along the lines of [25].Further, it can be shown that (with abuse of notation, we drop the average) the matrix A = δj δu ≡ ju u 2 −1 , where we take time independent averages at GGE [39].The matrix ju αβ = j α u β − j α u β , and u 2 αβ = u α u β − u α u β does not have any time dependence and is computed at equilibrium.This has a nice geometrical interpretation since the change of curvature of the differential equation is interpreted as a projection of the current to the conserved quantities in the system.A more rigorous procedure will be to follow the methods presented in [15,25].Before going to the next section, where we discuss the ballistic correlation, a small remark is that for a non-integrable system like the Fermi-Pasta-Ulam (FPU) chain with three conserved quantities, Eq. ( 6) consisting of 3 × 3 matrix corresponding to the stretch, density and energy conservation.In this case, truncating the local current to linear order is not enough and one needs to keep higher-order terms.This gives rise to evolution with non-ballistic scaling and is described by the non-linear fluctuating hydrodynamic theory.
Ballistic spread of correlations: We are interested in spatiotemporal correlation functions defined as C(x, t) = u(x, t)u T (0, 0) which follow the evolution as Eq. ( 6).It is easy to check that the total correlations dictated by this evolution is constant, i.e. ∂ t x C(x, t) = 0 as the system undergoes conservative dynamics.In general, the late-time evolution of the correlation functions follows a self-similar scaling from C(x, t) = 1 t γ f ( x t γ ).It is expected [although not obvious] that since our system is integrable, the transport will be ballistic, that is, γ = 1.Using this in Eq. ( 6), we see that the contribution of the diffusive term is subleading for ballistic scaling, and we can safely ignore the last term and are left with Solving the correlation equations: The evolution of the correlations is now a straightforward linear equation and using the Fourier transformation f The solution is then written as where C(k, 0) is the Fourier transform of the equal time correlations between the conserved quantities prepared in GGE.The matrix A is not symmetric, however the eigenvalues are always real as expected physically from the spatiotemporal symmetries of the system.This can be undertood as follows, when the equilibrium system is flipped spatially and its currents are time-reversed, it behaves identically to the original system.This requires that the eigenvalues of the matrix A are real.The matrix C is a symmetric matrix.The inverse Fourier transform then gives the time-dependent solution Taking the continuum limit and scaling the correlation functions, this can be rewritten in scale-free form as It is convenient to go to a diagonal basis of A, with RAR −1 = diag[λ], where R diagonalize the matrix.The eigenvalues of the matrix A are physically relevant quantities which can be interpreted as generalized velocities of ballistic propagating modes.In FPU-like non-integrable systems, with three conserved quantities, the matrix A is for example the 3 × 3 matrix: the eigenvalues give the sound speed of the system along with a zero eigenvalue, which corresponds to the stationary heat peak of the system.In this case, the Euler approximation does not reveal the structure of the peaks itself, and one must go to higher order expansions in the current to reveal the structure of the peaks.In an integrable system with ballistic scaling, the situation is interesting, as there is only one timescale.The system of equations gives rise to N nonzero velocities corresponding to the number of conserved quantities in the system.The eigenvalues are ordered in the sense that λ 1 ≤ λ 2 • • • ≤ λ N and form a linear space.The spectral gaps between the eigenvalues are smaller toward the edges and larger in the bulk.Taking these into account renormalizes the amplitude of the scale-free function at the k th site by λ k+1 − λ k .The scale-free equation Eq. ( 10) can be written in terms of these velocities as, The normalization can also be understood from a more physical point of view: Since the integrable dynamics of the system is conservative, normalization eliminates the gauge degree of freedom from the eigenvectors, and we have FIG.2: Comparison of scaled equilibrium correlations in a small finite sized Toda lattice with thermodynamic limit results.The circle and the stars are data from numerically solving the Euler scale hydrodynamic equation at Gibbs equilibrium as in Eq. 11 for small system sizes with N = 8 and N = 16.The noisy grey line is microscopic molecular dynamics simulation of the macroscopic number of N = 2048 particles with Toda interaction with an initial condition average of 10 7 initial conditions from the Gibbs distribution with T = 1, P = 1.The black dashed line is obtained by numerically solving the thermodynamic GHD equation.We see that the results of the finite short system already have an excellent overlap with the exact results in the thermodynamic limit.
the sum rule discussed above, z f (z)dz = C(0, 0), where the dz = λ k+1 − λ k are the discrete inverse density of the ordered eigenvalues and the right-hand side is the value of correlations computed at equilibrium.The density fluctuations are studied under the "unfolded" spectral gap of the eigenvalues of the integrable system akin to that used in random matrix theory.A discussion of this formula for the Toda lattice is shown in Sec.III.As we shall see, the above formula for a small system size, i.e. when N is of the order of O(1) reproduces the thermodynamic limit spatio-temporal correlations.Note that we have not taken any hydrodynamic or large system size limit in the above formula.This makes us wonder if the dynamics of small systems at equilibrium are sufficient to be described by hydrodynamics.

II. GENERALIZED HYDRODYNAMICS: GHD
Generalized hydrodynamics assumes that the dynamics of an integrable system can be described by a quasi-particle which undergoes ballistic diffusion and a drift with effective velocity.At thermal equilibrium, this is obtained by averaging the extensive number of conserved charges and their dynamics in a local GGE state.The natural hydrodynamic fields are the density of states (DOS) of the Lax matrix and the stretch, r, which we anticipated earlier.
For completeness, we first briefly review the basic concepts of the GHD in classical systems following the notation introduced in [15,26].The DOS is represented as rρ p (r), where the distribution ρ p encodes all information of the conserved quantities.Then the generalized hydrodynamic equations are where q 1 = r dw wρ p (w) is the average momentum of the system.The effective velocity is given by the solution of the linear integral equation of the form with the operator T being the integral of kernel (K) which takes into account the two-particle scattering shift.It is defined as an integral transform on a function ψ(w) as, where the kernel K depends on the microscopic model and can be computed by solving the large time asymptotic in two-particle scattering.
It is convenient to use a non-linear transformation to go to a coordinate frame, where the hydrodynamic equation becomes simplier.This is done by transformation to normal mode DOS ρ n , Key to this description is to introduce a dressing transformation that transforms the real valued function ψ as In this notation, the transformation between the normal mode and the effective velocity is given as and the normal mode hydrodynamic equation is It is convenient to denote the average over DOS as f p = f (w)ρ p (w)dw.In this notation, the conserved quantities are given as the integral identities r p = 1 and rw k p = I k .We are interested in the correlation function of the conserved quantities, Thermal states: For thermal states, the ρ n is the solution of the thermodynamic Bethe-Ansatz (TBA) like equation characterized by the chemical potential µ a function of the external pressure and temperature determined at equilibrium by condition n(w)dw = P , where V (w) = 1 2 βw 2 for the initial state in Gibbs equilibrium.Note that for a classical system there is no TBA, but this equation is obtained in Toda chain in a similar spirit using a nontrivial mapping to random matrix theory with the kernel given by K(x) = log(x), which is determined from the two-particle scattering.The solution to this highly non-linear root equation is often supplemented with solving the corresponding Fokker-Planck equation whose steady state solution satisfies the TBA equation above, where K denotes the derivative of the integral kernel K which is specific to the model.As a consistency check, the solution to Eq. ( 21) is inserted into the TBA Eq. ( 20), gives the chemical potential µ.At equilibrium, µ turns out to be exactly the free energy when we set w → 0 which acts as a cross-check for the results with the explicit formula available [20,30].At this point, it is convenient to introduce the quasi-energy parameterization of the normal mode obtained from the solution of the Fokker-Planck solution of TBA by (w) = − log[ρ n (w)].Numerically, this step is carried out by first obtaining from the Fokker-Planck equation and then plugging it into the TBA Eq. ( 20) and finding the roots of the equation.The average momentum in terms of dressed variable is The correlations in Eq. ( 19) in case of momentum correlation is conjectured as [Eq.3.17 in [15]], Other higher-order correlations can be computed appropriately by generalising this formula.

III. TODA CHAIN
We now focus on applying the methods described on the Toda chain [31] which is a model of classical particles with exponentially decaying nearest-neighbour interactions.This is an example of a classical interacting integrable particle system [4] which is solvable by Lax pairs as described before.A large Toda lattice is used as a discrete model for solitonic waves while a small chain of three particles can be mapped to the famous Hennon Heiles problem.The history of understanding dynamics in Toda lattice is decades old: The equilibrium dynamical properties in large Toda chains were extensively studied using various methods.Numerically, correlations have been explored in the literature in [20] using molecular dynamics simulations.In an integrable system of N particles, an equivalent number of conserved quantities is enough to describe the dynamics completely.[40] The Toda chain is defined as 2 + e −ri with r i = q i+1 − q i .The equations of motion (e.o.m.) are: These e.o.m. along with periodic boundary conditions r N +1 = r 1 and r 0 = r N can be cast to a Lax matrix as in Eq. ( 1), with λ = 0 and L = L + = L − .The symmetric matrix L is a tri-diagonal matrix with diagonal elements b i = pi 2 and off-diagonal elements a i = 1 2 e −ri/2 .The conserved quantities are trace of the powers of Lax matrix, denoted as N i=1 e i .In [24], explicit higher-order conserved quantities are remunerated.Importantly, in real space, the conserved quantities have increasing non-local support in space(see Fig. 1).However, this is not the only way to view the conserved quantities of the system; for example, the eigenvalues of the Lax matrix are also conserved under evolution.It has been suggested that the spacing between the eigenvalues of the Lax matrix provides a natural transition to normal mode frequencies in [32,33].These conserved quantities are an adiabatic invariant for small non-integrable perturbations [34].Prior attempts have been made to understand the dynamics of the Toda chain.The inverse scattering method is useful in an isolated system and provides a way to address not only the solitonic properties but also the large-scale hydrodynamic properties, as discussed earlier.Thermalization properties and the conundrum of the local observable are related to Gibbs or General Gibbs is addressed in [35].The majority of the works address single-point functions and finite temperature dynamical spatio-temporal properties although difficult to access has been attempted through non-interacting soliton gas analogy [36] and through quantumclassical correspondence with quantum Toda chain [37].Earlier attempts at dynamic properties of the Toda chain are reviewed in [38].After a gap in interest in Toda dynamics, recent development in the integrable hydrodynamic theory of Toda chain in [1,16], where the dynamics of single-point functions starting from domain wall initial conditions was addressed [26].Accurate numerical results for the dynamics of two-point functions were discussed in [20] and analytical results were only available in the limiting cases when the model was solvable.As a by product, this work fills the gap in theoretical description of the numerical results in spatio-temporal correlation functions described in [20] from the theory developed in [15,16].
Results: Fig. (2) summaries the main numerical results of computing the momentum spatio-temporal correlations at Gibbs thermal equilibrium.The system described by Eq. ( 24) is solved together with the initial conditions at equilibrium described by temperature (β −1 = T ) and pressure (P ) in a Gibbs ensemble with Z = i e −β(ei+P ri) .The noisy grey line is a microscopic numerical simulation of the macroscopic number of N = 2048 particles with Toda interaction with an average of 10 7 initial conditions sampled from the Gibbs distribution with T = 1, P = 1.The black dashed line is obtained by numerically solving the GHD equation in Eq. ( 23) with the input from the TBA equation for computing the effective velocity of the quasi-particle and analytically computed thermal expectation values.Exact numerical details of the relevant computation is given in [26], and a code is available at [28].The circle and the stars are data from numerically solving the Euler scale hydrodynamic equation at Gibbs equilibrium as in Eq. (??) for small system sizes with N = 8 and N = 16.As anticipated earlier in the introduction, surprisingly, even for the small system sizes, the spatiotemporal correlations have a good agreement with correlations in a large system sizes.This is indeed a surprising observation, and it make us wonder about the question of taking a thermodynamic limit, which serves as the motivation of the problem.Note that the positions for each of the dots corresponding to solving the Euler scale equation can be interpreted as "generalized" velocities of the system in the ballistic scaling limit.When the dimensions of the Lax matrix are of order N , there are N + 1 velocities, which overlaps with the thermodynamic solution of the GHD.It is interesting to note that these generalized velocities are neither symmetrical nor are they equally spaced.The peak of the curves are not captured well with this method, however, the possibility of the effects of taking of even-odd number of particles cannot be ruled out.The relative density of state of the distribution of the Lax eigenvalues are shown in Fig. (3), which shows that the eigenvalues are closer to the edges and more spread out in the bulk.
Numerical method for direct Euler scale dynamics: The numerical method used to simulate the hydrodynamics of the Euler scale is by Monte Carlo sampling of the Lax matrix.The trick here is to do the Monte Carlo sampling from the Lax matrix itself corresponding to the Gibbs ensemble of the system at a specific temperature and pressure.For the Toda system, it can be shown that in Gibbs distribution: i.e. when the system is prepared in P ({p i , r i }) = Z −1 e −β(H({pi,ri})+P ri) (Z is the normalization), the elements a i follow the chi-square distribution (Fig. 4 for numerical sampling), P a (x) = 2(4β) βP Γ(βP ) x 2βP −1 e −4βx 2 , and the momentum follow the Gaussian distribution.The eigenspectrum of the Lax matrix is used to compute the Euler-scale hydrodynamics along with the equilibrium correlations of conserved quantities.In our simulations, we used a 10 9 sample of the Lax matrix and checked the convergence of the eigenvalues of the matrix A. The equilibrium correlations of the conserved quantities are more difficult objects to deal with analytically.As we mentioned before, the higher-order conserved quantities have increasing non-local support in real space, and the convergence of equilibrium correlations becomes difficult.However, here we are interested in the system at Gibbs equilibrium, where the relevant conserved quantities that enter the distribution are energy and stretch, which are both extremely local in space.This makes the connected correlation of the higher-order conserved quantities negligible compared to the correlations of energy and the stretch variables.Interestingly, truncating the Lax matrix to include only the first three conserved quantities which appear in the Gibbs ensemble: the energy, momentum, and the stretch reproduce to 5% error on the motion of the peaks of the dynamical fluctuations [20] with the non-zero eigenvalues λ = ±0.8833and we conclude that the correlations of the higher non-local conserved quantities are important to capture the dynamical correlations of the more local conserved quantities.

Comparison and discussion and limitations of results
We reiterate that these results are strictly for correlations in the ballistic scaling limit.There are some intuitive results on the nature of diffusive corrections on the ballistic scale, which we will not discuss here [16].Although the method seems general, it is by no means more efficient numerically than solving the GHD equations directly.The GHD framework provides a general analytic framework for working with set generic initial conditions.On the other hand, the Euler scale hydrodynamics is a more brute-force approach to the problem using Monte Carlo sampling of the initial data and is limited to small system sizes.

IV. CONCLUSIONS
We have discussed the adequacy of using integrable hydrodynamics to describe small classical (and quantum) systems.We show that although the hydrodynamic theory is developed for a large number of particles in thermodynamic limit, practically, using an independent direct method discussed here, a relatively small system size in a classical interacting integrable system is already quite accurate in describing the results in thermodynamic limit.We have tested this for a short-ranged interacting system; however, questions related to the finite-size hydrodynamics in a long-ranged system (for example in Colegaro Sutherland model) remains to be explored.
Density of states:The relative density of states of the finite-dimensional Lax matrix for Toda used to calculate the correlation function is shown in Fig.3.

FIG. 3 : 4 .FIG. 4 :
FIG.3:The denisty of states 1/(λi+1 − λi) of the distribution of sorted eigenvalues of the Lax matrix for the Toda lattice in Gibbs equilibrium at temperature T = 1 and pressure P = 1.The sampling of the Lax matrix at equilibrium is achieved through transformation of the probability distribution.The diagonal elements of the Lax matrix are Gaussian to sample while the off-diagonal elements is shown in Fig.4.Monte-Carlo simulation of the Lax matrix: