Thermalization and hydrodynamics in an interacting integrable system: the case of hard rods

We consider the relaxation of an initial non-equilibrium state in a one-dimensional fluid of hard rods. Since it is an interacting integrable system, we expect it to reach the Generalized Gibbs Ensemble (GGE) at long times for generic initial conditions. Here we show that there exist initial conditions for which the system does not reach GGE even at very long times and in the thermodynamic limit. In particular, we consider an initial condition of uniformly distributed hard-rods in a box with the left half having particles with a singular velocity distribution (all moving with unit velocity) and the right half particles in thermal equilibrium. We find that the density profile for the singular component does not spread to the full extent of the box and keeps moving with a fixed effective speed at long times. We show that such density profiles can be well described by the solution of the Euler equations almost everywhere except at the location of the shocks, where we observe slight discrepancies due to dissipation arising from the initial fluctuations of the thermal background. To demonstrate this effect of dissipation analytically, we consider a second initial condition with a single particle at the origin with unit velocity in a thermal background. We find that the probability distribution of the position of the unit velocity quasi-particle has diffusive spreading which can be understood from the solution of the Navier-Stokes equation of the hard rods. Finally, we consider an initial condition with a spread in velocity distribution for which we show convergence to GGE. Our conclusions are based on molecular dynamics simulations supported by analytical arguments.


Introduction
A classical Hamiltonian many-body system will generally thermalize at long times in the sense that macroscopic obeservables can be described by the Gibbs Ensemble (GE).However, there may exist systems that do not thermalize to GE, because of the existence of macroscopic number of extra conservation laws which restrict their motion in the phase space.Such systems are often known as integrable manybody systems, which are believed to thermalize to the Generalized Gibbs Ensemble (GGE) [1,2,3,4].
They have been realized experimentally in one-dimensional trapped atoms [5,6].Their non-equilibrium dynamics close to local GGE is described by generalised hydrodynamics (GHD) [7,8,9,10,11].Integrable systems are very fine-tuned systems in the sense that the smallest of perturbations (which are always present in any experimental setup) can break integrability.However, in the presence of integrability breaking perturbations, it is expected that the system will still remain integrable for short times [5,6,12,13], and so integrable dynamics may still play an improtant role.Integrable systems are also important from the point of view of studying exact dynamics of systems far from equilibrium.Since they rekindle the hope of obtaining exact solutions to many-body systems out of equilibrium, it is possible to use them to study far from equilibrium states, which cannot be treated using hydrodynamics (for non-integrable case) or generalised hydrodynamics (for integrable case) because hydrodynamics (HD) can only handle states near local equilibrium (or local GGE).
It is useful to make a distinction between interacting and non-interacting integrable systems [14].
In non-interacting integrable systems, the quasiparticles move in straight lines at constant velocity.For example, in one dimensional hard point particle gas, the collisions happen at a point and two particles simply exchange velocity after colliding, and thus the system can be mapped to non-interacting one by interchanging the labels of the two colliding particles after the collision.In the mapped non-interacting problem, the new particles moving with a fixed velocity are called quasi-particles, and they move in straight lines at constant velocities, like in a non-interacting gas.This is not the case for interacting integrable systems.In the hard rod case, a quasiparticle will have a straight line motion interrupted by sudden jumps (of the size of rod length) owing to the collisions.This effectively leads to dissipation in the hydrodynamics of the hard rod gas, which can get manifested by the spreading of a tagged quasiparticle.
On the other hand, for the point particle case, the system has no dissipation term in its hydrodynamics and consequently no spreading in the position of a tagged quasiparticle.Such spreading was studied, from microscopic calculations, by Lebowitz, Percus and Sykes [15] and demonstrated the effect of dissipation.
Such dissipation terms appearing as Navier-Stokes (NS) corrections to the HD equation of hard rods was later established by Spohn [16], Boldrighini and Suhov [17] and recently discussed by Doyon and Spohn [18] and Ferrari and Olla [19].Due to the presence of the dissipation term, one generally expects that the hard rod gas would approach to a GGE state starting from a non-equilibrium initial condition.
The question of approach to the GGE state in integrable systems has been widely discussed in the quantum context [1,20,21,22] and the effect of dissipation was demonstrated in the context of evolution of a domain wall in the quantum Heisenberg spin chain [23].However, to the best of our knowledge, this has not been observed for classically integrable systems.Neither has the effect of the Navier-Stokes correction to the Euler GHD solutions been demonstrated in any study.In the context of the classical system of hard rods, the questions on evolution and effect of dissipation were addressed in [18] for the specific case of domain wall initial condition.This study demonstrated that the evolution from such initial condition can be very well accurately described by the solution of the Euler GHD equations.Although the corrections from the Navier-Stokes terms were discussed in [18], this could not be unambiguously established from the numerics.The aims of the present paper are: (i) to study the evolution of nonequilibrium initial states and see if they approach GGE at large times; (ii) to demonstrate the effect of dissipation in such an evolution.This paper is organized as follows.In Sec. 2 we define the model and the different initial conditions used in the study, and summarize the main results.We investigate the equilibration and the effect of dissipation in Sec.4.1 by comparing the predictions of hydrodynamics with those of MD simulations for different initial conditions.In Sec. 5 we provide discussions of our results and conclude.Some details of the calculations are provided in the appendix.

Model, observables and initial conditions
We consider a system of N hard rods each of length a and unit mass, moving inside a one dimensional box of size L. The rods move with constant velocity in between collisions.Two rods exchange their velocities at collisions with each other whereas at collisions with the walls at x = 0 and x = L the rods flip their velocities.This implies reflecting boundary conditions at x = 0 and x = L.This model with a ≠ 0 is an example of an interacting integrable system, while for a = 0, it becomes non-interacting integrable system.
The microscopic dynamics of hard rods can be mapped to that of hard point particles as follows [24,25,15] and v i , i = 1, 2, ..., N denote velocities of the rods.For each microscopic configuration {x i , v i } of hard rods, one can construct a configuration {x ′ i , v ′ i } of hard point particles by removing the inaccessible spaces between rods and, between rods and the walls.More precisely the mapping can be written as and consequently one has a set of hard point particles moving inside a box of length L ′ = L − N a.The dynamics of hard point particles can be further mapped to non-interacting point particles by the Jepsen mapping [26].This mapping has earlier been used to find several analytical results, such as quasiparticle distribution [15], free expansion problem [25] and sound and shock propagation [24].This mapping also allows one to simulate the hard rod dynamics efficiently and accurately.Throughout this paper, we represent configurations of the point particles by primed variables ({x ′ i , v ′ i }) and those of the rods by un-primed variables ({x i , v i }).
In this paper we study the evolution of the single particle phase space distribution, f (x, v, t), of the hard rods defined as where ⟨...⟩ denotes an average over the ensemble of initial conditions corresponding to fixed forms of the initial density profile and single particle velocity distribution.We investigate the possible approach to GGE and the effect of NS terms, for the following three different initial conditions: A The particles on each of the two halves, (0, L/2) and (L/2, L), are separately distributed uniformly at a finite density ρ 0 = N /L.We assign all particles on the left half [x ∈ (0, L/2)] with velocity v 0 (= 1), while the velocities of all particles on the right half [x ∈ (L/2, L)] are chosen from a Maxwell distribution with temperature T = 1 In this initial condition, one has two components of the gas -for the first component, each particle has velocity v 0 = 1, while in the second (which we call the background particles), the particles have velocity distributed according to h(v).
B Next we consider an initial condition where we first place a particle with a given velocity v 0 (=1) at the origin.The rest of the box is then filled with particles distributed uniformly in space with a density ρ 0 = N /L.These background particles have velocities chosen from the Maxwell velocity distribution h(v).This case is more analytically tractable than the previous case in (A) and was first studied by Lebowitz, Percus, Sykes (LPS) in [15].In this case, the initial single particle phase space C Finally we consider the set-up of free expansion from a half filled box.In this case the rods are uniformly distributed in the left half of the box at a constant density, 2ρ 0 , and the velocities are again chosen from the Maxwell distribution, h(v).The right half of the box is empty.This problem of free expansion was previously investigated in [25], where the evolution of various hydrodynamic variables was computed using a microscopic approach and with certain approximations that effectively amount to solving the Euler equations.For both the classical and quantum cases, the free expansion problem for point particles has recently been studied in the context of entropy growth [27,28,29].
For the three initial conditions mentioned above, we study the evolution of the density profile, ρ(x, t), and the velocity profile, u(x, t) (or equivalently, the momentum density profile p(x, t) = ρ(x, t)u(x, t)), defined as We investigate the effect of dissipation by comparing the profiles obtained from simulation with those predicted from the solution of Euler GHD.We also check if the system reaches GGE at long times.A simple test for GGE would be to check if the density and velocity profiles become time stationary, i.e., independent of time and uniform in space.Note that the velocity distribution is invariant under the integrable dynamics and specifies the GGE.We summarize here our main findings: -For initial condition (A), we find that the predictions for the evolution of the densities, from the solution of the Euler GHD, describe the profiles obtained from numerical simulations quite well almost everywhere except at the locations of the shocks where we observe clear discrepancies.These discrepancies appear due to the effect of dissipation described by the NS term.This leads to the width w(t) of the shock growing as w(t) ∼ √ t at early times, and saturating to a value w(t → ∞) ∼ √ N at large times.Thus the initial density of either of the components (v 0 = 1 and the thermal ones) never become homogeneous over the full system and each of the two components move inside the box with a constant effective speed v eff (see Eq. ( 38)) at large times.Hence the system never reaches GGE.
-This spreading is more prominent in the case of initial condition (B).For this case we provide an analytical understanding of the spreading based on the solution of the NS equation.In this case also the system never reaches GGE.
-For the case (C) of free expansion we find that the evolution of the density and momentum density profiles is completely described by Euler GHD and, since the discontinuity in the initial density profile disappears already at very early times, any effects of the NS corrections are too small to be observed.
In this case, the system at long times evolves to a state consistent with GGE.
3 Hydrodynamic equations for hard rods and solution of the Euler equation The hydrodynamic equation for the single particle phase space distribution f (x, v, t) for the hard rod gas is given by [18]: and ρ(x, t), u(x, t) are given in Eq. (4a) and (4b).The term ∂ x N represents the NS correction to the Euler equations [18,17].
We now discuss the solution of the Euler equation for general initial condition.As shown previously [25,24], the Euler equation for hard rods can be solved exactly for general initial conditions by mapping it to a non-interacting point particle problem.For completeness, we show below how the mapping to the non-interacting Euler equation can be obtained using the GHD approach.For this one defines a new function, where and B is the position of the left end of the container in which the hard rod fluid is contained.Note that F (x, t) is the cumulative density corresponding to ρ(x, t).We observe that which implies that f 0 (x ′ , v, t) is also a phase space distribution function.We now show that f 0 (x ′ , v, t) satisfies the Liouville equation of free ballistic particles and hence describes the single particle phase space distribution of the point particles.The first step towards this demonstration [10] is to define the function Using Eq. ( 6) and the relation ∂ t ρ = −∂ x (ρu) [for the fields defined in Eqs.(4a,4b)], it readily follows that Now, from Eq. 7, we have Taking the time derivative with respect to t on both sides of Eq. ( 12) and using ∂ t F (x, t) + ρu = 0 one finds On the other hand taking derivative with respect to x on both sides of Eq. ( 12) one has Inserting the forms from Eq. 13 and Eq. 14 in Eq. 11, one finds that the phase space distribution function f 0 (x ′ , v, t) satisfies the Liouville equation for the non-interacting particles This equation can be easily solved for arbitrary time and any initial condition f 0 (x ′ , v, 0).For example on the infinite line one has f 0 (x ′ , v, t) = f 0 (x ′ − vt, v, 0) while in the box one has to solve the single particle problem with repeated collisions with the walls [25,27].From f 0 (x ′ , v, t) one can find the solution for the phase space distribution f (x, v, t) of hard rods.To get the solution explicitly, we first note from Eq. ( 7) , where ρ 0 Hence inverting Eq. ( 7) and using Eq. ( 16), one finds The variable transformation x → x ′ can be inverted as using F (x, t) = F 0 (x ′ , t) which can be shown easily from Eq. ( 9).
While, as demonstrated above, the Euler equation can be solved exactly, it is difficult to solve the NS equation (5b) for arbitrary initial conditions.We expect that the difference between the solutions of the Euler and the NS equations are large at places where the spatial derivative of the Euler solution is large.
4 Results from numerical simulations for the three initial conditions

Initial condition A
In this case the initial condition can be written explicitly as where h(v) is given in Eq. (3) and with ρ 0 = N L−a and Θ(x) being Heaviside theta function.Note that we will be working in the thermody- As discussed in the previous section, solution to the Euler equation can be obtained by mapping to point particles.It is easy to show that initial phase space density f 0 (x ′ , v, 0) also has two components, the special component with velocity v = 1 and the background particles having velocity distributed according to Maxwell distribution.It is given explicitly as where L ′ = L − N a and The evolution of f 0 for any arbitrary initial distribution Since the Euler Equation ( 15) for point particles is linear, the distribution f 0 (x ′ , v, t) at time t still can be written as a sum of two components as The first term in this expression can be obtained by putting ϱ(y) = g 0 (y, 0) and p(u) = δ(u − 1) and performing the integration over v, yielding where, recall L ′ = L − N a.For the background component we set ϱ(y) = ρ 0 b (y, 0) and p(u) = h(u) to get: Using the Poisson resummation formula, this can be rewritten in the alternative series form: Note that shifting the origin to L ′ /2 (i.e., x ′ → z ′ = x ′ −L ′ /2) and taking L ′ → ∞, one obtains the solution of Euler GHD on the infinite line as where T = 1.The corresponding densities of the hard rods for the two components, respectively g(x, t) and ρ b (x, t), can be obtained using the inverse mapping Eq. 16 and Eq. ( 18) along with ρ 0 (x ′ , t) = g 0 (x ′ , t)+ρ 0 b (x ′ , t).We show in Fig. 1 the evolution of g(x, t) and ρ b (x, t) obtained from the solutions of the Euler GHD equation as well as the results from direct MD simulations.We see that even for very long time like t = 640, the profile obtained from MD does not relax to GGE, i.e., it does not become uniform.We also see that there is a discrepancy between MD and Euler solutions at the shock front due to dissipative effects.We have taken length of the box L = 2500, total number of particles N = 2000 and length of rod a = 1.0.We have performed ensemble averaging over 5000 realizations while doing MD.
For the solutions of the Euler GHD, we make the following observations: a.There is always a shock at the front of the density profiles for both the components.On the infinite line, the shocks for the two components move in opposite directions.Note that the density profiles g 0 (x ′ , t) and ρ 0 (x ′ , t) in the point particle gas evolve independently of each other.Consequently, g 0 (x ′ , t) will move with constant speed v 0 = 1 keeping the initial shape unchanged, i.e., with two discontinuities at L/2 separation.Hence the total density, ρ 0 (x ′ , t) = g 0 (x ′ , t) + ρ 0 b (x ′ , t) will also have discontinuities.Consequently, the density profiles g(x, t) and ρ(x, t) of the hard rods, obtained through the transformation in Eq. ( 16) also exhibits discontinuities, i.e., shocks.b.At early times the evolution of these density profiles correspond to that on an infinite line and can be described by g(x, t) and ρ b (x, t) obtained after transforming the solutions given in Eqs.(31) for the Euler equation of the point particles.c.At later times, each component of the gas gets reflected from the walls of the box which are described, in the point particle picture, by various terms in the series in Eq. ( 27) and Eq. ( 29).d.At the longest times both density profiles g(x, t) and ρ b (x, t) stop broadening further and settle to piece-wise flat profiles which move between the walls with some constant effective velocity v eff (see Fig. 2).The details of this solution will be discussed below.Since the density profile does not become time stationary even at the largest times, this indicates that for initial condition A, the hard rod system will never reach a GGE state ( which should be time stationary).Fig. 2: Lack of thermalization to GGE: Plot of g(x, t) vs. x at different (late) times.The dashed lines are obtained from molecular dynamics and the solid lines represents the solutions of the Euler equations.We observe that the profiles at t = 1000, 1100 and t = 1200 have moved by a displacement ∆x ≈ 300, implying v eff ≈ 3, in agreement with Eq. ( 38).The width of the pulse at the times t = 1000, 2000 are the same, thus indicating that it saturates and the whole profile does not become uniform, i.e., it does not relax to a GGE form.The inset shows a zoom of the shock at the two times t = 1000 and t = 2000, where we see that its width has saturated.We chose t = 1000 and t = 2000 as times for which the profiles coincided for the particular parameter values (i.e., a, N, L).Here N = 2000, L = 2500, a = 1 and ensemble averaging over 100 realizations were performed while doing MD.The values of g 1 , ρ b , ρb and v eff agree with the predictions in Sec.4.1).
While we see a very good overall agreement between Euler solution and MD simulations, there are clear differences.If we zoom near the shocks in Figs.1a and 1b, we notice that the simulation data for the hard rod density profile g(x, t) (dashed lines) shows a slight discrepancy with the Euler prediction.One observes similar discrepancy for ρ b (x, t) also.The simulated profiles display spreading at the locations of the shock in the Euler solutions.This is demonstrated in Figs.(3(a),4), where the density profiles are zoomed near the shock location after shifting appropriately so that the shock positions coincide.This spreading is a signature of the dissipation characterised by the Navier-Stokes term in Eq. ( 5).We observe that the width of the shock increases with time and scales as √ t as can be seen from Fig. 3 that a shock, of the Euler solution for g(x, t), encounters till time t.This fluctuations arise from the fluctuations in the initial conditions.For a given initial configuration of the positions and velocities of the rods, the shock remains sharp and does not widen.However, the place at which the shock appears at a given time fluctuates from one initial microstate to another (see Fig. 3(b)).This happens because the number of background rods that the special rods encounters is different for different initial microscopic configurations.Hence, on an average the shock widens.At small times, these fluctuations are independent as the rods have not realised the presence of boundaries of the box.The √ t growth at small times can be explained by considering the evolution of the density profile starting from initial condition A on an infinite line which is done in Appendix A.
The early time growth of the width of the shock stops after some time and saturates to a O( value as demonstrated in Fig. 4. As time progresses the rods move back and forth inside the box and consequently, the fluctuations in the number of background rods inside the region of the special rods (having velocity v 0 ) do not remain independent and get correlated.Consequently, the spreading of the shock cannot continue to grow as √ t and saturates to the observed O( √ N ) value.Thus even in the thermodynamic limit the pulse g(x, t) does not spread to the full extent of the system and remains in the shape of a rectangular pulse that keeps on moving back and forth inside the box.Consequently, the total density profile of the rods does not become homogeneous and stationary as one would expect in a GGE state.This implies that the a hard rod system inside a box, starting from initial condition A, does not reach the GGE state even in the thermodynamic limit.
Euler solution in t → ∞ limit: We now find the solution of the Euler equation in the t → ∞ limit.Recall that in the initial condition A, the rods are uniformly distributed in each half with density ρ 0 .The rods on the left half have velocity v = 1 and those on the right half have velocities distributed according to Eq. ( 3).Using Eq. 16, one maps this hard rod system to a point particle gas with uniform density ρ 0 (x, 0) = ρ0 1−aρ0 inside a smaller box of size L ′ = L − N a.The velocity distribution remains unchanged as in the hard rod gas, i.e., δ(v − 1) in the left half and h(v) in the right half.In the point particle gas the component with velocity v = 1 moves without changing its shape whereas the particles on the right half (called the background particles) perform free expansion, ignorant of the v = 1 particles since the gas is non-interacting.At long times, the background point particles will expand into the full box of length L ′ and become uniform with density half of their initial density, i.e., ρ 0 (x, t → ∞) = ρ 0 (x,0) 2 = ρ0 2(1−aρ0) .Thus at long times, one would observe the initial density pulse of the special point particles with velocity v = 1 moving in the uniform background of thermal particles (with Maxwell velocity distribution).Hence, at any instant, the total density profile has two regions: a uniform high density region where the v = 1 pulse is present (we call it the pulse region) and a uniform low density region in the remaining part of the box.Thus, the total density profile in the long time limit becomes piece-wise uniform which we now proceed to compute.The curves have been shifted so that the shock fronts for all the curves coincide.While doing MD, ensemble averaging over 5000 initial microstates were performed.We see that there is trend of increasing width with time while in the inset (which shows curves for t = 20, 40, 80), we see that there is a scaling collapse in the variable x/ √ t for short times when the v = 1 pulse does not know about the boundaries of the system and hence behaves like it is in an infinite system.We explain this √ t dependence in Appendix A. (b) This shows plot of g(x, t) for two different realizations (microstates) for initial condition A. We see that in a single realization the shock remains sharp, while the positions of the shock front in two realizations are different.Consequently, ensemble averaging will lead to smearing of the shock and thus is necessary to observe dissipation.Here we have chosen N = 8000, L = 10000.Fig. 4: System size dependence of the shock width for initial condition A: This shows the structure of the shock at late times (when the width saturates) for the initial condition A for different system sizes.For all the curves we have chosen t = 10000 which is much longer than the time at which the width of the pulse g(x, t) and that of the shock saturates.Even after this long time, the curve g(x, t) does not become uniform, i.e., it does not relaxes to GGE.The curve has been shifted so that the shock fronts for all the three curves coincide.We see that the shock broadens with system size, while in the inset, we see that there is a scaling collapse in the variable x/ √ N , thus showing that the shock broadens with the system size as √ N .In this case ensemble averaging over 500 realizations was performed.
Let us denote the value of the density of the v = 1 particles inside the pulse by g ′ 1 , in the point particle picture, and by g 1 in the hard rod picture.Similarly, we denote the density of the background particles inside the pulse region by ρ ′ b and ρ b , respectively, in the point particle and hard rod pictures.We also denote the density of background particles outside the pulse region by ρ′ b and ρb , once again, in the point particle and hard rod pictures, respectively.The total density of point particles in the pulse region is . Hence the total density there is 3ρ0 2(1−aρ0) .The density outside the pulse region is given by ρ′ b = ρ0 2(1−aρ0) .Now using the inverse mapping in Eq. ( 16) along with Eq. ( 18), one gets the late time densities in the hard rod picture.The density outside the pulse region is given by ρb and the total density inside the pulse region is given by To find individual values of g 1 and ρ b we use the conservation of the number of background particles where L 1 is the length of the pulse region at late times and L − L 1 is the length of the region outside the pulse.It is easy to see that L 1 = N 2g1 .Dividing both sides of Eq. ( 34) by N , we get Solving Eq. ( 33) and Eq. ( 35), we finally get The effective velocity with which the quasiparticles with v = 1 move at late times can be computed easily.The total density at late times in the pulse region is ρ = g 1 + ρ b = 3ρ0 2+aρ0 .The velocity field u in the pulse region at late times is given by u = g1 g1+ρ b = 2 3 .The effective velocity is thus: In our MD simulations in Figs.(1,2), we have taken ρ 0 = 4/5, a = 1.Plugging these values into the expressions above, we get ρb = 2 3 , ρ b = 2 7 , g 1 = 4 7 and v eff = 3.We have verified that these values match with our MD results at long times in Fig. 2. Note that v eff is the late-time speed of quasiparticles with bare velocity v = 1.

Initial condition B
In this case there is a special rod at the origin (middle of the box) with a fixed velocity v 0 = 1 and the two halves of the box on either side of the special particle are initially filled uniformly by hard rods.
The velocities of all rods, except the special one, are distributed according to the Maxwell distribution h(v) given in Eq. ( 3).This initial condition was studied by Lebowitz, Percus and Sykes (LPS) in [15].
The initial single particle phase space density is f (x, v, t = 0) = δ(v − v 0 )δ(x) + ρ 0 h(v).Since the initial distribution of the background rods are already in equilibrium, it does not change with time.However, the phase space distribution of the special rod (of velocity v 0 = 1) will change with time.At the Euler level the special rod moves ballistically with an effective velocity v eff = v0 1−aρ0 .Hence the Euler solution (for an infinite box) is given by f . However, by obtaining the exact microscopic solution of the problem in the thermodynamic limit and by performing an ensemble average, LPS showed that f (x, v 0 , t) spreads diffusively, along with a drift with velocity v eff [15], i.e, at long times one has the form f (x, v 0 , t) = δ(v − v 0 )δρ(x, t).They also obtained an explicit expression of the diffusion constant.The results of LPS were used in [18] to compute the current-current correlation and thus the Navier-Stokes (NS) term using the Green-Kubo formula.In Fig. 5a we present simulation results for δρ(x, t) which displays the spreading predicted by LPS.We observe that the spreading of the distribution increases with t and the data for different time collapse into a single function under scaling of space by √ t, as shown in Fig. 5b.This implies that the spreading grows with time as √ t [at late times it saturates in a finite box, due to the same reason for saturation in case (A)].
The origin of the growth of the width of the distribution at early times can be understood heuristically from a microscopic computation of the fluctuations of particle number as follows.Let N t be the number of particles in the interval [0, x t ], where x t is the position of the quasiparticle (special rod) with v = v 0 = 1 at time t.In the corresponding point particle picture the special particle, with velocity v = v 0 = 1, would move by a distance v 0 t in time t.Hence, the position of the rod with velocity v = v 0 = 1 is where N t , in the point particle picture, is the number of point particles that the special particle has crossed during its evolution, starting from the origin to the position v 0 t at time t.The number N t fluctuates from one realisation to another in an ensemble of initial conditions, and the fluctuation is proportional to The spread in f (x, v 0 , t) will also be proportional to the fluctuations, i.e., to √ ⟨N t ⟩.On an infinite line with uniform background of thermal particles, ⟨N t ⟩ grows linearly as t which thus leads to the √ t growth of the width in the distribution function.In a finite box, ⟨N t ⟩ cannot grow without bound, because the number of particles in the box is finite.On the hydrodynamic scale, the √ t spreading arises due to the Navier-Stokes terms in Eq. (5b) and we will now demonstrate this by obtaining a analytic solution of the Navier-Stokes equation (5b) on the infinite line.For this we make the ansatz: where h(v) is given in Eq. ( 3).This ansatz is motivated by the fact that the number of particles with a given velocity is conserved, and that the distribution of the background rods does not change with time.Plugging the ansatz into the Navier-Stokes equation, we get the following drift-diffusion equation for δρ(x, t) after ignoring the non linear terms proportional to (δρ) 2 : where The solution of this for the LPS-like initial condition is given by: which is exactly the solution that was obtained by LPS from a completely microscopic analysis [15].In Fig. 5b, we verify that the expression in Eq. ( 42) agrees with the MD simulation results.Our numerical results thus provide a direct demonstration of an observable effect of the NS terms in the hydrodynamic equations.We see that there is a good scaling collapse, and a nice agreement with the solution of NS equation.We have taken, N = 2 × 10 6 , L = 2.5 × 10 6 , a = 1.0, v 0 = 1.0 (and so µ(v 0 ) ≈ 1.0).For MD, ensemble averaging has been done over 10000 realizations.The times considered are much before the pulse hits the boundary of the box, hence the system is effectively infinite.

Euler vs MD for initial condition C
Finally we consider the free expansion set up in which the N hard rods are initially confined to the left half of the box of size L and distributed uniformly in space with density ρ 0 = N /L.The velocities of the rods are drawn from the Maxwell distribution h(v) in Eq. (3).As in the previous cases, we have hard reflecting walls at x = 0 and x = L.We now follow the same approach outlined in Sec.(4.1), to obtain a solution of the Euler equation for this initial condition, via the mapping to hard point gas.The solution in the point particle picture is similar to that obtained in [27], with the density given by: where ρ 0 = N /L and T = 1.In the a → 0 limit, the above expression of the density profile ρ 0 (x, t) matches with those obtained in [27].Using the inverse mapping in Eq. ( 16), the density profile ρ(x, t) of the rods can be found, where recall x = x ′ + aF 0 (x ′ , t) and the cumulative density profile, F 0 (x ′ , t), can be computed from ρ 0 (x ′ , t).In Fig. 6(a), we compare the theoretically computed profiles of the rods at different times with the density profiles obtained from MD simulation, and we observe excellent agreement.From this plot, we observe that with increasing time the density profile of the rods spreads to the right half of the box in a monotonic fashion and finally approaches a time-independent spatially uniform profile which is consistent with a GGE state.
In a similar way, the exact Euler expression for the momentum density field p(x, t) can be obtained.
First we compute the momentum field p 0 (x ′ , t) in the point particle picture and then transform to the momentum field for the hard rods using p(x, t) = p 0 (x ′ , t) We find where T = 1.In Fig. 6(b) we compare this with the results obtained from the MD simulations and we again see very good agreement.Here We observe that initially the momentum profile was zero everywhere.Once where T = 1 and 2ρ0 1−2aρ0 is the initial density for z ′ ∈ (−∞, 0).The early time plots (for t < 300) in Fig. ( 6) can be obtained by transforming the above simpler functions to ρ(x, t) and p(x, t) using transformations in Eq. ( 16), (44) and Eq.(18).The distortions of the densities of the point particles are appearing due to the non-linear transformations.We see that the blue lines lie at the edges of the red region and are curved in contrast to the straight lines for the point particle [see [27]].At late times, the distribution wraps around the allowed region multiple times and thus creates fine structures.Here, N = 10000, L = 4, a = 0.0001 Domain line for initial condition C: For the point particle case, f 0 (x ′ , v, t) has a discontinuity in x space (for a given v) for the free expansion problem.Since there is a mapping between the point particle Euler equation and the hard rod Euler equation, we expect that the Euler equation for hard rods will admit a similar discontinuity.We call the line of discontinuity of f (x, v, t) in the single-particle phase space as the "domain line".For the free expansion problem, the domain line can be found implicitly in the following manner.For times before the particles hit the right end of the container, the domain line for the point particle problem is given by x ′ = vt + L 2 − N a.For general times (including times after the particles hit the right end of the container), we can do an analysis similar to [27] to show that the single particle phase space distribution for the point particle problem for general times is given by: where n is an integer.From this, the domain line for the point particle problem can be computed as zeros of the argument of the theta functions appearing in the equation above.We know how x ′ maps to x (x = x ′ + aF 0 (x ′ , t)) within the framework of the Euler equation.Thus we can compute the domain line for the hard rod problem as predicted by the Euler equation.We computed the domain line for the hard rod problem as predicted by the Euler equation, and plotted it (blue line) along with the phase space plot for the hard rods (red dots) predicted by the MD simulation.The plots are shown in Fig. 7.We see that the blue line lies at the edge of the region occupied by the red dots.Thus the domain line predicted by Euler equation agrees with that predicted by MD simulation.We observe some key differences between hard rods (interacting integrable), alternate mass point-particle gas (non-integrable) and equal mass point particle gas (non-interacting integrable).In the hard rod gas, we see a sharp domain line which is not a straight line at early times (Fig. 7).In equal mass point particle gas discussed in [27], a sharp and straight domain line was observed.However, in the alternate mass point particle gas, no sharp domain line was observed [28].

Conclusion
In this paper we studied the macroscopic evolution of a collection of hard-rods in one dimension starting from three different initial conditions: (A) A uniformly filled box with an inhomogeneous velocity distribution -half of the box is in thermal equilibrium and the other half has particles with a fixed velocity v = 1, (B) One special particle with fixed velocity v 0 at the origin in the presence of a spatially uniform background of other particles having thermal velocity distribution and (C) Free expansion from half of the box filled uniformly with thermal velocity distribution.For initial conditions (A) and (C) we find that the molecular dynamics results agree very well with the solutions of the Euler equations.However, for (A) we observe shocks at all times and find discrepancies from the Euler solutions at the location of the shocks which can be attributed to the Navier-Stokes corrections to the Euler equations.For initial condition (B), the effect of the Navier-Stokes terms is more dramatic and here we show that the effect can be understood from the analytic solution of the Navier-Stokes equation.
Our second important finding is the absence of GGE for initial conditions (A) and (B), whereas for initial condition (C) the system at late time approaches to a GGE state.The absence of GGE in the initial conditions (A) and (B) are manifested by the fact the density profile remains time dependent at all times.
We find that the effect of the Navier-Stokes terms is very weak and to observe its effect one requires to have singular velocity distributions in the initial conditions such that a shock in the density profile survives for macroscopic time scale.At the location of the shock the large density gradient makes the contribution from the Navier-Stokes terms significant and consequently the solution near the shock (in Euler solution) becomes different from the Euler solutions.This is also what is observed in non-integrable systems [30,31,32].
Since the effect of dissipation is most noticeable near a shock, it is worth asking the question that for which initial conditions are shocks formed.It is easy to see there will be a shock only if the mapped point particle problem has a shock.This can be seen in the following way.Let δx 0 be the length scale over which the density is varying in the point particle problem.Then, using Percus' microscopic mapping, δx = δx 0 + aδN , where δx is the corresponding length scale in the hard rod problem, and δN is the number of point particles in the length scale δx 0 .If the point particle problem has a shock, then δx 0 will be small and δN ∼ O(1).Thus δx will be of the order of few rod lengths, and there will be a shock in the hard rod problem also.If there is no shock in the point particle problem, then both δx 0 and δN will be large, and hence δx will also be large.Thus there will not then be any shock in the hard rod problem.We find that shocks in density profile of the point particles ( and hence of the hard rods) persists with time if they start with singular velocity distributions [such as initial conditions (A) and (B)].On the other hand if the rods have smooth velocity distribution to start with, then even if there are discontinuities in the density profile initially, the profiles at later times becomes smooth [such as initial condition (C)].
For initial condition A, one may be curious what will happen if one chooses a Maxwellian distribution centered at v = 1 instead of a δ-function distribution [h(v) = δ(v − 1)] that we have choosen for the left half of the particles.If one does that, then the corresponding point particle problem will not have any shocks as both the halves will perform free expansion independently and we do not observe shocks in free expansion.If the spread of the Maxwellian is small enough, then one may observe a weak shock at small times, however the shock will weaken with time and eventually give rise to a smooth density profile.The time scale over which the shock will weaken will be inversely proportional to the spread of the Maxwellian, hence the shock for our initial condition A is infinitely long lived on the Euler level.

Fig. 1 :
Fig. 1: Comparing solution of Euler equation with MD simulation for initial condition A: Plot comparing the solution of the Euler equation with those of molecular dynamics for (a) the density of v = 1 particles, denoted by g(x) and (b) the density of background particles, denoted by ρ b (x).Dashed lines are MD simulations and solid lines are solutions of Euler equation.We have taken times t = 0 (dark blue), t = 40 (orange), t = 80 (green), t = 160 (red), t = 240 (violet), t = 320 (brown), t = 400 (pink), t = 480 (grey), t = 560 (mud green) and t = 640 (cyan).We see that even for very long time like t = 640, the profile obtained from MD does not relax to GGE, i.e., it does not become uniform.We also see that there is a discrepancy between MD and Euler solutions at the shock front due to dissipative effects.We have taken length of the box L = 2500, total number of particles N = 2000 and length of rod a = 1.0.We have performed ensemble averaging over 5000 realizations while doing MD.
(a) where profiles for different times collapse under the scaling of x by √ t.Microscopically the spreading originates from the fluctuations in the number of the background rods (having Maxwell velocity distribution)

Fig. 3 :
Fig. 3: (a) Time dependence of the width and (b) the fluctuation in the location of the shock: (a) This shows the structure of the shock for the initial condition A at different times for a given system size (N = 2000, L = 2500).The curves have been shifted so that the shock fronts for all the curves coincide.While doing MD, ensemble averaging over 5000 initial microstates were performed.We see that there is trend of increasing width with time while in the inset (which shows curves for t = 20, 40, 80), we see that there is a scaling collapse in the variable x/ √ t for short times when the v = 1 pulse does not know about the boundaries of the system and hence behaves like it is in an infinite system.We explain this √ t dependence in Appendix A. (b) This shows plot of g(x, t) for two different realizations (microstates) for initial condition A. We see that in a single realization the shock remains sharp, while the positions of the shock front in two realizations are different.Consequently, ensemble averaging will lead to smearing of the shock and thus is necessary to observe dissipation.Here we have chosen N = 8000, L = 10000.

Fig. 5 :
Fig. 5: Verifying NS equation for LPS-like initial condition: (a) This figure compares the results from MD simulation for the evolution of the density profile with those obtained from the solution of the NS equation, for the LPS-like initial condition.(b) We show the plot in terms of the scaling variables.We see that there is a good scaling collapse, and a nice agreement with the solution of NS equation.We have taken, N = 2 × 10 6 , L = 2.5 × 10 6 , a = 1.0, v 0 = 1.0 (and so µ(v 0 ) ≈ 1.0).For MD, ensemble averaging has been done over 10000 realizations.The times considered are much before the pulse hits the boundary of the box, hence the system is effectively infinite.

Fig. 6 :
Fig. 6: Comparing solution of Euler equation with MD simulation for initial condition C: Plots of the density and momentum profiles and comparison of the exact solution of Euler equation (solid lines) and the profiles obtained from MD (dashed lines) for the free expansion problem.We have shown for times t = 0 (dark blue), t = 10 (orange), t = 20 (green), t = 40 (red), t = 100 (violet), t = 150 (brown), 200 (pink), t = 300 (grey), t = 500 (muddy) and t = 1000 (light blue).We have taken N = 1000, L = 2500 and averaged over 100 realizations.
the gas is released, the rods with positive velocities near the middle of the box start moving to the right half.Thus the gas creates a positive momentum profile near the centre of the box.As time progresses, more particles move to the right half and consequently the momentum profile spreads on both halves of the box.After some time, of the order L/ √ T , finite size effects start showing, and some rods get reflected from the right walls.As a result, the motion of these rods start reducing the momentum field.At very late times, each rod has undergone several collisions with both the walls and the gas equilibrates.The time scale of equilibration is also of the order L/ √ T .At this stage, one has rods of opposite velocities with equal probabilities at any point of the box which leads again to a zero momentum profile everywhere.Note that shifting the origin to L/2 − N a (i.e., x ′ → z ′ = x ′ − L/2 + N a) in the point particle problem, and taking L, N → ∞ with N /L = ρ 0 , one obtains the solution of Euler GHD on the infinite line for times t << L √ T as:

Fig. 7 :
Fig. 7: Evolution of the domain lines for initial condition C: Plot of the phase space distribution of the rods (red dots) for the free expansion problem at times (a) t = 0.05 (b) t = 0.5, (c) t = 0.7 and (d) t = 4. Solid blue lines represent the domain lines obtained from the exact solution of the Euler equation.We see that the blue lines lie at the edges of the red region and are curved in contrast to the straight lines for the point particle [see[27]].At late times, the distribution wraps around the allowed region multiple times and thus creates fine structures.Here, N = 10000, L = 4, a = 0.0001