### Thermalization, Chaos and Hydrodynamics
in Classical Hamiltonian Systems

A thesis

Submitted to the

Tata Institute of Fundamental Research, Mumbai

for the degree of Doctor of Philosophy

in Physics

by

Santhosh Ganapa

International Centre for Theoretical Sciences

Tata Institute of Fundamental Research

Bangalore, India

September 13 2022

### Declaration

This thesis is a presentation of my original research work. Wherever
contributions of others are involved, every effort is made to indicate this
clearly, with due reference to the literature, and acknowledgement of
collaborative research and discussions.

The work was done under the guidance of Prof. Abhishek Dhar and Prof. Amit
Apte, at International Centre for Theoretical Sciences, Tata Institute of
Fundamental Research (ICTS-TIFR), Bangalore.

Santhosh Ganapa

In my capacity as supervisor of the candidate’s thesis, I certify that the
above statements are true to the best of my knowledge.

Abhishek Dhar Amit Apte

Date: September 13 2022

To my parents, sister, Dhidhi and Mamma.

### Acknowledgements

I thank God Almighty for giving me the strength and the endurance to pursue my
passion: physics. In particular, i would like to thank God for giving me an
opportunity to meet and work with my master Abhishek Dhar. I am very fortunate
to have him as my role model and my advisor. He must be credited for
converting a nobody to a somebody.

I would like to thank my co-advisor Amit Apte for introducing me to the
beautiful area of Hamiltonian chaos, which eventually lead me to work in
thermalization and hydrodynamics. I benefited from having Sriram Ramaswamy and
Samriddhi Sankar Ray in my thesis monitoring committee, who clarified my
confusions and also allowed me to express myself freely. I thank them for
their friendliness. I acknowledge Anupam Kundu and Vishal Vasan for teaching
me some of the tricks of the trade. I will forever be indebted to my almae
matres BITS Pilani and ICTS for helping me realize what i want to do in life.

I rejoice the time spent with my friends Kohinoor Ghosh and Hemanta Kumar, who
made my time during PhD even more memorable. I enjoyed discussions with Varun
Dubey and many of the other talented scholars of ICTS, from whom i learnt how
to do physics. My time in ICTS was made unforgettable by my kung fu and violin
teachers. I thank my football friends for helping me control my weight and
stay healthy during my PhD. Some of this work was done during the Covid-19
pandemic and i must credit the cooks, security guards, maintenance workers and
the cleaning ladies of ICTS for doing their part during the lockdown, which
made it easy for me to concentrate only on my work. I thank my friend Bhanu
Teja for always challenging me to achieve higher in every aspect of my life.

And finally, I dedicate this work to my family. I marvel at the hard work and
dedication of my parents. I also thank them for always giving me the freedom
to choose what i want to do in life, and my sister Pavani for always
counseling me after i fall into some sort of trouble by making some stupid
choice. I thank my grandmothers Dhidhi and Mamma for their sacrifices, and
looking after me when i needed them the most.

“The mystery of life is not a problem to be solved but a reality to be
experienced.”

…Alan Watts

### List of Publications

  1. 1.

Ganapa, S., Apte, A. and Dhar, A. Thermalization of Local Observables in the
$\alpha$-FPUT Chain. J Stat Phys 180 (1), 1010–1030 Springer (2020).
https://doi.org/10.1007/s10955-020-02576-2

  2. 2.

Ganapa, S., Chakraborti, S., Krapivsky, P., L., and Dhar, A. (2021). Blast in
the One-Dimensional Cold Gas: Comparison of Microscopic Simulations with
Hydrodynamic Predictions. Phys. Fluids 33 (8).
https://aip.scitation.org/doi/10.1063/5.0058152

  3. 3.

Chakraborti, S., Ganapa, S., Krapivsky, P., L., and Dhar, A. (2021). Blast in
a One-Dimensional Cold Gas: From Newtonian Dynamics to Hydrodynamics. Phy.
Rev. Lett. 126, 244503.
https://journals.aps.org/prl/pdf/10.1103/PhysRevLett.126.244503

### List of Collaborators

  1. –

Abhishek Dhar, Professor, ICTS-TIFR, Ind.

  2. –

Amit Apte, Professor, ICTS-TIFR and IISER Pune, Ind.

  3. –

Pavel. L. Krapivsky, Professor, Boston University, USA and Skolkovo Institute
of Science and Technology, Moscow, Rus.

  4. –

Subhadip Chakraborti, postdoctoral fellow, ICTS-TIFR (former) and
Friedrich–Alexander University Erlangen-Nuremberg, Deu (current).

### List of Abbreviations

  1. –

FPUT chain: Fermi-Pasta-Ulam-Tsingou chain (formerly Fermi-Pasta-Ulam chain).

  2. –

TvNS solution: Taylor-von Neumann-Sedov solution.

  3. –

1D: one dimensional.

  4. –

AHP gas: Alternate mass hard particle gas.

  5. –

NSF equations: Navier-Stokes-Fourier equations.

  6. –

SLE: Space localized excitations.

  7. –

NMLE: Normal mode localized excitations.

  8. –

BBGKY hierarchy: Bogoliubov–Born–Green–Kirkwood–Yvon hierarchy.

  9. –

KdV equation: Korteweg–de Vries equation.

  10. –

ODE: Ordinary differential equations.

  11. –

PDE: Partial differential equations.

  12. –

LE: Local equilibrium.

  13. –

MD simulations: Molecular dynamics simulations.

###### Contents

  1. 1 Introduction
    1. 1.1 A General Introduction to this Work
      1. 1.1.1 Thermalization and its Importance in Statistical Physics
      2. 1.1.2 Evolution Before Thermalization - Description in Terms of Conserved Fields
      3. 1.1.3 Structure of the Thesis
    2. 1.2 Problem Statements
      1. 1.2.1 Thermalization in the $\alpha-$FPUT Chain
      2. 1.2.2 Blast Wave in an Alternate Mass Hard Particle Gas
    3. 1.3 Methods Used in this Study
      1. 1.3.1 Symplectic Integrator
      2. 1.3.2 Equipartition Theorem
      3. 1.3.3 Averaging Methods
      4. 1.3.4 Estimation of Equilibration Time
      5. 1.3.5 Navier-Stokes-Fourier (NSF) Equations and the MacCormack Method
  2. 2 The Fermi-Pasta-Ulam-Tsingou Story
    1. 2.1 Continuum Limit and the KdV Equation
    2. 2.2 Overlap of Resonances
    3. 2.3 Role of Breather Solutions
    4. 2.4 Ideas from Wave Turbulence
  3. 3 Thermalization of Local Observables in the $\alpha-$FPUT Chain
    1. 3.1 Initial Conditions and the Observables that are Considered in this Study
    2. 3.2 Evolution of Local Observables from Space Localized Initial Conditions
    3. 3.3 Normal Mode Localized Initial Conditions
    4. 3.4 Dependence of Equilibration on the Averaging Procedure
    5. 3.5 Quantification of the Equilibration Time
    6. 3.6 Dependence of Equilibration on the Choice of Initial Conditions and Observables
    7. 3.7 Comparison with Other Integrable Models
    8. 3.8 Dependence of Equilibration on the Width of the Initial Distribution
    9. 3.9 Relation to Chaos
    10. 3.10 Summary
  4. 4 The Taylor$-$von Neumann$-$Sedov Blast Wave Solution
    1. 4.1 Euler Equations and its Self-similar Solutions
    2. 4.2 The TvNS Blast Wave Solution in One Dimension
    3. 4.3 Earlier Attempts to Observe TvNS Scaling in Microscopic Models
    4. 4.A Rankine-Hugoniot Conditions
    5. 4.B Other Scalings of the Euler Equations
  5. 5 Blast Wave in the One Dimensional Cold Gas
    1. 5.1 The Alternate Mass Hard Particle Gas
    2. 5.2 Earlier Studies done on the AHP gas
    3. 5.3 “Blast Wave” Initial Conditions for Microscopic Simulations
    4. 5.4 Microscopic Simulations and Comparison with the TvNS Solution
    5. 5.5 Various Aspects of Microscopic Simulations
      1. 5.5.1 Comparison with the Empirical Fields
      2. 5.5.2 Dependence on Initial Conditions
      3. 5.5.3 Dependence on the Initial Ensemble
      4. 5.5.4 Validity of Local Thermal Equilibrium
      5. 5.5.5 Structure of the Shock
    6. 5.6 Modelling by the Navier-Stokes-Fourier (NSF) Equations
    7. 5.7 Scaling near the Core of the Blast
    8. 5.8 Summary
    9. 5.A Derivation of Scaling Functions near the Core
  6. 6 Conclusion
    1. 6.1 Outcomes of the Study
    2. 6.2 Limitations of the Study
    3. 6.3 Future Directions
      1. 6.3.1 Thermalization and Chaos
      2. 6.3.2 Hydrodynamics
      3. 6.3.3 Concluding Remarks

##  Chapter 1 Introduction

Hamiltonian systems are ubiquitous in physics and attempts to understand their
behaviour in a better way have been made for a long time. Consider a one
dimensional system of $N$ particles with positions and momenta given by
$\\{q_{i},p_{i}\\}$, for $i=0,1,\ldots,N-1$. Let $H(\\{p_{i},q_{i}\\})$ be the
Hamiltonian of the system. The Hamilton’s equations of motion relate the
positions and momenta as:

$\begin{split}\dot{p}_{i}=-\frac{\partial H}{\partial q_{i}},\ \
\dot{q}_{i}=\frac{\partial H}{\partial p_{i}}.\end{split}$ (1.1)

Solving the above $2N$ first order differential equations for a given initial
condition describes the behaviour of the system at subsequent times. So,
describing a Hamiltonian system is now reduced to solving a set of
differential equations with some initial condition. A Hamiltonian system with
$2N$ degrees of freedom is said to be integrable if there exists $N$
independent constants of motion $I_{1},I_{2},...I_{N}$ which Poison commute
with each other: $\\{I_{i},I_{j}\\}=0\ $, $\forall\ i,j$. Then we have a class
of systems whose description is possible in principle, by solving the
equations of motion analytically. But integrable systems are rare, and we
encounter non-integrable systems that only have a few conserved quantities in
general, in which case it is not possible to solve the equations of motion Eq.
(1.1) analytically, and we have to resort to numerical methods. Now let’s say
we are interested in describing the behaviour of a macroscopic non-integrable
system. If we try to solve Eq. (1.1), then there would be a huge number of
equations to solve, which may hide the necessary information that would have
otherwise been sufficient to describe the system. So, we are not interested in
the behaviour of individual particles. We have to move on from solving the set
of Hamilton’s equations to something else. We use the methods of statistical
physics to describe emergent behaviour in physical systems when we are dealing
with such macroscopic systems that have a huge number of particles. We
describe the system as a whole now instead of its constituent particles.

###  1.1 A General Introduction to this Work

####  1.1.1 Thermalization and its Importance in Statistical Physics

There are many interesting problems in statistical physics that one can study
in the context of classical Hamiltonian systems. One of the fundamental
questions is to understand why statistical physics works, which has been a
long-standing problem. The theory of equilibrium statistical physics says that
the observed values of physical observables, for a system in thermal
equilibrium, can be obtained by taking an average over an appropriate
equilibrium statistical ensemble. A justification of why this works is usually
made by invoking the ergodic hypothesis, which says that generic Hamiltonian
systems would have few conserved quantities, for example the total energy, and
the motion of the system in phase space would fill the energy surface densely.
Given that Hamiltonian dynamics is area preserving, this means that the
infinite time average of a physical observable would equal the microcanonical
ensemble average. However, we note that typical measurements involve time
averages over a period which is minute compared to the time required to sample
the energy surface of a typical many-body system. In other words, we observe
thermal equilibrium at time scales much shorter than the time required for the
ergodic hypothesis to hold. This is especially true for macroscopic systems,
where the time taken for an initial point in the phase space to cover most of
it is astronomical. Hence, ergodic hypothesis cannot be a sufficient basis for
the foundations of statistical physics. Then the question is to understand
thermal equilibrium itself and to know what leads a system to thermal
equilibrium so that we are able to use statistical methods after much smaller
times. Some causes that can be explored are chaos, non-integrability,
nonlinearity and mixing. One can also study the role of coarse-graining, the
role of initial conditions and the choice of macrovariables used to estimate
equilibration. One of the first studies on this problem was done by Fermi,
Pasta, Ulam and Tsingou (FPUT) [1, 2]. They expected that nonlinearity should
be sufficient for a system to thermalize but their results did not indicate
that. They instead observed a quasiperiodic behaviour in the system. Since the
original work of FPUT, there have been a number of studies aimed at
understanding their results, which led to significant developments in
statistical physics, nonlinear dynamics and mathematical physics. To this day,
there remain some unanswered questions. Thus, understanding the requirements
for the validity of equilibrium statistical physics remains a challenge. At
the same time, with the drastic improvement of computational performance, we
are now in a better position to probe these questions even further.

####  1.1.2 Evolution Before Thermalization - Description in Terms of
Conserved Fields

Another interesting problem would be to study how a system thermalizes. This
question is different from what causes the system to thermalize, which is
addressed in the previous paragraph. Here, starting from a typical initial
condition, one aims to understand how the system evolves to the thermal state.
This can be done by studying emergent phenomena in the system - we would like
to find ways to actually simplify the description of a macroscopic Hamiltonian
system. One of the standard ways to do this is to consider a coarse-grained
description of the system in terms of the conserved fields, which are the slow
degrees of freedom. This is the purview of hydrodynamics. This is the realm of
nonequilibrium statistical physics, in contrast to equilibrium statistical
physics, where there are steady state values of the conserved quantities. One
of the interesting ways to study hydrodynamics is to initially excite the
system to a large amount of energy concentrated within a small region and then
study its evolution. This leads to the formation of shocks in the system. This
so called blast problem was first studied during the Second World War by
Taylor [3, 4], von Neumann [5] and Sedov [6, 7] (TvNS) in order to understand
atomic explosions. The study of the nonequilibrium regime of the system is now
reduced to studying the dynamics of shock propagation in the system. Because a
large amount of energy is initially excited in a small region, this is an
extreme example where we can test the predictions of hydrodynamics. One such
prediction is a self-similar structure of the shock wave (in time). It is then
possible to derive the exact solution of the evolution of conserved
quantities, which helps us understand the nonequilibrium regime. It will be
very interesting if one can observe these scaling solutions in the microscopic
simulations of a many body Hamiltonian system. However, an agreement between
molecular dynamics simulations of a macroscopic Hamiltonian system and
hydrodynamics remains an open question.

####  1.1.3 Structure of the Thesis

This work is aimed to address the problems motivated by the above-mentioned
line of thinking. As the title suggests, we study various aspects of
thermalization, chaos and hydrodynamics in one dimensional Hamiltonian
systems. We study two problems. In the first problem, we revisit the problem
of equilibration in the $\alpha-$Fermi-Pasta-Ulam-Tsingou (FPUT) system. In
this system we try to understand what causes a system to thermalize. The
second problem deals with the evolution of a blast wave in an alternate mass
hard particle (AHP) gas. Here we study the nonequilibrium regime in terms of
shock propagation in the system. We aim to describe the system using
hydrodynamic equations up to times the shock wave reaches the boundary of the
system. This thesis is organized as follows:

In this chapter we next go on to explain the problem statements and also the
methods that are used in this study. Chapters 2 and 3 discuss the first
problem mentioned in the previous paragraph. In chapter 2, we give a brief
literature review of the Fermi-Pasta-Ulam-Tsingou (FPUT) problem. There is a
lot of literature on this problem, some of which are mentioned in [8, 9, 10].
We emphasize only on those aspects which are relevant to our study. Our work
is motivated by the results of [11], and so we will explain the results of
that work in detail. Then in chapter 3, we present the results of our work
[12], and also compare with the results of [11]. Chapters 4 and 5 discuss the
second problem. In the chapter 4, we derive the Taylor-von Neumann-Sedov
(TvNS) solution of the Euler equations in one dimension, which will be used
later to compare our results [13, 14] of the molecular dynamics simulations of
the alternate mass hard particle (AHP) gas. We also give a literature review
of the works which aimed at observing the TvNS scaling in molecular dynamics
simulations. In chapter 5, we present our results of the molecular dynamics
simulations, which we compare with the exact TvNS scaling solution and also
with the numerical solution of the Navier-Stokes-Fourier (NSF) equations. In
chapter 6, we conclude our work and present some interesting extensions. For
completeness sake, we also mention the limitations of our study.

###  1.2 Problem Statements

####  1.2.1 Thermalization in the $\alpha-$FPUT Chain

Consider a harmonic chain of $N$ particles, each of mass $m$ connected by a
spring of spring constant $\mu$. The Hamiltonian of this system is then
described as a function of the position $q_{i}$ and momentum $p_{i}$ of each
particle as:

$H(\boldsymbol{p},\boldsymbol{q})=\sum_{i=0}^{N-1}\left[\frac{p_{i}^{2}}{2m}+\frac{\mu(q_{i+1}-q_{i})^{2}}{2}\right]~{},$
(1.2)

where we have considered periodic boundary conditions with $q_{N}\equiv q_{0}$
and $q_{-1}\equiv q_{N-1}$, and from now on we take $m=\mu=1$ unless otherwise
mentioned. One can first make a canonical transformation from $p,q$ variables
to the normal mode variables $P,Q$ defined by:

$Q_{k}=\frac{1}{N}\sum_{j=0}^{N-1}q_{j}e^{-i2\pi
kj/N},~{}~{}P_{k}=\frac{1}{N}\sum_{j=0}^{N-1}p_{j}e^{-i2\pi kj/N}.$ (1.3)

The Hamiltonian then becomes

$\displaystyle H=N\sum_{k=0}^{N-1}E_{k},~{}~{}{\rm
where}~{}E_{k}=\Bigg{[}\frac{|P_{k}|^{2}}{2}+\frac{\Omega_{k}^{2}|Q_{k}|^{2}}{2}\Bigg{]}~{}$
(1.4)

is the energy of each mode and $\Omega_{k}=2\sin(k\pi/N),k=0,1,2,...,N-1$ is
the normal mode frequency of that particular mode. The zero mode corresponding
to $k=0$ does not participate in the dynamics, and we only deal with initial
conditions that have $E_{0}=0$ (zero initial momentum). From Hamilton’s
equations of motion in these canonical coordinates we can observe that these
normal modes decouple. The evolution of each normal mode can then be found
independently of the others. So, energy exchange does not take place between
different normal modes. So, most of the phase space is not accessible to a
harmonic chain and hence the system does not thermalize. So, we would like the
normal modes to somehow interact with each other in order to expect
thermalization. For that we add a cubic interaction term to the Hamiltonian.
If we now ask whether the system thermalizes, then the answer is not very
obvious. We then get the $\alpha-$FPUT system described by the Hamiltonian:

$H(\boldsymbol{p},\boldsymbol{q})=\sum_{i=0}^{N-1}\left[\frac{p_{i}^{2}}{2m}+\frac{\mu(q_{i+1}-q_{i})^{2}}{2}+\frac{\alpha(q_{i+1}-q_{i})^{3}}{3}\right]~{}.$
(1.5)

The origin of this Hamiltonian can be understood as arising from a spring that
has restoring force $F$ dependent on the displacement of the particle $x$
connected to it as $F=-\mu x-\alpha x^{2}$. Thus, $\alpha$ is the origin of
nonlinearity in the system. A heuristic estimate of the strength of the
nonlinearity can be obtained by comparing the contribution of the nonlinear
interaction part to the total energy. Roughly, if $r$ is the average spacing
between particles, we expect $\mu r^{2}\sim E/(N-1)=e$ which gives a length
scale $r\sim\sqrt{e/\mu}$. An estimate of the ratio of the nonlinear and
harmonic energies is then given by the parameter

$\epsilon=\frac{\alpha r^{3}}{\mu r^{2}}=\frac{\alpha e^{1/2}}{\mu^{3/2}}~{}.$
(1.6)

This dimensionless number, $\epsilon$, and the system size, $N$, are the only
relevant parameters. In our subsequent discussions we assume that $\mu=1$ and
one can change $\epsilon=\alpha e^{1/2}$ by either changing the nonlinearity
strength $\alpha$ or equivalently, by changing the energy density $e$. Note
that the cubic potential of the $\alpha$-FPUT system implies that the system
stays bounded only if the total energy is sufficiently small and the precise
condition is $E<\mu^{3}/(6\alpha^{2})$, corresponding to all energy
contributing to the potential energy of a single particle. In practice this is
highly improbable and one can work with energies slightly higher than this
bound.

Because of the presence of the cubic term, the normal modes do not decouple
and energy exchange between different normal modes takes place and one might
expect that the entire constant energy surface
$H(\boldsymbol{p},\boldsymbol{q})=E$ becomes accessible. In the original
problem [1], the system was started in a highly atypical initial condition
(all energy in the first normal mode, $E_{1}=0.08$) with $\alpha=0.25$ (where
however the fixed boundary condition case was studied). It was asked if the
system evolves to a state where all the normal modes have the same energy. To
the surprise of the authors, they did not find such a state in the time scales
of their observations. Instead, they found quasi-periodic behaviour and near-
recurrences to the initial state as shown in Fig. 1.1.

Figure 1.1: $\alpha-$FPUT chain: Graph from the original paper [1] showing the
quasiperiodic evolution of the system with near recurrences to the initial
condition when the first mode is excited. The units for energy are arbitrary.

The original paper also tried with quartic nonlinearity (the $\beta-$FPUT
system) but no thermalization was observed. The problem of whether the FPUT
system thermalizes has been studied ever since and even today there remain
some open problems in this field. Recently, [11] was able to answer positively
that the $\alpha-$FPUT system thermalizes for an arbitrary (nonzero) value of
$\alpha$. Using the methods of wave turbulence [15, 16], they derived that the
time scale for equipartition $\tau_{eq}$ in the $\alpha-$FPUT system scaled
with the dimensionless nonlinearity parameter $\epsilon$ as: $\tau_{eq}\sim
1/\epsilon^{8}$. They also verified this scaling numerically for the system by
initially exciting the first mode and then checking that the energy indeed
gets divided between all the normal modes equally when some sort of averaging
is performed.

The current work, motivated by the results of [11], looks at some aspects of
thermalization in the $\alpha-$FPUT system from a somewhat different
perspective. The problems that are addressed, and the differences in the
approaches from the previous studies are listed below:

  1. –

Dependence of equilibration on the initial conditions: Most studies look at
initial conditions specified in the normal mode space - they excite the first
or the first few normal modes and then study the evolution of the
$\alpha-$FPUT chain. However, we study the evolution by using initial
conditions that are localized in the real space - we excite a few particles
initially and then study thermalization.

  2. –

Dependence on the choice of the observable used to measure equilibration: In
most of the literature, the analysis of equipartition was done for the energy
of normal modes of the corresponding integrable problem, the harmonic chain.
The analysis is in some sense “global”, since normal modes involve all the
particles in the chain. The present work attempts to analyse the problem
“locally” - by checking equipartition theorem at different sites in the chain
and tries to compute the time scale $\tau_{eq}$ the system needs to reach
equipartition. Our method would be especially relevant when we consider strong
nonlinearity, for which the normal mode picture becomes invalid. Several
earlier works have discussed this in the context of harmonic chains [17, 18,
19] and anharmonic chains [20]. We will then compare the dependence of
$\tau_{eq}$ on the nonlinearity parameter $\epsilon$ of both the observables.

  3. –

Role of averaging in determing equilibration: We also numerically investigate
the role played by averaging over an initial ensemble in determining whether
we observe thermalization in our system, and compare the results with the time
averaging protocol.

  4. –

Relation between thermalization and chaos in the system: Finally, we
investigate the question of dependence of the equilibration process on the
width of the initial distribution. By doing so, we aim to find a relation
between equilibration of local observables and sensitive dependence of the
$\alpha-$FPUT system on initial conditions. We also compare this relation of
the $\alpha-$FPUT system to the harmonic chain and the integrable Toda chain
(the potential $V$ depends on $r$ as $V(r)=g(e^{br}-br)/b$. The parameter
choice $b=2\alpha$ and $g=b^{-1}$ would then approximate to the $\alpha-$FPUT
potential upto leading nonlinearity).

####  1.2.2 Blast Wave in an Alternate Mass Hard Particle Gas

The evolution of a blast wave emanating from an intense explosion has drawn
much interest since the Second World War. The rapid release of a large amount
of energy in a localized region produces a shock wave propagating into the non
perturbed gas. The dynamics of the system can then be studied in terms of
shock propagation in the system. Taylor [3, 4], von Neumann [5] and Sedov [6,
7] (TvNS) understood that the compressible gas dynamics [21, 22] provides an
appropriate framework and described this problem using the evolution of
conserved quantities of the gas system: the density field $\rho(r,t)$, the
velocity field $v(r,t)$ and the energy field $e(r,t)$. They also found an
exact solution as long as the pressure in front of the shock wave can be
neglected compared to the pressure behind the shock wave. In this situation,
the hydrodynamic equations have a scaling solution and the hydrodynamic
variables acquire a self-similar form (in time) that are functions only of the
scaling variable. The hydrodynamic equations are thus converted from partial
differential equations (PDE) to ordinary differential equations (ODE)
depending only on this scaling variable, and they admit an exact solution.

It will be an interesting study to look for these scalings in the microscopic
simulations of a many-particle system. Such studies have been attempted only
recently, which can be partly attributed to the recent increase in
computational performance. As mentioned in the introduction, finding an
agreement between molecular dynamics simulations of a Hamiltonian system and
hydrodynamics remains an open question. Studies done in two and three
dimensions [23, 24, 25, 26] showed deviations from hydrodynamics. This was
never done in one dimensional systems.

This work attempts to understand the blast wave dynamics in a one dimensional
gas of hard point particles. In this system, the only interactions between the
particles is through point collisions between the nearest neighbors, which are
perfectly elastic and instantaneous. In between collisions the particles move
ballistically with constant velocities. Consider two colliding particles with
masses $(m_{1},m_{2})$ and initial velocities $(v_{1},v_{2})$. Their post-
collisional velocities are given by:

$\begin{split}v_{1}^{\prime}=\frac{2m_{2}v_{2}+(m_{1}-m_{2})v_{1}}{m_{1}+m_{2}},\
\
v_{2}^{\prime}=\frac{2m_{1}v_{1}+(m_{2}-m_{1})v_{2}}{m_{1}+m_{2}}.\end{split}$
(1.7)

These follow from the conservation of momentum and energy. If $m_{1}=m_{2}$,
then one can notice from the above equations that after the collision,
$v_{1}^{\prime}=v_{2},v_{2}^{\prime}=v_{1}$, i.e. the particles merely
exchange velocities. This means that the system effectively behaves as a non-
interacting system. This is an integrable system, whose evolution cannot be
described by the usual hydrodynamics. We choose a system that only has a few
conserved quantities (e.g. particle number, energy and momentum). This is the
case with typical non-integrable systems, in which case it is possible to
describe the evolution using the usual hydrodynamics. Thus, we take a non-
integrable system, where all odd numbered particles are taken to have mass
$m_{1}$ and all even numbered particles are taken to have mass $m_{2}$($\neq
m_{1}$). This system only has three conserved quantities, and so we expect the
usual hydrodynamics to be sufficient to describe it. We excite the system to a
sudden blast wave $-$ energy concentrated in a small region somewhere in the
middle of the system, and evolve the system only upto times the energy reaches
the boundary. The problems that this work addresses, are listed below:

  1. –

The TvNS solution of the 1D Euler equations for an ideal gas: We aim to derive
the expression for the position of the shock front in terms of the initial
total energy of the blast $E$ and the ambient density of the system
$\rho_{\infty}$. We also derive the expressions for scaling solution of the
conserved quantities. Together, this is the complete TvNS solution in one
dimension, which is an exact scaling solution of the Euler equations with an
ideal gas equation of state.

  2. –

Comparison of microscopic simulations of the 1D alternate mass hard particle
gas (AHP) and the TvNS solution: We first compare the results of molecular
dynamics of the alternate mass hard particle gas with the TvNS solution. We
also study the subtleties of the microscopic simulations like comparison of
ensemble averaging and spatial coarse graining protocols, and also the
comparison of the microcanonical and canonical ensemble averaging.

  3. –

Modelling using the NSF equations: In order to understand some of the
discrepancies between the microscopics and the TvNS solution, we then consider
the evolution of the conserved fields of AHP system with the Navier-Stokes
Fourier (NSF) equations (Euler equations with added dissipative terms due to
viscosity and heat conduction). We then move on to studying the various
scaling solutions that are possible for the hydrodynamic equations and show
TvNS scaling as one of the scaling solutions, and that it specifically results
from a blast wave initial condition.

  4. –

Study of deviations from the TvNS scaling and the cause of deviation: We
attribute the discrepancies between microscopic simulations and the TvNS
solution to the non-negligible contribution of the dissipation terms in the
NSF equations. We then perform a detailed analysis in this regime, where these
terms are not neglected and try to explain the deviation of the microscopics
from the TvNS scaling.

###  1.3 Methods Used in this Study

####  1.3.1 Symplectic Integrator

Simulations of the $\alpha-$FPUT system were done by using a $6^{th}$ order
symplectic integrator [27]. This is more accurate than the Runge-Kutta
methods. The reason for this is that the Runge-Kutta methods treat position
and momentum of each degree of freedom as independent of each other at each
instant of time. We note that an evolution in momentum at a certain time step
will automatically affect its corresponding canonical position through
$\dot{q}=p$/m. The Runge-Kutta methods do not take this interdependence into
consideration and thus introduce additional error in the solution, which is
avoided in symplectic integrators upto the desired level of accuracy. The
time-step size was taken as $0.01$ and we checked the relative energy change
in the system at the end of the computation to be of the order of $10^{-11}$.
Computations of the Lyapunov exponent were performed by solving the coupled
system of $2N+2N$ nonlinear and linearized equations using a fourth order
Runge-Kutta integrator. In this case the time-step size was taken as $0.001$
and the relative energy change in the system at the end of the computation is
around $10^{-9}$. We have checked that the results do not change significantly
on decreasing the time step by a factor of two.

####  1.3.2 Equipartition Theorem

We estimate equilibration time by checking Equipartition theorem at different
sites. According to the theorem, for a system in thermal equilibrium, the
following is true:

$\left\langle q_{i}\frac{\partial H}{\partial
q_{j}}\right\rangle_{eq}=\left\langle p_{i}\frac{\partial H}{\partial
p_{j}}\right\rangle_{eq}=c\delta_{ij},$ (1.8)

for all $i,j=0,1,\ldots,N-1$ and where $\left\langle...\right\rangle_{eq}$
represents an average over an equilibrium ensemble (the different types of
averaging are described below) and $c$ is a constant equal to $k_{B}(\partial
S/\partial E)^{-1}$ for the microcanonical ensemble and equal to $k_{B}T$ for
a canonical ensemble. Our method is in contrast to checking Equipartition
among normal modes because normal modes of the harmonic chain only
approximately describe the $\alpha-$FPUT system, while Eq. (1.8) is exact for
the latter system and can be used for a wide range of parameters and for
different systems. Choosing initial conditions localized in real space has the
distinct advantage that all the members of the ensemble used in the ensemble
averaging protocol lie on the same energy surface and hence, we will be able
to link thermalization and chaos. This cannot be done by initially exciting
the normal modes because all the members of the ensemble would then not lie on
the same energy surface - we cannot fix the normal mode energies and the total
energies at the same time.

####  1.3.3 Averaging Methods

One can use various averaging methods to check whether a system has reached
thermal equilibrium. It is clear that checking for equilibration during the
time evolution requires some kind of averaging, and we discuss three possible
approaches:

#####  Averaging over Initial Conditions

This approach involves performing an averaging over initial conditions [11],
with the expectation that such averages would represent the typical behaviour.
In our work, we perform an average over the initial conditions $({\bf
q}_{0},{\bf p}_{0})$ which are now chosen from a narrow distribution centred
around a specified point and that are still on the microcanonical surface of
constant energy $E$, momentum $P$ and number of particles $N$. Denoting the
initial distribution by $\rho_{I}({\bf q}_{0},{\bf p}_{0})$, we then obtain
the following average for any observable $A=A({\bf q},{\bf p})$:

$\langle A\rangle(t)=\int d{\bf q}_{0}d{\bf p}_{0}~{}A({\bf q}(t),{\bf
p}(t))\rho_{I}({\bf q}_{0},{\bf p}_{0})~{}.$ (1.9)

In our simulations we generate $R$ initial configurations from this
distribution and evolve each of them with Hamiltonian dynamics. The ensemble
average of a physical observable $A({\bf q},{\bf p})$ is estimated from Eq.
(1.9) as:

$\langle A\rangle(t)=\frac{\sum_{r=1}^{R}A_{r}(t)}{R},$

where the sum is over the $R$ members of the ensemble.

##### Temporal Coarse Graining

A second approach would be again to start with a fixed initial condition and
performing an averaging over time [28, 20, 29]. In this case, one can look at
microscopic variables, e.g the kinetic energy of individual particles, and ask
whether the time averaged value corresponds to the expected equilibrium value.
Starting the time evolution of the system from an arbitrary initial condition
$({\bf q}_{0},{\bf p}_{0})$, with energy $E$, we can define a time averaged
quantity for the observable $A=A({\bf q},{\bf p})$ as

$\bar{A}(t)=\frac{1}{t}\int_{0}^{t}dsA({\bf q}(s),{\bf p}(s))~{}.$ (1.10)

For a non-integrable system we would expect that for generic initial
conditions, at long times we should get thermal equilibration or
$\bar{A}(t\to\infty)=\langle A\rangle$. The time to reach the equilibrium
value should give a measure of equilibration time scales.

#####  Spatial Coarse Graining

We now discuss a protocol which is perhaps most relevant from the physical
point of view if one remembers that real observations are made on single
systems and on short time scales. It then makes sense to consider the system
being in a single microscopic configuration, but we now measure coarse grained
physical observables, involving some spatial averaging [30, 31]. In this
approach one should eventually look at large systems and the ideas related to
typicality and law of large numbers are relevant for understanding
equilibration. As an example, consider the problem of free expansion of a gas
on removal of a partition inside a box. We take as our observable to be the
number of particles in a tiny box, that is still large enough to contain
thousands of particles, and located in the initially empty half of the box. We
expect that the number of particles will increase with time to eventually
equilibrate around an average value with small fluctuations.

For the $\alpha-$FPUT system, we use the ensemble averaging protocol in order
to estimate the equilibration time. As mentioned before, doing so would give
us a way to relate equilibration of local observables to chaos. We briefly
compare the results of ensemble averaging with those of the time averaging
protocol. We also use the ensemble averaging protocol in order to study the
blast propagation in the alternate mass hard particle gas. Here, we compare
these results with those of spatial averaging. We will describe the spatial
coarse graining protocol in chapter $5$.

####  1.3.4 Estimation of Equilibration Time

As a measure of the level of equipartition that is achieved, we define the
following function, which has been referred to in the literature as entropy:

$S(t)=-\sum_{i=0}^{N-1}f_{i}(t)\ln f_{i}(t),$ (1.11)

where $f_{i}(t)=\nu_{i}(t)/\sum_{r=0}^{N-1}\nu_{r}(t)$ for $i=0,1,\dots,N-1$,
and $\nu_{i}$ could be either $\langle p_{i}\frac{\partial H}{\partial
p_{i}}\rangle$ or $\langle q_{i}\frac{\partial H}{\partial q_{i}}\rangle$ or
$\langle E_{i}\rangle$. The value of $S(t)$ is bounded between $0$,
corresponding to the highly nonequilibrium situation with all the energy in a
single degree of freedom, and $\ln N$, corresponding to the equilibrated
system with equipartition between all degrees with $\left\\{f_{i}\right\\}$
defining a uniform distribution over the set $\left\\{0,1,\dots,N-1\right\\}$.
Since it can theoretically take an infinite amount of time to reach $\ln N$,
we estimate the equilibration time as the time required for $S(t)$ to reach a
predetermined value of entropy which is close to the equilibration value. In
this work, we consider the following criteria to determine the equilibration
timescale:

$\left|\frac{S(t)-S_{max}}{S_{max}}\right|\leq 0.01~{}.$ (1.12)

The above threshold must be satisfied for two consecutive values of the time
that is sampled. The minimum value of $t$ for which the above criteria is
satisfied is termed as equilibration time (denoted by $\tau_{\rm eq}$). We use
this to quantify the equilibration time in the $\alpha-$FPUT system.

####  1.3.5 Navier-Stokes-Fourier (NSF) Equations and the MacCormack Method

The molecular dynamics simulations of the alternate mass hard particle gas
were done by an event-driven algorithm, which is much faster than the time-
driven algorithm. The only known conserved quantities of this system are the
total number of particles, the total momentum and the total energy of the
system, and we expect a hydrodynamic description in terms of the corresponding
conserved fields namely, the mass density field $\rho(x,t)$, the momentum
density field $\rho(x,t)v(x,t)$ ($v$ is the velocity field) and the energy
density field $E(x,t)$. We compare microscopic simulations with the results
obtained by solving the one-dimensional Navier-Stokes-Fourier (NSF) equations,

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

which are supplemented by the thermodynamic relations:
$T(x,t)=2\mu\,\epsilon(x,t)$ and
$p(x,t)=2\rho(x,t)\epsilon(x,t)=\mu^{-1}\rho(x,t)\,T(x,t).$ Here
$\epsilon(x,t)=e(x,t)-\frac{v(x,t)^{2}}{2}$ is the internal energy per unit
mass, $\mu=(m_{1}+m_{2})/2$, and we have set the Boltzmann’s constant
$k_{B}=1$. The dissipative effects are manifested in the bulk viscosity
$\zeta$ and the thermal conductivity $\kappa$ appearing in Eqs. (1.13).
Consistent with the Green-Kubo formalism and earlier studies, in our numerical
study we have 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 taken
to be unity.

The Navier-Stokes-Fourier (NSF) equations were solved using the MacCormack
method [32, 33], which is an explicit two-step, second order, finite
difference method. As an illustration of the method, consider the following
first order hyperbolic equation [34]:

$\frac{\partial u}{\partial t}+a\frac{\partial u}{\partial x}=0.$ (1.14)

The initial profile $u(x,0)$ is known, and we aim to find $u(x,t)$ by using
this method. The application of MacCormack method to the above equation
proceeds in two steps; a predictor step which is followed by a corrector step.
The predictor step is obtained by replacing the spatial and temporal
derivatives in Eq. (1.14) using forward differences. The output of this step
is the provisional value of $u$ at time level $n+1$ for each value of the
space level $i$ denoted by $\overline{u_{i}^{n+1}}$, which is given by:

$\overline{u_{i}^{n+1}}=u_{i}^{n}-a\frac{\Delta t}{\Delta
x}(u_{i+1}^{n}-u_{i}^{n}).$ (1.15)

In the corrector step, the predicted value $\overline{u_{i}^{n+1}}$ is
corrected by using backward finite difference approximation for the spatial
derivative. Also the time step used in the corrector step is $\Delta t/2$, in
contrast to the $\Delta t$ used in the predictor step. This results in the
desired output value $u_{i}^{n+1}$, which is given by:

$u_{i}^{n+1}=\frac{u_{i}^{n}+\overline{u_{i}^{n+1}}}{2}-a\frac{\Delta
t}{2\Delta x}(\overline{u_{i}^{n+1}}-\overline{u_{i-1}^{n+1}}).$ (1.16)

The order of differencing can be reversed for the time step (i.e.,
forward/backward followed by backward/forward). In our work, we chose the
former order, and we have taken $\Delta t=0.001$ and $\Delta x=0.1$.

##  Chapter 2 The Fermi-Pasta-Ulam-Tsingou Story

Ever since the original numerical experiment, there have been numerous studies
that tried to explain the recurrent behaviour in the Fermi-Pasta-Ulam-Tsingou
system. It all started with Fermi, Pasta, Ulam and Tsingou wanting to check if
nonlinearity is enough for a system to thermalize. But the failure to observe
thermal behaviour actually led to significant developments in statistical
physics, nonlinear dynamics and mathematical physics in the coming years.
There were some basic questions regarding the problem that needed answers. One
of them was to understand if there existed a stochasticity threshold in the
FPUT system - the critical value of the nonlinearity parameter below which the
system is not expected to thermalize. Maybe the parameter taken in the
original problem was lower than this threshold and hence thermalization was
not observed. May be the system needed to be simulated for much longer times
so that we could observe thermalization. Then there must be a theory which
explains what is that time scale. There is also a possibility of the existence
of breather solutions - time-periodic and space-localized stable solutions in
the FPUT system that could delay or even prevent thermalization.

In this chapter we discuss some of the existing literature on this problem.
FPUT problem has a vast literature, and here we only discuss some of the
literature and touch upon only those aspects which are relevant to the problem
that this work tries to address. We only refer to several of the review
articles on the FPUT problem [8, 9, 10, 28, 35]. Some of the recent studies
that focus on the specific aspects of the FPUT problem that are especially
relevant to this paper include the role of breathers [20, 36, 37, 38, 39],
wave-wave interaction theory [11, 40, 41], and the role of breakup of
invariant tori and the stochastic threshold [42, 43, 44, 45]. Our emphasis is
on [11], which explains the equilibration problem in terms of resonant
interactions between the nonlinear waves, and also tries to estimate the time
scale the system needs in order to attain equipartition.

###  2.1 Continuum Limit and the KdV Equation

Zabusky and Kruskal [46] showed that the continuum limit of the FPUT model
leads to the Korteweg–de Vries (KdV) equation, which is an integrable model
and has soliton solutions. Two or more solitons do not interact irreversibly,
and they just “pass through” each other virtually unaffected in size or shape.
Because the solitons are remarkably stable entities, preserving their identity
through numerous interactions, these authors expected the FPUT system to
exhibit thermalization (complete energy sharing among the corresponding linear
normal modes) only after extremely long times, if ever. One might then expect
the long wavelength initial condition used by FPUT to remain close to the
continuum description for long times and so there should not be any
thermalization. This is indeed consistent with what is observed in the
original problem, which was however done for N = 32. We will come back to the
KdV equation later when we discuss breather solutions.

###  2.2 Overlap of Resonances

As mentioned in the beginning of this chapter, one important aspect that was
not understood at the time the FPUT problem was first posed, was if there
existed a stochasticity threshold - critical value of the nonlinearity
parameter below which the system does not thermalize at all, and above which
there is stochasticity in the system, which would eventually lead to
thermalization. Izrailev and Chirikov [42] proposed that there exists a
stochasticity threshold value of the energy density ($e=E/N$) beyond which one
expects thermalization. This idea was developed using the criterion of
overlapping resonances. At the onset of stochasticity, the separation between
the resonances ($\Delta_{\omega}\equiv\omega_{k+1}-\omega_{k}$) is of the
order of resonance width of each resonance, arising due to the nonlinearity.
Here $\omega_{k}$ is the frequency of the $k^{th}$ mode of the harmonic
oscillator. It can be shown that for small values of $k\ (k<<N)$,
$\Delta_{\omega}\approx 2\pi/N$, while for larger values of $k\ (N-k<<N)$,
$\Delta_{\omega}\approx\pi^{2}/2N^{2}$, which is much less than that of low
frequency modes. So, the separation between the neighbouring frequencies is
smaller for larger values of k. If the system is initially excited in one
particular normal mode such that its frequency shift arising due to
nonlinearity is much less than $\Delta_{\omega}$, then one can neglect the
influence of neighbouring resonances, since there is no overlap. For smaller
values of k, the separation between the resonances is high and hence the
stochasticity threshold needed for equipartition is also high. The initial
condition taken by FPUT is below this stochastic border and hence there is a
recurrent behaviour. On the other hand, if a higher mode is excited initially,
then the separation between the resonances is relatively small. With the
stochasticity threshold relatively small, one can easily observe irregular
motion with strong energy sharing among a large number of high frequency
modes.

The stochastic threshold has also been studied numerically [43, 44, 45]. It is
numerically found that for generic initial conditions [43], the time evolution
of the maximum Lyapunov exponent $\Lambda$ is indistinguishable for the
$\alpha$ -FPUT and the Toda systems upto a characteristic time (called the
trapping time $\tau_{\rm tr}$) that increases with decreasing energy (at fixed
N). Then suddenly a bifurcation appears, which can be discussed in relation to
the breakup of regular, solitonlike structures. After this bifurcation,
$\Lambda$ of the $\alpha$ -FPUT model appears to approach a constant, while
$\Lambda$ of the Toda system appears to approach 0, consistent with its
integrability. This difference is attributed to the untrapping of the FPUT
system from its regular region in phase space and escaping to the chaotic
component of its phase space. They infer that once a trajectory enters the
chaotic component of phase space, equipartition will eventually be attained.
The authors computed the equilibration time $\tau_{\rm eq},\tau_{\rm tr}$ and
$\Lambda$ for system sizes $N=32,64,128$ at different energy densities and
found the existence of a threshold $e_{c}(N)$ such that for $e<e_{c}$,
$\tau_{\rm eq}$ and $\tau_{\rm tr}$ seemed to diverge while $\Lambda$
vanishes. The threshold decreases with system size as $e_{c}(N)\sim 1/N^{2}$.
Above $e_{c}$, power law dependences of the form $\tau_{\rm eq}\approx
1/e^{3},\tau_{\rm tr}\approx 1/e^{2.5}$ and $\Lambda\approx e^{2}$ were noted.
The FPUT parameters correspond to $e\ll e_{c}$.

###  2.3 Role of Breather Solutions

Breathers are time-periodic and space-localized solutions of nonlinear
dynamical systems [36]. By nature, they are non-thermal and one might expect
them to play a role in preventing thermalization. In a certain sense the idea
is similar to the one relating the presence of solitons in the KdV system to
the absence of equilibration in the FPUT system—the difference being that
breathers are stable solutions of the discrete system while the KdV is a
continuum approximation. The unexpected recurrences in the FPUT problem have
been linked to the choice of initial conditions used by FPUT, which are set
close to exact coherent time-periodic (or even quasiperiodic) trajectories,
e.g., q-breathers, which show exponential localization of energy in normal
mode space. Even if these trajectories have a measure zero in the phase space,
they might have a finite measure impact simply by being linearly stable [37].
Several other studies admit coherent time-periodic states localized in real
space, which are known as discrete breathers or intrinsic localized modes.

In [20], Danieli and others have shown how one can check for the presence of
breather solutions. They called trapping time, as the time the trajectory
spends off the equilibrium manifold before piercing it again. By measuring the
probability distribution functions of the trapping times in the FPUT chain,
they infer that the trajectory is with high probability getting trapped in
some parts of the phase space for long times. They also conjecture that these
trapping events are due to visiting regions of the phase space which are
substantially close to some regular orbits, which are supported by numerical
experiments.

###  2.4 Ideas from Wave Turbulence

The FPUT problem has recently been studied [11, 40, 41] using the approach of
wave turbulence [15, 16]. Based on requirement of resonance between sets of
normal modes, a detailed prediction has been made for the time-scale for
equilibration and it is argued that this is finite for any non-zero strength
of non-linearity, for finite sized systems. We now describe the work [11] in
some detail. Let us first describe the $\alpha-$FPUT system in terms of the
normal mode coordinates $a_{k}$:

$a_{k}=\frac{1}{\sqrt{2\omega_{k}}}(P_{k}-i\omega_{k}Q_{k}).$ (2.1)

If $\alpha\neq 0$, there will be coupling among the normal modes and let us
write the equations of motion. Writing in terms of the dimensionless
variables,

$a_{k}^{\prime}=\frac{(\mu/m)^{1/4}}{\sqrt{\sum_{k}\omega_{k}\mid
a_{k}(t=0)\mid^{2}}}a_{k},\ t^{\prime}=\sqrt{\frac{\mu}{m}}t,\
\omega_{k}^{\prime}=\sqrt{\frac{m}{\mu}}\omega_{k},$

we can write the equations of motion after removing the primes for brevity:

$i\frac{\partial a_{1}}{\partial
t}=\omega_{1}a_{1}+\epsilon\sum_{k2,k3}V_{1,2,3}(a_{2}a_{3}\delta_{1,2+3}+2a_{2}^{\star}a_{3}\delta_{1,3-2}+a_{2}^{\star}a_{3}^{\star}\delta_{1,-2-3})~{},$
(2.2)

where the matrix $V_{1,2,3}$ weights the transfer of energy between wave
numbers $k_{1}$, $k_{2}$ and $k_{3}$ and is given by:

$V_{1,2,3}=-\frac{1}{2\sqrt{2}}\frac{\sqrt{\omega_{1}\omega_{2}\omega_{3}}}{sign(sin(\frac{\pi
k_{1}}{N})sin(\frac{\pi k_{2}}{N})sin(\frac{\pi k_{3}}{N}))}~{}.$ (2.3)

The dimensionless parameter $\epsilon$ is given by:

$\epsilon=\frac{\alpha}{m}(\frac{\mu}{m})^{1/4}\sqrt{\sum_{k}\omega_{k}\mid
a_{k}(t=0)\mid^{2}}~{},$ (2.4)

where $a_{i}\equiv a(k_{1},t)$ and the delta function is 1 also if the
argument differs mod N. The variables $k_{2}$ and $k_{3}$ run from $-N/2+1$ to
$N/2$. It is easy to notice that this definition of $\epsilon$ reduces to the
same as that we have used in the previous chapter. Eq. (2.2) describes the
evolution of the normal mode variables in the $\alpha$-FPUT system and as
expected, it has coupling between the different modes. We are now interested
in the weak nonlinearity limit $\epsilon<<1$. The nonlinear behaviour of the
$\alpha$-FPUT system is due to three-wave interactions, with the strength of
coupling determined by $\epsilon$. First we try to find out if these three
waves interactions are resonant. The way to check for resonances is to first
do a canonical transformation and then try to remove the three wave
interactions. If a canonical transformation exists (does not diverge) then we
say that the three wave interactions are not resonant. Otherwise, we conclude
that there are resonances and these interactions cannot be removed by a
canonical transformation. For the $\alpha$-FPUT system, this canonical
transformation, in terms of $b$ variables, up to first order is given by:

$a_{1}=b_{1}+\epsilon\sum_{k2,k3}(A_{1,2,3}^{(1)}b_{2}b_{3}\delta_{1,2+3}+A_{1,2,3}^{(2)}b_{2}^{\star}b_{3}\delta_{1,3-2}+A_{1,2,3}^{(3)}b_{2}^{\star}b_{3}^{\star}\delta_{1,-2-3})+O(\epsilon^{2})~{}.$
(2.5)

Here,

$A_{1,2,3}^{(1)}=V_{1,2,3}/(\omega_{3}+\omega_{2}-\omega_{1}),$
$A_{1,2,3}^{(2)}=2V_{1,2,3}/(\omega_{3}-\omega_{2}-\omega_{1}),$
$A_{1,2,3}^{(3)}=V_{1,2,3}/(-\omega_{3}-\omega_{2}-\omega_{1}).$

It is easy to verify using the trigonometric identities that for the
frequencies $\\{\omega_{i}\\}$ of the harmonic chain the denominators in the
above transformations are never zero. Therefore, a canonical transformation
exists, the three wave interactions are not resonant and hence we observe the
quasiperiodic behaviour at short time scales. Substituting this transformation
in the equations of motion, we get:

$i\frac{\partial b_{1}}{\partial
t}=\omega_{1}b_{1}+\epsilon^{2}\sum_{k2,k3,k4}T_{1,2,3,4}b_{2}^{\star}b_{3}b_{4}\delta_{1+2,3+4}+O(\epsilon^{3})~{}.$
(2.6)

Now proceeding in a similar way, we try to remove these four wave interactions
using a suitable canonical transformation. This will be possible if these
interactions are not resonant. The above equation also has terms including the
Kronecker deltas $\delta_{1,3+4+2}$,$\delta_{1,-3-4-2}$ and
$\delta_{1,4-3-2}$. However, those terms are not resonant and can be removed
by higher order terms in the transformation, and are of no use in this
analysis. It is easy to verify that the resonances of Eq. (2.6) are described
by the equations:

$k_{1}+k_{2}-k_{3}-k_{4}=0\ mod\ N,\ \rm{and}\
\omega_{1}+\omega_{2}-\omega_{3}-\omega_{4}=0~{}.$ (2.7)

Now there are two kinds of solutions for this equation: (i) Trivial solutions:
These are obtained when $k_{1}=k_{3}$ and $k_{2}=k_{4}$, or $k_{1}=k_{4}$ and
$k_{2}=k_{3}$. These trivial solutions are responsible for a nonlinear
frequency shift and do not contribute to the energy transfer between the
modes. (ii) Nontrivial solutions: These are obtained when
$(k_{1},k_{2},k_{3},k_{4})=(k_{1},k_{2},-k_{1},-k_{2})$ or $k_{1}+k_{2}=mN/2$,
$m=0,\pm 1,\pm 2,...$. These solutions allow for a transfer of energy between
the the modes. However, these solutions are isolated from one another. For
example, two sets of solutions $(k_{1},k_{2},k_{3},k_{4})=(1,15,-1,-15)$ and
$(k_{1},k_{2},k_{3},k_{4})=(3,13,-3,-13)$ have no wave number in common. Hence
energy transfer can only take place between the modes that satisfy these
conditions and energy transfer does not take place across the spectrum. Thus
we say that these resonances are isolated and this would not lead to
equipartition of energy. Nevertheless, these resonances do exist and cannot be
removed by a higher order canonical transformation. Now, we move on to higher
order resonances. Doing another canonical transformation would lead to five
wave interactions, which can be shown to be non resonant. So, these can be
removed by an appropriate choice of a canonical transformation, again just
like a three wave interaction. So, we go to the next order interactions, which
are the six wave interactions. The evolution of the system is described in
terms of the following equation:

$i\frac{\partial b_{1}}{\partial t}=\omega_{1}b_{1}+\epsilon^{2}\sum
T_{1,2,3,4}b_{2}^{\star}b_{3}b_{4}\delta_{1+2,3+4}+\epsilon^{4}\sum
W_{1,2,3}^{4,5,6}b_{2}^{\star}b_{3}^{\star}b_{4}b_{5}b_{6}\delta_{1+2+3,4+5+6}+O(\epsilon^{5})~{}.$
(2.8)

In the above equation, the four wave interactions cannot be removed because
they are resonant and a canonical transformation suitable for removing those
modes would diverge. Resonances for this equation can be found by finding the
solutions of:

$k_{1}+k_{2}+k_{3}-k_{4}-k_{5}-k_{6}=0\ mod\ N,\ \rm{and}\
\omega_{1}+\omega_{2}+\omega_{3}-\omega_{4}-\omega_{5}-\omega_{6}=0~{}.$ (2.9)

There are three kinds of solutions for this equation now: (i) Trivial
solutions: As before, they would only lead to a nonlinear frequency shift.
(ii) Nontrivial symmetric resonances: These are given by
$(k_{1},k_{2},k_{3},k_{4},k_{5},k_{6})=(k_{1},k_{2},k_{3},-k_{1},-k_{2},-k_{3})$
or $k_{1}+k_{2}+k_{3}=mN/2$, $m=0,\pm 1,\pm 2,...$. For example,
$(k_{1},k_{2},k_{3},k_{4},k_{5},k_{6})=(1,2,13,-1,-2,-13)$ and
$(k_{1},k_{2},k_{3},k_{4},k_{5},k_{6})=(1,8,7,-1,-8,-7)$ have the wave number
$k=1$ in common. (iii) Nontrivial quasisymmetric resonances: These are given
by
$(k_{1},k_{2},k_{3},k_{4},k_{5},k_{6})=(k_{1},k_{2},k_{3},-k_{1},-k_{2},k_{3})$
or $k_{1}+k_{2}=mN/2$, $m=0,\pm 1,\pm 2,...$. For example,
$(k_{1},k_{2},k_{3},k_{4},k_{5},k_{6})=(1,15,13,-1,-15,13)$ is connected with
$(k_{1},k_{2},k_{3},k_{4},k_{5},k_{6})=(13,3,7,-13,-3,7)$ through $k=13$ mode.
One can also find a common wave number between nontrivial symmetric resonances
e.g., $(k_{1},k_{2},k_{3},k_{4},k_{5},k_{6})=(2,4,10,-2,-4,-10)$ and
nontrivial quasisymmetric resonances e.g.,
$(k_{1},k_{2},k_{3},k_{4},k_{5},k_{6})=(7,9,10,-7,-9,10)$. Thus, these
resonances are interconnected and they lead to the equipartition of energy. By
using BBGKY hierarchy (as given in [11]), one can estimate that the
equilibration time ($\tau_{\rm eq}$) of the $\alpha-$FPUT system depends on
the nonlinearity parameter as:

$\tau_{\rm eq}\propto\frac{1}{\epsilon^{8}}~{}.$ (2.10)

Had the four wave resonances been interconnected, as is the case in the
thermodynamic limit, the system would thermalize sooner. In [11] it is argued
that in the thermodynamic limit the dependence would instead have been:
$\tau_{\rm eq}\sim 1/\epsilon^{4}$. A more recent study [40] of the
$\beta$-FPUT chain suggests that wave-wave resonances lead to thermalization
at small $\epsilon$, where $\tau_{\rm eq}\sim 1/\epsilon^{4}$, while at larger
$\epsilon$ the overlap of resonances (Sec. 2.2) leads to $\tau_{\rm eq}\sim
1/\epsilon$.

##  Chapter 3 Thermalization of Local Observables in the $\alpha-$FPUT Chain

In this chapter, we study the problem of equilibration of the $\alpha-$FPUT
system using a different approach. As mentioned in the problem statement in
chapter 1, our work attempts to analyse the problem “locally”- by checking
equipartition theorem at different sites in the chain instead of checking the
normal mode energy which is “global”, since normal modes involve all the
particles in the chain. Also, most studies look at initial conditions
specified in the normal mode space - they excite the first or the first few
normal modes and then study the evolution of the FPUT chain. However, we study
the evolution by using initial conditions that are localized in the real space
- we excite a few particles initially and then study thermalization. We study
various aspects of thermalization like the dependence of the equilibration
time on the choice of initial conditions, choice of observables, choice of
averaging and the role of chaos in equilibration, and also compare our results
with the predictions of the wave turbulence theory [11].

###  3.1 Initial Conditions and the Observables that are Considered in this
Study

Let us first describe the initial conditions that we have used in our study.
Here we consider the case where the initial distribution lies on the constant
$(E,P,N)$ surface in $2N$ dimensional phase space and has a small spread, with
the size of the spread characterized by a dimensionless number $\gamma$. The
initial condition is also chosen to correspond to the initial energy of the
chain being localized in a small region called space localized excitations
(SLE). We set $q_{i}=0$, for $i=1,2,\ldots,N$ and $p_{i}=0$ for
$i=5,6\ldots,N-1,N$. The total energy $E$ is then distributed amongst the four
remaining particles in the following way:

$\displaystyle E_{1}$ $\displaystyle=$ $\displaystyle
E_{2}=(1-\gamma)E/4+\nu\gamma E/2,~{}$ (3.1) $\displaystyle E_{3}$
$\displaystyle=$ $\displaystyle E_{4}=(1-\gamma)E/4+(1-\nu)\gamma E/2,~{}$
(3.2)

where $0<\gamma<1$ is a number which specifies the “width” of the distribution
and $\nu$ is a uniformly distributed random number in the interval $(0,1)$
which basically gives us some randomness in the initial conditions. $\gamma=0$
corresponds to a fixed initial state while $\gamma=1$ corresponds to the
broadest distribution. Here we consider initial conditions that have zero
momentum. The first and the second particles are given velocities in opposite
directions, as are the third and the fourth. We shall refer to these initial
conditions as space localized excitations (SLE) as opposed to normal mode
localized excitations (NMLE) commonly used in most studies, where the first or
the first few normal modes are originally excited. In our simulations we
generate $R$ initial configurations from this distribution and evolve each
with Hamiltonian dynamics. We then consider the ensemble average of a physical
observable $A({\bf q},{\bf p})$ estimated as

$\langle A\rangle(t)=\frac{\sum_{r=1}^{R}A_{r}(t)}{R},$

where the sum is over the $R$ members of the ensemble. Instead of normal
modes, the observables that we are interested in are given by the equiparition
theorem. In the rest of this chapter we use the following notation:

$\displaystyle\left\langle T_{i}\right\rangle=\frac{1}{2}\left\langle
p_{i}\frac{\partial H}{\partial
p_{i}}\right\rangle\,\quad\textrm{and}\quad\left\langle
V_{i}\right\rangle=\left\langle q_{i}\frac{\partial H}{\partial
q_{i}}\right\rangle\,.$ (3.3)

$\left\langle T_{i}\right\rangle$ is nothing but the average kinetic energy of
the $i^{th}$ particle. Note however that $\left\langle V_{i}\right\rangle$ is
not the average potential energy in general.

###  3.2 Evolution of Local Observables from Space Localized Initial
Conditions

$(a)$$(b)$

Figure 3.1: $\alpha-$FPUT chain: Plot shows time evolution of $\left\langle
T_{j}\right\rangle$ (left panel) and $\left\langle V_{j}\right\rangle$ (right
panel) starting from space localized initial conditions for $j=0,1,2,...N-1$.
Parameter values for this plot are $N=32$, $E=31$, $\alpha=0.0848$
($\epsilon=0.0834$), $\gamma=0.9$ and $R=1000$. The solid black line
corresponds to the equipartition value $\langle T_{j}\rangle=0.484$ and
$\langle V_{j}\rangle=0.969$. The inset shows the time evolution at the
earliest times.

$(a)$$(b)$

Figure 3.2: $\alpha-$FPUT chain: Same as Fig. 3.1 but with $\gamma=10^{-8}$.
Note that the range of the time $t$ and the averages $\langle T_{j}\rangle$
and $\langle V_{j}\rangle$ are different from those in Fig. 3.1.

$(a)$$(b)$

Figure 3.3: $\alpha-$FPUT chain: Spatial profile of $\left\langle
T_{j}\right\rangle$ at different times for $\gamma=0.9$ (left panel) and
$\gamma=10^{-8}$ (right panel) and other parameters $N=32$, $E=31$,
$\alpha=0.0848$ and $R=1000$. We see that the initially localized energy
quickly spreads through the chain while equipartition is achieved at much
longer time scales.

In Fig. 3.1 we show the time evolution of $\left\langle T_{j}\right\rangle$
and $\left\langle V_{j}\right\rangle$ for $\gamma=0.9$. We see that there is a
long transient period and then we see equipartition at times $\sim 2\times
10^{4}$. At the earliest times, the inset in Fig. 3.1 shows near-recurrent
behaviour. The results for $\gamma=10^{-8}$ are plotted in Fig. 3.2, where we
now see that equipartition is achieved at somewhat longer times, around $t\sim
10^{5}$. We will return to the question of dependence of the equilibration
time on the width of the initial distribution $\gamma$ later, when we try to
relate thermalization and chaos. In Fig. 3.3 we plot the averaged kinetic
energy profile at different times. Here it can be seen that the energy spreads
quickly through the entire system, while equipartition is achieved at much
longer time scales.

$(a)$$(b)$

Figure 3.4: $\alpha-$FPUT chain: Here we examine how the fluctuations seen in
$\left\langle T_{j}\right\rangle$ depend on the number of realizations $R$. We
plot $\left\langle T_{j}\right\rangle$ for $R=10^{2}$ (left panel) and
$R=10^{4}$ (right panel). Other parameters were taken as $N=32$, $E=31$,
$\alpha=0.0848$, $\gamma=0.9$. The insets show zoom-ins at short and long
times. A vertical line is drawn at $t=1.2\times 10^{4}$, which is the
equilibration time. Up to this time, both the plots look nearly the same.

A closer examination of the plots in Figs. 3.1, 3.2 shows that even at late
times, the averaged quantities continue to fluctuate around their equilibrium
values. We show in Fig. 3.4 that these fluctuations in fact decrease with
increase in the number of realizations $R$ used to compute averages. We also
see that $\left\langle T_{j}\right\rangle$ at pre-thermalization times does
not depend significantly on $R$ and the plots are near identical for
$R=10^{2}$ and $R=10^{4}$ (up to the vertical line in Fig. 3.4).

###  3.3 Normal Mode Localized Initial Conditions

We now discuss and compare our results with those obtained in studies on
equipartition of normal mode energies using initial conditions which are
localized in the Fourier space [11]. The authors in [11] considered initial
conditions where the energy was distributed between the modes $k=1$ and $k=31$
with frequencies $\Omega_{1}=\Omega_{31}$. Averages were done over $1000$
initial conditions by choosing random phases for the modes. The time evolution
of the normal mode energies was monitored to check for equipartition. Here we
reproduce their numerical results and compare with the results in the previous
section. We consider again $N=32$ particles with total energy $E=31$. In Fig.
3.5 we show the time evolution of the energy of all the modes in the system.
Comparing with Figs. 3.1, 3.2 it is clear that equilibration now occurs at a
time scale ($\sim 10^{6}$) that is about an order of magnitude longer.

$(a)$$(b)$

Figure 3.5: $\alpha-$FPUT chain: Left panel shows time evolution of
$\left\langle E_{k}\right\rangle$ for $N=32,E=31,\alpha=0.0848$ and $R=1000$.
Only modes $k=1,31$ are initially excited. The inset shows quasiperiodic
behaviour at short time scale. Right panel shows the mode energy profile at
different times.

###  3.4 Dependence of Equilibration on the Averaging Procedure

$(a)$$(b)$

Figure 3.6: $\alpha-$FPUT chain: Plot shows time evolution of $T_{j}$ (left
panel) and $V_{j}$ (right panel) starting from a single initial condition.
Parameter values for this plot are $N=32$, $E=31$, $\alpha=0.0848$
($\epsilon=0.0834$). The insets show zoom-ins of the time evolution at early
times and at late times, and we see oscillatory behaviour in both cases. Thus
in this case, no equipartition is achieved. Figure 3.7: $\alpha-$FPUT chain:
Plot shows time evolution of $E_{k}$ starting from a single initial condition
with energy in two normal modes. Parameter values for this plot are $N=32$,
$E=31$, $\alpha=0.0848$. The insets show zoom-ins of the time evolution at
early times and at late times and we see oscillatory behaviour in both cases.
Thus in this case, no equipartition is achieved.

To demonstrate that the averaging procedure is crucial to the equilibration
process in this set-up, we show the evolution of $T_{j}$ and $V_{j}$ for a
single initial condition. In this case, we see in Fig. 3.6 that no
equilibration is achieved and the oscillatory behaviour persists up to
$t=5\times 10^{5}$. In Fig. 3.7 we see again that one does not see any signs
of equilibration in the evolution of $E_{k}$ for a single realization.

We can also discuss thermalization using a different protocol where one starts
with a single initial condition and then considers a time average of any given
observable given by Eq. (1.10). In Fig. 3.8 we show results obtained by using
this protocol for both the space-local and normal mode observables. The insets
in the figures show that thermalization time scales are completely different
from those obtained by the ensemble averaging protocol. This shows the
important role played by the choice of averaging protocol in determining how
fast the system equilibrates.

$(a)$$(b)$

Figure 3.8: $\alpha-$FPUT chain: Plot shows time evolution of running time
averages of $T_{j}$ for SLE (left panel) and $E_{k}$ for NMLE (right panel)
starting from a single initial condition. Parameter values for this plot are
$N=32$, $E=31$, $\alpha=0.0848$ ($\epsilon=0.0834$). The insets show
comparison between time averages and ensemble averages, and illustrates that
the latter procedure leads to faster thermalization.

###  3.5 Quantification of the Equilibration Time

$(a)$$(b)$

Figure 3.9: $\alpha-$FPUT chain: Plot of $S^{SL}(t)$ corresponding to
$\left\langle T_{j}\right\rangle$ for $\gamma=0.9$ (left panel) and
$\gamma=10^{-8}$ (right panel) with other parameters
$N=32,E=31,\alpha=0.0848$. The insets show the convergence of the
equilibration process as the number of realizations is increased.

To get a more systematic and quantitative estimate of the equilibration time
we now look at the entropy function defined in Eq. (1.11). This in some sense,
performs an average over all the degrees of freedom and attains its maximum
value $\ln N$ when all degrees have equilibrated. In Fig. 3.9, we plot the
evolution of entropy for the two parameter values $\gamma=0.9$ and
$\gamma=10^{-8}$. The insets show zoom-ins near the equilibrium value $\ln N$,
showing the approach to equilibration and its dependence on the number of
realizations. For higher number of realizations the fluctuations in the
entropy are lower and also the mean is closer to the equilibrium value.

We use the criterion of Eq. (1.12) to estimate the equilibration time from the
entropy $S^{SL}$ corresponding to $\langle T_{j}\rangle$ and find $\tau_{\rm
eq}\approx 100300$ and $12500$ for $\gamma=10^{-8}$ and $0.9$ respectively. In
general we find that as the “width” of the distribution $\gamma$ is increased,
thermalization happens faster. We will discuss this again later.

###  3.6 Dependence of Equilibration on the Choice of Initial Conditions and
Observables

Using our method, we compute the equilibration time for different values of
the dimensionless parameter $\epsilon$ for $N=32$. These results are plotted
in Fig. 3.10 where we find a power-law dependence of equilibration time,
$\tau_{\rm eq}$, on $\epsilon$ of the form

$\tau_{\rm eq}\propto\frac{1}{\epsilon^{a}}~{}.$ (3.4)

The value $a$ is found to depend on $\gamma$ and lies between 4 and 6. This is
significantly different from the form 1/$\epsilon^{8}$ obtained in [11], by
considering equilibration of normal modes. In Fig. 3.10 we also indicate the
relaxation time results for the normal modes which give $a\approx 7.7$.

It is to be expected that the equilibration time scale should depend not only
on the initial ensemble in which the system is prepared, but also on the
observable for which equipartition is being tested. We investigate this
question further by computing the entropy functions $S^{SL}$ and $S^{NML}$ for
both types of initial conditions, namely space localized (SLE) and normal mode
localized (NMLE). These results have been plotted in Fig. 3.11. We see clearly
that the relaxation of normal mode coordinates is slower than that of the
space localized observables, irrespective of initial distribution. This is
consistent with the higher value of the exponent $a$.

Figure 3.10: $\alpha-$FPUT chain: Graph showing $\tau_{\rm eq}$ of
$\left\langle T_{j}\right\rangle$ for $N=32,E=31$ as a function of $\epsilon$.
The slopes of the fitting lines give $a=5.9,6.0,4.0$ for $\gamma=10^{-8},0.01$
and $0.9$ respectively. We also show the equilibration times obtained from the
normal mode entropy function $S^{NML}$ for normal mode localized initial
conditions, which leads to an exponent $a\approx 7.7$.

$(a)$$(b)$

Figure 3.11: $\alpha-$FPUT chain: Graphs showing the evolution of the entropy
corresponding to space-local observables and normal modes, for $\gamma=0.9$
(left panel) and $\gamma=10^{-8}$ (right panel), with other parameters given
by $N=32,E=31,\alpha=0.0848$. Results for both space localized initial
conditions (SLE) and normal mode localized initial conditions (NMLE) are
shown.

###  3.7 Comparison with Other Integrable Models

To investigate the role of integrability, we now repeat the above computations
in two integrable models that are related to the $\alpha-$FPUT system in the
limit of weak nonlinearity. We consider the harmonic chain which is described
by the Hamiltonian in Eq. (1.2) and the Toda chain, described by the
Hamiltonian

$H(\boldsymbol{p},\boldsymbol{q})=\sum_{i=0}^{N-1}\left[\frac{p_{i}^{2}}{2}+\frac{g}{b}e^{b(q_{i+1}-q_{i})}\right]~{}.$
(3.5)

The Toda system is known to be integrable [47, 48] and has been much studied
as the integrable limit of the FPUT chain [28, 43, 49, 50]. The parameter
choice $b=2\alpha$ and $g=b^{-1}$ would then approximate the $\alpha-$FPUT
potential to leading nonlinearity. Starting with the same space-localized
initial conditions as in the previous sections, we now check equipartition of
kinetic energy $T_{i}$. In Fig. 3.12 we see that no equilibration is achieved
for the harmonic chain. On the other hand, somewhat surprisingly, we see in
Fig. 3.13 that the Toda chain does equilibrate, provided we start with a wider
initial distribution ($\gamma=0.9$). In Fig. 3.14 we plot the entropy function
$S^{SL}$ and using the criterion in Eq. (1.12), estimate the equilibration
time and find $\tau_{\rm eq}\approx 66000$ for $\gamma=0.9$.

$(a)$$(b)$

Figure 3.12: Harmonic chain: Plot of $\left\langle T_{j}\right\rangle$ as a
function of time, for $N=32,E=31,R=1000$ and $\gamma=10^{-8}$ (left panel) and
$\gamma=0.9$ (right panel). In both cases there is no sign of thermalization.

$(a)$$(b)$

Figure 3.13: Toda chain: Plot of $\left\langle T_{j}\right\rangle$ as a
function of time, for $N=32,E=31,\alpha=0.0848,b=2\alpha,g=b^{-1},R=1000$ and
$\gamma=10^{-8}$ (left panel), and $\gamma=0.9$ (right panel). Now we observe
thermalization in the right panel.

$(a)$$(b)$

Figure 3.14: Toda chain: Plot of the entropy $S^{SL}(t)$ for $\left\langle
T_{j}\right\rangle$ of the Toda chain with same parameters as Fig. 3.13 and
$\gamma=10^{-8}$ (left panel), and $\gamma=0.9$ (right panel). The second case
shows clear equilibration.

###  3.8 Dependence of Equilibration on the Width of the Initial Distribution

We now discuss the dependence of the equilibration process in the
$\alpha-$FPUT system on the width of the initial distribution $\gamma$ used to
compute the ensemble averages. For the $\alpha-$FPUT system, we now take
$N=32,E=31,R=1000$ and $\alpha=0.0848$. We compute the equilibration time
$\tau_{\rm eq}$ for different values of $\gamma$. We find that broader the
initial distribution (higher the $\gamma$), faster the equilibration. More
precisely, we observe that $\tau_{\rm eq}\propto log(\gamma)$. In Fig. 3.15a,
the solid line describes an ensemble of initial conditions described by Eqs.
(3.1)-(3.2), referred to as zero volume ensemble, since the latter occupies
zero volume in the phase space. The magnitude of the slope of this line on a
semilog plot is found to be $4.5\times 10^{3}$. Its inverse is $2.2\times
10^{-4}$, close to the Lyapunov exponent of the system $\Lambda\approx
4.1\times 10^{-4}$. We have also studied an ensemble of initial conditions
that has randomness in all the degrees of freedom (while maintaining momentum
conservation), thereby occupying a finite volume in the phase space,
quantified by the number $\gamma$. The results for this are shown by the
dashed line in Fig. 3.15a where again we see the logarithmic dependence on
$\gamma$ and in fact find a closer agreement between the magnitude of the
slope (the inverse of which is $3.3\times 10^{-4}$) and the Lyapunov exponent.
For the Toda system however, we observe a different behaviour. In Fig. 3.15b
we show the dependence of $\tau_{\rm eq}$ on $\gamma$ for the Toda chain,
again with $\alpha=0.0848$. The slope on a log-log plot is close to $-1$
suggesting $\tau_{\rm eq}\propto 1/\gamma$. In the next section, we show how
this dependence links equilibration of local observables to the growth of
perturbations of initial conditions, and hence to chaos. This idea then also
explains the absence of equilibration in the harmonic chain and slow
equilibration in the Toda chain.

$(a)$$(b)$

Figure 3.15: Graphs showing the dependence of equilibration time of
$\left\langle T_{j}\right\rangle$ on $\gamma$ for the $\alpha$-FPUT chain
(left panel) and the Toda chain for $\alpha=0.0848$ (right panel). Other
parameters are $N=32$, $E=31$, and $R=1000$. For the $\alpha$-FPUT chain,
plotted on a linear-log scale, there are two different ensembles. The solid
line represents an ensemble of initial conditions occupying zero volume in the
phase space. It is described by Eqs. (3.1)-(3.2). The slope of this line is
$-4.5\times 10^{3}$. The dashed line represents an ensemble of initial
conditions occupying a finite volume in the phase space. The slope of this
line is $-3.0\times 10^{3}$. The slope for the Toda chain, plotted on a log-
log scale, is $-0.81$. Figure 3.16: The time-evolution of the ensemble of
initial conditions ($\gamma=0.9$) obtained using Eqs. (3.1)-(3.2) is shown for
the FPUT chain with $\alpha=0.0848$ (left column), the corresponding Toda
chain (middle column), and the harmonic chain (right column), in the
$(p_{13},q_{13})$ and $(p_{0},q_{0})$ planes (top and bottom rows,
respectively). The different lines for harmonic chain are at increased times,
with increasing marker size.

###  3.9 Relation to Chaos

We first quantify the growth of perturbations in the system. Let us consider
an infinitesimal perturbation $\\{\delta q_{i}(0),\delta p_{i}(0)\\}$ of an
initial condition $\\{q_{i}(0),p_{i}(0)\\}$. We compute the quantity

$\displaystyle Z(t)$ $\displaystyle=\frac{\sum_{i=1}^{N}\left[\delta
q_{i}^{2}(t)+\delta p_{i}^{2}(t)\right]}{\sum_{i=1}^{N}\left[\delta
q_{i}^{2}(0)+\delta p_{i}^{2}(0)\right]}\,.$ (3.6)

We then compute the ensemble averaged time-dependent Lyapunov exponent
$\lambda(t)$ defined as

$\displaystyle\lambda(t)=\frac{1}{2t}\langle\ln Z(t)\rangle~{},$ (3.7)

where $\langle...\rangle$ denotes an average over initial conditions
$\\{q_{i}(0),p_{i}(0)\\}$ chosen from the distribution $\rho_{0}({\bf q},{\bf
p})$. As described earlier, the numerical integration of
$\\{q_{i}(t),p_{i}(t),\delta q_{i}(t),\delta p_{i}(t)\\}$ is done by solving
$2N+2N$ nonlinear and linearized equations. The largest Lyapunov exponent
$\Lambda$ is then given by $\Lambda=\lim_{t\to\infty}\lambda(t)$. Since the
harmonic chain and the Toda chain are integrable, $\Lambda=0$ for both of
them. So, let us first understand why the Toda chain thermalizes for a broad
initial distribution and the harmonic chain never thermalizes. Since they are
both integrable, they can be written in terms of action-angle variables. The
$N$ action variables (I) are constants of motions, while the $N$ angle
variables ($\boldsymbol{\theta}$) evolve linearly with time as:

$\theta_{i}(t)=\theta_{i}(0)+\omega_{i}(\textbf{I})t~{}.$ (3.8)

The separation between any two initial conditions
($\boldsymbol{I_{1}},\boldsymbol{\theta_{1}}$) and
($\boldsymbol{I_{2}},\boldsymbol{\theta_{2}}$) evolves as:

$\theta_{1i}(t)-\theta_{2i}(t)=\theta_{1i}(0)-\theta_{2i}(0)+(\omega_{1i}(\boldsymbol{I_{1}})-\omega_{2i}(\boldsymbol{I_{2}}))t~{}.$
(3.9)

for $i=0,1,2,...N-1$. In case of the harmonic chain, the $\omega$’s are not
dependent on the action variables and are constants for different initial
conditions. Thus, the separation between the two initial conditions is
bounded. For the Toda chain however, the $\omega$’s are dependent on the
action variables and are thus different for different initial conditions. Here
the separation between the two initial conditions grows linearly with time.
Hence we expect the Toda system to thermalize for a broad enough initial
distribution.

The basic picture that illustrates the difference between the three models is
shown in Fig. 3.16, where we show the time-evolution of the ensemble of
initial conditions obtained using Eqs. (3.1)-(3.2). The FPUT chain shows a
fast growth in phase space because of its positive Lyapunov exponent, while it
takes much longer for the Toda chain because of its integrability (and linear
temporal growth of perturbations). There is no spread in the harmonic chain.

In Fig. 3.17 we plot $\langle\log Z(t)\rangle$ for the FPUT chain as well for
the corresponding harmonic chain and Toda chain (with $b=2\alpha$,
$g=b^{-1}$). We confirm the expected exponential growth of $Z(t)$ for the FPUT
chain at large times, it’s linear growth for the Toda chain and the lack of
growth in case of the harmonic chain. In the inset of the right panel, we also
show the line $\langle\log Z(t)\rangle=10.3$, which corresponds to the
equilibration time of the $\alpha-$FPUT chain ($\approx 12000$, which we found
earlier in this chapter). It turns out that the point of intersection of this
horizontal line with the Toda chain is very close to its equilibration time
($\approx 66000$). This gives us a means to relate the thermalization
properties of both the systems to the growth of their perturbations.

$(a)$$(b)$

Figure 3.17: Left panel shows $Z(t)$ of $\alpha-$FPUT chain, Toda chain and
the harmonic chain at short times. It can be seen that the growth of
perturbations is bounded for the harmonic chain and linear for the Toda chain.
Right panel shows $Z(t)$ of $\alpha-$FPUT chain and the Toda chain for long
times. The growth of perturbations in the $\alpha-$FPUT chain though linear at
short times, becomes exponential at long times. It can be seen that the growth
of perturbations is linear for the FPU chain and logarithmic for the Toda
chain in the logarithmic scale. A horizontal line drawn at $\langle\log
Z(t)\rangle=10.3$ intersects both the curves. Its significance is explained in
the text.

Let us define $D(t)=\left[\sum_{i=1}^{N}\left(\Delta q_{i}^{2}(t)+\Delta
p_{i}^{2}(t)\right)\right]^{1/2}$ to be the distance at time $t$, of points
that are initially separated by a small but _finite_ separation $D(0)$, due to
a finite perturbation of initial conditions. In particular, we will interpret
$D(t)$ as the “spread” of an ensemble of trajectories, for example, using the
ensemble of initial conditions as described in Eqs. (3.1)-(3.2). At early
times we expect that the growth can be described by

$\displaystyle D(t)\sim\begin{cases}D(0)e^{\Lambda t}&\textrm{for FPUT,}\\\
D(0)t&\textrm{for Toda.}\end{cases}$ (3.10)

We say early times because $\\{\Delta{\bf q}(t),\Delta{\bf p}(t)\\}$ cannot
grow forever since $\\{{\bf q}(t),{\bf p}(t)\\}$ are constrained to be on the
constant energy surface and are themselves bounded. The initial width $D(0)$
should be proportional to $\gamma$, which characterizes the width of our
initial phase-space distribution. Thus we have, taking $\log$ on both sides:

$\displaystyle\log(D(t))\sim\begin{cases}\Lambda t+\log(\gamma)&\textrm{for
FPUT,}\\\ \log(t)+\log(\gamma)&\textrm{for Toda.}\end{cases}$ (3.11)

Thus $D(t)$ is bounded - we expect it to attain the maximum value after
thermalization has been reached. Therefore, if we plot $\tau_{\rm eq}$ as a
function of $\gamma$ for the $\alpha-$FPUT system then we should expect a
straight line on a log-linear plot, which is what is observed in Fig. 3.15a.
As reported earlier, the magnitude of the slope of the line is close to the
inverse of the Lyapunov exponent. Thus, we have made a strong point for our
claim regarding the relation between the equilibration of local observables
and the sensitive dependence of the system on initial conditions. Similarly,
if we plot the $\tau_{\rm eq}$ as a function of $\gamma$ for the Toda system
then we should expect a straight line on a log-log plot of slope $-1$, which
is again, what is observed in Fig. 3.15b. Thus, we have explained the results
of Fig. 3.15.

The system is translationally invariant. So, the results are identical if we
initially distribute the energy to a different set of four successive
particles, while maintaining the order of the initial distribution. If the
energies are distributed in a different permutation then we find that while
the precise equilibration times are different, the dependence of $\tau_{\rm
eq}$ on $\epsilon$ and $\gamma$ are still the same.

For the $\alpha-$FPUT system, as mentioned in the previous section, 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.
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. Nevertheless, this method gives
us a way to link chaos and thermalization in the $\alpha-$FPUT system, and we
expect this method to also work for other chaotic systems.

###  3.10 Summary

We studied the time scale of thermalization of local variables in the
$\alpha$-FPUT chain and its two limiting integrable versions, namely the
harmonic chain and the Toda chain. Considering systems with $N=32$ particles
and total energy $E=31$, we estimated the thermalization time $\tau_{\rm eq}$
by measuring $\langle z_{i}\partial H/\partial z_{i}\rangle$ ($z_{i}$
indicating phase space coordinates) and finding the time to attain
equipartition. The averaging is done over initial conditions chosen from a
distribution whose width is characterized by the parameter $0\leq\gamma\leq
1$, with $\gamma=1$ corresponding to the broadest distribution and $\gamma=0$
corresponding to a fixed initial condition. The initial distribution is taken
to be one where energy is localized initially in real space instead of normal
mode space. The system is described by a single dimensionless parameter
$\epsilon=\alpha(E/N)^{1/2}$ characterizing the effective nonlinearity. Some
of our main findings are as follows:

  * –

For the $\alpha$-FPUT chain we find $\tau_{\rm eq}\propto{1}/{\epsilon^{a}}$
with $a$ between 4 and 6, and $\gamma$ dependent, in contrast to normal mode
equilibration times [11], where one finds $a\approx 8$.

  * –

We find that local variables equilibrate at much shorter time scales than
normal modes for both normal mode localized excitations (NMLE) and phase space
localized excitations (SLE).

  * –

For the $\alpha$-FPUT chain we observe the time scales for equilibration to be
completely different when we perform ensemble averaging and time averaging
protocols.

  * –

The thermalization time depends on the initial ensemble and we find $\tau_{\rm
eq}\propto\ln(\gamma)$, with the proportionality constant being close to the
inverse of the maximal Lyapunov exponent of the system, thus quantifying the
relation between thermalization and chaos for the $\alpha$-FPUT system.

  * –

Surprisingly, we find that the Toda chain equilibrated on very long time
scales if the width of the initial distribution is broad enough. In fact we
obtain $\tau_{\rm eq}\sim 1/\gamma$. On the other hand, the harmonic chain
never equilibrates.

  * –

We provide a simple geometric understanding of these results – the
equilibration time is simply related to the time it takes for an ensemble of
initial conditions in the $2N$ dimensional phase space to spread over the
microcanonical energy surface. For the FPUT chain, for energies such that the
system is chaotic with a positive Lyapunov exponent, a fast exponential (in
time) spreading occurs. For the Toda chain the growth is linear and so
thermalization takes more time. We provide numerical evidence to support this
picture.

##  Chapter 4 The Taylor$-$von Neumann$-$Sedov Blast Wave Solution

Up until now, we have studied what causes a system to thermalize and we
explored the roles of initial conditions, averaging, choice of observables and
chaos in thermalization. Now we move on to a different problem. For a one
dimensional system, instead of describing the evolution by solving $2N$
equations, can we somehow use significantly less number of equations to
describe the evolution? As usual, we are only looking at a coarse-grained
description of the system when we are trying to describe the evolution. We
therefore enter the realm of hydrodynamics, which describes the evolution of
conserved fields, which are the slow degrees of freedom. We say slow because
these fields evolve slower than the evolution of individual particles. So, we
are in the non-equilibrium regime where, we study the evolution of conserved
quantities by using the Euler equations.

In this chapter, we study the blast problem. In Sec. 4.1 we introduce the
blast problem and the Euler equations in $d$ dimensions. Now solving Euler
equations, which are PDEs is not always be easy. For the blast problem, where
the energy is initially given only to a small region, the Euler equations
however admit solutions that are self-similar (in time) up to the times the
perturbation reaches the boundary. In this case the Euler equations can be
converted into ODEs (the scaling variable would then be the only independent
variable) and we can get an exact solution. This was first noticed in three
dimensions by Taylor [3, 4], von Neumann [5] and Sedov [6, 7] (TvNS), who
studied the self-similar solutions in the context of atomic explosions. For
our purposes, we derive the exact scaling solution in one dimension in Sec.
4.2, which will be valid upto the times the shock reaches the boundary of the
system. We then cite some of the literature which tried to find an agreement
between hydrodynamics and molecular dynamics simulations of a macroscopic
Hamiltonian system in Sec. 4.3. Some technical details like derivation of the
Rankine-Hugoniot conditions and exploring what are the various scaling regimes
of the Euler equations are presented in the appendices.

###  4.1 Euler Equations and its Self-similar Solutions

Consider a system of interacting particles in which the only known conserved
quantities are the total number of particles, the total momentum and the total
energy of the system. We begin by discussing the general case of propagation
of a blast wave in a $d$ dimensional system for which, because of the radial
symmetry of the problem there are three Euler equations for the hydrodynamic
fields which depend only on the radial coordinate $r$. We therefore expect a
hydrodynamic description in terms of the corresponding conserved fields
namely, the mass density field $\rho(r,t)$, the momentum density field
$\rho(r,t)v(r,t)$ ($v$ is the velocity field) and the energy density field
$E(r,t)=\rho(r,t)e(r,t)$.

The blast is caused by the instantaneous release of energy $E$ in a very small
region that creates a spherical shock wave propagating through the quiescent
gas. The blast is infinitely strong if the pressure behind the shock wave can
be neglected. In normal conditions, this is valid up to a certain time; for
the gas at zero temperature, it remains valid forever. Below we consider the
infinitely strong blast. Dimensional analysis [51] alone gives the position of
the shock wave $R=R(t)$ in terms of time $t$, the released energy $E$, and the
background density $\rho_{\infty}$:

$R(t)=\left(\frac{Et^{2}}{A\rho_{\infty}}\right)^{\frac{1}{d+2}}.$ (4.1)

The dimensionless amplitude $A$ must be determined somehow, but the dependence
of the position of the shock front $R(t)$ on the basic parameters
automatically follows from dimensional considerations. Interestingly, whether
the shock propagation in the system is diffusive, sub-diffusive or super-
diffusive depends on the dimension $d$ the system is set up in, as can be seen
from the above equation.

The density is uniform and equals to $\rho_{\infty}$ everywhere in front of
the shock wave, i.e., for $r>R(t)$. The velocity of the shock wave is

$U=\frac{dR}{dt}=\delta\,\frac{R}{t},\qquad\delta\equiv\frac{2}{d+2}\,.$ (4.2)

Behind the shock wave $0\leq r<R(t)$, the radial velocity $v(r,t)$, density
$\rho(r,t)$ and pressure $p(r,t)$ satisfy:

$\displaystyle\partial_{t}\rho+\partial_{r}(\rho v)+\frac{d-1}{r}\,\rho v=0,$
(4.3a)
$\displaystyle(\partial_{t}+v\partial_{r})v+\frac{1}{\rho}\,\partial_{r}p=0,$
(4.3b)
$\displaystyle(\partial_{t}+v\partial_{r})\ln\frac{p}{\rho^{\gamma}}=0,$
(4.3c)

where $\gamma$ is the adiabatic index. These equations are the Euler equations
for conserved fields and we used the entropy form for the ideal gas. Instead
of energy we have used pressure, which is related to energy by the equation of
state. The Rankine-Hugoniot conditions [21], describing the jump between the
states on both sides of the shock wave are given by:

$\frac{\rho(R)}{\rho_{\infty}}=\frac{\gamma+1}{\gamma-1}\,,~{}~{}\frac{v(R)}{U}=\frac{2}{\gamma+1}\,,~{}~{}\frac{p(R)}{\rho_{\infty}U^{2}}=\frac{2}{\gamma+1}$
(4.4)

in the case of the infinitely strong blast. The pressure $p_{\infty}$ in front
of the shock wave can be neglected as long as $p(R)\gg p_{\infty}$. Using Eqs.
(4.1)–(4.2) and (4.4) one finds that this is valid in the time range

$t\ll t_{*},\quad
t_{*}\sim\left(\frac{E}{p_{\infty}}\right)^{\frac{1}{d}}\sqrt{\frac{\rho_{\infty}}{p_{\infty}}}.$
(4.5)

In the case when the particles surrounding the blast are initially at rest,
$p_{\infty}=0$, so $t_{*}=\infty$ and the blast forever remains infinitely
strong. Now we convert the PDEs given by Eqs. (4.3) to ODEs.

Instead of pressure, we can use the temperature field given by (for the ideal
gas) $T=\mu p/\rho$. Dimensional analysis ensures that the hydrodynamic
variables acquire a self-similar form

$\rho=\rho_{\infty}G(\xi),\quad v=\delta\,\frac{r}{t}\,V(\xi),\quad
T=\frac{\mu\delta^{2}}{\gamma}\,\frac{r^{2}}{t^{2}}\,Z(\xi).$ (4.6)

The factors $\delta,\gamma$ are inserted for convenience; e.g., from Eq. (4.2)
we see that the velocity of the shock wave is $\delta{R}/{t}$ and this
suggests the use of the factor $\delta$ for $v$. The fields now depend on the
single dimensionless variable

$\xi=\frac{r}{R}.$ (4.7)

One seeks to find the scaling functions $G(\xi),V(\xi)$ and $Z(\xi)$ behind
the shock wave, $0\leq\xi\leq 1$. The Rankine-Hugoniot conditions Eq. (4.4)
(also derived in Appendix 4.A) become:

$\displaystyle G(1)=\frac{\gamma+1}{\gamma-1},$ (4.8a) $\displaystyle
V(1)=\frac{2}{\gamma+1},$ (4.8b) $\displaystyle
Z(1)=\frac{2\gamma(\gamma-1)}{(\gamma+1)^{2}}.$ (4.8c)

The conservation of energy allows to express $Z$ through the scaled velocity
$V$:

$Z=\frac{\gamma(\gamma-1)(1-V)V^{2}}{2(\gamma V-1)}.$ (4.9)

This integral of motion is usually established [21] for $d=3$, but the same
derivation works in arbitrary dimension and yields the universal result Eq.
(4.9).

Plugging the ansatz Eqs. (4.6)–(4.7) into Eq. (4.3a) we obtain

$\frac{dV}{d\ell}+(V-1)\,\frac{d\ln G}{d\ell}=-dV,$ (4.10)

where $\ell=\ln\xi$. Similarly we transform Eq. (4.3c) into

$\frac{d\ln Z}{d\ell}-(\gamma-1)\,\frac{d\ln G}{d\ell}=\frac{d+2-2V}{V-1}.$
(4.11)

Equations (4.10)–(4.11) are solvable for arbitrary $d$ and $\gamma>1$. One
must solve these equations even if one merely wants to determine the amplitude
$A=A(d,\gamma)$ in Eq. (4.1). The energy conservation gives

$\displaystyle E$ $\displaystyle=$
$\displaystyle\int_{0}^{R}dr\,\,\Omega_{d}\,r^{d-1}\,\rho\left[\frac{v^{2}}{2}+\frac{c^{2}}{\gamma(\gamma-1)}\right]$
(4.12) $\displaystyle=$
$\displaystyle\rho_{\infty}\,\Omega_{d}\,\delta^{2}\,\frac{R^{d+2}}{t^{2}}\int_{0}^{1}d\xi\,\xi^{d+1}\,\frac{(\gamma-1)V^{3}}{2(\gamma
V-1)}\,G,$

where $c=\sqrt{\gamma p/\rho}$ is the sound speed, $\Omega_{d}$ is the surface
area of the unit sphere and we have used Eq. (4.9). Combining Eq. (4.12) with
Eq. (4.1) we obtain

$A=\Omega_{d}\,\delta^{2}\,\frac{\gamma-1}{2}\int_{0}^{1}d\xi\,\xi^{d+1}\,\frac{V^{3}}{\gamma
V-1}\,G.$ (4.13)

###  4.2 The TvNS Blast Wave Solution in One Dimension

In the classical literature, the blast problem is studied in three dimensions;
the two-dimensional solutions are mentioned in [7]. Here we present the
derivation of the solution in $d=1$. The problem is solvable for an arbitrary
adiabatic index $\gamma>1$, but we use $\gamma=3$ following the general
prediction, $\gamma=1+2/d$, of kinetic theory for monoatomic gases [52]. The
Rankine-Hugoniot conditions Eqs. (4.8) become

$G(1)=2,\quad V(1)=\tfrac{1}{2},\quad Z(1)=\tfrac{3}{4}.$ (4.14)

We insert Eq. (4.9) into Eq. (4.11) and find

$\left[\frac{2}{V}+\frac{1}{V-1}-\frac{3}{3V-1}\right]\frac{dV}{d\ell}={2}\,\frac{d\ln
G}{d\ell}+\frac{3-2V}{V-1}.$ (4.15)

Using Eqs. (4.10) and (4.15), we express the derivatives of $V$ and $\ln G$
through $V$:

$\displaystyle 2\,\frac{dV}{d\ell}=-V\,\frac{3-13V+12V^{2}}{1-4V+6V^{2}},$
(4.16a) $\displaystyle 2\,\frac{d\ln
G}{d\ell}=V\,\frac{5V-1}{1-5V+10V^{2}-6V^{3}}.$ (4.16b)

Dividing Eq. (4.16b) by Eq. (4.16a) yields

$\frac{d\ln G}{dV}=\frac{5V-1}{(3V-1)(3-4V)(1-V)}$ (4.17)

which is integrated to give

$G=2^{\frac{16}{5}}\,(1-V)^{2}\,(3V-1)^{\frac{1}{5}}\,(3-4V)^{-\frac{11}{5}}.$
(4.18a) The amplitude is fixed by the boundary conditions Eq. (4.14).
Similarly integrating Eq. (4.16a) one gets
$\xi^{5}=2^{-\frac{4}{3}}\,(3V-1)^{2}\,V^{-\frac{10}{3}}\,(3-4V)^{-\frac{11}{3}}$
(4.18b) that implicitly determines $V=V(\xi)$. Finally $Z$ is given by Eq.
(4.9) which when $\gamma=3$ becomes $Z=\frac{3(1-V)V^{2}}{3V-1}.$ (4.18c)

Equations (4.18) constitute the exact solution. Now all that is remaining is
to find the value of the dimensionless parameter $A$ in Eq. (4.1). When $d=1$,
we have $\Omega_{1}=2,~{}\delta=2/3$ and $\gamma=3$. Therefore Eq. (4.13)
reduces to

$A=\left(\frac{2}{3}\right)^{2}\int_{0}^{1}d\xi\,\xi^{2}\,\frac{2V^{3}}{3V-1}\,G.$
(4.19)

Using Eqs. (4.18a)–(4.18b) one can reduce the integral in Eq. (4.19) to a
rather complicated integral over $V$ which is computed to give

$A=\frac{152}{1071}.$ (4.20)

Our basic analytical prediction about the position of the shock wave becomes

$R(t)=\left(\frac{1071}{152}\,\frac{Et^{2}}{\rho_{\infty}}\right)^{\frac{1}{3}}.$
(4.21)

Equations (4.18) and (4.21) provide the complete TvNS solution in one
dimension. In the next chapter we compare this solution with results from
direct simulations of our microscopic model.

###  4.3 Earlier Attempts to Observe TvNS Scaling in Microscopic Models

The scaling law for $R(t)$ (Eq. (4.1)) is obeyed in nuclear explosions [3, 4].
This growth law and its generalization arise in numerous applications, e.g. in
astrophysical blast waves [53, 54, 55, 56, 57, 58], laser-driven blast waves
in gas jets [59], plasma [60], granular materials [61, 62] and relativistic
blast problems [63, 64, 65, 66, 67]. It is interesting to compare the scaling
law with the microscopic simulations of the underlying many-particle system.
Such studies have been attempted only recently, which can be partly attributed
to the recent increase in computational performance. An earlier analysis [68]
of the blast in a homogeneous elastic hard-sphere gas in one and two
dimensions veriﬁed the growth law Eq. (4.1). The authors also studied the
dependence of the number of moving particles at time $t$ and the number of
collisions up to time $t$. These growth laws satisfied the predictions of the
hydrodynamic theory for the shock wave emanating from an explosion. More
detailed studies of the scaling functions, defined in Eqs. (4.1), (4.6)
through molecular dynamic simulations have recently been performed for the
hard-sphere gas in various dimensions [69, 70, 71] both with elastic and
inelastic collisional dynamics [23, 24, 25, 26]. For elastic systems in two
dimensions, a reasonable agreement between the TvNS solution and microscopic
simulations has been reported [23] when the ambient density is not too large.
More extensive simulations of the hard-sphere gas in three dimensions [24] and
two dimensions [25] however found significant departure of the simulation
results from the TvNS solution both near the center of the explosion and at
the shock front. The agreement between molecular dynamics simulations of a
Hamiltonian system and hydrodynamics still remains an open question. While the
scaling form seems robust, the scaling functions $G,V$ and $Z$ appear to be
different from the TvNS prediction. Recently, [24] tried to find a good match
between the microscopic simulations and the hydrodynamic predictions by
replacing the ideal gas equation of state with the equation of state that
included up to tenth order virial coefﬁcients. The hydrodynamics would then
predict a modified scaling solution with the same scaling exponents. It was
then observed that while there is an improvement, they still failed to get a
good match. This deviation was attributed to the violation of the assumption
of local equilibrium and the neglect of heat conduction, which is not
accounted for by the Euler equations. If the later is indeed the reason, then
we have to use a different model to describe the hydrodynamic behaviour. One
way is to use some additional dissipation terms in the Euler equations in
order to capture this effect, which is the theme of our work, to be described
in the next chapter.

##  Appendix 4.A Rankine-Hugoniot Conditions

In this section we derive the Rankine-Hugoniot conditions. The
Rankine–Hugoniot conditions, also referred to as Rankine–Hugoniot jump
conditions or Rankine–Hugoniot relations, describe the relationship between
the states on both sides of a shock wave. They are needed as boundary
conditions of the scaling solutions so that we can construct the scaling
solutions of the Euler equations. We start with the Euler equations in $d$
dimensions. Because of the radial symmetry of the problem, these three
equations for the hydrodynamic ﬁelds depend only on the radial coordinate $r$.
Thus we have the three one dimensional equations:

$\displaystyle\partial_{t}\rho+\partial_{r}(\rho v)+\frac{d-1}{r}\,\rho v=0,$
(4.22a)
$\displaystyle(\partial_{t}+v\partial_{r})v+\frac{1}{\rho}\,\partial_{r}p=0,$
(4.22b)
$\displaystyle(\partial_{t}+v\partial_{r})\ln\frac{p}{\rho^{\gamma}}=0,$
(4.22c)

where $r$ is the radial coordinate. Now consider Eq. (4.22a). Let the solution
exhibit a shock at $r=R(t)$. We now integrate this equation over a small
region, $r\in(r_{1},r_{2})$. Using the Leibnitz rule, the first term becomes:

$\int_{r_{1}}^{r_{2}}\partial_{t}\rho\
dr=\partial_{t}\int_{r_{1}}^{r_{2}}\rho\
dr+\rho(r_{1},t)\dot{r}_{1}-\rho(r_{2},t)\dot{r}_{2}~{}.$

Then Eq. (4.22a) can be written as:

$\partial_{t}\int_{r_{1}}^{r_{2}}\rho\
dr+\rho(r_{1},t)\dot{r}_{1}-\rho(r_{2},t)\dot{r}_{2}+\rho(r_{2},t)v(r_{2},t)-\rho(r_{1},t)v(r_{1},t)+\int_{r_{1}}^{r_{2}}dr\
\frac{d-1}{r}\rho v=0.$ (4.23)

We now go to the limit where $r_{1}$ is just before the shock front and
$r_{2}$ is just after the shock front. We take $r_{1}=R(t)-\epsilon$ and
$r_{2}=R(t)+\epsilon$, with $\epsilon\rightarrow 0$. Then
$\dot{r}_{1}=\dot{r}_{2}=\dot{R}(t)$, the velocity of the shock front. Also,
the first term and the last term in Eq. (4.23) vanish because $r_{1}$ and
$r_{2}$ are getting close to each other in this limit. Considering that the
shock has not reached $r=r_{2}$ at time $t$ yet, we have $v(r_{2},t)=0$. This
leads to the following equation:

$\frac{\rho(R)v(R)}{\rho(R)-\rho_{\infty}}=\dot{R}.$ (4.24)

Similarly, we get the following equations by integrating Eq. (4.22b) and Eq.
(4.22c):

$\displaystyle\frac{\rho(R)v^{2}(R)+P(R)}{\rho(R)v(R)}=\dot{R},$ (4.25)
$\displaystyle\frac{\rho(R)v(R)e(R)+P(R)v(R)}{\rho(R)e(R)}=\dot{R}.$ (4.26)

Now we use the forms of the scaling solutions as given in Eq. (4.6):

$\rho=\rho_{\infty}G(\xi),\quad v=\delta\,\frac{r}{t}\,V(\xi),\quad
T=\frac{\mu\delta^{2}}{\gamma}\,\frac{r^{2}}{t^{2}}\,Z(\xi).$

The position of the shock front corresponds to $\xi=1$. Also, the velocity of
the shock front is related to its position by $\dot{R}=\delta R/t$, as given
by Eqs. (4.1), (4.2). Substituting these forms at $\xi=1$ in Eqs. (4.24),
(4.25), (4.26) we get the Rankine-Hugoniot conditions. They are given in terms
of $G(1),V(1)$ and $Z(1)$ by:

$\displaystyle\frac{GV}{G-1}$ $\displaystyle=1,$ (4.27a)
$\displaystyle\frac{V^{2}+Z/\gamma}{V}$ $\displaystyle=1,$ (4.27b)
$\displaystyle\frac{V^{3}/2+ZV/(\gamma-1)}{V^{2}/2+Z/(\gamma(\gamma-1))}$
$\displaystyle=1.$ (4.27c)

Solving these three equations would finally lead us to the Rankine-Hugoniot
conditions:

$\displaystyle G(1)=\frac{\gamma+1}{\gamma-1},$ (4.28a) $\displaystyle
V(1)=\frac{2}{\gamma+1},$ (4.28b) $\displaystyle
Z(1)=\frac{2\gamma(\gamma-1)}{(\gamma+1)^{2}}.$ (4.28c)

These boundary conditions can then be used to fix the integration constants
and we thus derive the TvNS solution.

##  Appendix 4.B Other Scalings of the Euler Equations

In this chapter we have derived the complete scaling solution of the Euler
equations in one dimension. In this section, we emphasize the role played by
the blast wave initial conditions and re-derive the form of the TvNS scaling
in one dimension using a different method. We also explore if there are other
scaling solutions that are possible in the Euler equations. We start with the
Euler equations in one dimension:

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

where $p(x,t)$ is the pressure and we need the equation of state that relates
pressure to the conserved fields. We consider one dimensional systems that can
be described by ideal gas thermodynamics:

$\begin{split}T(x,t)&=2\mu\,\epsilon(x,t),\\\
p(x,t)&=2\rho(x,t)\epsilon(x,t)=\mu^{-1}\rho(x,t)\,T(x,t),\end{split}$ (4.30)

where $\epsilon=e-v^{2}/2$ is the internal energy per unit mass and
$\mu=(m_{1}+m_{2})/2$. Using these the Euler equations now become:

$\displaystyle\partial_{t}\rho+\partial_{x}(\rho v)=0,$ (4.31a)
$\displaystyle\partial_{t}(\rho v)+\partial_{x}(2\rho e)=0,$ (4.31b)
$\displaystyle\partial_{t}(\rho e)+\partial_{x}(3\rho ev-\rho v^{3})=0.$
(4.31c)

Let us assume the solution to the above equations be of the form:

$\displaystyle\rho(x,t)$
$\displaystyle=\rho_{\infty}+t^{b}\tilde{\rho}(xt^{a}),~{}$ (4.32a)
$\displaystyle v(x,t)$ $\displaystyle=t^{c}\tilde{v}(xt^{a}),~{}$ (4.32b)
$\displaystyle e(x,t)$ $\displaystyle=t^{d}\tilde{e}(xt^{a}).~{}$ (4.32c)

Here $\rho_{\infty}$ is the background density into which the blast is
propagating. $\tilde{\rho},\tilde{v}$ and $\tilde{e}$ are the scaling
functions of the conserved fields (which are related to $G,\ V$ and $Z$
dervied in this chapter). We aim to find the values of $a,b,c$ and $d$ here.
Now substituting these scaling forms into Eqs. (4.31a)-(4.31c), and defining
the scaling variable $z=xt^{a}$, we arrive at the following equations:

$\displaystyle
b\tilde{\rho}+az\tilde{\rho}^{\prime}+t^{a+c+1}(\tilde{\rho}\tilde{v})^{\prime}+\rho_{\infty}t^{a-b+c+1}\tilde{v}^{\prime}=0,~{}$
(4.33)
$\displaystyle(b+c)\tilde{\rho}\tilde{v}+az(\tilde{\rho}\tilde{v})^{\prime}+\rho_{\infty}ct^{-b}\tilde{v}+\rho_{\infty}azt^{-b}\tilde{v}^{\prime}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}+2t^{a-c+d+1}(\tilde{\rho}\tilde{e})^{\prime}+2\rho_{\infty}t^{a+d-c-b+1}\tilde{e}^{\prime}=0,~{}$
(4.34)
$\displaystyle(b+d)\tilde{\rho}\tilde{e}+az(\tilde{\rho}\tilde{e})^{\prime}+\rho_{\infty}dt^{-b}\tilde{e}+\rho_{\infty}azt^{-b}\tilde{e}^{\prime}+3t^{a+c+1}(\tilde{\rho}\tilde{e}\tilde{v})^{\prime}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}+3\rho_{\infty}t^{a+c-b+1}(\tilde{e}\tilde{v})^{\prime}-t^{a+3c-d+1}(\tilde{\rho}\tilde{v}^{3})^{\prime}-\rho_{\infty}t^{a+3c-b-d+1}(\tilde{v}^{3})^{\prime}=0.~{}$
(4.35)

Here prime denotes derivative wrt $z$. In order to get the scaling solution
all we need is that there should not be $x$ and $t$ anywhere in the equations
except in the form $z=xt^{a}$. Requiring now that there exists a scaling
solution, we set all powers of the time variable to zero, which then gives us
the conditions $b=0$, $a+c+1=0$ and $d=2c$. We still need one more condition
to determine the exponents completely and, as first noticed by Taylor [3, 4],
the energy conservation condition turns out to be sufficient. The constancy of
energy means that the integral

$\int_{-\infty}^{\infty}dx\rho(x,t)e(x,t)=E_{0},$ (4.36)

should be a constant, independent of time. Plugging in our scaling forms for
$\rho$ and $e$, and making a change of variables, the above equation gives

$t^{d-a}\int_{-\infty}^{\infty}dz\tilde{\rho}(z)\tilde{e}(z)=E_{0}.$ (4.37)

Demanding time-independence then requires $d=a$ and together with the earlier
equations we finally get:

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

This is the TvNS scaling in one dimension. It is easy to see that with this
choice of the exponents, Eqs. (4.32) reduce to Eqs. (4.6). For instance,
$\widetilde{\rho}=\rho_{\infty}G-\rho_{\infty}$ and
$\widetilde{v}=t^{1/3}v(x,t)=(2\alpha/3)\xi V$, where
$\alpha=[E/(A\rho_{\infty})]^{1/3}$.

Note that we got the TvNS scaling for a blast wave initial condition - energy
is initially excited in a small region and there is no energy in the ambient
medium. If however, there is energy in the surrounding medium initially then
we get a different scaling. For example, the following scaling solution leads
to a ballistic scaling and we get $a=-1,b=0,c=0$ and $d=0$:

$\displaystyle\rho(x,t)$
$\displaystyle=\rho_{\infty}+t^{b}\tilde{\rho}(xt^{a}),~{}$ (4.39a)
$\displaystyle v(x,t)$ $\displaystyle=t^{c}\tilde{v}(xt^{a}),~{}$ (4.39b)
$\displaystyle e(x,t)$ $\displaystyle=e_{0}+t^{d}\tilde{e}(xt^{a}).~{}$
(4.39c)

The above scaling will be observed if we are interested in studying the
propagation of a step discontinuity in the system (Riemann problem), which
results in a ballistic scaling.

##  Chapter 5 Blast Wave in the One Dimensional Cold Gas

In the previous chapter we derived the TvNS solution in one dimension, which
is the scaling solution of the Euler equations for a blast wave initial
condition. In this chapter we try to show an agreement between the TvNS
scaling and microscopic simulations. We consider the well-studied alternate
mass hard particle (AHP) gas system [68, 72, 73, 74, 75, 76, 77, 78]. We
describe the relation between hydrodynamic fields and the microscopic
variables of this system in Sec. 5.1. We then go on to brief some of the
earlier works which studied various aspects of the AHP gas in Sec. 5.2. We
study the blast problem in the AHP gas system. We excite the system to a
“blast wave” initial condition, which will be described in Sec. 5.3. This will
result in the propagation of a shock wave in the system. We first compare the
results of microscopic simulations with the TvNS solution in Sec. 5.4. We find
a very good agreement between them everywhere except in a small region near
the core of the blast. We attribute this deviation near the core to the non-
negligible contribution of heat conduction, which is not captured by the Euler
equations. So, we model our system by adding viscosity and the heat
conductivity terms to the Euler equations, which then become the Navier-
Stokes-Fourier (NSF) equations [79], which are described in Sec. 5.6. Now we
observe a very good agreement between the scaling solutions of the microscopic
simulations and the NSF equations. We thus propose that the NSF equations
capture the hydrodynamics of the one dimensional AHP gas more accurately than
the Euler equations. Then in Sec. 5.7 we analyse the role of the conductivity
term in more detail and study its role in the evolution of the conserved
quantities in the system near the core of the blast.

###  5.1 The Alternate Mass Hard Particle Gas

As mentioned in chapter 1, the reason for considering alternate masses is that
we need some sort of integrability breaking. Having equal masses would
essentially mean a system of non-interacting particles, which is integrable
and leads to a macroscopic number of conservation quantities. Such systems are
described by generalized hydrodynamics. We are however, interested in
describing the hydrodynamics of a system that only has a few conserved
quantities. We choose odd-numbered particles to have mass $m_{1}$, so even-
numbered particles have mass $m_{2}$ ($\neq m_{1}$). In this case, the usual
hydrodynamics is sufficient. The only known conserved quantities of this
system are the total number of particles, the total momentum and the total
energy, and we expect a hydrodynamic description in terms of the corresponding
conserved fields. The hydrodynamic fields $(\rho,v,E)$, are related to the
microscopic variables:

$\displaystyle\begin{split}&\rho(x,t)=\sum_{j}\left\langle
m_{j}\delta(q_{j}(t)-x)\right\rangle,\\\ &\rho(x,t)v(x,t)=\sum_{j}\left\langle
m_{j}v_{j}\delta(q_{j}(t)-x)\right\rangle,\\\
&E(x,t)=\rho(x,t)e(x,t)=\tfrac{1}{2}\sum_{j}\left\langle
m_{j}v_{j}^{2}\delta(q_{j}(t)-x)\right\rangle,\end{split}$ (5.1)

where $\langle...\rangle$ indicates an average over an initial distribution of
micro-states that correspond to the same initial macro-state.

We study the blast problem in this system. We consider a gas with density
$\rho_{\infty}$ and particles initially at rest (zero pressure and
temperature). The pressure in front of the shock wave remains zero throughout
the evolution if it is zero initially. We suddenly inject, at time $t=0$,
energy $E$ into a localized region. The details of injection are not
important, the explosion quickly becomes radially symmetric and keeps [80, 81]
its spherical shape, so the hydrodynamic functions depend on the distance $r$
from the center of the explosion and time $t$. Since the temperature is zero
in the unperturbed region, local equilibrium and hydrodynamic behaviors may
emerge only sufficiently far behind the shock, namely in the region where
particles have already undergone a few collisions. Hence the blast problem
with particles initially at rest is an excellent laboratory to study
deviations from hydrodynamics. Also, at any time only a finite number of
particles move thereby allowing us to simulate an infinite-particle system.

###  5.2 Earlier Studies done on the AHP gas

Here we discuss some works which studied various properties of the alternate
mass hard particle (AHP) gas. The AHP gas has also been extensively studied in
the context of the breakdown of Fourier’s law of heat conduction in one
dimension [72, 73, 74, 75, 76, 77, 78, 82, 83]. Remarkably, the thermal
conductivity $\kappa$ diverges with system size $L$ as $\kappa\sim L^{1/3}$.
Some theoretical understanding of this has been achieved using nonlinear
fluctuating hydrodynamics [84, 85, 86, 87], which makes very detailed
predictions about the form of equilibrium dynamical correlation functions.
These predictions have been verified for two microscopic models, the AHP gas
[86, 88] and the FPUT chain [89]. As far as the evolution of macroscopic
profiles is considered, the AHP gas was used recently [90] to study the so-
called Riemann problem [91], where one considers an infinite system with
domain wall initial conditions. A comparison was made between the simulations
and the predictions from the Euler equations. Note that the blast-wave initial
conditions lead to a different class of self-similar solutions than those that
one gets for the Riemann problem [ballistic scaling with $R(t)\sim t$]. We
have already seen this in Appendix (4.B), where we derived the ballistic
scaling for the hydrodynamic equations for a different set of initial
conditions. The one-dimensional hard-particle gas has also been used to probe
the breakdown of the hydrodynamic description [92, 93]. It was observed [93]
that the field profiles in front of a rapidly moving piston differ from
hydrodynamic predictions. The splash problem in the 1D AHP gas has also been
studied recently [94].

###  5.3 “Blast Wave” Initial Conditions for Microscopic Simulations

We consider an initial macro-state where the gas has a finite uniform density
$\rho_{\infty}$, zero flow velocity $v$, and is at zero temperature everywhere
except in a region of width $\sigma$ centered at $x=0$. This is the region of
the blast and has a specified finite temperature profile. In numerical
simulations, two different profiles for the initial energy have been used:

$\displaystyle{\rm(i)~{}Gaussian~{}profile}~{}~{}E(x,0)=\frac{E}{\sqrt{2\pi\sigma^{2}}}e^{-x^{2}/{2\sigma^{2}}},$
(5.2)
$\displaystyle{\rm(ii)~{}Box~{}profile}~{}~{}E(x,0)=\begin{cases}\frac{E}{2\sigma},~{}-\sigma<x<\sigma,\\\
0,~{}~{}|x|>\sigma.\end{cases}~{}$ (5.3)

In both cases, the total energy of the blast is $E$. We now describe how one
can realize these macro-states in the microscopic simulations of the AHP gas.
Different choices of the initial ensemble of micro-states can lead to the same
average macroscopic profile. Here we discuss two possible ensemble choices.

(A) Ensemble where energy and momentum are fixed on average: In this case the
two initial macro-states can be realized as follows:

(i) Gaussian temperature profile: We distribute $N$ particles labelled
$i=-N/2+1,-N/2+2,\ldots,N/2$ uniformly in the interval $(-L/2,L/2)$. The
number density is $n_{0}=N/L=\rho_{\infty}/\mu$, where $\rho_{\infty}$ is the
ambient mass density of the gas. For the $N_{c}$ particles at the center, with
$i=-N_{c}/2+1,\ldots,N_{c}/2$, the velocities are chosen from the Maxwell
distribution with temperature $T=2E/N_{c}$:
$\text{Prob}(v_{i})=\sqrt{m_{i}/(2\pi T)}e^{-m_{i}v_{i}^{2}/(2T)}$. The
velocities of other $N-N_{c}$ particles are set to zero. We note that
$E(x,0)=\langle\sum m_{i}v_{i}(0)^{2}\delta(x-q_{i}(0))/2\rangle$. For large
$N$, each particle’s position is a Gaussian with mean $\bar{q}_{i}=i/n_{0}$
and variance $\sigma^{2}=L/(4n_{0})$. Defining the length scale
$s=N_{c}/n_{0}$, one can express

$\displaystyle E(x,0)$
$\displaystyle=\sum_{i=-N_{c}/2+1}^{N_{c}/2}\frac{T}{2\sqrt{2\pi\sigma^{2}}}e^{-(x-\bar{q}_{i})^{2}/(2\sigma^{2})}$
$\displaystyle\approx\frac{n_{0}T}{2}\int_{-s/2}^{s/2}dy\frac{e^{-(x-y)^{2}/(2\sigma^{2})}}{\sqrt{2\pi\sigma^{2}}}$
$\displaystyle=\frac{n_{0}T}{4}\left[\operatorname{erf}{\left(\frac{s-2x}{2\sqrt{2}\sigma}\right)}+\operatorname{erf}{\left(\frac{s+2x}{2\sqrt{2}\sigma}\right)}\right]$
(5.4)

via the error function
$\operatorname{erf}(z)=\frac{2}{\sqrt{\pi}}\int_{0}^{z}du\,e^{-u^{2}}$. When
$s\ll\sigma$, Eq. (5.4) simplifies to

$\displaystyle E(x,0)\approx
E\,\frac{e^{-x^{2}/(2\sigma^{2})}}{\sqrt{2\pi\sigma^{2}}}\,.$ (5.5)

(ii) Box temperature profile: In this case we distribute $(N-N_{c})/2$
particles uniformly in $(-L/2,-s/2)$, and $(N-N_{c})/2$ particles in
$(s/2,L/2)$, with $s=N_{c}/n_{0}$. The remaining $N_{c}$ particles are
distributed uniformly in $(-s/2,s/2)$. The mean number density is again
$n_{0}=N/L$. The velocities of the $N_{c}$ particles in the center are chosen
from the Maxwell distribution with temperature $T=2E/N_{c}$, viz.
$\text{Prob}(v_{i})=\sqrt{m_{i}/(2\pi T)}e^{-m_{i}v_{i}^{2}/(2T)}$; the
initial velocities of all other particles are set to zero.

(B) Ensemble where initial energy is fixed exactly to $E$ and momentum exactly
to $0$: In this case we distribute particles in space using the same protocol
as before. We now choose the velocities of only $N_{c}/2$ particles
($-N_{c}/2+1,-N_{c}/2+2,\ldots,0$) from the Maxwell distribution,
$\text{Prob}(v_{i})=\sqrt{m_{i}/(2\pi T)}e^{-m_{i}v_{i}^{2}/(2T)}$, where
$T=2E/N_{c}$. For every realization, we rescale these velocities by a constant
factor such that the total energy of the $N_{c}/2$ particles is exactly $E/2$.
The remaining particles ($i=1,2,\ldots,N_{c}/2$) are assigned mirror-image
momenta: $p_{i}=p_{-i+1}$. Thus we have an ensemble of initial micro-states
which all have total energy exactly equal to $E$ and total momentum exactly
zero.

Simulation details: In all our simulations we took $m_{1}=1,\ m_{2}=2$ (so
that $\mu=(m_{1}+m_{2})/2=1.5$), $\rho_{\infty}=1.5$, $E=32$ and $N_{c}=32$.
In largest simulations, we took $N=24000$, $L=24000$ and averaged over an
ensemble of ${\cal{R}}=10^{5}$ initial conditions. In all cases, simulations
were carried out up to times such that the shocks have not reached the
boundaries, so the results correspond to an infinite system. The molecular
dynamics simulations for this system can be done efficiently using an event-
driven algorithm that updates the system between successive collisions.

###  5.4 Microscopic Simulations and Comparison with the TvNS Solution

We now present the results of the microscopic simulations. We mainly discuss
the results of the Gaussian initial temperature profile and with the averaging
over an initial ensemble of micro-states where the total energy and total
momentum are fixed exactly (ensemble (B) in the previous section). We also
show a comparison with the empirical profiles (obtained from spatial coarse-
graining of single realization) and outline results obtained for the box
profile and for the ensemble choice where energy and momentum are fixed on
average.

In Fig. 5.1 we illustrate at the microscopic level a typical blast in one
dimension, by showing the trajectories of the moving particles. Around $240$
particles are shown in this simulation; other parameters are $N_{c}=32$,
$\rho(x,0)=\rho_{\infty}=1.5$ and the initial energy profile was taken to be a
Gaussian (see Sec. 5.3). The expected shock positions $\pm R$, with $R(t)$
given by Eq. (4.1), are shown by the black dashed lines.

In Fig. 5.2 we show the spatio-temporal evolution of $\rho(x,t)$, $v(x,t)$,
$T(x,t)$ in a gas initially at rest perturbed by a localized Gaussian
temperature profile. We see a clear sub-ballistic evolution of the fields. The
mass density and velocity are peaked around the shock front, while the
temperature has an additional peak at the center.

Figure 5.1: The trajectories of about $240$ moving particles in a typical
blast, starting from initial conditions corresponding to a Gaussian initial
temperature profile. The expected shock positions, $\pm R(t)$, are shown as
dashed lines. The stationary particles are not shown. In this realization, the
left-most and right-most moving particles are quite a bit ahead of where they
are expected to be according to the deterministic hydrodynamics.

In Fig. 5.3 we plot the early time evolution of the hydrodynamic fields
$\rho(x,t)$, $v(x,t)$, and $T(x,t)$ as functions of position $x$ at different
times $t$. In Fig. 5.4 we plot the evolution of the blast wave at long times.
The lower panel in Fig. 5.4 shows
$\widetilde{\rho}=\rho_{\infty}G-\rho_{\infty},~{}\widetilde{v}=t^{1/3}v(x,t)=(2\alpha/3)\xi
V$ and $\widetilde{T}=t^{2/3}T(x,t)=(4\mu\alpha^{2}/27)\xi^{2}Z$, where
$\alpha=[E/(A\rho_{\infty})]^{1/3}$, as functions of the scaling variable
$x/t^{2/3}$. We find a very good scaling collapse of the data confirming the
expected TvNS scaling. For comparison, we plot the exact predictions Eq.
(4.18). We find an almost perfect agreement between the simulation results and
the TvNS solution except for deviations near the center of the blast. As we
discuss later in Sec. 5.7, the TvNS solution predicts a divergence of the
temperature field with $\widetilde{T}(\xi)\sim\xi^{-1/2}$ and this implies
that the effect of the heat conduction term in the hydrodynamic equation for
energy conservation cannot be neglected. Thus one has to deal with the full
Navier-Stokes-Fourier equations instead of the Euler equations.

Figure 5.2: Heat maps of the spatio-temporal evolution of (a) density, (b)
velocity, and (c) temperature starting from initial conditions corresponding
to a Gaussian initial temperature profile and an ensemble average over
$10^{4}$ initial conditions.

(a)(b)(c)

Figure 5.3: Early time evolution of (a) density, (b) velocity and (c)
temperature obtained by ensemble averaging starting from initial conditions
corresponding to a Gaussian initial temperature profile and an average over
$10^{5}$ initial conditions.

(a)(b)(c)

(d)(e)(f)

Figure 5.4:  (a,b,c) Late time evolution of 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 gives a good collapse of the data at the longest times
everywhere except near the origin. The exact scaling solution of the
hydrodynamic Euler equations is shown by the black dashed line.

(a)(b)(c)

Figure 5.5: Evolution of the empirical density, velocity, and temperature
obtained from a single realization, and after spatial coarse-graining
according to Eqs. (5.6-5.8) and using sliding bins. The black dashed lines
indicate the ensemble-averaged fields, and we see that apart from the
fluctuations seen in the empirical fields, there is a good agreement between
the two. The coarse-graining length scale, $\Delta=50$, and same initial
conditions as in Fig. 5.4 are used.

###  5.5 Various Aspects of Microscopic Simulations

####  5.5.1 Comparison with the Empirical Fields

Ideally, we would like to see a comparison of the hydrodynamic theory with
simulations for the empirical fields that correspond to the physical situation
where we observe a single realization, and there is only a spatial averaging
of the microscopic degrees. We divide space into boxes of size $\Delta$ that
are small compared to the length of the region with moving particles but still
typically contain many particles. The empirical fields are defined as:

$\displaystyle\bar{\rho}(x,t)=\frac{1}{\Delta}\sum m_{j}\delta[q_{j}(t),x],$
(5.6) $\displaystyle\bar{\rho}(x,t)\bar{v}(x,t)=\frac{1}{\Delta}\sum
m_{j}v_{j}\delta[q_{j}(t),x],$ (5.7)
$\displaystyle\bar{E}(x,t)=\bar{\rho}(x,t)\bar{e}(x,t)=\frac{1}{2\Delta}\sum
m_{j}v_{j}^{2}\delta[q_{j}(t),x],$ (5.8)

where $\delta[q_{j},x]$ is an indicator function taking value $1$ if
$x-\Delta/2\leq q_{j}\leq x+\Delta/2$ and zero otherwise. The system is
initially excited in the same manner as the previous subsection with the same