initial energy $E$ and zero momentum, but we now take only _one realization_
of the microscopic process. The comparison of the empirical density profiles
with the ensemble-averaged profiles is shown in Fig. 5.5, and we see a close
agreement between the two. Here we used a bin size $\Delta=50$ and, in order
to get more data points, we used sliding bins. The fluctuations seen in the
empirical profile are expected to vanish in the thermodynamic limit.

####  5.5.2 Dependence on Initial Conditions

The long-time scaling solution should be independent of the details of the
initial conditions. It may depend on two parameters, the total blast energy
$E$ and the background density $\rho_{\infty}$. We verify this in simulations
by starting with the box initial condition Eq. (5.3), with the same values of
$E=32$ and $\rho_{\infty}=1.5$. From the results plotted in Fig. 5.6, it is
clear that the profiles of all the three fields from the two different initial
conditions quickly converge and become almost indistinguishable.

(a)(b)(c)

Figure 5.6: The comparison of the evolution of three field profiles starting
from two different initial conditions corresponding to either a Gaussian
temperature profile (solid lines) or a box-profile (dotted lines). The total
energy was fixed at $E=32$, the background density was taken as
$\rho_{\infty}=1.5$, and the average was taken over $10^{5}$ realizations.

####  5.5.3 Dependence on the Initial Ensemble

In Fig. 5.7 we repeat the simulations with initial conditions chosen from the
ensemble (A) discussed in Sec. 5.3 where the energy $E$ and total zero initial
momentum value are fixed only on average. This means that there are
realizations where the energy can be larger or smaller than $E$ while the
momentum can have a non-zero value. As a result, this leads to tails
developing in the shock front and we lose agreement with the hydrodynamic
theory. The correct physically relevant initial condition is the fixed energy,
fixed momentum ensemble and so we will not discuss further the case of
ensemble (A).

(a)(b)(c)

(d)(e)(f)

Figure 5.7: Initial ensemble with fluctuations in initial total energy and
total momentum, with mean energy fixed at $E$ and mean momentum at 0 (a,b,c)
Evolution of the density, velocity, and temperature starting from the initial
conditions corresponding to a Gaussian initial temperature profile and an
average over $10^{4}$ initial conditions. (d,e,f) The $x/t^{2/3}$ scaling
demonstrates data collapse in the long time limit everywhere except at the
center. Comparing with Fig. 5.4 we see that the agreement with the TvNS
scaling solution, shown by dashed lines, is now much poorer with significant
deviations both at the blast center and near the shock fronts.

####  5.5.4 Validity of Local Thermal Equilibrium

The local equilibrium (LE) assumption is an important requirement for the
validity of the hydrodynamic description. In case of equilibrium, the velocity
distribution is Gaussian. Earlier works in two and three dimensions [24, 25]
find strong deviations from LE, through observations of the skewness and
kurtosis of the local velocity distributions, which indicate non-Gaussianity.
These could be a cause for the observed disagreement with the hydrodynamic
predictions. Here we show results for these quantities obtained from MD
simulations of our 1D system. We define the local skewness
$S(x,t)=\frac{1}{n(x,t)}\left\langle\sum[v_{i}-v(x,t)]^{3}\delta(q_{i}-x)\right\rangle$
and the local kurtosis
$\chi(x,t)=\frac{1}{n(x,t)}\left\langle\sum
m_{i}^{2}[v_{i}-v(x,t)]^{4}\delta(q_{i}-x)\right\rangle$, where the number
density field $n(x,t)=\left\langle\sum\delta(q_{i}-x)\right\rangle$. In 1D,
the Gaussian distribution of velocity fluctuation implies that the skewness
$S(x,t)=0$ and the kurtosis $\chi(x,t)=3T^{2}$.

Figure 5.8: The skewness $S$ and the kurtosis $\chi$, corresponding to the
local velocity distribution at different points inside the blast, plotted as a
function of the scaled variable $\xi=x/R(t)$. The skewness shown in (a) is
scaled by $[{\rm thermal~{}velocity}]^{3}$ while the kurtosis shown in (b) is
scaled by the expected thermal value of $3T^{2}(x,t)$. The results at two
different times are shown.

In Fig. 5.8, we plot the dimensionless scaled skewness and kurtosis with the
scaled variable $\xi=x/R(t)$ at two different times and find that the LE
assumption is satisfied quite accurately, and significantly more than what is
observed in higher dimensions [24, 25]. We note that there are small
deviations from LE which interestingly, peak in the core and shock regions.
While this could provide a possible explanation for the disagreement between
TvNS and simulations in the core region, we note however that we continue to
have good agreement between these theories near the shock front. We now
explore the structure of the shock in some more detail.

####  5.5.5 Structure of the Shock

At any moment, we have a finite number of moving particles surrounded by an
infinite sea of stationary particles. One can then define the shock positions
as that of the extreme moving particles. Denote by $X_{L}(t)$ and $X_{R}(t)$
the positions of the left-most and right-most moving particles. These
quantities fluctuate between different microscopic realizations. Fluctuations
are outside the scope of deterministic hydrodynamics. We probe the basic
features of fluctuations by measuring the mean front positions,
$\bar{X}_{L}(t)=\langle X_{L}(t)\rangle$ and $\bar{X}_{R}(t)=\langle
X_{R}(t)\rangle$, and the variances, $\sigma^{2}_{L}(t)=\langle
X^{2}_{L}(t)\rangle-\langle X_{L}(t)\rangle^{2}$ and
$\sigma^{2}_{R}(t)=\langle X^{2}_{R}(t)\rangle-\langle X_{R}(t)\rangle^{2}$.
Fig. 5.9a shows that the mean positions $\bar{X}_{L}(t)$ and $\bar{X}_{R}(t)$
closely follow the hydrodynamic predictions $\pm R(t)$. Fig. 5.9b shows that
the fluctuations remain bounded around the equilibrium value $\sigma_{\rm
eq}^{2}=L^{2}/(4N)$, which for our parameters gives $\sigma_{\rm
eq}^{2}\approx 45$. Thus we do not observe a broadening of the shock front
with increasing time. This behaviour is in sharp contrast to the moving piston
case studied in [93] where a diffusive broadening of the shock front was
observed.

Figure 5.9: Time evolution of (a) the mean shock positions
$\bar{X}_{L}(t),~{}\bar{X}_{R}(t)$, and (b) their fluctuations
$\sigma_{L}(t),~{}\sigma_{R}(t)$. The dashed line $R(t)$ is the TvNS
prediction. The MD simulations were done for $N=L=8000$ with an averaging over
$10^{4}$ realizations.

The hydrodynamic description assumes local equilibrium (LE). The LE fails
inside the shock where the hydrodynamic fields undergo large changes over a
small length scale (comparable to the mean free path, i.e., to the gap between
adjacent particles in one dimension). The Euler description, where shock
fronts are singular surfaces ignores the shock wave structure, while the NSF
description is just an uncontrolled approximation. The Boltzmann equation
approach provides an adequate framework for analyzing shock waves.
Unfortunately, it is notoriously difficult for theoretical analysis even for
infinitely strong shock waves [95, 96, 97].

###  5.6 Modelling by the Navier-Stokes-Fourier (NSF) Equations

We now compare microscopic simulations with results obtained by solving the
one-dimensional Navier-Stokes-Fourier (NSF) equations:

$\displaystyle\partial_{t}\rho+\partial_{x}(\rho v)=0,$ (5.9a)
$\displaystyle\partial_{t}(\rho v)+\partial_{x}(\rho
v^{2}+p)=\partial_{x}(\zeta\partial_{x}v),$ (5.9b)
$\displaystyle\partial_{t}(\rho e)+\partial_{x}(\rho ev+pv)=\partial_{x}(\zeta
v\partial_{x}v+\kappa\partial_{x}T),$ (5.9c)

supplemented with thermodynamic relations Eqs. (4.30). We aim to show that NSF
equations capture the “essence” of the AHP gas system better than the Euler
equations.

The dissipative effects are manifested by the bulk viscosity $\zeta$ and the
thermal conductivity $\kappa$ appearing in Eqs. (5.9). These transport
coefficients can depend on the fields. The Green-Kubo relations and kinetic
theory arguments [52, 98] lead to the temperature dependence $\zeta\sim
T^{1/2}$ and $\kappa\sim T^{1/2}$ characterizing hard-sphere gases in an
arbitrary spatial dimension. A recent numerical study [77] suggests the
density dependence $\kappa\sim\rho^{1/3}$ specific to one dimension. In our
numerical study we have thus used the forms $\zeta=D_{1}T^{1/2}$ and
$\kappa=D_{2}\rho^{1/3}T^{1/2}$, where $D_{1}$ and $D_{2}$ are constants. We
solve the NSF equations Eqs. (5.9) numerically subject to the same initial
conditions as in the microscopic simulations:
$\rho(x,0)=\rho_{\infty},v(x,0)=0$ and $E(x,0)$ given by Eq. (5.4). We use the
MacCormack method [32, 33] which is second order in both space and time, with
discretization $dx=0.1$ and $dt=0.001$. We set the amplitudes in the
dissipative terms to unity: $D_{1}=D_{2}=1$. In our numerical solution we took
system size L=4000 and evolved up to the time ($t\approx 6400$) at which the
shock front reached the boundaries.

The numerical solution at early times plotted in Fig. 5.10 shows that sharp
fronts quickly develop for all three fields. The late time numerical results
are plotted in Fig. 5.11 (a,b,c) where we show the profiles of the fields,
$\rho(x,t)$, $v(x,t)$ and $T(x,t)$ at different times. In Fig. 5.11 (d,e,f) we
verify the expected scaling forms
$\rho(x,t)=\rho_{\infty}+\widetilde{\rho}(x/t^{2/3})$,
$v(x,t)=t^{-1/3}\widetilde{v}(x/t^{2/3})$ and
$T(x,t)=t^{-2/3}\widetilde{T}(x/t^{2/3})$, but observe that there is again a
core region where we do not get a good data collapse for the density and
temperature fields. For comparison, the exact TvNS scaling solution is also
plotted (dark dashed lines) and we see now that the singularities at the
origin have now disappeared and the profiles are closer to those observed in
simulations confirming thus the important role of dissipation. In Fig. 5.12 we
show a comparison of the numerical predictions of the dissipative
hydrodynamics (at the longest available time) with MD simulations and the
exact solution of Euler equations. Simulation results and the numerical
solution of the NSF equations show that the TvNS scaling breaks down in the
core region of the blast and we should therefore seek a different solution in
this region. In the next section, we analyze the NSF equations to answer this
question.

(a)(b)(c)

Figure 5.10: Early time evolution of the hydrodynamic fields
$\rho(x,t),v(x,t),T(x,t)$ as obtained from a numerical solution of Eqs. (5.9),
starting from the initial conditions with Gaussian temperature profile as in
the simulations for Fig. 5.4. We see the formation of a sharp front.

(a)(b)(c)

(d)(e)(f)

Figure 5.11: (a,b,c) Evolution of the hydrodynamic fields
$\rho(x,t),v(x,t),T(x,t)$ as obtained from a numerical solution of Eqs. (5.9)
starting from the initial conditions with Gaussian temperature profile used in
the simulations for Fig. 5.4. (d,e,f) Plots of the scaling functions
$\widetilde{\rho},\widetilde{v},\widetilde{T}$ and a comparison with the TvNS
solution (black dashed lines).

(a)(b)(c)

Figure 5.12: Comparison of the scaled fields
$\widetilde{\rho},\widetilde{v},\widetilde{T}$ plotted as a function of the
scaled variable $x/t^{2/3}$, from microscopic simulation (at time $t=80000$)
given by red dash dotted lines, from dissipative hydrodynamics (at $t=6400$)
given by blue solid lines and from the exact TvNS solution of the Euler
equations given by black dashed lines.

###  5.7 Scaling near the Core of the Blast

We start with the exact TvNS solution in one dimension, which was derived in
the previous chapter. The scaling functions $G(\xi),V(\xi)$ and $Z(\xi)$ are
implicitly given by:

$\displaystyle
G=2^{\frac{16}{5}}\,(1-V)^{2}\,(3V-1)^{\frac{1}{5}}\,(3-4V)^{-\frac{11}{5}},$
(5.10a)
$\displaystyle\xi^{5}=2^{-\frac{4}{3}}\,(3V-1)^{2}\,V^{-\frac{10}{3}}\,(3-4V)^{-\frac{11}{3}},$
(5.10b) $\displaystyle Z=\frac{3(1-V)V^{2}}{3V-1}.$ (5.10c)

Near the center of the explosion, the scaled hydrodynamic quantities exhibit
the following behaviours:

$\displaystyle G$ $\displaystyle\simeq B\,|\xi|^{\frac{1}{2}}\,,~{}\quad
B=2^{\frac{16}{3}}\,3^{-\frac{1}{2}}\,5^{-\frac{11}{6}},$ (5.11a)
$\displaystyle 3V-1$
$\displaystyle\simeq\beta\,|\xi|^{\frac{5}{2}}\,,\quad~{}~{}\beta=2^{\frac{2}{3}}\,3^{-\frac{7}{2}}\,5^{\frac{11}{6}},$
(5.11b) $\displaystyle Z$ $\displaystyle\simeq C\,|\xi|^{-\frac{5}{2}}\,,\quad
C=2^{\frac{1}{3}}\,3^{\frac{3}{2}}\,5^{-\frac{11}{6}},$ (5.11c)

where $\xi=x/R(t)$. Equation (5.11b) follows from Eq. (5.10b). Substituting
Eq. (5.11b) into Eq. (5.10a) [resp. Eq. (5.10c)] one derives Eq. (5.11a)
[resp. Eq. (5.11c)].

In one dimension, the scaling form in Eq. (4.6) gives

$\rho=\rho_{\infty}\,G(\xi),~{}~{}v=\frac{2x}{3t}\,V(\xi),~{}~{}T=\frac{4\mu
x^{2}}{27t^{2}}\,Z(\xi).$ (5.12)

These equations are valid for all $x\in[-R(t),R(t)]$. The scaled variable is
still given by $\xi=|x|/R$ and varies in the $0\leq\xi\leq 1$ range. Equations
(5.11a)–(5.11c) and Eq. (5.12) determine the behavior of the original
hydrodynamic quantities near the center of explosion. Near the center of
explosion, the gas moves with velocity varying linearly with distance from the
center of explosion

$v\simeq\frac{2x}{9t}.$ (5.13)

The velocity vanishes at the origin, while the temperature diverges. Indeed,
using Eqs. (5.11c) and Eq. (5.12) together with $T=\mu c^{2}/\gamma=\mu
c^{2}/3$ one gets

$T\simeq\frac{4\mu C}{27}\,\frac{x^{2}}{t^{2}}\,\xi^{-\frac{5}{2}}.$ (5.14)

Since $\xi\sim|x|/t^{2/3}$, the temperature diverges as

$T\sim|x|^{-\frac{1}{2}}\,t^{-\frac{1}{3}}.$ (5.15)

The divergence of the temperature is unphysical. Simulation results (Sec. 5.4)
show that the temperature at the center of the explosion remains finite and
decreases with time. The divergence is the consequence of an idealized
modeling using the Euler equations that ignores dissipation. The divergence
indicates that the heat conduction term becomes non-negligible near the
origin.

Let us now analyze the source of this problem from a slightly different point
of view. In the previous chapter in Appendix 4.B, we derived the TvNS
exponents for the scaling solutions of the Euler equations as:

$a=-2/3,~{}b=0,c=-1/3,~{}d=-2/3.$ (5.16)

It is easy to verify that if we substitute this scaling in Eq. (5.9), the
disspation terms will have a factor of $t^{-1/3}$ and we therefore say that
they do not contribute in the long time limit. Now we say in the previous
paragraph that the heat conduction term becomes non-negligible near the core
of the blast. That can only happen if the TvNS scaling is not valid near the
core of the blast. Not just the scaling functions, but also the values of the
exponents $a,b,c$ and $d$. We should therefore expect different values of
$a,b,c$ and $d$ which lead to a finite contribution of the conduction term
near the core and also a different scaling solution. Therefore, we have to
rectify our analysis.

To rectify the analysis we now take into account the dissipative processes and
consider an analysis of the full hydrodynamic equations Eq. (5.9). As shown in
Appendix 5.A, the terms corresponding to viscous dissipation drop off in the
scaling limit and we effectively get the following equations for the fields
$\rho,v,T$:

$\displaystyle\partial_{t}\rho+\partial_{x}(\rho v)=0,$ (5.17a)
$\displaystyle\rho(\partial_{t}+v\partial_{x})v+\mu^{-1}\partial_{x}(\rho
T)=0,$ (5.17b)
$\displaystyle\rho^{3}(\partial_{t}+v\partial_{x})\left[\frac{T}{\rho^{2}}\right]=2\mu\partial_{x}\big{(}D_{2}T^{1/2}\,\rho^{1/3}\partial_{x}T\big{)}.$
(5.17c)

These are the basic governing equations for the core region. Using Eq. (5.17c)
we deduce an estimate

$\rho\,\frac{T}{t}\sim\frac{T^{3/2}}{x^{2}}\,\rho^{1/3},$ (5.18)

which is combined with $\rho\sim\sqrt{\xi}\sim|x|^{1/2}/t^{1/3}$ to give an
estimate of the temperature in the core:

$T\sim|x|^{\frac{14}{3}}\,t^{-\frac{22}{9}},$ (5.19)

where thermal conductivity plays a significant role. The estimates Eqs. (5.15)
and (5.19) are comparable at $|x|=X$ determined by

$X^{\frac{14}{3}}\,t^{-\frac{22}{9}}\sim
X^{-\frac{1}{2}}\,t^{-\frac{1}{3}}\,.$

Thus the size $X$ of the core grows as

$X\sim t^{\frac{38}{93}}\,.$ (5.20)

The relative size of the core compared to the size of the excited region
approaches to zero as $X\,t^{-\frac{2}{3}}\sim t^{-\frac{8}{31}}$.

The temperature at the center of the explosion is estimated from Eq. (5.15)
after setting $x=X$, with $X$ given by Eq. (5.20). Similarly, the density at
the center of the explosion is estimated from $\rho_{0}\sim X^{1/2}/t^{1/3}$
and Eq. (5.20). We thus arrive at the following asymptotic decay laws

$T_{0}\sim t^{-\frac{50}{93}}\,,\qquad\rho_{0}\sim t^{-\frac{4}{31}}.$ (5.21)

Using Eqs. (5.13), (5.20) we estimate the velocity near the core: $v_{0}\sim
X/t\sim t^{-\frac{55}{93}}$. The decay laws Eq. (5.21) agree better with
microscopic simulations than the TvNS solution as shown in Fig. 5.13.

The hydrodynamic variables in the hot core where dissipative effects play the
dominant role exhibit scaling behaviors, but instead of $\xi=x/R$ we should
use the scaled spatial coordinate $\eta=x/X$. Thus the hydrodynamic variables
have a self-similar form

$\rho=t^{-\frac{4}{31}}\widetilde{G}(\eta),\quad
v=t^{-\frac{55}{93}}\widetilde{V}(\eta),\quad
T=t^{-\frac{50}{93}}\widetilde{Z}(\eta)$ (5.22)

in the scaling region

$x\to\infty,\quad t\to\infty,\quad\eta\sim
x\,t^{-\frac{38}{93}}\to\text{finite}.$ (5.23)

The scaling functions $\widetilde{G},\widetilde{V}$ and $\widetilde{Z}$ are
now different from $G,V$ and $Z$. The data from the NSF numerics and the
microscopic simulations satisfy these scaling forms and give a better collapse
than the TvNS scaling, see Fig. 5.14. All that is remaining now is to find
these scaling functions $\widetilde{G},\widetilde{V}$ and $\widetilde{Z}$ near
the core.

We now combine the scaling ansatz Eqs. (5.22)–(5.23) with the governing
equations (5.17) and we set $\mu=1.5,\ D_{2}=1$. Eq. (5.17b) yields the most
simple outcome $-$ it becomes $(\widetilde{G}\widetilde{Z})^{\prime}=0$ in the
leading order. (Hereinafter prime denotes differentiating with respect to
$\eta$.) Thus $\widetilde{G}=1/\widetilde{Z}$, where for the moment we have
set the integration constant to one (we will discuss this point again later).
Eq. (5.17a) becomes

$\widetilde{V}^{\prime}-\tfrac{4}{31}=\big{(}\widetilde{V}-\tfrac{38}{93}\eta\big{)}(\ln\widetilde{Z})^{\prime},$
(5.24a) while Eq. (5.17c) becomes
$\big{(}\widetilde{Z}^{\frac{1}{6}}\,\widetilde{Z}^{\prime}\big{)}^{\prime}+\tfrac{13}{93}=\tfrac{3}{2}\big{(}\widetilde{V}-\tfrac{38}{93}\eta\big{)}(\ln\widetilde{Z})^{\prime}.$
(5.24b)

Comparing Eqs. (5.24a) and (5.24b) we obtain

$\big{(}\widetilde{Z}^{\frac{1}{6}}\,\widetilde{Z}^{\prime}\big{)}^{\prime}+\tfrac{13}{93}=\tfrac{3}{2}\big{(}\widetilde{V}^{\prime}-\tfrac{4}{31}\big{)},$

which is integrated to give

$\widetilde{V}=\tfrac{2}{9}\eta+\tfrac{2}{3}\widetilde{Z}^{\frac{1}{6}}\,\widetilde{Z}^{\prime}.$
(5.25)

The integration constant vanishes since $\widetilde{V}(0)=0$ and
$\widetilde{Z}^{\prime}(0)=0$ due to symmetry. Inserting Eq. (5.25) into Eq.
(5.24b) we obtain

$\big{(}\widetilde{Z}^{\frac{1}{6}}\,\widetilde{Z}^{\prime}\big{)}^{\prime}+\tfrac{13}{93}=\big{(}-\tfrac{26}{93}\eta+\widetilde{Z}^{\frac{1}{6}}\,\widetilde{Z}^{\prime}\big{)}(\ln\widetilde{Z})^{\prime}.$
(5.26)

Solving Eq. (5.26) numerically [equivalently Eqs. (5.24b, 5.25)], with
$\widetilde{Z}(0)$ fixed from the NSF numerical data at $t=6400$ and
$\widetilde{Z}^{\prime}(0)=0$, we find the scaling functions
$\widetilde{G},\widetilde{V},\widetilde{Z}$. In Fig. 5.14, we plot these
functions and compare with the data from the NSF. We find qualitative
agreement in the form of the scaling functions but not a complete agreement.
See Appendix 5.A for a more detailed derivation of the equations for
$\widetilde{G},\widetilde{V},\widetilde{Z}$, where we also introduce an
additional scale factor $\lambda$ to make $\eta$ dimensionless. This gives
another fitting parameter but does not lead to a significant improvement in
comparison with data.

(a)(b)

Figure 5.13: Time evolution of (a) density and (b) temperature at the center
of the blast. The data is taken from microscopic simulations and the NSF
numerics. The scaling pre-factors correspond to the predictions from TvNS and
the core-scaling analysis. For the density, it is difficult to differentiate
between the two scaling forms, while for temperature, the scaling prediction
$T_{0}\sim t^{-50/93}$ in the core describes the data better than the TvNS
prediction $T_{0}\sim t^{-2/3}$.

(a)(b)(c)

Figure 5.14: A numerical verification of the core scaling forms Eqs. (5.22)
for data for the fields $\rho,~{}v$ and $T$ from both the microscopic
simulations and the NSF hydrodynamic numerics. The black dashed lines show the
analytic scaling functions obtained from a numerical solution of Eqs. (5.24b,
5.25).

The asymptotic behaviors of the velocity field are

$\displaystyle\widetilde{V}\simeq\tfrac{4}{31}\eta\quad\qquad\qquad\text{when}\quad|\eta|\ll
1,$ (5.27a)
$\displaystyle\tfrac{2}{9}\eta-\widetilde{V}\sim\eta^{-\frac{19}{12}}\qquad\text{when}\quad\eta\gg
1.$ (5.27b)

The asymptotic of Eq. (5.27a) follows from Eq. (5.24a). The slope predicted by
Eq. (5.27a) is a bit smaller than the slope $\frac{2}{9}$ in the blast region.
The asymptotic Eq. (5.27b) follows from Eq. (5.25). Further,

$\widetilde{Z}\sim\eta^{-1/2}\qquad\text{when}\quad\eta\gg 1$ (5.28)

follows from Eq. (5.26). Indeed,
$\widetilde{Z}^{\frac{1}{6}}\,\widetilde{Z}^{\prime}$ is asymptotically
negligible, so Eq. (5.26) simplifies to
$(\ln\widetilde{Z})^{\prime}=(2\eta)^{-1}$ leading to Eq. (5.28).

Matching of the inner and outer solutions: The asymptotic behavior Eq. (5.28)
yields $\widetilde{G}=\widetilde{Z}^{-1}\sim\eta^{1/2}$ matching the small
$\xi$ behavior Eq. (5.11a). For the temperature in the outer region we have

$\displaystyle
T\sim\left(\frac{x}{t}\right)^{2}Z(\xi)\sim\left(\frac{x}{t}\right)^{2}\left(\frac{x}{t^{2/3}}\right)^{-\frac{5}{2}}\sim
x^{-1/2}t^{-1/3},$ (5.29)

while in the inner region

$\displaystyle T=\frac{\widetilde{Z}(\eta)}{t^{50/93}}\sim
t^{-50/93}\left(\frac{x}{t^{38/93}}\right)^{-1/2}\sim x^{-1/2}t^{-1/3}.$
(5.30)

Thus we have perfect matching in the overlap of the inner and outer regions.
Similarly for the density

$\displaystyle\rho={G}(\xi)\sim(x/t^{2/3})^{1/2}\sim x^{1/2}t^{-1/3}$ (5.31)

in the outer region, while in the inner region

$\displaystyle\rho=t^{-4/31}\widetilde{G}(\eta)\sim
t^{-4/31}(x/t^{38/93})^{1/2}\sim x^{1/2}t^{-1/3},$ (5.32)

where we used $\widetilde{G}=1/\widetilde{Z}$. Again we get perfect matching.

###  5.8 Summary

We analyzed the connection between the atomistic description of matter with
particles obeying Newton’s equations and the continuum hydrodynamic
descriptions relying on the Euler equations (i.e., neglecting dissipative
effects) and Navier-Stokes-Fourier equations. Our system is a blast wave in a
one-dimensional gas initially at rest (zero temperature). The hydrodynamic
description is provided by $\rho(x,t),v(x,t),T(x,t)$ describing the mass
density, flow velocity, and temperature. We compared the exact TvNS solution
of the hydrodynamic Euler equations with results of direct molecular dynamics
simulations of a gas composed of hard point particles with alternating masses.
Our main results are:

  1. –

We obtained the exact TvNS scaling solution of the 1D Euler equation for the
ideal gas. The position of the shock is
$R(t)=[1071Et^{2}/(152\rho_{\infty})]^{1/3}$. The hydrodynamic fields are
expressed through the scaling variable $\xi=x/R(t)$.

  2. –

For the AHP gas, we find a remarkably close agreement between the simulation
results and the TvNS predictions. This includes a matching of the location of
the front $R(t)$ and of the scaling functions over most of the range of the
scaling variable $\xi$, except for a core region near the origin whose size
shrinks (when measured in $\xi$) at larger times.

  3. –

We obtained results both for the ensemble averaged fields and the empirical
fields and find agreement between them. We also showed that it is important to
use a microcanonical ensemble — otherwise the energy fluctuations in a
canonical ensemble lead to a broadening of the shock front.

  4. –

The deviations from the TvNS scaling at the core occur due to the contribution
of heat conduction that becomes non-negligible in this region. Euler’s
equations do not take these terms into account.

  5. –

We thus model the hydrodynamics of our system using the NSF equations, which
provided us a detailed understanding of the distinct scaling form in the core
region, whose size scales as $X\sim t^{38/93}$. The divergence of the
temperature field at the origin and the vanishing of the density field, which
are unphysical predictions of the TvNS solution, are avoided in the solutions
of the dissipative equations.

##  Appendix 5.A Derivation of Scaling Functions near the Core

Consider the heat equation,

$\rho\partial_{t}(T/2\mu)=\partial_{x}(\kappa\partial_{x}T),$ (5.33)

with $\kappa=D_{2}\rho^{1/3}T^{1/2}$. Near the origin, the ratio of the right-
hand side to the left-hand side of Eq. (5.33) is estimated from the TvNS
solution to yield

$\displaystyle\frac{T^{1/2}t}{\rho^{2/3}x^{2}}\sim\frac{1}{\xi^{5/4}\xi^{1/3}x}\sim\frac{t^{19/18}}{x^{31/12}}.$
(5.34)

We see that this becomes large in the region $x\lesssim t^{38/93}$. Hence we
expect a different scaling solution in the region $x<t^{-a}$ where $a=-38/93$.
As argued in Sec. 5.7, the expected scaling form in the core reads

$\displaystyle\rho=t^{b}\widetilde{G}(xt^{a}),~{}~{}v=t^{c}\widetilde{V}(xt^{a}),~{}~{}T=t^{d}\widetilde{Z}(xt^{a}),$
(5.35)

with $a=-\frac{38}{93},~{}b=-\frac{4}{31},~{}c=-\frac{55}{93}$ and
$d=-\frac{50}{93}$. We plug this scaling ansatz into the governing NSF
equations:

$\displaystyle\partial_{t}\rho+\partial_{x}(\rho v)=0,$ (5.36)
$\displaystyle\rho(\partial_{t}+v\partial_{x})v+\partial_{x}p=\partial_{x}(\eta\partial_{x}v),$
(5.37)
$\displaystyle\rho(\partial_{t}+v\partial_{x})\epsilon+p\partial_{x}v=\partial_{x}(\eta
v\partial_{x}v)+\partial_{x}(\kappa\partial_{x}T),$ (5.38)

with $\epsilon=T/(2\mu)$ and $p=\rho T/\mu$, and arrive at

$\displaystyle
b\widetilde{G}+a\eta\widetilde{G}^{\prime}+t^{a+c+1}\lambda(\widetilde{G}\widetilde{V})^{\prime}=0,$
(5.39) $\displaystyle
t^{c-1}(c\widetilde{G}\widetilde{V}+a\eta\widetilde{G}\widetilde{V}^{\prime})+t^{2c+a}\lambda\widetilde{G}\widetilde{V}\widetilde{V}^{\prime}+t^{a+d}\lambda\frac{(\widetilde{G}\widetilde{Z})^{\prime}}{\mu}$
$\displaystyle=t^{d/2+c+2a-b}\lambda^{2}(\widetilde{Z}^{1/2}\widetilde{V}^{\prime})^{\prime},$
(5.40) $\displaystyle
t^{b+d-2a-1}\,\frac{d\widetilde{G}\widetilde{Z}+a\eta\widetilde{G}\widetilde{Z}^{\prime}}{2\mu\lambda^{2}}+t^{b+c+d-a}\,\frac{\lambda\widetilde{G}\big{[}\widetilde{V}\widetilde{Z}^{\prime}+2\widetilde{Z}\widetilde{V}^{\prime}\big{]}}{2\mu\lambda^{2}}$
$\displaystyle=t^{d/2+2c}(\widetilde{Z}^{1/2}\widetilde{V}\widetilde{V}^{\prime})^{\prime}+t^{b/3+3d/2}(\widetilde{G}^{1/3}\widetilde{Z}^{1/2}\widetilde{Z}^{\prime})^{\prime}.$
(5.41)

Since $c=-1-a=-\frac{55}{93}$ and $b=3d/4+3a+3/2=-\frac{4}{31}$, the leading
term in Eq. (5.40) gives $(\widetilde{G}\widetilde{Z})^{\prime}=0$ implying
$\widetilde{G}=k/\widetilde{Z}$, where $k=\widetilde{G}(0)\widetilde{Z}(0)$ is
a constant. Using this then leads to the following equations for
$\widetilde{V}$ and $\widetilde{Z}$:

$\displaystyle\widetilde{V}^{\prime}-\frac{4}{31}+\left(\frac{38\eta}{93}-\widetilde{V}\right)\frac{\widetilde{Z}^{\prime}}{\widetilde{Z}}=0,$
(5.42)
$\displaystyle\widetilde{V}^{\prime}-\frac{25}{93}-\frac{1}{2}\left(\frac{38\eta}{93}-\widetilde{V}\right)\frac{\widetilde{Z}^{\prime}}{\widetilde{Z}}=\nu\lambda^{2}(\widetilde{Z}^{1/6}\widetilde{Z}^{\prime})^{\prime},$
(5.43)

where $\nu=\mu/k^{2/3}$. From these two equations we get

$\displaystyle
9\widetilde{V}^{\prime}-2=6\nu\lambda^{2}(\widetilde{Z}^{1/6}\widetilde{Z}^{\prime})^{\prime},$
(5.44)

which is integrated to give

$\displaystyle\widetilde{V}=2\eta/9+(2\nu\lambda^{2}/3)\widetilde{Z}^{1/6}\widetilde{Z}^{\prime}.$
(5.45)

Substituting this into Eq. (5.42) yields a non-linear ordinary differential
equation for $\widetilde{Z}$:

$\displaystyle\lambda^{2}\nu(\widetilde{Z}^{1/6}\widetilde{Z}^{\prime})^{\prime}+\left(\frac{26\eta}{93}-\lambda^{2}\nu\widetilde{Z}^{1/6}\widetilde{Z}^{\prime}\right)\frac{\widetilde{Z}^{\prime}}{\widetilde{Z}}+\frac{13}{93}=0.$

The constants $\widetilde{G}(0),\widetilde{V}(0),\widetilde{Z}(0)$ and
$\lambda$ should be fixed by matching the $\eta\gg 1$ behavior with the TvNS
solution at $\xi\ll 1$, which however did not lead to a perfect agreement as
shown in Fig. 5.14.

##  Chapter 6 Conclusion

###  6.1 Outcomes of the Study

In this work we studied various aspects of thermalization, chaos and
hydrodynamics in one dimensional Hamiltonian systems. Our motivation was two
folds. First, to understand what leads a system to thermalization and we
studied the role of initial conditions, choice of observables, choice of
averaging and the role of chaos. For this, we studied the $\alpha-$FPUT system
through the equipartition of energy among the local observables. And second,
we studied the evolution of the alternate mass hard particle (AHP) gas from a
nonequilibrium initial state and tried to understand the evolution of
conserved fields via hydrodynamics.

For the $\alpha-$FPUT system, we found that the equilibration time $\tau_{eq}$
of the local observables scaled with the dimensionless nonlinearity parameter
$\epsilon$ as $\epsilon^{-a}$, with the value of $a$ found to be between $4$
and $6$, dependent on the width of the initial distribution $\gamma$ used to
compute the ensemble averages. This was in contrast to the equilibration of
normal modes, where $a\approx 8$. We found that local observables equilibrate
on much shorter time scales when compared to the normal modes. For a given
$\epsilon$, we also found that $\tau_{eq}$ depended on the initial ensemble
and infact we observe $\tau_{eq}\propto\log\gamma$ with the slope being close
to the inverse of the maximal Lyapunov exponent of the system, thus
quantifying the relation between thermalization and chaos in the $\alpha-$FPUT
system. We also found that the Toda chain thermalizes for a broad enough
initial distribution. We tried to explain these differences on the basis of
their phase space dynamics.

We then studied the evolution of the AHP gas from “blast-wave” initial
condition until the time the energy reaches its boundary. This leads to the
formation of shock waves in the system and the evolution of the conserved
fields become self-similar in time. The scaling solution of the conserved
fields is given by the Taylor, von Neumann and Sedov (TvNS) solution. We find
that in our system, the evolution of conserved quantities matches perfectly
with the TvNS predictions everywhere except in a small region near the core of
the blast whose size shrinks at large times. We explain this deviation as
arising due to the non-negligible contribution of the conductivity term near
the core of the blast, which is not accounted for by the Euler equations. We
thus, model our system by using NSF equations - Euler equations with
additional dissipation terms due to viscosity and conductivity.

###  6.2 Limitations of the Study

For the $\alpha-$FPUT system we found that $\tau_{\rm eq}\propto\log(\gamma)$.
The inverse of the proportionality constant $3.3\times 10^{-4}$, is close to
the Lyapunov exponent $\Lambda\approx 4.1\times 10^{-4}$ of the system. We
believe that the properties of the initial conditions, such as its symmetries
and its vicinity from breather solutions can affect the exponential growth of
perturbations of the initial conditions, which would lead to a retardation of
thermalization and higher values of $\tau_{\rm eq}$. This could explain why we
did not get a closer agreement. Far away from breathers we can expect there is
no such effect. This needs to be investigated further and is beyond the scope
of this work.

For the alternate mass hard particle gas system, our understanding of the core
region remains incomplete. We don’t get a perfect agreement in the form of
scaling solutions near the core of the blast. A comprehensive agreement
between microscopic simulations and hydrodynamics is difficult to claim due to
the lack of knowledge of the dissipative terms in one dimension. Anomalous
transport in one dimensional systems may also have a say in this. This has to
be studied further.

###  6.3 Future Directions

####  6.3.1 Thermalization and Chaos

For the $\alpha-$FPUT system, our results indicate that the equilibration time
for local observables $\tau_{eq}$ scales with $\epsilon$ as $\epsilon^{-a}$,
with the value of $a$ found to be between $4$ and $6$, in contrast to the
equipartition of normal modes, where $a=8$. One problem is to understand the
source of this deviation. Why is this different? One way to study this is to
extend the formalism of wave turbulence to the case of space localized initial
conditions and observables. Another problem would be to study thermalization
in the AHP gas (a non-chaotic system), by studying the evolution of conserved
quantities for much longer times and compare results of the dependence of
$\tau_{eq}$ on $\gamma$ with the FPUT system, which can give us more clues
regarding the role of chaos.

Thermalization in finite Hamiltonian systems has usually been studied by
either considering a time averaging protocol or an ensemble averaging
protocol. As we illustrate, they can lead to very different estimates for the
time scale of equilibration. In our example, the time averaging protocol gives
a thermalization time that is several orders of magnitude larger than that
obtained from ensemble averaging. We believe that the ensemble averaging
protocol is relevant for understanding aspects of the classical-quantum
correspondence in the context of thermalization in finite systems. For quantum
systems several studies show, e.g [99], that a finite quantum system prepared
in a pure initial state and evolving under unitary dynamics can exhibit
thermalization (without requiring any time averaging). A corresponding
statement for the classical system is difficult if (corresponding to the
quantum pure state) one considers a single point in phase space. However if we
smear out the point into an initial blob in phase space, as is done in our
study, then we can get an equivalent classical statement. The smearing out
process can be thought of as arising from the uncertainty principle. As an
example we point out recent work [100] on the quantum-classical correspondence
in Floquet systems where such an averaging protocol is essential.

We note that another important question is that of thermalization of
macroscopic systems. There we expect that thermalization requires neither a
temporal nor an ensemble averaging protocol, but arises from the fact that
physical observables are macro-variables and their measurement typically
involves an averaging over many degrees of freedom [30, 31]. An interesting
problem would be to compare the results of the spatial averaging protocol with
that of time and ensemble averaging protocols in the thermodynamic limit. This
would give us a better understanding thermalization of macroscopic systems.
Studying the evolution of conserved quantities in the FPUT system through
hydrodynamics is challenging because the equation of state is not known
analytically and can only be computed numerically, nevertheless this would be
an interesting study.

####  6.3.2 Hydrodynamics

Earlier studies [23, 24, 25, 26, 70, 71] of hard-sphere gases in two and three
dimensions were unable to obtain a clear and convincing agreement between
simulations and the TvNS solution even in the non-core region. Some possible
reasons for this are: (i) for our system, an ideal gas, the equilibrium
properties including the equation of state are known exactly and hence one can
obtain an exact solution; (ii) simulations in the 1D AHP gas can be done very
accurately, the system has good ergodicity properties, and the asymptotic
behavior is easy to study. Following our work, an agreement between molecular
dynamics simulations and hydrodynamics equations was reported in higher
dimensions [101], once dissipative terms are included.

One of the most interesting extensions of the present work would be to
consider the effect of potential energy terms in the Hamiltonian, e.g to study
systems such as Lennard-Jones gases or lattice models like the Fermi-Pasta-
Ulam-Tsingou (FPUT) chain. For the FPUT model, the Euler equations are known
[90], however the equation of state is not known in closed form and this makes
a detailed comparison between simulations and hydrodynamic predictions quite
nontrivial. Recently, the Euler equations along with Navier-Stokes
corrections, for integrable Hamiltonian systems such as the Toda lattice and
the hard rod gas have been written and the evolution of domain wall initial
conditions investigated [102, 103, 104]. The effect of integrability would be
another interesting direction to explore. For integrable systems, we use
generalized hydrodynamics and for non-integrable systems we use the usual
hydrodynamics. It will be very interesting to consider what happens with near-
integrable systems, and if generalized hydrodynamics is valid in these
systems. The evolution of blast-wave initial conditions in these systems, and
in particular the nature of self-similar solutions of the conserved fields are
fascinating open questions.

The emergence of the mean-field hydrodynamic behavior is questionable in low
dimensional systems [98] at low temperatures. Contrary to these generic
expectations, our extremal one-dimensional zero-temperature system disturbed
by sudden local release of energy demonstrates the hydrodynamic behaviour
surprisingly well. Our understanding of the core region remains incomplete. A
comprehensive agreement between microscopic simulations and hydrodynamics is
difficult to claim due to the lack of knowledge of the dissipative terms in
one dimension. The well-known phenomenon of anomalous heat transport in one-
dimensional systems (with mass, momentum, and energy conserved) means that the
thermal conductivity $\kappa$ is not given by the standard Green-Kubo formula,
see [105] for a recent resolution.

####  6.3.3 Concluding Remarks

I would like to conclude this thesis by pointing out that thermalization and
hydrodynamics in Hamiltonian systems are far from being completely understood.
On the one hand, we use statistical physics to understand the emergent
behaviour in macroscopic systems. On the other, there is a lot to understand
about finite sized systems from the point of view of phase space dynamics. It
is possible that Hamiltonian chaos [106, 107, 108] could give us a deeper
understanding regarding the transition from finite sized systems to
macroscopic systems. However, there are some studies [109, 110, 111, 112, 113]
that provide evidence of chaos being misleading, misguided or unnecessary in
the context of discussing thermalization of macroscopic systems. It would
nevertheless be very interesting to probe this question further and actually
evaluate the exact role of chaos. And one wonders if Bob Dylan, an American
singer is correct when he said “I accept chaos, I’m not sure whether it
accepts me.”

## Bibliography

  * [1] E. Fermi, J. Pasta, and S. Ulam.  Studies of the nonlinear problems, (1955).  Los Alamos National Laboratory Report (LA1940).
  * [2] T. Dauxois.  Fermi, Pasta, Ulam, and a mysterious lady.  Physics today, 61(1), (2008).
  * [3] G. I. Taylor.  The formation of a blast wave by a very intense explosion - I. Theoretical discussion.  Proc. Roy. Soc. London. Series A, 201(1065), (1950).
  * [4] G. I. Taylor.  The formation of a blast wave by a very intense explosion - II. The atomic explosion of 1945.  Proc. Roy. Soc. London. Series A, 201(1065), (1950).
  * [5] J. Von Neumann.  John Von Neumann: Collected Works. Theory of games, astrophysics, hydrodynamics and meteorology.  Number 6. Pergamon Press, (1963).
  * [6] L. I. Sedov.  Propagation of strong shock waves.  J. Appl. Math. Mech., 10, (1946).
  * [7] L. I. Sedov.  Similarity and Dimensional Methods in Mechanics.  Elsevier Science, (2014).
  * [8] T. P. Weissert.  The genesis of simulation in dynamics: pursuing the Fermi-Pasta-Ulam problem.  Springer, (2012).
  * [9] G. P. Berman and F. M. Izrailev.  The Fermi-Pasta-Ulam problem: Fifty years of progress.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 15(1), (2005).
  * [10] G. Gallavotti.  Introduction to FPU.  In Giovanni Gallavotti, editor, The Fermi-Pasta-Ulam Problem: A Status Report. Springer Berlin Heidelberg, (2008).
  * [11] M. Onorato, L. Vozella, D. Proment, and Y.V. Lvov.  Route to thermalization in the $\alpha$-Fermi–Pasta–Ulam system.  PNAS, 112(14):4208–4213, (2015).
  * [12] S. Ganapa, A. Apte, and A. Dhar.  Thermalization of local observables in the $\alpha$-FPUT chain.  Journal of Statistical Physics, 180(1):1010–1030, (2020).
  * [13] S. Ganapa, S. Chakraborti, P.L. Krapivsky, and A. Dhar.  Blast in the one-dimensional cold gas: Comparison of microscopic simulations with hydrodynamic predictions.  Physics of Fluids, 33(8):087113, (2021).
  * [14] S. Chakraborti, S. Ganapa, P.L. Krapivsky, and A. Dhar.  Blast in a one-dimensional cold gas: From Newtonian dynamics to hydrodynamics.  Physical Review Letters, 126(24):244503, (2021).
  * [15] V. E. Zakharov V. S. L'vov and G. Falkovich.  Statistical description of weak wave turbulence.  In Kolmogorov Spectra of Turbulence I: Wave Turbulence, pages 63–82. Springer Berlin Heidelberg, (1992).
  * [16] S. Nazarenko.  Wave Turbulence.  Springer Berlin Heidelberg, (2011).
  * [17] P Mazur and E Montroll.  Poincaré cycles, ergodicity, and irreversibility in assemblies of coupled harmonic oscillators.  Journal of mathematical physics, 1(1):70–84, (1960).
  * [18] EM Titulaer.  Ergodic features of harmonic-oscillator systems. I.  Physica, 70(2):257–275, (1973).
  * [19] UM Titulaer.  Ergodic features of harmonic-oscillator systems. II. asymptotic ergodicity.  Physica, 70(2):276–296, (1973).
  * [20] C. Danieli, D. K. Campbell, and S. Flach.  Intermittent many-body dynamics at equilibrium.  Phys. Rev. E, 95(6):060202, (2017).
  * [21] L. D. Landau and E. M. Lifshitz.  Fluid Mechanics.  Pergamon Press, New York, (1987).
  * [22] Y. B. Zel’dovich and Y. P. Raizer.  Physics of Shock Waves and High-Temperature Hydrodynamic Phenomena.  Academic Press, New York, (1967).
  * [23] M. Barbier, D. Villamaina, and E. Trizac.  Microscopic origin of self-similarity in granular blast waves.  Phys. Fluids, 28:083302, (2016).
  * [24] J. P. Joy, S. N. Pathak, and R. Rajesh.  Shock propagation following an intense explosion: comparison between hydrodynamics and simulations.  J. Stat. Phys., 182:1–22, (2021).
  * [25] J. P. Joy and R. Rajesh.  Shock propagation in the hard sphere gas in two dimensions: Comparison between simulations and hydrodynamics.  Journal of Statistical Physics, 184(3):1–22, (2021).
  * [26] J. P. Joy, S. N. Pathak, D. Das, and R. Rajesh.  Shock propagation in locally driven granular systems.  Phys. Rev. E, 96:032908, (2017).
  * [27] H. Yoshida.  Construction of higher order symplectic integrators.  Physics Letters A, 150(5-7):262–268, (1990).
  * [28] G. Benettin, H. Christodoulidi, and A. Ponno.  The Fermi-Pasta-Ulam Problem and Its Underlying Integrable Dynamics.  Journal of Statistical Physics, 152:195–212, (2013).
  * [29] G. Benettin, S. Pasquali, and A. Ponno.  The Fermi–Pasta–Ulam problem and its underlying integrable dynamics: An approach through Lyapunov exponents.  Journal of Statistical Physics, 171(4):521–542, Mar (2018).
  * [30] J. L. Lebowitz.  Statistical mechanics: A selective review of two central issues.  Reviews of Modern Physics, 71(2):S346–S357, Mar (1999).
  * [31] S. Goldstein, D. A. Huse, J. L. Lebowitz, and R. Tumulka.  Macroscopic and microscopic thermal equilibrium.  Annalen der Physik, 529(7):1600301, Feb (2017).
  * [32] R. W. MacCormack.  The effect of viscosity in hypervelocity impact cratering.  Journal of spacecraft and rockets, 40(5):757–763, (2003).
  * [33] R. W. MacCormack.  A numerical method for solving the equations of compressible viscous flow.  AIAA journal, 20:1275–1281, (1982).
  * [34] Wikipedia.  MacCormack method - Wikipedia, the free encyclopedia.  http://en.wikipedia.org/w/index.php?title=MacCormack%20method&oldid=1049809004, (2022).  [Online; accessed 20-April-2022].
  * [35] J. Ford.  The Fermi-Pasta-Ulam problem: Paradox turns discovery.  Physics Reports, 213(5):271 – 310, (1992).
  * [36] J. Marin and S. Aubry.  Breathers in nonlinear lattices: numerical calculation from the anticontinuous limit.  Nonlinearity, 9:1501–1528, (1996).
  * [37] S. Flach, M. V. Ivanchenko, and O. I. Kanakov.  $q$-Breathers and the Fermi-Pasta-Ulam problem.  Phys. Rev. Lett., 95(6):064102, (2005).
  * [38] S. Flach, M. V. Ivanchenko, and O. I. Kanakov.  $q$-breathers in Fermi-Pasta-Ulam chains: Existence, localization, and stability.  Phys. Rev. E, 73(3):036618, (2006).
  * [39] H. Christodoulidi, C. Efthymiopoulos, and T. Bountis.  Energy localization on $q$-tori, long-term stability, and the interpretation of Fermi-Pasta-Ulam recurrences.  Phys. Rev. E, 81(1):016210, (2010).
  * [40] Y. V. Lvov and M. Onorato.  Double scaling in the relaxation time in the $\beta$-Fermi-Pasta-Ulam-Tsingou model.  Phys. Rev. Lett., 120:144301, (2018).
  * [41] L. Pistone, S. Chibbaro, M. Bustamante, Y. L’vov, and M. Onorato.  Universal route to thermalization in weakly-nonlinear one-dimensional chains.  Mathematics in Engineering, 1(4):672, (2019).
  * [42] F. M. Izrailev and B. V. Chirikov.  Statistical properties of a nonlinear string.  Dokl. Akad. Nauk SSSR, 166(1):57–59, (1966).
  * [43] L. Casetti, M. Cerruti-Sola, M. Pettini, and E. Cohen.  The Fermi-Pasta-Ulam problem revisited: Stochasticity thresholds in nonlinear Hamiltonian systems.  Phys. Rev. E, 55(6), (1996).
  * [44] J. Deluca, A. J. Lichtenberg, and S. Ruffo.  Energy transitions and time scales to equipartition in the Fermi-Pasta-Ulam oscillator chain.  Phys. Rev. E, 51(4):2877–2885, (1995).
  * [45] R. Livi, M. Pettini, S. Ruffo, M. Sparpaglione, and A. Vulpiani.  Equipartition threshold in nonlinear large Hamiltonian systems: The Fermi-Pasta-Ulam model.  Phys. Rev. A, 31(2):1039–1045, (1985).
  * [46] N. J. Zabusky and M. D. Kruskal.  Interaction of “Solitons” in a Collisionless Plasma and the Recurrence of Initial States.  Phys. Rev. Lett., 15:240–243, (1965).
  * [47] M. Toda.  Vibration of a chain with nonlinear interaction.  Journal of the Physical Society of Japan, 22(2):431–436, (1967).
  * [48] M. Toda.  Studies of a non-linear lattice.  Physics Reports, 18(1):1 – 123, (1975).
  * [49] T. Goldfriend and J. Kurchan.  Equilibration of quasi-integrable systems.  Phys. Rev. E, 99(2), Feb (2019).
  * [50] W. Fu, Y. Zhang, and H. Zhao.  Universal law of thermalization for one-dimensional perturbed Toda lattices.  New Journal of Physics, 21(4):043009, Apr (2019).
  * [51] G. I. Barenblatt.  Scaling, Self-similarity, and Intermediate Asymptotics: Dimensional Analysis and Intermediate Asymptotics.  Cambridge Texts in Applied Mathematics. Cambridge University Press, (1996).
  * [52] P. Resibois and M. De Leener.  Classical Kinetic Theory of Fluids.  Wiley, (1977).
  * [53] J. P. Ostriker and C. F. McKee.  Astrophysical blastwaves.  Rev. Mod. Phys., 60, (1988).
  * [54] C. F. McKee and J. K. Truelove.  Explosions in the interstellar medium.  Phys. Reports, 256:157–172, (1995).
  * [55] E. Waxman and R. Shvarts.  Asymptotic self-similar solutions with a characteristic timescale.  Astrophys. J., 721:1928, (2010).
  * [56] D. Kushnir and B. Katz.  Early hydrodynamic evolution of a stellar collision.  Astrophys. J., 785:124, (2014).
  * [57] D. Branch and J. C. Wheeler.  Supernova Explosions.  Springer, (2017).
  * [58] E. Modelevsky and R. Sari.  Revisiting the strong shock problem: Converging and diverging shocks in different geometries.  Phys. Fluids, 33:056105, (2021).
  * [59] M. J. Edwards, A. J. MacKinnon, J. Zweiback, K. Shigemori, D. Ryutov, A. M. Rubenchik, K. A. Keilty, E. Liang, B. A. Remington, and T. Ditmire.  Investigation of ultrafast laser-driven radiative blast waves.  Phys. Rev. Lett., 87:085004, (2001).
  * [60] A. S. Moore, D. R. Symes, and R. A. Smith.  Tailored blast wave formation: Developing experiments pertinent to laboratory astrophysics.  Phys. Plasmas, 12:052707, (2005).
  * [61] J. F. Boudet, J. Cassagne, and H. Kellay.  Blast shocks in quasi-two-dimensional supersonic granular flows.  Phys. Rev. Lett., 103:224501, (2009).
  * [62] J. F. Boudet and H. Kellay.  Unstable blast shocks in dilute granular flows.  Phys. Rev. E, 87:052202, (2013).
  * [63] R. D. Blandford and C. F. McKee.  Fluid dynamics of relativistic blast waves.  Phys. Fluids, 19(8):1130–1138, (1976).
  * [64] P. Best and R. Sari.  Second-type self-similar solutions to the ultrarelativistic strong explosion problem.  Phys. Fluids, 12(11):3029–3035, (2000).
  * [65] R. Sari.  First and second type self-similar solutions of implosions and explosions containing ultrarelativistic shocks.  Phys. Fluids, 18(2):027106, (2006).
  * [66] J. Tian.  Revisiting first type self-similar solutions of explosions containing ultrarelativistic shocks.  Phys. Fluids, 30(5):056102, (2018).
  * [67] E. Coughlin.  Energy-conserving, relativistic corrections to strong shock propagation.  Astrophys. J., 880:108, (2019).
  * [68] T. Antal, P. L. Krapivsky, and S. Redner.  Exciting hard spheres.  Phys. Rev. E, 78:030301, (2008).
  * [69] Z. Jabeen, R. Rajesh, and P. Ray.  Universal scaling dynamics in a perturbed granular gas.  EPL, 89:34001, (2010).
  * [70] M. Barbier, D. Villamaina, and E. Trizac.  Blast dynamics in a dissipative gas.  Phys. Rev. Lett., 115:214301, (2015).
  * [71] M. Barbier.  Kinetics of blast waves in one-dimensional conservative and dissipative gases.  JSTAT, (2015):P11019, (2015).
  * [72] A Dhar.  Heat conduction in a one-dimensional gas of elastically colliding particles of unequal masses.  Phys. Rev. Lett., 86:3554–3557, (2001).
  * [73] P Grassberger, W Nadler, and L Yang.  Heat conduction and entropy production in a one-dimensional hard-particle gas.  Phys. Rev. Lett., 89:180601, (2002).
  * [74] G. Casati and T. Prosen.  Anomalous heat conduction in a one-dimensional ideal gas.  Phys. Rev. E, 67:015203, (2003).
  * [75] P Cipriani, S Denisov, and A Politi.  From anomalous energy diffusion to levy walks and heat conductivity in one-dimensional systems.  Phys. Rev. Lett., 94:244301, (2005).
  * [76] S. Chen, J. Wang, G. Casati, and G. Benenti.  Nonintegrability and the Fourier heat conduction law.  Phys. Rev. E, 90:032134, (2014).
  * [77] P. I. Hurtado and P. L. Garrido.  A violation of universality in anomalous Fourier’s law.  Scientific reports, 6:1–10, (2016).
  * [78] S. Lepri, R. Livi, and A. Politi.  Too close to integrable: Crossover from normal to anomalous heat diffusion.  Phys. Rev. Lett., 125:040604, (2020).
  * [79] R. K. Zeytounian.  Navier-Stokes-Fourier equations: a rational asymptotic modelling point of view.  Springer Science & Business Media, (2012).
  * [80] D. Ryu and E. T. Vishniac.  The growth of linear perturbations of adiabatic shock waves.  Ap. J., 313:820–841, (1987).
  * [81] J. Sanz, S. E. Bouquet, C. Michaut, and J. Miniere.  The spectrum of the Sedov–Taylor point explosion linear stability.  Phys. Plasmas, 23:062114, (2016).
  * [82] P. L. Garrido, P. I. Hurtado, and B. Nadrowski.  Simple one-dimensional model of heat conduction which obeys Fourier law.  Phys. Rev. Lett., 86:5486, (2001).
  * [83] H. Zhao and W. Wang.  Fourier heat conduction as a strong kinetic effect in one-dimensional hard-core gases.  Phys. Rev. E, 97:010103, (2018).
  * [84] O. Narayan and S. Ramaswamy.  Anomalous heat conduction in one-dimensional momentum-conserving systems.  Phys. Rev. Lett., 89:200601, (2002).
  * [85] H. Van Beijeren.  Exact results for anomalous transport in one-dimensional hamiltonian systems.  Phys. Rev. Lett., 108:180601, (2012).
  * [86] C. B. Mendl and H. Spohn.  Dynamic correlators of Fermi-Pasta-Ulam chains and nonlinear fluctuating hydrodynamics.  Phys. Rev. Lett., 111:230601, (2013).
  * [87] H. Spohn.  Nonlinear fluctuating hydrodynamics for anharmonic chains.  J. Stat. Phys., 154:1191–1227, (2014).
  * [88] C. B. Mendl and H. Spohn.  Equilibrium time-correlation functions for one-dimensional hard-point systems.  Phys. Rev. E, 90:012147, (2014).
  * [89] S. G. Das, A. Dhar, K. Saito, C. B. Mendl, and H. Spohn.  Numerical test of hydrodynamic fluctuation theory in the Fermi-Pasta-Ulam chain.  Phys. Rev. E, 90:012124, (2014).
  * [90] C. B. Mendl and H. Spohn.  Shocks, rarefaction waves, and current fluctuations for anharmonic chains.  J. Stat. Phys., 166:841, (2017).
  * [91] E. F. Toro.  Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction.  Springer, Berlin, (2009).
  * [92] Y. Du, H. Li, and L. P. Kadanoff.  Breakdown of hydrodynamics in a one-dimensional system of inelastic particles.  Phys. Rev. Lett., 74:1268–1271, (1995).
  * [93] P. I. Hurtado.  Breakdown of hydrodynamics in a simple one-dimensional fluid.  Phys. Rev. Lett., 96:010601, (2006).
  * [94] S. Chakraborti, A. Dhar, and P.L. Krapivsky.  A splash in a one-dimensional cold gas.  arXiv:2203.11097, (2022).
  * [95] C. Cercignani, A. Frezzotti, and P. Grosfils.  The structure of an infinitely strong shock wave.  Phys. Fluids, 11(9):2757–2764, (1999).
  * [96] S. Takata, K. Aoki, and C. Cercignani.  The velocity distribution function in an infinitely strong shock wave.  Phys. Fluids, 12(8):2116–2127, (2000).
  * [97] M. Colangeli.  From kinetic models to hydrodynamics: some novel results.  Springer, (2013).
  * [98] P. L. Krapivsky, S. Redner, and E. Ben-Naim.  A kinetic view of statistical physics.  Cambridge University Press, (2010).
  * [99] M Rigol, V Dunjko, and M Olshanii.  Thermalization and its mechanism for generic isolated quantum systems.  Nature, 452(7189):854–858, (2008).
  * [100] O Howell, P Weinberg, D Sels, A Polkovnikov, and M Bukov.  Asymptotic prethermalization in periodically driven classical spin chains.  Phys. Rev. Lett., 122(1):010602, (2019).
  * [101] A. Kumar and R. Rajesh.  Blast waves in two and three dimensions: Euler versus Navier-Stokes equations.  arxiv, 2111.09213, (2021).
  * [102] H. Spohn.  Hydrodynamic equations for the Toda lattice.  arXiv:2101.06528, (2021).
  * [103] B. Doyon.  Generalized hydrodynamics of the classical Toda system.  Journal of Mathematical Physics, 60(7):073302, (2019).
  * [104] C. B. Mendl and H. Spohn.  High-low pressure domain wall for the classical toda lattice.  arXiv preprint arXiv:2011.11008, (2020).
  * [105] K. Saito, M. Hongo, A. Dhar, and S. Sasa.  Microscopic theory of the fluctuating hydrodynamics in nonlinear lattices.  Phys. Rev. Lett., 127:010601, (2021).
  * [106] V. I. Arnol’d.  Mathematical methods of classical mechanics, volume 60.  Springer Science & Business Media, (2013).
  * [107] B. V. Chirikov.  A universal instability of many-dimensional oscillator systems.  Physics reports, 52(5):263–379, (1979).
  * [108] A. J. Lichtenberg and M. A. Lieberman.  Regular and chaotic dynamics, volume 38.  Springer Science and Business Media, (2013).
  * [109] J. L. Lebowitz.  Macroscopic laws, microscopic dynamics, time’s arrow and Boltzmann’s entropy.  Physica A: Statistical Mechanics and its Applications, 194(1-4):1–27, (1993).
  * [110] J. L. Lebowitz.  Boltzmann’s entropy and time’s arrow.  Physics today, 46:32, (1993).
  * [111] J. L. Lebowitz.  From time-symmetric microscopic dynamics to time-asymmetric macroscopic behavior: An overview.  Boltzmann’s legacy, pages 63–88, (2007).
  * [112] S. Goldstein, J. L. Lebowitz, R. Tumulka, and N. Zanghì.  Gibbs and Boltzmann entropy in classical and quantum mechanics.  In Statistical Mechanics and Scientific Explanation: Determinism, Indeterminism and Laws of Nature, pages 519–581. World Scientific, (2020).
  * [113] S. Chakraborti, A. Dhar, S. Goldstein, A. Kundu, and J. L. Lebowitz.  Entropy growth during free expansion of an ideal gas.  arXiv:2109.07742, (2021).