# Navier-Stokes equations for nearly integrable quantum gases

Maciej Łebek<EMAIL_ADDRESS>Miłosz Panfil<EMAIL_ADDRESS>Faculty of Physics, University of Warsaw, Pasteura 5, 02-093 Warsaw, Poland

(August 28, 2024)

###### Abstract

The Navier-Stokes equations are paradigmatic equations describing
hydrodynamics of an interacting system with microscopic interactions encoded
in transport coefficients. In this work we show how the Navier-Stokes
equations arise from the microscopic dynamics of nearly integrable $1d$
quantum many-body systems. We build upon the recently developed hydrodynamics
of integrable models to study the effective Boltzmann equation with collision
integral taking into account the non-integrable interactions. We compute the
transport coefficients and find that the resulting Navier-Stokes equations
have two regimes, which differ in the viscous properties of the resulting
fluid. We illustrate the method by computing the transport coefficients for an
experimentally relevant case of coupled $1d$ cold-atomic gases.

Hydrodynamics is a universal theory describing phenomena of transport and
applicable to systems of sizes ranging from nuclear through biological up to
cosmological scales [1, 2, 3]. When viewed as an emergent theory, it assumes a
huge reduction of the degrees of freedom present at the microscopic scale. The
high-energy degrees of freedom are effectively integrated out which leads to
dissipation effects present in the Navier-Stokes (NS) equations and captured
by the transport coefficients. These are non-universal quantities depending on
the details of a microscopic theory. In the present work we show how the NS
equations arise from quantum-many body dynamics for systems that are nearly
integrable. We also derive expressions for the transport coefficients which
treat the integrable part of the interactions in an exact and non-perturbative
way.

Recent years have witnessed important developments in the field of non-
equilibrium dynamics of isolated quantum-many body systems. The program was
especially successful in $1+1$ dimensions when methods of quantum
integrability can be applied. This was fuelled by both experimental progress
in creating and manipulating cold-atomic systems [4, 5, 6, 7, 8, 9, 10, 11]
and by theoretical developments [12, 13, 14, 15, 16]. It resulted in an
understanding that typically quantum integrable systems relax to a generalized
Gibbs ensemble (GGE) [17, 18, 19] which includes, beside the energy, also the
additional conserved charges present in such theories.

The presence of additional conservation laws affects also the hydrodynamics of
integrable models. The resulting generalized hydrodynamics (GHD) [20, 21, 22,
14, 23] takes a form of a coupled continuity equations for an extensive number
of conserved charge densities. The dissipative processes are also taken into
account a as diffusion intrinsic to the integrable dynamics [24, 25, 26]. As
the hydrodynamic picture assumes that locally the system is in equilibrium, it
requires a separation of of time-scales between small $t_{\rm GGE}$ of the
local thermalization to the GGE and large $t_{\rm GHD}$ associated with the
flow created by spatial inhomogeneities.

Figure 1: top panel: Typically quantum integrable systems locally relax to the
GGE with some characteristic time-scale $t_{\rm GGE}$. When the system is
globally inhomogeneous further evolution is described by the GHD and occurs at
time-scale associated with dynamics of spatial gradients of conserved
densities. Finally, in the presence of integrability breaking terms and at
larger time-scales, the energy dissipates from higher conserved charges to the
$3$ elementary ones resulting in the NS equations. bottom panel: Our method
separates the full Hamiltonian into integrable Hamiltonian and weak
integrability breaking term $V_{2}$. The transport coefficients have then two
contributions, one from the interactions present in the integrable theory
(integrable collisions) and second from the integrability breaking term (non-
integrable collisions).

This formalism can be extended beyond purely integrable models by capturing
the additional interactions through the Boltzmann collision integral [27, 28,
29, 7, 8, 30, 31, 32, 33, 34, 35]. Building on that results one finds that
there is an additional time scale, that we call $t_{\rm NS}$, associated with
the integrability breaking terms. In this work we show that for models with
Galilean invariance and with the new interaction term which conserves only the
particle number, momentum and energy this results in the NS equations. We also
find that there are two contributions to the transport coefficients. They
origin from two types of collisions present in the system. First are
integrable collisions captured by the diffusion term of the GHD. Second are
the integrability breaking collisions captured by the Boltzmann collision
integral. This is summarized in fig. 1. Our approach to the problem
generalizes the classic kinetic theory derivations of hydrodynamics through
Chapman-Enskog expansion [36].

Generalized hydrodynamics: The GHD is a theory of hydrodynamics for integrable
models, such as integrable spin chains or 1d Bose gas, where an extensive
number of conservation laws are present. Whereas the theory is fairly new [20,
21] (with earlier developments in [37]), it has been studied quite extensively
(e.g. see the collections of reviews [22]) and was verified numerically [38,
24, 39] and experimentally [5, 40, 28]. Recently it was also shown how to
extend its applicability, by the logic of Boltzmann collision integral, to the
cases when the integrability is weakly broken [41, 27, 42, 28, 29, 43, 44, 45,
46, 47].

The GHD assumes a local relaxation of the system to the GGE state described by
the rapidity distribution $\rho_{\rm p}\equiv\rho_{\rm p}(\lambda;x,t)$ [48].
The GHD equation is then

$\partial_{t}\rho_{\rm p}+\partial_{x}\left(v_{\rho}\rho_{\rm
p}\right)+\mathfrak{f}\partial_{\lambda}\rho_{\rm
p}=\frac{1}{2}\partial_{x}\left(\mathfrak{D}_{\rho}\partial_{x}\rho_{\rm
p}\right)+\mathcal{I}[\rho_{\rm p}],$ (1)

with effective velocity $v_{\rho}$, the diffusion kernel $\mathfrak{D}_{\rho}$
and the Boltzmann collision integral $\mathcal{I}[\rho_{\rm p}]$. We also
consider external potential $U(x)$ and associated force term
$\mathfrak{f}=-\partial_{x}U$ [49, 50, 51, 52]. The index $\rho$ signifies
that the effective velocity and diffusion kernel are functionals of $\rho_{\rm
p}$ and are space-time dependent. This renders the GHD equations highly
nonlinear. The Boltzmann integral can be evaluated from the matrix elements of
the perturbing operator in the spirit of the Fermi’s golden rule [27, 29].

The expressions for the effective velocity and the diffusion kernel are quite
involved and as their explicit form for the present discussion is not
important we relegate the details to [53]. Similarly, for the Boltzmann
collision integral, the only important aspect are their collision invariants
which we now discuss.

In what follows we quite generally assume that there are $3$ collision
invariants corresponding to the total particle number, total momentum and
energy such that

$\int{\rm d}\lambda\,\mathscr{h}_{j}(\lambda)\,\mathcal{I}[\rho_{\rm
p}](\lambda)=0,\quad j=0,1,2.$ (2)

where $\mathscr{h}_{j}(\lambda)$ are single-particle eigenvalues of the
conserved charges of the underlying integrable model,

$\mathscr{q}_{j}=\langle\mathscr{h}_{j}\rangle\equiv\int{\rm
d}\lambda\,\mathscr{h}_{j}(\lambda)\rho_{\rm p}(\lambda).$ (3)

Functions $\mathscr{h}_{j}(\lambda)$ for a Galilean invariant theory are
$\mathscr{h}_{j}(\lambda)=\lambda^{j}/j!$. The $3$ collisions invariants
determine the nature of the stationary state of (1) to be a homogeneous
thermal boosted state.

Macroscopic conservation laws: The $3$ collision invariants lead to
macroscopic conservation laws. These are conveniently formulated using
hydrodynamic fields

$u(x,t)=\frac{\mathscr{q}_{1}(x,t)}{\mathscr{q}_{0}(x,t)},\quad
e(x,t)=\frac{\mathscr{q}_{2}(x,t)}{\mathscr{q}_{0}(x,t)}-\frac{\mathscr{q}_{1}(x,t)^{2}}{2\mathscr{q}_{0}(x,t)},$
(4)

and we additionally rename $\varrho(x,t)=\mathscr{q}_{0}(x,t)$. To this end,
we multiply the GHD eq. (1) by $\mathscr{h}_{j}(\lambda)$, integrate over
$\lambda$ and use the property (2). Following this procedure for $j=0,1,2$ we
obtain

$\displaystyle\partial_{t}\varrho$ $\displaystyle=-\partial_{x}(\varrho
u),\qquad\partial_{t}(\varrho u)=-\partial_{x}(\varrho
u^{2}+\mathcal{P})+\mathfrak{f}\varrho,$ (5)
$\displaystyle\partial_{t}(\varrho e)$ $\displaystyle=-\partial_{x}(u\varrho
e+\mathcal{J})-\mathcal{P}\partial_{x}u,$

with $\mathcal{P}(x,t)$ local pressure and $\mathcal{J}(x,t)$ local heat
current. The pressure and the heat current have two contributions. One
intrinsic to the diffusion term present in the GHD and given by
$\mathcal{P}_{\mathfrak{D}}=\mathscr{h}_{1}\mathfrak{D}_{\rho}\partial_{x}\rho_{\rm
p}/2$ and
$\mathcal{J}_{\mathfrak{D}}=\mathscr{h}_{2}\mathfrak{D}_{\rho}\partial_{x}\rho_{\rm
p}/2$. The second one resembles the standard expressions from the kinetic
theory,

$\displaystyle\mathcal{P}_{v}=\langle(\lambda-u)v(\lambda)\rangle,\quad\mathcal{J}_{v}=\frac{1}{2}\langle(\lambda-u)^{2}(v(\lambda)-u)\rangle.$
(6)

The derivation of the conservation laws together with the formulas for
$\mathcal{P}$ and $\mathcal{J}$ is elementary and presented in [53]. We only
note that in deriving the first conservation law the identity $\langle
v\rangle=\langle\mathscr{h}_{1}\rangle$, a consequence of the Galiliean
invariance, is operatorial. This property establishes the velocity field $u$ a
relevant hydrodynamic degree of freedom.

The equations (5) are not closed because the pressure and heat currents are
functions of the full distribution $\rho_{\rm p}$ rather than just its $3$
moments encoded through fields $\varrho$, $u$ and $e$. To circumvent this
problem, in the phenomenological treatment one assumes the following
expressions

$\displaystyle\mathcal{P}=P-\zeta\partial_{x}u,\qquad\mathcal{J}=-\kappa\partial_{x}T,$
(7)

in terms of the hydrodynamic fields. Here $P=P(\varrho,e)$ is the hydrostatic
pressure, $\zeta$ is the bulk viscosity and $\kappa$ is the thermal
conductivity, respectively. The temperature field and the energy field are
related through the thermodynamic identity $T^{-1}=(\partial s/\partial
e)_{\varrho}$. This leaves the theory with two phenomenological constants
$\zeta$ and $\kappa$ called transport coefficients. In this way we obtain the
one-dimensional version of NS equations [54, 55, 56], for the fields
$\varrho(x,t)$, $u(x,t)$ and $e(x,t)$. In the following we confirm this
phenomenological picture and provide explicit formulas for the transport
coefficients.

Chapman-Enskog method: The standard procedure to justify the phenomenological
picture and to turn the hydrodynamic equations into a closed set is the
Chapman-Enskog (ChE) method. In the process it also provides a systematic
computation of the transport coefficients from the microscopic theory.
Originally it was formulated in the context of classical kinetic theory [36,
1]. Here we adapt it to the GHD case. The ChE method assumes the time
evolution of $\rho_{\rm p}$ is determined by the hydrodynamic fields, namely
$\rho_{\rm p}(\lambda,x|\varrho_{\alpha}(x,t))$, called normal solutions. With
this assumption the GHD equations (1) can be rewritten as

$\sum_{\alpha}\frac{\partial\rho_{\rm
p}}{\partial\varrho_{\alpha}}\partial_{t}\varrho_{\alpha}+\partial_{x}\left(v_{\rho}\rho_{\rm
p}\right)=\frac{1}{2}\partial_{x}\left(\mathfrak{D}_{\rho}\partial_{x}\rho_{\rm
p}\right)+\mathcal{I}[\rho_{\rm p}],$ (8)

with $\\{\varrho_{\alpha}\\}={\varrho,u,e}$. As the external potential does
not influence the determination of the transport coefficients we set it here
to zero. In what follows it will be useful to have a compact notation for the
conservation laws. We write
$\partial_{t}\varrho_{\alpha}=D_{\alpha}[\\{\varrho_{\beta}\\}]$ where
$D_{\alpha}[\\{\varrho_{\beta}\\}]$ can be read out from eq. (5). With this we
find the ChE equation

$\sum_{\alpha}\frac{\partial\rho_{\rm
p}}{\partial\varrho_{\alpha}}D_{\alpha}+\partial_{x}\left(v_{\rho}\rho_{\rm
p}\right)=\frac{\delta_{1}}{2}\partial_{x}\left(\mathfrak{D}_{\rho}\partial_{x}\rho_{\rm
p}\right)+\frac{1}{\delta_{2}}\mathcal{I}[\rho_{\rm p}],$ (9)

which is an equation for $\rho_{\rm p}$ given the $\\{\varrho_{\alpha}\\}$ and
contains only spatial derivatives. The effective velocity is now determined
fully from the knowledge of $\\{\varrho_{\alpha}\\}$. In writing eq. (9) we
have introduced two small parameters $\delta_{1}$ and $\delta_{2}$. Their
inverses define typical distances (in units of a microscopic length) at which
effects of GHD diffusion and integrability breaking occur respectively.
Analyzing the problem order by order $\delta_{i}$, we find that the leading
order requires that $\mathcal{I}[\rho_{\rm p}^{(0)}]=0$. This is the Euler
scale hydrodynamics. The solution is the locally varying boosted thermal
state. Specifically, explicit computations show that the hydrodynamic pressure
$\mathcal{P}$ reduces then to the hydrostatic pressure $P(\varrho,e)$ [57, 53,
15], and both the bulk viscosity and the heat current vanish identically. The
hydrodynamic equations at the Euler scale are

$\displaystyle\partial_{t}\varrho$ $\displaystyle=-\partial_{x}(\varrho
u),\qquad\partial_{t}(\varrho u)=-\partial_{x}(\varrho u^{2}+P),$ (10)
$\displaystyle\partial_{t}(\varrho e)$ $\displaystyle=-\partial_{x}(u\varrho
e)-P\partial_{x}u,$

and are fully determined from the thermodynamics of the system encoded in the
equation of state $P=P(\varrho,e)$. We thus partially confirm the
phenomenological picture of eq. (7) in which $\mathcal{P}$ and $\mathcal{J}$,
beside the hydrostatic pressure $P$, contain derivative terms and do not enter
at the Euler scale. This renders the Euler scale dynamic reversible as can be
also witnessed by the conservation of entropy. Indeed, simple computations
show that the thermodynamic entropy density $\mathsf{s}(x,t)$ [57, 48, 58]
obeys the continuity equation
$\partial_{t}\mathsf{s}+\partial_{x}(v\mathsf{s})=0$ [53].

Expanding the ChE Eq. (9) to the first orders in both parameters
$\delta_{1},\delta_{2}$ allows us to find the transport coefficients. As
announced earlier there are two contributions. The one from the GHD diffusion
amounts to evaluating $\mathcal{P}_{\mathfrak{D}}$ and
$\mathcal{J}_{\mathfrak{D}}$ on a boosted thermal state. The second
contribution is due to the collision integral and is more technically
involved. However both contributions have a feature known from the standard
kinetic theory [36, 59]. Namely, that the transport the coefficients are a
simple generalization of transport coefficients found in a linearized
hydrodynamics. Therefore, in the present work we discuss the linearized
theory, which also has the advantage of revealing the physical ingredients and
defer the further computations within the ChE method to a longer publication
[60].

The transport coefficients that we find are
$\zeta=\zeta_{\mathcal{I}}+\zeta_{\mathfrak{D}}$ and
$\kappa=\kappa_{\mathcal{I}}+\kappa_{\mathfrak{D}}$ with

$\displaystyle\zeta_{\mathcal{I}}$
$\displaystyle=\mathscr{h}_{1}\mathcal{B}\phi_{\zeta},\quad\zeta_{\mathfrak{D}}=\beta\mathscr{h}_{1}\mathfrak{D}\mathcal{C}\mathscr{h}_{1}/2,$
(11) $\displaystyle\kappa_{\mathcal{I}}$
$\displaystyle=\mathscr{h}_{2}\mathcal{B}\phi_{\kappa},\quad\kappa_{\mathfrak{D}}=\beta^{2}\mathscr{h}_{2}\mathfrak{D}\mathcal{C}\mathscr{h}_{2}/2,$

where $\mathcal{B}$ and $\mathcal{C}$ are current and charge susceptibility
matrices of the integrable model [61, 62]. Additionally, $\beta\equiv 1/T$ and
functions $\phi_{\kappa,\zeta}(\lambda)$ are found from the following integral
equations

$\Gamma\phi_{\kappa,\zeta}=\chi_{\kappa,\zeta},$ (12)

where $\Gamma=-(\delta\mathcal{I}/\delta\rho_{\rm p})\mathcal{C}$ evaluated in
a thermal state, and

$\displaystyle\chi_{\zeta}=\sqrt{\beta\varrho}\,\mathcal{B}h_{1}-\sqrt{\beta/\varkappa_{T}}\,\mathcal{C}(h_{0}+\sqrt{c_{P}/c_{V}-1}\,h_{2}),$
(13) $\displaystyle\chi_{\kappa}=\beta\sqrt{\varrho
c_{V}}\,\mathcal{B}(h_{2}-\sqrt{c_{P}/c_{V}-1}\,h_{0}).$

Isothermal compressibility $\varkappa_{T}$ and specific heats $c_{P},c_{V}$
are those of the unperturbed, integrable model. As we show in [53],
thermodynamic quantities are represented as explicit functionals of the
underlying thermal state with density $\varrho$ and temperature $T$. Functions
$h_{n}(\lambda)$ are charges constructed with hydrodynamic scalar product
$(\cdot|\cdot)$, discussed later in the text. The integral equations (12) are
supplemented with additional conditions
$(\mathscr{h}_{i}|\phi_{\kappa,\zeta})=0,\,\,i=0,1,2$, which render their
solutions unique. Equations (12) naturally generalize ChE integral equations
[36] to interacting, nearly integrable models. In the process we confirm that
there are no other gradient terms in $\mathcal{P}$ and $\mathcal{J}$ beside
the two given in the phenomenological expression (7). We will provide now a
justification of these results from the linearized hydrodynamics.

Linearized hydrodynamics: We linearize the GHD equations (1) around the
homogeneous thermal state such that $\rho_{\rm p}(x,t,\lambda)=\rho_{\rm
th}(\lambda)+\delta\rho_{\rm p}(x,t,\lambda)$,

$\partial_{t}\delta\rho_{\rm p}+\partial_{x}\left(v_{\rm th}\delta\rho_{\rm
p}\right)=\frac{1}{2}\mathfrak{D}_{\rm th}\partial_{x}^{2}\delta\rho_{\rm
p}+\frac{\delta\mathcal{I}}{\delta\rho_{\rm p}}\delta\rho_{\rm p}.$ (14)

The linearized GHD equation (14) leads to an equation for perturbations of
conserved quantities. The structure of the resulting equations is simpler if
we work in the orthonormal basis of the conserved charges with respect to the
hydrodynamic scalar product $(g|h)\equiv g\mathcal{C}h$. We denote the
resulting charges by $q_{n}$ and the corresponding single-particle eigenvalues
$h_{n}(\lambda)$, c.f. (3).

Standard computations using the methods of Thermodynamic Bethe Ansatz yield,
see [53],

$\partial_{t}\delta q_{n}+A_{nm}\partial_{x}\delta
q_{m}-D_{nm}\partial_{x}^{2}\delta q_{m}=-\Gamma_{nm}\delta q_{m},$ (15)

where we assume a summation over the repeated indices and $\delta
q_{n}(x,t)=\int{\rm d}\lambda\,h_{n}(\lambda)\delta\rho_{\rm p}(\lambda;x,t)$.
Moreover, $A_{nm}$, $D_{nm}$ and $\Gamma_{nm}$ are now interpreted as matrix
elements of operators associated with velocity, diffusion and collision
integral evaluated in the orthonormal basis. We observe that $D_{0n}=D_{n0}$
which is a consequence of Markovianity of the GHD diffusion operator [24, 53],
while $\Gamma_{jn}=0$ for $j=0,1,2$ which reflects the collision invariants
(2).

The linearized GHD equations (15) can be now solved perturbatively in
gradients of $\delta q_{n}$. This is the easiest to formulate by going to the
Fourier space and considering a time evolution of a single mode $\delta
q_{n}(x,t)\sim\exp(ikx-\Lambda(k)t)\delta q_{n}(k)$. This leads to the
eigenvalue problem

$\left(\Gamma_{nm}+ikA_{nm}+k^{2}D_{nm}\right)\delta q_{m}(k)=\Lambda(k)\delta
q_{n}(k),$ (16)

with $\Lambda(k)$ the dispersion relations for different modes. Among them
there are $3$ gapless modes related to the $3$ conservation laws respected by
the collision integral. For small $k$ a standard perturbation theory methods
apply and we find two sound modes and a heat mode (thermal mode) with
dispersions (we neglect $\mathcal{O}(k^{3})$ terms)

$\displaystyle\Lambda_{\pm}(k)$ $\displaystyle=\pm
iv_{s}k+\frac{1}{2\varrho}\bigg{[}\left(\frac{1}{c_{V}}-\frac{1}{c_{P}}\right)\kappa+\zeta\bigg{]}k^{2},$
(17) $\displaystyle\Lambda_{\rm th}(k)$ $\displaystyle=\frac{\kappa}{\varrho
c_{P}}k^{2},$

which agree with dispersion relations of linearized NS equations. The sound
velocity is $v_{s}=\left(A_{1,0}^{2}+A_{2,1}^{2}\right)^{1/2}$ and is equal to
thermodynamic formula $v_{s}=1/\sqrt{\varrho\varkappa_{S}}$, where
$\varkappa_{S}$ is the adiabatic compressibility of the integrable model
calculated for the thermal background state $\rho_{\rm th}$. The transport
coefficients $\zeta$ and $\kappa$ are given by (11) computed with the same
state $\rho_{\rm th}$. The details are referred to [53].

The corresponding eigenmodes are related to the hydrodynamic fields $\varrho$,
$u$ and $T$

$\displaystyle\delta\varrho=\varrho\sqrt{T\varkappa_{T}}\,\delta
q_{0},\,\,\,\,\delta u=\sqrt{T/\varrho}\,\delta q_{1},\,\,\,\,\delta
T=T/\sqrt{\varrho c_{V}}\,\delta q_{2},$ (18)

where we took $\varrho(x,t)=\varrho+\delta\varrho(x,t)$, $u(x,t)=0+\delta
u(x,t)$ and $T(x,t)=T+\delta T(x,t)$. Therefore the first $3$ orthonormalized
charges are equivalent to the hydrodynamic degrees of freedom. Expressing heat
and sound modes in terms of $\delta\varrho,\delta u,\delta T$ we recover
precisely the modes of NS equations [53].

Application to the Lieb-Liniger model: We now apply the introduced formalism
to a situation relevant to experiments with cold atomic gases. We consider a
setup of two 1d Bose gases, described by the integrable Lieb-Liniger (LL)
model [63], coupled by a long-range interaction which breaks the
integrability. The integrable part of the Hamiltonian is $H_{\rm
int}=H_{1}+H_{2}$ with

$H_{i}=\int{\rm
d}x\,\psi_{i}^{\dagger}(x)\left(-\frac{1}{2}\partial_{x}^{2}+\frac{c}{2}\,\psi_{i}^{\dagger}(x)\psi_{i}(x)\right)\psi_{i}(x),$
(19)

where we work in the units $\hbar=m=1$ and $\psi_{i}(x),\psi_{i}^{\dagger}(x)$
are canonical bosonic fields. The integrability breaking coupling between the
tubes is given by $V_{2}=\int{\rm d}x{\rm
d}yV(x-y)\hat{\rho}_{1}(x)\hat{\rho}_{2}(y)$, with
$\hat{\rho}_{i}(x)=\psi_{i}^{\dagger}(x)\psi_{i}(x)$ and $V(x-y)$ interaction
potential. In an experimental realizations with cold-atomic gases it
originates from dipole-dipole interactions and its strength is tunable [6,
10].

Of special interest is the case when the two tubes are initially in the same
state. Then they stay the same throughout the dynamics and the system can be
described by a single equation (1). The collision integral based on form-
factor approach [64, 65, 66] and the resulting relaxation for the homogeneous
initial state was studied in [29]. For the computation of the transport
coefficients we need $\mathcal{I}$ linearized around a thermal state,
c.f.(12). The derivation is presented in [53]. The result for the matrix
elements of $\Gamma$ in the basis $\\{h_{j}\\}$ is

$\displaystyle\Gamma_{nm}\\!=\frac{1}{\tau}\\!\int{\rm
d}\boldsymbol{\lambda}_{1}{\rm
d}\boldsymbol{\lambda}_{2}\,g_{nm}(\lambda_{1},\lambda_{2})\left(\mathcal{T}^{\rm
dr}(\lambda_{1},\lambda_{2})\right)^{2}|v_{12}|,$ (20)

with the integration measure ${\rm d}\boldsymbol{\lambda}=\rho_{\rm
p}(\lambda)(1-n(\lambda)){\rm d}\lambda$ with $n(\lambda)$ the filling
function, $v_{i}=v(\lambda_{i})$, $v_{12}=v_{1}-v_{2}$, and

$g_{nm}(\lambda_{1},\lambda_{2})=\int{\rm
d}\boldsymbol{\lambda}\,(k^{\prime}(\lambda))^{2}\frac{\mathcal{S}[h_{n}^{\rm
Dr}]\mathcal{S}[h_{m}^{\rm dr}]}{(v-v_{1})^{2}(v-v_{2})^{2}}.$ (21)

Here $\mathcal{T}$ is the differential scattering kernel, for the LL model
$\pi\mathcal{T}(\lambda-\lambda^{\prime})=c/(c^{2}+(\lambda-\lambda^{\prime})^{2})$
and ${\rm dr}$ and ${\rm Dr}$ are dressing operations [67], $k(\lambda)$ is
the Dressed momentum while $S[\varphi]\equiv
S[\varphi](\lambda,\lambda_{1},\lambda_{2})$ is

$S[\varphi]=\frac{\varphi^{\prime}}{k^{\prime}(v_{1}-v_{2})}+\frac{\varphi_{1}^{\prime}}{k_{1}^{\prime}(v_{2}-v)}+\frac{\varphi_{2}^{\prime}}{k_{2}^{\prime}(v-v_{1})}.$
(22)

The timescale $\tau$ is related to the strength of inter-tube coupling
$\tau^{-1}=\tilde{V}^{2}m/(4\hbar^{2})$, where $\tilde{V}^{2}=\int
dk\tilde{V}(k)k^{2}$ and $\tilde{V}(k)$ is Fourier transform of $V(x-y)$.

The transport coefficients are then determined from (11) and the numerical
solution is shown in Fig. 2. We focus on the strong coupling regime of the LL
model in the nomenclature of [68] where both interaction and thermal
fluctuations are important. We find that there are two competing behaviors.
When intratube interactions dominate the resulting fluid has relatively small
viscosity and large heat current. Instead, when intertube interactions
dominate then both viscosity and heat current are of the same order. In the
limit of noninteracting integrable model $c,c^{-1}\to 0$ the bulk viscosity
vanishes and thermal conductivity diverges [53] as expected from the classical
kinetic theory of weakly interacting models [54].

Figure 2: Transport coefficients in the log-log scale for the coupled LL
models as a function of intratube interaction coupling $c$. In the main plot
we show separate contributions to the heat current and viscosity from
integrable and non-integrable collisions. In the inset we show the total heat
current and viscosity for $\tau=10^{-2}$. The viscosity is non-zero across
different values of temperature and couplings.

Outlook: In this work we have shown that the NS equations universally arise
for $1d$ quantum systems with weakly broken integrability and with Galilean
invariance, assuming the GHD equation supplemented with the collision
integral. We have provided an explicit equations (12) for the transport
coefficients linking the effective hydrodynamic description with the exact
microscopic thermodynamics of the integrable model. Unlike the standard
approach in the kinetic theory, we treat an important part of the interactions
non-perturbatively. One consequence of this is a non-zero viscosity, which in
the derivation of the NS equations based on the free theory, is identically
zero [36, 55]. The resulting fluid has two regimes. One dominated by the
integrability breaking collisions with relatively small viscosity, and the
second, dominated by the integrable collisions with heat current and viscosity
of the same order.

It is well-known that transport in one dimension with conserved momentum
becomes anomalous upon introducing noise to the system [69, 70, 71, 72]. The
linearized framework with $U(x)=0$ presented here together with initial
results for GHD with noise [35] open possibility for further explorations in
this direction.

Recently, the NS equations were also found for the low temperature dynamics of
the GHD without an explicit integrability breaking terms [73]. The effective
projection to the three lowest conserved densities is attributed then to the
low temperatures and holds at intermediate timescales. The transport
coefficients are then of the form $\zeta_{\mathfrak{D}}$ and
$\kappa_{\mathfrak{D}}$ taken in the $T\to 0$ limit, namely they origin from
the integrable collisions. Our results apply also to classical integrable
Galilean-invariant models, whose dynamics is described by (1). In such cases
collision integrals may follow from BBGKY equations or be constructed
phenomenologically, as in the relaxation time approximation [30].

The authors thank Romain Vasseur for his inspiration and encouragement to
tackle this problem. We also thank Alvise Bastianello and Sarang
Gopalakrishnan for discussions. The authors acknowledge support by the
National Science Centre (NCN), Poland via projects 2018/31/D/ST3/03588 and
2022/47/B/ST2/03334.

## References

  * Dorfman  _et al._ [2021] J. R. Dorfman, H. van Beijeren, and T. R. Kirkpatrick,  _Contemporary Kinetic Theory of Matter_ (Cambridge University Press, 2021).
  * Rezzolla and Zanotti [2013] L. Rezzolla and O. Zanotti,  _Relativistic Hydrodynamics_ (Oxford University Press, 2013).
  * Vogel [1996] S. Vogel,  _Life in Moving Fluids: The Physical Biology of Flow_ , Princeton paperbacks (Princeton University Press, 1996).
  * Kinoshita  _et al._ [2006] T. Kinoshita, T. Wenger, and D. S. Weiss, Nature 440, 900 (2006).
  * Schemmer  _et al._ [2019] M. Schemmer, I. Bouchoule, B. Doyon, and J. Dubail, Phys. Rev. Lett. 122, 090601 (2019), arXiv:1810.07170 [cond-mat.quant-gas] .
  * Tang  _et al._ [2018] Y. Tang, W. Kao, K.-Y. Li, S. Seo, K. Mallayya, M. Rigol, S. Gopalakrishnan, and B. L. Lev, Physical Review X 8, 021030 (2018).
  * Cataldini  _et al._ [2022] F. Cataldini, F. Møller, M. Tajik, J. Sabino, S.-C. Ji, I. Mazets, T. Schweigler, B. Rauer, and J. Schmiedmayer, Physical Review X 12, 041032 (2022).
  * Møller  _et al._ [2021] F. Møller, C. Li, I. Mazets, H.-P. Stimming, T. Zhou, Z. Zhu, X. Chen, and J. Schmiedmayer, Physical Review Letters 126, 090602 (2021).
  * Li  _et al._ [2020] C. Li, T. Zhou, I. Mazets, H.-P. Stimming, F. S. Møller, Z. Zhu, Y. Zhai, W. Xiong, X. Zhou, X. Chen, and J. Schmiedmayer, SciPost Physics 9, 058 (2020).
  * Li  _et al._ [2023] K.-Y. Li, Y. Zhang, K. Yang, K.-Y. Lin, S. Gopalakrishnan, M. Rigol, and B. L. Lev, Physical Review A 107, L061302 (2023).
  * Wilson  _et al._ [2020] J. M. Wilson, N. Malvania, Y. Le, Y. Zhang, M. Rigol, and D. S. Weiss, Science 367, 1461 (2020).
  * D’Alessio  _et al._ [2016] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, Advances in Physics 65, 239 (2016).
  * Mori  _et al._ [2018] T. Mori, T. N. Ikeda, E. Kaminishi, and M. Ueda, Journal of Physics B: Atomic, Molecular and Optical Physics 51, 112001 (2018).
  * Doyon  _et al._ [2023] B. Doyon, S. Gopalakrishnan, F. Møller, J. Schmiedmayer, and R. Vasseur, (2023), arXiv:2311.03438 [cond-mat.stat-mech] .
  * Bouchoule and Dubail [2022] I. Bouchoule and J. Dubail, Journal of Statistical Mechanics: Theory and Experiment 2022, 014003 (2022).
  * Doyon  _et al._ [2017] B. Doyon, J. Dubail, R. Konik, and T. Yoshimura, Physical Review Letters 119, 195301 (2017).
  * Rigol  _et al._ [2008] M. Rigol, V. Dunjko, and M. Olshanii, Nature 452, 854 (2008).
  * Ilievski  _et al._ [2015] E. Ilievski, J. De Nardis, B. Wouters, J.-S. Caux, F. H. L. Essler, and T. Prosen, Phys. Rev. Lett. 115, 157201 (2015).
  * Langen  _et al._ [2015] T. Langen, S. Erne, R. Geiger, B. Rauer, T. Schweigler, M. Kuhnert, W. Rohringer, I. E. Mazets, T. Gasenzer, and J. Schmiedmayer, Science 348, 207 (2015).
  * Bertini  _et al._ [2016a] B. Bertini, M. Collura, J. De Nardis, and M. Fagotti, Phys. Rev. Lett. 117, 207201 (2016a).
  * Castro-Alvaredo  _et al._ [2016] O. A. Castro-Alvaredo, B. Doyon, and T. Yoshimura, Phys. Rev. X 6, 041065 (2016).
  * Bastianello  _et al._ [2022] A. Bastianello, B. Bertini, B. Doyon, and R. Vasseur, Journal of Statistical Mechanics: Theory and Experiment 2022, 014001 (2022).
  * Doyon [2020] B. Doyon, SciPost Phys. Lect. Notes , 18 (2020).
  * De Nardis  _et al._ [2018] J. De Nardis, D. Bernard, and B. Doyon, Phys. Rev. Lett. 121, 160603 (2018).
  * Gopalakrishnan  _et al._ [2018] S. Gopalakrishnan, D. A. Huse, V. Khemani, and R. Vasseur, Phys. Rev. B 98, 220303 (2018).
  * Nardis  _et al._ [2019] J. D. Nardis, D. Bernard, and B. Doyon, SciPost Physics 6, 049 (2019).
  * Durnin  _et al._ [2021a] J. Durnin, M. J. Bhaseen, and B. Doyon, Phys. Rev. Lett. 127, 130601 (2021a).
  * Møller  _et al._ [2021] F. Møller, C. Li, I. Mazets, H.-P. Stimming, T. Zhou, Z. Zhu, X. Chen, and J. Schmiedmayer, Phys. Rev. Lett. 126, 090602 (2021).
  * Panfil  _et al._ [2023] M. Panfil, S. Gopalakrishnan, and R. M. Konik, Phys. Rev. Lett. 130, 030401 (2023).
  * Lopez-Piqueres  _et al._ [2021] J. Lopez-Piqueres, B. Ware, S. Gopalakrishnan, and R. Vasseur, Physical Review B 103, L060302 (2021).
  * Lopez-Piqueres  _et al._ [2022] J. Lopez-Piqueres, S. Gopalakrishnan, and R. Vasseur, Journal of Physics A: Mathematical and Theoretical 55, 234005 (2022).
  * Lopez-Piqueres and Vasseur [2023] J. Lopez-Piqueres and R. Vasseur, Physical Review Letters 130, 247101 (2023).
  * Caux  _et al._ [2019] J.-S. Caux, B. Doyon, J. Dubail, R. Konik, and T. Yoshimura, SciPost Physics 6, 070 (2019).
  * Bouchoule  _et al._ [2020] I. Bouchoule, B. Doyon, and J. Dubail, SciPost Physics 9, 044 (2020).
  * Bastianello  _et al._ [2020a] A. Bastianello, J. D. Nardis, and A. D. Luca, Physical Review B 102, 161110 (2020a).
  * Résibois and De Leener [1977] P. Résibois and M. De Leener,  _Classical Kinetic Theory of Fluids_ (John Wiley & Sons, 1977).
  * Bonnes  _et al._ [2014] L. Bonnes, F. H. L. Essler, and A. M. Läuchli, Phys. Rev. Lett. 113, 187203 (2014).
  * Doyon and Spohn [2017a] B. Doyon and H. Spohn, Journal of Statistical Mechanics: Theory and Experiment 2017, 073210 (2017a).
  * Biagetti  _et al._ [2023] L. Biagetti, G. Cecile, and J. D. Nardis, (2023), arXiv:2307.05379 [cond-mat.stat-mech] .
  * Malvania  _et al._ [2021] N. Malvania, Y. Zhang, Y. Le, J. Dubail, M. Rigol, and D. S. Weiss, Science 373, 1129 (2021).
  * Friedman  _et al._ [2020] A. J. Friedman, S. Gopalakrishnan, and R. Vasseur, Phys. Rev. B 101, 180302 (2020).
  * Bastianello  _et al._ [2021] A. Bastianello, A. D. Luca, and R. Vasseur, Journal of Statistical Mechanics: Theory and Experiment 2021, 114003 (2021).
  * Pawłowski  _et al._ [2024] J. Pawłowski, M. Panfil, J. Herbrych, and M. Mierzejewski, Physical Review B 109, L161109 (2024).
  * Łebek  _et al._ [2024] M. Łebek, M. Panfil, and R. M. Konik, (2024), arXiv:2303.12490 [cond-mat.quant-gas] .
  * Bertini  _et al._ [2015] B. Bertini, F. H. Essler, S. Groha, and N. J. Robinson, Physical Review Letters 115, 180601 (2015).
  * Bertini  _et al._ [2016b] B. Bertini, F. H. L. Essler, S. Groha, and N. J. Robinson, Physical Review B 94, 245117 (2016b).
  * del Vecchio del Vecchio  _et al._ [2022] G. del Vecchio del Vecchio, A. D. Luca, and A. Bastianello, SciPost Physics 12, 060 (2022).
  * Mossel and Caux [2012] J. Mossel and J.-S. Caux, Journal of Physics A: Mathematical and Theoretical 45, 255001 (2012).
  * Doyon and Yoshimura [2017] B. Doyon and T. Yoshimura, SciPost Physics 2, 014 (2017).
  * Bastianello  _et al._ [2020b] A. Bastianello, A. D. Luca, B. Doyon, and J. D. Nardis, Physical Review Letters 125, 240604 (2020b).
  * Durnin  _et al._ [2021b] J. Durnin, A. D. Luca, J. D. Nardis, and B. Doyon, Journal of Physics A: Mathematical and Theoretical 54, 494001 (2021b).
  * Cao  _et al._ [2018] X. Cao, V. B. Bulchandani, and J. E. Moore, Physical Review Letters 120, 164101 (2018).
  * [53] See the Supplementary Material.
  * Balescu [1975] R. Balescu,  _Equilibrium and nonequilibrium statistical mechanics_ (John Wiley & Sons, 1975).
  * Spohn [2012] H. Spohn,  _Large Scale Dynamics of Interacting Particles_ , Theoretical and Mathematical Physics (Springer Berlin Heidelberg, 2012).
  * Chakraborti  _et al._ [2021] S. Chakraborti, S. Ganapa, P. Krapivsky, and A. Dhar, Physical Review Letters 126, 244503 (2021).
  * Yang and Yang [1969] C. N. Yang and C. P. Yang, Journal of Mathematical Physics 10, 1115 (1969).
  * Takahashi [1999] M. Takahashi,  _Thermodynamics of One-Dimensional Solvable Models_ (Cambridge University Press, 1999).
  * Resibois [1970] P. Resibois, Journal of Statistical Physics 2, 21 (1970).
  * [60] M. Łebek and M. Panfil, Chapman-Enskog theory for nearly integrable quantum gases, in preparation .
  * Doyon and Spohn [2017b] B. Doyon and H. Spohn, SciPost Phys. 3, 039 (2017b).
  * De Nardis  _et al._ [2022] J. De Nardis, B. Doyon, M. Medenjak, and M. Panfil, Journal of Statistical Mechanics: Theory and Experiment 2022, 014002 (2022).
  * Lieb and Liniger [1963] E. H. Lieb and W. Liniger, Physical Review 130, 1605 (1963).
  * Nardis and Panfil [2015] J. D. Nardis and M. Panfil, Journal of Statistical Mechanics: Theory and Experiment 2015, P02019 (2015).
  * Nardis and Panfil [2018] J. D. Nardis and M. Panfil, Journal of Statistical Mechanics: Theory and Experiment 2018, 033102 (2018).
  * Panfil [2021] M. Panfil, Journal of Statistical Mechanics: Theory and Experiment 2021, 013108 (2021).
  * Ruggiero  _et al._ [2020] P. Ruggiero, P. Calabrese, B. Doyon, and J. Dubail, Phys. Rev. Lett. 124, 140603 (2020).
  * Kheruntsyan  _et al._ [2003] K. V. Kheruntsyan, D. M. Gangardt, P. D. Drummond, and G. V. Shlyapnikov, Phys. Rev. Lett. 91, 040403 (2003).
  * van Beijeren [2012] H. van Beijeren, Physical Review Letters 108, 180601 (2012).
  * Narayan and Ramaswamy [2002] O. Narayan and S. Ramaswamy, Physical Review Letters 89, 200601 (2002).
  * Forster  _et al._ [1977] D. Forster, D. R. Nelson, and M. J. Stephen, Physical Review A 16, 732 (1977).
  * Spohn [2014] H. Spohn, Journal of Statistical Physics 154, 1191 (2014).
  * Urichuk  _et al._ [2023] A. Urichuk, S. Scopa, and J. D. Nardis, (2023), arXiv:2309.14476 [cond-mat.str-el] .

Supplemental Material:
Navier-Stokes equations for nearly integrable quantum gases
Maciej Łebek and Miłosz Panfil

Faculty of Physics, University of Warsaw, Pasteura 5, 02-093 Warsaw, Poland

In the Supplemental Material we provide computations supporting the results
presented in the main text.

In Sec. I we introduce the basic ingredients of the GHD equation: effective
velocity, diffusion kernel and susceptibility matrices along with properties
of these objects relevant in the present context. In Sec. II we derive
hydrodynamic conservation laws for the three conserved charges and introduce
dynamic pressure and heat current. Next, in Sec. III we analyze the basic
properties of thermal boosted state of integrable models. We show that for
such states the heat current vanishes and dynamic pressure becomes the
hydrostatic pressure. In Sec. IV we derive expressions for thermodynamic
quantities such as sound velocity and specific heats in terms of hydrodynamic
matrices.

The following sections concern the linearized hydrodynamics and computation of
transport coefficients. First, in Sec. V we present a derivation of a
linearized transport equation for local charge densities. In Sec. VI we
linearize Navier-Stokes equations and find the dispersion relations of heat
and sound modes. We use these results in Sec. VII where we treat the
eigenproblem perturbatively and identify expressions for bulk viscosity and
thermal conductivity. Additionally, we give numerical results for the
transport coefficients. Later in Sec. VIII we show that the results found in
the perturbation theory can be formulated as generalized Chapman-Enskog
integral equations, which are stated in the main text.

Finally, in Sec. IX we analyze the collision integral for two Lieb-Liniger
gases coupled by the density-density interactions.

###### Contents

  1. I Effective velocity, diffusion operator and susceptibility matrices
  2. II Hydrodynamic conservation laws
  3. III Thermal boosted state
  4. IV Thermodynamics and hydrodynamic matrices
    1. IV.1 Sound velocity
    2. IV.2 Specific heats, isothermal compressibility and additional identities
  5. V Derivation of the linearized equation for charge perturbations
  6. VI Linearized Navier-Stokes equations and transport coefficients
  7. VII Perturbation theory
    1. VII.1 Perturbation theory at order $k$
    2. VII.2 Perturbation theory at order $k^{2}$
    3. VII.3 Transport coefficients
  8. VIII Integral equations for $\kappa^{\rm coll}$ and $\zeta^{\rm coll}$ transport coefficients
  9. IX Collision integral for the coupled Lieb-Liniger models

##  I Effective velocity, diffusion operator and susceptibility matrices

We consider a Galilean-invariant model with energy
$E(\lambda)=\frac{1}{2}\lambda^{2}$ and momentum $p(\lambda)=\lambda$. We set
the mass to unity. As in the main text we introduce functions
$\mathscr{h}_{j}(\lambda)=\lambda^{j}/j!$ related to ultra-local charges. The
effective velocity [23] is given by $v(\lambda)=(E^{\prime})^{\rm
dr}(\lambda)/(p^{\prime})^{\rm dr}(\lambda)$ where the dressing operation of
an arbitrary function $g(\lambda)$ is implemented by the following integral
equation

$g^{\rm dr}(\lambda)=g(\lambda)+\int{\rm
d}\lambda^{\prime}\mathcal{T}(\lambda-\lambda^{\prime})n(\lambda^{\prime})g^{\rm
dr}(\lambda^{\prime}),$ (I.1)

and $\mathcal{T}(\lambda-\lambda^{\prime})$ is a symmetric scattering kernel,
which depends on the model. For instance, in the case of Lieb-Liniger model
[63] it reads $\mathcal{T}(\lambda)=\frac{c}{\pi}\frac{1}{c^{2}+\lambda^{2}}$,
where $c$ is the coupling constant. We also defined the occupation function
$n(\lambda)=\rho_{\rm p}(\lambda)/\rho_{\rm tot}(\lambda)$, where total
density of states $\rho_{\rm tot}(\lambda)$ is

$\rho_{\rm
tot}(\lambda)=\frac{1}{2\pi}+\int\differential\lambda^{\prime}\mathcal{T}(\lambda-\lambda^{\prime})\rho_{\rm
p}(\lambda).$ (I.2)

It can be shown that the effective velocity fulfils the following integral
equation [23]

$v(\lambda)=\mathscr{h}_{1}(\lambda)+2\pi\int\differential\lambda^{\prime}\mathcal{T}(\lambda-\lambda^{\prime})\rho_{\rm
p}(\lambda^{\prime})\left(v(\lambda)-v(\lambda^{\prime})\right).$ (I.3)

With this, it is straightforward to check the identity used in the main text

$\int\differential\lambda v(\lambda)\rho_{\rm
p}(\lambda)=\int\differential\lambda\mathscr{h}_{1}(\lambda)\rho_{\rm
p}(\lambda).$ (I.4)

Indeed, we have

$\int\differential\lambda v(\lambda)\rho_{\rm p}(\lambda)=\int
d\lambda\mathscr{h}_{1}(\lambda)\rho_{\rm
p}(\lambda)+2\pi\int\differential\lambda\differential\lambda^{\prime}\mathcal{T}(\lambda-\lambda^{\prime})\rho_{\rm
p}(\lambda)\rho_{\rm
p}(\lambda^{\prime})\left(v(\lambda)-v(\lambda^{\prime})\right)=\int\differential\lambda\mathscr{h}_{1}(\lambda)\rho_{\rm
p}(\lambda),$ (I.5)

where the integral
$\int\differential\lambda\differential\lambda^{\prime}(\ldots)$ vanished due
to the symmetry of the kernel $\mathcal{T}(-\lambda)=\mathcal{T}(\lambda)$.

In addition to the Euler-scale term involving the effective velocity, there is
also the diffusion term [24]
$\frac{1}{2}\partial_{x}\left(\mathfrak{D}\partial_{x}\rho_{\rm p}\right)$
with $(\mathfrak{D}\partial_{x}\rho_{\rm
p})(\theta)=\int\differential\alpha\mathfrak{D}(\theta,\alpha)\partial_{x}\rho_{\rm
p}(\alpha)$. The diffusion kernel is

$\mathfrak{D}=(1-n\mathcal{T})^{-1}\rho_{\rm
tot}^{-1}\tilde{\mathfrak{D}}\rho_{\rm
tot}^{-1}(1-n\mathcal{T}),\qquad\tilde{\mathfrak{D}}(\theta,\alpha)=\delta(\theta-\alpha)w(\theta)-W(\theta,\alpha),$
(I.6)

and

$\displaystyle W(\theta,\alpha)=\rho_{\rm
p}(\theta)f(\theta)\left(\mathcal{T}^{\rm
dr}(\theta,\alpha)\right)^{2}|v(\theta)-v(\alpha)|,\qquad w(\theta)=\int{\rm
d}\alpha W(\alpha,\theta),$ (I.7)

where $f(\theta)$ denotes statistical factor [23] defined in (IV.5), which
depends on the particles’ statistics. The dressed scattering kernel fulfills

$\mathcal{T}^{\rm
dr}(\lambda,\lambda^{\prime})=\mathcal{T}(\lambda-\lambda^{\prime})+\int\differential\lambda^{\prime\prime}\mathcal{T}(\lambda-\lambda^{\prime\prime})n(\lambda^{\prime\prime})\mathcal{T}^{\rm
dr}(\lambda^{\prime\prime},\lambda^{\prime}).$ (I.8)

Note also that our definition of $\tilde{\mathfrak{D}}$ differs from the
standard one [24] by a factor of $\rho_{\rm tot}^{2}(\theta)$.

From the Galilean invariance follows that diffusion kernel has left
eigenvector with zero eigenvalue [24],

$\int d\theta\mathfrak{D}(\theta,\alpha)=0.$ (I.9)

This implies that the local density of particles is conserved by the diffusion
term. Moreover from Markovianity of the diffusion kernel, there is a
corresponding right eigenvector with zero eigenvalue, namely

$\int
d\alpha\mathfrak{D}(\theta,\alpha)(\mathcal{C}\mathscr{h}_{0})(\alpha)=0,$
(I.10)

where $\mathcal{C}$ is the charge susceptibility matrix defined as [23, 61]

$\mathcal{C}=(1-n\mathcal{T})^{-1}\rho_{\rm p}f(1-\mathcal{T}n)^{-1}.$ (I.11)

We have used here the fact [23] that $\rho_{\rm
tot}(\lambda)=\mathscr{h}_{0}^{\rm dr}(\lambda)/(2\pi)$. For later use, we
also introduce here the current susceptibility matrix $\mathcal{B}$ [23, 61]

$\mathcal{B}=(1-n\mathcal{T})^{-1}\rho_{\rm p}fv(1-\mathcal{T}n)^{-1},$ (I.12)

we will discuss matrix elements of the operators $\mathcal{B}$ and
$\mathcal{C}$ in more detail in Sec. IV.

##  II Hydrodynamic conservation laws

We derive here the conservation laws for the 3 hydrodynamic modes. We adopt
notation from the main text for the expectation values of charges

$\mathscr{q}_{n}=\int\differential\lambda\mathscr{h}_{n}(\lambda)\rho_{\rm
p}(\lambda).$ (II.1)

We start with the GHD equation including the diffusion, force term
$\mathfrak{f}=-\partial_{x}U(x)$ and the Boltzmann scattering integral

$\partial_{t}\rho_{\rm p}+\partial_{x}\left(v_{\rho}\rho_{\rm
p}\right)+\mathfrak{f}\partial_{\lambda}\rho_{\rm
p}=\frac{1}{2}\partial_{x}\left(\mathfrak{D}_{\rho}\partial_{x}\rho_{\rm
p}\right)+\mathcal{I}[\rho_{\rm p}].$ (II.2)

To obtain the conservation laws we multiply the GHD equation by
$\mathscr{h}_{j}(\lambda)$ and integrate over $\lambda$.

  * •

Mass conservation: for $j=0$ we obtain

$\partial_{t}\langle\mathscr{q}_{0}\rangle+\partial_{x}\langle\mathscr{q}_{1}\rangle=0,$
(II.3)

where we used that both the diffusion and the collision term conserve the
particles’ density. This equation reflects the conservation of particles’
density. We have also used that

$\int{\rm d}\lambda\,v(\lambda)\rho_{\rm p}(\lambda)=\int{\rm
d}\lambda\,\mathscr{h}_{1}(\lambda)\rho_{\rm p}(\lambda),$ (II.4)

  * •

Momentum conservation: for $j=1$ we obtain

$\partial_{t}\langle\mathscr{q}_{1}\rangle+\partial_{x}\langle\mathscr{q}_{1}v\rangle=-\partial_{x}\mathcal{P}_{\mathfrak{D}}+\mathfrak{f}\varrho,$
(II.5)

where we defined

$\mathcal{P}_{\mathfrak{D}}=-\frac{1}{2}\int{\rm d}\lambda{\rm
d}\mu\,\mathscr{h}_{1}(\lambda)\mathfrak{D}_{\rho}(\lambda,\mu)\partial_{x}\rho_{\rm
p}(\mu).$ (II.6)

The average in the second term can be decomposed into two contributions: the
convective flow and the pressure term. To this end we define the velocity
field

$u=\frac{\langle\mathscr{q}_{1}\rangle}{\langle\mathscr{q}_{0}\rangle},$
(II.7)

and introduce rapidity variable $\xi$ defined with respect to $u$,
$\lambda=u+\xi$. From the definition of $u$ it then follows that

$\int{\rm d}\xi\,\xi\rho_{\rm p}(u+\xi)=0,$ (II.8)

which is easily shown by transforming back to rapidity $\lambda$. We then have

$\displaystyle\langle\mathscr{q}_{1}v\rangle=\int{\rm d}\lambda
v(\lambda)\lambda\rho_{\rm p}(\lambda)=\int{\rm d}\xi(u+\xi)v(u+\xi)\rho_{\rm
p}(u+\xi)=\frac{\langle\mathscr{q}_{1}\rangle^{2}}{\langle\mathscr{q}_{0}\rangle}+\mathcal{P}_{v},$
(II.9)

where we defined the pressure by

$\mathcal{P}_{v}=\int{\rm d}\xi\,\xi v(u+\xi)\rho_{\rm p}(u+\xi)=\int{\rm
d}\lambda\,(\lambda-u)v(\lambda)\rho_{\rm p}(\lambda).$ (II.10)

We will later see that for a system locally in the thermal equilibrium state
this reduces to the hydrostatic pressure $P$.

Summing up, the continuity equation reads

$\partial_{t}\langle\mathscr{q}_{1}\rangle+\partial_{x}\left(\frac{\langle\mathscr{q}_{1}\rangle^{2}}{\langle\mathscr{q}_{0}\rangle}+\mathcal{P}\right)=\mathfrak{f}\varrho,$
(II.11)

with $\mathcal{P}=\mathcal{P}_{v}+\mathcal{P}_{\mathfrak{D}}$.

  * •

Energy conservation: Finally, for $j=2$ we find

$\partial_{t}\langle\mathscr{q}_{2}\rangle+\partial_{x}\langle\mathscr{q}_{2}v\rangle=-\partial_{x}\left(u\mathcal{P}_{\mathfrak{D}}+\mathcal{J}_{\mathfrak{D}}\right)+\mathfrak{f}\varrho
u,$ (II.12)

where we have defined

$u\mathcal{P}_{\mathfrak{D}}+\mathcal{J}_{\mathfrak{D}}=-\frac{1}{2}\int{\rm
d}\lambda{\rm
d}\mu\,\mathscr{h}_{2}(\lambda)\mathfrak{D}_{\rho}(\lambda,\mu)\partial_{x}\rho_{\rm
p}(\mu).$ (II.13)

The average on the left hand side of (II.12) can be again divided into various
contributions. Transforming to $\xi$ variable we have

$\langle\mathscr{q}_{2}v\rangle=\frac{1}{2}\int{\rm}d\xi(u+\xi)^{2}v(u+\xi)\rho_{\rm
p}(u+\xi)=\frac{1}{2}u^{2}\langle\mathscr{q}_{1}\rangle+u\mathcal{P}_{v}+\frac{1}{2}\int{\rm}d\xi\xi^{2}v(u+\xi)\rho_{\rm
p}(u+\xi).$ (II.14)

The remaining integral can be still divided into two contributions according
to

$\frac{1}{2}\int{\rm}d\xi\xi^{2}v(u+\xi)\rho_{\rm
p}(u+\xi)=\frac{1}{2}u\int{\rm}d\xi\,\xi^{2}\rho_{\rm
p}(u+\xi)+\frac{1}{2}\int{\rm}d\xi\,\xi^{2}(v(u+\xi)-u)\rho_{\rm p}(u+\xi).$
(II.15)

The first integral, by going back to $\lambda$ variable, becomes

$\frac{1}{2}\int{\rm}d\xi\,\xi^{2}\rho_{\rm
p}(u+\xi)=\langle\mathscr{q}_{2}\rangle-\frac{1}{2}\frac{\langle\mathscr{q}_{1}\rangle^{2}}{\langle\mathscr{q}_{0}\rangle},$
(II.16)

whereas the second defines the $\mathcal{J}_{v}$ contribution to the heat
current

$\mathcal{J}_{v}=\frac{1}{2}\int{\rm d}\xi\xi^{2}(v(u+\xi)-u)\rho_{\rm
p}(u+\xi)=\frac{1}{2}\int{\rm d}\lambda(\lambda-u)^{2}(v(\lambda)-u)\rho_{\rm
p}(\lambda).$ (II.17)

The total heat current is
$\mathcal{J}=\mathcal{J}_{v}+\mathcal{J}_{\mathfrak{D}}$. Finally we obtain

$\partial_{t}\langle\mathscr{q}_{2}\rangle+\partial_{x}\left(\frac{\langle\mathscr{q}_{1}\rangle\langle\mathscr{q}_{2}\rangle}{\langle\mathscr{q}_{0}\rangle}+\frac{\langle\mathscr{q}_{1}\rangle}{\langle\mathscr{q}_{0}\rangle}\mathcal{P}+\mathcal{J}\right)=\mathfrak{f}\varrho
u.$ (II.18)

Summarizing we obtain the following set of equations

$\displaystyle\partial_{t}\langle\mathscr{q}_{0}\rangle+\partial_{x}\langle\mathscr{q}_{1}\rangle=0,\qquad\partial_{t}\langle\mathscr{q}_{1}\rangle+\partial_{x}\left(\frac{\langle\mathscr{q}_{1}\rangle^{2}}{\langle\mathscr{q}_{0}\rangle}+\mathcal{P}\right)=\mathfrak{f}\varrho,$
(II.19)
$\displaystyle\partial_{t}\langle\mathscr{q}_{2}\rangle+\partial_{x}\left(\frac{\langle\mathscr{q}_{1}\rangle\langle\mathscr{q}_{2}\rangle}{\langle\mathscr{q}_{0}\rangle}+\frac{\langle\mathscr{q}_{1}\rangle}{\langle\mathscr{q}_{0}\rangle}\mathcal{P}+\mathcal{J}\right)=\mathfrak{f}\varrho
u.$ (II.20)

The continuity equations written above are for moments of the particles
distribution. More standard form is found if we express them in terms of the
velocity field $u$ and internal energy field $e$ defined by

$e=\frac{\langle\mathscr{q}_{2}\rangle}{\langle\mathscr{q}_{0}\rangle}-\frac{1}{2}\frac{\langle\mathscr{q}_{1}\rangle^{2}}{\langle\mathscr{q}_{0}\rangle}.$
(II.21)

Introducing also $\varrho(x,t)=\langle\mathscr{q}_{0}\rangle$ we find

$\displaystyle\partial_{t}\varrho+\partial_{x}(\varrho
u)=0,\qquad\partial_{t}(\varrho u)+\partial_{x}\left(\varrho
u^{2}+\mathcal{P}\right)=\mathfrak{f}\varrho,$ (II.22)
$\displaystyle\partial_{t}\left(\varrho
e+\frac{1}{2}u^{2}\varrho\right)+\partial_{x}\left(u\varrho
e+\frac{1}{2}u^{3}\varrho+u\mathcal{P}+\mathcal{J}\right)=\mathfrak{f}\varrho
u.$ (II.23)

This can be still simplified using that

$\frac{1}{2}\partial_{t}\left(u^{2}\varrho\right)=\frac{1}{2}u^{2}\partial_{x}\left(\varrho
u\right)-u\partial_{x}\left(\varrho
u^{2}+\mathcal{P}\right)+\mathfrak{f}\varrho u.$ (II.24)

Substituting this into the energy conservation law we obtain the final set of
equations

$\displaystyle\partial_{t}\varrho+\partial_{x}(\varrho
u)=0,\qquad\partial_{t}(\varrho u)+\partial_{x}\left(\varrho
u^{2}+\mathcal{P}\right)=\mathfrak{f}\varrho,$ (II.25)
$\displaystyle\partial_{t}\left(\varrho e\right)+\partial_{x}\left(u\varrho
e+\mathcal{J}\right)+\mathcal{P}\partial_{x}u=0,$ (II.26)

reported in the main text.

##  III Thermal boosted state

In this Section we study thermal boosted states. We show that a thermal
boosted state characterized by chemical potentials $\beta_{0}$, $\beta_{1}$
and $\beta_{2}$ can be related to a thermal state with zero momentum and
appropriately chosen $\bar{\beta}_{0}$, $\bar{\beta}_{2}$. Moreover, we will
show that heat current $\mathcal{J}_{v}$ vanishes in such case and that the
pressure $\mathcal{P}_{v}$ of a thermal boosted state is equal to the
hydrostatic pressure $P$.

We start with the free energy density [23]

$\mathsf{f}=\frac{1}{2\pi}\int\differential\lambda F[\epsilon(\lambda)],$
(III.1)

where function $F[\epsilon(\lambda)]$ depends on the particle’ statistics. For
example, for fermionic case relevant for Lieb-Liniger model [57]
$F[\epsilon(\lambda)]=-\log(1+e^{-\epsilon})$. The pseudoenergy
$\epsilon(\lambda)$ is determined from equation

$\epsilon(\lambda)=\epsilon_{0}(\lambda)+\int\differential\lambda^{\prime}\mathcal{T}(\lambda^{\prime}-\lambda)F[\epsilon(\lambda^{\prime})],$
(III.2)

where for thermal boosted state the bare pseudoenergy reads

$\epsilon_{0}(\lambda)=\beta_{0}\mathscr{h}_{0}(\lambda)+\beta_{1}\mathscr{h}_{1}(\lambda)+\beta_{2}\mathscr{h}_{2}(\lambda).$
(III.3)

In a thermal state, on the other hand
$\beta_{0}=-\mu/T,\,\beta_{1}=0,\,\beta_{2}=1/T$ with $\mu,T$ equal to
chemical potential and temperature, respectively (we set $k_{B}=1$).

We move now to discuss properties of a thermal boosted state, let $\rho_{\rm
p}(\lambda)$ be its distribution function. Then its first moment is non-zero,

$u=\frac{\int{\rm d}\lambda\,\lambda\rho_{\rm p}(\lambda)}{\int{\rm
d}\lambda\,\rho_{\rm p}(\lambda)}.$ (III.4)

Define $\bar{\rho}_{\rm p}(\lambda)=\rho_{\rm p}(\lambda+u)$. We will show
that $\bar{\rho}_{\rm p}(\lambda)$ is a thermal state and find expressions for
$\bar{\beta}_{0}$ and $\bar{\beta}_{2}$ in terms of the chemical potentials of
the boosted state. This can be argued on the basis of the Galilean invariance
of the theory or computed explicitly. We follow the latter path.

Consider the total density $\rho_{\rm tot}(\lambda)$, it obeys the integral
relation (I.2). Using the definition of $\bar{\rho}_{\rm p}(\lambda)$ it is
easy to show that the corresponding total density is given by $\bar{\rho}_{\rm
tot}(\lambda)=\rho_{\rm tot}(\lambda+u)$. From this follows that the
occupation functions and dressed pseudo-energies $\epsilon(\lambda)$ and
$\bar{\epsilon}(\lambda)$ obey the same relation
$\bar{\epsilon}(\lambda)=\epsilon(\lambda+u)$. Inspecting equation (III.2) we
finally find that also $\bar{\epsilon}_{0}(\lambda)=\epsilon_{0}(\lambda+u)$.
Writing this relation explicitly we have

$\bar{\beta}_{0}+\frac{1}{2}\bar{\beta}_{2}\lambda^{2}=\beta_{0}+\beta_{1}(\lambda+u)+\frac{1}{2}\beta_{2}(\lambda+u)^{2}.$
(III.5)

Matching the terms with the same power of $\lambda$ we find

$\bar{\beta}_{0}=\beta_{0}+u\beta_{1}+\frac{1}{2}u^{2}\beta_{2},\quad\bar{\beta}_{2}=\beta_{2},\quad\beta_{1}+u\beta_{2}=0.$
(III.6)

The third equality is always fulfilled and gives $u=-\beta_{1}/\beta_{2}$.

The last property follows from the observation that for $\lambda=u$ function
$\rho_{\rm p}(\lambda)$ has the maximum. On the other hand

$\frac{\partial\rho_{\rm
p}(\lambda)}{\partial\lambda}\sim\frac{\partial\epsilon(\lambda)}{\partial\lambda}.$
(III.7)

Function $\epsilon^{\prime}(\lambda)$ obeys a linear integral equation

$\epsilon^{\prime}(\lambda)=\epsilon_{0}^{\prime}(\lambda)+\int{\rm
d}\mu\mathcal{T}(\lambda-\mu)n(\mu)\epsilon^{\prime}(\mu).$ (III.8)

which follows from (III.2) and the identity [23]

$n=\frac{dF}{d\epsilon}$ (III.9)

and because of its linearity it maps zero to zero hence we need to solve
$\epsilon_{0}^{\prime}(\lambda)=0$. Using
$\epsilon_{0}(\lambda)=\beta_{0}+\beta_{1}\lambda+\frac{1}{2}\beta_{2}\lambda^{2}$
we readily find the relation between $u$ and the chemical potentials.

We can also show using (I.3) that the effective velocity transforms as
$\bar{v}(\lambda)+u=v(\lambda+u)$. This has a consequence that $v(u)=u$ where
we used that $\bar{v}(0)=0$.

With these expressions we can show that the hydrodynamic pressure of a thermal
boosted state reduces to the hydrostatic pressure $P$, see [15] and the end of
the present section. Hence the bulk viscosity vanishes is such case. We have

$\mathcal{P}_{v}=\int{\rm d}\xi\,\xi v(\xi+u)\rho_{\rm p}(\xi+u)=\int{\rm
d}\xi\,\xi(\bar{v}(\xi)+u)\bar{\rho}_{\rm p}(\xi)=\int{\rm
d}\xi\,\xi\bar{v}(\xi)\bar{\rho}_{\rm p}(\xi)=P.$ (III.10)

For the heat current we have

$\mathcal{J}_{v}=\frac{1}{2}\int{\rm d}\xi\xi^{2}(v(u+\xi)-u)\rho_{\rm
p}(u+\xi)=\frac{1}{2}\int{\rm d}\xi\xi^{2}\bar{v}(\xi)\bar{\rho}_{p}(\xi)=0,$
(III.11)

because $\bar{v}(\xi)$ is an odd function while $\bar{\rho}_{p}(\xi)$ is even.

The entropy of the boosted state is equal to the entropy of the unboosted
state. Recall the formula for the entropy density

$\mathsf{s}=\int{\rm d}\lambda\rho_{\rm tot}(\lambda)g(n(\lambda)),$ (III.12)

where function $g(n)$ depends on particles’ statistics [23, 24]. For the
fermionic models we have $g[n]=n\log n+(1-n)\log(1-n)$. In general
$g=\epsilon(\lambda)n(\lambda)-F[\epsilon(\lambda)]$. Using now the relations
between the particle distribution and filling function of the boosted and
unboosted state we see that the entropies are indeed equal.

### Expression for the hydrostatic pressure in terms of $v(\lambda)$

Let us close this section by further validating equation (III.10). Our aim is
to show that indeed

$P=\int\differential\lambda\lambda v(\lambda)\rho_{\rm p}(\lambda)$ (III.13)

where $\rho_{\rm p}(\lambda)$ is thermal state, is the hydrostatic pressure as
defined in thermodynamics, see eq. (IV.10) of the next section. To this end we
generalize the computations for the Lieb-Liniger model presented in [15].

We start by noting that $\rho_{\rm
tot}(\lambda)=\frac{1}{2\pi}\mathscr{h}_{0}^{\rm dr}(\lambda)$ and
$v(\lambda)=\mathscr{h}_{1}^{\rm dr}(\lambda)/\mathscr{h}_{0}^{\rm
dr}(\lambda)$, hence

$P=\frac{1}{2\pi}\int\differential\lambda\,\lambda\,n(\lambda)\mathscr{h}_{1}^{\rm
dr}(\lambda).$ (III.14)

Note that we have also shown at this point that $P=\mathscr{j}_{1}$, where
$\mathscr{j}_{1}$ is the momentum current (compare with formula for current
expectation value in [23]). To get $\mathscr{h}_{1}^{\rm dr}(\lambda)$ we look
at (III.8). For thermal state it reads

$\epsilon^{\prime}(\lambda)=\frac{1}{T}\mathscr{h}_{1}(\lambda)+\int\differential\lambda^{\prime}\mathcal{T}(\lambda-\lambda^{\prime})n(\lambda^{\prime})\epsilon^{\prime}(\lambda^{\prime}).$
(III.15)

Comparing with the definition of dressing operation (I.1) we conclude that
$\mathscr{h}_{1}^{\rm dr}(\lambda)=T\epsilon^{\prime}(\lambda)$. Moreover,
using (III.9) and integrating by parts we get

$P=-\frac{T}{2\pi}\int\differential\lambda F[\epsilon(\lambda)]=-T\mathsf{f}$
(III.16)

in agreement with (IV.10).

##  IV Thermodynamics and hydrodynamic matrices

In this section we will relate the matrix elements of hydrodynamic matrices
[61] $A$, $B$ and $C$ to the thermodynamics of the integrable model. This will
be very helpful in linking results from linearized GHD with thermodynamic
quantities, which appear naturally in the Navier-Stokes equations.

We start by defining the free energy flux [23], in addition to the free energy
density $\mathsf{f}$ defined in (III.1),

$\mathsf{g}=\frac{1}{2\pi}\int\differential\lambda
E^{\prime}(\lambda)F[\epsilon(\lambda)],$ (IV.1)

where the pseudoenergy $\epsilon(\lambda)$ satisfies (III.2).

For a generic GGE state, the bare pseudoenergy includes higher conserved
charges [48]. The simplest choice being the ultra-local charges
$\mathscr{q}_{j}$,

$\epsilon_{0}(\lambda)=\sum_{j=0}^{\infty}\beta_{j}\mathscr{h}_{j}(\lambda).$
(IV.2)

The expectation values of the charges in the GGE state follow then from
derivatives of the free energy density with respect to the generalized
temperatures. In a similar fashion, the expectation values of the currents are
given by the derivatives of the free energy flux,

$\mathscr{q}_{i}=\frac{\partial\mathsf{f}}{\partial\beta_{i}},\qquad\mathscr{j}_{i}=\frac{\partial\mathsf{g}}{\partial\beta_{i}}.$
(IV.3)

The hydrodynamic matrices $\mathcal{B}$ and $\mathcal{C}$ are susceptibilities
of the currents and charges

$\mathcal{B}_{i,j}=-\frac{\partial\mathscr{j}_{i}}{\partial\beta_{j}}=\int{\rm
d}\lambda\,\rho_{\rm p}(\lambda)f(\lambda)v(\lambda)\mathscr{h}_{i}^{\rm
dr}(\lambda)\mathscr{h}_{j}^{\rm
dr}(\lambda),\qquad\mathcal{C}_{i,j}=-\frac{\partial\mathscr{q}_{i}}{\partial\beta_{j}}=\int
d\lambda\,\rho_{\rm p}(\lambda)f(\lambda)\mathscr{h}_{i}^{\rm
dr}(\lambda)\mathscr{h}_{j}^{\rm dr}(\lambda).$ (IV.4)

The matrices defined above are clearly symmetric. Above we defined

$f(\lambda)=-\frac{d^{2}F/d\epsilon^{2}}{dF/d\epsilon},$ (IV.5)

which is called the statistical factor [23]. For instance, in the case of
fermionic statistics one finds $f(\theta)=1-n(\theta)$. We also note that for
Galilean invariant systems there is an useful relation
$\mathcal{B}_{1,0}=\mathcal{C}_{1,1}$.

The hydrodynamic matrices depend on the choice of the basis for the ultra-
local charges. For the thermodynamics it is natural to work with matrix
elements of the ultralocal charges $\mathscr{h}_{n}(\lambda)$ as defined
above. Instead, for the hydrodynamics, it is more convenient to work with the
orthonormalized basis obtained from the ultra-local charges through the Gram-
Schmidt method with the inner product (V.6). For $h_{n}(\lambda)$ denoting the
single-particle eigenvalues of the resulting charges, the hydrodynamic
matrices are

$B_{i,j}=-\frac{\partial j_{i}}{\partial\beta_{j}}=\int{\rm
d}\lambda\,\rho_{\rm p}(\lambda)f(\lambda)v(\lambda)h_{i}^{\rm
dr}(\lambda)h_{j}^{\rm dr}(\lambda),\qquad C_{i,j}=-\frac{\partial
q_{i}}{\partial\beta_{j}}=\int d\lambda\,\rho_{\rm
p}(\lambda)f(\lambda)h_{i}^{\rm dr}(\lambda)h_{j}^{\rm
dr}(\lambda)=\delta_{i,j}.$ (IV.6)

In the remainder of this section we will relate matrix elements of
$\mathcal{B}$ and $\mathcal{C}$ to thermodynamic quantities. We focus on the
Gibbs state with a thermal $\epsilon_{0}(\lambda)$ function,

$\epsilon_{0}(\lambda)=-\frac{\mu}{T}\mathscr{h}_{0}(\lambda)+\frac{1}{T}\mathscr{h}_{2}(\lambda).$
(IV.7)

We also introduce the notation

$\varrho=\mathscr{q}_{0}=\int\differential\lambda\mathscr{h}_{0}(\lambda)\rho_{\rm
p}(\lambda),\qquad\mathsf{e}=\mathscr{q}_{2}=\int\differential\lambda\mathscr{h}_{2}(\lambda)\rho_{\rm
p}(\lambda),\qquad\mathsf{s}=\int\differential\lambda\rho_{\rm
tot}(\lambda)g[n(\lambda)]$ (IV.8)

for particle, energy and entropy density, respectively. Moreover, by

$s=\varrho^{-1}\mathsf{s},\qquad e=\varrho^{-1}\mathsf{e}$ (IV.9)

we denote entropy and energy per particle. For Gibbs state the free energy
$\mathsf{f}$ takes the form of the grand potential per unit length divided by
temperature and can be written as

$\mathsf{f}=-\frac{\mu}{T}\varrho+\frac{1}{T}\mathsf{e}-\mathsf{s}=-\frac{P}{T}.$
(IV.10)

With this, we compute various thermodynamic quantities. We start with the
sound velocity.

###  IV.1 Sound velocity

The sound velocity is related to the adiabatic compressibility $\varkappa_{S}$
and defined as (we work with unit mass)

$v_{s}^{2}=\left(\frac{\partial
P}{\partial\varrho}\right)_{s}=\frac{1}{\varrho\varkappa_{S}},\qquad\varkappa_{S}=\frac{1}{\varrho}\left(\frac{\partial\varrho}{\partial
P}\right)_{s}$ (IV.11)

To compute the relevant derivative, we realize that $P,\varrho$ and $s$ depend
only on two thermodynamic parameters, which specify the state: $\mu,T$.
Equivalently, they depend on $\beta_{0},\beta_{2}$. We may thus write

$\differential s=\left(\frac{\partial
s}{\partial\beta_{0}}\right)_{\beta_{2}}\differential\beta_{0}+\left(\frac{\partial
s}{\partial\beta_{2}}\right)_{\beta_{0}}\differential\beta_{2},$ (IV.12)

and similarly for $\differential\varrho$ and $\differential P$. To keep the
entropy fixed $\differential\beta_{0}$ and $\differential\beta_{2}$ must be
related through

$\differential\beta_{0}=-\frac{\left(\frac{\partial
s}{\partial\beta_{2}}\right)_{\beta_{0}}}{\left(\frac{\partial
s}{\partial\beta_{0}}\right)_{\beta_{2}}}\differential\beta_{2},$ (IV.13)

which leads to the following formula

$\left(\frac{\partial\varrho}{\partial
P}\right)_{s}=\frac{\left(\frac{\partial\varrho}{\partial\beta_{0}}\right)_{\beta_{2}}\left(\frac{\partial
s}{\partial\beta_{2}}\right)_{\beta_{0}}-\left(\frac{\partial\varrho}{\partial\beta_{2}}\right)_{\beta_{0}}\left(\frac{\partial
s}{\partial\beta_{0}}\right)_{\beta_{2}}}{\left(\frac{\partial
P}{\partial\beta_{0}}\right)_{\beta_{2}}\left(\frac{\partial
s}{\partial\beta_{2}}\right)_{\beta_{0}}-\left(\frac{\partial
P}{\partial\beta_{2}}\right)_{\beta_{0}}\left(\frac{\partial
s}{\partial\beta_{0}}\right)_{\beta_{2}}}.$ (IV.14)

The derivatives of the entropy can be computed from the expression for free
energy (IV.10). The entropy per particle equals to

$s=\frac{\beta_{0}\varrho+\beta_{2}\mathsf{e}-\mathsf{f}}{\varrho},$ (IV.15)

and the two required derivatives are

$\left(\frac{\partial
s}{\partial\beta_{0}}\right)_{\beta_{2}}=\frac{\mathcal{C}_{0,0}(\beta_{2}\mathsf{e}-\mathsf{f})-\mathcal{C}_{2,0}\beta_{2}\varrho}{\varrho^{2}},\qquad\left(\frac{\partial
s}{\partial\beta_{2}}\right)_{\beta_{0}}=\frac{\mathcal{C}_{2,0}(\beta_{2}\mathsf{e}-\mathsf{f})-\mathcal{C}_{2,2}\beta_{2}\varrho}{\varrho^{2}}.$
(IV.16)

We use now the relation $\mathscr{j}_{1}=P$, valid for Galilean-invariant
models, see (III.14). Together with $P=-T\mathsf{f}=-\mathsf{f}/\beta_{2}$
they yield

$\mathcal{B}_{1,0}=-\left(\frac{\partial
P}{\partial\beta_{0}}\right)_{\beta_{2}}=\frac{\varrho}{\beta_{2}}=T\varrho,\qquad\mathcal{B}_{2,1}=-\left(\frac{\partial
P}{\partial\beta_{2}}\right)_{\beta_{0}}=\frac{\beta_{2}\mathsf{e}-\mathsf{f}}{\beta_{2}^{2}},$
(IV.17)

which gives

$\left(\frac{\partial
s}{\partial\beta_{0}}\right)_{\beta_{2}}=\frac{1}{\mathcal{B}_{1,0}^{2}}\left(\mathcal{C}_{0,0}\mathcal{B}_{2,1}-\mathcal{C}_{2,0}\mathcal{B}_{1,0}\right)\qquad\left(\frac{\partial
s}{\partial\beta_{2}}\right)_{\beta_{0}}=\frac{1}{\mathcal{B}_{1,0}^{2}}\left(\mathcal{C}_{2,0}\mathcal{B}_{2,1}-\mathcal{C}_{2,2}\mathcal{B}_{1,0}\right)$
(IV.18)

We now go back to (IV.14) and (IV.11) to find the final expression for the
adiabatic compressibilty

$\varkappa_{S}=T\frac{\mathcal{C}_{2,2}\mathcal{C}_{0,0}-\mathcal{C}_{2,0}^{2}}{\mathcal{B}_{1,0}^{2}\mathcal{C}_{2,2}+\mathcal{B}_{2,1}^{2}\mathcal{C}_{0,0}-2\mathcal{B}_{1,0}\mathcal{B}_{2,1}\mathcal{C}_{2,0}},$
(IV.19)

from which the final expression for the sound velocity follows

$v_{s}=\sqrt{\frac{1}{\mathcal{B}_{1,0}}\frac{\mathcal{B}_{1,0}^{2}\mathcal{C}_{2,2}+\mathcal{B}_{2,1}^{2}\mathcal{C}_{0,0}-2\mathcal{B}_{1,0}\mathcal{B}_{2,1}\mathcal{C}_{2,0}}{\mathcal{C}_{2,2}\mathcal{C}_{0,0}-\mathcal{C}_{2,0}^{2}}}.$
(IV.20)

We have analytically confirmed that this general formula reduces to known
formulas for free fermions $v_{s}^{\rm
Fermi}=\sqrt{3T\frac{f_{3/2}(e^{\mu/T})}{f_{1/2}(e^{\mu/T})}}$ where
$f_{\alpha}(x)$ are Fermi functions and for hard-rods [1] $v_{s}^{\rm
HR}=\sqrt{3T/(1-a\varrho)^{2}}$, where $a$ is rod’s length. For the Lieb-
Liniger model we confirmed with derivatives computed numerically from TBA.

###  IV.2 Specific heats, isothermal compressibility and additional identities

We calculate in a similar fashion several more thermodynamic quantities, which
are useful for purposes of this work. Isothermal compressibility: The
isothermal compressibility is defined as

$\varkappa_{T}=\frac{1}{\varrho}\left(\frac{\partial\varrho}{\partial
P}\right)_{T},$ (IV.21)

and standard manipulations give

$\varkappa_{T}=\frac{1}{\varrho}\frac{\left(\frac{\partial\varrho}{\partial\beta_{0}}\right)_{\beta_{2}}}{\left(\frac{\partial
P}{\partial\beta_{0}}\right)_{\beta_{2}}}=\frac{1}{\varrho}\frac{\mathcal{C}_{0,0}}{\mathcal{B}_{1,0}}.$
(IV.22)

Specific heat at constant volume: The definition is

$c_{V}=\left(\frac{\partial e}{\partial
T}\right)_{\varrho}=\frac{1}{\varrho}\left(\frac{\partial\mathsf{e}}{\partial
T}\right)_{\varrho}.$ (IV.23)

Expressing $\differential T$ with $\differential\beta_{0}$ and
$\differential\beta_{2}$ we obtain

$c_{V}=\frac{1}{\varrho}\frac{-\left(\frac{\partial\mathsf{e}}{\partial\beta_{0}}\right)_{\beta_{2}}\left(\frac{\partial\varrho}{\partial\beta_{2}}\right)_{\beta_{0}}+\left(\frac{\partial\mathsf{e}}{\partial\beta_{2}}\right)_{\beta_{0}}\left(\frac{\partial\varrho}{\partial\beta_{0}}\right)_{\beta_{2}}}{-T^{2}\left(\frac{\partial\varrho}{\partial\beta_{0}}\right)_{\beta_{2}}},$
(IV.24)

which can now be expressed with the hydrodynamic matrices as

$c_{V}=\frac{1}{T}\frac{\mathcal{C}_{2,2}\mathcal{C}_{0,0}-\mathcal{C}_{2,0}^{2}}{\mathcal{C}_{0,0}\mathcal{B}_{1,0}}.$
(IV.25)

Specific heat at constant pressure: The definition is

$c_{P}=T\left(\frac{\partial s}{\partial T}\right)_{P}.$ (IV.26)

Proceeding in an analogous way as in the previous cases we find

$c_{P}=\frac{1}{T}\frac{\mathcal{B}_{1,0}^{2}\mathcal{C}_{2,2}+\mathcal{B}_{2,1}^{2}\mathcal{C}_{0,0}-2\mathcal{B}_{1,0}\mathcal{B}_{2,1}\mathcal{C}_{2,0}}{\mathcal{B}_{1,0}^{3}}.$
(IV.27)

As a check for our expressions, we also observe that the two compressibilities
and two specific heats are not independent from each other. Instead, they obey
the thermodynamic identity

$\varkappa_{S}=\frac{c_{V}}{c_{P}}\varkappa_{T}.$ (IV.28)

One can easily verify that this identity holds for the expressions given in
(IV.19), (IV.22), (IV.25) and (IV.27).

Expressions for the matrix elements of $A$: The sound velocity may be
rewritten with $A$ matrix elements in GS-orthonormalized basis. These matrix
elements are defined in (V.14) and appear explicitly in the equation (VII.2),
from which dispersion relations are determined. It is thus useful to write
thermodynamic quantities in terms of $A$ matrix in this basis

$A_{1,0}=\sqrt{\frac{\mathcal{B}_{1,0}}{\mathcal{C}_{0,0}}}=\frac{1}{\sqrt{\varrho\varkappa_{T}}},\qquad
A_{2,1}=\frac{\mathcal{B}_{2,1}\mathcal{C}_{0,0}-\mathcal{B}_{1,0}\mathcal{C}_{2,0}}{\sqrt{\mathcal{C}_{0,0}\mathcal{B}_{1,0}(\mathcal{C}_{2,2}\mathcal{C}_{0,0}-\mathcal{C}_{2,0}^{2})}}=\frac{T}{\varrho^{2}c_{V}}\left(\frac{\partial
P}{\partial
T}\right)_{\varrho}=\frac{1}{\sqrt{\varrho\varkappa_{T}}}\sqrt{c_{P}/c_{V}-1}.$
(IV.29)

Above we have used a general thermodynamic identity [36]

$c_{P}-c_{V}=\frac{T}{\varrho^{2}}\frac{\left(\frac{\partial P}{\partial
T}\right)_{\varrho}^{2}}{\left(\frac{\partial
P}{\partial\varrho}\right)_{T}}.$ (IV.30)

In the perturbation theory, we will find that speed of sound becomes

$v_{s}=\sqrt{A_{1,0}^{2}+A_{2,1}^{2}},$ (IV.31)

which agrees with (IV.20). Moreover, we observe that

$\frac{A_{2,1}^{2}}{A_{1,0}^{2}}=c_{P}\left(\frac{1}{c_{V}}-\frac{1}{c_{P}}\right)=\frac{c_{P}}{c_{V}}-1.$
(IV.32)

Additional identities: For analysis of the structure of heat and sound modes
the following identities are useful

$\left(\frac{\partial P}{\partial
T}\right)_{\varrho}=\frac{1}{T^{2}}\frac{\mathcal{B}_{2,1}\mathcal{C}_{0,0}-\mathcal{B}_{1,0}\mathcal{C}_{2,0}}{\mathcal{C}_{0,0}},\qquad\frac{1}{\varrho}\left(\mathsf{e}+P-T\left(\frac{\partial
P}{\partial
T}\right)_{\varrho}\right)=\frac{\mathcal{C}_{2,0}}{\mathcal{C}_{0,0}}.$
(IV.33)

They can be derived with the methods presented above. Lastly, let us comment
that the method presented in this section gives thermodynamic quantities as
explicit functionals (that is, combinations of matrix elements of hydrodynamic
matrices) of the thermal rapidity distribution $\rho_{\rm p}(\lambda)$. This
representation is convenient and may be useful for applications which go
beyond the scope of this work.

##  V Derivation of the linearized equation for charge perturbations

We again invoke the GHD equation supplemented with the collision integral
term, this time with $U(x)=0$,

$\partial_{t}\rho_{\rm p}+\partial_{x}\left(v_{\rho}\rho_{\rm
p}\right)=\frac{1}{2}\partial_{x}\left(\mathfrak{D}_{\rho}\partial_{x}\rho_{\rm
p}\right)+\mathcal{I}[\rho_{\rm p}].$ (V.1)

Our aim here is to first linearize this equation and then rewrite it as a
continuity equation for the densities of local charges. To this end we
linearize it around a stationary uniform thermal state $\rho(x,t)=\rho_{\rm
th}+\delta\rho_{\rm p}(x,t)$ getting equation

$\partial_{t}\delta\rho_{\rm
p}+\partial_{x}\left((1-n\mathcal{T})^{-1}v_{\rho}(1-n\mathcal{T})\delta\rho_{\rm
p}\right)=\frac{1}{2}\mathfrak{D}_{\rm th}\partial_{x}^{2}\delta\rho_{\rm
p}+\left(\frac{\delta\mathcal{I}}{\delta\rho_{\rm p}}\right)_{\rm
th}\delta\rho_{\rm p}.$ (V.2)

The second term in the equation above (linearized expression for effective
velocity term) was derived in [61].

Motivated by the structure of $\Gamma$ matrix that we consider (see Sec. IX),
we wish to change the degrees of freedom from a perturbation of quasiparticle
distribution to a perturbation of the bare pseudoenergy $\epsilon_{0}$. To
establish a relation between them, we start with

$\rho_{\rm p}=n\rho_{\rm tot},\qquad\delta\rho_{\rm p}=\delta n\rho_{\rm
tot}+n\delta\rho_{\rm tot}.$ (V.3)

The perturbation occupation function may be related to the perturbation of
pseudoenergy via (III.9), hence $\delta n=-fn\delta\epsilon$. Moreover from
(I.2) we have $\delta\rho_{\rm tot}=\mathcal{T}\delta\rho_{\rm p}$ and from
(III.2) we get $\delta\epsilon=(1-\mathcal{T}n)^{-1}\delta\epsilon_{0}$.
Combining, we find the following relation

$\delta\rho_{\rm p}=-(1-n\mathcal{T})^{-1}\rho_{\rm
p}f\delta\epsilon=-(1-n\mathcal{T})^{-1}\rho_{\rm
p}f(1-\mathcal{T}n)^{-1}\delta\epsilon_{0}=-\mathcal{C}\delta\epsilon_{0}.$
(V.4)

where the charge susceptibility $\mathcal{C}$ operator was introduced in
(I.11) and the operator $(1-\mathcal{T}n)^{-1}$ implements the dressing (I.1).
Function $\delta\epsilon_{0}(\lambda)$ can be expanded as

$\delta\epsilon_{0}(\lambda)=\sum_{m}(h_{m}|\delta\epsilon_{0})h_{m}(\lambda),$
(V.5)

where $h_{m}(\lambda)$ form an orthonormal set of functions with respect to
the hydrodynamic scalar product $(\cdot|\cdot)$, see also [23],

$(h|g)=\int{\rm d}\lambda\,\rho_{\rm p}(\lambda)f(\lambda)h^{\rm
dr}(\lambda)g^{\rm dr}(\lambda).$ (V.6)

In practice, $h_{m}$ are constructed by the Gram-Schmidt orthonormalization of
the ultra-local charges with single particle eigenvalue
$\mathscr{h}_{m}(\lambda)=\lambda^{m}/m!$. The most important first three
functions are

$h_{0}(\lambda)=\frac{1}{\sqrt{\mathcal{C}_{0,0}}}\mathscr{h}_{0}(\lambda),\qquad
h_{1}(\lambda)=\frac{1}{\sqrt{\mathcal{C}_{1,1}}}\mathscr{h}_{1}(\lambda),\qquad
h_{2}(\lambda)=\sqrt{\frac{\mathcal{C}_{0,0}}{\mathcal{C}_{2,2}\mathcal{C}_{0,0}-\mathcal{C}_{2,0}^{2}}}\left(\mathscr{h}_{2}(\lambda)-\frac{\mathcal{C}_{2,0}}{\mathcal{C}_{0,0}}\mathscr{h}_{0}(\lambda)\right).$
(V.7)

The expansion coefficients are directly related to the change in the
expectation value of the $m$-th conserved charge,
$(h_{m}|\delta\epsilon_{0})=-\delta q_{m}$, where

$\delta q_{m}=\int{\rm d}\lambda\,h_{m}(\lambda)\delta\rho_{\rm p}(\lambda).$
(V.8)

This relation follows from the following identity [23] (here dot-product is
$h\cdot g=\int{\rm d}\lambda h(\lambda)g(\lambda)$)

$(1-n\mathcal{T})^{-1}ng\cdot h=g\cdot n(1-\mathcal{T}n)^{-1}h=g\cdot nh^{\rm
dr},$ (V.9)

which follows from the symmetry of the kernel $\mathcal{T}$. With the help of
this identity we find

$\displaystyle\delta q_{m}$ $\displaystyle=-\int{\rm
d}\lambda\,h_{m}(\lambda)(1-n\mathcal{T})^{-1}\rho_{\rm
p}f(1-\mathcal{T}n)^{-1}\delta\epsilon_{0}$ $\displaystyle=-\int{\rm
d}\lambda\,\rho_{\rm p}(\lambda)f(\lambda)h_{m}^{\rm
dr}(\lambda)\delta\epsilon_{0}^{\rm dr}(\lambda)=-(h_{m}|\delta\epsilon_{0})$
(V.10)

where we used in the second step with $h=h_{m}$ and $g=\rho_{\rm
p}f\delta\epsilon_{0}^{\rm dr}/n$. This allows us to write

$\delta\rho_{\rm p}=(1-n\mathcal{T})^{-1}\rho_{\rm
p}f(1-\mathcal{T}n)^{-1}\sum_{m}h_{m}\delta q_{m}.$ (V.11)

Functions $h_{n}$ allow us to reformulate the GHD equations as equations for
the charge perturbations. We multiply the GHD equation (V.2) by
$h_{n}(\lambda)$, plug (V.11) and integrate over $\lambda$ taking advantage of
the orthonormality of functions $h_{n}$. Moreover, we introduce $\Gamma$
operator

$\Gamma=\frac{\delta\mathcal{I}}{\delta\epsilon_{0}}=-\frac{\delta\mathcal{I}}{\delta\rho_{\rm
p}}\mathcal{C}.$ (V.12)

The result is

$\partial_{t}\delta q_{n}+\sum_{m}A_{nm}\partial_{x}\delta
q_{m}=\sum_{m}D_{nm}\partial_{x}^{2}\delta q_{m}-\sum_{m}\Gamma_{nm}\delta
q_{m}$ (V.13)

with

$\displaystyle A_{nm}$ $\displaystyle=B_{nm}=\int{\rm d}\lambda\,\rho_{\rm
p}(\lambda)f(\lambda)v(\lambda)h_{n}^{\rm dr}(\lambda)h_{m}^{\rm
dr}(\lambda),$ (V.14) $\displaystyle D_{nm}$
$\displaystyle=\frac{1}{2}\int{\rm d}\lambda\,{\rm d}\mu\,h_{n}^{\rm
dr}(\lambda)\rho_{\rm tot}^{-1}(\lambda)\tilde{\mathfrak{D}}_{\rm
th}(\lambda,\mu)\rho_{\rm tot}^{-1}(\mu)\rho_{\rm p}(\mu)f(\mu)h_{m}^{\rm
dr}(\mu),$ (V.15) $\displaystyle\Gamma_{nm}$ $\displaystyle=\int{\rm
d}\lambda\,{\rm d}\mu\,h_{n}(\lambda)\Gamma(\lambda,\mu)h_{m}^{\rm dr}(\mu).$
(V.16)

Furthermore, we can take into account that $\Gamma$ operator hides the
Dressing operation [27, 29, 67]. The Dressing operation is defined through
(where $F(\mu|\lambda)$ is backflow function),

$f^{\rm Dr}(\lambda)=f(\lambda)-\int{\rm
d}\mu\,n(\mu)F(\mu|\lambda)f^{\prime}(\mu)$ (V.17)

such that the collision integral is

$\Gamma(\lambda,\mu)=\int{\rm
d}\nu\,\Gamma_{0}(\nu,\mu)\mathbb{F}(\lambda,\nu),\qquad\mathbb{F}(\lambda,\nu)=\delta(\lambda-\nu)+\partial_{\lambda}(n(\lambda)F(\lambda|\nu)),$
(V.18)

where $\Gamma_{0}$ is the bare collision integral. This Dressing operation may
be shifted to $h_{n}(\lambda)$ so finally we get

$\Gamma_{nm}=\int{\rm d}\lambda\,{\rm d}\mu\,h_{n}^{\rm
Dr}(\lambda)\Gamma_{0}(\lambda,\mu)h_{m}^{\rm dr}(\mu).$ (V.19)

Let us comment on general properties of the matrix elements listed above.
Matrix elements $A_{nm}$ are non-zero where $n,m$ are of different parity,
instead $\Gamma_{mn}$ and $D_{nm}$ are non-zero only when $n,m$ have the same
parity. We also have $\Gamma_{0,n}=\Gamma_{1,n}=\Gamma_{2,n}=0$ which reflects
the $3$ collision invariants, whereas $D_{0,n}=D_{n,0}=0$ which follows from
the Markovianity of the diffusion kernel discussed in Section I.

For the spectrum of $\Gamma$ operator we assume that it is non-negative and
gapped. The eigenvectors correspond to either even or odd functions of
rapidities. We confirm these properties for the collision integral analyzed in
Sec. IX.

##  VI Linearized Navier-Stokes equations and transport coefficients

We linearize the Navier-Stokes equations, which can be understood as (II.25)
with the phenomenological expressions

$\displaystyle\mathcal{P}=P-\zeta\partial_{x}u,\qquad\mathcal{J}=-\kappa\partial_{x}T,$
(VI.1)

for the dynamic pressure and heat current. In this section, we do not specify
$\zeta$ and $\kappa$. Our aim is to focus on general structure of N-S
equations and find relation between the spectrum of the hydrodynamic modes and
the transport coefficients.

We start by introducing the temperature field to the Navier-Stokes equations.
For this, we note the following general thermodynamic equality, which can be
derived from basic principles [36]

$\delta T=\frac{1}{c_{V}}\delta
e-\frac{1}{c_{V}}\left(\frac{P}{\varrho^{2}}-\frac{T}{\varrho^{2}}\left(\frac{\partial
P}{\partial T}\right)_{\varrho}\right)\delta\varrho.$ (VI.2)

One-dimensional version of Navier-Stokes equations [54, 56] for the fields
$\varrho(x,t)$, $u(x,t)$ and $T(x,t)$ reads then

$\displaystyle\partial_{t}\varrho$ $\displaystyle=-\partial_{x}(\varrho u),$
(VI.3) $\displaystyle\partial_{t}u$
$\displaystyle=-u\partial_{x}u-\frac{1}{\varrho}\left(\frac{\partial
P}{\partial\varrho}\right)_{T}\partial_{x}\varrho-\frac{1}{\varrho}\left(\frac{\partial
P}{\partial
T}\right)_{\varrho}\partial_{x}T+\frac{\zeta}{\varrho}\partial_{x}^{2}u,$
$\displaystyle\partial_{t}T$ $\displaystyle=-\frac{T}{\varrho
c_{V}}\left(\frac{\partial P}{\partial
T}\right)_{\varrho}\partial_{x}u-u\partial_{x}T+\frac{\kappa}{\varrho
c_{V}}\partial_{x}^{2}T+\frac{\zeta}{\varrho c_{V}}(\partial_{x}u)^{2},$

where $\kappa,\zeta$ are thermal conductivity and bulk viscosity,
respectively. Linearizing around a homogenous thermal state

$\varrho(x,t)=\varrho+\delta\varrho(x,t),\qquad u(x,t)=0+\delta u(x,t),\qquad
T(x,t)=T+\delta T(x,t),$ (VI.4)

one gets

$\displaystyle\partial_{t}\delta\varrho$
$\displaystyle=-\varrho\partial_{x}\delta u,$ (VI.5)
$\displaystyle\partial_{t}\delta u$
$\displaystyle=-\frac{1}{\varrho}\left(\frac{\partial
P}{\partial\varrho}\right)_{T}\partial_{x}\delta\varrho-\frac{1}{\varrho}\left(\frac{\partial
P}{\partial T}\right)_{\varrho}\partial_{x}\delta
T+\frac{\zeta}{\varrho}\partial_{x}^{2}\delta u,$
$\displaystyle\partial_{t}\delta T$ $\displaystyle=-\frac{T}{\varrho
c_{V}}\left(\frac{\partial P}{\partial T}\right)_{\varrho}\partial_{x}\delta
u+\frac{\kappa}{\varrho c_{V}}\partial_{x}^{2}\delta T.$

Moving to the $k$-space and anticipating the time-dependence as

$\delta\varrho(x,t)=\int\differential ke^{ikx-\Lambda
t}\delta\varrho(k),\qquad\delta u(x,t)=\int\differential ke^{ikx-\Lambda
t}\delta u(k),\qquad\delta T(x,t)=\int\differential ke^{ikx-\Lambda t}\delta
T(k),$ (VI.6)

we get an eigenproblem with non-Hermitian matrix. As the matrix is non-
Hermitean the left- and right-eigenvectors are different. This problem be
cured by changing the variables as

$\delta\varrho=\sqrt{T\varrho/(\partial
P/\partial\varrho)_{T}}\,\delta\tilde{\varrho}=\varrho\sqrt{T\varkappa_{T}}\,\delta\tilde{\varrho},\qquad\delta
u=\sqrt{T/\varrho}\,\delta\tilde{u},\qquad\delta T=T/\sqrt{\varrho
c_{V}}\,\delta\tilde{T}.$ (VI.7)

In these variables, the relevant matrix becomes symmetric. Explicitly, when
looking for eigenvalues $\Lambda$ we will deal with

$\text{det}\begin{pmatrix}-\Lambda&i\frac{1}{\sqrt{\varrho\varkappa_{T}}}k&0\\\
i\frac{1}{\sqrt{\varrho\varkappa_{T}}}k&-\Lambda+\frac{\zeta}{\varrho}k^{2}&i\sqrt{\frac{T}{\varrho^{2}c_{V}}}\left(\frac{\partial
P}{\partial T}\right)_{\varrho}k\\\
0&i\sqrt{\frac{T}{\varrho^{2}c_{V}}}\left(\frac{\partial P}{\partial
T}\right)_{\varrho}k&-\Lambda+\frac{\kappa}{\varrho
c_{V}}k^{2}\end{pmatrix}=0.$ (VI.8)

We use now the general thermodynamic identities (IV.28), (IV.30) and
$\frac{c_{P}}{c_{V}}\left(\frac{\partial
P}{\partial\varrho}\right)_{T}=v_{s}^{2}$ [36] and find the following
dispersion relations

$\displaystyle\Lambda_{\pm}(k)=\pm
iv_{s}k+\frac{1}{2\varrho}\bigg{[}\left(\frac{1}{c_{V}}-\frac{1}{c_{P}}\right)\kappa+\zeta\bigg{]}k^{2}+O(k^{3}),\qquad\Lambda_{\rm
th}(k)=\frac{\kappa}{\varrho c_{P}}k^{2}+O(k^{3}).$ (VI.9)

These dispersion relations correspond to two sound modes and one heat mode,
respectively. We observe that in the case of vanishing bulk viscosity
$\zeta=0$, the ratio of the quadratic coefficients is $(c_{P}/c_{V}-1)/2$ and
therefore is universal. Depends only on the thermodynamics of the underlying
gas but not on the integrability breaking terms. Vanishing bulk viscosity is
an expected result for an ideal gas perturbed with weak interactions [54]. We
will see that in our setup the bulk viscosity vanishes in the Tonks-Girardeau
($c\rightarrow\infty$) and in $c\to 0$ limits, see Sec. VIII.

In addition to the eigenvalues, we can find eigenvectors. In the limit of
$k\to 0$ we have

$\delta_{\pm}(k\to
0)=\frac{1}{\sqrt{2}}\left(\sqrt{\frac{c_{V}}{c_{P}}}\delta\tilde{\varrho}\pm\delta\tilde{u}+\sqrt{\frac{c_{P}-c_{V}}{c_{P}}}\delta\tilde{T}\right),\qquad\delta_{\rm
th}(k\to
0)=-\sqrt{\frac{c_{P}-c_{V}}{c_{P}}}\delta\tilde{\varrho}+\sqrt{\frac{c_{V}}{c_{P}}}\delta\tilde{T}.$
(VI.10)

The general structure revealed in this section will be our reference point for
perturbative treatment of (V.13). It will allow us to identify Navier-Stokes
hydrodynamics appearing in GHD equation supplemented with a collision term.

##  VII Perturbation theory

As discussed in the main text the linearized problem (V.13) is solved by
Fourier transform followed by perturbative solution in $k$. The dynamics of a
single mode is

$\delta q_{n}(x,t)\sim e^{ikx-\Lambda(k)t}\delta q_{n}(k).$ (VII.1)

Substituting this expression to the linearized GHD equations we obtain the
following spectral problem

$\left(\Gamma_{nm}+ikA_{nm}+k^{2}D_{nm}\right)\delta q_{m}(k)=\Lambda(k)\delta
q_{m}(k),$ (VII.2)

where summation over the repeated indices is assumed. We approach the problem
perturbatively assuming that $k$ is small. The spectrum consists then of
gapped and gapless modes. The value of the gap is determined by the smallest
non-zero eigenvalue of $\Gamma$. The gapped modes decay in time and become
irrelevant at sufficiently large time scales. On the other hand, there are
also gapless modes that are related to the collision invariants of $\Gamma$.
In the logic of the perturbation theory in $k$ the spectrum of the gapless
modes has then the following structure

$\Lambda(k)=ak+\omega k^{2}+\mathcal{O}(k^{3}),$ (VII.3)

and importantly the linear coefficients depend only on $A$ and therefore are
universal, whereas the second order coefficients depend on both the collision
integral and the GHD diffusion. The dependence on $\Gamma$ enters through the
second-order perturbation theory from coupling of the zero-modes with the
dispersive ones. On the other hand, the dependence on $D$, comes directly from
the first order perturbation theory because of the $k^{2}$ factor standing in
front of it.

For the purpose of perturbation theory calculations it is convenient to
introduce bra-ket notation such that $\Gamma_{nm}=\langle
h_{n}|\Gamma|h_{m}\rangle$ and similarly for the other operators. More
concretely,

$\langle f|g\rangle=\int{\rm d}\lambda f(\lambda)g(\lambda),\qquad
K|f\rangle=\int{\rm d}\lambda\,K(\cdot,\lambda)f(\lambda),$ (VII.4)

with matrix elements understood in this notation as

$A_{nm}=\langle h_{n}|A|h_{m}\rangle=\langle
h_{n}|\mathcal{A}\mathcal{C}|h_{m}\rangle,\qquad D_{nm}=\langle
h_{n}|D|h_{m}\rangle=\frac{1}{2}\langle
h_{n}|\mathfrak{D}\mathcal{C}|h_{m}\rangle,\qquad\Gamma_{nm}=\langle
h_{n}|\Gamma|h_{m}\rangle.$ (VII.5)

Where $\mathcal{A}=(1-n\mathcal{T})^{-1}v(1-n\mathcal{T})$ [23, 61] and
operators $\mathcal{C},\mathfrak{D}$ are discussed in Sec. I. Operator
$\Gamma$ is discussed in Sec. V with a concrete example in IX. Moreover, for
the basis $h_{n}$ introduced in Sec. V, the orthonormality condition is
$\langle h_{n}|\mathcal{C}|h_{m}\rangle=\delta_{n,m}$.

With this notation, $|h_{i}\rangle$ with $i=0,1,2$ are eigenstates of $\Gamma$
with zero eigenvalues. We denote by $|g_{\beta}\rangle$ the eigenstates with
non-zero eigenvalue $\gamma_{\beta}\neq 0$,
$\Gamma|g_{\beta}\rangle=\gamma_{\beta}|g_{\beta}\rangle$. Whenever summation
over $\beta$ appears, it means summation over eigenvectors with non-zero
eigenvalue. On the other hand, index $\alpha$ always corresponds to thermal
and sound modes.

###  VII.1 Perturbation theory at order $k$

At this order we can neglect the GHD diffusion. The linear part of the
dispersion relation of the hydrodynamic modes is then found by solving the
eigenproblem in the degenerate zero subspace of the collision integral. The
problem reduces to diagonalizing the following matrix

$ikA_{\rm trunc}=ik\begin{pmatrix}0&A_{1,0}&0\\\ A_{1,0}&0&A_{2,1}\\\
0&A_{2,1}&0\end{pmatrix}.$ (VII.6)

We find in total three modes with eigenvalues $\Lambda_{\pm}=\pm iv_{s}k$
(sound modes), and $\Lambda_{\rm th}=0$ (heat mode). The sound velocity given
by the following formula, in agreement with thermodynamic sound velocity, see
Sec. IV,

$v_{s}=\sqrt{A_{1,0}^{2}+A_{2,1}^{2}}=1/\sqrt{\varrho\varkappa_{S}},$ (VII.7)

and associated eigenvectors are (normalization factor is
$\mathcal{N}=\sqrt{1+(A_{2,1}/A_{1,0})^{2}}$)

$|h_{\pm}\rangle=\frac{1}{\sqrt{2}\mathcal{N}}\bigg{(}\frac{A_{1,0}}{A_{2,1}}|h_{0}\rangle\pm\mathcal{N}|h_{1}\rangle+|h_{2}\rangle\bigg{)},\qquad|h_{\rm
th}\rangle=\frac{1}{\mathcal{N}}\bigg{(}-|h_{0}\rangle+\frac{A_{1,0}}{A_{2,1}}|h_{2}\rangle\bigg{)}.$
(VII.8)

Therefore, in the first order of the perturbation theory, the degeneracy has
been completely lifted. It is worth emphasising here that the found
eigenvalues (yielding the speed of sound) are not eigenvalues of the full
matrix $A$. These are eigenvalues of matrix $A_{\rm trunc}$ truncated to
subspace of charges ${q_{0},q_{1},q_{2}}$ conserved by the perturbation
$\Gamma$. Observe that eigenvectors (VII.8) may be written entirely in terms
of the specific heats, using the relations between hydrodynamic matrices and
thermodynamics derived in Sec. IV,

$|h_{\pm}\rangle=\frac{1}{\sqrt{2}}\bigg{(}\sqrt{\frac{c_{V}}{c_{P}}}|h_{0}\rangle\pm|h_{1}\rangle+\sqrt{\frac{c_{P}-c_{V}}{c_{P}}}|h_{2}\rangle\bigg{)},\qquad|h_{\rm
th}\rangle=-\sqrt{\frac{c_{P}-c_{V}}{c_{P}}}|h_{0}\rangle+\sqrt{\frac{c_{V}}{c_{P}}}|h_{2}\rangle.$
(VII.9)

This result can be compared with known formulas for the modes in the classical
three-dimensional ideal gas. For the ideal gas $c_{P}/c_{V}=5/3$ and our
formulas above reproduce results for sound and heat modes, when these values
are substituted, see [54]. Furthermore using (VII.9), we also find useful
relations for the sound modes in terms of the heat mode

$|h_{\pm}\rangle=\sqrt{\frac{c_{P}-c_{V}}{2c_{V}}}|h_{\rm
th}\rangle+\sqrt{\frac{c_{P}}{2c_{V}}}|h_{0}\rangle\pm\frac{1}{\sqrt{2}}|h_{1}\rangle,$
(VII.10)

and heat mode in terms of ultra-local basis functions

$|h_{\rm th}\rangle=\frac{1}{\sqrt{\varrho
T^{2}c_{P}}}\left(|\mathscr{h}_{2}\rangle-\frac{\mathcal{B}_{2,1}}{\mathcal{B}_{1,0}}|\mathscr{h}_{0}\rangle\right).$
(VII.11)

#### Gapless modes are equal to the modes of Navier-Stokes equations

We complete considerations at this level of perturbation theory by noting that
the found heat and sound modes coincide with the modes of Navier-Stokes
equations (VI.10). For the purpose of this part, the bra-ket notation (which
simplifies perturbation theory calculations) is not the most convenient. We
therefore explicitly write sound and heat modes as

$\delta q_{\pm}=\frac{1}{\sqrt{2}}\bigg{(}\sqrt{\frac{c_{V}}{c_{P}}}\delta
q_{0}\pm\delta q_{1}+\sqrt{\frac{c_{P}-c_{V}}{c_{P}}}\delta
q_{2}\bigg{)},\quad\delta q_{\rm th}=-\sqrt{\frac{c_{P}-c_{V}}{c_{P}}}\delta
q_{0}+\sqrt{\frac{c_{V}}{c_{P}}}\delta q_{2}.$ (VII.12)

Already here we can see that these are exactly the modes of Navier-Stokes
equation (VI.10), after the following identification is made

$\delta q_{0}=\delta\tilde{\varrho},\qquad\delta
q_{1}=\delta\tilde{u},\qquad\delta q_{2}=\delta\tilde{T}.$ (VII.13)

We can also observe the truncated matrix (VII.6) expressed with thermodynamic
quantities (IV.29) coincides with the relevant matrix (VI.8) for linearized
Navier-Stokes equations, when terms of order $k^{2}$ are neglected (these
terms do not influence the eigenvectors in $k\to 0$ limit)

$ikA_{\rm trunc}=ik\begin{pmatrix}0&\frac{1}{\sqrt{\varrho\varkappa_{T}}}&0\\\
\frac{1}{\sqrt{\varrho\varkappa_{T}}}&0&\sqrt{\frac{T}{\varrho^{2}c_{V}}}\left(\frac{\partial
P}{\partial T}\right)_{\varrho}\\\
0&\sqrt{\frac{T}{\varrho^{2}c_{V}}}\left(\frac{\partial P}{\partial
T}\right)_{\varrho}&0\end{pmatrix}.$ (VII.14)

The relation (VII.13) can also be justified from the very definition of the
charges $q_{0},q_{1},q_{2}$. To show this, we start with perturbations to
ultra-local charges
$\delta\mathscr{q}_{0},\delta\mathscr{q}_{1},\delta\mathscr{q}_{2}$, see (V.7)
and find

$\delta\mathscr{q}_{0}=\sqrt{\mathcal{C}_{0,0}}\delta
q_{0},\qquad\delta\mathscr{q}_{1}=\sqrt{\mathcal{C}_{1,1}}\delta
q_{1},\qquad\delta\mathscr{q}_{2}=\sqrt{\frac{\mathcal{C}_{2,2}\mathcal{C}_{0,0}-\mathcal{C}_{2,0}^{2}}{\mathcal{C}_{0,0}}}\delta
q_{2}+\frac{\mathcal{C}_{2,0}}{\sqrt{C_{0,0}}}\delta q_{0}.$ (VII.15)

Perturbations of ultra-local charges can be easily connected to perturbations
of the fields $\delta\varrho,\delta u,\delta e$ via relations (II.7) and
(II.21)

$\delta\varrho=\delta\mathscr{q}_{0}\qquad\delta
u=\frac{1}{\varrho}\delta\mathscr{q}_{1}\qquad\delta
e=\frac{1}{\varrho}\delta\mathscr{q}_{2}-\frac{\mathscr{q}_{2}}{\varrho^{2}}\delta\mathscr{q}_{0}.$
(VII.16)

Lastly, to go to $\delta\varrho,\delta u,\delta T$ we need (VI.2). Combining
all three transformations together with thermodynamic formulas from Sec.IV (in
particular, (IV.33)), we find

$\delta\varrho=\sqrt{T\varrho\left(\frac{\partial
P}{\partial\varrho}\right)_{T}^{-1}}\,\delta
q_{0}=\varrho\sqrt{T\varkappa_{T}}\,\delta q_{0},\qquad\delta
u=\sqrt{T/\varrho}\,\delta q_{1},\qquad\delta T=T/\sqrt{\varrho c_{V}}\,\delta
q_{2}.$ (VII.17)

These relations are obviously equivalent to (VII.13). Summing up, we have
shown that perturbations of charges $\delta q_{0},\delta q_{1},\delta_{2}$ are
equal to perturbations of rescaled hydrodynamic fields.

###  VII.2 Perturbation theory at order $k^{2}$

We move now to calculations at the $k^{2}$ order. As already seen in (VI.9),
at this level we should expect transport coefficients $\zeta,\kappa$ appearing
in corrections to the heat and sound modes. The perturbation theory will
provide us with explicit expressions for transport coefficients.

The corrections of the order $k^{2}$ have two origins. There are second-order
corrections coming from the coupling between the hydrodynamic modes and
dispersive modes through the collision integral; we denote them
$\omega_{\alpha}^{\rm coll}$. There are also corrections coming from the GHD
diffusion denoted $\omega_{\alpha}^{\rm diff}$. Here $\alpha=\pm,{\rm th}$
specifies one of the hydrodynamic modes. Thus, the general structure is

$\displaystyle\Lambda_{\rm th}(k)$ $\displaystyle=\left(\omega_{\rm th}^{\rm
diff}+\omega_{\rm th}^{\rm coll}\right)k^{2}+\mathcal{O}(k^{3}),$ (VII.18)
$\displaystyle\Lambda_{\rm\pm}(k)$ $\displaystyle=\pm
iv_{s}k+\left(\omega_{\rm\pm}^{\rm diff}+\omega_{\rm\pm}^{\rm
coll}\right)k^{2}+\mathcal{O}(k^{3}).$

The corrections due to the diffusion take a simple form. They are given by the
first order perturbation theory in the diffusion operator
$\omega_{\alpha}^{\rm diff}=\langle h_{\alpha}|D|h_{\alpha}\rangle$. Because
the diffusion operator is zero between eigenstates of different parity and
$|h_{0}\rangle$ is its zero eigenvector (both left and right) there are only
$2$ non-zero matrix elements $D_{1,1}$ and $D_{2,2}$. In terms of them we find

$\displaystyle\omega_{\rm th}^{\rm diff}$
$\displaystyle=\frac{c_{V}}{c_{P}}D_{2,2},$ (VII.19)
$\displaystyle\omega_{\pm}^{\rm diff}$
$\displaystyle=\frac{1}{2}D_{1,1}+\frac{1}{2}\frac{c_{P}-c_{V}}{c_{P}}D_{2,2}.$

Using the decomposition (VII.10) we also find

$\omega_{\pm}^{\rm diff}=\frac{c_{P}-c_{V}}{2c_{V}}\omega_{\rm th}^{\rm
diff}+\frac{1}{2}D_{1,1}.$ (VII.20)

The second-order correction to the eigenvalues of the hydrodynamic modes come
from the standard formula from the perturbation theory

$\omega_{\alpha}^{\rm coll}=\sum_{\beta}\frac{\langle
h_{\alpha}|A|g_{\beta}\rangle\langle
g_{\beta}|A|h_{\alpha}\rangle}{\gamma_{\beta}},\qquad\alpha=\pm,{\rm th},$
(VII.21)

where the sum does not include the kernel of operator $\Gamma$. Inspecting the
formula we find $\omega_{+}^{\rm coll}=\omega_{-}^{\rm coll}$. This follows
from the fact, that $\Gamma$ operator has eigenfunctions which are either
symmetric or antisymmetric functions of rapidities.

Note that action of $A$ on $|h_{0}\rangle$ gives a mode proportional to
$|h_{1}\rangle$ (this can be verified directly at the level of matrix elements
(V.14)) and therefore the middle term from the decomposition (VII.10) does not
contribute to $\omega_{\pm}^{\rm coll}$ because $\langle
g_{\beta}|h_{1}\rangle=0$ from the definition of $|g_{\beta}\rangle$. In the
expression for $\omega_{\pm}^{\rm coll}$ there appear also mixed terms of the
form $\langle h_{\rm th}|A|g_{\beta}\rangle\langle g_{\beta}|A|h_{1}\rangle$.
Such terms vanish because $A$ is non-zero only for states of different parity.
This allows us to write

$\omega_{\pm}^{\rm coll}=\frac{c_{P}-c_{V}}{2c_{V}}\omega_{{\rm th}}^{\rm
coll}+\frac{1}{2}\omega_{1}^{\rm coll},\qquad\omega_{1}^{\rm
coll}=\sum_{\beta}\frac{(\langle
g_{\beta}|A|h_{1}\rangle)^{2}}{\gamma_{\beta}}.$ (VII.22)

Figure 1: We present numerical solutions of (VII.2) for the coupled Lieb-
Liniger model analyzed in the main text. The presented results are for $T=2$,
$n=1$ and $c=2$. Left: small-$k$ part of the spectrum of dispersion relations.
For sufficiently small $k$, numerical results (points) agree with universal
Navier-Stokes dispersion relations (VI.9) with transport coefficients (VII.23)
(lines). At larger momenta we see deviations caused by the higher order
corrections in $k$. We have multiplied the real parts of the dispersions by
factor $10$ to make them of the same order as imaginary parts. Right: a
numerical solution requires truncating the infinite system of equations to a
finite one. Truncating to first $12$ (squares), $18$ (diamonds) and $24$
(solid lines) charges yields results which differ on a scale below the
resolution of the plot.

###  VII.3 Transport coefficients

We can finally read off the transport coefficients by comparing spectra (VI.9)
of the linearized Navier-Stokes with the results of the perturbation theory.
We find

$\displaystyle\kappa^{\rm coll}$ $\displaystyle=\varrho c_{P}\omega_{\rm
th}^{\rm coll},\qquad\zeta^{\rm coll}=\varrho\omega_{1}^{\rm coll},$ (VII.23)
$\displaystyle\kappa^{\rm diff}$ $\displaystyle=\varrho
c_{V}D_{2,2},\qquad\zeta^{\rm diff}=\varrho D_{1,1}.$

with $\kappa=\kappa^{\rm coll}+\kappa^{\rm diff}$ and $\zeta=\zeta^{\rm
coll}+\zeta^{\rm diff}$.

We confirm these formulas on Fig.1, which presents small-$k$ part of the
spectrum found in numerical diagonalization of (VII.2). Our results agree with
numerics in the samll-$k$ regime, as expected. The dispersion relations for
larger number of modes and for higher values of $k$ are presented in Fig.2.
For small $k$ the spectrum consists of Navier-Stokes modes and gapped modes.
For high enough $k$, the spectrum is those of unperturbed GHD. This
corresponds to operator $\Gamma$ being negligible in the eigenproblem (VII.2),
as compared to $A$ and $D$. In between, there is a region which cannot be
described neither by Navier-Stokes hydrodynamics, nor the pure diffusive GHD.
Therefore, by numerically solving (VII.2) we are able to quantatively capture
the crossover from Navier-Stokes hydrodynamics to GHD.

Figure 2: The real part of the dispersion relations for the lowest lying modes
from exact numerical solution of (VII.2). At small momenta, the spectrum is
characterized by the $3$ gapless modes (the orange curve is doubly degenerate)
of the Navier-Stokes regime whereas at (relatively) large momenta by the modes
of the GHD. The dashed lines are the gapless modes of the GHD without the
integrability breaking terms. The joining and splitting of the modes defines
an intermediate regime (shaded) where the dynamics cannot be reduced to either
Navier-Stokes or unperturbed GHD.

Equations (VII.23) together with (VII.22) provide also an efficient way to
evaluate the transport coefficients numerically. The procedure is the
following. For a given thermal equilibrum state, we construct an orthonormal
basis $\\{|h_{i}\rangle\\}$ using the hydrodynamic scalar product $(f,g)$.
This allows us to compute the matrix elements $A_{nm}$, $\Gamma_{mn}$ and
$D_{nm}$. The matrix $\Gamma_{mn}$ can be now diagonalized numerically. We
obtain the spectrum $\gamma_{\beta}$ but also expressions for
$|g_{\beta}\rangle$ in terms of linear combinations of $|h_{i}\rangle$’s. This
allows us to evaluate the matrix elements $\langle
g_{\beta}|A|h_{\alpha}\rangle$ and $\omega_{\alpha}^{\rm coll}$. The remaining
ingredient is the heat capacity $c_{P}$ which can be computed from the
hydrodynamic matrices $\mathcal{B}$ and $\mathcal{C}$ according to (IV.27).
The numerical solution requires truncating the infinite basis to a finite one.
In Fig. 1 we show that it is sufficient to take around $20$ basis functions.

Let us note that contributions $\kappa^{\rm diff}$ and $\zeta^{\rm diff}$ can
be also written as matrix elements in the ultra-local basis, as in the main
text. Using relations (V.7), properties of the diffusion kernel and adopting
notation of matrix elements of diffusion operator (V.14) as in the main text
we find that indeed

$\kappa^{\rm
diff}=\frac{1}{2T^{2}}\mathscr{h}_{2}\mathfrak{D}\mathcal{C}\mathscr{h}_{2},\qquad\zeta^{\rm
diff}=\frac{1}{2T}\mathscr{h}_{1}\mathfrak{D}\mathcal{C}\mathscr{h}_{1}.$
(VII.24)

Importantly $\kappa^{\rm diff},\zeta^{\rm diff}\geq 0$, due to the positivity
of operator $\mathfrak{D}\mathcal{C}$ [24]. Moreover, it is easy to observe
that $\omega_{1}^{\rm coll},\omega^{\rm coll}_{\rm th}\geq 0$ and hence
$\kappa^{\rm coll},\zeta^{\rm coll}\geq 0$. We thus conclude that transport
coefficients found by us are non-negative, which assures that entropy
production in the Navier-Stokes dynamics is positive [54].

In the next section we will confirm that $\kappa^{\rm coll}$ and $\zeta^{\rm
coll}$ found in the second-order perturbation theory (VII.23) can be written
as generalized Chapman-Enskog expressions given in the main text.

##  VIII Integral equations for $\kappa^{\rm coll}$ and $\zeta^{\rm coll}$
transport coefficients

Expressions for the collisional transport coefficients involve sums over all
non-zero eigenstates of the $\Gamma$ operator. In the following, we will show
how they can be rewritten as integral equations, which are identical to the
integral equations found in the Chapman-Enskog procedure. We start with the
thermal conductivity related to $\omega_{\rm th}^{\rm coll}$. As mentioned
earlier for any $\beta$ we have $\langle h_{0}|A|g_{\beta}\rangle=0$. Using
explicit form (VII.11) of $|h_{\rm th}\rangle$ in ultra-local basis we find

$\langle h_{\rm th}|A|g_{\beta}\rangle=\frac{1}{\sqrt{\varrho
T^{2}c_{P}}}\langle\mathscr{h}_{2}|A|g_{\beta}\rangle,$ (VIII.1)

and thus we have

$\omega_{\rm th}^{\rm coll}=\frac{1}{\sqrt{\varrho
T^{2}c_{P}}}\sum_{\beta}\frac{\langle\mathscr{h}_{2}|A|g_{\beta}\rangle\langle
g_{\beta}|A|h_{\rm th}\rangle}{\gamma_{\beta}}.$ (VIII.2)

For the heat mode the sum can be easily extended to the whole space by noting
that $A|h_{\rm th}\rangle$ is orthogonal to the kernel of $\Gamma$. For
$|h_{0}\rangle,|h_{2}\rangle$ it is clear from symmetry properties of $A$. For
$|h_{1}\rangle$ we can show $\langle h_{1}|A|h_{\rm th}\rangle=0$, which be
seen from explicit computation using (VII.11). This allows us to write

$\omega_{\rm th}^{\rm coll}=\lim_{\epsilon\to 0}\frac{1}{\sqrt{\varrho
T^{2}c_{P}}}\langle\mathscr{h}_{2}|A\Gamma^{-1}_{\epsilon}A|h_{\rm
th}\rangle,$ (VIII.3)

where we used the completeness of the basis which implies that

$\lim_{\varepsilon\to
0}\sum_{m}|g_{m}\rangle\frac{1}{\gamma_{m}-\varepsilon}\langle
g_{m}|=\lim_{\varepsilon\to
0}\frac{1}{\Gamma-\varepsilon\mathcal{C}}\equiv\lim_{\epsilon\to
0}\Gamma^{-1}_{\epsilon},$ (VIII.4)

where $m$ extends now over all the eigenstates. Now, we can formally introduce
vector $|\phi_{\kappa}\rangle$ such that,

$\Gamma_{\epsilon}^{-1}A|h_{\rm th}\rangle=\frac{T}{\sqrt{\varrho
c_{P}}}|\phi_{\kappa}\rangle,\qquad\omega_{\rm th}^{\rm coll}=\frac{1}{\varrho
c_{P}}\langle\mathscr{h}_{2}|A|\phi_{\kappa}\rangle.$ (VIII.5)

To find $|\phi_{\kappa}\rangle$ we can go to quasiparticle representation. We
look for a function $\eta_{\kappa}(\lambda)$ which fulfills the following
integral equation

$\Gamma_{\epsilon}\phi_{\kappa}=(1-n\mathcal{T})^{-1}\rho_{\rm
p}f\eta_{\kappa}$ (VIII.6)

with function $\eta_{\kappa}(\lambda)$ related to heat mode and given by

$\eta_{\kappa}(\lambda)=\frac{\sqrt{\varrho c_{P}}}{T}v(\lambda)h_{\rm
th}^{\rm dr}(\lambda)=\frac{\sqrt{\varrho
c_{V}}}{T}\left(h_{2}-\sqrt{\frac{c_{P}-c_{V}}{c_{V}}}h_{0}\right)^{\rm
dr}(\lambda).$ (VIII.7)

For $\epsilon\neq 0$, the operator $\Gamma_{\epsilon}$ is invertible and the
solution to (VIII.6) is unique. However, we would like to explicitly take the
limit $\epsilon\to 0$, which will invalidate uniqueness of the solution. This
is clear from the existence of collision invariants of $\Gamma$ \- for any
particular solution $\bar{\phi}_{\kappa}$ to the equation (VIII.6) with
$\epsilon=0$ the function

$\bar{\phi}_{\kappa}+\sum_{i=0}^{2}a_{i}\mathscr{h}_{i}$ (VIII.8)

is an another solution (coefficients $a_{i}$ are arbitrary). Therefore, before
we set $\epsilon=0$, let us multiply (VIII.6) with collision invariants
$\mathscr{h}_{i}(\lambda)$ where $i=0,1,2$ and integrate over
$\differential\lambda$. This leads to the conditions

$\langle\mathscr{h}_{i}|\mathcal{C}|\phi_{\kappa}\rangle=(\mathscr{h}_{i}|\phi_{\kappa})=0,$
(VIII.9)

which we keep after taking $\epsilon\to 0$. The conditions above make the
solution unique by fixing the constants $a_{i}$ in the general solution.

In the main text we shorten the notation introducing
$\chi_{\kappa}=(1-n\mathcal{T})^{-1}\rho_{\rm p}f\eta_{\kappa}$, which
explicitly reads

$\chi_{\kappa}=\frac{\sqrt{\varrho
c_{V}}}{T}\mathcal{B}(h_{2}-\sqrt{c_{P}/c_{V}-1}\,h_{0}).$ (VIII.10)

The integral equation can be then written as
$\Gamma\phi_{\kappa}=\chi_{\kappa}$, as in the main text. Finally, thermal
conductivity $\kappa^{\rm coll}$ is given by

$\kappa^{\rm coll}=\int\differential\lambda\rho_{\rm
p}(\lambda)f(\lambda)v(\lambda)\mathscr{h}_{2}^{\rm
dr}(\lambda)\phi_{\kappa}^{\rm
dr}(\lambda)=\mathscr{h}_{2}\mathcal{B}\phi_{\kappa}.$ (VIII.11)

We focus now on $\zeta^{\rm coll}$. We proceed similarly, with the exception
that now $\langle g_{\beta}|A|h_{1}\rangle\neq 0$, therefore we have to
properly regularize matrix elements before we introduce representation of
inverse $\Gamma$. We start with expression

$\langle g_{\beta}|A|h_{1}\rangle,$ (VIII.12)

which enters the sum for $\zeta^{\rm coll}$. We can add to it two zero terms,
namely:

$\langle g_{\beta}|A|h_{1}\rangle+a\langle g_{\beta}|h_{0}\rangle+b\langle
g_{\beta}|h_{2}\rangle.$ (VIII.13)

We want $a,b$ to be such that the vector
$|w_{1}\rangle=A|h_{1}\rangle+a|h_{0}\rangle+b|h_{2}\rangle$ is orthogonal to
$|h_{0}\rangle,|h_{1}\rangle,|h_{2}\rangle$. We find $a=-A_{1,0},b=-A_{2,1}.$
Thus, we write (we have used here that $|h_{1}\rangle=1/\sqrt{\varrho
T}|\mathscr{h}_{1}\rangle$, see (V.7) and Sec.IV)

$\omega_{1}^{\rm coll}=\frac{1}{\sqrt{\varrho T}}\lim_{\epsilon\rightarrow
0}\langle\mathscr{h}_{1}|A\Gamma_{\epsilon}^{-1}|w_{1}\rangle.$ (VIII.14)

As previously, we introduce a vector $|\phi_{\zeta}\rangle$ such that
$\Gamma_{\epsilon}^{-1}|w_{1}\rangle=\sqrt{T/\varrho}|\phi_{\zeta}\rangle$. In
quasiparticle representation finding such vector corresponds to solving the
following integral equation

$\Gamma_{\epsilon}\phi_{\zeta}=(1-n\mathcal{T})^{-1}\rho_{\rm
p}f\eta_{\zeta},$ (VIII.15)

with function (we have used here (IV.29))

$\eta_{\zeta}(\lambda)=\sqrt{\frac{\varrho}{T}}\left[v(\lambda)h_{1}^{\rm
dr}(\lambda)-\frac{1}{\sqrt{\varrho\varkappa_{T}}}\left(h_{0}+\sqrt{\frac{c_{P}-c_{V}}{c_{V}}}h_{2}\right)^{\rm
dr}(\lambda)\right].$ (VIII.16)

Similarly as for thermal conductivity we set $\epsilon\to 0$ by imposing
additional conditions $(\mathscr{h}_{i}|\phi_{\zeta})=0,\,\,i=0,1,2$. We
shorten the notation and introduce
$\chi_{\zeta}=(1-n\mathcal{T})^{-1}\rho_{\rm p}f\eta_{\zeta}$, which
explicitly reads

$\chi_{\zeta}=\sqrt{\frac{\varrho}{T}}\mathcal{B}h_{1}-\frac{1}{\sqrt{T\varkappa_{T}}}\,\mathcal{C}(h_{0}+\sqrt{c_{P}/c_{V}-1}\,h_{2}).$
(VIII.17)

The integral equation is then $\Gamma\phi_{\zeta}=\chi_{\zeta}$ and the bulk
viscosity $\zeta^{\rm coll}$ finally is

$\zeta^{\rm coll}=\int\differential\lambda\rho_{\rm
p}(\lambda)f(\lambda)v(\lambda)\mathscr{h}_{1}^{\rm
dr}(\lambda)\phi_{\zeta}^{\rm
dr}(\lambda)=\mathscr{h}_{1}\mathcal{B}\phi_{\zeta}.$ (VIII.18)

In this way, we have shown that results from perturbation theory can be
reformulated as generalized Chapman-Enskog integral equations, given in the
main text. One of the advantages of integral equation approach is simple
analysis of the non-interacting limit of the underlying integrable model. This
is done in the next section.

### Non-interacting limit of transport coefficients

In this section we look more closely on the limit of non-interacting
integrable model with trivial dressings. Transport coefficients (VII.23) have
contributions from collision integral and from diffusion operator. In the non-
interacting limit, diffusion operator vanishes together with the associated
contributions to the transport coefficients. We thus turn to contributions
from collision integral, which in principle may be non-zero even in the limit
of free integrable model. Obviously, we still consider non-zero integrability
breaking term $\mathcal{I}$.

In the case of trivial dressings, the following symmetries appear on the level
of hydrodynamic matrices

$\mathcal{B}_{1,0}=2C_{2,0},\qquad\mathcal{B}_{2,1}=2\mathcal{C}_{2,2}.$
(VIII.19)

Let us look how these affect expressions for $\eta_{\kappa,\zeta}(\lambda)$
given by (VIII.7), (VIII.16). It is convenient to represent these functions in
the ultra-local basis first, we find

$\displaystyle\eta_{\kappa}(\lambda)=$
$\displaystyle\frac{1}{T^{2}}v(\lambda)\left(\mathscr{h}_{2}^{\rm
dr}(\lambda)-\frac{\mathcal{B}_{2,1}}{\mathcal{B}_{1,0}}\mathscr{h}_{0}^{\rm
dr}(\lambda)\right),$ (VIII.20) $\displaystyle\eta_{\zeta}(\lambda)=$
$\displaystyle\frac{1}{T}\left(v(\lambda)\mathscr{h}^{\rm
dr}_{1}(\lambda)+\frac{\mathcal{B}_{2,1}\mathcal{C}_{2,0}-\mathcal{B}_{1,0}\mathcal{C}_{2,2}}{\mathcal{C}_{2,2}\mathcal{C}_{0,0}-\mathcal{C}_{2,0}^{2}}\mathscr{h}_{0}^{\rm
dr}(\lambda)-\frac{\mathcal{B}_{2,1}\mathcal{C}_{0,0}-\mathcal{B}_{1,0}\mathcal{C}_{2,0}}{\mathcal{C}_{2,2}\mathcal{C}_{0,0}-\mathcal{C}_{2,0}^{2}}\mathscr{h}_{2}^{\rm
dr}(\lambda)\right),$ (VIII.21)

where we have expressed thermodynamic quantities with hydrodynamic matrices,
see Sec. IV. Moreover, we have represented GS-orthonormalized vectors (V.7)
with ultra-local counterparts. Let us consider a limit, where dressing can be
neglected and hence $v(\lambda)=\mathscr{h}_{1}(\lambda)$. Taking into account
also the symmetries (VIII.19), we observe

$\displaystyle\eta_{\kappa}(\lambda)\to$
$\displaystyle\,\frac{1}{T^{2}}\mathscr{h}_{1}(\lambda)\left(\mathscr{h}_{2}(\lambda)-\frac{\mathcal{B}_{2,1}}{\mathcal{B}_{1,0}}\mathscr{h}_{0}(\lambda)\right),$
(VIII.22) $\displaystyle\eta_{\zeta}(\lambda)\to$ $\displaystyle\,0.$
(VIII.23)

The rate on which $\eta_{\zeta}(\lambda)$ approaches zero with vanishing
coupling depends on the model. This might be important in cases where $\Gamma$
vanishes in the free limit as well. When $\eta_{\zeta}(\lambda)$ approaches
zero faster than $\Gamma$ with the interaction coupling, the bulk viscosity
tends to 0. This is the case for collision integral studied by us. On the
other hand, as $\eta_{\kappa}(\lambda)$ remains finite in the non-interacting
limit, the $\Gamma$ operator, which vanishes in the free limit leads then to
infinite thermal conductivity. This is a known result for non-interacting gas
[54] and we observe it also for collision integral considered by us (see Fig.
2 in the main text).

##  IX Collision integral for the coupled Lieb-Liniger models

In this supplementary material we derive the linearized collision integral for
the two coupled Lieb-Liniger models. It has all the necessary properties
which, we have assumed about the $\Gamma$ operator and will allow us to
illustrate our general formalism in an experimentally relevant situation.

Generally speaking there are two contributions to the collision integral
$\mathcal{I}$. The direct contributions comes from the very scattering events
caused by the perturbing potential. But there is also an indirect
contribution, that arises because a modification to the distribution of
rapidities in one place (in the rapidity space) reshuffles rapidities
elsewhere due to the interaction present in the Lieb-Liniger model. The two
contributions together appear as a dressing of the bare scattering integral
$\mathcal{I}_{0}$

$\mathcal{I}[\rho_{\rm{p},1},\rho_{\rm{p},2}](\lambda)=\int\differential\mu\mathbb{F}(\lambda,\mu)\mathcal{I}_{0}[\rho_{\rm{p},1},\rho_{\rm{p},2}](\mu),$
(IX.1)

where $\mathbb{F}$ is backflow operator introduced in (V.18).

As discussed in [29, 44] in the homogeneous case when two tubes are in the
same state the leading contribution to the collision integral comes from
processes of type $(1,2)$ and $(2,1)$. These are processes in which $1$
particle-hole pair is created in the first tube and $2$ particle-hole pair is
created in the second tube or vice verse,
$\mathcal{I}[\rho_{\rm{p}}]=\mathcal{I}^{(1,2)}[\rho_{\rm{p}}]+\mathcal{I}^{(2,1)}[\rho_{\rm{p}}]+(\dots)$,
where $(\dots)$ stands for processes of higher order in momentum transferred
between the tubes.

The bare collision integral for these processes is [29, 44]

$\displaystyle\mathcal{I}_{0}[\rho_{\rm p}](\lambda)=$
$\displaystyle\frac{(2\pi)^{2}}{(2!)^{2}}\int\differential
p_{1,1}\differential h_{1,1}\differential p_{1,2}\differential
h_{1,2}\differential p_{2}\differential h_{2}\,\left(\delta(\lambda-
p_{1,1})+\delta(\lambda-p_{1,2})+\delta(\lambda-
p_{2})\right)\,\tilde{A}^{2}\big{(}k_{2\rm
ph}(\mathbf{p}_{1},\mathbf{h}_{1})\big{)}\times$ (IX.2) $\displaystyle|F_{2\rm
ph}(\mathbf{p}_{1},\mathbf{h}_{1})|^{2}|F_{1\rm
ph}(\mathbf{p}_{2},\mathbf{h}_{2})|^{2}\,\delta\big{(}k_{2\rm
ph}(\mathbf{p}_{1},\mathbf{h}_{1})+k_{1\rm
ph}(\mathbf{p}_{2},\mathbf{h}_{2})\big{)}\,\delta\big{(}\omega_{2\rm
ph}(\mathbf{p}_{1},\mathbf{h}_{1})+\omega_{1\rm
ph}(\mathbf{p}_{2},\mathbf{h}_{2})\big{)}\times$ $\displaystyle\Big{(}J_{2\rm
ph}(\mathbf{p}_{1},\mathbf{h}_{1})\,J_{1\rm
ph}(\mathbf{p}_{2},\mathbf{h}_{2})-(\mathbf{h}\leftrightarrow\mathbf{p})\Big{)},$

where $\tilde{A}(k)$ is Fourier transform of the interaction potential with
$\tilde{A}(k)=\tilde{A}(-k)$, $F_{n\rm ph}(\mathbf{p},\mathbf{h})$ are density
operator form factors and

$k_{n\rm
ph}(\mathbf{p},\mathbf{h})=\sum_{i=1}^{n}\left(k(p_{i})-k(h_{i})\right)\qquad\omega_{n\rm
ph}(\mathbf{p},\mathbf{h})=\sum_{i=1}^{n}\left(\omega(p_{i})-\omega(h_{i})\right).$
(IX.3)

where $k(\lambda)=(\lambda)^{\rm Dr}$ and $\omega(\lambda)=(\lambda^{2})^{\rm
Dr}$. Finally, $J$ factors are

$J_{n\rm ph}(\mathbf{p},\mathbf{h})=\prod_{i=1}^{n}\rho_{\rm
p}(h_{i})\rho_{\rm h}(p_{i}).$ (IX.4)

We linearize now the collision integral around a thermal state. Any thermal
state is a stationary state fulfilling the condition

$J_{2\rm ph}(\mathbf{p}_{1},\mathbf{h}_{1})\,J_{1\rm
ph}(\mathbf{p}_{2},\mathbf{h}_{2})-J_{2\rm
ph}(\mathbf{h}_{1},\mathbf{p}_{1})\,J_{1\rm
ph}(\mathbf{h}_{2},\mathbf{p}_{2})=0,$ (IX.5)

on a $4$-dimensional manifold given by $\omega_{2\rm
ph}(\mathbf{p}_{1},\mathbf{h}_{1})+\omega_{1\rm
ph}(\mathbf{p}_{2},\mathbf{h}_{2})=0$ and $k_{2\rm
ph}(\mathbf{p}_{1},\mathbf{h}_{1})+k_{1\rm
ph}(\mathbf{p}_{2},\mathbf{h}_{2})=0$. To expand the collision integral around
the thermal state it is useful to switch from the distribution functions to
the pseudoenergy $\epsilon(\lambda)$ for which we have

$\rho_{\rm p}(\lambda)=\frac{\rho_{\rm
tot}(\lambda)}{1+e^{\epsilon(\lambda)}},\qquad\rho_{\rm
h}(\lambda)=\frac{\rho_{\rm tot}(\lambda)}{1+e^{-\epsilon(\lambda)}}.$ (IX.6)

The ratio of the $J$-factors has then a simple form

$\frac{J_{2\rm ph}(\mathbf{p}_{1},\mathbf{h}_{1})\,J_{1\rm
ph}(\mathbf{p}_{2},\mathbf{h}_{2})}{J_{2\rm
ph}(\mathbf{h}_{1},\mathbf{p}_{1})\,J_{1\rm
ph}(\mathbf{h}_{2},\mathbf{p}_{2})}=\exp\left(\epsilon(p_{1,1})-\epsilon(h_{1,1})+\epsilon(p_{1,2})-\epsilon(h_{1,2})+\epsilon(p_{2,1})-\epsilon(h_{2,1})\right),$
(IX.7)

and the argument of the exponential vanishes for thermal states for which
$\epsilon(\lambda)=\beta_{2}\omega(\lambda)+\beta_{0}$.

To linearize the collision integral we consider a deviation from the thermal
state $\epsilon\to\epsilon+\delta\epsilon$. The leading correction arrises
from the $J$-factors, the other terms factors contribute at higher orders in
$\delta\epsilon$. The result is

$\mathcal{I}_{0}[\rho]=\Gamma_{0}[\rho_{\rm
th}]\delta\epsilon+\mathcal{O}((\delta\epsilon)^{2}),$ (IX.8)

where

$\displaystyle\Gamma_{0}(\lambda,\mu)=$
$\displaystyle\frac{(2\pi)^{2}}{8}\int\differential p_{1,1}\differential
h_{1,1}\differential p_{1,2}\differential h_{1,2}\differential
p_{2}\differential
h_{2}\,G(\mathbf{p}_{1},\mathbf{h}_{1},\mathbf{p}_{2},\mathbf{h}_{2})\times$
(IX.9) $\displaystyle\Big{(}\delta(\lambda-p_{1,1})-\delta(\lambda-
h_{1,1})+\delta(\lambda-p_{1,2})-\delta(\lambda-h_{1,2})+\delta(\lambda-
p_{2})-\delta(\lambda-h_{2})\Big{)}\times$ $\displaystyle\Big{(}\delta(\mu-
p_{1,1})-\delta(\mu-h_{1,1})+\delta(\mu-p_{1,2})-\delta(\mu-
h_{1,2})+\delta(\mu-p_{2})-\delta(\mu-h_{2})\Big{)},$

with

$\displaystyle
G(\mathbf{p}_{1},\mathbf{h}_{1},\mathbf{p}_{2},\mathbf{h}_{2})=$
$\displaystyle\tilde{A}^{2}\big{(}k_{2\rm
ph}(\mathbf{p}_{1},\mathbf{h}_{1})\big{)}|F_{2\rm
ph}(\mathbf{p}_{1},\mathbf{h}_{1})|^{2}|F_{1\rm
ph}(\mathbf{p}_{2},\mathbf{h}_{2})|^{2}J_{2\rm
ph}(\mathbf{p}_{1},\mathbf{h}_{1})\,J_{1\rm
ph}(\mathbf{p}_{2},\mathbf{h}_{2})\times$ $\displaystyle\delta\big{(}k_{2\rm
ph}(\mathbf{p}_{1},\mathbf{h}_{1})+k_{1\rm
ph}(\mathbf{p}_{2},\mathbf{h}_{2})\big{)}\,\delta\big{(}\omega_{2\rm
ph}(\mathbf{p}_{1},\mathbf{h}_{1})+\omega_{1\rm
ph}(\mathbf{p}_{2},\mathbf{h}_{2})\big{)}.$ (IX.10)

Note that
$G(\mathbf{p}_{1},\mathbf{h}_{1},\mathbf{p}_{2},\mathbf{h}_{2})=G(\mathbf{h}_{1},\mathbf{p}_{1},\mathbf{h}_{2},\mathbf{p}_{2})$,
namely function $G$ is $\rm ph$-symmetric. This implies operator $\Gamma_{0}$
is symmetric, $\Gamma_{0}(\lambda,\mu)=\Gamma_{0}(\mu,\lambda)$. Moreover,
$\Gamma_{0}$ has an additional symmetry. By performing change of variables in
the integrals (changing sign in each one) we get
$\Gamma_{0}(\lambda,\mu)=\Gamma_{0}(-\lambda,-\mu)$. This property leads to
certain selection rules once we consider the action of $\Gamma$ in the basis
of ultra-local conserved charges.

The linearized collision integral has also three left eigenvectors with zero
eigenvalues [44]. They reflect the fact that the total number of particles,
the total momentum and the total energy are conserved by the perturbation,

$\int{\rm d}\lambda(1)^{\rm Dr}(\lambda)\Gamma_{0}(\lambda,\mu)=\int{\rm
d}\lambda(\lambda)^{\rm Dr}(\lambda)\Gamma_{0}(\lambda,\mu)=\int{\rm
d}\lambda(\lambda^{2})^{\rm Dr}(\lambda)\Gamma_{0}(\lambda,\mu)=0,$ (IX.11)

where to derive the last two properties we used conservation laws under the
integral.

We consider now matrix elements of $\Gamma$ in the basis of the ultra-local
conserved charges. The matrix elements are defined according to eq.(V.14),

$\Gamma_{nm}=\int\differential\lambda\differential\mu(h_{n})^{\rm
Dr}(\lambda)\Gamma_{0}(\lambda,\mu)(h_{m})^{\rm dr}(\mu).$ (IX.12)

The dressing on $h_{m}$ appears because $\Gamma_{0}$ acts naturally on
$\delta\epsilon$ (IX.8) and not on $\delta\epsilon_{0}$, which in the end is
expanded, see Sec.V. Therefore, we have to use relation
$\delta\epsilon=\delta\epsilon_{0}^{\rm dr}$, discussed in Sec.V. Moreover, to
make connection with notation in the main text let us note that, the $\Gamma$
operator, which appears in generalized Chapman-Enskog equations can be written
as $\Gamma=\mathbb{F}\,\Gamma_{0}(1-\mathcal{T}n)^{-1}$.

From the construction of our basis (see Sec.V), functions $h_{n}$ are either
symmetric or antisymmetric (and both dressings do not change this property).
This implies that $\Gamma_{nm}=0$ whenever functions $h_{n}$ and $h_{m}$ have
different parity. According to the above discussion on the zero eigenvectors
we also have $\Gamma_{nm}=0$ for $n=0,1,2$. Let us note here, that when
dressing is neglected (large coupling expansion), one finds that
$\Gamma_{nm}=0$ when $n=0,1,2$ but also when $m=0,1,2$.

### Small momentum limit

The matrix elements $\Gamma_{nm}$ read explicitly

$\displaystyle\Gamma_{nm}=$
$\displaystyle\frac{(2\pi)^{2}}{8}\int\differential p_{1,1}\differential
h_{1,1}\differential p_{1,2}\differential h_{1,2}\differential
p_{2}\differential
h_{2}\,G(\mathbf{p}_{1},\mathbf{h}_{1},\mathbf{p}_{2},\mathbf{h}_{2})\times$
(IX.13) $\displaystyle\Big{(}h_{n}^{\rm Dr}(p_{1,1})-h_{n}^{\rm
Dr}(h_{1,1})+h_{n}^{\rm Dr}(p_{1,2})-h_{n}^{\rm Dr}(h_{1,2})+h_{n}^{\rm
Dr}(p_{2})-h_{n}^{\rm Dr}(h_{2})\Big{)}\times$ $\displaystyle\Big{(}h_{m}^{\rm
dr}(p_{1,1})-h_{m}^{\rm dr}(h_{1,1})+h_{m}^{\rm dr}(p_{1,2})-h_{m}^{\rm
dr}(h_{1,2})+h_{m}^{\rm dr}(p_{2})-h_{m}^{\rm dr}(h_{2})\Big{)}.$

We will find now a much simpler expression for them in the small momentum
limit. To this end we express the integrand in terms of center of mass
rapidities and the momentum transferred between the tubes.

First, we change the integration variables to the center-of-mass rapidities.
These are defined through $\lambda=(p+h)/2$ and $\alpha=p-h$ for each
particle-hole excitation. The Jacobian of the transformation from $(p,h)$ to
$(\lambda,\alpha)$ is $1$ hence this amounts to simply replacing the
integration measure by
$\differential\lambda_{1,1}\differential\lambda_{1,2}\differential\lambda_{2}\differential\alpha_{1,1}\differential\alpha_{1,2}\differential\alpha_{2}$.
We then assume that the main contribution to the integral comes from the small
momentum excitations. Therefore, the integrand can be expanded in small
$\alpha$’s. Specifically, the energy-momentum constraints become

$k^{\prime}(\lambda_{1,1})\alpha_{1,1}+k^{\prime}(\lambda_{1,2})\alpha_{1,2}+k^{\prime}(\lambda_{2})\alpha_{2}=0,$
(IX.14a)
$\omega^{\prime}(\lambda_{1,1})\alpha_{1,1}+\omega^{\prime}(\lambda_{1,2})\alpha_{1,2}+\omega^{\prime}(\lambda_{2})\alpha_{2}=0.$
(IX.14b)

Solving them for $\alpha_{1,1}$ and $\alpha_{1,2}$ gives

$\displaystyle\bar{\alpha}_{1,1}=-\alpha_{2}\frac{k^{\prime}(\lambda_{2})}{k^{\prime}(\lambda_{1,1})}\frac{v(\lambda_{2})-v(\lambda_{1,2})}{v(\lambda_{1,1})-v(\lambda_{1,2})},\qquad\bar{\alpha}_{1,2}=-\alpha_{2}\frac{k^{\prime}(\lambda_{2})}{k^{\prime}(\lambda_{1,2})}\frac{v(\lambda_{2})-v(\lambda_{1,1})}{v(\lambda_{1,2})-v(\lambda_{1,1})}.$
(IX.15)

We can now use the following property of $\delta$-functions

$\delta(f(x,y))\delta(g(x,y))=\frac{\delta(x-\bar{x})\delta(y-\bar{y})}{|\det
J|},\qquad J=\begin{pmatrix}\partial_{x}f(x,y)&\partial_{y}f(x,y)\\\
\partial_{x}g(x,y)&\partial_{y}g(x,y)\end{pmatrix},$ (IX.16)

where $(\bar{x},\bar{y})$ is a zero of $(f(x,y),g(x,y))$. This gives for the
$\delta$-functions implementing the conservation of energy-momentum the
following expression

$\delta\big{(}k_{2\rm ph}(\mathbf{p}_{1},\mathbf{h}_{1})+k_{1\rm
ph}(\mathbf{p}_{2},\mathbf{h}_{2})\big{)}\,\delta\big{(}\omega_{2\rm
ph}(\mathbf{p}_{1},\mathbf{h}_{1})+\omega_{1\rm
ph}(\mathbf{p}_{2},\mathbf{h}_{2})\big{)}=\frac{\delta(\alpha_{1,1}-\bar{\alpha}_{1,1})\delta(\alpha_{1,2}-\bar{\alpha}_{1,2})}{k^{\prime}(\lambda_{1,1})k^{\prime}(\lambda_{1,2})|v(\lambda_{1,1})-v(\lambda_{1,2})|}.$
(IX.17)

After these preparational steps we evaluate the integrals over $\alpha_{1,1}$
and $\alpha_{1,2}$ with the help of the Dirac $\delta$-functions. This yields

$\displaystyle\Gamma_{nm}=$
$\displaystyle\frac{(2\pi)^{2}}{8}\int\differential\lambda_{1,1}\differential\lambda_{1,2}\differential\lambda_{2}\differential\alpha_{2}\,\frac{G(\mathbf{p}_{1},\mathbf{h}_{1},\mathbf{p}_{2},\mathbf{h}_{2})}{k^{\prime}(\lambda_{1,1})k^{\prime}(\lambda_{1,2})|v(\lambda_{1,1})-v(\lambda_{1,2})|}\times$
(IX.18) $\displaystyle\Big{(}(h_{n}^{\rm
Dr})^{\prime}(\lambda_{1,1})\bar{\alpha}_{1,1}+(h_{n}^{\rm
Dr})^{\prime}(\lambda_{1,2})\bar{\alpha}_{1,2}+(h_{n}^{\rm
Dr})^{\prime}(\lambda_{2})\alpha_{2}\Big{)}\times$
$\displaystyle\Big{(}(h_{m}^{\rm
dr})^{\prime}(\lambda_{1,1})\bar{\alpha}_{1,1}+(h_{m}^{\rm
dr})^{\prime}(\lambda_{1,2})\bar{\alpha}_{1,2}+(h_{m}^{\rm
dr})^{\prime}(\lambda_{2})\alpha_{2}\Big{)}.$

The second step is to evaluate the integral over $\alpha_{2}$ which is
proportional to the momentum transferred between the tubes. To this end we
observe that the combinations of functions appearing in the last two lines of
$\Gamma_{nm}$ can be rewritten according to

$\bar{\alpha}_{1,1}\varphi(\lambda_{1,1})+\bar{\alpha}_{1,2}\varphi(\lambda_{1,2})+\alpha_{2}\varphi(\lambda_{2})=\frac{\alpha_{2}k^{\prime}(\lambda_{2})}{v(\lambda_{1,1})-v(\lambda_{1,2})}\mathcal{S}[\varphi],$
(IX.19)

where we introduced symbol $\mathcal{S}[\varphi]$ defined by

$\mathcal{S}[\varphi]=\frac{\varphi^{\prime}(\lambda_{1,1})}{k^{\prime}(\lambda_{1,1})}(v(\lambda_{1,2})-v(\lambda_{2}))+\frac{\varphi^{\prime}(\lambda_{1,2})}{k^{\prime}(\lambda_{1,2})}(v(\lambda_{2})-v(\lambda_{1,1}))+\frac{\varphi^{\prime}(\lambda_{2})}{k^{\prime}(\lambda_{2})}(v(\lambda_{1,1})-v(\lambda_{1,2})).$
(IX.20)

We note that $\mathcal{S}[\varphi]$ is zero whenever any two of the $3$
variables $\lambda_{1,1}$, $\lambda_{1,2}$ or $\lambda_{2}$ are equal.
Moreover, it also vanishes identically for constant $\varphi$ or
$\varphi(\lambda)=k(\lambda),\omega(\lambda)$, and these functions naturally
correspond to the collision invariants. In the limit of small momentum
transfer between the tubes we can also simplify function $G$. The density
factors and the form-factors are independent of $\alpha_{2}$ in this limit.
The density factors are $J_{m{\rm ph}}(\mathbf{p},\mathbf{h})=J_{m{\rm
ph}}(\boldsymbol{\lambda},\boldsymbol{\lambda})$, while the form-factors yield
[64, 65, 66]

$\displaystyle F_{1\rm ph}(\mathbf{p}_{2},\mathbf{h}_{2})$
$\displaystyle=k^{\prime}(\lambda_{2}),$ (IX.21) $\displaystyle F_{2\rm
ph}(\mathbf{p}_{1},\mathbf{h}_{1})$ $\displaystyle=\frac{k^{2}\mathcal{T}^{\rm
dr}(\lambda_{1,1},\lambda_{1,2})}{k^{\prime}(\lambda_{1,1})k^{\prime}(\lambda_{1,2})\bar{\alpha}_{1,1}\bar{\alpha}_{1,2}}=-\mathcal{T}^{\rm
dr}(\lambda_{1,1},\lambda_{1,2})\frac{(v(\lambda_{1,1})-v(\lambda_{1,2}))^{2}}{(v(\lambda_{2})-v(\lambda_{1,2}))(v(\lambda_{2})-v(\lambda_{1,1}))},$
(IX.22)

with
$k=k^{\prime}(\lambda_{1,1})\bar{\alpha}_{1,1}+k^{\prime}(\lambda_{1,2})\bar{\alpha}_{1,2}$.
Here, $\mathcal{T}^{\rm dr}$ is the dressed integral kernel defined in (I.8).
With this we find

$\displaystyle\Gamma_{nm}=$
$\displaystyle\frac{(2\pi)^{2}}{8}\int\differential\lambda_{1,1}\differential\lambda_{1,2}\differential\lambda_{2}\differential\alpha_{2}\,\tilde{A}^{2}\big{(}k^{\prime}(\lambda_{2})\alpha_{2}\big{)}\,J_{2\rm
ph}(\boldsymbol{\lambda}_{1},\boldsymbol{\lambda}_{1})\,J_{1\rm
ph}(\lambda_{2},\lambda_{2})\times$ (IX.23)
$\displaystyle\frac{\alpha_{2}^{2}k^{\prime}(\lambda_{2})^{4}\mathcal{T}^{\rm
dr}(\lambda_{1,1},\lambda_{1,2})^{2}}{k^{\prime}(\lambda_{1,1})k^{\prime}(\lambda_{1,2})}\frac{|v(\lambda_{1,1})-v(\lambda_{1,2})|}{(v(\lambda_{2})-v(\lambda_{1,2}))^{2}(v(\lambda_{2})-v(\lambda_{1,1}))^{2}}\mathcal{S}[h_{n}^{\rm
Dr}]\mathcal{S}[h_{m}^{\rm dr}].$

The integral over $\alpha_{2}$ can be now performed. The result is the
following

$\displaystyle\Gamma_{nm}$
$\displaystyle=\frac{(2\pi)^{2}\tilde{A}^{2}}{8}\int\differential\lambda_{1,1}\differential\lambda_{1,2}\differential\lambda_{2}\frac{k^{\prime}(\lambda_{2})\mathcal{T}^{\rm
dr}(\lambda_{1,1},\lambda_{1,2})^{2}J_{2\rm
ph}(\boldsymbol{\lambda}_{1},\boldsymbol{\lambda}_{1})\,J_{1\rm
ph}(\lambda_{2},\lambda_{2}))}{k^{\prime}(\lambda_{1,1})k^{\prime}(\lambda_{1,2})}\times$
(IX.24)
$\displaystyle\frac{|v(\lambda_{1,1})-v(\lambda_{1,2})|}{(v(\lambda_{2})-v(\lambda_{1,2}))^{2}(v(\lambda_{2})-v(\lambda_{1,1}))^{2}}\mathcal{S}[h_{n}^{\rm
Dr}]\mathcal{S}[h_{m}^{\rm dr}],$

where $\tilde{A}^{2}=\int{\rm d}k\tilde{A}^{2}\big{(}k\big{)}k^{2}$. It is
convenient to rewrite $\Gamma_{nm}$ in the following form, which appears in
the main text

$\Gamma_{nm}=\frac{\tilde{A}^{2}}{8}\int{\rm d}\boldsymbol{\lambda}_{1}{\rm
d}\boldsymbol{\lambda}_{2}\,\mathcal{T}^{\rm
dr}(\lambda_{1},\lambda_{2})^{2}|v(\lambda_{1})-v(\lambda_{2})|g_{nm}(\lambda_{1},\lambda_{2}),$
(IX.25)

where the measure is $\differential\boldsymbol{\lambda}=2\pi\rho_{\rm
p}(\lambda)\rho_{\rm
h}(\lambda)/k^{\prime}(\lambda)\differential\lambda=\rho_{\rm
p}(\lambda)f(\lambda)\differential\lambda$ and

$g_{nm}(\lambda_{1},\lambda_{2})=\frac{1}{2\pi}\int{\rm
d}\boldsymbol{\lambda}(k^{\prime}(\lambda))^{2}\frac{\mathcal{S}[h_{n}^{\rm
Dr}]\mathcal{S}[h_{m}^{\rm
dr}]}{(v(\lambda)-v(\lambda_{1}))^{2}(v(\lambda)-v(\lambda_{2}))^{2}}.$
(IX.26)

The apparent poles in $\Gamma_{mn}$ at $\lambda_{1,1}=\lambda$ and
$\lambda_{1,2}=\lambda$ are regularized by the expressions in the brackets
which vanish in those limits. The double limit when all the $\lambda$’s are
equal is additionally regularized (and vanishes) because of the factor
$|v(\lambda_{1,1})-v(\lambda_{1,2})|$ in the numerator.

Explicitly, we check (we shorten the notation and introduce
$\tilde{\varphi}=\varphi^{\prime}/k^{\prime}$, moreover $v=v(\lambda)$ and
$v_{i}=v(\lambda_{i})$, similarly for $\tilde{\varphi}$)

$\displaystyle\lim_{v_{2}\rightarrow v}\frac{S[\varphi]}{v_{2}-v}$
$\displaystyle=\tilde{\varphi}_{1}+\lim_{v_{2}\rightarrow
v}\frac{\tilde{\varphi}_{2}(v-v_{1})+\tilde{\varphi}(v_{1}-v_{2})}{v_{2}-v}$
(IX.27) $\displaystyle=\tilde{\varphi}_{1}-v_{1}\lim_{v_{2}\rightarrow
v}\frac{\tilde{\varphi}_{2}-\tilde{\varphi}}{v_{2}-v}+\lim_{v_{2}\rightarrow
v}\frac{\tilde{\varphi}_{2}v-\tilde{\varphi}v_{2}}{v_{2}-v}$
$\displaystyle=\tilde{\varphi}_{1}-\tilde{\varphi}-(v_{1}-v)\tilde{\varphi}^{\prime}.$

The other apparent pole at $v=v_{1}$ is regularized in the same way. This a