11institutetext: 1 International Centre for Theoretical Sciences, Tata
Institute of Fundamental Research, Bangalore 560089, India

# Generalized hydrodynamics and approach to Generalized Gibbs equilibrium for
a classical harmonic chain

Saurav Pandey1 Abhishek Dhar1 and Anupam Kundu1

###### Abstract

We study the evolution of a classical harmonic chain with nearest-neighbor
interactions starting from domain wall initial conditions. The initial state
is taken to be either a product of two Gibbs Ensembles (GEs) with unequal
temperatures on the two halves of the chain or a product of two Generalized
Gibbs Ensembles (GGEs) with different parameters in the two halves. For this
system, we construct the Wigner function and demonstrate that its evolution
defines the Generalized Hydrodynamics (GHD) describing the evolution of the
conserved quantities. We solve the GHD for both finite and infinite chains and
compute the evolution of conserved densities and currents. For a finite chain
with fixed boundaries, we show that these quantities relax as $\sim
1/\sqrt{t}$ to their respective steady-state values given by the final
expected GE or GGE state, depending on the initial conditions. Exact
expressions for the Lagrange multipliers of the final expected GGE state are
obtained in terms of the steady state densities. In the case of an infinite
chain, we find that the conserved densities and currents at any finite time
exhibit ballistic scaling while, at infinite time, any finite segment of the
system can be described by a current-carrying non-equilibrium steady state
(NESS). We compute the scaling functions analytically and show that the
relaxation to the NESS occurs as $\sim 1/t$ for the densities and as $\sim
1/t^{2}$ for the currents. We compare the analytic results from hydrodynamics
with those from exact microscopic numerics and find excellent agreement.

###### Contents

  1. 1 Introduction
  2. 2 Microscopic description of the system
  3. 3 Wigner function, conserved quantities and connection to GHD
    1. 3.1 Wigner function: its relation to the correlation matrix and evolution
    2. 3.2 Conserved quantities and connection to GHD
  4. 4 Comparison between microscopic and hydrodynamic evolution from domain wall initial condition
    1. 4.1 Domain wall initial condition composed of two GE states
    2. 4.2 Domain wall initial condition composed of two GGE states
  5. 5 Conclusions
  6. 6 Acknowledgements
  7. A GGE partition function
  8. B Density evolution corresponding to the initial condition composed of two GGEs

##  1 Introduction

Thermalization or the absence of it in isolated Hamiltonian systems has been
an active area of research, ever since the first numerical study of the FPUT
model fput1955 ; dauxois2008  was carried out. For generic non-integrable
systems, energy is typically the only conserved quantity. Such a system
usually thermalizes to a Gibbs Ensemble (GE) in the sense that, starting from
a non-equilibrium initial state, the macroscopic properties of the system
relax to the values determined by the GE. On the other hand, integrable
systems have an extensive number of conserved quantities and thus do not
thermalize to a GE. However, they are expected to relax to a steady state
described by a Generalized Gibbs Ensemble (GGE). The relaxation and approach
to a steady state has been shown in non-interacting integrable models such as
the hard particle gas (HPG) chakraborti2022  and the harmonic chain
vulpiani2022  and in interacting integrable models such as the alternate mass
hard particle gas (AMHPG) chakraborti2023  and the Toda chain ganapa2020 .
Thermalization for a low-density free fermion lattice model is shown
rigorously in tasaki2023  while relaxation of ideal quantum gases is studied
in pandey2023 . The approach to GGE however has been demonstrated mainly in
quantum integrable models rigol2007 ; rigol2011 ; calabrese2011 ; essler2022
and recently for the classical hard rods gas sahil2023 .

An effective understanding of a system’s macroscopic evolution and relaxation
to equilibrium can be achieved via the hydrodynamic approach. The hydrodynamic
theory is based on identifying the conserved quantities for writing down the
continuity equations involving densities and currents and the assumption of
local equilibrium for expressing currents in terms of the densities. For non-
integrable systems, one writes the hydrodynamic equations for the few
conserved quantities such as mass, momentum, and energy. In contrast,
integrable systems with an extensive number of conservation laws are described
by Generalized Hydrodynamics (GHD) which is written in terms of the quasi-
particle density. One needs to distinguish between non-interacting and
interacting integrable systems. For the former case, the hydrodynamic
equations are easier to write since the quasi-particles do not undergo phase
shifts upon collisions. Common examples include the harmonic chain for which
hydrodynamic equations were derived both at the Euler level dobrushin1986
(see also mielke2006 ) and with higher derivative corrections dobrushin1988 ,
and the HPG chakraborti2022 . On the other hand for the interacting case,
quasi-particles undergo phase shifts upon collision. The simplest example of
such systems is a collection of hard rods for which the hydrodynamic equations
were written long back percus1969exact ; boldrighini1983euler ;
boldrighini1997navierstokes . For the case of interacting integrable quantum
systems, GHD was derived more recently castro2016emergent ;
bertini2016transport  and quite remarkably they have the same structure as the
classical systems doyon2020 . A number of studies have established GHD for
several interacting classical and quantum integrable systems such as Toda
chain spohn2021toda ; cao2019gge ; mazzuca2023equilibrium ;
kundu2023integrable , the $\delta$-Bose gas ruggiero2020  and hard rods
doyon2017 ; sahil2023 . GHD has been successfully used to understand unusual
equilibration of trapped integrable systems caux2019hydrodynamics ;
schemmer2019generalized ; cao2018incomplete  and also the predictions of GHD
have been experimentally verified malvania2021generalized .

The harmonic chain provides us with an example of an integrable system that is
analytically tractable and allows macroscopic description in terms of
hydrodynamics dobrushin1986 ; dobrushin1988 . For this system, a detailed
comparison of the microscopic and hydrodynamic calculations is possible.
Approach to GGE starting from a domain wall initial condition with unequal
temperatures in the two halves has been studied for the case of an infinite
quantum harmonic chain in eisler2014  where the Lagrange multipliers of the
final GGE state were evaluated explicitly. Convergence of the phase space
distribution to a Gaussian stationary measure, starting from an initial state
which is also Gaussian, is studied in dudnikova2003convergence . In
boldrighini1983conservedQ , two families of conserved quantities for infinite
harmonic lattices were shown and the convergence of the covariance matrix to
its stationary form was studied. Hydrodynamic approaches have also been
applied to harmonic systems. In dobrushin1986 , starting from a family of
initial states that are slowly varying in space, the hydrodynamic equation of
an infinite classical harmonic chain was derived at the Euler level in terms
of what is referred to as the spectral density matrix function (SDMF) in the
paper, which is essentially the Fourier transform of the correlation matrix
defined in a coarse-grained cell around a macroscopic point. The study was
extended in dobrushin1988  where the next-order correction (in space-time
scaling parameter) to the Euler equation was derived. Alternatively,
macroscopic evolution in the harmonic chain has been studied using Wigner
functions spohn2006boltzmann ; dudnikova2005local ; bernardin2019 . The
harmonic chain being an integrable system, one expects a generalized
hydrodynamic description for the macro-evolution. However, the connections
between the correlation matrix approach dobrushin1986 ; dobrushin1988  and the
Wigner function method spohn2006boltzmann  have not been elaborated in the
physics literature. Furthermore, and somewhat surprisingly, their relationship
to GHD remains unexplored. In this paper, we aim to address this gap and
elucidate these connections which we use to understand non-equilibrium
evolution and approach to the stationary state.

In this paper, we consider a classical harmonic chain with nearest-neighbor
interactions either of finite length (with fixed boundary conditions) or of
infinite extent. Initially, we prepare the system in a domain wall
configuration by taking the two halves of the chain either in thermal
equilibrium with unequal temperatures or in different GGEs with unequal
parameters (Lagrange multipliers). Starting from this initial condition, we
study the evolution of the chain and the approach to the stationary state in
the long time limit (for the finite case). Following is a summary of our main
results:

  * $\bullet$

We clarify how the Wigner function method spohn2006boltzmann  is related to
the correlation matrix approach in dobrushin1986 . In particular, we write an
explicit expression for the Wigner function in terms of the Fourier transform
of the local correlation function — this allows us to show the equivalence
between the Euler equation in dobrushin1986  and the equation satisfied by the
Wigner function [see Sec. (3.1)].

  * $\bullet$

We express all the local conserved densities and currents in terms of the
Wigner function which allows us to write the hydrodynamic equations for all
the conserved quantities. The equation for the coarse-grained Wigner function
can thus be considered as the GHD equation for the harmonic chain. This
approach is similar to that used for deriving the GHD equation for free
fermions essler2022 .

  * $\bullet$

Using the non-interacting structure, we solve the GHD equations analytically
and obtain explicit expressions for the time evolution of the density and
current profiles starting from non-equilibrium initial conditions. We compare
these solutions with microscopic computations and find remarkable agreement.

  * $\bullet$

On the infinite line, the domain wall spreads ballistically and from our
solution of GHD equations we can compute the exact scaling functions
corresponding to both the density and current profiles. At any fixed location,
the density and the current corresponding to any conserved quantity approach a
stationary finite value in the $t\to\infty$ limit. The approach to the steady
state value is $\sim 1/t$ and $\sim 1/t^{2}$, for densities and currents
respectively. Therefore the stationary state inside any finite segment should
be describable by a GGE which, in addition to the conserved densities, also
includes all the associated currents that are conserved in the infinite
geometry eisler2014 .

  * $\bullet$

For the domain wall initial condition on a finite chain, we observe that the
chain relaxes to a final stationary state that is consistent with a GGE state
in which all the conserved densities become homogeneous in space and
stationary in time while the corresponding currents vanish. Note that in this
case, current conservation is violated by the reflections at the boundaries.
The conserved densities and the associated currents approach their stationary
values as $1/\sqrt{t}$ in the $t\to\infty$ limit.

  * $\bullet$

The final expected GGE state in the finite chain is characterized by chemical
potentials $\\{\lambda_{n}\\}$ that are completely determined by the values of
the conserved quantities. We provide explicit expressions of these chemical
potentials in terms of the conserved densities. For the particular case of the
domain wall initial condition composed of two Gibbs states at different
temperatures, the final stationary state is also a Gibbs state. For the other
case studied in this work in which the initial domain wall state is composed
of two GGEs where only the first two Lagrange parameters are non-zero in
either half of the chain, we find that the final state is also expected to be
a GGE with Lagrange parameters that decay exponentially in strength for large
$n$.

The rest of the paper is organized as follows. In Sec. (2), we introduce our
system of a classical harmonic chain with nearest-neighbor interactions and
study its microscopic evolution averaged over initial states drawn from a
distribution that is a product of two Generalized Gibbs Ensembles with
different values for the chemical potentials in the two halves of the chain.
The Wigner function for harmonic lattices is introduced in Sec. (3.1) and its
connection to the correlation matrix approach used in the literature is made
explicit. In Sec. (3.2), we construct the conserved quantities for the
harmonic chain and express the corresponding densities and currents in terms
of the Wigner function. In Sec. (4), we derive the hydrodynamic evolution for
the conserved quantities using the Wigner function approach and compare them
with the exact microscopic results. The case where the two halves are
initially described by a GE with different temperatures is studied in Sec.
(4.1) whereas Sec. (4.2) deals with the case where the two halves are
initially described by a GGE with only the first two parameters being nonzero
and unequal in the left and right parts of the chain. We conclude our findings
in Sec. (5).

##  2 Microscopic description of the system

…..…..$\\{\lambda_{n}^{L}\\}$$\\{\lambda_{n}^{R}\\}$$t=0$$1$$2$N$-$1NN$+$1N$+$22N$-$12N

…..…..$t>0$

Figure 1: Plot showing a schematic of the harmonic chain system of length
$2N$. Initially, at $t=0$, the left and right halves are disconnected and in
general, described by two GGEs with unequal chemical potentials
$\lambda_{n}^{L}\neq\lambda_{n}^{R}$. At $t>0$, the two halves interact and
the system evolves as a whole towards a new stationary state.

In this section, we define the microscopic model of our classical harmonic
chain and state the exact expressions for various correlation functions that
can be obtained from the solution of the microscopic dynamics. We consider a
chain of $2N$ particles of unit masses whose positions and momenta are denoted
by $\\{q_{j},p_{j}\\},~{}j=1,\ldots 2N$. Assuming fixed boundary conditions:
$q_{0}=q_{2N+1}=0$, the Hamiltonian is given by

$H=\sum\limits_{j=1}^{2N}\frac{p_{j}^{2}}{2}+\sum\limits_{j=0}^{2N}\frac{(q_{j+1}-q_{j})^{2}}{2}=\displaystyle{\frac{1}{2}}\left(\mathbf{p}^{T}\mathbf{p}+\mathbf{q}^{T}\Phi\mathbf{q}\right),$
(1)

where $\mathbf{q}=\\{q_{1},\ldots,q_{2N}\\}$ and
$\mathbf{p}=\\{p_{1},\ldots,p_{2N}\\}$ are vectors containing the particle
positions and momenta respectively and $\Phi$ denotes the force matrix with
elements
$\Phi_{j,j^{\prime}}=2\delta_{j,j^{\prime}}-\delta_{j+1,j^{\prime}}-\delta_{j-1,j^{\prime}}$
for $j,j^{\prime}=1,2,...,2N$, where $\delta_{jj^{\prime}}$ is Kronecker
delta. This gives us the equations of motion (EOM)

$\displaystyle\ddot{\mathbf{q}}$ $\displaystyle=-\Phi\mathbf{q}.$ (2)

Let $U$ be the matrix of eigenfunctions that diagonalises $\Phi$ i.e.,
$U^{T}\Phi U=\Omega^{2}$ where
$\Omega^{2}_{k,k^{\prime}}=\omega_{k}^{2}\delta_{k,k^{\prime}}$. The
eigenvalues and eigenfunctions are given by

$\displaystyle\omega_{k_{\ell}}^{2}=2(1-\cos k_{\ell}),$ (3a)
$\displaystyle\psi_{k_{\ell}}(j)=U_{jk_{\ell}}=\displaystyle{\sqrt{\frac{2}{2N+1}}}\sin
jk_{\ell},$ (3b)
$\displaystyle\text{where}~{}k_{\ell}=\displaystyle{\frac{\pi\ell}{2N+1}},~{}\ell=1,\ldots,2N.$
(3c)

We can write down a general solution to Eq. (2) in matrix notation as

$\displaystyle\mathbf{q}(t)$
$\displaystyle=\dot{A}(t)\mathbf{q}(0)+A(t)\mathbf{p}(0),$ (4a)
$\displaystyle\mathbf{p}(t)$
$\displaystyle=\ddot{A}(t)\mathbf{q}(0)+\dot{A}(t)\mathbf{p}(0),$ (4b)

where we have defined

$\displaystyle A(t)=U\displaystyle{\frac{\sin(\Omega t)}{\Omega}}U^{T}.$ (5)

For a harmonic chain of $N$ particles, one can construct $N$ number of local
conserved quantities which we denote by $Q^{(n)}({\bf
p,q}),~{}n=0,1,2,...,N-1$. These conserved quantities are of the form:

$\displaystyle Q^{(n)}({\bf
p,q})=\frac{1}{2}\left(\mathbf{p}^{T}B^{(n)}\mathbf{p}+\mathbf{q}^{T}M^{(n)}\mathbf{q}\right),$
(6)

where $B^{(n)}$ and $M^{(n)}$ are $N\times N$ symmetric matrices [see Sec.
(3.2) for a detailed discussion]. The general stationary state of the chain is
therefore described by a GGE of the form

$\displaystyle P_{\rm
GGE}({\bf{p,q}})=\frac{1}{Z(\\{\lambda_{n}\\})}\exp\left(-\sum\limits_{n=0}^{N-1}\lambda_{n}Q^{(n)}\right),$
(7)

where $\lambda_{n}$ are the corresponding Lagrange multipliers that determine
the average values of the conserved quantities. Here, $Z(\\{\lambda_{n}\\})$
is the GGE partition function for a chain of $N$ particles, which is
explicitly given by [see Appendix (A)]

$\displaystyle
Z(\\{\lambda_{n}\\})=(2\pi)^{N}\prod_{\ell=1}^{N}\frac{1}{\omega_{k_{\ell}}\sum_{n=0}^{N-1}\lambda_{n}\cos\left(nk_{\ell}\right)},$
(8)

where $\omega_{k}$ and $k_{\ell}$ are given in Eqs. (3a, 3c), respectively.

Initial condition: As shown in the schematic in Fig. (1), our system comprises
of two such $N$-particle chains that are initially in a domain wall
configuration described by a product of two GGE states

$\displaystyle P_{\rm
in}=\frac{e^{-\sum_{n}\lambda_{n}^{L}Q^{(n)}_{L}}}{Z^{L}}\times\frac{e^{-\sum_{n}\lambda_{n}^{R}Q^{(n)}_{R}}}{Z^{R}},$
(9)

with different values of the Lagrange multipliers in the left and right
halves, denoted respectively by $\\{\lambda^{L}_{n}\\}$ and
$\\{\lambda^{R}_{n}\\}$, corresponding to the conserved quantities
$Q^{(n)}_{L}$ and $Q^{(n)}_{R}$. Substituting the expression for $Q^{(n)}$
from Eq. (6) into Eq. (9), we can express $P_{\rm in}$ in the following
Gaussian form

$\displaystyle P_{\rm
in}=\frac{1}{Z^{L}Z^{R}}\exp\left(-\displaystyle{\frac{1}{2}}[\mathbf{p}^{\text{T}}\mathcal{B}\mathbf{p}+\mathbf{q}^{\text{T}}\mathcal{M}\mathbf{q}]\right),$
(10)

where $\mathcal{B}$ and $\mathcal{M}$ are $2N\times 2N$ symmetric matrices
defined as

$\displaystyle\mathcal{B}=\begin{pmatrix}\sum\limits_{n}\lambda_{n}^{L}B^{(n)}&0\\\
0&\sum\limits_{n}\lambda_{n}^{R}B^{(n)}\end{pmatrix},~{}~{}~{}~{}\mathcal{M}=\begin{pmatrix}\sum\limits_{n}\lambda_{n}^{L}M^{(n)}&0\\\
0&\sum\limits_{n}\lambda_{n}^{R}M^{(n)}\end{pmatrix},$ (11)

where the Lagrange multipliers are chosen such that $\mathcal{B}$ and
$\mathcal{M}$ are positive definite matrices.

The initial correlations can be immediately inferred from Eq. (10) and are
given by $\expectationvalue{\mathbf{q}(0)\mathbf{q}(0)^{T}}=\mathcal{M}^{-1}$,
$\expectationvalue{\mathbf{p}(0)\mathbf{p}(0)^{T}}=\mathcal{B}^{-1}$, and
$\expectationvalue{\mathbf{p}(0)\mathbf{q}(0)^{T}}=0$. Using Eq. (4), the
equal time two-point correlations can then be easily obtained and we find

$\displaystyle\expectationvalue{\mathbf{q}(t)\mathbf{q}(t)^{T}}$
$\displaystyle=A(t)\mathcal{B}^{-1}A(t)+\dot{A}(t)\mathcal{M}^{-1}\dot{A}(t),$
(12a) $\displaystyle\expectationvalue{\mathbf{p}(t)\mathbf{p}(t)^{T}}$
$\displaystyle=\dot{A}(t)\mathcal{B}^{-1}\dot{A}(t)+\ddot{A}(t)\mathcal{M}^{-1}\ddot{A}(t),$
(12b) $\displaystyle\expectationvalue{\mathbf{q}(t)\mathbf{p}(t)^{T}}$
$\displaystyle=A(t)\mathcal{B}^{-1}\dot{A}(t)+\dot{A}(t)\mathcal{M}^{-1}\ddot{A}(t),$
(12c)

where $A(t)$ is given in Eq. (5).

##  3 Wigner function, conserved quantities and connection to GHD

In this section, we discuss the Generalized Hydrodynamics (GHD) for the
harmonic chain through its connection to the Wigner function formalism. This
is similar to the approach followed in essler2022  for non-interacting
fermions on a one-dimensional lattice. Using GHD, we obtain an analytical
understanding of the correlations and conserved quantities on a macroscopic
scale. Unless stated otherwise, we will assume the chain to be infinite at
both ends for which the eigen-spectrum in Eq. (3a) becomes continuous:
$\omega^{2}(k)=2(1-\cos k)$ where $k\in[-\pi,\pi]$.

### 3.1 Wigner function: its relation to the correlation matrix and evolution

Although the Wigner function was originally introduced as a quantum analogue
of the classical phase space distribution wigner1932 ; hillery1997 , it has
also been used widely in the context of transport and hydrodynamics of
classical harmonic crystals. Following spohn2006boltzmann , we construct the
Wigner function for the classical harmonic chain defined in Eq. (1). For this,
we first define the quantity

$\displaystyle
a(k,\tau)=\displaystyle{\frac{1}{\sqrt{2}}}\left[\sqrt{\omega(k)}q(k,\tau)+\frac{i}{\sqrt{\omega(k)}}p(k,\tau)\right],$
(13)

where the Fourier and inverse Fourier transforms are defined as

$\displaystyle q(k,\tau)$
$\displaystyle=\sum_{j\in\mathbb{Z}}e^{-ijk}q_{j}(\tau);~{}~{}~{}q_{j}(\tau)=\int_{-\pi}^{\pi}\displaystyle{\frac{dk}{2\pi}}e^{ijk}q(k,\tau),$
(14a) $\displaystyle p(k,\tau)$
$\displaystyle=\sum_{j\in\mathbb{Z}}e^{-ijk}p_{j}(\tau);~{}~{}~{}p_{j}(\tau)=\int_{-\pi}^{\pi}\displaystyle{\frac{dk}{2\pi}}e^{ijk}p(k,\tau),$
(14b) $\displaystyle a(k,\tau)$
$\displaystyle=\sum_{j\in\mathbb{Z}}e^{-ijk}a_{j}(\tau);~{}~{}~{}a_{j}(\tau)=\int_{-\pi}^{\pi}\displaystyle{\frac{dk}{2\pi}}e^{ijk}a(k,\tau),$
(14c)

where $a_{j}(\tau)$ is defined as the Fouier transform of $a(k,\tau)$. The
quantity $a(k,\tau)$ has a simple evolution:
$a(k,\tau)=e^{-i\omega(k)\tau}a(k,0)$. We then define the Wigner function in
terms of $a(k,\tau)$ as

$\displaystyle W^{1}(y,k,\tau)$
$\displaystyle=\int_{-2\pi}^{2\pi}\displaystyle{\frac{d\xi}{4\pi}}e^{i\xi
y}\expectationvalue{a^{\star}(k-\xi/2,\tau)a(k+\xi/2,\tau)}=\sum_{j\in\mathbb{Z}}e^{-ijk}\expectationvalue{a^{\star}_{y-j/2}(\tau)a_{y+j/2}(\tau)},$
(15)

where $y\in\mathbb{Z}/2$. The sum on $j$ runs over even/odd integers if $y$ is
integer/half-integer and the average is over an initial ensemble. Substituting
Eq. (13) into Eq. (15), we obtain

$\displaystyle\begin{split}W^{1}(y,k,\tau)&=\displaystyle{\frac{1}{2}}\int_{-2\pi}^{2\pi}\displaystyle{\frac{d\xi}{4\pi}}\sum_{j,j^{\prime}\in\mathbb{Z}}e^{i\xi(y-(j+j^{\prime})/2)}e^{i(j-j^{\prime})k}\\\
&\times~{}\left[\sqrt{\omega(k-\xi/2)\omega(k+\xi/2)}\expectationvalue{q_{j}q_{j^{\prime}}}(\tau)+i\sqrt{\frac{\omega(k-\xi/2)}{\omega(k+\xi/2)}}\expectationvalue{q_{j}p_{j^{\prime}}}(\tau)\right.\\\
&~{}~{}\left.-i\sqrt{\frac{\omega(k+\xi/2)}{\omega(k-\xi/2)}}\expectationvalue{p_{j}q_{j^{\prime}}}(\tau)+\frac{1}{\sqrt{\omega(k-\xi/2)\omega(k+\xi/2)}}\expectationvalue{p_{j}p_{j^{\prime}}}(\tau)\right].\end{split}$
(16)

Note that the above expression readily provides connection between the Wigner
function and the correlations. If we assume the correlations to be homogeneous
in space, then we obtain a steady-state Wigner function

$\displaystyle
W_{ss}(k)=\displaystyle{\frac{1}{2}}\left[\omega(k)F_{11}(k)-iF_{12}(k)+iF_{21}(k)+\displaystyle{\frac{1}{\omega(k)}}F_{22}(k)\right],$
(17)

where $F_{mn}(k)$ are elements of the $2\times 2$ correlation matrix

$\displaystyle
F(k)=\sum_{r\in\mathbb{Z}}e^{ikr}\bar{F}(r),~{}~{}\text{where}~{}~{}\bar{F}({j-j^{\prime}})=\begin{pmatrix}\expectationvalue{q_{j}q_{j^{\prime}}}&\expectationvalue{q_{j}p_{j^{\prime}}}\\\\[5.0pt]
\expectationvalue{p_{j}q_{j^{\prime}}}&\expectationvalue{p_{j}p_{j^{\prime}}}\end{pmatrix}.$
(18)

We point out that $F_{11}^{*}(k)=F_{11}(k),~{}F_{22}^{*}(k)=F_{22}(k)$ and
$F_{12}^{*}(k)=F_{21}(k)$ so that the Wigner function is always real. In
thermal equilibrium at an inverse temperature $\beta$, we simply get

$\displaystyle W_{\beta}(k)=\displaystyle{\frac{1}{\beta\omega(k)}}.$ (19)

In general, for a chain of $N$ particles in GGE described by the distribution
in Eq. (7), the Wigner function can be obtained by computing the correlations
$F_{j-j^{\prime}}$ in Eq. (18) and using them in Eq. (17). One finds the
following explicit expression

$\displaystyle
W_{GGE}(k)=\displaystyle{\frac{1}{\omega(k)\sum\limits_{n=0}^{N-1}\lambda_{n}\cos
nk}}.$ (20)

On the other hand, when the correlation functions are inhomogeneous and evolve
over a large space-time scale, one can assume the system to be locally in a
steady state characterized by the local values of the correlations.
Transforming to the scaled variables: $x=\epsilon y,t=\epsilon\tau$, and
taking the $\epsilon\to 0$ limit in Eq. (16), the relation in Eq. (17) between
the Wigner function and the correlation matrix gets modified to

$\displaystyle
W(x,k,t)=\displaystyle{\frac{1}{2}}\left[\omega(k)F_{11}(x,k,t)-iF_{12}(x,k,t)+iF_{21}(x,k,t)+\displaystyle{\frac{1}{\omega(k)}}F_{22}(x,k,t)\right],$
(21)

where

$\displaystyle W(x,k,t)$ $\displaystyle:=\lim\limits_{\epsilon\to
0}W^{1}(x/\epsilon,k,t/\epsilon)=\lim\limits_{\epsilon\to
0}\epsilon\int_{-2\pi/\epsilon}^{2\pi/\epsilon}\displaystyle{\frac{d\eta}{4\pi}}e^{i\eta
x}\expectationvalue{a^{\star}(k-\epsilon\eta/2,t/\epsilon)a(k+\epsilon\eta/2,t/\epsilon)},$
(22)

and

$\displaystyle F_{mn}(x,k,t)=\lim\limits_{\epsilon\to
0}\sum\limits_{r=-\infty}^{\infty}e^{-irk}\expectationvalue{y^{m}_{[\epsilon^{-1}x-r/2]}~{}y^{n}_{[\epsilon^{-1}x+r/2]}}(\epsilon^{-1}t).$
(23)

Here $y^{1}_{j}=q_{j}$ is the displacement and $y^{2}_{j}=p_{j}$ is the
momentum of the $j^{\text{th}}$ oscillator. The local correlation functions
$F_{mn}(x,k,t)$ were shown, in Ref. dobrushin1986 , to satisfy the Euler
equations, using which one can show that the Wigner function $W(x,k,t)$
satisfies the following transport equation:

$\displaystyle\partial_{t}W(x,k,t)+\omega^{\prime}(k)\partial_{x}W(x,k,t)=0.$
(24)

A more direct way to derive this equation is to start with $W^{1}(y,k,\tau)$
defined in Eq. (15) which satisfies

$\displaystyle\partial_{\tau}W^{1}(y,k,\tau)=i\int_{-2\pi}^{2\pi}\displaystyle{\frac{d\xi}{4\pi}}e^{i\xi
y}[\omega(k-\xi/2)-\omega(k+\xi/2)]\expectationvalue{a^{\star}(k-\xi/2,\tau)a(k+\xi/2,\tau)}.$
(25)

We then transform to the scaled variables $(x=\epsilon y,t=\epsilon\tau)$ on
both sides of Eq. (25) and expand the right-hand side in $\epsilon$. At the
leading order in $\epsilon$, one finds the evolution equation (24) for
$W(x,k,t)$ defined in Eq. (22).

On an infinite line, the solution of Eq. (24) at a later time can be simply
obtained by boosting the initial function $W(x,k,0)$ with the phonon group
velocity $\omega^{\prime}(k)$:

$\displaystyle W(x,k,t)=W(x-\omega^{\prime}(k)t,k,0).$ (26)

Note that in thermal equilibrium the average internal energy at a site can be
obtained from $W_{\beta}(k)$ [given in Eq. (19)] as
$\int_{-\pi}^{\pi}dk/(2\pi)\omega(k)W_{eq}(k)=1/\beta$. Generalizing this
relation to the inhomogeneous case we find that the average internal energy
density can be expressed in terms of the Wigner function as

$\displaystyle
e(x,t)=\int_{-\pi}^{\pi}\displaystyle{\frac{dk}{2\pi}}\omega(k)W(x,k,t).$ (27)

It is easy to see that $e(x,t)$ satisfies a continuity equation
$\partial_{t}e(x,t)+\partial_{x}\mathfrak{j}(x,t)=0$ where the average energy
current density is given by

$\displaystyle\mathfrak{j}(x,t)=\int_{-\pi}^{\pi}\displaystyle{\frac{dk}{2\pi}}\omega(k)\omega^{\prime}(k)W(x,k,t).$
(28)

In the following, we discuss other conserved densities and the associated
currents and show how one can express them in terms of the Wigner function.

### 3.2 Conserved quantities and connection to GHD

The Hamiltonian in Eq. (1) for $N\to\infty$ can be written in Fourier space as

$\displaystyle
H=\bigintsss\limits_{-\pi}^{\pi}\displaystyle{\frac{dk}{2\pi}}\mathcal{E}(k)~{}~{}\text{with}~{}~{}\mathcal{E}(k)=\displaystyle{\frac{\omega(k)}{2}}[a^{\star}(k,\tau)a(k,\tau)+a^{\star}(-k,\tau)a(-k,\tau)],$
(29)

where $a(k,\tau)$ is defined in Eq. (13). Note that
$\mathcal{E}(k)={\frac{|p(k,\tau)|^{2}}{2}}+\omega(k)^{2}\frac{|q(k,\tau)|^{2}}{2}$
represents the energy of a single mode and is a conserved quantity. This
allows us to define a set of new conserved quantities as the following linear
combinations:

$\displaystyle
Q^{(n)}:=\int\limits_{-\pi}^{\pi}\displaystyle{\frac{dk}{2\pi}}\mathcal{E}(k)e^{ink}=\int\limits_{-\pi}^{\pi}\displaystyle{\frac{dk}{2\pi}}\cos
nk~{}\omega(k)a^{\star}(k,\tau)a(k,\tau)~{}~{}\text{for}~{}n=0,1,2,...$ (30)

We now show that these conserved quantities are local in nature, i.e., they
can be expressed as a sum over local terms. To see this, we insert Eq. (13)
into Eq. (30) and after some simplifications, obtain

$\displaystyle
Q^{(n)}=\sum_{\ell\in\mathbb{Z}}\varrho^{(n)}_{\ell}(\tau)~{}~{}\text{where}~{}~{}\varrho^{(n)}_{\ell}(\tau)=\displaystyle{\frac{1}{2}}\left[p_{\ell}p_{\ell-n}+2q_{\ell}q_{\ell-n}-q_{\ell}q_{\ell-n+1}-q_{\ell}q_{\ell-n-1}\right].$
(31)

It is easy to see that the microscopic densities $\varrho^{(n)}_{\ell}(\tau)$
obey a discrete continuity equation

$\displaystyle\partial_{\tau}\varrho^{(n)}_{\ell}(\tau)=\mathfrak{j}^{(n)}_{\ell-1}(\tau)-\mathfrak{j}^{(n)}_{\ell}(\tau),$
(32)

where the microscopic current $\mathfrak{j}^{(n)}_{\ell}(\tau)$ from site
$\ell$ to $\ell+1$ is given by

$\displaystyle\mathfrak{j}^{(n)}_{\ell}(\tau)=\displaystyle{\frac{1}{2}}\left(p_{\ell-n+1}q_{\ell}-p_{\ell-n}q_{\ell+1}\right).$
(33)

Note that $Q^{(n)}$ in Eq. (31) can also be written in matrix form as

$\displaystyle
Q^{(n)}=\displaystyle{\frac{1}{2}}\left(\mathbf{p}^{T}B^{(n)}\mathbf{p}+\mathbf{q}^{T}M^{(n)}\mathbf{q}\right),$
(34)

where $B^{(n)}$ and $M^{(n)}$ are symmetric matrices defined as

$\displaystyle B^{(n)}_{\ell,\ell^{\prime}}$
$\displaystyle=\frac{1}{2}\left(\delta_{\ell-n,\ell^{\prime}}+\delta_{\ell+n,\ell^{\prime}}\right),$
(35a) $\displaystyle M^{(n)}_{\ell,\ell^{\prime}}$
$\displaystyle=\delta_{\ell-n,\ell^{\prime}}+\delta_{\ell+n,\ell^{\prime}}-\frac{1}{2}\left(\delta_{\ell-n+1,\ell^{\prime}}+\delta_{\ell+n-1,\ell^{\prime}}+\delta_{\ell-n-1,\ell^{\prime}}+\delta_{\ell+n+1,\ell^{\prime}}\right).$
(35b)

These matrix forms are particularly useful for numerical microscopic
computations. For initial conditions chosen from a distribution such as in Eq.
(9), one is interested in the average values of these conserved densities
$\rho_{\ell}^{(n)}:=\expectationvalue{\varrho_{\ell}^{(n)}}$, which, on a
macroscopic scale can be expressed in terms of the coarse-grained Wigner
function $W(x,k,t)$. Noting
$\expectationvalue{a^{*}(k,\tau)a(k,\tau)}=\sum_{y}W^{1}(y,k,\tau)$ from Eq.
(15) and substituting it in Eq. (30) gives

$\displaystyle\expectationvalue{Q^{(n)}}=\sum_{y}\int\limits_{-\pi}^{\pi}\displaystyle{\frac{dk}{2\pi}}\cos
nk~{}\omega(k)W^{1}(y,k,\tau).$ (36)

Defining $\rho^{(n)}(x,t)=\lim_{\epsilon\to
0}\epsilon^{-1}\rho^{(n)}_{[\epsilon^{-1}x]}(\epsilon^{-1}t)$, we can write
$\expectationvalue{Q^{(n)}}=\int dx\rho^{(n)}(x,t)$, where

$\displaystyle\rho^{(n)}(x,t)=\int\limits_{-\pi}^{\pi}\displaystyle{\frac{dk}{2\pi}}\cos
nk~{}\omega(k)W(x,k,t).$ (37)

Using Eq. (24), it is easy to see that the densities in Eq. (37) obey the
continuity equation

$\displaystyle\partial_{t}\rho^{(n)}(x,t)+\partial_{x}\mathfrak{j}^{(n)}(x,t)=0,$
(38)

where the macroscopic currents $\mathfrak{j}^{(n)}(x,t)$ are given by

$\displaystyle\mathfrak{j}^{(n)}(x,t)=\int\limits_{-\pi}^{\pi}\displaystyle{\frac{dk}{2\pi}}\cos
nk~{}\omega(k)\omega^{\prime}(k)W(x,k,t).$ (39)

We now comment on the connection to GHD. We see that the coarse-grained Wigner
function in Eq. (22) indeed has the structure of GHD equations if we identify
$W(x,k,t)$ as a phase space distribution of non-interacting quasi-particles
(basically the phonons) with position $x$, momentum $k$, and velocity
$\omega^{\prime}(k)$. This identification is natural as was pointed out in
spohn2006boltzmann . The definition of conserved densities and currents in
Eqs. (37, 39) is consistent with this identification. One question is
regarding the positivity of $W$. At the microscopic level, it can be negative
but coarse-graining can lead to a positive-definite quantity. Various
possibilities for coarse-graining exist, for example through the construction
of the Husimi function (see mielke2006 ). Another more physical construction
was discussed in pandey2023  in the context of a non-interacting quantum gas
where the coarse-grained Wigner function has the direct interpretation as the
number of quasi-particles in a coarse-grained cell. Note that at a physical
level, the GHD equation for the harmonic chain is the same as the Peierls-
Boltzmann equation for a weakly anharmonic chain with collision terms
neglected (see simoncelli2022 ).

##  4 Comparison between microscopic and hydrodynamic evolution from domain
wall initial condition

In this section, we apply the hydrodynamic formulation developed in Sec. (3)
to the study of equilibration of the harmonic chain starting from a domain
wall initial condition as depicted in Fig. (1). We study both a chain of
infinite extent and a chain of finite length with fixed boundaries. For the
hydrodynamics of a finite-sized chain, we assume its length to be $2L$ with
$2N$ particles and lattice constant $\epsilon$ such that $N\to\infty$ and
$\epsilon\to 0$ with $L=N\epsilon$ held fixed. The microscopic numerical
calculations are done for finite $N$ with unit lattice constant and to compare
the results, we make the identification $L=N$. For both cases, we consider two
types of domain wall initial conditions. In the first choice, only the
(inverse) temperatures in the two halves of the chain are different, namely
$\beta_{1}$ and $\beta_{2}$, while all other Lagrange multipliers are zero
i.e., $\lambda^{L}_{n}=\lambda^{R}_{n}=0,~{}\forall~{}n\geq 1$. For a finite
chain, this implies that the only conserved quantity with a non-zero average
value initially is the energy. As a consequence of conservation laws, the
values of higher conserved quantities remain zero in the final stationary
state which is thus again expected to be described by a GE with the final
temperature being the mean of the initial temperatures in the two halves. In
Sec. (4.1), by studying the evolution of the conserved densities we show how
the system relaxes to the final expected GE state. In Sec. (4.2), we consider
the second choice of the initial condition in which the first two parameters
(Lagrange multipliers) are non-zero and take different values in the two
halves while other parameters are zero on both sides. In this case, all the
conserved quantities have non-zero average values initially, and consequently,
we expect the system to finally reach a GGE stationary state. Once again, by
solving the GHD equation in Eq. (24), we study the evolution of the conserved
densities and the associated currents and show how the system in this case
approaches the final expected GGE state. Note that for both types of initial
conditions, the infinite chain is expected to go to a GGE stationary state
which is current carrying. This is due to the presence of additional conserved
quantities (see boldrighini1983conservedQ ).

### 4.1 Domain wall initial condition composed of two GE states

As described above, initially, there is a domain wall at the center of the
chain such that the left and the right halves are described by two GE states
with unequal inverse temperatures $\beta_{1}$ and $\beta_{2}$ respectively. We
first consider the simpler case of an infinite chain which will be followed by
a discussion on finite-size chain.

Infinite chain: We assume the domain wall to be initially located at $x=0$.
Using Eq. (19), we can write down the Wigner function at $t=0$ as

$\displaystyle
W(x,k,0)=\displaystyle{\frac{1}{\omega(k)}}\left[\frac{1}{\beta_{1}}-\left(\frac{1}{\beta_{1}}-\frac{1}{\beta_{2}}\right)\Theta\left(x\right)\right].$
(40)

From Eq. (26) the Wigner function at time $t$ can straightforwardly be written
as

$\displaystyle
W(x,k,t)=\displaystyle{\frac{1}{\omega(k)}}\left[\frac{1}{\beta_{1}}-\left(\frac{1}{\beta_{1}}-\frac{1}{\beta_{2}}\right)\Theta\left(x-\omega^{\prime}(k)t\right)\right].$
(41)

Substituting Eq. (41) into Eq. (27) and performing the integral, one can
compute the average internal energy density $e(x,t)$. Since we assume the
lattice spacing to be unity, the particle density is unity everywhere. Hence,
$e(x,t)$ is essentially the temperature profile which is explicitly given by

$\displaystyle
T(x,t):=e(x,t)=\begin{cases}\displaystyle{\frac{1}{\beta_{1}}},~{}x<-t\\\ \\\
\displaystyle{\frac{1}{2}}\left(\displaystyle{\frac{1}{\beta_{1}}}+\displaystyle{\frac{1}{\beta_{2}}}\right)-\displaystyle{\frac{1}{\pi}}\left(\displaystyle{\frac{1}{\beta_{1}}}-\displaystyle{\frac{1}{\beta_{2}}}\right)\sin^{-1}\left(\frac{x}{t}\right),~{}|x|<t\\\
\\\ \displaystyle{\frac{1}{\beta_{2}}},~{}x>t\end{cases}$ (42)

The expressions for other conserved quantities $\rho^{(n)}(x,t)$ for $n>0$ can
be calculated similarly and are given by

$\displaystyle\rho^{(n)}(x,t)=\displaystyle{\frac{1}{2n\pi}}\left(\displaystyle{\frac{1}{\beta_{1}}}-\displaystyle{\frac{1}{\beta_{2}}}\right)\sin\left\\{2n\cos^{-1}\left(\frac{x}{t}\right)\right\\}\Theta\left(1-\frac{|x|}{t}\right).$
(43)

Using the continuity equation in Eq. (38) for $n=0$, we can also write down
the energy current density

$\displaystyle\mathfrak{j}^{(0)}(x,t)=\displaystyle{\frac{1}{\pi}}\left(\displaystyle{\frac{1}{\beta_{1}}}-\displaystyle{\frac{1}{\beta_{2}}}\right)\sqrt{1-\left(\frac{x}{t}\right)^{2}}\Theta\left(1-\frac{|x|}{t}\right).$
(44)

The higher currents $\mathfrak{j}^{(n)}(x,t)$ for $n>0$ can be similarly
evaluated and one finds the following expression

$\displaystyle\mathfrak{j}^{(n)}(x,t)=\displaystyle{\frac{(-1)^{n+1}}{\pi(4n^{2}-1)}}\left(\displaystyle{\frac{1}{\beta_{1}}}-\displaystyle{\frac{1}{\beta_{2}}}\right)\left[2n\displaystyle{\frac{x}{t}}\sin\left(2n\mu\left(\frac{x}{t}\right)\right)+\sqrt{1-\left(\frac{x}{t}\right)^{2}}\cos\left(2n\mu\left(\frac{x}{t}\right)\right)\right]\Theta\left(1-\frac{|x|}{t}\right),$
(45)

where $\mu(z)=\sin^{-1}\left(z\right)$.

To compare the numerical microscopic computations with the theoretical results
from hydrodynamics, we coarse-grain the microscopic densities (and currents)
by averaging over $\Delta$ number of consecutive sites:

$\displaystyle\bar{\rho}^{(n)}_{y}(t)=\displaystyle{\frac{1}{\Delta}}\sum\limits_{j=(\ell-1)\Delta+1}^{\ell\Delta}\rho^{(n)}_{j}(t)~{}~{}\text{where}~{}~{}y=\ell\displaystyle{\frac{\Delta}{2}}.$
(46)

Note that we have expressed the coarse-grained densities as functions of
macroscopic time variable $t$, because, as mentioned earlier, in our case
$\epsilon=1$ which is small compared to $L=N$. In Fig. (2a) and Fig. (2b), we
plot the coarse-grained temperature and the corresponding current density
obtained numerically from the microscopic computations in a chain of length
$2N$. In the insets of these figures, we show that these profiles have
ballistic scaling and we compare the scaling functions with those (dashed
black lines) obtained analytically in Eqs. (42) and (44) for temperature and
the associated current respectively. The microscopic calculations are
performed by evaluating averages of the corresponding microscopic quantities
given in Eq. (31) and Eq. (33) and for that, we use the correlations from Eq.
(12). We see that the hydrodynamic results for both temperature and energy
current density match with the corresponding finite-size microscopic
computations as long as $t<N$. This can be understood since $N$ is the time
taken by the step initial condition to spread and reach the boundaries.

It is interesting to note that in the $t\to\infty$ limit, the currents at any
fixed point $x$ in the system approach stationary nonzero values given by

$\displaystyle\mathfrak{j}^{(n)}_{ss}=\displaystyle{\frac{(-1)^{n+1}}{\pi(4n^{2}-1)}}\left(\displaystyle{\frac{1}{\beta_{1}}}-\displaystyle{\frac{1}{\beta_{2}}}\right),$
(47)

and the relaxation happens as $\sim 1/t^{2}$. On the other hand, at any fixed
location in the $t\to\infty$ limit, all the densities of the conserved
quantities approach zero except temperature which approaches the mean
temperature. The relaxation to these stationary values occurs as $\sim 1/t$.
Therefore any finite segment of the infinite chain reaches a non-equilibrium
steady state (NESS) which cannot described by a GE but possibly by a current
carrying GGE instead, as shown for the quantum harmonic chain in eisler2014 .

As we will see below, this is different from the case of the finite chain
where the steady state is indeed expected to be described by a GE since the
currents, decay to zero due to reflections at the boundaries.

Figure 2: Plot showing the early time behavior of the coarse-grained
microscopic (a) temperature profile and (b) energy current density profile at
different times (in the scaled units $\tilde{t}=t/(2N)$). The insets in both
figures show the collapse of the profiles upon $x/t$ scaling of the x-axis
where $x=y-N$. The dashed black lines in the insets correspond to the analytic
scaling functions obtained from Eq. (42) and Eq. (44) respectively. The
coarse-graining was done over $\Delta=8$ consecutive sites. The parameters
used are $\beta_{1}=0.25,\beta_{2}=1.0$ and $2N=512$.

Finite chain: We now turn our attention to the case of a finite chain. Let the
length of the chain be $2L$ and the initial condition be that of a domain wall
at $x=L$ as given in Eq. (9) with
$\lambda_{0}^{L}=\beta_{1},\lambda_{0}^{R}=\beta_{2}$ and all other
$\lambda_{n}^{L,R}=0$. Using Eq. (19), we can write down the Wigner function
for the finite chain at $t=0$ as

$\displaystyle
W(x,k,0)=\displaystyle{\frac{1}{\omega(k)}}\left[\frac{1}{\beta_{2}}+\left(\frac{1}{\beta_{1}}-\frac{1}{\beta_{2}}\right)\left[\Theta\left(x+L\right)-\Theta\left(x-L\right)\right]\right],$
(48)

where $x\in[0,2L]$. Using Eq. (26) and taking care of reflections at the
boundaries (following a similar procedure as done in Appendix A of Ref.
chakraborti2022 ), the Wigner function at a later time $t$ is given by

$\displaystyle
W(x,k,t)=\displaystyle{\frac{1}{\omega(k)}}\left[\frac{1}{\beta_{2}}+\left(\frac{1}{\beta_{1}}-\frac{1}{\beta_{2}}\right)\sum_{s=-\infty}^{\infty}\left[\Theta\left(x-\omega^{\prime}(k)t+4sL+L\right)-\Theta\left(x-\omega^{\prime}(k)t+4sL-L\right)\right]\right].$
(49)

Substituting Eq. (49) into Eq. (37) with $n=0$, we obtain the following series
representation for the temperature

$\displaystyle
T(x,t)=\displaystyle{\frac{1}{\beta_{2}}}+\left(\frac{1}{\beta_{1}}-\frac{1}{\beta_{2}}\right)\sum\limits_{s=-\infty}^{\infty}\left[I\left(\frac{x+4sL+L}{t}\right)-I\left(\frac{x+4sL-L}{t}\right)\right],$
(50)

where we have defined the integral

$\displaystyle
I(z):=\int\limits_{-1}^{1}\displaystyle{\frac{du}{\pi}}\frac{\Theta(z-u)}{\sqrt{1-u^{2}}}=\begin{cases}0,~{}z<-1\\\
\\\ \displaystyle{\frac{1}{2}}+\frac{1}{\pi}\sin^{-1}z,~{}|z|<1\\\ \\\
1,~{}z>1\end{cases}$ (51)

Note that Eq. (49) can alternatively be written in a compact form

$\displaystyle
W(x,k,t)=\displaystyle{\frac{1}{\omega(k)}}\left[\frac{1}{\beta_{1}}-\left(\frac{1}{\beta_{1}}-\frac{1}{\beta_{2}}\right)\Theta\left(\left\\{\displaystyle{\frac{x-\omega^{\prime}(k)t+L}{4L}}\right\\}-\frac{1}{2}\right)\right],$
(52)

where $\\{x\\}=x-\lfloor x\rfloor$ is the fractional part function. Defining
$z=\left\\{\displaystyle{\frac{x-\omega^{\prime}(k)t+L}{4L}}\right\\}$, we
write the Fourier series representation,
$\Theta\left(z-\displaystyle{\frac{1}{2}}\right)=\displaystyle{\frac{1}{2}}-\sum\limits_{\ell\neq
0}\frac{(e^{\pi i\ell}-1)e^{-2\pi i\ell z}}{2\pi i\ell}$. Putting this in Eq.
(52) gives

$\displaystyle
W(x,k,t)=\displaystyle{\frac{1}{2\omega(k)}\left(\frac{1}{\beta_{1}}+\frac{1}{\beta_{2}}\right)}+\frac{1}{\omega(k)}\left(\frac{1}{\beta_{1}}-\frac{1}{\beta_{2}}\right)\sum\limits_{\ell\neq
0}\frac{\sin(\pi\ell/2)}{(\pi\ell)}\exp\left(-2\pi
i\ell\frac{x-\omega^{\prime}(k)t}{4L}\right).$ (53)

Now substituting Eq. (53) into Eq. (37) after setting $n=0$ and integrating
over $k$ yields an alternate series representation for the temperature profile

$\displaystyle
T(x,t)=\rho^{(0)}(x,t)=\displaystyle{\frac{1}{2}}\left(\frac{1}{\beta_{1}}+\frac{1}{\beta_{2}}\right)+\left(\frac{1}{\beta_{1}}-\frac{1}{\beta_{2}}\right)\sum\limits_{\ell=1}^{\infty}\frac{\sin(\pi\ell/2)}{(\pi\ell/2)}\cos\frac{\pi\ell
x}{2L}J_{0}\left(\frac{\pi\ell t}{2L}\right),$ (54)

where $J_{\nu}(z)$ is the Bessel function. From Eq. (53), we can also compute
all other densities $\rho^{(n)}(x,t)$ for $n>0$ and we find the following
evolution

$\displaystyle\rho^{(n)}(x,t)=(-1)^{n}\left(\frac{1}{\beta_{1}}-\frac{1}{\beta_{2}}\right)\sum\limits_{\ell=1}^{\infty}\frac{\sin(\pi\ell/2)}{(\pi\ell/2)}\cos\frac{\pi\ell
x}{2L}J_{2n}\left(\frac{\pi\ell t}{2L}\right).$ (55)

Similarly putting Eq. (53) into Eq. (39) and integrating over $k$ gives us the
current densities

$\displaystyle\mathfrak{j}^{(0)}(x,t)$
$\displaystyle=\left(\frac{1}{\beta_{1}}-\frac{1}{\beta_{2}}\right)\sum\limits_{\ell=1}^{\infty}\frac{\sin(\pi\ell/2)}{(\pi\ell/2)}\sin\frac{\pi\ell
x}{2L}J_{1}\left(\frac{\pi\ell t}{2L}\right),$ (56a)
$\displaystyle\mathfrak{j}^{(n)}(x,t)$
$\displaystyle=\displaystyle{\frac{(-1)^{n}}{2}}\left(\frac{1}{\beta_{1}}-\frac{1}{\beta_{2}}\right)\sum\limits_{\ell=1}^{\infty}\frac{\sin(\pi\ell/2)}{(\pi\ell/2)}\sin\frac{\pi\ell
x}{2L}\left[J_{2n+1}\left(\frac{\pi\ell
t}{2L}\right)-J_{2n-1}\left(\frac{\pi\ell t}{2L}\right)\right],~{}n>0.$ (56b)

The Bessel function $J_{\nu}(z)$ decays as $1/\sqrt{z}$ for large $z$.
Therefore, the temperature goes to a steady state value of
$(1/2)(\beta_{1}^{-1}+\beta_{2}^{-1})$ while all other quantities in Eqs. (55,
56) go to zero. All quantities relax as $1/\sqrt{t}$ in the large $t$ limit.

In Figs. (3-5), we show a comparison of the analytic results from
hydrodynamics, for temperature, energy current density, and $Q^{(1)}$ density,
with those obtained numerically from the microscopics by computing averages of
the corresponding microscopic quantities in Eqs. (31, 33) and using the
correlations from Eq. (12).

In Fig. (3a) we show the temperature profiles at different scaled times
$\tilde{t}=t/2N$. The initial front spreads freely as long as $t<N$. At $t=N$
the front reaches the boundaries for the first time and gets reflected back
into the bulk. These reflections continue to even out the temperature in the
two halves as can be observed from the profile corresponding to the largest
time $\tilde{t}=9.68$ shown in (3a). The temperature eventually becomes flat
and equal to the mean of the initial temperatures. In Fig. (3b), we show how
the temperature evolves with time at a coarse-grained point
$y\approx(2N)\times 0.46$ that is slightly to the left of the center. The
temperature remains unchanged until $\tilde{t}=(N-y)/(2N)\approx 0.04$ which
is the time taken by the edge of the front to reach that point (not visible in
the scale of the plot). After that it starts decaying as $\sim 1/t$ [see inset
of Fig. (3b)] similar to the infinite chain case and continues to do so until
an edge of the reflected front [in Fig. (3a)] crosses the particular coarse-
grained point for the first time. The temperature at any given point continues
to show an overall decay interspersed with the presence of isolated kinks that
appear corresponding to the passing of an edge of the reflected front.
Eventually, the temperature at a fixed coarse-grained point approaches the
final value as $\sim 1/\sqrt{t}$ (as also expected from our hydrodynamic
calculations). In both figures, points represent numerical data, and the solid
lines are analytical results from Eq. (54) and we observe an excellent
agreement between the two, thus demonstrating the validity of the GHD results
for the case of finite-size chains.

In Fig. (4a), we plot the energy current density profiles at a few different
times. The current profile is initially zero everywhere and starts to grow
from the center. As long as $t<N$, the current evolves according to Eq. (44)
after which the reflections from the boundaries become relevant. Fig. (4b)
shows the evolution of energy current density at a coarse-grained point that
lies near the middle of the left half. The current stays zero until
$\tilde{t}=(N-y)/(2N)\approx 0.27$ which is the time taken by the edge of the
front to reach that point. Here also we observe kinks that correspond to the
passing of an edge of the reflected front of the current profile in Fig. (4a).
The current eventually approaches zero as $\sim 1/\sqrt{t}$ [see inset of Fig.
(4b)]. We again see good agreement between the numerical microscopic
computations and the analytical results from hydrodynamics.

(a)(b)

Figure 3: (a) Plot showing the coarse-grained temperature profile at different
times $\tilde{t}=t/(2N)$. (b) shows the evolution of the temperature with time
at a fixed coarse-grained point $y/(2N)\approx 0.46$. The inset in (b) shows
the same curve on a log-log scale but with the steady state value subtracted.
We observe a $\sim 1/t$ decay at small times which, after many reflections
from the boundaries, changes to a $\sim 1/\sqrt{t}$ decay at large times as
expected from the analytical expressions. In both figures, the points are
calculated by coarse-graining the microscopic temperature profile whereas the
continuous lines are obtained from the analytic GHD solutions. We find an
excellent agreement between the coarse-grained microscopic and the GHD
results. The coarse-graining is done over $\Delta=8$ consecutive sites. The
parameters used are $\beta_{1}=0.25,\beta_{2}=1.0$ and $2N=512$.

(a)(b)

Figure 4: (a) Plot showing the coarse-grained energy current density profile
at different times $\tilde{t}=t/(2N)$. (b) shows the evolution of the energy
current density with time at a fixed coarse-grained point $y/(2N)\approx
0.23$. The inset in (b) shows the same curve on a log-log scale. We observe a
$\sim 1/\sqrt{t}$ decay at large times as expected from the analytical
expressions. In both figures, the points are calculated by coarse-graining the
microscopic energy current density profile whereas the continuous lines are
obtained from the analytic GHD expressions. We find an excellent agreement
between the coarse-grained microscopic and the GHD results. The coarse-
graining is done over $\Delta=8$ consecutive sites. The parameters used are
$\beta_{1}=0.25,\beta_{2}=1.0$ and $2N=512$.

(a)(b)

Figure 5: (a) Plot showing the coarse-grained $Q^{(1)}$ density profile at
different times $\tilde{t}=t/(2N)$. (b) shows the evolution of the $Q^{(1)}$
density with time at a fixed coarse-grained point $y/(2N)\approx 0.30$. The
inset in (b) shows the same curve on a log-log scale. We again observe a $\sim
1/\sqrt{t}$ decay at large times as expected from the analytical expressions.
In both figures, the points are calculated by coarse-graining the microscopic
$Q^{(1)}$ density profile whereas the continuous lines are obtained from the
analytic GHD expressions. We find an excellent agreement between the coarse-
grained microscopic and the GHD results. The coarse-graining is done over
$\Delta=8$ consecutive sites. The parameters used are
$\beta_{1}=0.25,\beta_{2}=1.0$ and $2N=512$.

### 4.2 Domain wall initial condition composed of two GGE states

We now consider a more general initial condition where the two halves are
described by a product of two GGEs, as given in Eq. (9). The initial Wigner
functions in the left and right parts are given by Eq. (20) with the Lagrange
parameters $\\{\lambda_{n}^{L}\\}$ and $\\{\lambda_{n}^{R}\\}$ respectively.
In the following, we focus on the particular case where only the first two
parameters are non-zero and the rest are zero on both sides of the chain.

Infinite chain: Once again we assume that the domain wall is initially located
at $x=0$ and the first two Lagrange multipliers have values
$\lambda_{0}^{L}=\beta_{1},~{}\lambda_{1}^{L}=-\beta_{1}\gamma_{1},~{}\lambda_{0}^{R}=\beta_{2}$,
and $\lambda_{1}^{R}=-\beta_{2}\gamma_{2}$. The initial Wigner function in
this case is

$\displaystyle
W(x,k,0)=\displaystyle{\frac{1}{\omega(k)}}\left[\frac{1}{\beta_{1}(1-\gamma_{1}\cos
k)}-\left(\frac{1}{\beta_{1}(1-\gamma_{1}\cos
k)}-\frac{1}{\beta_{2}(1-\gamma_{2}\cos
k)}\right)\Theta\left(x\right)\right],$ (57)

which according to Eq. (26) takes the following form at time $t$

$\displaystyle
W(x,k,t)=\displaystyle{\frac{1}{\omega(k)}}\left[\frac{1}{\beta_{1}(1-\gamma_{1}\cos
k)}-\left(\frac{1}{\beta_{1}(1-\gamma_{1}\cos
k)}-\frac{1}{\beta_{2}(1-\gamma_{2}\cos
k)}\right)\Theta\left(x-\omega^{\prime}(k)t\right)\right].$ (58)

Using this solution in Eq. (37) and simplifying we find that the density
profiles associated with the conserved quantities are given by

$\displaystyle\rho^{(n)}(x,t)=\begin{cases}g^{(n)}(\beta_{1},\gamma_{1}),~{}x<-t\\\
\\\
\displaystyle{\frac{1}{2}}[g^{(n)}(\beta_{1},\gamma_{1})+g^{(n)}(\beta_{2},\gamma_{2})]\\\
-\bigintss\limits_{0}^{\sin^{-1}(x/t)}\displaystyle{\frac{d\theta}{\pi}}\cos\\{2n\cos^{-1}|\sin\theta|\\}\left(\frac{1}{\beta_{1}(1+\gamma_{1}\cos
2\theta)}-\frac{1}{\beta_{2}(1+\gamma_{2}\cos 2\theta)}\right),~{}|x|<t\\\ \\\
g^{(n)}(\beta_{2},\gamma_{2}),~{}x>t\end{cases}$ (59)

where the quantity $g^{(n)}(\beta,\gamma)$ is given by

$\displaystyle
g^{(n)}(\beta,\gamma)=\displaystyle{\frac{1}{\beta(1+\gamma)}}{}_{3}F^{\text{reg}}_{2}\left(\\{1/2,1,1\\};\\{1-n,1+n\\};\displaystyle{\frac{2\gamma}{1+\gamma}}\right),$
(60)

with
${}_{p}F^{\text{reg}}_{q}(\\{a_{1},\ldots,a_{p}\\};\\{b_{1},\ldots,b_{q}\\};z)$
being the regularized generalized Hypergeometric function. The details of the
calculation are given in Appendix (B). Putting $n=0$ in Eq. (59) and
performing the integrals we find that the temperature profile is given by

$\displaystyle
T(x,t)=\rho^{(0)}(x,t)=\begin{cases}g^{(0)}(\beta_{1},\gamma_{1}),&x<-t\\\ \\\
\displaystyle{\frac{1}{2}}\left(g^{(0)}(\beta_{1},\gamma_{1})+g^{(0)}(\beta_{2},\gamma_{2})\right)-\displaystyle{\frac{1}{\pi}}\left[g^{(0)}(\beta_{1},\gamma_{1})\tan^{-1}\left\\{u(\gamma_{1})s(x/t)\right\\}\right.\\\
\left.-g^{(0)}(\beta_{2},\gamma_{2})\tan^{-1}\left\\{u(\gamma_{2})s(x/t)\right\\}\right],&|x|<t\\\
\\\ g^{(0)}(\beta_{2},\gamma_{2}),&x>t\end{cases}$ (61)

where $g^{(0)}(\beta,\gamma)=\frac{1}{\beta\sqrt{1-\gamma^{2}}}$ and $s(z)$
and $u(\gamma)$ are defined as

$\displaystyle
s(z)=\frac{z}{\sqrt{1-z^{2}}},~{}u(\gamma)=\sqrt{\displaystyle{\frac{1-\gamma}{1+\gamma}}}.$
(62)

The expression for the density corresponding to $Q^{(1)}$ is similarly given
by

$\displaystyle\rho^{(1)}(x,t)=\begin{cases}g^{(1)}(\beta_{1},\gamma_{1}),&x<-t\\\
\\\
\displaystyle{\frac{1}{2}}\left(g^{(1)}(\beta_{1},\gamma_{1})+g^{(1)}(\beta_{2},\gamma_{2})\right)\\\
-\displaystyle{\frac{1}{\pi}}\left[g^{(1)}(\beta_{1},\gamma_{1})\tan^{-1}\left\\{u(\gamma_{1})s(x/t)\right\\}-g^{(1)}(\beta_{2},\gamma_{2})\tan^{-1}\left\\{u(\gamma_{2})s(x/t)\right\\}\right]&|x|<t\\\
+\displaystyle{\frac{1}{\pi}}\left[\left(\displaystyle{\frac{1}{\beta_{1}\gamma_{1}}}\tan^{-1}\left\\{\displaystyle{\frac{[1-u(\gamma_{1})]s(x/t)}{1+u(\gamma_{1})s(x/t)^{2}}}\right\\}-\frac{1}{\beta_{2}\gamma_{2}}\tan^{-1}\left\\{\displaystyle{\frac{[1-u(\gamma_{2})]s(x/t)}{1+u(\gamma_{2})s(x/t)^{2}}}\right\\}\right)\right],\\\
\\\ g^{(1)}(\beta_{2},\gamma_{2}),&x>t\end{cases}$ (63)

where
$g^{(1)}(\beta,\gamma)=\frac{1}{\beta\gamma}\left(\frac{1}{\sqrt{1-\gamma^{2}}}-1\right)$
and $s(z)$ and $u(\gamma)$ are defined in Eq. (62).

In Fig. (6), we show a comparison of the analytic results from hydrodynamics
with those obtained microscopically for a chain of length $2N$. The averages
of the corresponding microscopic quantities in Eq. (31) are evaluated by using
the correlations from Eq. (12). We see that the hydrodynamic results in Eqs.
(61) and (63) for temperature and $Q^{(1)}$ density, respectively, match with
the corresponding finite-size microscopic computations as long as $t<N$.

Figure 6: Plot showing the early time behavior of the coarse-grained
microscopic (a) temperature profile and (b) $Q^{(1)}$ density profile at
different times (in the scaled units $\tilde{t}=t/(2N)$). The insets in both
figures show the collapse of the profiles upon $x/t$ scaling of the x-axis
where $x=y-N$. The dashed black lines in the insets correspond to the analytic
scaling functions obtained from Eq. (61) and Eq. (63) respectively. The
coarse-graining was done over $\Delta=8$ consecutive sites. The parameters
used are $\beta_{1}=0.25,\beta_{2}=1.0,\gamma_{1}=0.4,\gamma_{2}=0.8$ and
$2N=512$.

Finite chain: We now turn our attention to the case of a finite chain. Let the
length of the chain be $2L$ and the initial condition be that of a domain wall
at $x=L$ as given in Eq. (9) with
$\lambda_{0}^{L}=\beta_{1},\lambda_{1}^{L}=-\beta_{1}\gamma_{1},\lambda_{0}^{R}=\beta_{2},\lambda_{1}^{R}=-\beta_{2}\gamma_{2}$
and all other $\lambda_{n}^{L,R}=0$. Using Eq. (20) with these parameters in
the two halves of the chain, we can write down the Wigner function at $t=0$ as

$\displaystyle
W(x,k,0)=\displaystyle{\frac{1}{\omega(k)}}\left[\frac{1}{\beta_{2}(1-\gamma_{2}\cos
k)}+\left(\frac{1}{\beta_{1}(1-\gamma_{1}\cos
k)}-\frac{1}{\beta_{2}(1-\gamma_{2}\cos
k)}\right)\left[\Theta\left(x+L\right)-\Theta\left(x-L\right)\right]\right],$
(64)

where $x\in[0,2L]$. As was argued in Sec. (4.1) for the finite chain, we can
similarly write down the Wigner function at time $t$ as

$\displaystyle\begin{split}W(x,k,t)=\displaystyle{\frac{1}{\omega(k)}}&\left[\frac{1}{\beta_{2}(1-\gamma_{2}\cos
k)}+\left(\frac{1}{\beta_{1}(1-\gamma_{1}\cos
k)}-\frac{1}{\beta_{2}(1-\gamma_{2}\cos k)}\right)\right.\\\
&\left.\times\sum_{s=-\infty}^{\infty}\left[\Theta\left(x-\omega^{\prime}(k)t+4sL+L\right)-\Theta\left(x-\omega^{\prime}(k)t+4sL-L\right)\right]\right].\end{split}$
(65)

Substituting Eq. (65) into Eq. (37) with $n=0$, we obtain the following series
representation for the temperature

$\displaystyle
T(x,t)=\displaystyle{\frac{1}{\beta_{2}\sqrt{1-\gamma_{1}^{2}}}}+\sum\limits_{s=-\infty}^{\infty}\left[\frac{1}{\beta_{1}}I_{\gamma_{1}}\left(\frac{x+4sL+L}{t}\right)-\frac{1}{\beta_{2}}I_{\gamma_{2}}\left(\frac{x+4sL-L}{t}\right)\right],$
(66)

where $I_{\gamma}(z)$ is defined as the integral

$\displaystyle
I_{\gamma}(z):=\int\limits_{-1}^{1}\displaystyle{\frac{du}{\pi}}\frac{\Theta(z-u)}{\sqrt{1-u^{2}}[1-\gamma(2u^{2}-1)]}.$
(67)

Note that, as before, Eq. (65) can alternatively be written in a compact form

$\displaystyle
W(x,k,t)=\displaystyle{\frac{1}{\omega(k)}}\left[\frac{1}{\beta_{1}(1-\gamma_{1}\cos
k)}-\left(\frac{1}{\beta_{1}(1-\gamma_{1}\cos
k)}-\frac{1}{\beta_{2}(1-\gamma_{2}\cos
k)}\right)\Theta\left(\left\\{\displaystyle{\frac{x-\omega^{\prime}(k)t+L}{4L}}\right\\}-\frac{1}{2}\right)\right],$
(68)

As done in Sec. (4.1) for the finite case, we define
$z=\left\\{\displaystyle{\frac{x-\omega^{\prime}(k)t+L}{4L}}\right\\}$ and
write the Fourier series representation,
$\Theta\left(z-\displaystyle{\frac{1}{2}}\right)=\displaystyle{\frac{1}{2}}-\sum\limits_{\ell\neq
0}\frac{(e^{\pi i\ell}-1)e^{-2\pi i\ell z}}{2\pi i\ell}$. Putting this in Eq.
(68) gives

$\displaystyle\begin{split}W(x,k,t)&=\displaystyle{\frac{1}{2\omega(k)}\left(\frac{1}{\beta_{1}(1-\gamma_{1}\cos
k)}+\frac{1}{\beta_{2}(1-\gamma_{2}\cos k)}\right)}\\\
&+\frac{1}{\omega(k)}\left(\frac{1}{\beta_{1}(1-\gamma_{1}\cos
k)}-\frac{1}{\beta_{2}(1-\gamma_{2}\cos k)}\right)\sum\limits_{\ell\neq
0}\frac{\sin(\pi\ell/2)\exp\left(-2\pi
i\ell\frac{x-\omega^{\prime}(k)t}{4L}\right)}{\pi\ell}.\end{split}$ (69)

Substituting Eq. (69) into Eq. (37), the evolution for $\rho^{(n)}(x,t)$ can
be formally written down as

$\displaystyle\begin{split}\rho^{(n)}(x,t)&=\rho^{(n)}_{GGE}+\sum\limits_{\ell=1}^{\infty}\frac{\sin(\pi\ell/2)}{(\pi\ell/2)}\cos\frac{\pi\ell
x}{2L}\left[\displaystyle{\frac{1}{\beta_{1}}}C^{(n)}_{\ell}(\gamma_{1},t)-\displaystyle{\frac{1}{\beta_{2}}}C^{(n)}_{\ell}(\gamma_{2},t)\right],\end{split}$
(70)

where the stationary state density $\rho^{(n)}_{GGE}$ is given by

$\displaystyle\rho^{(n)}_{GGE}=\displaystyle{\frac{1}{2}}\left[g^{(n)}(\beta_{1},\gamma_{1})+g^{(n)}(\beta_{2},\gamma_{2})\right],$
(71)

with $g^{(n)}(\beta,\gamma)$ defined in Eq. (60) and
$C^{(n)}_{\ell}(\gamma,t)$ given by the integral

$\displaystyle
C^{(n)}_{\ell}(\gamma,t)=\int\limits_{0}^{\pi}\displaystyle{\frac{dk}{2\pi}}\frac{2\cos
nk}{1-\gamma\cos k}\cos\left(\frac{\pi\ell t}{2L}\cos\frac{k}{2}\right).$ (72)

Similarly, we can write down the formal expression for the currents
$\mathfrak{j}^{(n)}(x,t)$

$\displaystyle\mathfrak{j}^{(n)}(x,t)=\sum\limits_{\ell=1}^{\infty}\frac{\sin(\pi\ell/2)}{(\pi\ell/2)}\sin\frac{\pi\ell
x}{2L}\left[\displaystyle{\frac{1}{\beta_{1}}}D^{(n)}_{\ell}(\gamma_{1},t)-\displaystyle{\frac{1}{\beta_{2}}}D^{(n)}_{\ell}(\gamma_{2},t)\right],$
(73)

where $D^{(n)}_{\ell}(\gamma,t)$ is given by the integral

$\displaystyle
D^{(n)}_{\ell}(\gamma,t)=\int\limits_{0}^{\pi}\displaystyle{\frac{dk}{2\pi}}\frac{2\cos
nk\cos(k/2)}{1-\gamma\cos k}\sin\left(\frac{\pi\ell
t}{2L}\cos\frac{k}{2}\right).$ (74)

It can be shown that, in the large time limit, the densities and currents in
Eq. (70) and Eq. (73) respectively relax to their equilibrium values as $\sim
1/\sqrt{t}$, similar to the finite-size case in Sec. (4.1). In Figs. (7-10),
we show a comparison of the analytic results from hydrodynamics, for
temperature, energy current density, $Q^{(1)}$ density, and $Q^{(1)}$ current
density with those obtained numerically from the microscopics by computing
averages of the corresponding microscopic quantities in Eqs. (31, 33) and
using the correlations in Eq. (12). The qualitative features of the evolution
of these quantities are similar to those for the finite chain in Sec. (4.1).
We again find a good agreement between hydrodynamics and exact numerics.

(a)(b)

Figure 7: (a) Plot showing the coarse-grained temperature profile at different
times $\tilde{t}=t/(2N)$. (b) shows the evolution of the temperature with time
at a fixed coarse-grained point $y/(2N)\approx 0.30$. The inset in (b) shows
the same curve on a log-log scale but with the steady state value subtracted.
We observe a $\sim 1/t$ decay at small times which, after many reflections
from the boundaries, changes to a $\sim 1/\sqrt{t}$ decay at large times as
expected from the analytical expressions. In both figures, the points are
calculated by coarse-graining the microscopic temperature profile whereas the
continuous lines are obtained from the analytic GHD expressions. We find an
excellent agreement between the coarse-grained microscopic and the GHD
results. The coarse-graining is done over $\Delta=8$ consecutive sites. The
parameters used are
$\beta_{1}=0.25,\beta_{2}=1.0,\gamma_{1}=0.4,\gamma_{2}=0.8$ and $2N=512$.

(a)(b)

Figure 8: (a) Plot showing the coarse-grained energy current density profile
at different times $\tilde{t}=t/(2N)$. (b) shows the evolution of the energy
current density with time at a fixed coarse-grained point $y/(2N)\approx
0.23$. The inset in (b) shows the same curve on a log-log scale. We observe a
$\sim 1/\sqrt{t}$ decay at large times as expected from the analytical
expressions. In both figures, the points are calculated by coarse-graining the
microscopic energy current density profile whereas the continuous lines are
obtained from the analytic GHD expressions. We find an excellent agreement
between the coarse-grained microscopic and the GHD results. The coarse-
graining is done over $\Delta=8$ consecutive sites. The parameters used are
$\beta_{1}=0.25,\beta_{2}=1.0,\gamma_{1}=0.4,\gamma_{2}=0.8$ and $2N=512$.

(a)(b)

Figure 9: (a) Plot showing the coarse-grained $Q^{(1)}$ density profile at
different times $\tilde{t}=t/(2N)$. (b) shows the evolution of the $Q^{(1)}$
density with time at a fixed coarse-grained point $y/(2N)\approx 0.30$. The
inset in (b) shows the same curve on a log-log scale but with the steady state
value subtracted. We observe a $\sim 1/\sqrt{t}$ decay at large times as
expected from the analytical expressions. In both figures, the points are
calculated by coarse-graining the microscopic $Q^{(1)}$ density profile
whereas the continuous lines are obtained from the analytic GHD expressions.
We find an excellent agreement between the coarse-grained microscopic and the
GHD results. The coarse-graining is done over $\Delta=8$ consecutive sites.
The parameters used are
$\beta_{1}=0.25,\beta_{2}=1.0,\gamma_{1}=0.4,\gamma_{2}=0.8$ and $2N=512$.

(a)(b)

Figure 10: (a) Plot showing the coarse-grained $Q^{(1)}$ current density
profile at different times $\tilde{t}=t/(2N)$. (b) shows the evolution of the
$Q^{(1)}$ current density with time at a fixed coarse-grained point
$y/(2N)\approx 0.23$. The inset in (b) shows the same curve on a log-log
scale. We observe a $\sim 1/\sqrt{t}$ decay at large times as expected from
the analytical expressions. In both figures, the points are calculated by
coarse-graining the microscopic $Q^{(1)}$ current density profile whereas the
continuous lines are obtained from the analytic GHD expressions. We find an
excellent agreement between the coarse-grained microscopic and the GHD
results. The coarse-graining is done over $\Delta=8$ consecutive sites. The
parameters used are
$\beta_{1}=0.25,\beta_{2}=1.0,\gamma_{1}=0.4,\gamma_{2}=0.8$ and $2N=512$.

We expect the final state to be a GGE state characterised by Lagrange
multipliers $\\{\lambda_{n}\\}$ which are related to the average densities as

$\displaystyle\rho^{(n)}_{GGE}=\int\limits_{-\pi}^{\pi}\displaystyle{\frac{dk}{2\pi}}\displaystyle{\frac{\cos
nk}{\sum\limits_{m=0}^{\infty}\lambda_{m}\cos mk}}.$ (75)

The expression in Eq. (75) can be easily inverted using the Fourier cosine
series and we obtain the following explit expressions for the Lagrange
multipliers of the finite chain

$\displaystyle\lambda_{0}$
$\displaystyle=\displaystyle{\frac{1}{2\pi}}\int\limits_{0}^{\pi}dk~{}\frac{1}{f(k)},$
(76a) $\displaystyle\lambda_{n}$
$\displaystyle=\displaystyle{\frac{1}{\pi}}\int\limits_{0}^{\pi}dk~{}\frac{\cos
nk}{f(k)}~{}~{}\text{for}~{}~{}n>1,$ (76b)

where $f(k)$ is given by

$\displaystyle
f(k)=\displaystyle{\frac{1}{2}}\rho^{(0)}_{GGE}+\sum\limits_{n=1}^{\infty}\rho^{(n)}_{GGE}\cos
nk.$ (77)

Putting the values of the conserved densities from Eq. (71) into Eq. (76), we
can compute the $\\{\lambda_{n}\\}$ in the final state. In Fig. (11), we plot
the stationary state density and the absolute values of the Lagrange
multipliers on a log-linear scale. Unlike the case where the initial state is
composed of two GEs, here we find that all densities and Lagrange multipliers
are non-zero in the final state. However, both quantities decay exponentially
for large $n$. Note that Eqs. (75-77) are valid even if the initial state is
described by a product of two GGEs with an arbitrary number of nonzero
parameters that have unequal values in the two halves of the chain.

Figure 11: Plot showing exponential decay of the densities $\rho^{(n)}_{GGE}$
and the absolute values of the Lagrange multipliers $\lambda_{n}$ with
increasing $n$. The points in red are calculated using Eq. (71) whereas the
blue points are evaluated using Eq. (76). The dashed lines show an exponential
fit to the data points. The parameters used are
$\beta_{1}=0.25,\beta_{2}=1.0,\gamma_{1}=0.4,\gamma_{2}=0.8$.

##  5 Conclusions

In this work, we first provided a physically motivated derivation of the
generalized hydrodynamics of a classical harmonic chain, based on the approach
of the Wigner function. Our work thus makes connections between earlier work
on harmonic chains dobrushin1986 ; dobrushin1988 ; boldrighini1983euler ;
boldrighini1997navierstokes  with the recent developments in GHD. As an
application, we studied the problem of thermalization of a harmonic chain
starting from domain wall initial conditions.

Specifically, we studied the evolution of a classical harmonic chain with
nearest-neighbor interactions starting from a domain wall initial state
composed of either two GEs or two GGEs with unequal chemical potentials in the
left and right halves of the chain. For each type of the initial distribution,
we study both finite and infinite chains. In all cases, we found excellent
agreement, at all times, between results from hydrodynamics with those
obtained from exact numerics. We observe that a finite chain relaxes as $\sim
1/\sqrt{t}$ to a steady state that is expected to be described by a GE or a
GGE depending on whether we use the GE or GGE type of initial condition. For
the GGE initial condition, the final state is also expected to be a GGE
characterized by a set of Lagrange multipliers which we evaluate exactly in
terms of the steady state densities. We find that this GGE steady state has an
infinite number of Lagrange parameters with exponentially decaying strengths.
All currents decay to zero as $\sim 1/\sqrt{t}$ as we approach the steady
state due to reflections at the boundaries. For the case of an infinite chain,
we find that even for the initial condition composed of two GEs, any finite
region of the system tends, in the $t\to\infty$ limit, towards a non-
equilibrium steady state (NESS) in which the currents have nonzero but
constant values. The NESS in the infinite chain is therefore expected to be a
current-carrying GGE. The relaxation to this NESS happens as $\sim 1/t$ for
the densities and as $\sim 1/t^{2}$ for the currents.

We note that the Wigner GHD equation is applicable for generic harmonic
systems incorporating more complex force matrices in arbitrary dimensions and
quantum statistics. These can be interesting extensions of the present work.
For the quantum case, the computation of entanglement entropy in free
fermionic systems using semiclassical hydrodynamics has generated a lot of
interest bertini2018GHD ; saha2024page  and it would be of great interest to
explore this in the context of quantum harmonic crystals. Another interesting
open question would be to prove or demonstrate equilibration in the harmonic
chain for a single realization that starts from an arbitrary initial
condition. This has recently been proved tasaki2024macroscopic  in the quantum
case for non-interacting fermions while for a classical ideal gas, it has been
established for initial conditions chosen from continuous distributions
deBievre2017Kac .

##  6 Acknowledgements

We thank Herbert Spohn for very useful discussions. We acknowledge support
from the Department of Atomic Energy, Government of India, under Project No.
RTI4001. AD acknowledges the J.C. Bose Fellowship (JCB/2022/000014) of the
Science and Engineering Research Board of the Department of Science and
Technology, Government of India. A. K. would like to acknowledge the support
of DST, Government of India Grant under Project No. ECR/2017/000634 and the
MATRICS grant MTR/2021/000350 from the SERB, DST, Government of India.

##  Appendix A GGE partition function

The GGE partition function in Eq. (7) is given by the integral

$\displaystyle\begin{split}Z(\\{\lambda_{n}\\})&:=\int_{-\infty}^{\infty}\prod_{j=1}^{N}dp_{j}dq_{j}\exp\left(-\sum\limits_{n=0}^{N-1}\lambda_{n}Q^{(n)}({\bf
p},{\bf q})\right)\\\
&=\int_{-\infty}^{\infty}\prod_{j=1}^{N}dp_{j}dq_{j}\exp\left(-\frac{1}{2}\left[{\bf{p}}^{T}\sum_{n}\lambda_{n}B^{(n)}{\bf{p}}+{\bf{q}}^{T}\sum_{n}\lambda_{n}M^{(n)}{\bf{q}}\right]\right)\end{split}$
(78)

where we have used the expression for $Q^{(n)}$ in Eq. (34) and the Lagrange
multipliers are chosen such that $\sum_{n}\lambda_{n}B_{n}$ and
$\sum_{n}\lambda_{n}M_{n}$ are positive definite matrices. The Gaussian
integration can be easily performed once we know the eigenvalues of $B^{(n)}$
and $M^{(n)}$ matrices which are respectively given by $\cos(nk_{\ell})$ and
$\omega_{k_{\ell}}^{2}\cos(nk_{\ell})$ with
$k_{\ell}=\displaystyle{\frac{\pi\ell}{N+1}},~{}\ell=1,\ldots,N$ and
$\omega_{k}^{2}=2(1-\cos k)$. We therefore get

$\displaystyle
Z(\\{\lambda_{n}\\})=\sqrt{\displaystyle{\frac{(2\pi)^{N}}{\prod_{\ell=1}^{N}\left(\sum_{n=0}^{N-1}\lambda_{n}\cos(nk_{\ell})\right)}}}\times\sqrt{\displaystyle{\frac{(2\pi)^{N}}{\prod_{\ell=1}^{N}\omega_{k_{\ell}}^{2}\left(\sum_{n=0}^{N-1}\lambda_{n}\cos(nk_{\ell})\right)}}},$
(79)

which can be simplified to yield Eq. (8).

##  Appendix B Density evolution corresponding to the initial condition
composed of two GGEs

Here we show the steps leading to the density profiles in Eq. (59) starting
from the expression for the Wigner function in Eq. (58). Substituting the
Wigner function from Eq. (58) into Eq. (37), we obtain the following
expression for $\rho^{(n)}(x,t)$

$\displaystyle\rho^{(n)}(x,t)=2\int\limits_{0}^{\pi/2}\displaystyle{\frac{d\theta}{\pi}}\cos
2n\theta\frac{1}{\beta_{1}(1-\gamma_{1}\cos
2\theta)}-I^{(n)}(x,t)=g^{(n)}(\beta_{1},\gamma_{1})-I^{(n)}(x,t),$ (80)

where $g^{(n)}(\beta,\gamma)$ is defined in Eq. (60) and $I^{(n)}(x,t)$ is
given by

$\displaystyle
I^{(n)}(x,t)=\int\limits_{0}^{\pi/2}\displaystyle{\frac{d\theta}{\pi}}\cos
2n\theta\left(\frac{1}{\beta_{1}(1-\gamma_{1}\cos
2\theta)}-\frac{1}{\beta_{2}(1-\gamma_{2}\cos
2\theta)}\right)\left[\Theta(x-t\cos\theta)+\Theta(x+t\cos\theta)\right].$
(81)

Differentiating Eq. (81) w.r.t $x$ on both sides and simplifying yields

$\displaystyle\displaystyle{\frac{\partial I^{(n)}}{\partial
x}}=\displaystyle{\frac{\cos 2n\theta^{\star}}{\pi
t\sin\theta^{\star}}}\left(\frac{1}{\beta_{1}(1-\gamma_{1}\cos
2\theta^{\star})}-\frac{1}{\beta_{2}(1-\gamma_{2}\cos
2\theta^{\star})}\right),~{}|x|<t$ (82)

where $\theta^{\star}=\cos^{-1}(|x|/t)$. Integrating Eq. (82) w.r.t $x$ gives

$\displaystyle I^{(n)}(x,t)=\begin{cases}0,~{}x<-t\\\ \\\
\bigintss\limits_{0}^{\pi/2}\displaystyle{\frac{d\theta}{\pi}}\cos
2n\theta\left(\frac{1}{\beta_{1}(1-\gamma_{1}\cos
2\theta)}-\frac{1}{\beta_{2}(1-\gamma_{2}\cos 2\theta)}\right)\\\
+\bigintss\limits_{0}^{\sin^{-1}(x/t)}\displaystyle{\frac{d\theta}{\pi}}\cos\\{2n\cos^{-1}|\sin\theta|\\}\left(\frac{1}{\beta_{1}(1+\gamma_{1}\cos
2\theta)}-\frac{1}{\beta_{2}(1+\gamma_{2}\cos 2\theta)}\right),~{}|x|<t\\\ \\\
2\bigintss\limits_{0}^{\pi/2}\displaystyle{\frac{d\theta}{\pi}}\cos
2n\theta\left[\left(\frac{1}{\beta_{1}(1-\gamma_{1}\cos
2\theta)}-\frac{1}{\beta_{2}(1-\gamma_{2}\cos
2\theta)}\right)\right],~{}x>t\end{cases}$ (83)

Using Eq. (60), we can rewrite Eq. (83) as

$\displaystyle I^{(n)}(x,t)=\begin{cases}0,~{}x<-t\\\ \\\
\displaystyle{\frac{1}{2}}[g^{(n)}(\beta_{1},\gamma_{1})-g^{(n)}(\beta_{2},\gamma_{2})]\\\
+\bigintss\limits_{0}^{\sin^{-1}(x/t)}\displaystyle{\frac{d\theta}{\pi}}\cos\\{2n\cos^{-1}|\sin\theta|\\}\left(\frac{1}{\beta_{1}(1+\gamma_{1}\cos
2\theta)}-\frac{1}{\beta_{2}(1+\gamma_{2}\cos 2\theta)}\right),~{}|x|<t\\\ \\\
g^{(n)}(\beta_{1},\gamma_{1})-g^{(n)}(\beta_{2},\gamma_{2}),~{}x>t\end{cases}$
(84)

Substituting Eq. (84) into Eq. (80) gives us the required expression in Eq.
(59).

## References

  * (1) E. Fermi, P. Pasta, S. Ulam, and M. Tsingou, “Studies of the nonlinear problems,” tech. rep., Los Alamos National Lab.(LANL), Los Alamos, NM (United States), 1955.
  * (2) T. Dauxois, “Fermi, Pasta, Ulam, and a mysterious lady,” Physics Today, vol. 61, no. 1, pp. 55–57, 2008.
  * (3) S. Chakraborti, A. Dhar, S. Goldstein, A. Kundu, and J. L. Lebowitz, “Entropy growth during free expansion of an ideal gas,” Journal of Physics A: Mathematical and Theoretical, vol. 55, no. 39, p. 394002, 2022.
  * (4) N. Cocciaglia, A. Vulpiani, and G. Gradenigo, “Thermalization without chaos in harmonic systems,” Physica A: Statistical Mechanics and its Applications, vol. 601, p. 127581, 2022.
  * (5) S. Chakraborti, A. Dhar, and A. Kundu, “Boltzmann’s entropy during free expansion of an interacting gas,” Journal of Statistical Physics, vol. 190, no. 4, p. 74, 2023.
  * (6) S. Ganapa, A. Apte, and A. Dhar, “Thermalization of local observables in the $\alpha$-FPUT chain,” Journal of Statistical Physics, vol. 180, no. 1, pp. 1010–1030, 2020.
  * (7) N. Shiraishi and H. Tasaki, “Nature abhors a vacuum: A simple rigorous example of thermalization in an isolated macroscopic quantum system,” arXiv preprint arXiv:2310.18880, 2023.
  * (8) S. Pandey, J. M. Bhat, A. Dhar, S. Goldstein, D. A. Huse, M. Kulkarni, A. Kundu, and J. L. Lebowitz, “Boltzmann entropy of a freely expanding quantum ideal gas,” Journal of Statistical Physics, vol. 190, no. 8, p. 142, 2023.
  * (9) M. Rigol, V. Dunjko, V. Yurovsky, and M. Olshanii, “Relaxation in a completely integrable many-body quantum system: an ab initio study of the dynamics of the highly excited states of 1d lattice hard-core bosons,” Physical Review Letters, vol. 98, no. 5, p. 050405, 2007.
  * (10) A. C. Cassidy, C. W. Clark, and M. Rigol, “Generalized thermalization in an integrable lattice system,” Physical Review Letters, vol. 106, no. 14, p. 140405, 2011.
  * (11) P. Calabrese, F. H. Essler, and M. Fagotti, “Quantum quench in the transverse-field Ising chain,” Physical Review Letters, vol. 106, no. 22, p. 227203, 2011.
  * (12) F. H. Essler, “A short introduction to generalized hydrodynamics,” Physica A: Statistical Mechanics and its Applications, p. 127572, 2022.
  * (13) S. K. Singh, A. Dhar, H. Spohn, and A. Kundu, “Thermalization and hydrodynamics in an interacting integrable system: the case of hard rods,” arXiv preprint arXiv:2310.18684, 2023.
  * (14) R. Dobrushin, A. Pellegrinotti, Y. M. Suhov, and L. Triolo, “One-dimensional harmonic lattice caricature of hydrodynamics,” Journal of Statistical Physics, vol. 43, pp. 571–607, 1986.
  * (15) A. Mielke, “Macroscopic behavior of microscopic oscillations in harmonic lattices via Wigner-Husimi transforms,” Archive for Rational Mechanics and Analysis, vol. 181, pp. 401–448, 2006.
  * (16) R. Dobrushin, A. Pellegrinotti, Y. M. Suhov, and L. Triolo, “One-dimensional harmonic lattice caricature of hydrodynamics: second approximation,” Journal of Statistical Physics, vol. 52, pp. 423–439, 1988.
  * (17) J. Percus, “Exact solution of kinetics of a model classical fluid,” The Physics of Fluids, vol. 12, no. 8, pp. 1560–1563, 1969.
  * (18) C. Boldrighini, R. Dobrushin, and Y. M. Sukhov, “One-dimensional hard rod caricature of hydrodynamics,” Journal of Statistical Physics, vol. 31, pp. 577–616, 1983.
  * (19) C. Boldrighini and Y. M. Suhov, “One-dimensional hard-rod caricature of hydrodynamics:“Navier–Stokes correction” for local equilibrium initial states,” Communications in Mathematical Physics, vol. 189, no. 2, pp. 577–590, 1997.
  * (20) O. A. Castro-Alvaredo, B. Doyon, and T. Yoshimura, “Emergent hydrodynamics in integrable quantum systems out of equilibrium,” Physical Review X, vol. 6, no. 4, p. 041065, 2016.
  * (21) B. Bertini, M. Collura, J. De Nardis, and M. Fagotti, “Transport in out-of-equilibrium XXZ chains: Exact profiles of charges and currents,” Physical Review Letters, vol. 117, no. 20, p. 207201, 2016.
  * (22) B. Doyon, “Lecture notes on generalized hydrodynamics,” SciPost Physics Lecture Notes, p. 018, 2020.
  * (23) H. Spohn, “Hydrodynamic equations for the Toda lattice,” arXiv preprint arXiv:2101.06528, 2021.
  * (24) X. Cao, V. B. Bulchandani, and H. Spohn, “The GGE averaged currents of the classical Toda chain,” Journal of Physics A: Mathematical and Theoretical, vol. 52, no. 49, p. 495003, 2019.
  * (25) G. Mazzuca, T. Grava, T. Kriecherbauer, K. T.-R. McLaughlin, C. B. Mendl, and H. Spohn, “Equilibrium spacetime correlations of the Toda lattice on the hydrodynamic scale,” Journal of Statistical Physics, vol. 190, no. 8, p. 149, 2023.
  * (26) A. Kundu, “Integrable hydrodynamics of Toda chain: case of small systems,” The European Physical Journal Special Topics, pp. 1–10, 2023.
  * (27) P. Ruggiero, P. Calabrese, B. Doyon, and J. Dubail, “Quantum generalized hydrodynamics,” Physical Review Letters, vol. 124, no. 14, p. 140603, 2020.
  * (28) B. Doyon and H. Spohn, “Dynamics of hard rods with initial domain wall state,” Journal of Statistical Mechanics: Theory and Experiment, vol. 2017, no. 7, p. 073210, 2017.
  * (29) J.-S. Caux, B. Doyon, J. Dubail, R. Konik, and T. Yoshimura, “Hydrodynamics of the interacting Bose gas in the Quantum Newton Cradle setup,” SciPost Physics, vol. 6, no. 6, p. 070, 2019.
  * (30) M. Schemmer, I. Bouchoule, B. Doyon, and J. Dubail, “Generalized hydrodynamics on an atom chip,” Physical Review Letters, vol. 122, no. 9, p. 090601, 2019.
  * (31) X. Cao, V. B. Bulchandani, and J. E. Moore, “Incomplete thermalization from trap-induced integrability breaking: Lessons from classical hard rods,” Physical Review Letters, vol. 120, no. 16, p. 164101, 2018.
  * (32) N. Malvania, Y. Zhang, Y. Le, J. Dubail, M. Rigol, and D. S. Weiss, “Generalized hydrodynamics in strongly interacting 1d Bose gases,” Science, vol. 373, no. 6559, pp. 1129–1133, 2021.
  * (33) V. Eisler and Z. Zimborás, “Entanglement negativity in the harmonic chain out of equilibrium,” New Journal of Physics, vol. 16, no. 12, p. 123020, 2014.
  * (34) T. Dudnikova, A. Komech, and H. Spohn, “On the convergence to statistical equilibrium for harmonic crystals,” Journal of Mathematical Physics, vol. 44, no. 6, pp. 2596–2620, 2003.
  * (35) C. Boldrighini, A. Pellegrinotti, and L. Triolo, “Convergence to stationary states for infinite harmonic systems,” Journal of Statistical Physics, vol. 30, pp. 123–155, 1983.
  * (36) H. Spohn, “The phonon boltzmann equation, properties and link to weakly anharmonic lattice dynamics,” Journal of Statistical Physics, vol. 124, pp. 1041–1104, 2006.
  * (37) T. Dudnikova and H. Spohn, “Local stationarity for lattice dynamics in the harmonic approximation,” arXiv preprint math-ph/0505031, 2005.
  * (38) C. Bernardin, F. Huveneers, and S. Olla, “Hydrodynamic limit for a disordered harmonic chain,” Communications in Mathematical Physics, vol. 365, pp. 215–237, 2019.
  * (39) E. Wigner, “On the quantum correction for thermodynamic equilibrium,” Physical Review, vol. 40, no. 5, p. 749, 1932.
  * (40) M. Hillery, R. F. O’Connell, M. O. Scully, and E. P. Wigner, Distribution Functions in Physics: Fundamentals, pp. 273–317.  Berlin, Heidelberg: Springer Berlin Heidelberg, 1997.
  * (41) M. Simoncelli, N. Marzari, and F. Mauri, “Wigner formulation of thermal transport in solids,” Physical Review X, vol. 12, no. 4, p. 041011, 2022.
  * (42) B. Bertini, M. Fagotti, L. Piroli, and P. Calabrese, “Entanglement evolution and generalized hydrodynamics: noninteracting systems,” Journal of Physics A: Mathematical and Theoretical, vol. 51, no. 39, p. 39LT01, 2018.
  * (43) M. Saha, M. Kulkarni, and A. Dhar, “Generalised Hydrodynamics description of the Page curve-like dynamics of a freely expanding fermionic gas,” arXiv preprint arXiv:2402.18422, 2024.
  * (44) H. Tasaki, “Macroscopic Irreversibility in Quantum Systems: ETH and Equilibration in a Free Fermion Chain,” arXiv preprint arXiv:2401.15263, 2024.
  * (45) S. De Bievre and P. E. Parris, “A rigourous demonstration of the validity of Boltzmann’s scenario for the spatial homogenization of a freely expanding gas and the equilibration of the Kac ring,” Journal of Statistical Physics, vol. 168, pp. 772–793, 2017.