Blast in a One-Dimensional Cold Gas: From Newtonian Dynamics to Hydrodynamics
Subhadip Chakraborti, Santhosh Ganapa, P. L. Krapivsky, Abhishek Dhar
BBlast in the one-dimensional cold gas: From Newton to Euler and Navier-Stokes
Subhadip Chakraborti , Santhosh Ganapa , P. L. Krapivsky and Abhishek Dhar International Centre for Theoretical Sciences, Tata Institute of Fundamental Research, Bengaluru 560089, India and Department of Physics, Boston University, Boston, MA 02215, USA (Dated: February 17, 2021)It is known that a gas composed of a large number of atoms following Newtonian dynamics can bedescribed by the continuum laws of hydrodynamics. Proving this rigorously is one of the outstandingopen problems in physics and mathematics. Surprisingly, precise numerical demonstrations of theequivalence of the hydrodynamic and microscopic descriptions are rare. We test this equivalencein the context of the classic problem of the evolution of a blast-wave, a problem that is expectedto be at the limits where hydrodynamics would work. We study a one-dimensional gas for whichthe hydrodynamic Euler equations for the conserved fields of density, momentum and energy areknown to have self-similar scaling solutions. Our microscopic model consists of hard point particleswith alternate masses, which is a non-integrable system with strong mixing dynamics. Our extensivemicroscopic simulations find a remarkable agreement with hydrodynamics, with deviations in a smallcore region that are understood as arising due to heat conduction.
The Navier-Stokes-Fourier (NSF) equations [1] de-scribing the evolution of the density, velocity and temper-ature apply to an enormous range of phenomena, e.g., toatmospheric flows. At the same time we know that thereis a different description for air where we view it as acollection of molecules which follow Newton’s equationsof motion. How accurate is the hydrodynamic descrip-tion and what are the limits of its applicabilty? Herewe address this question, which is of fundamental as wellas practical importance. A rigorous derivation of thecontinuum hydrodynamics description from the atom-istic one is an open problem [2–4]. Much progress hasbeen made towards a partial solution where one startsfrom the Boltzmann kinetic theory equation for the singleparticle distribution function and derives the equationsof hydrodynamics as a systematic expansion in a smallparameter [5–7]. This still leaves one with the problemof deriving the Boltzmann equation from Newton’s equa-tion which has been achieved only for the case of a dilutegas. For the case where one adds a weak noise to theNewtonian dynamics (still satisfying the same conserva-tion laws) a rigorous derivation of the Euler equations forthe hydrodynamic fields has been achieved [8].Somewhat surprisingly, there appears to be no directnumerical verification of the fact that the equations ofhydrodynamics accurately reproduce the evolution of theconserved fields following from the microscopic dynamics.The present work provides such a detailed comparison inthe context of the classic blast-wave problem for whicha self-similar scaling solution of the Euler equations wasobtained more than sixty years back, independently byTaylor [9, 10], von Neumann [11] and Sedov [12, 13], andis referred to as the TvNS solution. The evolution of ablast wave emanating from an intense explosion was firststudied to understand the mechanical effect of bombs.The rapid release of large amount of energy in a local-ized region produces a surface of discontinuity beyondwhich the quantities concerned like density, velocity andtemperature fields change discontinuously [14, 15]. Theblast-wave problem thus presents an extreme case to testthe validity of the hydrodynamic description. From the point of view of microscopic models, the hard sphere sys-tem would be the natural candidate since much is knownanalytically (equation of state in terms of virial expan-sions and kinetic theory via the Chapman-Enskog expan-sion). Secondly, this system can be simulated very effi-ciently using event driven simulations. However, largescale molecular dynamics (MD) simulations [16–23] haveso far not found clear and complete agreement with theTvNS solution from hydrodynamics. It was suggestedthat possible reasons for the differences could be the lackof local equilibration or due to the contribution of viscos-ity and heat conduction in the hydrodynamic equations(not included in the TvNS analysis).In this Letter, we address this question by studyinghard point particles with binary mass distribution —in one dimension, particles with equal masses just ex-change velocities, so mass dispersion is necessary for re-laxation, and the binary mass distribution is thus thesimplest setting where relaxation is possible. Further-more, we assume that adjacent particles have differentmasses (say m and m ). This alternating hard particle(AHP) gas has been extensively investigated in the con-text of the breakdown of Fourier’s law of heat conductionin 1D [24–32]. The hard particle (and hard rod) systemhas been investigated earlier in the context of breakdownof the hydrodynamic description in 1D [33, 34] and morerecently it has also been used to study the so-called Rie-mann problem [35] and thermalization [36]. Compared tothe hard sphere system in higher dimensions, the 1D gashas several advantages — the equation of state is knownexactly and is that of an ideal gas and simulations arefaster since collisions occur only with nearest neighbors.Note that while the equilibrium physics is that of an idealgas, the dynamics is non-integrable and known to havegood ergodic properties (while being non-chaotic) [37].We provide here a detailed study of the evolution of theblast wave initial condition in a 1D gas with an ideal gasequation of state. From extensive molecular dynamicssimulations of the AHP gas we compute the evolutionof the profiles of density, velocity and temperature fieldsand thereby extract the scaling forms obtained in the a r X i v : . [ c ond - m a t . s t a t - m ec h ] F e b FIG. 1.
Heat maps showing the spatio-temporal evolution of the density, velocity and temperature fields, starting from ini-tial conditions corresponding to a Gaussian initial temperature profile and ρ ( x,
0) = ρ ∞ = 1 . , v ( x,
0) = 0 . The simulationparameters were N = L = 24000 , E = 32 , µ = 1 . and an ensemble average was taken over initial conditions. long time limit. We make comparisons with the scalingsolution predicted from the TvNS solution of the Eulerequations for which we provide an exact solution. Wefind that a complete explanation of the simulation of theblast requires us to go beyond the Euler equation andinclude the effect of heat conduction. We thus discussthe NSF equations, for which we present results from anumerical solution as well as a scaling analysis. The Taylor-von Neumann-Sedov solution —
We con-sider a 1D ideal gas at zero temperature, ambient density ρ ∞ and we inject an amount of energy E in a localizedregion of extent σ . The hydrodynamic Euler equationsfor the density field, ρ ( r, t ), velocity field, v ( r, t ), and thetemperature field, T ( r, t ) are given by ∂ t ρ + ∂ x ( ρv ) = 0 , (1a) ∂ t ( ρv ) + ∂ x ( ρv + P ) = 0 , (1b) ∂ t ( ρe ) + ∂ x ( ρve + P v ) = 0 , (1c)where for an ideal gas we have e = v / T / (2 µ ), P = ρk B T /µ while µ = ( m + m ) / k B = 1. Fromdimensional analysis one can show [14] that at long timesthe solution of these equations develop discontinuities ata distance R ( t ) = (cid:18) Et Aρ ∞ (cid:19) , (2)where A is a dimensionless constant factor which, quiteremarkably, can be determined exactly [38]. From dimen-sional analysis one can further show [14] that the fieldswill take the following self-similar scaling form ρ ( x, t ) = ρ ∞ G ( ξ ) , (3a) v ( x, t ) = 23 xt V ( ξ ) = 2 α t / ξV ( ξ ) , (3b) T ( x, t ) = 4 µ x t Z ( ξ ) = 4 µα t / ξ Z ( ξ ) , (3c)where ξ = x/R ( t ) is the scaled position, α =[ E/ ( Aρ ∞ )] / and the scaling functions G, V, Z need tobe determined. The factors 2 / , /
27 are inserted for convenience, e.g. from Eq. (2) one finds that the velocityof the shock wave is (2 / R/t and this suggests the use ofthe factor 2 /
3. Plugging these into Eqs. (8a-8c) we findthat the scaling functions satisfy a set of coupled firstorder ordinary differential equations in the variable ξ .Using the condition of conservation of energy and the so-called Rankine-Hugoniot conditions, specifying the fielddiscontinuities at the shock front, allows one to obtaina complete closed-form solution of the problem, i.e. thefunctions G, V, Z and the constant A [39]. In [38] and[40] we describe some details of the solution. Microscopic dynamics and initial conditions —
Oursystem is a 1D gas of N hard point particles movinginside a box ( − L/ , L/ m andall odd ones m . In this case the only known conservedquantities are particle number, total momentum and en-ergy and we expect a hydrodynamic description in termsof the corresponding conserved fields. The hydrodynamicfields are obtained from the microscopic variables usingthe following standard relations: ρ ( x, t ) = N (cid:88) j =1 m j (cid:104) δ ( q j ( t ) − x ) (cid:105) , (4a) p ( x, t ) = ρ ( x, t ) v ( x, t ) = N (cid:88) j =1 m j (cid:104) v j δ ( q j ( t ) − x ) (cid:105) , (4b) E ( x, t ) = ρ ( x, t ) e ( x, t ) = N (cid:88) j =1 m j (cid:10) v j δ ( q j ( t ) − x ) (cid:11) . (4c)Here (cid:104) ... (cid:105) indicates an average over an initial distribu-tion of microstates that correspond to the same initial (a) (d)(b) (e)(c) (f ) FIG. 2.
Results from Molecular DynamicsSimulations : ( a, b, c ) Late time evolution of density, ve-locity and temperature fields, starting from the initialconditions corresponding to a Gaussian initial temperatureprofile and ρ ( x,
0) = ρ ∞ = 1 . , v ( x,
0) = 0 . The simulationparameters were N = L = 24000 , E = 32 , µ = 1 . and anaverage was taken over initial conditions. ( d, e, f ) Thisshows the x/t / scaling of the data. We observe a very goodcollapse of the data at the longest times and an agreementwith the TvNS scaling solution (dashed line). macrostate. For a non-integrable system, it is expectedthat the evolving system is in local thermal equilibriumand the three fields should contain the local thermody-namic information at any space-time point x, t . Usingthermodynamics we obtain the other local fields in termsof the basic fields ρ, v, e . Thus the internal energy perunit mass is given by (cid:15) ( x, t ) = e − v / T ( x, t ) = 2 (cid:15) ( x, t ) , P ( x, t ) =2 ρ ( x, t ) (cid:15) ( x, t ).We consider an initial macrostate where the gas has afinite uniform density ρ ∞ , zero flow velocity v , and is atzero temperature everywhere except in a region of width σ centred at x = 0. This is the region of the blast and wetake a smooth Gaussian profile E ( x,
0) = E √ πσ e − x / σ .We have verified that our long-time results are indepen-dent of the initial profile as long as it is localized and thetotal energy E is fixed [40]. The procedure to realize thismacrostate in the microscopic simulations of the AHPgas and other details of our numerics are given in [38]. Comparison of simulations with the TvNS solution —
We now present the results of the microscopic simula-tions for the evolution of the hydrodynamic fields ρ ( x, t ), v ( x, t ), and T ( x, t ) starting from the blast-wave initial conditions. In Fig. (6) we show the spatio-temporalevolution of the three fields (for individual particle tra-jectories see [38]). We see a sharp shock-front thatevolves sub-ballistically. The mass density and flow ve-locity of the gas are peaked around the blast front whilethe temperature has an additional peak at the centre.In Fig. (2a,b,c) we plot the evolution of the blast waveat long times. In Figs. (2d,e,f) we show the scaledfields (cid:101) ρ = ρ ∞ G − ρ ∞ , (cid:101) v = t / v ( x, t ) = (2 α/ ξV and (cid:101) T = t / T ( x, t ) = (4 µα / ξ Z as functions of the scal-ing variable ξ . We find an excellent collapse of the dataeverywhere except near the blast centre. In the regionwhere there is a collapse of data, we find a perfect agree-ment with the exact TvNS scaling solution (plotted asblack dashed lines). Close to the origin, the TvNS so-lution predicts the singular forms [40] G ( ξ ) ∼ ξ / and Z ( ξ ) ∼ ξ − / which implies that near the origin the den-sity vanishes as ρ ∼ x / and the temperature diverges as T ∼ x − / . These are unphysical and clearly the simula-tions do not observe this. The reason for this deviationis that the effect of dissipation terms, specifically heatconduction, become important in the core region and weneed to use the NSF equations. We now discuss this. Comparison of simulations with numerical solution ofthe NSF equations. —
We now compare our simulationresults and the TvNS solution with those from the fulldissipative hydrodynamic equations. We re-write the Eu-ler part from Eqs. (8a-8c) in a slightly different form and,after including dissipation, get the NSF equations: ∂ t ρ + ∂ x ( ρv ) = 0 , (5a) ρ ( ∂ t + v∂ x ) v + ∂ x ( ρT /µ ) = ∂ x ( ζ∂ x v ) , (5b) ρ µ ( ∂ t + v∂ x ) (cid:18) Tρ (cid:19) = ∂ x ( ζv∂ x v + κ∂ x T ) , (5c)where ζ denotes the bulk viscosity and κ is the thermalconductivity of the system. These transport coefficientscan depend on the fields and, based on the Green-Kuborelations, we expect their temperature dependence to be ζ ∼ T / and κ ∼ T / . A recent numerical study [30]suggests the density dependence κ ∼ ρ / . In our numer-ical study we have thus used the forms ζ = D T / and κ = D ρ / T / , where D and D are constants.We solve Eqs. (5a-5c) mumerically [38] for the sameinitial conditions as considered in the microscopic simula-tions, namely ρ ( x,
0) = ρ ∞ , v ( x,
0) = 0 and E ( x,
0) givenby the Gaussian form with total energy E . The numer-ical results are plotted in Figs. (refnoscaleHydnewa,b,c)which shows the profiles of the fields, ρ ( x, t ), v ( x, t ) and T ( x, t ) at different times. In Figs. (3d,e,f) we plot thescaled fields (cid:101) ρ, (cid:101) v, (cid:101) T and verify the expected agreementwith the TvNS solution everywhere except in the core. Inthe core region the NSF numerical solution does not haveany singularities, unlike the TvNS solution. In Fig. (4)we plot together the long time microscopic simulationresults and the NSF results and it can be seen that inthe core region, the NSF solution is in better qualitativeagreement with simulations as compared to the TvNSsolution. (a) (d)(b) (e)(c) (f ) FIG. 3.
Results from solution of the Navier-Stokes-Fourier equations : ( a, b, c ) Evolution of the hydrodynamicfields ρ ( x, t ) , v ( x, t ) , T ( x, t ) as obtained from a numerical so-lution of the NSF equations (5a) - (5c) , starting from the sameinitial conditions as used in the simulations for Fig. (2) . Theother parameters used in the numerics are D = 1 , D = 1 and L = 4000 . ( d, e, f ) Scaling plot and comparison with theTvNS solution. (a) (b) FIG. 4.
Comparison of the scaled fields (cid:101) ρ, (cid:101) T from micro-scopic simulations (at time t = 80000 ), from Navier-Stokes-Fourier equations (at t = 6400 ) and from the TvNS scaling solution . We proceed to estimate the size of the core and thecorresponding scaling solution. Using Eq. (5c) we seethat the heat conduction term becomes important at alength scale X such that ρ Tt ∼ T / X ρ / , where we ne-glect constant factors. Using the ξ → ρ ∼ √ ξ ∼ | x | / /t / and T ∼ ξ − / ∼ t / / | x | / , leadsto the estimate X ∼ t (6)for the size X of the hot core where heat conduction is (a) (b) FIG. 5.
Verification of the scaling form Eqs. (7) , for ρ and T , in the core of the blast. The comparison with data bothfrom simulations and from the NSF solution are shown. Thedata are the same as in Figs. (2,3). The value η = 10 x/t / corresponds to ξ ≈ . . important. The outer core solution (TvNS) and the innercore solution are comparable at | x | = X and this allowsus to determine the inner core scaling laws. Thus, thetemperature at the centre of the explosion is estimatedfrom T ∼ X − t − and Eq. (6). Similarly, the densityat the centre of the explosion can be estimated from ρ ∼ X / /t / and Eq. (6). We then arrive at the asymptoticdecay laws ρ ∼ t − , T ∼ t − while the velocityin the hot core scales as X/t ∼ t − / .The hydrodynamic fields in the hot core where dis-sipative effects play the dominant role should thereforeexhibit scaling behaviors in terms of the scaled spatialcoordinate η = x/X . For large x, t , in a region with η ∼ O (1), we expect the self-similar forms ρ = t − (cid:101) G ( η ) , v = t − (cid:101) V ( η ) , T = t − (cid:101) Z ( η ) . (7)In Fig. (5) we show that the data from the numericalsolution of the NSF equations approximately satisfy thisscaling form, though it appears that the convergence issomewhat slow. For comparison we also plot the simula-tion data at the last two times under the same scaling.A more detailed discussion of the inner-core scaling solu-tion and its “matching” with the outer-core solution canbe found in Ref. [40]. Conclusions. —
A detailed comparison was made ofthe predictions of hydrodynamics for a compressible one-dimensional gas with results from microscopic dynamicsof a many body interacting gas, for the physical exam-ple of energy expanding into a cold environment. Asour main results, the front position R ( t ) ∼ t / and thescaling functions were obtained in closed form (TvNSsolution) and a remarkable agreement was observed be-tween these and microscopic simulations. Deviationswere seen in a core region whose width vanishes if mea-sured in units of the TvNS scaling variable ξ ∼ x/t / .Thus we numerically establish that the Euler equationscompletely describing the hydrodynamic behavior of theAHP gas. Earlier attempts in two- and three-dimensionalsystems [18–23] have not been able to obtain the unam-biguous agreement as reported here— this is because the1D hard point gas is an ideal model for this study. It’sequilibrium properties including the equation of state areknown exactly which allows us to obtain the TvNS scal-ing functions exactly. In contrast, in higher dimensionsone relies on a virial expansion to get the equation ofstate and then the TvNS solution has to be numericallyfound. Secondly our system is non-integrable and knownto have good ergodicity properties and can be simulatedvery efficiently for large systems and long times.The deviation at the core was understood as arisingfrom the contribution of thermal conduction terms inthe energy conservation equation. This then leads us togo beyond Euler and examine carefully the full Navier-Stokes-Fourier equation for our 1D gas. For an accuratecomparison with the simulation results, we need the pre-cise form of the thermal conductivity, κ , and in particularits dependence on temperature and density. For the AHPgas it is known that κ ∼ T / and earlier studies suggest κ ∼ ρ / . Our analysis of the resulting NSF equationsgives us the estimate | X | ∼ t / for the core size as wellas scaling forms for the fields which we verified in thenumerical solution of the NSF equations. The inclusionof the dissipative terms leads to a qualitative agreementbetween the microscopic simulation results and hydrody-namics even in the core region, in particular it cures oneof the chief problem of the TvNS solution — the diver-gence of the temperature at the centre of explosion. Anopen and challenging problem that remains is to obtainquantitative agreement between simulations and hydro-dynamics at the core centre which would however requiremore information on the precise form of the thermal con-ductivity of our 1D gas. Although we focused on 1Dwhere we can provide detailed numerical comparison, ourapproach can be extended to higher dimensions where wefind a similar growing core region where thermal conduc-tivity plays a role. The application of hydrodynamics toquantum systems is also of much recent interest [41, 42]and exploring this would be another interesting direction.We thank Anupam Kundu, R. Rajesh, Sriram Ra-maswamy, Samriddhi Sankar Ray and Vishal Vasan foruseful discussions. We acknowledge support of the De-partment of Atomic Energy, Government of India, underproject no. RTI4001. Supplemental material for ‘Blast in theone-dimensional cold gas: From Newton to Eulerand Navier-Stokes’A. Exact TvNS scaling solution of Euler equations
We start with the equations ∂ t ρ + ∂ x ( ρv ) = 0 , (8a) ∂ t ( ρv ) + ∂ x ( ρv + P ) = 0 , (8b) ∂ t ( ρe ) + ∂ x ( ρve + P v ) = 0 , (8c)where we recall that e = v / T / (2 µ ) and P = ρT /µ and look for a self-similar scaling solution of these equa-tions with a shock at the location R ( t ). At the shockthe condition of weak solutions lead to the the followingRankine-Hugoniot boundary conditions ρ ( R ) v ( R ) ρ ( R ) − ρ ∞ = U, (9) ρ ( R ) v ( R ) + P ( R ) ρ ( R ) v ( R ) = U, (10) ρ ( R ) v ( R ) e ( R ) + P ( R ) v ( R ) ρ ( R ) e ( R ) = U, (11)where U = ˙ R . In addition, we need to use the conditionof energy conservation: (cid:90) R dxρ ( v + T /µ ) = E. (12)We look for scaling solutions of the following form ρ = ρ ∞ G ( ξ ) , v = 2 x t V ( ξ ) , T = 4 µ x t Z ( ξ ) . (13)where the fields now depend on the single dimensionlessvariable ξ = xR , where R ( t ) = (cid:18) Et Aρ ∞ (cid:19) / . (14)Using the scaling form and after some simplification, theRankine-Hugoniot boundary conditions lead to G (1) = 2 , V (1) = , Z (1) = . (15)Plugging the scaling forms (13)–(14) into Eq. (8a) andEq. (8b) and, defining (cid:96) = ln ξ , we obtain dVd(cid:96) + ( V − d ln Gd(cid:96) = − V, (16) d ln Zd(cid:96) − d ln Gd(cid:96) = 3 − VV − . (17)Using the fact that energy is conserved in the region( − ξR ( t ) < x < ξR ( t )), one can obtain the following re-lation between Z and V as [14, 40] Z = γ ( γ − − V ) V γV − . (18)The total energy conservation condition in Eq. (12) gives E = ρ ∞ R t (cid:90) dξ ξ V V − G, (19)where we have used (18). Combining (19) with (14) weobtain A = 49 (cid:90) dξ ξ V V − G. (20)As shown in [40] it is possible to obtain an explicit so-lution of the set of equations (16,17,18) with the bound-ary conditions in Eq. (15) and the value of A given byEq. (20).We find [40] A = 152 / V is obtained bysolving the implicit equation ξ = 2 − (3 V − V − (3 − V ) − . (21) Z is then given by Eq. (18), and finally, G = 2 (1 − V ) (3 V − (3 − V ) − . (22)Thus we have a completely closed form expression for theTvNS solution for the 1D ideal gas. B. Specification of microscopic dynamics andinitial conditions
Let q j , p j , and m j be respectively the position, momen-tum, and mass of the j th particle with j = 1 , , . . . , N .In between collisions the particles move ballistically withconstant speeds. On collision between particles i and i + 1, their post-collisional velocities follow from conser-vation of energy and momentum and are given by v (cid:48) i = [2 m i +1 v i +1 + v i ( m i − m i +1 )]( m i + m i +1 ) (23) v (cid:48) i +1 = [2 m i v i + v i +1 ( m i +1 − m i )]( m i + m i +1 ) . (24)In our numerical simulations we consider the followingGaussian profile for the initial energy: E ( x,
0) = E √ πσ e − x / σ . (25)The total energy of the blast is E .To obtain this in the particle model, we first dis-tribute N ordered particles numbered say as i = − N/ , − N/ , . . . , N/ x = − L/ L/
2, so that the number density is n ∞ = N/L = ρ ∞ /µ ,where µ = ( m + m ) / ρ ∞ the am-bient mass density of the gas. For the N c centre particleswith i = − N c / , − N c / , . . . N c /
2, we choose theirvelocities from the Maxwell distribution,
P rob ( v i ) = (cid:112) m i / (2 πT i ) e − m i v i / (2 T i ) , where T i = 2 E/N c . Set thevelocities of all other particles to zero. The size of theinitial blast is thus approximately s = N c /n ∞ . t x FIG. 6.
This plot shows the trajectories of N = 1000 particles,starting from initial conditions corresponding to a Gaussianinitial temperature profile and ρ ( x,
0) = ρ ∞ = 1 . , v ( x,
0) = 0 .The expected shock front position R ( t ) is shown as a dashedline. The stationary particles are not shown. We note that E ( x,
0) = (cid:104) (cid:80) Ni =1 m i v i (0) / δ ( x − q i (0)) (cid:105) .Since, for large N , each particle’s position is a Gaussianwith mean ¯ q i = i/n ∞ and variance σ = L/ (4 n ∞ ), wethen get: E ( x,
0) = N c / (cid:88) i = − N c / T i √ πσ e − ( x − ¯ q i ) / (2 σ ) ≈ n ∞ T i (cid:90) s/ − s/ dy e − ( x − y ) / (2 σ ) √ πσ = n ∞ T i (cid:20) erf (cid:18) s − x √ σ (cid:19) + erf (cid:18) s + 2 x √ σ (cid:19)(cid:21) . (26)When s (cid:28) σ this simplifies to E ( x, ≈ E e − x / (2 σ ) √ πσ . (27) C. Details of numerical techniques
Molecular dynamics simulations:
In all our sim-ulations we took m = 1 , m = 2 (so that µ = ( m + m ) / . ρ ∞ = 1 . E = 32 , N c = 32. In ourlargest simulations we took N = 24000, L = 24000 andaveraged over an ensemble of R = 10 initial conditions.For each microscopic initial condition, we evolved thesystem with the Hamiltonian dynamics. The moleculardynamics for this system was done using an event-drivenalgorithm which updates the system between successivecollisions. Solution of Euler equations with Navier-Stokesterms : The numerical solution used the MacCormackmethod [43], which is second order in both space andtime and used a discretization dx = 0 . dt = 0 . D = D = 1 and the system size was taken as L = 4000. We have evolved up to time t = 6400 whichis before the energy reaches the boundary. [1] R. K. Zeytounian. Navier-Stokes-Fourier equations: arational asymptotic modelling point of view . Springer Sci-ence & Business Media, 2012.[2] R. Esposito and M. Pulvirenti. From particles to flu-ids. In S. Friedlander and D. Serre, editors,
Handbook ofmathematical fluid dynamics , volume 3, chapter 1, pages1–82. Elsevier North-Holland, Amsterdam, 2004.[3] I. Gallagher. From newton to navier–stokes, or how toconnect fluid mechanics equations from microscopic tomacroscopic scales.
Bull. Amer. Math. Soc. , 56(1):65–85, 2019.[4] A. N. Gorban. Hilbert’s sixth problem: the endless roadto rigour.
Phil. Trans. R. Soc. A. , 376(2118):20170238,2018.[5] S. Chapman, T. G. Cowling, and D. Burnett.
The math-ematical theory of non-uniform gases: an account of thekinetic theory of viscosity, thermal conduction and diffu-sion in gases . Cambridge university press, 1990.[6] A. de Masi, R. Esposito, and J. L. Lebowitz. Incom-pressible navier-stokes and euler limits of the boltzmannequation.
Commun. Pure Appl. Math. , 42(8):1189–1214,1989.[7] L. Saint-Raymond.
Hydrodynamic limits of the Boltz-mann equation . Springer, Berlin, 2009.[8] S. Olla, S. R. S. Varadhan, and H-T. Yau. Hydrody-namical limit for a hamiltonian system with weak noise.
Commun. Math. Phys. , 155(3):523–560, 1993.[9] G.I. Taylor. The formation of a blast wave by a veryintense explosion I. Theoretical discussion.
Proc. R. Soc.A , 201(1065):159–174, 1950.[10] G.I. Taylor. The formation of a blast wave by a very in-tense explosion. - II. The atomic explosion of 1945.
Proc.R. Soc. A , 201(1065):175–186, 1950.[11] J. Von Neumann.
John Von Neumann: Collected Works.Theory of games, astrophysics, hydrodynamics and me-teorology . Number 6. Pergamon Press, 1963.[12] L. I. Sedov. Propagation of strong shock waves.
J. Appl.Math. Mech. , 10:241–250, 1946.[13] L. I. Sedov.
Similarity and Dimensional Methods in Me-chanics . Academic Press, New York, 1959.[14] L. D. Landau and E. M. Lifshitz.
Fluid Mechanics . Perg-amon Press, New York, 1987.[15] Y. B. Zel’dovich and Y. P. Raizer.
Physics of ShockWaves and High-Temperature Hydrodynamic Phenom-ena . Academic Press, New York, 1967.[16] T. Antal, P. L. Krapivsky, and S. Redner. Exciting hardspheres.
Phys. Rev. E , 78:030301, Sep 2008.[17] Z Jabeen, R. Rajesh, and P. Ray. Universal scaling dy-namics in a perturbed granular gas.
Europhys. Lett. ,89(3):34001, Feb 2010.[18] M. Barbier, D. Villamaina, and E. Trizac. Blast dynamicsin a dissipative gas.
Phys. Rev. Lett. , 115:214301, Nov2015.[19] M. Barbier. Kinetics of blast waves in one-dimensionalconservative and dissipative gases.
J. Stat. Mech. ,2015(11):P11019, Nov 2015.[20] M. Barbier, D. Villamaina, and E. Trizac. Microscopicorigin of self-similarity in granular blast waves.
Phys. Fluids , 28(8):083302, 2016.[21] J. P. Joy, S. N. Pathak, D. Das, and R. Rajesh. Shockpropagation in locally driven granular systems.
Phys.Rev. E , 96:032908, Sep 2017.[22] J. P. Joy, S. N. Pathak, and R. Rajesh. Shockpropagation following an intense explosion: compari-son between hydrodynamics and simulations. preprintarXiv:1812.03638 , 2018.[23] J. P. Joy and R. Rajesh. Shock propagation in the hardsphere gas in two dimensions: comparison between sim-ulations and hydrodynamics. preprint arXiv:1907.03416 ,2019.[24] P. L. Garrido, P. I. Hurtado, and B. Nadrowski. Simpleone-dimensional model of heat conduction which obeysfourier’s law.
Phys. Rev. Lett. , 86(24):5486, 2001.[25] A Dhar. Heat conduction in a one-dimensional gas ofelastically colliding particles of unequal masses.
Phys.Rev. Lett. , 86:3554–3557, Apr 2001.[26] P Grassberger, W Nadler, and L Yang. Heat conduc-tion and entropy production in a one-dimensional hard-particle gas.
Phys. Rev. Lett. , 89:180601, Oct 2002.[27] G. Casati and T. Prosen. Anomalous heat conduction ina one-dimensional ideal gas.
Phys. Rev. E , 67(1):015203,2003.[28] P Cipriani, S Denisov, and A Politi. From anomalous en-ergy diffusion to levy walks and heat conductivity in one-dimensional systems.
Phys. Rev. Lett. , 94(24):244301,2005.[29] S. Chen, J. Wang, G. Casati, and G. Benenti. Nonin-tegrability and the Fourier heat conduction law.
Phys.Rev. E , 90(3):032134, 2014.[30] P. I. Hurtado and P. L. Garrido. A violation of univer-sality in anomalous Fourier’s law.
Sci. Rep. , 6(1):38823,2016.[31] H. Zhao and W. Wang. Fourier heat conduction as astrong kinetic effect in one-dimensional hard-core gases.
Phys. Rev. E , 97:010103, Jan 2018.[32] S. Lepri, R. Livi, and A. Politi. Too close to inte-grable: Crossover from normal to anomalous heat dif-fusion.
Phys. Rev. Lett. , 125:040604, Jul 2020.[33] Y. Du, H. Li, and L. P. Kadanoff. Breakdown of hydrody-namics in a one-dimensional system of inelastic particles.
Phys. Rev. Lett. , 74:1268–1271, Feb 1995.[34] P. I. Hurtado. Breakdown of hydrodynamics in a simpleone-dimensional fluid.
Phys. Rev. Lett. , 96(1):010601,2006.[35] C. B. Mendl and H. Spohn. Shocks, rarefaction waves,and current fluctuations for anharmonic chains.
J. Stat.Phys. , 166(3-4):841–875, Oct 2016.[36] X. Cao, V. B. Bulchandani, and J. E. Moore. In-complete thermalization from trap-induced integrabilitybreaking: Lessons from classical hard rods.
Phys. Rev.Lett. , 120:164101, Apr 2018.[37] G Casati and T Prosen. Mixing property of triangularbilliards.
Phys. Rev. Lett. , 83(23):4729, 1999.[38] See supplemental material.[39] G. I. Barenblatt.
Scaling, Self-similarity, and Intermedi-ate Asymptotics: Dimensional Analysis and Intermediate
Asymptotics . Cambridge Texts in Applied Mathematics.Cambridge University Press, 1996.[40] S. Ganapa, S. Chakraborti, P. L. Krapivsky, and A. Dhar.The Taylor - von Neumann - Sedov blast-wave solution:comparisons with microscopic simulations of a one di-mensional gas. preprint arXiv:2010.15868 , 2020.[41] J. A. Joseph, J. E. Thomas, M. Kulkarni, and A. G.Abanov. Observation of shock waves in a strongly inter- acting fermi gas.
Phys. Rev. Lett. , 106(15):150401, 2011.[42] V. B. Bulchandani, C. Karrasch, and J. E. Moore. Su-perdiffusive transport of energy in one-dimensional met-als.
PNAS , 117(23):12713–12718, 2020.[43] R. W. MacCormack. A numerical method for solving theequations of compressible viscous flow.