11institutetext: 1 International Centre for Theoretical Sciences, Bengaluru,
India-560089
2 Department of Physics, Boston University, Boston, MA 02215, USA
3 Santa Fe Institute, 1399 Hyde Park Road, Santa Fe, NM 87501, USA

# Blast waves in the zero temperature hard sphere gas: double scaling
structure

Sahil Kumar Singh1 Subhadip Chakraborti1 Abhishek Dhar1 and P. L. Krapivsky2,3

(Received: date / Accepted: date)

###### Abstract

We study the blast generated by sudden localized release of energy in a cold
gas. Specifically, we consider one-dimensional hard-rod gas and two-
dimensional hard disc gas. For this problem, the Taylor-von Neumann-Sedov
(TvNS) solution of Euler equations has a self-similar form. The shock wave
remains infinitely strong for the zero-temperature gas, so the solution
applies indefinitely. The TvNS solution ignores dissipation, however. We show
that this is erroneous in the core region which, in two dimensions, expands as
$t^{2/5}$ while the shock wave propagates as $t^{1/2}$. A new self-similar
solution depending on the scaling variable $r/t^{2/5}$ describes the core,
while the TvNS solution describes the bulk. We demonstrate this from a
numerical solution of the Navier-Stokes (NS) equations and from molecular
dynamics simulations for a gas of hard rods in two dimensions and hard rods in
one dimension. In both cases, the shock front position predicted by NS
equations and by the TvNS solution agrees with that predicted by molecular
dynamics simulations. However, the NS equations fail to describe the near-core
form of the scaling functions.

###### Contents

  1. 1 Introduction
  2. 2 Blast in the 2D hard-disc gas
    1. 2.1 TvNS solution
    2. 2.2 Effect of dissipation on scaling in the blast core
    3. 2.3 Numerical results
  3. 3 Blast in 2D ideal gas
    1. 3.1 TvNS solution
    2. 3.2 Core scaling solution
    3. 3.3 Numerical results
  4. 4 Blast in the hard rod gas
    1. 4.1 TvNS solution
    2. 4.2 Core scaling solution
    3. 4.3 Numerical results
  5. 5 Hard point gas in $d-$dimensions: Inner and outer solutions
    1. 5.1 Outer solution
      1. 5.1.1 Euler equations in $d-$dimensions
      2. 5.1.2 Self-similar solution
      3. 5.1.3 Singular behavior
    2. 5.2 Inner solution
      1. 5.2.1 Navier-Stokes equations
      2. 5.2.2 Heuristic analysis for core scaling
  6. 6 Discussion

##  1 Introduction

The non-equilibrium state of a fluid is usually described by hydrodynamics.
However, on a microscopic level, the individual particles follow Newton’s laws
of motion. How does hydrodynamics emerge from microscopic Newtonian dynamics?
A rigorous derivation of hydrodynamics starting from microscopic Newtonian
dynamics is lacking. Phenomenological textbook derivations involve applying
conservation of mass, momentum and energy to a parcel of fluid, and there are
also derivations based on the Boltzmann equation landau9 . The latter assume
that the hydrodynamic limit is achieved. More precisely, the spatial
variations of the coarse-grained variables are sufficiently slow, i.e., local
equilibrium is achieved, and the deviations from local equilibrium are small.
Thus to zeroth order, one assumes local equilibrium and gets the Euler
equations which do not have dissipation. To next order, one includes
deviations from local equilibrium to get Navier-Stokes equations which include
dissipation as higher order derivative corrections to the Euler equations.

In some cases the hydrodynamic limit may not be achieved. For example, when
there is shock formation, there is a sharp variation in the hydrodynamic
fields, and one may expect hydrodynamics to break down. One situation where
there can be shock formation is when the energy released from a blast is left
to propagate freely in an otherwise cold gas. Also, the center of the blast is
a low density region where collisions are rare and local equilibrium is
difficult to reach. Thus the blast problem is a classic problem to test the
validity of hydrodynamics. It was first studied in the context of atomic
explosions and underlies the behavior of many astrophysical systems, see
Ostriker ; Goodman ; Gal-Yam ; Chevalier ; Burrows . An exploding star sends a
blast wave into the stellar medium and later into the interstellar medium. The
released energy is tremendous in astrophysical applications, and hence the
assumption of a zero temperature, zero pressure ambient gas, in which the
explosion propagates, is an excellent approximation.

One of the central problems is to determine how the shock wave advances and to
compute the hydrodynamic fields behind the shock. The blast problem is
traditionally studied in the framework of an ideal compressible gas neglecting
the effect of dissipation. Dimensional analysis sedov ; landau  allows one to
express the radius $R(t)$ of the shock wave through the time $t$ counted from
the moment of explosion, the released energy $E_{0}$ and the density
$\rho_{\infty}$ in front of the shock

$R(t)=\left(\frac{E_{0}t^{2}}{A_{d}\,\rho_{\infty}}\right)^{\frac{1}{d+2}}.$
(1)

Here $d$ is the spatial dimension. The amplitude $A_{d}$ is a priori unknown,
fixing it is a part of the solution. It was shown by Taylor, von Neumann, and
Sedov taylor1 ; taylor1950 ; vN ; sedov1 ; sedov  that at long times the
hydrodynamic fields behind the shock front have a self-similar scaling form
that can be obtained from a solution of the Euler equations. This solution
will be referred to as the Taylor-von Neumann-Sedov (TvNS) solution. One of
the main issues that we address here is the effect of dissipation on the TvNS
predictions.

By writing the Navier-Stokes (NS) equations, one can easily see that
dissipation in the flow behind the shock is indeed negligible, apart from the
region near the center of the explosion where the temperature obtained from
ideal hydrodynamics, i.e., Euler equations, diverges. This is physically
impossible since the energy injected is finite. Dissipation is thus expected
to affect the flow in the region surrounding the origin which we will refer to
as the core region. Earlier studies of the influence of heat conduction Bethe
; zel ; Oppenheim ; Steiner  and viscosity Johnny ; Brode ; Latter ; Plooster
on the blast were performed in the realm of atomic explosions.

Recently, a comparison was made between the results of Newtonian dynamics from
molecular dynamics (MD) simulations and that of hydrodynamics in the one-
dimensional (1D) alternating mass hard particle (AHP) gas for a blast-like
chakraborti ; ganapa  and splash chakraborti2022  initial condition, and it
was found that the TvNS solution describes, surprisingly accurately, the
results of MD simulations except in the core of the blast. As expected, due to
the effect of dissipation, the TvNS solution fails in the core. It was noted
in chakraborti ; ganapa  that the core has size $X(t)\sim t^{38/93}$ which, at
large times, is much smaller than the size of the blast given from Eq. (1),
for $d=1,$ by $R(t)\sim t^{2/3}$. The core region does not have a sharp
boundary, unlike the shock wave front. An analysis of the hydrodynamic
equations including dissipative terms yields a different scaling solution in
the core which was compared with results from MD. While dissipation kills the
unphysical divergences of the Euler solution, the detailed agreement between
simulations and dissipative hydrodynamics is not very good in the core region.
One reason could be the presence of anomalous heat transport in $d=1$ LLP2003
; dhar2008 . Fourier’s law in one dimension breaks down, and it is replaced by
a non-local law kundu1 ; kundu2 . Thus, we expect hydrodynamics to break down
in one dimension, and such a breakdown has indeed been observed in hurtado .
Thus, the near-perfect agreement found between the TvNS solution and molecular
dynamics simulations in chakraborti ; ganapa  is quite surprising. It can only
be speculated that since the dissipative region is small in the 1D blast
problem, much of the perturbed region in the gas is described by dissipation-
less Euler equations.

A natural question is to investigate the effect of dissipation in higher
dimensions. The first study Trizac2  comparing the predictions of molecular
dynamics with those of Euler equations in two dimensions (2D) found a
reasonable agreement between the two. Other studies in two and three-
dimensional hard-sphere gases joy1 ; joy2 ; amit  found disagreement between
the TvNS prediction and MD simulations (even away from the core region). It
turns out private  that the reason for the observed discrepancy in joy1  arose
from an incorrect estimation of the value of initial energy of the blast used
in the simulation. We have verified this in our MD simulations where we find
perfect agreement with TvNS in the outer region (see Fig. (4) .

In this work we first consider the blast in a 2D gas of hard discs and in a
toy 2D gas of point particles. We also consider a 1D gas of hard rods with
alternating masses (to ensure non-integrability) to see if the finite rod size
changes the remarkable agreement between MD results and hydrodynamics seen in
chakraborti ; ganapa . Our main results are:

(i) For the 2D system, we show that, as in 1D, the dissipative terms in the NS
equations introduce a new growing length scale which we denote by $X(t)$. We
find a new self-similar solution in the core, different from the TvNS scaling
solution. This core scaling solution is verified from the numerical solution
of the NS equations. Thus there are two scaling regions: the bulk region
described by TvNS scaling and the core region described by the dissipative
scaling. We further show that the two regions are connected by the rules of
asymptotic matching.

(ii) We find that NS results agree perfectly with TvNS in the bulk region and
also agree with the MD simulation results. The MD results satisfy TvNS scaling
in the bulk region and also the core scaling in the central region. However,
the scaling functions in the core region are different from the predictions of
hydrodynamics.

(iii) For the gas of hard rods in 1D we find an excellent agreement between
TvNS, NS and MD simulations in the bulk region, including for the position of
the shock front and for the scaling functions. In the blast core, the NS
scaling is observed in the MD simulations, but the scaling functions do not
agree with those obtained from NS.

We summarize here the plan of the paper. In Sec. (2.1), we review the TvNS
solution. In Sec. (2.2), we discuss the effect of dissipation and the new
scaling solution that exists in the core. In Sec. (2.3), we present numerical
results of the NS equations which verify the core scaling forms. We also show
that MD simulations satisfy this scaling and point out the agreement between
MD and NS. In Sec. (3), we present results on the 2D gas when virial
corrections in the equation of state are ignored. In Sec. (4), we present
results on hard rods. Sec. (5) contain a brief discussion on the form of the
inner and outer solutions in arbitrary dimensions. We conclude in Sec. (6)
with a discussion.

##  2 Blast in the 2D hard-disc gas

We first consider a gas of hard discs of diameter $a$, and mass $m$, whose
only interaction is via elastic collisions that conserve momentum and energy.
Particles move with constant velocities between collisions. The finite size of
the particles implies that the equation of state is no longer given by the
ideal gas equation. This feature complicates the solution of the hydrodynamic
equations. At a microscopic level, the blast initial condition consists of
taking an infinite gas of discs that are distributed uniformly in the 2D plane
with mass density, $\rho_{\infty}$, and all at rest. Particles inside a
localized region are then excited such that their total energy is $E_{0}$, and
the total momentum is zero. The excitation then evolves in a radially
symmetric way and we are interested, in particular, in how the three conserved
fields of mass density, $\rho(r)$, momentum density, $p(r)$, and energy
density, $e(r)$, evolve with time. We first review the TvNS solution as
presented in joy1 , and then we discuss the effect of dissipation by
considering the NS equations.

###  2.1 TvNS solution

The classic TvNS solution was first obtained for a 3D ideal gas using Euler
equations and ignoring any virial corrections caused by the finite cross
sectional area of the colliding molecules. A similar analysis was done in joy1
for hard-disc gas, but by including the virial corrections in the equation of
state.

The first step in finding TvNS solution is to find how the shock front $R(t)$
grows as a function of time. As discussed above this can be found from purely
dimensional analysis. The only variables that the shock front can depend on at
large times are time $t$, energy injected $E_{0}$, and the ambient density
$\rho_{\infty}$. Applying dimensional analysis, we get
$R(t)=[E_{0}t^{2}/(A_{d}\rho_{\infty})]^{\frac{1}{d+2}}$, where $A_{d}$ is an
unknown dimensionless constant, that only depends on the volume fraction
$(\rho_{\infty}/m)V_{d}a^{d}$, where $V_{d}=2^{-d}\pi^{d/2}/\Gamma(d/2+1)$ is
the volume of a $d$-dimensional sphere of unit diameter. The Euler equations
in any dimension are given by:

$\displaystyle\partial_{t}\rho+\partial_{i}(\rho v_{i})=0,$ (2a)
$\displaystyle\partial_{t}(\rho v_{i})+\partial_{j}(\rho
v_{i}v_{j})+\partial_{i}P=0,$ (2b) $\displaystyle\partial_{t}(\rho
e)+\partial_{i}(\rho ev_{i})+\partial_{i}(Pv_{i})=0,$ (2c)

where $P$ is the pressure, which can be related to the other hydrodynamic
fields via equation of state, and $e$ is the energy per unit mass,
$e={v^{2}}/{2}+{dT}/{2}$. Henceforth, in most of the discussions here we set
$a=k_{B}=m=1$, unless otherwise specified.

For the 2D hard disc gas, the equation of state at low densities is given by
the virial equation of state:

$P=\rho T\Big{(}1+\sum_{n=2}^{\infty}B_{n}\rho^{n-1}\Big{)}.$ (3)

Following joy1 , we only take terms up to $n=10$. The values of the constants
$B_{n}$ are given in joy1 . According to TvNS, the blast wave initial
condition evolves, at long times to a self-similar scaling form. We define the
following scaling variable:

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

In terms of this dimensionless variable, one again uses dimensional analysis
to make the following ansatz for the long time solution of Eqs. (2) for the
density, velocity and temperature fields:

$\displaystyle\rho$ $\displaystyle=\rho_{\infty}G(\xi),$ (5a) $\displaystyle
v$ $\displaystyle=\frac{r}{t}V(\xi),$ (5b) $\displaystyle T$
$\displaystyle=\frac{r^{2}}{t^{2}}Z(\xi).$ (5c)

Putting these scaling forms into the Euler equations, we get the following
ODEs for the scaling functions:

$\displaystyle\xi\Big{(}V-\frac{1}{2}\Big{)}G\frac{dV}{d\xi}+\xi\frac{d}{d\xi}\left[ZGB(G)\right]-GV+GV^{2}+2GZB(G)=0,$
(6a) $\displaystyle\xi\Big{(}V-\frac{1}{2}\Big{)}\frac{dG}{d\xi}+\xi
G\frac{dV}{d\xi}+2GV=0,$ (6b) $\displaystyle-\frac{\xi
B(G)}{G}\frac{dG}{d\xi}+\frac{\xi}{Z}\frac{dZ}{d\xi}+2\frac{(V-1)}{(V-1/2)}=0,$
(6c)

where

$B(G)=1+\sum_{n=2}^{\infty}B_{n}\rho_{\infty}^{n-1}G^{n-1}.$ (7)

Under the dynamics of the Euler equations, energy is conserved, and this gives
us the further constraint:

$\int_{0}^{R(t)}dr\leavevmode\nobreak\ 2\pi r\rho e=E_{0},$ (8)

which, in terms of the scaling functions becomes:

$2\pi\int_{0}^{1}d\xi\leavevmode\nobreak\
\xi^{3}G(\xi)\Big{(}\frac{V^{2}(\xi)}{2}+Z(\xi)\Big{)}=A_{2}.$ (9)

The boundary conditions required to solve these ODEs are provided by the
Rankine-Hugoniot conditions which specify the discontinuity of the fields
across the shock:

$\displaystyle\frac{\rho(R)v(R)}{\rho(R)-\rho_{\infty}}=\dot{R},$ (10a)
$\displaystyle\frac{\rho(R)v^{2}(R)+P(R)}{\rho(R)v(R)}=\dot{R},$ (10b)
$\displaystyle\frac{\rho(R)v(R)e(R)+P(R)v(R)}{\rho(R)e(R)}=\dot{R}.$ (10c)

These, in terms of the scaling functions are given by:

$\displaystyle\frac{1}{G(1)}\left(1+\frac{2}{B[G(1)]}\right)=1,$ (11a)
$\displaystyle V(1)=\frac{1}{G(1)\leavevmode\nobreak\ B[G(1)]},$ (11b)
$\displaystyle Z(1)=\frac{1}{2}V^{2}(1).$ (11c)

These equations can be solved numerically to find the values of the scaling
functions at the shock front, which can be used as boundary values to solve
the ODEs. The dimensionless constant $A_{2}$ can be found using Eq. (9).

###  2.2 Effect of dissipation on scaling in the blast core

In the Euler framework, it is assumed that the hydrodynamic fields satisfy the
Euler equations throughout the blast and that dissipative terms can be
neglected. This assumption is valid only if the dissipative terms
corresponding to the TvNS solution are small compared to the Euler terms. The
TvNS solution predicts a divergence $Z(\xi)\sim\xi^{-4}$ for small $\xi$ which
means that the temperature field will have a steep spatial variation. Since
the spatial derivatives in the dissipative terms in the NS equations are of
one order higher than the Euler terms, this signals that the dissipative terms
corresponding to the TvNS solution are much larger than the Euler terms, thus
suggesting that dissipation becomes important near the core. However, far from
the core, the slopes of the various scaling functions are not steep, so
dissipation is not important in those regions. Since dissipation terms far
from the core are not important, we expect the ODEs to still hold in the bulk
region.

We can estimate the size of the core where dissipation becomes important. The
Navier-Stokes equations in $d$-dimensions are given by:

$\displaystyle\partial_{t}\rho+\partial_{i}(\rho v_{i})=0,$ (12a)
$\displaystyle\partial_{t}(\rho v_{i})+\partial_{j}(\rho
v_{i}v_{j})+\partial_{j}\sigma_{ij}=0,$ (12b) $\displaystyle\partial_{t}(\rho
e)+\partial_{i}(\rho
ev_{i})+\partial_{i}(\sigma_{ij}v_{j})+\partial_{i}(-\kappa\partial_{i}T)=0,$
(12c)

with the stress tensor $\sigma$ given by:

$\sigma_{ij}=P\delta_{ij}-\mu[\partial_{i}v_{j}+\partial_{j}v_{i}-\frac{2}{d}(\nabla\cdot\boldsymbol{v})\delta_{ij}]-\zeta(\nabla\cdot\boldsymbol{v})\delta_{ij},$
(13)

where $\mu$ is the shear viscosity, $\zeta$ is the bulk viscosity and $\kappa$
is the heat conductivity. Kinetic theory predicts that $\mu=D_{\mu}\sqrt{T}$,
$\kappa=D_{\kappa}\sqrt{T}$ and $\zeta=0$ for a monoatomic gas in the low
density limit landau9 ; beijeren , where $D_{\kappa}$ and $D_{\mu}$ are some
constants. However, we emphasize that these results are only first order terms
in a virial expansion landau9 . For a radially symmetric blast, all the
hydrodynamic fields are radially symmetric. Using this fact, the NS equations
become (for $d=2$):

$\displaystyle\partial_{t}\rho+\frac{\rho v}{r}+\partial_{r}(\rho v)=0,$ (14a)
$\displaystyle\rho(\partial_{t}v+v\partial_{r}v)+\partial_{r}P=\frac{1}{r^{2}}\partial_{r}[r^{2}\mu(\partial_{r}v-r^{-1}v)],$
(14b)
$\displaystyle\rho(\partial_{t}e+v\partial_{r}e)+\frac{1}{r}\partial_{r}(rPv)=\mu(\partial_{r}v-\frac{v}{r})^{2}+\frac{v}{r^{2}}\partial_{r}[r^{2}\mu(\partial_{r}v-\frac{v}{r})]+\frac{1}{r}\partial_{r}(r\kappa\partial_{r}T).$
(14c)

It is known joy1 , that the TvNS solution predicts (in 2D) that for small
$\xi$, $G\sim\xi^{2}$, $V-1/4\sim\xi^{4}$ and $Z\sim\xi^{-4}$ [see also Sec.
(3)]. Using this small $\xi$ behaviour predicted by TvNS solution, we can
estimate if the dissipation terms are indeed small, e.g, by looking at the
ratio of the terms on the right and left hand sides of Eq. (14b). This shows
that dissipation is small for
$r\gtrsim{D_{\mu}^{1/5}E_{0}^{1/5}}{\rho_{\infty}^{-2/5}}t^{2/5}$, thus
verifying the assumption leading to the TvNS solution. However, for
$r\lesssim{D_{\mu}^{1/5}E_{0}^{1/5}}{\rho_{\infty}^{-2/5}}t^{2/5}$,
dissipation terms are larger than Euler terms, and thus the assumption leading
to TvNS solution breaks down. Hence, we expect TvNS solution not to be valid
for $r\lesssim{D_{\mu}^{1/5}E_{0}^{1/5}}{\rho_{\infty}^{-2/5}}t^{2/5}$. This
introduces a new length scale, which is an estimate of the size of the core,

$X(t)=\frac{D_{\mu}^{1/5}E_{0}^{1/5}}{\rho_{\infty}^{2/5}}t^{2/5},$ (15)

and which grows with time as $t^{2/5}$. From this length scale, we can find an
estimate for the density, velocity and temperature in the core, which we
denote respectively by $\rho^{*}$, $v^{*}$ and $T^{*}$. This will give us new
density, velocity and temperature scales (here $\xi_{*}={X(t)}/{R(t)}$ is the
scaled size of the core):

$\displaystyle\rho^{*}=\rho_{\infty}\xi_{*}^{2}=b_{\rho}t^{-1/5},$ (16a)
$\displaystyle v^{*}=\frac{X(t)}{t}=b_{v}t^{-3/5},$ (16b) $\displaystyle
T^{*}=\frac{X^{2}(t)}{t^{2}}\xi_{*}^{-4}\sim b_{T}t^{-4/5},$ (16c)

where $b_{\rho}={A^{1/2}D_{\mu}^{2/5}\rho_{\infty}^{7/10}}{E_{0}^{-1/10}}$,
$b_{v}=D_{\mu}^{1/5}{\rho_{\infty}^{-2/5}}E_{0}^{1/5}$ and
$b_{T}={A^{-1}D_{\mu}^{-2/5}\rho_{\infty}^{-1/5}}{E_{0}^{3/5}}$. From these
scales, we can construct a new scaling form for the various hydrodynamic
fields:

$\displaystyle\rho=t^{-\frac{1}{5}}\tilde{G}(\eta),$ (17a) $\displaystyle
v=t^{-\frac{3}{5}}\tilde{V}(\eta),$ (17b) $\displaystyle
T=t^{-\frac{4}{5}}\tilde{Z}(\eta),$ (17c)

where we have defined the new scaling variable

$\eta=\frac{r}{t^{2/5}},$ (18)

and the constants have been absorbed into the scaling functions $\tilde{G}$,
$\tilde{V}$ and $\tilde{Z}$.

To get the ODEs satisfied by the core scaling functions $\tilde{G}(\eta)$,
$\tilde{V}(\eta)$ and $\tilde{Z}(\eta)$, we plug the scaling ansatz (Eqs.
(17)) into the NS equations and take the long time limit. We get the following
coupled ODEs:

$\displaystyle-\frac{\tilde{G}}{5}-\frac{2\eta\tilde{G}^{\prime}}{5}+\frac{(\eta\tilde{G}\tilde{V})^{\prime}}{\eta}=0,$
(19a) $\displaystyle(\tilde{Z}\tilde{G})^{\prime}=0,$ (19b)
$\displaystyle-\frac{4}{5}\tilde{Z}\tilde{G}-\frac{2\eta\tilde{G}\tilde{Z}^{\prime}}{5}+\tilde{G}\tilde{V}\tilde{Z}^{\prime}+\frac{\tilde{Z}\tilde{G}\tilde{V}}{\eta}+(\tilde{G}\tilde{V}\tilde{Z})^{\prime}=\frac{D_{\kappa}\tilde{Z}^{1/2}\tilde{Z}^{\prime}}{\eta}+\frac{2}{3}(D_{\kappa}\tilde{Z}^{3/2})^{\prime\prime}.$
(19c)

We note the following: (a) The virial corrections do not appear in these
equations because they vanish in the long time limit; (b) only the thermal
conductivity term makes an appearance [this is true also if we consider the
form of transport coefficients in Eq. (29)]. These ODEs are first order in
$\tilde{G}$ and $\tilde{V}$ and second order in $\tilde{Z}$. Thus, we need 4
boundary values: two for each $\tilde{G}$ and $\tilde{V}$, and two for
$\tilde{Z}$. Two of the boundary conditions can be determined using the radial
symmetry of the problem. Because of radial symmetry, there is no flow velocity
and heat flow in the origin. This gives $\tilde{V}(\eta=0)=0$ and
$\tilde{Z}^{\prime}(\eta=0)=0$. The rest two boundary conditions can be
determined from the asymptotic matching condition, which is given by:

$h(\eta\rightarrow\infty)=h(\xi\rightarrow 0),$ (20)

where $h$ is some hydrodynamic field.

We call the solution obtained from the ODEs in Eqs. (6) for the TvNS scaling
functions as the outer solution, and those obtained from the ODEs in Eqs. (19)
for the core scaling functions as the inner solution. We have not succeeded in
finding analytical expressions for the inner solution. However, approximate
solutions can be found in the small $\eta$ region, as we now show. Since
$\tilde{V}(0)=0$ and $\tilde{Z}^{\prime}(0)=0$ we can take for small $\eta$:

$\tilde{V}(\eta)\approx V_{0}\eta^{\alpha}\text{, $\alpha>0$},$ (21)
$\tilde{Z}^{\prime}(\eta)\approx F_{0}\eta^{\beta}\text{, $\beta>0$}.$ (22)

Thus,

$\tilde{Z}(\eta)=F_{0}\frac{\eta^{1+\beta}}{1+\beta}+C_{1},\text{ for small
$\eta$}.$ (23)

Using these and Eq. (19b), we have

$\tilde{G}=\frac{C_{2}}{\tilde{Z}}\approx\frac{C_{2}}{C_{1}}\Big{(}1-\frac{F_{0}\eta^{1+\beta}}{C_{1}(1+\beta)}\Big{)},\text{
for small $\eta$}.$ (24)

Putting these into Eq. (19a) and comparing powers and taking small $\eta$
limit, we get $\alpha=1$, $V_{0}=1/10$. Using these values for $\alpha$ and
$V_{0}$ and Eq. (19c), we get:

(a)

(b)

(c)

(d)

(e)

(f)

Figure 1: MD simulations of hard discs: The three hydrodynamic fields,
obtained from MD simulation data of Ref. joy1 , plotted in (a,b,c) according
to the expected TvNS scaling form [Eqs. (5)] and in (d,e,f) according to the
expected core scaling form [Eqs. (17)]. We can see that in (a,b,c) there is a
very good collapse of data everywhere except near the core, where the scaling
is not very good, especially for the temperature field. On the other hand in
(d,e,f) we see excellent data collapse near the core. Thus MD simulations
exhibit the expected core and bulk scaling forms.

$\displaystyle-\frac{3C_{2}}{5}-\frac{2\eta
C_{2}F_{0}}{5C_{1}}\Big{(}\eta^{\beta}-\frac{F_{0}\eta^{1+2\beta}}{C_{1}(1+\beta)}\Big{)}+\frac{C_{2}F_{0}}{10C_{1}}\Big{(}\eta^{1+\beta}-\frac{F_{0}\eta^{2+2\beta}}{C_{1}(1+\beta)}\Big{)}$
$\displaystyle-
F_{0}D_{\kappa}\sqrt{C_{1}}\Big{(}\eta^{\beta-1}+\frac{F_{0}\eta^{2\beta}}{2C_{1}(1+\beta)}\Big{)}-D_{\kappa}\sqrt{C_{1}}F_{0}\beta\eta^{\beta-1}=0.$
(25)

Since $\beta>0$, the equation above can only be satisfied when $\beta=1$ and

$F_{0}=-\frac{3C_{2}}{10D_{\kappa}\sqrt{C_{1}}}.$ (26)

The only unknown constants are $C_{1}$ and $C_{2}$, which can be fixed from
boundary values determined from either asymptotic matching or the data from
the numerical solution of the full NS Eqs. (14).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 2: NS solution for hard discs: The three hydrodynamic fields, obtained
from numerical solution of NS equations in Eq. (14), are plotted in (a,b,c)
according to the expected TvNS scaling form [Eqs. (5)] and in (d,e,f)
according to the expected core scaling form [Eqs. (17)]. We can see that in
(a,b,c) there is a very good collapse of data everywhere except near the core.
On the other hand in (d,e,f) we see excellent data collapse near the core.
Thus the NS solutions exhibit the expected core and bulk scaling forms. In
(a,b,c) we also see a good fit of the scaled data to the TvNS scaling
functions, $(G,V,Z)$, also obtained from a numerical solution of Eqs. (6). In
(d,e,f) we see a reasonable fit in the core region to the core scaling
functions, $(\tilde{G},\tilde{V},\tilde{Z}$), obtained from a numerical
solution of Eqs. (19) with boundary conditions obtained from the solution of
the NS equations in Eqs. (14). The parameter values were taken as $E_{0}=4.0$,
$\rho_{\infty}=0.382$, and the Henderson equation of state was used.

Thus to first non-zero order in $\eta$, we have:

$\displaystyle\tilde{G}(\eta)=\frac{C_{2}}{C_{1}}\Big{(}1+\frac{3C_{2}\eta^{2}}{20D_{\kappa}C_{1}^{3/2}}+\ldots\Big{)},$
(27a) $\displaystyle\tilde{V}(\eta)=\frac{\eta}{10}+\ldots,$ (27b)
$\displaystyle\tilde{Z}(\eta)=C_{1}-\frac{3C_{2}\eta^{2}}{20D_{\kappa}\sqrt{C_{1}}}+\ldots.$
(27c)

###  2.3 Numerical results

From the discussions in the previous two subsections, we expect two different
scaling form for the hydrodynamic fields, namely the TvNS scaling form given
by Eqs. (5) in the bulk region $X(t)\lesssim r<R(t)$, and the core scaling
form given by Eqs. (17) in the region $0<r\lesssim X(t)$). We now check if
these predictions for the bulk and core scaling forms can be verified in data
obtained from MD simulations and from numerical solutions of the NS equations.

We use the MD simulation results from Ref. joy1 . They considered a gas of
hard discs with diameter $a=1$, $E_{0}=2.0$ and $\rho_{\infty}=0.382$. The
initial condition was chosen with a Gaussian temperature profile,
$v(\boldsymbol{r},t=0)=0$ and $\rho(\boldsymbol{r},t=0)=\rho_{\infty}$. The
parameters of the Gaussian temperature profile were chosen such that the
initial total energy is $E_{0}$. We note here that the long-time scaling
solution does not depend on the width of the Gaussian chosen as long as it
corresponds to energy $E_{0}$. For the numerical solution of the NS equations
in Eqs. (14), we considered initial conditions that correspond to the ones
used in the MD simulations. In all our computations involving hard-disc gas,
we use the Henderson equation of state (EOS) henderson  instead of the
truncated virial expansion for the ease of computation. The Henderson EOS is
given by:

$\displaystyle p=\rho T\frac{128+\pi^{2}\rho^{2}a^{4}}{8(4-\pi\rho
a^{2})^{2}}.$ (28)

The MacCormack method maccormack  was used to solve the PDEs, with
discretization given by $dx=0.05$, $dt=10^{-5}$. The initial temperature
profile was taken to be a Gaussian:
$T(r,0)=({T_{0}}/{\sqrt{2\pi\sigma^{2}}})e^{-r^{2}/2\sigma^{2}}$, with
$T_{0}={E_{0}}/{(\sqrt{2\pi}\rho_{\infty}\sigma)}$ and $\sigma=0.5$. We note
that the lowest order non-zero term of the Enskog expansion for the transport
coefficients gives gass1971 ; garcia2006 ; beijeren ; kremer2010

$\displaystyle\kappa=\frac{2}{a{\pi}^{1/2}}\leavevmode\nobreak\
T^{1/2},\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\mu=\frac{1}{2a{\pi}^{1/2}}\leavevmode\nobreak\ T^{1/2},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\zeta=\frac{8\phi^{2}}{\pi^{3/2}}T^{1/2}=\frac{\pi^{1/2}a^{3}}{2}\leavevmode\nobreak\
\rho^{2}T^{1/2},$ (29)

where $\phi=\pi a^{2}\rho/4$ is the volume fraction. An important general
property of monoatomic gases is the vanishing of the bulk viscosity in the
dilute regime, $\phi\to 0$. This, together with relation $\gamma=1+2/d$, is
the crucial input coming from kinetic theory landau9 ; beijeren ; kremer2010 ;
phenomenological hydrodynamics is compatible with arbitrary $\gamma>1$ and
$\zeta\geq 0$. When $\phi>0$, the bulk viscosity becomes positive, and Eq.
(29) gives its leading asymptotic for small $\phi$. We note here that the
$\rho$-dependence of the bulk viscosity will not change the core scaling, as
the core is a low density region where bulk viscosity is negligible.

In our NS numerics we followed Ref. amit  and chose (with $a=1$) the values
$\kappa=\mu={(\sqrt{\pi}}/{8})\sqrt{T}$ and $\zeta=0$, which are different
from those in Eqs. (29). We have verified that our main conclusions remain
unchanged on using Eqs. (29) for the transport coefficients.

(a)

(b)

(c)

Figure 3: Asymptotic matching in hard discs: Figure comparing the inner
solution and outer solution with the full solution of NS equations. The
boundary values for finding the inner solution were taken from the full
solution of the NS equations. Also shown in the middle plot is the line with
height $0.25$ which is the near core value of $V(\xi)$ predicted by the TvNS
solution. We see that the near core asymptotic of the outer solution matches
with the far core asymptotic of the inner solution, which is nothing but the
rule of asymptotic matching.

MD simulations: In Fig. (1) we show the results of MD simulation data and
check if they satisfy the predicted scaling. In Figs. (1a,1b,1c), following
Eqs. (5), we plot the scaled fields $(\rho,t^{1/2}v,tT)$ as a function of
$r/t^{1/2}$, at different times, and find that this gives a very good data
collapse in the region far from the core. However, as expected, the data
collapse near the core is not quite good. In Figs. (1d,1e,1f), we follow Eqs.
(17) and plot the scaled fields $(t^{1/5}\rho,t^{3/5}v,t^{4/5}T)$ as a
function of $r/t^{2/5}$ and see that this scaling gives a much better data
collapse in this region. This establishes that there are two different regions
with different scaling.

NS equations: We now discuss the results obtained from the numerical solution
of NS equations. In Fig. (2), we plot the hydrodynamic fields obtained from
the numerical solution of the NS equations. In Figs. (2a,2b,2c) we verify that
these satisfy the TvNS scaling form in the bulk region and agree very well
with the TvNS scaling functions. However, in the core we see a clear departure
from TvNS scaling. In Figs. (2d,2e,2f) we see that the core scaling form is
satisfied in the core region. We also plot the solution of Eqs. (19) and we
see that they agree with the numerical solution of NS equations in the core
region.

Asymptotic matching: In Fig. (3), we verify the rules of asymptotic matching
by using the inner and outer solutions. Once we have the inner solution
$\tilde{G}$, $\tilde{V}$ and $\tilde{Z}$, we can find the corresponding
hydrodynamic fields for density, velocity and temperature. We can then compare
how they match with the outer solution. In Fig. (3), we show this comparison
and find that the $\xi\rightarrow 0$ limit of the outer solution has the same
asymptotic as the $\xi\rightarrow\infty$ of the inner solution. Thus we verify
the rule of asymptotic matching.

Finally in Fig. (4), we compare MD data of joy1  (with the same value of
$E_{0}$) with NS solution and TvNS solution. We see that NS and TvNS solutions
agree with MD data. The only discrepancies are in the near core behavior for
the scaling functions.

(a)

(b)

(c)

Figure 4: 2D hard discs: Comparison of the NS solution and TvNS with MD
simulations. We see that the MD data matches with both the NS or TvNS solution
in the bulk region.

##  3 Blast in 2D ideal gas

In this section, we consider the hydrodynamics of a toy model of a 2D gas,
where we ignore the virial corrections in the equation of state, and take it
to be that of an ideal gas. This amounts to considering point particles which
in two dimensions would mean a non-interacting gas. So naively we would not
expect any evolution. However, one can work in the Boltzmann-Grad limit, that
is, $a\to 0$ and $\rho\to\infty$ keeping $\rho a$ finite (in 2D). In this
situation, one expects finite transport coefficients predicted by kinetic
theory, while still preserving the ideal gas equation of state (since the
volume fraction $\pi\rho a^{2}/4\to 0$). This system would be very difficult
to simulate but we can still analyze the NS equations. One advantage of this
toy model is that we can find an exact TvNS solution.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 5: NS solution for 2D hard point gas: The three hydrodynamic fields,
obtained from numerical solution of NS equations in Eq. (14) with ideal gas
equation of state, are plotted in (a,b,c) according to the expected TvNS
scaling form in Eqs. (30) and in (d,e,f) according to the expected core
scaling form [Eqs. (17)]. We can see that in (a,b,c) there is a very good
collapse of data everywhere except near the core. On the other hand in (d,e,f)
we see excellent data collapse near the core. Thus the NS solutions exhibit
the expected core and bulk scaling forms. In (a,b,c) we also see a good fit of
the scaled data to the TvNS scaling functions $(G,V,Z)$, obtained from the
exact solution [ Eqs. (32,34). In (d,e,f) we see a reasonable fit in the core
region to the core scaling functions $\tilde{G},\tilde{V},\tilde{Z}$, obtained
from a numerical solution of Eqs. (19) with boundary conditions obtained from
the solution of the hard point NS equations (PDEs). The parameter values were
taken as $E_{0}=4.0$, $\rho_{\infty}=0.382$.

###  3.1 TvNS solution

For this toy model, our first aim is to find the TvNS solution. To this end we
take the following scaling ansatz for the form of the hydrodynamic fields at
long times:

$\displaystyle\rho=\rho_{\infty}G(\xi),$ (30a) $\displaystyle
v=\frac{r}{t}V(\xi),$ (30b) $\displaystyle T=\frac{r^{2}}{t^{2}}Z(\xi),$ (30c)

where $\xi={r}/{R}$, $R(t)=({E_{0}t^{2}}/{(A_{2}\rho_{\infty})})^{{1}/{4}}$
and $A_{2}$ is some dimensionless constant that we will determine later. Using
these scaling ansatz, the Euler equations can be reduced to the following
equations

$\displaystyle\Big{(}V-\frac{1}{2}\Big{)}G\xi\frac{dV}{d\xi}+\xi\frac{d}{d\xi}(ZG)-GV+GV^{2}+2GZ=0,$
(31a) $\displaystyle\Big{(}V-\frac{1}{2}\Big{)}\xi\frac{dG}{d\xi}+\xi
G\frac{dV}{d\xi}+2GV=0,$ (31b)
$\displaystyle-\Big{(}V-\frac{1}{2}\Big{)}\frac{\xi}{G}\frac{dG}{d\xi}+\Big{(}V-\frac{1}{2}\Big{)}\frac{\xi}{Z}\frac{dZ}{d\xi}+2(V-1)=0.$
(31c)

As a result of the scaling form we see that the total energy in the region
$0<r<\xi R(t)$ is conserved. This leads to the following integral of motion
(see also Ref. landau  for an alternate derivation):

$Z=\frac{V^{2}(1-2V)}{2(4V-1)}.$ (32)

The Rankine-Hugoniot boundary conditions give us

$\displaystyle G(1)=3,\leavevmode\nobreak\ \leavevmode\nobreak\
V(1)=\frac{1}{3},\leavevmode\nobreak\ \leavevmode\nobreak\ Z(1)=\frac{1}{18}.$
(33a)

These equations can be solved to give

$\displaystyle\xi^{4}$
$\displaystyle=\frac{|4V-1|}{108V^{2}\Big{(}V-\frac{1}{2}\Big{)}^{2}},$ (34a)
$\displaystyle G$
$\displaystyle=\frac{3\sqrt{3}\sqrt{|4V-1|}e^{-\frac{1}{2V-1}}}{e^{3}}.$ (34b)

Together with Eq. (32), these provide a complete solution of Eqs. (31), except
that the value of the constant $A_{2}$ is still undetermined. This can be
found easily by using the conservation of energy, exactly as it was done in
chakraborti ; ganapa  for the case of 1D ideal gas:

$A_{2}=2\pi\int_{0}^{1}G(\xi)\Big{(}\frac{V^{2}(\xi)}{2}+Z(\xi)\Big{)}\xi^{3}d\xi.$
(35)

Plugging the exact solution for $G$, $V$, $Z$, we will obtain an expression
for $A_{2}$ as an integral which cannot be evaluated exactly, but can be
evaluated numerically. After doing the full calculation we obtain
$A_{2}=0.3519359068$.

Using the exact solution, one finds that the following $\xi\rightarrow 0$
behaviour of the various scaling functions:

$\displaystyle V-\frac{1}{4}\sim\xi^{4},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ G\sim\xi^{2},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ Z\sim\xi^{-4}.$ (36)

###  3.2 Core scaling solution

Just like we observe a new scaling in the core region for hard discs, we
expect a similar new scaling in the core region for our toy model. Following
the same argument as in Sec. (2.2), and using the fact that the $\xi\to 0$
behaviour of the TvNS solution, Eq. (36), remains the same as for hard discs,
we see immediately that again the size of the dissipative core grows with time
as $t^{2/5}$. It is clear also that the core scaling form will again be
exactly as in Eq. (17). Finally, we note that the higher order terms, in the
virial expansion for pressure in the hard disc gas, do not appear in the
equations for the scaling functions
$\tilde{G}(\eta),\tilde{V}(\eta),\tilde{Z}(\eta)$ in Eq. (19). Hence we have
the same ODE equations in the core as for our 2D point particle gas. As before
we need two boundary conditions to be determined either from NS data, or from
the asymptotic matching condition, which is given by Eqs. (20). The boundary
conditions are different for the hard point gas and will lead to a different
solution for the scaling functions.

###  3.3 Numerical results

In Fig. (5), we compare results from NS solutions with the predicted TvNS and
core scaling forms. In Fig. (5)(a,b,c), we see that the NS shock front matches
with the TvNS prediction. We also see a clear TvNS scaling in the bulk region
but significant departures in the core. Plotting the exact TvNS scaling
functions, we notice that they agree with the bulk NS results. In Fig.
(5)(d,e,f), we verify that the NS solutions satisfy the core scaling. The
comparison between the numerical solution of the core scaling ODEs in Eqs.
(19) and that of the full NS equations is also good in the collapsed region.
For the solution of Eqs. (19), we took the boundary conditions from the full
solution of the NS equations.

In Fig. (6) we plot the exact TvNS solution, the core scaling functions and
the solution from the full NS equations. We see a very good verification of
the rules of asymptotic matching.

(a)

(b)

(c)

Figure 6: 2D hard point gas: In this figure, the inner core scaling solution,
the exact TvNS solution and the solution obtained from the full NS equations
are plotted together. The boundary values for finding the inner solution were
taken from the full solution of the NS equations. We see that the near core
asymptotics of outer solution matches with the far core asymptotics of inner
solution, which is nothing but the rule of asymptotic matching.

##  4 Blast in the hard rod gas

In this section, we consider the blast problem in the 1D hard rod gas, and
compare the predictions of MD simulations with those of NS equations and TvNS
solution.

We consider a gas of hard rods of length $a$. If all the hard rods are of
equal masses, the particles simply exchange velocities during collisions,
hence the system would be integrable and never reach local equilibrium. Thus,
as in chakraborti ; ganapa , we consider the alternate mass hard rod gas where
successive rods on the line have masses $m$ and $M>m$ with mean mass
$\bar{m}={(m+M)}/2=1$. In units where $k_{B}=1$, the equation of state for
this hard rod gas is given by:

$P=\frac{\rho T}{1-\rho a}.$ (37)

###  4.1 TvNS solution

The shock front in 1D grows with time as
$R(t)=\Big{(}\frac{E_{0}t^{2}}{A_{1}\rho_{\infty}}\Big{)}^{\frac{1}{3}}$,
where $\rho_{\infty}$ is the background density. As usual, we choose this
growing length scale to define the TvNS scaling variable $\xi=x/R(t)$ and the
following scaling forms for the hydrodynamic fields:

$\displaystyle\rho=\rho_{\infty}G(\xi),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
v=\frac{x}{t}V(\xi),\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ T=\frac{x^{2}}{t^{2}}Z(\xi).$ (38)

Plugging these into the Euler equations, we get the following ODEs for the
scaling functions:

$\displaystyle-\frac{2}{3}\xi G^{\prime}+\xi(GV)^{\prime}+GV=0,$ (39a)
$\displaystyle-V+V^{2}-\frac{2\xi V^{\prime}}{3}+\xi VV^{\prime}+\frac{\xi
ZG^{\prime}}{G(1-\phi G)^{2}}+\frac{2Z}{1-\phi G}+\frac{\xi Z^{\prime}}{1-\phi
G}=0,$ (39b) $\displaystyle-2Z-\frac{2\xi Z^{\prime}}{3}+2ZV+\xi
VZ^{\prime}+\frac{2ZV}{1-\phi G}+\frac{2\xi ZV^{\prime}}{1-\phi G}=0,$ (39c)

where $\phi=\rho_{\infty}a$ is the background packing fraction. As a result of
the scaling form we see that the total energy in the region $0<r<\xi R(t)$ is
conserved. This leads to the following integral of motion landau :

$\displaystyle Z=\frac{V^{2}(1-3V/2)}{3V/2+3V(1-\phi G)^{-1}-1}.$ (40)

The Rankine-Hugoniot conditions give the following boundary values needed to
solve Eqs. (39):

$\displaystyle G(1)=\frac{2}{1+\phi},\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ V(1)=\frac{1-\phi}{3},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ Z(1)=\frac{(1-\phi)^{2}}{9}.$ (41)

These equations have to be solved along with the equation which specifies the
total energy, which will fix the constant $A_{1}$. For point particles, Eq.
(40) makes it possible to solve the TvNS equations in any dimensions. However,
for the case of hard rods we have not been able to find a closed form
solution. We have verified numerically that the numerical solution of the ODEs
above have the same behaviour near the core as the TvNS solution of ideal gas
obtained in chakraborti ; ganapa , i.e.,

$\displaystyle G(\xi)\sim|\xi|^{1/2},\leavevmode\nobreak\ \leavevmode\nobreak\
V(\xi)=2/9,\leavevmode\nobreak\ \leavevmode\nobreak\ Z(\xi)\sim|\xi|^{-5/2},$
(42)

for small $\xi$. This is expected as the small $\xi$ region is a region of low
density making the virial corrections almost negligible and the gas ideal.

###  4.2 Core scaling solution

Following Refs. chakraborti ; ganapa , we take for the heat conductivity in NS
to be of the form $\kappa=D_{\kappa}\rho^{1/3}T^{1/2}$. An analysis similar to
the one we did in Sec. (2.2) leads us then to a core growing as $X(t)\sim
t^{38/93}$ and the following forms for the scaling near the core:

$\displaystyle\rho=t^{-\frac{4}{31}}\tilde{G}(\eta),$ (43a) $\displaystyle
v=t^{-\frac{55}{93}}\tilde{V}(\eta),$ (43b) $\displaystyle
T=t^{-\frac{50}{93}}\tilde{Z}(\eta),$ (43c)

where $\eta=xt^{-38/93}$.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 7:  Hard rods: Plot checking the TvNS scaling in (a,b,c) and the new
scaling in the core (obtained by including the dissipative effects) in (d,e,f)
from the results of MD simulation. We can see that there is a nice collapse of
data at two different times in the bulk region in (a,b,c) and in the core
region in (d,e,f). We have taken $\rho_{\infty}=0.4$, $E_{0}=0.4$.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 8:  Hard rods: Plot checking the TvNS scaling in (a,b,c) and the new
scaling in the core (obtained by including the dissipative effects) in (d,e,f)
from numerical solution of the NS equations. Also shown are the outer solution
in (a,b,c) and the inner solution in (d,e,f). We can see that there is a nice
collapse of data at two different times in the bulk region in (a,b,c) and in
the core region in (d,e,f). We can also see that the outer solution matches
with the NS data in the bulk region in (a,b,c) and the inner solution matches
with the NS data in the core region in (d,e,f). We have taken
$\rho_{\infty}=0.4$, $E_{0}=0.4$.

We can solve for the core scaling functions similar to the way we did it for
hard discs. Since the virial corrections to the equation of state are
negligible in the core and the gas is almost ideal, we directly take the
equations for the core scaling functions from ganapa :

$\displaystyle-\frac{4}{31}\tilde{G}-\frac{38}{93}\eta\tilde{G}^{\prime}+(\tilde{G}\tilde{V})^{\prime}=0,$
(44a) $\displaystyle(\tilde{G}\tilde{Z})^{\prime}=0,$ (44b)
$\displaystyle-\frac{25}{93}\tilde{G}\tilde{Z}-\frac{19}{93}\eta\tilde{G}\tilde{Z}^{\prime}+\frac{\tilde{G}\tilde{V}\tilde{Z}^{\prime}}{2}+\tilde{G}\tilde{Z}\tilde{V}^{\prime}=(\tilde{G}^{1/3}\tilde{Z}^{1/2}\tilde{Z}^{\prime})^{\prime}.$
(44c)

As for the 2D case in Eqs. (19) we observe again that the change in the
equation of state does not affect these equations and only the thermal
conductivity term makes an appearance in the core scaling equations.

###  4.3 Numerical results

In our simulations, we have taken $m=4/5$ and $M=6/5$. We verify the core
scaling and bulk scaling from MD data in Figs. (7).

We solved the TvNS ODEs numerically, and we solved the NS equations following
the method of chakraborti ; ganapa . The plots comparing the two solutions are
shown in Figs. (8a, 8b, 8c). We can see that the shock fronts predicted by the
two solutions matches. Also there is a nice data collapse in the bulk for the
NS solutions at two different times.

We solve Eqs. (44) (with boundary conditions taken from the NS data) and
compare the resulting solution with the full numerical solution of the NS
equations in Figs. (8d, 8e, 8f). We see that there is agreement in the core,
and that there is data collapse in the core of the NS solutions at two
different times.

While solving the NS equations, we take (following chakraborti ; ganapa ) heat
conductivity $\kappa=D_{\kappa}\rho^{1/3}T^{1/2}$, with $D_{\kappa}=1$ and
also a finite bulk viscosity $\zeta=D_{\zeta}T^{1/2}$, with $D_{\zeta}=1$
(setting this to zero does not change any of our conclusions). We also
compared our NS numerical solution and TvNS solution with MD simulations (Fig.
(9)) and find that the shock front position agrees quite well with that
obtained from MD simulations. Near the core, there are some discrepancies,
possibly due to anomalous heat conduction in 1D, or due to large deviations
from local equilibrium (also observed for the case for point particles in
chakraborti ; ganapa ).

(a)

(b)

(c)

Figure 9: Hard rods Plot comparing the NS solution and TvNS solution with the
MD simulations. We see that the shock front matches. We have taken
$\rho_{\infty}=0.4$, $E_{0}=0.4$.

##  5 Hard point gas in $d-$dimensions: Inner and outer solutions

In this section, we briefly discuss some features of the inner and outer
solutions for the point particle gas in arbitrary dimensions. As in Sec. (3),
we consider the limit of $\rho_{\infty}\to\infty$ and $a\to 0$ such that one
has an ideal gas with finite transport coefficients. We first discuss the
outer solution given by the TvNS-type analysis of the Euler equations. Then we
discuss the inner solution where one has to write the full NS equations
including dissipation terms.

###  5.1 Outer solution

####  5.1.1 Euler equations in $d-$dimensions

From dimensional analysis we have in this case the shock front position
$R(t)=\left(\frac{E_{0}t^{2}}{A_{d}\rho_{\infty}}\right)^{\frac{1}{d+2}}$.
Behind the shock, $0\leq r<R(t)$, the radial velocity $v(r,t)$, density
$\rho(r,t)$ and pressure $P(r,t)$ satisfy the Euler equations

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

where $\gamma$ is the adiabatic index.

The velocity of the shock wave is

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

The Rankine-Hugoniot conditions landau  describing the jump of the
hydrodynamic variables on both sides of the shock wave read

$\frac{P(R)}{\rho_{\infty}U^{2}}=\frac{2}{\gamma+1}\,,\leavevmode\nobreak\
\leavevmode\nobreak\
\frac{\rho(R)}{\rho_{\infty}}=\frac{\gamma+1}{\gamma-1}\,,\leavevmode\nobreak\
\leavevmode\nobreak\ \frac{v(R)}{U}=\frac{2}{\gamma+1},$ (47)

in the case of the infinitely strong blast.

####  5.1.2 Self-similar solution

Dimensional analysis assures that the hydrodynamic variables acquire a self-
similar form

$v=\delta\,\frac{r}{t}\,V,\quad\rho=\rho_{\infty}\,G,\quad
c^{2}=\delta^{2}\,\frac{r^{2}}{t^{2}}\,Z.$ (48)

Here $c^{2}=\gamma P/\rho=\gamma T$ is the square of the speed of sound and we
have introduced constant factors of $\delta$ and $\delta^{2}$ in the
definition of the scaling functions to simplify subsequent computations. The
dimensionless quantities $V,G,Z$ depend on the single variable

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

One seeks the behavior of $V(\xi),G(\xi)$ and $Z(\xi)$ behind the shock wave,
$0\leq\xi\leq 1$. The Rankine-Hugoniot conditions Eqs. (47) become

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

Using energy conservation, one expresses $Z$ through the scaled velocity $V$
to give landau :

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

Plugging the ansatz Eqs. (48–49) into Eq. (45a) gives

$\frac{dV}{d\ell}+(V-1)\,\frac{d\ln G}{d\ell}=-dV,$ (52a) where $\ell=\ln\xi$.
Similarly we transform Eq. (45b) into
$\frac{d\ln(Z/G^{\gamma-1})}{d\ell}=\frac{d+2-2V}{V-1}.$ (52b)

Equations (52a–52b) can be solved for arbitrary $d$ and $\gamma>1$. Having
solved the problem, one can compute the energy

$\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]$
$\displaystyle=$ $\displaystyle
n_{\infty}\,\Omega_{d}\,\delta^{2}\,\frac{R^{d+2-a}}{t^{2}}\int_{0}^{1}d\xi\,\xi^{d+1}\,\frac{(\gamma-1)V^{3}}{2(\gamma
V-1)}\,G.$

Here $\Omega_{d}=2\pi^{d/2}/\Gamma(d/2)$ is the surface area of the unit
sphere and we have used Eq. (51). Substituting Eq. (1) into the above
expression for the energy we fix the dimensionless constant $A_{d}$:

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

In the following we limit ourselves to a monoatomic gas where $\gamma=1+2/d$.
The adiabatic index generally depends only on the spatial dimension and
independent of the interaction between particles fluid ; Cercignani . The
monoatomic gases are relevant in astrophysics where even diatomic molecules
are extremely rare. For the monoatomic gas, the Rankine-Hugoniot conditions
Eqs. (50a–50c) become

$V(1)=\frac{d}{d+1}\,,\quad G(1)=d+1,\quad Z(1)=\frac{d+2}{(d+1)^{2}},$ (54)

and Eq. (51) reduces to

$Z=\frac{(1+2/d)(1-V)V^{2}}{(d+2)V-d}.$ (55)

Using Eq. (52a) and Eq. (52b) with $\gamma=1+2/d$ and $Z$ given by Eq. (55) we
express $dV/d\ell$ and $d(\ln G)/d\ell$ through the scaled velocity $V$:

$\displaystyle\frac{dV}{d\ell}=V\,\frac{X}{2D},$ (56a)
$\displaystyle\frac{d\ln G}{d\ell}=\frac{V}{V-1}\,\frac{Y}{2D},$ (56b)

where we write in short

$\begin{split}X&=[(d+2)V-d](d+2-4V),\\\ D&=d-2(d+1)V+(d+1)(1+2/d)V^{2},\\\
Y&=(2-d)d+(d-2)(2+3d)V+2(d+2)(1-d)V^{2}.\end{split}$

Integrating Eq. (56a) yields an implicit solution for the scaled velocity

$\xi=c_{d}[(d+2)V-d]^{\frac{2}{d^{2}+4}}\,V^{-\frac{2}{d+2}}\,(d+2-4V)^{-\lambda_{d}},$
(57)

with

$\begin{split}c_{d}&=\left(\frac{d}{d+1}\right)^{\frac{2}{d+2}-\frac{2}{d^{2}+4}}\,\left(\frac{d^{2}-d+2}{d+1}\right)^{\lambda_{d}},\\\
\lambda_{d}&=\frac{8+4d+10d^{2}-d^{3}+d^{4}}{2d(2+d)(4+d^{2})}.\end{split}$
(58)

We then divide Eq. (56b) by Eq. (56a) and integrate in $V$ to give

$G=C_{d}\,[(d+2)V-d]^{\frac{d^{2}}{d^{2}+4}}\,(1-V)^{-\frac{2d}{d-2}}\,(d+2-4V)^{\Lambda_{d}},$
(59)

with

$\begin{split}C_{d}&=(d+1)^{\frac{d^{2}}{d^{2}+4}-\frac{2d}{d-2}}\,\,d^{-\frac{d^{2}}{d^{2}+4}}\,\left(\frac{d^{2}-d+2}{d+1}\right)^{-\Lambda_{d}},\\\
\Lambda_{d}&=\frac{8+4d+10d^{2}-d^{3}+d^{4}}{2(d-2)(4+d^{2})}.\end{split}$
(60)

Finally, $Z=Z(V)$ is given by Eq. (55).

####  5.1.3 Singular behavior

Using Eq. (57) and Eq. (59), one extracts the asymptotic behaviors near the
center of the explosion ($\xi\ll 1$):

$G\sim\xi^{\frac{d^{2}}{2}}\,,\quad Z^{-1}\sim
V-\tfrac{d}{d+2}\sim\xi^{\frac{d^{2}+4}{2}}.$ (61)

In terms of the original coordinates

$\rho\sim r^{\frac{d^{2}}{2}}\,t^{-\frac{d^{2}}{d+2}}\,,\qquad T\sim
r^{-\frac{d^{2}}{2}}\,t^{\frac{d(d-2)}{d+2}}.$ (62)

This physically dubious behavior (the density vanishes, while the temperature
diverges at $r=0$) indicates that the framework we have employed so far
becomes incomplete near the center of the explosion. We rectify it by
including dissipation.

###  5.2 Inner solution

####  5.2.1 Navier-Stokes equations

Here we are considering the hard-point gas in the Boltzmann-Grad limit. Notice
however that, even if we were to take the model of hard sphere gas to ensure
that there are collisions, our core analysis would not change. This would mean
virial corrections in the equation of state and corrections to the transport
coefficients. However, since the core is a region of low density, these
corrections are negligible, and the core equations would be the same as that
of an ideal gas.

The continuity equation is not affected by dissipative processes landau . Thus
for our radial flow Eq. (45a) remains valid. The momentum equation is
generally affected by viscosity—adding an extra term
$\nabla\cdot\boldsymbol{\sigma}$ to the Euler equation yields the Navier-
Stokes equation. The viscous stress tensor is given by

$\boldsymbol{\sigma}=\mu\left[\nabla{\bf v}+(\nabla{\bf
v})^{T}-\tfrac{2}{d}(\nabla\cdot{\bf v}){\bf I}\right]+\zeta(\nabla\cdot{\bf
v}){\bf I},$

where ${\bf I}$ is the unit tensor and the coefficients of shear and bulk
viscosity are denoted by $\mu$ and $\zeta$. Little is known about bulk
viscosity. Fortunately, for dilute monoatomic gases the coefficient of bulk
viscosity vanishes fluid . Thus for the dilute hard sphere (HS) gas

$\boldsymbol{\sigma}=\mu\left[\nabla{\bf v}+(\nabla{\bf
v})^{T}-\tfrac{2}{d}(\nabla\cdot{\bf v}){\bf I}\right].$ (63)

According to kinetic theory fluid ; Cercignani , transport coefficients
exhibit similar behavior. For the dilute HS gas, these coefficients are
proportional to $\sqrt{T}$. This is the classical prediction of kinetic
theory. More precisely, the coefficient of shear viscosity for the dilute hard
sphere (HS) gas reads

$\mu=M_{d}\,a^{-(d-1)}\sqrt{mT},$ (64)

where $m$ and $a$ denote spheres’ masses and radii. The dimensionless
amplitude $M_{3}$ is known only approximately fluid ; Cercignani . The
coefficient of thermal conductivity for the dilute HS gas is given by a
similar formula

$\kappa=K_{d}\,a^{-(d-1)}\sqrt{mT},$ (65)

which differs from the expression Eq. (64) only by the dimensionless
amplitude, $K_{d}$ instead of $M_{d}$.

The entropy equation is affected by heat conduction and viscous dissipation.
To estimate the size of the core region where dissipative effects matter it
suffices to keep heat conduction and ignore viscous dissipation. In this
situation Eq. (45b) is replaced by equation

$\rho
T(\partial_{t}+v\partial_{r})s=r^{-(d-1)}\,\partial_{r}(r^{d-1}\kappa\partial_{r}T),$
(66)

for the entropy per unit mass

$s=\frac{d}{2}\,\ln\\!\left(\frac{T}{\rho^{2/d}}\right).$ (67)

Using Eqs. (67–65) we recast Eq. (66) into

$\rho^{1+2/d}\,(\partial_{t}+v\partial_{r})\frac{T}{\rho^{2/d}}=\frac{1}{r^{d-1}}\,\partial_{r}\big{(}r^{d-1}T^{1/2}\,\partial_{r}T\big{)}.$
(68)

We have set $K_{d}a^{-(d-1)}\sqrt{m}=1$ in Eq. (68) to avoid cluttering the
formulas; this factor can be always restored in the final results on
dimensional grounds.

####  5.2.2 Heuristic analysis for core scaling

Using 68 we deduce an estimate

$\rho\,\frac{T}{t}\sim\frac{T^{3/2}}{r^{2}},$ (69)

which is combined with 62 to give an estimate of the growth of the radius $X$
of the core where the heat transfer plays significant role

$X\propto R^{h_{d}}\,,\quad h_{d}=\frac{4+3d^{2}}{8+3d^{2}}.$ (70)

According to Eq. (62), the density vanishes and the temperature diverges at
the center of the explosion. Heat conduction rectifies these predictions.
Indeed, instead of using Eq. (62) at $r=0$, one should substitute $r=X$ into
Eq. (62). This allows us to estimate the density, temperature and pressure in
the core region including the center of the explosion

$\rho^{*}\propto R^{-\nu_{d}}\,,\quad T^{*}\propto R^{-(d-\nu_{d})}\,,\quad
p^{*}\propto R^{-d},$ (71)

with $\nu_{d}=2d^{2}/(8+3d^{2})$.

##  6 Discussion

In this work, we presented the results of numerical solutions of the Navier-
Stokes (NS) equations of dissipative hydrodynamics for a 2D hard disc gas, a
2D hard point gas, and a 1D hard rod gas. We revealed that the long-time
solution has a double scaling form — it consists of an outer solution
described by the well-known TvNS scaling form, and an inner solution described
by different scaling functions. For the 2D gas, the inner solution is valid in
the region $0<r\lesssim t^{2/5}$ while the outer solution is valid in the
region $t^{2/5}\lesssim r<R(t)$, with $R(t)\sim t^{1/2}$ known quite
precisely. The inner solution arises due to the dominance of dissipative terms
(namely, the viscosity and thermal conductivity terms in the NS equations) in
the core of the blast.

We also pointed out that molecular dynamics (MD) simulation data for both the
1D and 2D gases agree with the results from the solution of NS equations,
except in the core. The MD data does satisfy the scaling forms corresponding
to the inner and outer scaling solutions, however, the scaling functions
differ. This mismatch in the core occurs also for 1D point particles
chakraborti ; ganapa  and a possible reason could be that transport
coefficients are anomalous in low dimensional systems.

Note that the ratio of the size of the core to the size of the blast
$X(t)/R(t)$ in 1D (${t^{38/93}}/{t^{2/3}}\sim t^{-24/93}$ chakraborti ; ganapa
) is much smaller than in 2D (${t^{2/5}}/{t^{1/2}}\sim t^{-1/10}$). This ratio
is even larger in 3D (${t^{62/175}}/{t^{2/5}}$). Despite the ratio being
larger in 2D than in 1D, we find that dissipative corrections in the core do
not affect the position of the shock front in 2D. In general $d$-dimensions
($d\geq 2$), the shock front position, $R(t)$, and the core size, $X(t)$, grow
as [see Sec. (5) for the details of the derivation]:

$R(t)\sim t^{\frac{2}{(d+2)}},\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ X(t)\sim
t^{\frac{2(4+3d^{2})}{(d+2)(8+3d^{2})}}.$ (72)

We present below a table showing the size of the shock front and the core, the
hydrodynamic fields in the core and the ratio of the core energy ($E^{*}$) to
the total energy in 2D and 3D [the details of the derivation are given in Sec.
(5)].

$d$ | $R$ | $X$ | $\rho^{*}$ | $T^{*}$ | $p^{*}$ | $v^{*}$ | $E^{*}/E_{0}$
---|---|---|---|---|---|---|---
2 |  $t^{1/2}$ | $t^{2/5}$ | $t^{-1/5}$ | $t^{-4/5}$ | $t^{-1}$ | $t^{-3/5}$ | $t^{-1/5}$
3 |  $t^{2/5}$ |  $t^{62/175}$ | $t^{-36/175}$ | $t^{-174/175}$ | $t^{-6/5}$ | $t^{-113/175}$ | $t^{-24/175}$

The core is the region where deviations from local equilibrium are large (see
also joy1 ), and treating the deviation as a first order approximation (which
is what is done in the NS framework) may not suffice. Apart from this, there
is the issue of long time tails of correlation functions and the divergence of
transport coefficients in dimensions $d\leq 2$ alder1970 ; beijeren ; dhar2008
; lepri2016 . This can be also be a reason for the mismatch between
hydrodynamics and microscopic simulations observed in the core. Finally,
another reason for the failure of hydrodynamics in the core could be the fact
that the Knudsen number (which is the ratio of the mean free path to the
length scale of variation of the hydrodynamic fields) in the core for the 2D
case has a very slow decay than that in the bulk. In the core, the Knudsen
number is $\sim{\rho^{*}}^{-1}t^{-2/5}\sim t^{-1/5}$, while that in the bulk
is $\sim\rho_{\infty}^{-1}t^{-1/2}\sim t^{-1/2}$. This reason is connected to
the reason mentioned before, as the Knudsen number is a measure of the
deviation from local equilibrium and how well the hydrodynamic limit is
reached.

Interesting open problems would be to resolve the observed disagreement
between hydrodynamics and microscopic dynamics in the core for the blast
problem, and to see if this disagreement disappears in higher dimensions.

## Conflict of interest statement

The authors declare that there is no conflict of interest.

## Data Availability

Data will be made available on reasonable request.

## Acknowledgments

We thank Amit Kumar and R. Rajesh for discussions and for permission to use
data from Ref. joy1 . S.K.S and A. D. acknowledge support from the Department
of Atomic Energy, Government of India, under project no. 19P1112R&D. AD and
PLK acknowledge the support of the Erwin Schrodinger Institute where several
discussions related to the work were held during the program ‘Large
Deviations, Extremes and Anomalous Transport in Non-equilibrium Systems’. The
numerical computations were done on ICTS clusters Contra, Tetris and Mario.

## References

  * [1] E. M. Lifshitz and L. P. Pitaevskii.  Physical Kinetics. Course of Theoretical Physics, Vol. 10.  Butterworth-Heinemann, 1981.
  * [2] Jeremiah P. Ostriker and Christopher F. McKee.  Astrophysical blastwaves.  Rev. Mod. Phys., 60:1–68, Jan 1988.
  * [3] Jeremy Goodman.  Convective instability of hollow sedov-taylor blast waves.  Astrophys. J., 358:214, jul 1990.
  * [4] A. Gal-Yam, D. B. Fox, P. A. Price, E. O. Ofek, M. R. Davis, D. C. Leonard, A. M. Soderberg, B. P. Schmidt, K. M. Lewis, B. A. Peterson, S. R. Kulkarni, E. Berger, S. B. Cenko, R. Sari, K. Sharon, D. Frail, D. S. Moon, P. J. Brown, A. Cucchiara, F. Harrison, T. Piran, S. E. Persson, P. J. McCarthy, B. E. Penprase, R. A. Chevalier, and A. I. MacFadyen.  A novel explosive process is required for the $\gamma$-ray burst GRB 060614.  Nature, 444(7122):1053–1055, 2006.
  * [5] Xiaping Tang and Roger A. Chevalier.  Shock evolution in non-radiative supernova remnants.  Monthly Notices of the Royal Astronomical Society, 465(4):3793–3802, 11 2016.
  * [6] Adam Burrows.  Colloquium: Perspectives on core-collapse supernova theory.  Rev. Mod. Phys., 85:245–261, Feb 2013.
  * [7] L. I. Sedov.  Similarity and dimensional methods in mechanics.  Academic Press, New York, 1959.
  * [8] L. D. Landau and E. M. Lifshitz.  Fluid Mechanics.  Pergamon Press, New York, 1987.
  * [9] G. I. Taylor.  The formation of a blast wave by a very intense explosion I. Theoretical discussion.  Proc. R. Soc. A, 201(1065):159–174, 1950.
  * [10] Geoffrey Ingram Taylor.  The formation of a blast wave by a very intense explosion.-ii. the atomic explosion of 1945.  Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 201(1065):175–186, 1950.
  * [11] J. von Neumann.  Collected Works.  Pergamon Press,Oxford, 1963.
  * [12] Leonid Ivanovich Sedov.  Propagation of strong shock waves.  Journal of Applied Mathematics and Mechanics, 10:241–250, January 1946.
  * [13] H A Bethe, K Fuchs, J von Neuman, R Peierls, W G Penney, and notes written by J.O. Hirschfelder.  Shock hydrodynamics and blast waves.  10 1944.
  * [14] J. von Neumann.  Physics of shock waves and high-temperature hydrodynamic phenomena, Vol. 1-2.  Academic Press, New York, 1966-67.
  * [15] A. F. Ghoniem, M. M. Kamel, S. A. Berger, and A. K. Oppenheim.  Effects of internal heat transfer on the structure of self-similar blast waves.  Journal of Fluid Mechanics, 117:473–491, 1982.
  * [16] H. Steiner and W. Gretler.  The propagation of spherical and cylindrical shock waves in real gases.  Physics of Fluids, 6(6):2154–2164, 1994.
  * [17] J. VonNeumann and R. D. Richtmyer.  A method for the numerical calculation of hydrodynamic shocks.  Journal of Applied Physics, 21(3):232–237, 1950.
  * [18] Harold L. Brode.  Numerical solutions of spherical blast waves.  Journal of Applied Physics, 26(6):766–775, 1955.
  * [19] Richard Latter.  Similarity solution for a spherical shock wave.  Journal of Applied Physics, 26(8):954–960, 1955.
  * [20] Myron N. Plooster.  Shock waves from line sources. numerical solutions and experimental measurements.  The Physics of Fluids, 13(11):2665–2675, 1970.
  * [21] Subhadip Chakraborti, Santhosh Ganapa, P. L. Krapivsky, and Abhishek Dhar.  Blast in a one-dimensional cold gas: From newtonian dynamics to hydrodynamics.  Phys. Rev. Lett., 126:244503, Jun 2021.
  * [22] Santhosh Ganapa, Subhadip Chakraborti, P. L. Krapivsky, and Abhishek Dhar.  Blast in the one-dimensional cold gas: Comparison of microscopic simulations with hydrodynamic predictions.  Physics of Fluids, 33:087113, 2021.
  * [23] Subhadip Chakraborti, Abhishek Dhar, and Paul Krapivsky.  A splash in a one-dimensional cold gas.  SciPost Physics, 13(3):074, 2022.
  * [24] Stefano Lepri, Roberto Livi, and Antonio Politi.  Thermal conduction in classical low-dimensional lattices.  Physics reports, 377(1):1–80, 2003.
  * [25] Abhishek Dhar.  Heat transport in low-dimensional systems.  Advances in Physics, 57(5):457–537, 2008.
  * [26] Aritra Kundu, Cédric Bernardin, Keji Saito, Anupam Kundu, and Abhishek Dhar.  Fractional equation description of an open anomalous heat conduction set-up.  Journal of Statistical Mechanics: Theory and Experiment, 2019(1):013205, jan 2019.
  * [27] Abhishek Dhar, Anupam Kundu, and Aritra Kundu.  Anomalous heat transport in one dimensional systems: A description using non-local fractional-type diffusion equation.  Frontiers in Physics, 7, 2019.
  * [28] Pablo I. Hurtado.  Breakdown of hydrodynamics in a simple one-dimensional fluid.  Phys. Rev. Lett., 96:010601, Jan 2006.
  * [29] M. Barbier, D. Villamaina, and E. Trizac.  Microscopic origin of self-similarity in granular blast waves.  Physics of Fluids, 28:083302, 2016.
  * [30] Jilmy P. Joy and R. Rajesh.  Shock Propagation in the Hard Sphere Gas in Two Dimensions: Comparison Between Simulations and Hydrodynamics.  Journal of Statistical Physics, 184(1):3, July 2021.
  * [31] Jilmy P. Joy, Sudhir N. Pathak, and R. Rajesh.  Shock Propagation Following an Intense Explosion: Comparison Between Hydrodynamics and Simulations.  Journal of Statistical Physics, 182(2):34, February 2021.
  * [32] A. Kumar and R. Rajesh.  Blast waves in two and three dimensions: Euler versus navier–stokes equations.  Journal of Statistical Physics, 188(2), jun 2022.
  * [33] A. Kumar.  Private communication.
  * [34] J. R. Dorfman, Henk van Beijeren, and T. R. Kirkpatrick.  Contemporary Kinetic Theory of Matter.  Cambridge University Press, 2021.
  * [35] Alain Barrat and Emmanuel Trizac.  Molecular dynamics simulations of vibrated granular gases.  Phys. Rev. E, 66:051303, Nov 2002.
  * [36] R. W. MacCormack.  A numerical method for solving the equations of compressible viscous flow.  AIAA Journal, 20(9):1275–1281, 1982.
  * [37] David M Gass.  Enskog theory for a rigid disk fluid.  The Journal of Chemical Physics, 54(5):1898–1902, 1971.
  * [38] Ramón García-Rojo, Stefan Luding, and J Javier Brey.  Transport coefficients for dense hard-disk systems.  Physical Review E, 74(6):061305, 2006.
  * [39] Gilberto M Kremer.  An introduction to the Boltzmann equation and transport processes in gases.  Springer Science & Business Media, 2010.
  * [40] P. Resibois and M. De Leener.  Classical Kinetic Theory of Fluids.  Wiley, New York, 1977.
  * [41] C. Cercignani.  The Boltzmann equation and its applications.  Springer, New York, 1988.
  * [42] BJ Alder and TE Wainwright.  Decay of the velocity autocorrelation function.  Physical review A, 1(1):18, 1970.
  * [43] Stefano Lepri.  Thermal transport in low dimensions: from statistical physics to nanoscale heat transfer, volume 921.  Springer, 2016.


$B_{3}\subseteq\mathrm{interior}(V)$, $K:=\max_{B_{3}}\left\lvert
f\right\rvert$ exists and $f$ is $(\frac{1}{\delta}\max_{B_{3}}\left\lvert
f\right\rvert)$-Lipschitz on $B_{1}$. A $\delta$ satisfying the first two
conditions exists since $V$ is open and bcause $f$ is locally bounded and
$B_{3}$ is compact. We will show that for all such $\delta$, the third
condition already holds.

For that, fix $y,z\in B_{1}$ and call $r:=\frac{d(y,z)}{2\delta}\in[0,1]$. We
want to show that $\left\lvert
f(y)-f(z)\right\rvert\leq\frac{K}{\delta}d(y,z)=2Kr$. Wlog $y\neq z$,
otherwise there is nothing to show.

So $r\neq 0$ and we can consider $u:=\frac{1+r}{r}z-\frac{1}{r}y$,
$v:=\frac{1}{r}y-\frac{r-1}{r}z$. We have
$u,v\in\overline{B_{2\delta}(z)}=:B_{2}$. Indeed, $d(u,z)=\lVert
u-z\rVert=\lVert\frac{z}{r}+z-\frac{y}{r}-z\rVert=\frac{\lVert
z-y\rVert}{r}=2\delta$ and similarly $d(v,z)=2\delta$. Geometrically, those
are actually the intersections between $B_{2}$ and the line generated by $y$
and $z$, see Figure 4. Now we have the convex combinations
$z=\frac{1}{1+r}y+\frac{r}{1+r}u$ and $y=(1-r)z+rv$, so the concavity of $f$
entails on one hand:
$f(z)\geq\frac{1}{1+r}f(y)+\frac{r}{1+r}f(u)\geq\frac{f(y)}{1+r}-\frac{rK}{1+r}$,
i.e. $f(y)-f(z)\leq r(K+f(z))\leq 2rK$, and on the other hand:
$f(y)\geq(1-r)f(z)+rf(v)\geq f(z)-r(f(z)+K)$, i.e. $f(z)-f(y)\leq
r(f(z)+K)\leq 2rK$. In the previous inequalities we have used that
$f(u),f(v)\geq-K$. This follows because $u,v\in B_{2}\subseteq B_{3}$, as it
can be immediately checked, thus $\left\lvert f(u)\right\rvert,\left\lvert
f(v)\right\rvert\leq K$. Putting the final inequalities together, we have
$\left\lvert f(y)-f(z)\right\rvert\leq 2rK$, i.e. the thesis. ∎

Therefore, since $(0,\infty)^{X}$ is open and convex in
$(\mathbb{R}^{X},\lVert\cdot\rVert)$ and all tps are non-negative, we
immediately have Theorem 6.3.

##  Appendix D Proofs from Section 7: Lipschitz meets Taylor

###  D.1 Proof of Theorem 7.1

Theorem 7.1 states the validity of Taylor expansion in
$\mathbb{L}\mathsf{Rel}_{!}$. We must check for it the following equation,
given $f\in\mathbb{L}\mathsf{Rel}_{!}(C,B^{A})$ and
$g\in\mathbb{L}\mathsf{Rel}_{!}(C,A)$:

$\mathsf{ev}\circ\langle
f,g\rangle=\inf_{n\in\mathbb{N}}\left\\{((\cdots(\Lambda^{-}(f)\underbrace{\star
g)\cdots)\star g}_{n\text{
times}})\circ\langle\mathrm{id},\infty\rangle\right\\}$

where:

  1. 1.

$\mathsf{ev}\in\mathbb{L}\mathsf{Rel}_{!}(B^{A}+A,B)$ is the canonical
_evaluation_ morphism;

  2. 2.

$\Lambda^{-}(\\_):=\mathsf{ev}\circ(\\_\times\mathrm{id})$ is the _uncurry_
operator;

  3. 3.

given $f\in\mathbb{L}\mathsf{Rel}_{!}(C+A,B)$ and
$g\in\mathbb{L}\mathsf{Rel}_{!}(C,A)$, $f\star
g\in\mathbb{L}\mathsf{Rel}_{!}(C+A,B)$ is the morphism obtained by
differentiating $f$ in its second component and applying $g$ in that
component, i.e.

$f\star
g=D(f)\circ\langle\langle\infty,g\circ\pi_{1}\rangle,\mathrm{id}_{C+A}\rangle.$

We do it in the following $4$ steps.

  1. 1.

Let us compute the morphism $\mathsf{ev}$ explicitly:
$\mathsf{ev}\in\mathbb{L}^{\mathcal{M}_{\mathsf{fin}}((\mathcal{M}_{\mathsf{fin}}(A)\times
B)+A)\times B}$ is given by

$\mathsf{ev}_{\mu,y}=\begin{cases}0&\text{ if
}\mu=[\langle\rho,y\rangle]\oplus\rho\\\ \infty&\text{ otherwise}\end{cases}$

and observe that, given $f\in\mathbb{L}\mathsf{Rel}_{!}(C,B^{A})$ and
$g\in\mathbb{L}\mathsf{Rel}_{!}(C,A)$,

$\big{(}\mathsf{ev}\circ\langle
f,g\rangle\big{)}_{\chi,y}=\inf\Big{\\{}\sum_{i=1}^{m}g_{\chi_{i},x_{i}}+f_{\chi^{\prime},\langle[x_{1},\dots,x_{m}],y\rangle}\mid
x_{1},\dots,x_{m}\in A,\chi=\chi^{\prime}+\sum_{i=1}^{m}\chi_{i},\Big{\\}}$

  2. 2.

Let us compute the morphism $\Lambda^{-}$ explicitly: given
$g\in\mathbb{L}\mathsf{Rel}_{!}(C,B^{A})$,
$\Lambda^{-}(g)\in\mathbb{L}\mathsf{Rel}_{!}(C+A,B)$ is given by

$\big{(}\Lambda^{-}(g)\big{)}_{\rho\oplus\mu,y}=g_{\rho,\langle\mu,y\rangle}$

  3. 3.

Let us compute the morphism $\star$ explicitly: $f\star g$ is given by

$(f\star
g)_{\rho\oplus\mu,y}=\inf\Big{\\{}g_{\rho^{\prime},x}+f_{\rho^{\prime\prime}\oplus(\mu+x)}\mid
x\in A,\rho=\rho^{\prime}+\rho^{\prime\prime}\Big{\\}}$

  4. 4.

We can now conclude: given the definition of $\mathsf{ev}\circ\langle
f,g\rangle$, to check the Taylor equation it is enough to check that, for all
$N\in\mathbb{N}$,

$\left(((\cdots(\Lambda^{-}(f)\underbrace{\star g)\cdots)\star g}_{N\text{
times}})\circ\langle\mathrm{id},\infty\rangle\right)_{\chi,y}=\inf\Big{\\{}\sum_{i=1}^{N}g_{\chi_{i},x_{i}}+f_{\chi^{\prime},\langle[x_{1},\dots,x_{N}],y\rangle}\mid\begin{matrix}x_{1},\dots,x_{N}\in
A,\\\ \chi=\chi^{\prime}+\sum_{i=1}^{N}\chi_{i}\end{matrix}\Big{\\}}$

Let us show, by induction on $N$, the following equality, from which the
desired equality easily descends:

$\big{(}(\cdots(\Lambda^{-}(f)\underbrace{\star g)\cdots)\star g}_{N\text{
times}}\big{)}_{\chi\oplus\mu,y}=\inf\Big{\\{}\sum_{i=1}^{N}g_{\chi_{i},x_{i}}+f_{\chi^{\prime},\langle\mu+[x_{1},\dots,x_{N}],y\rangle}\mid\begin{matrix}x_{1},\dots,x_{N}\in
A,\\\ \chi=\chi^{\prime}+\sum_{i=1}^{N}\chi_{i}\end{matrix}\Big{\\}}$

     * •

if $N=0$, the right-hand term reduces to
$f_{\chi,\langle\mu,y\rangle}=(\Lambda^{-}(f))_{\chi\oplus\mu,y}$;

     * •

otherwise, let $F=((\cdots(\Lambda^{-}(f)\underbrace{\star g)\cdots)\star
g}_{N-1\text{ times}})$, so that by I.H. we have

$F_{\chi\oplus\mu,y}=\inf\Big{\\{}\sum_{i=1}^{N-1}g_{\chi_{i},x_{i}}+f_{\chi^{\prime},\langle\mu+[x_{1},\dots,x_{N-1}],y\rangle}\mid\begin{matrix}x_{1},\dots,x_{N-1}\in
A,\\\ \chi=\chi^{\prime}+\sum_{i=1}^{N-1}\chi_{i}\end{matrix}\Big{\\}}$

Then we have

$\displaystyle\big{(}F\star g\big{)}_{\chi\oplus\mu,y}$
$\displaystyle=\inf\left\\{g_{\chi^{\prime},x}+F_{\chi^{\prime\prime}\oplus(\mu+x)}\mid
x\in A,\chi=\chi^{\prime}+\chi^{\prime\prime}\right\\}$
$\displaystyle=\inf\left\\{g_{\chi^{\prime},x}+\inf\left\\{\sum_{i=1}^{N-1}g_{\chi_{i},x_{i}}+f_{\chi^{\prime\prime},\langle\mu+[x_{1},\dots,x_{N-1}]+x,y\rangle}\mid\begin{matrix}x_{1},\dots,x_{N-1}\in
A,\\\
\chi^{*}=\chi^{\prime\prime}+\sum_{i=1}^{N-1}\chi_{i}\end{matrix}\right\\}\
\Big{|}\ \begin{matrix}x\in A,\\\
\chi=\chi^{\prime}+\chi^{*}\end{matrix}\right\\}$
$\displaystyle=\inf\Big{\\{}g_{\chi^{\prime},x}+\sum_{i=1}^{N-1}g_{\chi_{i},x_{i}}+f_{\chi^{\prime\prime},\langle\mu+[x_{1},\dots,x_{N-1}]+x,y\rangle}\mid\begin{matrix}x,x_{1},\dots,x_{N-1}\in
A,\\\
\chi=\chi^{\prime}+\chi^{\prime\prime}+\sum_{i=1}^{N-1}\chi_{i}\end{matrix}\Big{\\}}$
$\displaystyle=\inf\Big{\\{}\sum_{i=1}^{N}g_{\chi_{i},x_{i}}+f_{\chi^{\prime},\langle\mu+[x_{1},\dots,x_{N}],y\rangle}\mid\begin{matrix}x_{1},\dots,x_{N}\in
A,\\\ \chi=\chi^{\prime}+\sum_{i=1}^{N}\chi_{i}\end{matrix}\Big{\\}}.$

##  Appendix E Proofs of Section 8: $\mathbb{L}$-modules and Generalized
Metric Spaces

###  E.1 Complete $\mathbb{L}$-categories and their $\mathbb{L}$-module
structure

In this subsection we quickly recall the notion of complete
$\mathbb{L}$-category and its associated $\mathbb{L}$-module structure.

Functors of shape $\Phi:X\times Y^{\mathrm{op}}\to\mathbb{L}$ are called
_distributors_ and usually noted
$\Phi:Y\mathrel{\ooalign{\hfil$\mapstochar\mkern 5.0mu$\hfil\cr$\to$\cr}}X$.

######  Definition 12 (weighted colimits).

Let $X,Y,Z$ be $\mathbb{L}$-categories,
$\Phi:Z\mathrel{\ooalign{\hfil$\mapstochar\mkern 5.0mu$\hfil\cr$\to$\cr}}Y$ be
a distributor and $f:Y\to X$ be a functor. A functor $g:Z\to X$ is the _$\Phi$
-weighted colimit of $f$ over $X$_, noted $\mathrm{colim}(\Phi,f)$, if for all
$z\in Z$ and $x\in X$

$\displaystyle X(g(z),x)=\sup_{y\in
Y}\left\\{X(f(y),x)\dotdiv\Phi(y,z)\right\\}$ (19)

A functor $f:X\to Y$ is _continuous_ if it commutes with all existing weighted
colimits in $X$, i.e. $f(\mathrm{colim}(\Phi,g))=\mathrm{colim}(\Phi,f\circ
g)$. A $\mathbb{L}$-enriched category $X$ is said _categorically complete_ (or
just _complete_) if all weighted colimits over $X$ exist.

An important example of colimit is the following:

######  Definition 13 (tensors).

Let $X$ be a $\mathbb{L}$-category, $x\in X$ and $\epsilon\in\mathbb{L}$. The
_tensor of $x$ and $\epsilon$_, if it exists, is the colimit $\epsilon\otimes
x:=\mathrm{colim}([\epsilon],\Delta x)$, where
$[\epsilon]:\\{\star\\}\mathrel{\ooalign{\hfil$\mapstochar\mkern
5.0mu$\hfil\cr$\to$\cr}}\\{\star\\}$ is the constantly $\epsilon$ distributor
and $\Delta x:\\{\star\\}\to X$ is the constant functor.

The $\mathbb{L}$-module structure of a complete $\mathbb{L}$-category has
order given by $x\preceq_{X}y$ iff $X(y,x)=0$, and action given by tensors
$\epsilon\otimes x$.

To conclude our correspondence between $\mathbb{L}$-modules and complete
$\mathbb{L}$-categories, it remains to observe that the two constructions
leading from one structure to the other are one the inverse of the other: for
any $\mathbb{L}$-module $(M,\preceq,\star)$, $x\preceq_{M}y$ iff $M(y,x)=0$
iff $x\preceq y=0\star y$, and, from $M(\epsilon\star
x,y)=M(x,y)\dotdiv\epsilon$, we deduce $\epsilon\otimes x=\epsilon\star x$.
Conversely, for any complete $\mathbb{L}$-category $X$ and $x,y\in X$, one can
check that $X(y,x)=\inf\\{\epsilon\mid X(\epsilon\otimes y,x)=0\\}$.

###  E.2 Exponential and Differential Structure of
$\mathbb{L}\mathsf{Mod}\simeq\mathbb{L}\mathsf{CCat}$

In this subsection we show that the category
$\mathbb{L}\mathsf{Mod}\simeq\mathbb{L}\mathsf{CCat}$ can be endowed with an
exponential modality $!$ so that that the coKleisli category
$\mathbb{L}\mathsf{Mod}_{!}$ is a model of the differential $\lambda$-calculus
extending the category $\mathbb{L}\mathsf{Rel}_{!}$.

First, we need to define a Lafont exponential $!$ over
$\mathbb{L}\mathsf{Mod}$. Since $\mathbb{L}\mathsf{Mod}$ is a SMCC with
biproducts, where the latter commute with tensors (see e.g. [70, Theorem
4.7.11]), we can apply a well-known recipe from [61, 51], which yields $!$ as
the _free exponential modality_ (i.e. such that $!X$ can be given the
structure of the cofree commutative comonoid over $X$).

First, we define the symmetric algebra $!_{n}M:=\mathrm{Sym}_{n}(M)$ as the
equalizer of all permutative actions on $n$-tensors $M\otimes\dots\otimes M$.
Notice that each element of $!_{n}M$ can be described as a join of “multisets”
$[x_{1},\dots,x_{n}]$, where the latter is the equivalence class of the tensor
$x_{1}\otimes\dots\otimes x_{n}\in M^{\otimes_{n}}$ under the permutative
actions. Moreover, the $\mathbb{L}$-module $!_{n}M$ is a complete
$\mathbb{L}$-category with distance function defined on basic “multisets” as
follows:

$\displaystyle(!_{n}M)(\alpha,\beta)=\sup_{\sigma\in\mathfrak{S}_{n}}\inf_{\tau\in\mathfrak{S}_{n}}\sum_{i=1}^{n}X(x_{\sigma(i)},y_{\tau(i)})$
(20)

where $\alpha=[x_{1},\dots,x_{n}]$ and $\beta=[y_{1},\dots,y_{n}]$, and
extended to arbitrary elements $\alpha=\bigvee_{i}\alpha_{i}$ and
$\beta=\bigvee_{j}\beta_{j}$ by
$(!_{n}M)(\alpha,\beta)=\sup_{i}\inf_{j}(!_{n}M)(\alpha_{i},\beta_{j})$.

Finally, we define $!M$ as the infinite biproduct $\prod_{n}!_{n}M$, yielding
the cofree commutative comonoid over $M$ (cf. [61, Proposition 1]).

The construction of $!$ for $\mathbb{L}\mathsf{Mod}$ generalizes the one for
$\mathbb{L}\mathsf{Rel}$:

######  Proposition 6.

$!(\mathbb{L}^{X})\simeq\mathbb{L}^{\mathcal{M}_{\mathrm{fin}}(X)}$. In
particular,
$\mathbb{L}\mathsf{Mod}_{!}(\mathbb{L}^{X},\mathbb{L}^{Y})\simeq\mathbb{L}\mathsf{Rel}_{!}(X,Y)$.

###### Proof.

Let us show that the morphism
$h:\mathbb{L}^{\mathcal{M}_{n}(S)}\to\mathbb{L}^{S\times\dots\times S}$
defined by $h(f)(\langle s_{1},\dots,s_{n}\rangle)=h([s_{1},\dots,s_{n}])$ is
the equalizer of the diagram
${\mathbb{L}^{\mathcal{M}_{n}(S)}}$${\mathbb{L}^{S\times\dots\times
S}}$${\mathbb{L}^{S\times\dots\times
S}}$$\scriptstyle{h}$$\scriptstyle{[\sigma]}$, where $[\sigma](x)(\langle
s_{1},\dots,s_{n}\rangle)=x(\langle
x_{\sigma(1)},\dots,x_{\sigma(n)}\rangle)$, with $\sigma$ varying over
$\mathfrak{S}_{n}$.

It is immediate that $h\circ[\sigma]=h\circ[\tau]$, for all
$\sigma,\tau\in\mathfrak{S}_{n}$. Let now
$k:C\to\mathbb{L}^{S\times\dots\times S}$ satisfy
$k\circ[\sigma]=k\circ[\tau]$: then for all $c\in C$, $k(c)(\langle
s_{1},\dots,s_{n}\rangle)=k(c)(\langle
s_{\sigma(1)},\dots,s_{\sigma(n)}\rangle)$, so $k(c)$ actually defines a
unique element of $\mathbb{L}^{\mathcal{M}_{n}(S)}$, and thus $k$ splits in a
unique way as $C\stackrel{{\scriptstyle
k^{\prime}}}{{\to}}\mathbb{L}^{\mathcal{M}_{n}(S)}\stackrel{{\scriptstyle
h}}{{\to}}\mathbb{L}^{S\times\dots\times S}$.

Now, observe that $\mathbb{L}^{S\times\dots\times
S}\simeq(\mathbb{L}^{S})^{\otimes_{n}}$ (cf. [70, Corollary 4.7.12 ($iii$)]),
and then, since equalizers are unique up to a unique isomorphism, we obtain an
isomorphism $\mathbb{L}^{\mathcal{M}_{n}(S)}\simeq!_{n}\mathbb{L}^{S}$. From
this we obtain the claim via
$!(\mathbb{L}^{X})=\prod_{n}!_{n}(\mathbb{L}^{X})\simeq\prod_{n}\mathbb{L}^{\mathcal{M}_{n}(X)}\simeq\mathbb{L}^{\coprod_{n}\mathcal{M}_{n}(X)}\simeq\mathbb{L}^{\mathcal{M}_{\mathrm{fin}}(X)}$.
∎

At this point, [54, Theorem 21], which states that an additive Lafont category
with free exponential modality and finite biproducts is a _differential
category_ [15], yields:

######  Theorem E.1.

$\mathbb{L}\mathsf{Mod}$ (equivalently $\mathbb{L}\mathsf{CCat}$) is a
differential category.

Finally, from Theorem E.1 we can conclude that $\mathbb{L}\mathsf{Mod}_{!}$
can be endowed with a differential operator $E$ making it a CC$\partial$C [18,
Proposition 3.2.1], i.e. Theorem 8.1 is proved.

To conclude, we make the definition of the differential operator $E$ of
$\mathbb{L}\mathsf{Mod}_{!}$ explicit: for $f:!M\to N$, we let

$\displaystyle Ef(\alpha)=\bigvee\left\\{f(\beta\cup[x])\ \Big{|}\
\iota_{n}(\beta)\otimes\iota_{1}(x)\leq S(\alpha)\right\\}$ (21)

where $\iota_{k}:M_{k}\to\prod_{i\in I}M_{i}$ is the injection morphism given
by $\iota_{k}(x)(k)=x$ and $\iota_{k}(x)(i\neq k)=\infty$, and $S:!(M\times
N)\to!M\otimes!N$ is the Seely isomorphism [61], and $E$ satisfies all
required axioms.

One can easily check that, when
$f\in\mathbb{L}\mathsf{Mod}_{!}(\mathbb{L}^{X},\mathbb{L}^{Y})\simeq\mathbb{L}\mathsf{Rel}_{!}(X,Y)$,
its derivative $Ef$ coincides with the derivative $D_{!}f$ defined for tps in
Section 5.


* [34] R. T. Rockafellar.  Conjugate Duality and Optimization.  SIAM, 1974.
  * [35] R. T. Rockafellar.  Extensions of subgradient calculus with applications to optimization.  Nonlinear Analysis: Theory, Methods & Applications, 9:665–698, 1985.
  * [36] R. T. Rockafellar and R. J-B Wets.  Variational Analysis.  Springer, 3rd printing-2009 edition, 1998.
  * [37] J. O. Royset.  Stability and error analysis for optimization and generalized equations.  SIAM J. Optimization, 30(1):752–780, 2020.
  * [38] J. O. Royset.  Good and bad optimization models: Insights from Rockafellians.  In J. G. Carlsson, editor, Tutorials in Operations Research: Emerging Optimization Methods and Modeling Techniques with Applications, pages 131–160. INFORMS, Cantonsville, MD, 2021.
  * [39] J. O. Royset.  Consistent approximations in composite optimization.  Mathematical Programming, 201:339–372, 2023.
  * [40] J. O. Royset, L. L. Chen, and E. Eckstrand.  Rockafellian relaxation and stochastic optimization under perturbations.  Preprint arXiv:2204.04762, 2022.
  * [41] J. O. Royset and R. J-B Wets.  Variational theory for optimization under stochastic ambiguity.  SIAM J. Optimization, 27(2):1118–1149, 2017.
  * [42] J. O. Royset and R. J-B Wets.  Lopsided convergence: an extension and its quantification.  Mathematical Programming, 177(1):395–423, 2019.
  * [43] J. O. Royset and R. J-B Wets.  An Optimization Primer.  Springer, 2021.
  * [44] M. Sion.  On general minimax theorems.  Pacific J. Mathematics, 8(1):171–176, 1958.