# The dissipative Generalized Hydrodynamic equations and their numerical
solution

F. Møller1, N. Besse2, I.E. Mazets1,3, H.P. Stimming3, N.J. Mauser3

(  1 Vienna Center for Quantum Science and Technology (VCQ), Atominstitut, TU
Wien, 1020 Vienna, Austria 2 Observatoire de la Côte d’Azur, Univ. Côte
d’Azur, 06300 Nice, Cedex 4, France 3 Research platform MMM ”Mathematics-
Magnetism-Materials” c/o Fak. Mathematik, Univ. Wien, 1090 Vienna, Austria )

###### Abstract

“Generalized Hydrodynamics” (GHD) stands for a model that describes one-
dimensional integrable systems in quantum physics, such as ultra-cold atoms or
spin chains. Mathematically, GHD corresponds to nonlinear equations of kinetic
type, where the main unknown, a statistical distribution function
$f(t,z,\theta)$, lives in a phase space which is constituted by a one-
dimensional position variable $z$, and a one-dimensional ”kinetic” variable
$\theta$, actually a wave-vector, called “rapidity”. Two key features of GHD
equations are first a non-local and nonlinear coupling in the advection term,
and second an infinite set of conserved quantities, which prevent the system
from thermalizing. To go beyond this, we consider the dissipative GHD
equations, which are obtained by supplementing the right-hand side of the GHD
equations with a non-local and nonlinear diffusion operator or a Boltzmann-
type collision integral. In this paper, we deal with new high-order numerical
methods to efficiently solve these kinetic equations. In particular, we devise
novel backward semi-Lagrangian methods for solving the advective part (the so-
called Vlasov equation) by using a high-order time-Taylor series expansion for
the advection fields, whose successive time derivatives are obtained by a
recursive procedure. This high-order temporal approximation of the advection
fields are used to design new implicit/explicit Runge–Kutta semi-Lagrangian
methods, which are compared to Adams–Moulton semi-Lagrangian schemes. For
solving the source terms, constituted by the diffusion and collision
operators, we use and compare different numerical methods of the literature.

##  1 Introduction

Over the last decades, significant advances in experimentally realizing and
manipulating quantum many-body systems have been made. In particular, the
platform of ultracold atoms has demonstrated great success, owing to its
versatility; by controlling the trapping geometry or interactions of the
atoms, the effective Hamiltonian of the system can be tuned to follow various
condensed matter or field theory models [16, 28, 52]. Hence, ultracold gases
are sometimes referred to as ’analog quantum simulators’ [15, 55].

However, the behavior of such interacting quantum many-body systems is
notoriously complex, thus leading to a high demand for theoretical methods
capable of simulating their dynamics. Particularly successful theoretical
developments include, among others, the Density Matrix Renormalization Group
[34, 94, 95] and quantum Monte Carlo methods [1, 47, 86]. Alternative to
treating microscopic interactions, emergent properties of interacting quantum
many-body system can often be described in terms of quasi-particles
representing collective degrees of freedom. The concept of quasi-particles was
originally developed to study the transport of electrons in a solid [65], but
has since been employed to describe a number of phenomena, such as colossal
magnetoresistance [75] and superconductivity [67].

In integrable systems, quasi-particles descriptions are especially powerful,
as they provide an exact solution following the Bethe Ansatz [13, 69, 50, 63].
Such systems include, but are not limited to, repulsively interacting Bose
[71, 70] and Fermi [103, 49] gases and the anisotropic XXZ Heisenberg spin
chain [83]. Integrable systems share a number of key properties; they are one-
dimensional, all associated multi-particle scattering matrices factorize into
two-particle matrices, and they obey an infinite number of conservation laws.
Further, by virtue of their restricted geometry, the systems are intrinsically
strongly correlated and exhibit phenomena not found in higher-dimensional
system, such as the absence of thermalization [88, 27, 66]. The resulting
dynamics and relaxation properties are thus rather particular (see the special
issue [23]), as is perhaps best demonstrated in the seminal quantum Newton’s
cradle experiment [62].

Recently, the development of the theory of ”Generalized Hydrodynamics” (GHD)
[24, 7] has drastically simplified the study of quantum many-body dynamics in
integrable system by formulating the conservation laws as a kinetic equation
for the quasi-particles. The theory derives it name from the generalized Gibbs
ensemble [90, 89]; these statistical ensembles describe the non-thermal,
maximum-entropy steady states of integrable models. Following its inception, a
large number of extensions have been made to GHD (see the special issue [4]),
treating the spreading of entanglement and correlations [2], transport
properties [82, 19] and perturbations to integrability [78, 6]. The validity
of GHD and some of its extensions have been demonstrated following
observations from several experiments with (quasi) one-dimensional Bose gases
[93, 74, 78, 25] (see also the review [18]).

However, despite its significance and impact [80, 79, 42, 101, 40, 92, 77], no
systematic approach to develop and compare numerical methods has been done in
relation to GHD, with most studies employing only low-order schemes [20, 3,
21, 72, 5]. For advancing the use of GHD, especially for numerical modeling of
experiments, it is very useful to further develop and test numerical
techniques for solving the associated kinetic equations. In this work, we
develop a set of novel backward semi-Lagrangian (BSL) schemes of high-order
(implicit/explicit Runge–Kutta semi-Lagrangian schemes), which give very
accurate results proved by comparison with other high-order semi-Lagrangian
schemes (Adams–Moulton). Our results represent the highest numerical precision
solution of the GHD equation so far. Further, we present several versions of
BSL hybrid schemes for solving the GHD equation in the presence of diffusion
(higher order corrections to the typical Euler scale formulation of GHD) or
non-integrable perturbations (accounted for by a Boltzmann-type collision
integral). Of the latter we here treat the effect of highly energetic
collisions in a Bose gas, which can lead to 3D excitations of the atoms in the
trapping potential used to confine the gas. However, the methods presented can
be generalized to a number of other relevant perturbations [6, 43] in a
straightforward manner.

This paper is organized as follows: In Section 2, we introduce the theory of
Generalized Hydrodynamics (GHD) and present the key equations, namely the
advective and dissipative GHD equations. Of the latter, we focus on two
particular dissipative source terms: the diffusion operator and the Boltzmann-
type collision integral. Next, in Section 3, we derive exact time-derivatives
of the GHD advection fields, while in Section 4 we discuss how to solve the
dressing equation ubiquitous in GHD. In Section 5, we recall the backward
semi-Lagrangian method, followed, in Section 6, by the presentation of various
high-order time discretization schemes. In particular, we develop novel
implicit/explicit Runge–Kutta methods. In Section 7, we conduct a number of
numerical benchmarks, testing the accuracy of BSL schemes for solving the
advective GHD equation and BSL-hybrid schemes for solving the dissipative GHD
equation. Finally, in Section 8, we present some concluding remarks.

##  2 The GHD equations

Here, we present the GHD equations in two steps. In Section 2.1, we first
describe the GHD advection equation, which can be seen as a collisionless
kinetic equation of Vlasov type with non-local and nonlinear advection fields.
Then, in Section 2.2, we present two source terms for the right-hand side of
the GHD equation, resulting in a advective-dissipative equation. The first
source term is a diffusion operator (see Section 2.2.1), while the second one
is a Boltzmann-type collision operator (see Section 2.2.2). Both operators are
are non-local and nonlinear, but their dissipation mechanisms and mathematical
structures are very different.

###  2.1 The GHD advection equation

By virtue of the infinite number of conservation laws found in integrable
systems, their emergent hydrodynamics consist, in principle, of an infinite
set of continuity equations. Generalized Hydrodynamics employs the so-called
”thermodynamic Bethe Ansatz” [99] to formulate all the conservation laws as a
single kinetic equation. Following the thermodynamic Bethe Ansatz, the local
equilibrium state of an integrable systems is fully characterized by a density
of quasi-particles $\rho_{\mathrm{p}}$. The quasi-particles represent
collective excitations of the system and are uniquely labelled by their
wavevector $\theta$, the rapidity, which in the thermodynamic limit becomes a
continuous parameter. Note that with $\theta$ having dimension 1/lenght,
$p=\hbar\theta$ is a momentum. Assuming small variations between neighbouring
mesoscopic regions of the system, a single position-dependent quasi-particle
density describes the macrostate of the entire system. Large scale dynamics
can then be viewed as the propagation of quasi-particles between different
regions of the system. The mean propagation velocity of a given quasi-
particle, dubbed the effective velocity, is highly dependent on interactions
with other quasi-particles at the same point in space. At the lowest order
(Euler scale), where all currents depend only on the local densities of
conserved quantities, the propagation of the quasi-particles is purely
ballistic. Expressions for the currents were first presented in [24, 7] and
later more rigorously proven in [33, 17, 45, 22]. Including first-order
positional derivatives of the densities in the current functions results in
diffusive corrections to the quasi-particle propagation [36, 37].

Equivalently to the quasi-particle density, the local equilibrium state can be
fully characterized by the so-called filling function $f$, which is non-
negative and can be seen as a statistical distribution function of particle
occupation rate in the phase space. Indeed, the probability density function
$f$ encodes the fraction of possible momentum states occupied by the
particles. Furthermore, each rapidity component of the filling function can be
identified as a fluid mode [24, 39]. Following GHD, the infinite set of
advection equations for the fluid modes can be written as a single, two-
dimensional advection equation for the filling function $f=f(t,z,\theta)$,
namely

$\partial_{t}f+v^{\mathrm{eff}}\>\partial_{z}f+a^{\mathrm{eff}}\>\partial_{\theta}f=0,$
(1)

where $t\in{\mathbb{R}}_{+}$, $z\in{\mathbb{R}}$, and $\theta\in{\mathbb{R}}$,
represent, respectively, the time, position and momentum (rapidity) variables
[24, 7]. The advection equation (1), which is the main equation of the GHD
equations, must be supplemented with an initial condition,
$f(0,z,\theta)=f_{0}(z,\theta)$. The effective velocity
$v^{\mathrm{eff}}=v^{\mathrm{eff}}(t,z,\theta)$ is a non-local and nonlinear
expression in terms of $f$. Similarly, the effective acceleration
$a^{\mathrm{eff}}=a^{\mathrm{eff}}(t,z,\theta)$, accounting for inhomogeneous
Hamiltonians [3], is also in general a non-local and nonlinear expression in
terms of $f$. The effective velocity $v^{\mathrm{eff}}$ is defined by

$v^{\mathrm{eff}}:=\frac{(\partial_{\theta}\epsilon)^{\mathrm{dr}}}{(\partial_{\theta}p)^{\mathrm{dr}}},$
(2)

and the effective acceleration $a^{\mathrm{eff}}$ is defined by

$a^{\mathrm{eff}}:=\frac{(-\partial_{z}\epsilon)^{\mathrm{dr}}}{(\partial_{\theta}p)^{\mathrm{dr}}},$
(3)

where $\epsilon=\epsilon(z,\theta)$ and $p=p(\theta)$ are the single-particle
energy and momentum, respectively. Note that $a^{eff}$ has the dimension of an
acceleration if you divide by $\hbar m$ that is 1 in our scaling. Similar for
$v^{eff}$. The exact expressions for $\epsilon$ and $p$ depend on the
integrable model in question. Meanwhile, the superscript dr denotes the
dressing operation, which describes the modification of single-particle
properties following interactions with other particles. Given a state $f$, the
dressing operation is defined as

$g^{\mathrm{dr}}(\theta)=g(\theta)+\int_{{\mathbb{R}}}\mathrm{d}\theta^{\prime}\>T(\theta-\theta^{\prime})f(\theta^{\prime})g^{\mathrm{dr}}(\theta^{\prime}).$
(4)

Here, $T(\theta)$ is the scattering kernel given by

$T(\theta)=\frac{\partial_{\theta}\varphi(\theta)}{2\pi},$ (5)

with $\varphi(\theta)$ being the two-body scattering phase, which also depends
on the specific model. Since the effective velocity (and in general the
acceleration) are functionals of the distribution function $f$ itself, it
turns out that the GHD advection equation (1) is a nonlinear kinetic equation
with a self-consistent coupling.

Note, the dressing operation defines the function $g^{\mathrm{dr}}$ as the
solution of an integral equation in the rapidity variable $\theta$, with the
given source term $g$, and where the other variables, i.e. $(t,z)$, are fixed.
Introducing the following shorthand notation for the action of an integral
operator $\widehat{K}$ with kernel $K$ on any function $g$,

$\big{(}\widehat{K}g\big{)}(\theta)\coloneqq\int_{{\mathbb{R}}}\mathrm{d}\theta^{\prime}\>K(\theta,\theta^{\prime})g(\theta^{\prime}),$
(6)

the dressing operation rewrites as

$g^{\mathrm{dr}}=g+\widehat{T}fg^{\mathrm{dr}},$ (7)

and formally the dressing operator is $(1-\widehat{T}f)^{-1}$, i.e.
$g^{\mathrm{dr}}=(1-\widehat{T}f)^{-1}g$, for any function $g$.

As for collisionless kinetic equations of Vlasov type, a conservative form of
the GHD equations, written in terms of the quasi-particle density
$\rho_{\mathrm{p}}=\rho_{\mathrm{p}}(t,z,\theta)$, exists and reads as

$\partial_{t}\rho_{\mathrm{p}}+\partial_{z}\big{(}v^{\mathrm{eff}}\>\rho_{\mathrm{p}}\big{)}+\partial_{\theta}\big{(}a^{\mathrm{eff}}\>\rho_{\mathrm{p}}\big{)}=0.$
(8)

Likewise, the effective velocity and acceleration can be expressed in terms of
the quasi-particle density $\rho_{\mathrm{p}}$ yielding the equations

$v^{\mathrm{eff}}(t,z,\theta)=\partial_{\theta}\epsilon(z,\theta)+2\pi\int_{{\mathbb{R}}}\mathrm{d}\theta^{\prime}\>T(\theta-\theta^{\prime})\rho_{\mathrm{p}}(t,z,\theta^{\prime})\big{(}v^{\mathrm{eff}}(t,z,\theta^{\prime})-v^{\mathrm{eff}}(t,z,\theta)\big{)},$
(9)

and

$a^{\mathrm{eff}}(t,z,\theta)=-\partial_{z}\epsilon(z,\theta)+2\pi\int_{{\mathbb{R}}}\mathrm{d}\theta^{\prime}\>T(\theta-\theta^{\prime})\rho_{\mathrm{p}}(t,z,\theta^{\prime})\big{(}a^{\mathrm{eff}}(t,z,\theta^{\prime})-a^{\mathrm{eff}}(t,z,\theta)\big{)}.$
(10)

In fact, the conservative form of the GHD equations (8)-(10) is equivalent to
the advective form of the GHD equations, which are constituted by equations
(1), (2) and (3). Indeed, using the density of states
$\rho_{\mathrm{s}}=\rho_{\mathrm{s}}(t,z,\theta)$, defined by

$\rho_{\mathrm{s}}:=\frac{1^{\mathrm{dr}}}{2\pi},$ (11)

and the fundamental definition,

$f:=\frac{\rho_{\mathrm{p}}}{\rho_{\mathrm{s}}},$ (12)

we obtain, from the definitions of the dressing operation (4) and the
scattering kernel $T$ given by (5), that the equations (2)-(3) for the
advection field $\mathbf{F}\equiv(v^{\mathrm{eff}},\,a^{\mathrm{eff}})^{T}$ is
the same as (9)-(10). To pass from the advective form (1) to the conservative
form (8), we use the relation (12), and the fact that the density of states
$\rho_{\mathrm{s}}$, defined by (11), satisfies the same conservation law as
$\rho_{\mathrm{p}}$, namely,

$\partial_{t}\rho_{\mathrm{s}}+\partial_{z}\big{(}v^{\mathrm{eff}}\,\rho_{\mathrm{s}}\big{)}+\partial_{\theta}\big{(}a^{\mathrm{eff}}\,\rho_{\mathrm{s}}\big{)}=0.$
(13)

In Appendix A we show that the conservation law (13) for $\rho_{\mathrm{s}}$
can be recovered from the conservation law (8) for $\rho_{\mathrm{p}}$. In
addition, in Appendix A we show that the GHD equations preserve in time some
integral invariants such as the total number of particles $N$, the total
density of states $J$, the total energy $E$ and an infinite number of entropy
$S$. More precisely, for all $t\in{\mathbb{R}}_{+}$, we have,

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}N(t)=0,\quad\mbox{ with }\ \ \
N(t)=\int_{\mathbb{R}}\mathrm{d}z\int_{\mathbb{R}}\mathrm{d}\theta\,\rho_{\mathrm{p}},$
(14) $\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}J(t)=0,\quad\ \mbox{ with }\
\ \
J(t)=\int_{\mathbb{R}}\mathrm{d}z\int_{\mathbb{R}}\mathrm{d}\theta\,\rho_{\mathrm{s}},$
(15) $\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}E(t)=0,\quad\ \mbox{ with }\
\ \
E(t)=\int_{\mathbb{R}}\mathrm{d}z\int_{\mathbb{R}}\mathrm{d}\theta\,\rho_{\mathrm{p}}\epsilon,$
(16) $\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}S(t)=0,\quad\ \mbox{ with }\
\ \
S(t)=\int_{\mathbb{R}}\mathrm{d}z\int_{\mathbb{R}}\mathrm{d}\theta\,\rho_{\mathrm{p}}\Psi(f),$
(17)

where the function $\Psi:{\mathbb{R}}_{+}\mapsto{\mathbb{R}}$ is any smooth
function.

The Lieb–Liniger model [71, 70], which describes a one-dimensional Bose gas
with repulsive point-like interaction, is among the most common models treated
with GHD. Therefore, within the scope of this paper, we will concern ourselves
only with this particular model. However, by virtue of the universal
formulation of GHD, all results and methods presented can be generalized in a
straightforward manner to other integrable models [39]. Given a gas of $N$
bosons, the Lieb–Liniger Hamiltonian reads, after a rescaling s.t.
$\frac{\hbar^{2}}{2m}=1$

$\mathcal{H}=-\sum_{i}^{N}\frac{\partial^{2}}{\partial
z_{i}^{2}}+2c\sum_{i<j}^{N}\delta\left(z_{i}-z_{j}\right),$ (18)

where $c>0$ is the coupling constant of the atoms and $z_{i}$ the position of
the $i$’th boson. For the Lieb–Liniger model, we have
$\partial_{\theta}\varphi(\theta)=2c/(c^{2}+\theta^{2})$. Further, the single-
particle momentum and energy read

$p=p(\theta):=\theta,$ (19)

and

$\epsilon=\epsilon(z,\theta):=\theta^{2}+W(z,\theta),$ (20)

respectively, where $W=W(z,\theta)$ is a general potential. Most commonly the
potential $W$ is an external trapping potential $W(z,\theta)=V(z)$ used to
confine the atoms to a particular region in space. In this case, the effective
acceleration simplifies to
$a^{\mathrm{eff}}(t,z,\theta)=a^{\mathrm{eff}}(z)=-\partial_{z}V(z)$ [41].
Note that for the Lieb-Liniger model we have
$(\partial_{\theta}p)^{\mathrm{dr}}=1^{\mathrm{dr}}$.

###  2.2 The dissipative GHD equation

In the presence of diffusion or perturbations to integrability, the GHD
evolution of the filling function $f$ takes the general form of a advective-
dissipative equation

$\partial_{t}f+v^{\mathrm{eff}}\>\partial_{z}f+a^{\mathrm{eff}}\>\partial_{\theta}f=\mathcal{N}[f],$
(21)

where the source term $\mathcal{N}[f]$ is a nonlinear functional of the
filling function $f$. Unlike the GHD advection equation (1), the dissipative
GHD equation (21) no longer preserves the infinite set of conserved
quantities, and thus enables thermalization. We will be studying two different
nonlinear source terms $\mathcal{N}[f]$ in the following. First, the diffusion
operator $\mathcal{N}=\mathcal{D}$, which accounts for high-order corrections
to the GHD advection equation [36, 37, 5], and second the Boltzmann-type
collision integral $\mathcal{N}=\mathcal{I}$, which accounts for non-
integrable perturbations to the system [6, 43]. Here, we will be treating an
integrability-breaking perturbation in the form of 3D excitations of atoms in
the Bose gas [68, 78]. Generally, in the presence of either of these nonlinear
source terms, the total number of particles $N$ and the total energy $E$
remain conserved, while entropies $S$ decrease or increase monotonically,
according to the choice of the sign convention (see [37, 5, 78]). However, due
to simplifications made in the construction of the collision integral for the
3D excitations, the total energy is not fully conserved for all interaction
strengths of the Lieb–Liniger model. Lastly, an important difference between
these two dissipation operators is that the diffusion operator acts in the
variable $z$ locally (with a second-order derivative in this variable) and
non-locally in the variable $\theta$ (with some integral operators in this
variable), while the Boltzmann-type collision operator only acts in a non-
local way in the variable $\theta$ (without any derivative in the variables
$z$ and $\theta$).

####  2.2.1 Diffusion operator

The GHD model can be extended by adding a diffusive term : While considering
excitations with only one particle-hole pair above the reference state is
sufficient to derive the GHD on the Euler scale, the diffusion arises due to
simultaneous excitation of multiple (at least 2) particle-hole pairs. The
resonance conditions (i.e., the momentum and energy conservation) for the
scattering of two pairs are fulfilled in a general case, which leads to the
violation of the genuinely ballistic propagation of quasiparticles. More
detailed analysis [37] shows that the subleading to ballistic expansion order
in this propagation is diffusive (its contribution to the two-point
correlation functions of conserved charge grows linearly in time).

For the advective form of the extended GHD equation (21), the diffusion
operator reads [37, 5],

$\mathcal{D}_{f}[f]=\frac{1}{2\rho_{\mathrm{s}}}(1-f\widehat{T})\>\partial_{z}\big{(}(1-f\widehat{T})^{-1}\rho_{\mathrm{s}}\>\widehat{D}\>\partial_{z}f\big{)},$
(22)

which, notably, depends on the $z$-derivative of the filling function. Here,
the kernel of the integral operator $\widehat{D}$ is given by

$D(\theta,\alpha)=\frac{1}{\rho_{\mathrm{s}}^{2}(\theta)}\Big{(}\delta(\theta-\alpha)\int_{\mathbb{R}}\mathrm{d}\gamma\>W(\gamma,\alpha)-W(\theta,\alpha)\Big{)},$
(23)

where

$W(\theta,\alpha)=\rho_{\mathrm{p}}(\theta)\big{(}1-f(\theta)\big{)}\big{[}T^{\mathrm{dr}}(\theta,\alpha)\big{]}^{2}\big{|}v^{\mathrm{eff}}(\theta)-v^{\mathrm{eff}}(\alpha)\big{|}.$
(24)

In (24) the dressed scattering kernel $T^{\mathrm{dr}}(\theta,\alpha)$ is the
kernel of the operator $(1-\widehat{T}f)^{-1}\widehat{T}$, i.e.
$T^{\mathrm{dr}}(\theta,\alpha)=((1-\widehat{T}f)^{-1}T(\cdot-\alpha))(\theta)$,
or in other words the dressing of the function $T(\theta-\alpha)$ as a
function of its first argument $\theta$. Even if the variables $(t,z)$ are not
explicitly written in (23)-(24), the kernels $D(\theta,\alpha)$ and
$W(\theta,\alpha)$ depend on the time variable $t$ and the position variable
$z$, through the functions $\rho_{\mathrm{s}}$, $\rho_{\mathrm{p}}$, and $f$.
In the presence of diffusion and an external potential $V(z)$, the total
number of quasi-particles $N$, defined by (14), and the total energy $E$,
defined by (16), remain conserved throughout evolution [5] (see Appendix B for
a proof), while entropies $S$, defined by (17), decrease monotonically [37].
For this, it is more convenient to work with the conservative form of the GHD
equations, namely

$\partial_{t}\rho_{\mathrm{p}}+\partial_{z}\big{(}v^{\mathrm{eff}}\>\rho_{\mathrm{p}}\big{)}+\partial_{\theta}\big{(}a^{\mathrm{eff}}\>\rho_{\mathrm{p}}\big{)}=\mathcal{D}_{\rho_{\mathrm{p}}}[\rho_{\mathrm{p}}],$
(25)

where the corresponding diffusion operator
$\mathcal{D}_{\rho_{\mathrm{p}}}[\rho_{\mathrm{p}}]$ reads

$\mathcal{D}_{\rho_{\mathrm{p}}}[\rho_{\mathrm{p}}]=\frac{1}{2}\partial_{z}\big{(}(1-f\widehat{T})^{-1}\rho_{\mathrm{s}}\>\widehat{D}\>\rho_{\mathrm{s}}^{-1}(1-f\widehat{T})\>\partial_{z}\rho_{\mathrm{p}}\big{)}.$
(26)

Passing from the conservative form (25)-(26) to the advective form (21), with
$\mathcal{N}[f]=\mathcal{D}_{f}[f]$ and $\mathcal{D}_{f}[f]$ given by (22), is
detailed in [37].

####  2.2.2 Boltzmann-type collision integral

In an experimental setting, the Lieb-Liniger model can be realized by tightly
confining an ultracold Bose gas along two of its axes. Assuming the transverse
confinement to be harmonic with trapping frequency $\omega_{\perp}$ and width
$l_{\perp}=\sqrt{\hbar/m\omega_{\perp}}$, the dynamics of the gas is
effectively restricted to a single position axis if the transverse level
spacing $\hbar\omega_{\perp}$ exceeds all internal energy scales of the gas
[53, 54, 61, 64]. However, real systems are often only quasi one-dimensional,
as excited states of the transverse confinement can be populated through high
energy collisions of atoms. Such collisions break the integrability of the
system and eventually lead to thermalization [68]. To account for non-
integrable scattering processes within the framework of GHD, one can construct
a Boltzmann-type collision integral [43, 6]. The construction of the collision
integral accounting for transverse state-changing collisions is summarized in
Appendix C; the resulting collision integral is given by

$\mathcal{I}[g]=\sum_{n=1}^{2}\frac{1}{2}\big{(}\big{\\{}\mathcal{I}_{\mathrm{h}}^{+}[g]-\mathcal{I}_{\mathrm{p}}^{-}[g]\big{\\}}+\big{\\{}\mathcal{I}_{\mathrm{h}}^{-}[g]-\mathcal{I}_{\mathrm{p}}^{+}[g]\big{\\}}\nu_{n}^{\beta_{n}}\big{)},$
(27)

where $g$ is respectively the filling function $f$ or the density of quasi-
particles $\rho_{\mathrm{p}}$, depending on whether the advective or
conservative form of the extended GHD equations is considered [78, 25]. Here,
$\nu_{n}$ is the probability for an atom to be in the $n$’th transverse
excited state, while $\beta_{1}=2$ and $\beta_{2}=1$ are the number of atoms
changing state via the collisions. Note that the probability
$\nu_{n}=\nu_{n}(t)$ is a function of the time variable $t$. In addition, for
solving the dissipative GHD equation with the Boltzmann-type collision
integral (27), one must also solve the rate equations for the excitation
probabilities $\nu_{n}(t)$, namely

$\frac{\mathrm{d}\nu_{n}}{\mathrm{d}t}=\frac{1}{2}\beta_{n}\big{\\{}\Gamma_{\mathrm{h}}^{+}-\Gamma_{\mathrm{p}}^{+}\,\nu_{n}^{\beta_{n}}\big{\\}},$
(28)

where, for $\alpha\in\\{\mathrm{p},\,\mathrm{h}\\}$,

$\Gamma_{\alpha}^{+}=(2N)^{-1}\int_{\mathbb{R}}\mathrm{d}z\int_{\mathbb{R}}\mathrm{d}\theta\,\rho_{\mathrm{s}}\,\mathcal{I}_{\alpha}^{+}[f]=(2N)^{-1}\int_{\mathbb{R}}\mathrm{d}z\int_{\mathbb{R}}\mathrm{d}\theta\,\mathcal{I}_{\alpha}^{+}[\rho_{\mathrm{p}}].$
(29)

The terms $\mathcal{I}_{\alpha}^{\pm}$ of eq. (27) are collision integrals of
individual scattering processes, which for the advective form of the
dissipative GHD equation (21) can be written in the form

$\mathcal{I}_{\alpha}^{\pm}[f](\theta)=\frac{(2\pi)^{2}\hbar}{m}\int_{\mathcal{R}_{\pm}}\mathrm{d}\theta^{\prime}\;\lvert\theta-\theta^{\prime}\rvert\>P_{\updownarrow}(\lvert\theta-\theta^{\prime}\rvert,\,\lvert\theta_{\pm}-\theta^{\prime}_{\pm}\rvert)\>f_{\alpha}(\theta)f_{\alpha}(\theta^{\prime})f_{\bar{\alpha}}(\theta_{\pm})f_{\bar{\alpha}}(\theta^{\prime}_{\pm})\>\rho_{\mathrm{s}}(\theta^{\prime})\rho_{\mathrm{s}}(\theta_{\pm})\rho_{\mathrm{s}}(\theta_{\pm}^{\prime}).$
(30)

Here, $\bar{\alpha}={\mathrm{h}}$, for $\alpha={\mathrm{p}}$ and vice versa,
while $f_{\mathrm{p}}=f$ and $f_{\mathrm{h}}=1-f$. Further,
$P_{\updownarrow}(\theta_{1},\,\theta_{2})=4c^{2}\theta_{1}\theta_{2}/[\theta_{1}^{2}\theta_{2}^{2}+c^{2}(\theta_{1}+\theta_{2})^{2}]$
is the scattering probability, while
$\theta_{\pm}=\frac{1}{2}(\theta+\theta^{\prime})+\frac{1}{2}(\theta-\theta^{\prime})\sqrt{1\pm
8/[(\theta-\theta^{\prime})l_{\perp}]^{2}}$ and
$\theta_{\pm}^{\prime}=\frac{1}{2}(\theta+\theta^{\prime})-\frac{1}{2}(\theta-\theta^{\prime})\sqrt{1\pm
8/[(\theta-\theta^{\prime})l_{\perp}]^{2}}$ are the rapidities after a
collision leading to excitation or de-excitation of the transverse states,
respectively. The integration ranges in equation (30) are the following:
$\mathcal{R}_{+}$ is the whole real axis, and $\mathcal{R}_{-}$ is comprised
of those real values of $\theta^{\prime}$, which yield real $\theta_{-}$ and
$\theta_{-}^{\prime}$, i.e.
$\mathcal{R}_{-}=\\{\theta^{\prime}:\theta^{\prime}<\theta-\sqrt{8}/l_{\perp}\\}\cup\\{\theta^{\prime}:\theta^{\prime}>\theta+\sqrt{8}/l_{\perp}\\}$.
Even if the variables $(t,z)$ are not explicitly written in eqs. (27)-(30),
the collision integrals $\mathcal{I}_{\alpha}^{\pm}[f](\theta)$ depend on the
time variable $t$ and position variable $z$ through the functions
$\rho_{\mathrm{s}}$ and $f$.

Boltzmann-type collision integrals are explicitly constructed to preserve the
total number of particles $N$. Generally, this is also the case for the total
energy $E$, however, in the derivation of the collision integral (30) the
backflow, describing a shift in all local rapidities following an excitation
[70], was neglected. Hence, the associated energy shift is not accounted for
and the total energy is not conserved, although for most experimentally
realizable parameters of the Lieb–Liniger model the error is not too severe
[78]. For more details, see Appendix C.

Again, the conservation of particle number is most obvious when examining its
equivalent conservative form given by

$\partial_{t}\rho_{\mathrm{p}}+\partial_{z}\big{(}v^{\mathrm{eff}}\>\rho_{\mathrm{p}}\big{)}+\partial_{\theta}\big{(}a^{\mathrm{eff}}\>\rho_{\mathrm{p}}\big{)}=\mathcal{I}[\rho_{\mathrm{p}}].$
(31)

The form of the collision integral $\mathcal{I}[\rho_{\mathrm{p}}]$, which is
originally constructed in [78], has the same form as equation (27), however,
the collision integrals of the individual collision channels now read

$\mathcal{I}_{\alpha}^{\pm}[\rho_{\mathrm{p}}](\theta)=\frac{(2\pi)^{2}\hbar}{m}\int_{\mathcal{R}_{\pm}}\mathrm{d}\theta^{\prime}\,|\theta-\theta^{\prime}|P_{\updownarrow}(|\theta-\theta^{\prime}|,\,|\theta_{\pm}-\theta^{\prime}_{\pm}|)\rho_{\alpha}(\theta)\rho_{\alpha}(\theta^{\prime})\rho_{\bar{\alpha}}(\theta_{\pm})\rho_{\bar{\alpha}}(\theta^{\prime}_{\pm}),$
(32)

where $\bar{\alpha}=\mathrm{h}$, for $\alpha=\mathrm{p}$ and vice versa, and
with $\rho_{\bar{\alpha}}=\rho_{\mathrm{s}}-\rho_{\alpha}$.

##  3 Analytic time derivatives of the advection fields

In this section, we derive analytic time derivatives of the advection fields
$v^{\mathrm{eff}}$ and $a^{\mathrm{eff}}$, which will be used to high-order
schemes in time for solving the GHD advection equation (1). For this, we
define auxillary advection fields $w^{\mathrm{eff}}$ and $b^{\mathrm{eff}}$ as
follows,

$w^{\mathrm{eff}}\coloneqq 1^{\mathrm{dr}}v^{\mathrm{eff}},\quad\mbox{ and
}\quad b^{\mathrm{eff}}\coloneqq 1^{\mathrm{dr}}a^{\mathrm{eff}}.$

Given the definition of the dressing equation (4), $w^{\mathrm{eff}}$ and
$b^{\mathrm{eff}}$ solve the following equations,

$\displaystyle\begin{split}w^{\mathrm{eff}}(\theta)=\partial_{\theta}\epsilon+\int_{\mathbb{R}}\mathrm{d}\theta^{\prime}\>T(\theta-\theta^{\prime})f(\theta^{\prime})w^{\mathrm{eff}}(\theta^{\prime}),\\\
b^{\mathrm{eff}}(\theta)=-\partial_{z}\epsilon+\int_{\mathbb{R}}\mathrm{d}\theta^{\prime}\>T(\theta-\theta^{\prime})f(\theta^{\prime})b^{\mathrm{eff}}(\theta^{\prime}),\end{split}$

where $\epsilon=\epsilon(\theta,z)=\theta^{2}+V(z)$ is the one-particle
Hamiltonian. Note that the filling function $f$ is also a function of time and
position, i.e. $f=f(t,z,\theta)$; although these arguments will be omitted in
the following whenever they appear the same on either side of an equality.
Differentiating $w^{\mathrm{eff}}$ in time, we obtain

$\partial_{t}w^{\mathrm{eff}}(\theta)=\int_{\mathbb{R}}\mathrm{d}\theta^{\prime}\>T(\theta-\theta^{\prime})\partial_{t}f(\theta^{\prime})w^{\mathrm{eff}}(\theta^{\prime})+\int_{\mathbb{R}}\mathrm{d}\theta^{\prime}\>T(\theta-\theta^{\prime})f(\theta^{\prime})\partial_{t}w^{\mathrm{eff}}(\theta^{\prime}).$

Here, $\partial_{t}f(\theta^{\prime})$ is given through the GHD advection
equation (1), implying that the time derivative of $f$ is expressed in terms
of the position and rapidity derivatives of $f$, and of other known
quantities, such as the advection fields
$(v^{\mathrm{eff}},\,a^{\mathrm{eff}})$. Defining the source term $S_{w}$ as

$S_{w}(\theta)\coloneqq\int_{\mathbb{R}}\mathrm{d}\theta^{\prime}\>T(\theta-\theta^{\prime})\partial_{t}f(\theta^{\prime})w^{\mathrm{eff}}(\theta^{\prime}),$

we see that

$\partial_{t}w^{\mathrm{eff}}(\theta)=S_{w}^{\mathrm{dr}}(\theta).$

Likewise, a source term for $b^{\mathrm{eff}}$ can be defined such that

$\partial_{t}b^{\mathrm{eff}}(\theta)=S_{b}^{\mathrm{dr}}(\theta),$

where

$S_{b}(\theta)=\int_{\mathbb{R}}\mathrm{d}\theta^{\prime}\>T(\theta-\theta^{\prime})\partial_{t}f(\theta^{\prime})b^{\mathrm{eff}}(\theta^{\prime}).$

Knowing the derivatives $\partial_{t}w^{\mathrm{eff}}$ and
$\partial_{t}b^{\mathrm{eff}}$, the time derivatives of the advection fields
can be computed following

$\displaystyle\begin{split}\partial_{t}v^{\mathrm{eff}}=\frac{\partial_{t}w^{\mathrm{eff}}-\partial_{t}1^{\mathrm{dr}}v^{\mathrm{eff}}}{1^{\mathrm{dr}}},\\\
\partial_{t}a^{\mathrm{eff}}=\frac{\partial_{t}b^{\mathrm{eff}}-\partial_{t}1^{\mathrm{dr}}a^{\mathrm{eff}}}{1^{\mathrm{dr}}},\end{split}$
(33)

where the unit-cell $1^{\mathrm{dr}}$ obeys the conservation law

$\partial_{t}1^{\mathrm{dr}}+\partial_{z}\big{(}v^{\mathrm{eff}}\>1^{\mathrm{dr}}\big{)}+\partial_{\theta}\big{(}a^{\mathrm{eff}}\>1^{\mathrm{dr}}\big{)}=0.$
(34)

Thereby, knowing both $1^{\mathrm{dr}}$ and the advection fields
$(v^{\mathrm{eff}},\,a^{\mathrm{eff}})$, we can compute
$\partial_{t}1^{\mathrm{dr}}$ from (34) and in turn the time derivatives of
the advection fields $(v^{\mathrm{eff}},\,a^{\mathrm{eff}})$ from (33).
Higher-order time derivatives of the advection fields can be derived in
similar fashion. If we let $S_{w}^{(1)}\equiv S_{w}$, we can construct higher-
order source terms fulfilling

$\partial_{t}^{(k)}w^{\mathrm{eff}}(\theta)=\big{(}S_{w}^{(k)}\big{)}^{\mathrm{dr}}(\theta),$

where $S_{w}^{(k)}(\theta)$ depends on all previous source terms
$S_{w}^{(j)}(\theta)$ with $j<k$. Of course, we have similar expressions for
$\partial_{t}^{(k)}b^{\mathrm{eff}}(\theta)$ in terms of
$S_{b}^{(k)}(\theta)$. For example, the second-order source term
$S_{w}^{(2)}(\theta)$ reads

$S_{w}^{(2)}(\theta)=\int_{\mathbb{R}}\mathrm{d}\theta^{\prime}\>T(\theta-\theta^{\prime})\partial_{t}^{(2)}f(\theta^{\prime})w^{\mathrm{eff}}(\theta^{\prime})+2\int_{\mathbb{R}}\mathrm{d}\theta^{\prime}\>T(\theta-\theta^{\prime})\partial_{t}f(\theta^{\prime})\partial_{t}w^{\mathrm{eff}}(\theta^{\prime}).$
(35)

In order to obtain high-order time derivatives of $f$, such as
$\partial_{t}^{(2)}f$ appearing in (35), we use the successive time
derivatives of the GHD advection equation (1), where for any mixed derivative
of $f$, containing some time derivatives, each time derivative of $f$ is
replaced by the $(z,\theta)$-derivative of $f$ through the GHD advection
equation (1). In such a way, high-order time derivatives of $f$ can be
expressed only in terms of the $(z,\theta)$-derivatives of $f$ and of other
known quantities, such as the previously calculated time derivatives of the
advection fields $(v^{\mathrm{eff}},\,a^{\mathrm{eff}})$. In the same spirit,
we use the successive time derivatives of equations (33)-(34) to obtain the
high-order time derivatives of $(v^{\mathrm{eff}},\,a^{\mathrm{eff}})$ from
those of $(w^{\mathrm{eff}},\,b^{\mathrm{eff}})$. Finally, knowing the
advection fields $(v^{\mathrm{eff}},\,a^{\mathrm{eff}})$ at time $\tau$ and
their time derivatives up to the order $k$, we can obtain a high-order
approximation of the advection fields at any time $t$ following the time-
Taylor series expansion,

$\displaystyle\begin{split}v^{\mathrm{eff}}(t)&=v^{\mathrm{eff}}(\tau)+(t-\tau)\>\partial_{t}v^{\mathrm{eff}}(\tau)+\ldots+\frac{(t-\tau)^{k}}{k!}\>\partial_{t}^{(k)}v^{\mathrm{eff}}(\tau),\\\
a^{\mathrm{eff}}(t)&=a^{\mathrm{eff}}(\tau)+(t-\tau)\>\partial_{t}a^{\mathrm{eff}}(\tau)+\ldots+\frac{(t-\tau)^{k}}{k!}\>\partial_{t}^{(k)}a^{\mathrm{eff}}(\tau).\end{split}$
(36)

##  4 Solving the dressing equation

The dressing operation, given by equation (4), is ubiquitous in GHD, as it
describes how a single particle, or bare function is modified in the presence
of collective interactions. Importantly, a dressed quantity becomes a function
of time and position (and of course rapidity) through the dependence on the
filling function $f(t,z,\theta)$. Therefore, at every step of time
propagation, the dressing operation presents a new set of equations to be
solved.

To solve the dressing equation, we first discretize position and rapidity
creating the grids $z_{i}\in[z_{1},z_{2},\ldots,z_{N_{z}}]$ and
$\theta_{j}\in[\theta_{1},\theta_{2},\ldots,\theta_{N_{\theta}}]$. Although
the rapidity integral in equation (4) is taken over the whole real axis, in
practice it is sufficient to integrate over the interval
$[\theta_{\mathrm{min}},\theta_{\mathrm{max}}]$ containing the support in
$\theta$ of the filling function $f$. Since the dressing operation is
completely local in position, the equation can be solved independently for
every point in position $z_{i}$. Then, in discretized form, the dressing
equation for a single point in position reads

$g^{\mathrm{dr}}(\theta_{i})=g(\theta_{i})+\sum_{j=1}^{N_{\theta}}\mathrm{w}_{j}\,T(\theta_{i}-\theta_{j})\,f(\theta_{j})g^{\mathrm{dr}}(\theta_{j}),$
(37)

where $\mathrm{w}_{j}$ are the quadrature weights for the corresponding
discretization scheme of the rapidity axis. Note that the position and time
arguments have been omitted entirely for a more compact notation, as they are
the same on either side of the equality sign. The discretized form of the
dressing equation (37) can be rewritten as

$g(\theta_{i})=\sum_{j=1}^{N_{\theta}}\big{[}\delta_{ij}-\mathrm{w}_{j}\,T(\theta_{i}-\theta_{j})\,f(\theta_{j})\big{]}g^{\mathrm{dr}}(\theta_{j}),$
(38)

where $\delta_{ij}$ is the standard Kronecker symbol. Hence, equation (38) can
be written in compact form as the matrix product $g=\mathbf{U}\cdot
g^{\mathrm{dr}}$, where the elements of the matrix $\mathbf{U}$ are
$\mathbf{U}_{ij}=\delta_{ij}-\mathrm{w}_{j}\,T(\theta_{i}-\theta_{j})\,f(\theta_{j})$.
Applying the dressing operation numerically thus amounts to solving the set of
linear equations of equation (38) either through inverting the matrix
$\mathbf{U}$ or, as in the case of the iFluid framework [80], using a
subroutine such as MATLAB mldivide function [76], capable of leveraging
symmetries of the matrix.

Figure 1:  Comparison of effective velocity (2) calculated for the test state
(39) using the trapezoidal (blue) and Simpson’s (red) rule to discretize the
integral over rapidities featuring in the dressing equation (4). The effective
velocity is calculated and plotted for increasingly fine discretization of the
rapidity axis ($N_{\theta}$ is the number of grid points). Below, the
difference between the two results are plotted.

In general, when solving the dressing equation it is most convenient to employ
the same discretization scheme used in the backward semi-Lagrangian method.
Otherwise interpolation of the filling function $f$ between the two grids must
be performed at every evaluation of the dressing equation. Hence, for the
entirety of this work, a uniform discretization of both the position and
rapidity axis will be employed. For such discretization several quadrature
rules exist to approximate this integral. Here, we test two of the most
standard: the trapezoidal rule and the Simpson’s rule. Following the
trapezoidal rule, the quadrature weights $\mathrm{w}_{j}$ read

$\\{\mathrm{w}_{1},\mathrm{w}_{2},\mathrm{w}_{3},\ldots,\mathrm{w}_{N_{\theta}-1},\mathrm{w}_{N_{\theta}}\\}=\frac{(\theta_{\mathrm{max}}-\theta_{\mathrm{min}})}{2(N_{\theta}-1)}\\{1,2,2,\ldots,2,1\\},$

while for the Simpson’s rule they read

$\\{\mathrm{w}_{1},\mathrm{w}_{2},\mathrm{w}_{3},\mathrm{w}_{4},\mathrm{w}_{5},\ldots,\mathrm{w}_{N_{\theta}-2},\mathrm{w}_{N_{\theta}-1},\mathrm{w}_{N_{\theta}}\\}=\frac{(\theta_{\mathrm{max}}-\theta_{\mathrm{min}})}{3(N_{\theta}-1)}\\{1,4,2,4,2\ldots,2,4,1\\}.$

To compare the two quadrature rules, we compute the effective velocity (2) at
different grid resolutions of the reference state

$f_{0}(\theta)=0.9\left[\exp\left(-\frac{\left(\theta-\theta_{\mathrm{B}}\right)^{2}}{2\sigma}\right)+\exp\left(-\frac{\left(\theta+\theta_{\mathrm{B}}\right)^{2}}{2\sigma}\right)\right],$
(39)

where $\theta_{\mathrm{B}}=2$ and $\sigma=1/\sqrt{2}$. The results can be seen
in Fig. 1. Although a clear systematic difference between the two quadrature
rules is apparent, the difference is nonetheless almost vanishing for
$N_{\theta}>100$. Then, in order to keeps matters simple, for the entirety of
this work the trapezoidal rule will be employed.

##  5 The backward semi-Lagrangian (BSL) method

The characteristic equations associated with the first-order PDE (1) are

$\displaystyle\frac{\mathrm{d}Z(t)}{\mathrm{d}t}=v^{\mathrm{eff}}(t,Z(t),\Theta(t))\quad,\quad
Z(t;t,z,\theta)=z,$ (40a)
$\displaystyle\frac{\mathrm{d}\Theta(t)}{\mathrm{d}t}=a^{\mathrm{eff}}(t,Z(t),\Theta(t))\quad,\quad\Theta(t;t,z,\theta)=\theta,$
(40b)

while the filling function $f$ is constant along the characteristics, i.e.,

$f(t,z,\theta)=f(s,Z(s;t,z,\theta),\Theta(s;t,z,\theta)),\quad\forall
s\in{\mathbb{R}}.$ (41)

We recall that the advection field
$\mathbf{F}:=(v^{\mathrm{eff}},\,a^{\mathrm{eff}})^{T}$ is a functional of the
filling function $f$, through the integral equations (2)-(3).

Several methods for numerically solving equation (1) exist. In purely
Lagrangian methods, the implicit solution (41) is utilized by considering the
initial state $f(0,z,\theta)=f_{0}(z,\theta)$ as a distribution of fluid
packages, whose trajectories, following time evolution, are encoded in the
characteristics defined by equations (40). However, following dynamics, the
fluid packages will be scattered around the entire phase space, which may
result in the mesh needed for evaluating the effective velocity becoming
heavily distorted and irregular [102]. Backward Semi-Lagrangian (BSL) schemes
[12] circumvent this issue by assuming the arrival point of the
characteristics at each time step to coincide with a pre-defined grid,
effectively re-meshing at every time step [98]. The BSL schemes then proceed
in the following way. First, we discretize position and rapidity, creating the
grids $z_{i}\in\\{z_{1},z_{2},\ldots,z_{N_{z}}\\}$ and
$\theta_{j}\in\\{\theta_{1},\theta_{2},\ldots,\theta_{N_{\theta}}\\}$. For a
single time step $\Delta t$, the GHD advection equation (1), integrated in
time, writes as

$f(t^{n+1},z_{i},\theta_{j})=f(t^{n},Z(t^{n};t^{n+1},z_{i},\theta_{j}),\Theta(t^{n};t^{n+1},z_{i},\theta_{j})),$
(42)

where $t^{n}=n\Delta t$. To compute $f(t^{n+1})$ at the grid points
$(z_{i},\theta_{j})$, we consider the characteristics which arrive at time
$t^{n+1}$ on the grid points $(z_{i},\theta_{j})$. Then, the position and
rapidity couple
$\big{(}Z(t^{n};t^{n+1},z_{i},\theta_{j}),\Theta(t^{n};t^{n+1},z_{i},\theta_{j})\big{)}$
denotes the departure point of the fluid package at time $t^{n}$, which at the
next time $t^{n+1}$ will arrive at the grid point $(z_{i},\theta_{j})$. For a
more compact notation we write the phase-space coordinate as
$\mathbf{x}\equiv(z,\theta)^{T}$, the characteristics as
$\mathbf{X}(t^{n})\equiv\big{(}Z(t^{n};t^{n+1}),\Theta(t^{n};t^{n+1})\big{)}^{T}$,
and the advection fields as
$\mathbf{F}\equiv(v^{\mathrm{eff}},\,a^{\mathrm{eff}})^{T}$. Hence, a solution
to equation (1) is obtained from (42) and by following the characteristics X
backward in time, starting from the grid points $\mathbf{x}_{ij}$, i.e. by
solving,

$\begin{dcases}\frac{\mathrm{d}\mathbf{X}(t)}{\mathrm{d}t}&=\mathbf{F}(t,\mathrm{X}(t))\\\
\mathbf{X}(t^{n+1})&\equiv\mathbf{X}^{n+1}\coloneqq\mathbf{x}_{ij}.\end{dcases}$
(43)

Note that the characteristics in general do not coincide with the grid points.
So we must employ interpolation in order to obtain values of the filling
function $f$ or the advection field F at any point in the phase space. With
this in mind, a single propagation step following a BSL scheme can be divided
into two sub-steps: First, the origin of the characteristics
$\mathbf{X}(t^{n})$ is computed. An illustration of this can be seen in Fig.
2. Next, the filling function at the origin $f(t^{n},\textbf{X}(t^{n}))$ is
reconstructed using interpolation from the known values on the mesh. Since the
value of the filling function is constant along the characteristics, the state
at time $t^{n+1}$ is obtained as follows,

$f^{n+1}(\mathbf{x}_{ij})=\Pi_{h}f^{n}(\mathbf{X}^{n}),$

where $\Pi_{h}$ denotes the interpolation operator. The choice of $\Pi_{h}$ is
crucial for the accuracy and stability of the scheme [8, 9, 11]. In the
following, we will employ an interpolation scheme based on a cubic spline
using not-a-knot end conditions [35, 96].

Figure 2:  Illustration of a single propagation step using a backward semi-
Lagrangian scheme. At time $t^{n}$ the values of the distribution $f$ are
known on the mesh (dotted lines). Over the time interval $[t^{n},t^{n+1}]$ the
distribution propagates following trajectories, along which it is constant.
The arrival point $\textbf{X}(t^{n+1})$ (marked as a circle) is assumed to
coincide with the mesh, while the departure point $\textbf{X}(t^{n})$ (marked
as a square) is unknown. The solution of the corresponding ODEs (43), plotted
as a solid line, can be approximated using various numerical schemes (dashed
line).

##  6 High-order schemes in time to solve the BSL method

In this section we present a number of high-order schemes for solving the
advective and dissipative GHD equations. First, in Section 6.1, we focus on
the purely advective case, for which we develop a set of novel
implicit/explicit Runge–Kutta semi-Lagrangian methods. These methods employ
the exact time-derivatives of the GHD advection fields derived in Section 3 to
estimate the advection fields at intermediate time steps. For comparison, we
introduce both Adams–Moulton schemes in Section 6.2 and leap-frog schemes in
Section 6.3. Finally, in Section 6.4 we present a number of BSL-hybrid schemes
for solving the dissipative GHD equation.

###  6.1 Implicit/explicit Runge–Kutta semi-Lagrangian schemes

Runge–Kutta methods are a family of methods that employ a number of
intermediate time steps in order to obtain a high-order approximation for the
solution of ordinary differential equations, such as the one of equation (43).
Suppose that one knows the advection field $\mathbf{F}(t,\cdot)$ at any time
$t$. For the purpose of solving the GHD advection equation (1), and following
the recursive method of Section 3, the advection field $\mathbf{F}(t,\cdot)$
can be predicted by computing successively the time derivatives
$\partial_{t}^{(k)}\mathbf{F}(t^{n},\cdot)$, for $k\geq 0$, and then using the
time-Taylor series expansion (36) to reconstruct $\mathbf{F}(t,\cdot)$ at any
time $t\in[t^{n},t^{n+1}]$ (see Section 3 for more details). Then, using the
Runge–Kutta schemes of order $s$, that we formally write as
$\mathrm{RK}_{\mathbf{F}}^{(s)}(\Delta
t;t^{n},t^{n+c_{2}},t^{n+c_{3}},\ldots,t^{n+c_{s}};\mathbf{X}(t^{n}))$, with
$0<c_{i}<1$ and $i\in\\{2,\ldots,s\\}$, we obtain

$\displaystyle\mathbf{X}(t^{n+1})-\mathbf{X}(t^{n})$
$\displaystyle=\int_{t^{n}}^{t^{n+1}}\mathbf{F}(s,\mathbf{X}(s))\,\mathrm{d}s$
(44) $\displaystyle\simeq\Delta t\>\mathrm{RK}_{\mathbf{F}}^{(s)}\big{(}\Delta
t;t^{n},t^{n+c_{2}},t^{n+c_{3}},\ldots,t^{n+c_{s}};\mathbf{X}(t^{n})\big{)}.$

The function $\mathrm{RK}_{\mathbf{F}}^{(s)}(\Delta
t;t^{n},t^{n+c_{2}},t^{n+c_{3}},\ldots,t^{n+c_{s}};\cdot)$ is known explicitly
in terms of the functions $\mathbf{F}(t^{n},\cdot)$,
$\mathbf{F}(t^{n+c_{2}},\cdot)$, $\ldots$, $\mathbf{F}(t^{n+c_{s}},\cdot)$.
Then, to obtain the departure points $\mathbf{X}(t^{n})$ from equation (44),
one must solve the following fixed-point problem corresponding to finding the
zeroes of the following nonlinear function,

$\mathbf{X}\to G(\mathbf{X})\coloneqq\mathbf{x}_{ij}-\mathbf{X}-\Delta
t\>\mathrm{RK}_{\mathbf{F}}^{(s)}\big{(}\Delta
t;t^{n},t^{n+c_{2}},t^{n+c_{3}},\ldots,t^{n+c_{s}};\mathbf{X}\big{)}.$ (45)

The problem, $G(\mathbf{X})=0$, can be solved by using the Picard iteration
(first order) or the Newton iteration (second order). Here, we employ the
former, which usually converges in less than five iterations. While the fixed-
point problem of equation (45) is implicit, one can also construct explicit
Runge–Kutta schemes such as

$\mathbf{X}(t^{n})\simeq\mathbf{x}_{ij}-\Delta
t\>\mathrm{RK}_{\mathbf{F}}^{(s)}\big{(}\Delta
t;t^{n},t^{n+c_{2}},t^{n+c_{3}},\ldots,t^{n+c_{s}};\mathbf{x}_{ij}\big{)}.$

Note that the functions $\mathrm{RK}_{\mathbf{F}}^{(s)}$ used in the implicit
and explicit schemes are usually not the same. Explicit schemes rely only on
known quantities and do not require to solve any fixed-point problem by
successive iterations, whereby they often are faster than implicit schemes.
However, given two schemes at the same order, the implicit scheme is often
more accurate, allowing for overall larger time steps.

####  6.1.1 First-order schemes (RK1)

Here, we consider a first order scheme, just to demonstrate how to compute the
origin of the characteristics. For most practical purposes first-order schemes
are not used, as they become too diffusive and are not accurate enough. At
first order, the implicit Runge–Kutta method reads

$\mathrm{RK}_{\mathbf{F}}^{(1)}(\Delta
t;\mathbf{X})=\mathbf{F}(t^{n},\mathbf{X}).$

To evaluate $\mathbf{F}(t^{n},\mathbf{X})$ at any points $(z,\theta)$ in the
phase-space, we need to interpolate them from the known grid values
$\mathbf{F}_{ij}(t^{n})$. Historically, this scheme was first employed in [20]
to solve equation (1) in the absence of any external potential. Similarly, we
can construct an explicit first-order Runge–Kutta scheme such as

$\mathrm{RK}_{\mathbf{F}}^{(1)}(\Delta
t;\mathbf{x}_{ij})=\mathbf{F}(t^{n+1},\mathbf{x}_{ij}),$

where $\mathbf{F}(t^{n+1},\mathbf{x}_{ij})\equiv\mathbf{F}_{ij}(t^{n+1})$ is
estimated via equation (36).

####  6.1.2 Second-order schemes (RK2)

At second order, the Runge–Kutta schemes employ the midpoint advection field
$\mathbf{F}(t^{n+1/2})$ to obtain a more accurate estimate of the origin of
the characteristics. To estimate the midpoint advection field on the grid we
employ the time-Taylor series expansion (36) and the recursive procedure of
Section 3 to compute the time derivatives of the advection field
$\\{\partial_{t}^{(k)}\mathbf{F}(t)\\}_{k\geq 0}$. Next, we use cubic-spline
interpolation to evaluate $\mathbf{F}(t^{n+1/2},\mathbf{X})$ at any points
$(z,\theta)$ in the phase-space. Integrating the ODEs between $t^{n}$ and
$t^{n+1}$ with a second-order quadrature in time, we obtain the following
implicit scheme,

$\mathrm{RK}_{\mathbf{F}}^{(2)}(\Delta
t;\mathbf{X})=\mathbf{F}\bigg{(}t^{n+1/2},\frac{\mathbf{X}+\mathbf{X}(t^{n+1})}{2}\bigg{)},$

where
$(\mathbf{X}(t^{n})+\mathbf{X}(t^{n+1}))/2\simeq\mathbf{X}(t^{n+1/2})+\mathcal{O}(\Delta
t^{2})$ is the second-order approximation of the midpoint characteristics. The
corresponding explicit scheme reads

$\mathrm{RK}_{\mathbf{F}}^{(2)}(\Delta
t;\mathbf{x}_{ij})=\mathbf{F}\Big{(}t^{n+1/2},\mathbf{x}_{ij}-\frac{\Delta
t}{2}\mathbf{F}_{ij}(t^{n})\Big{)}.$

Similar time discretizations already appeared in [102].

####  6.1.3 Fourth-order schemes (RK4)

Assume that for all times $t\in[t^{n},\,t^{n+1}]$, the approximate advection
field $\mathbf{F}(t)$ is known. This can be obtained by using the recursive
method of Section 3. We recall that this method allows us to compute
recursively the time derivatives of the advection field
$\\{\partial_{t}^{(k)}\mathbf{F}(t)\\}_{k\geq 0}$, in order to construct the
time-Taylor series expansion (36) of the advection field $\mathbf{F}(t)$. At
fourth order the Runge–Kutta scheme then reads

$\mathrm{RK}_{\mathbf{F}}^{(4)}(\Delta
t;\mathbf{X})=\frac{1}{6}(K_{1}+2K_{2}+2K_{3}+K_{4}),$

where the implicit advection estimates are given by

$\displaystyle K_{1}=\mathbf{F}(t^{n},\mathbf{X})$ (46a) $\displaystyle
K_{2}=\mathbf{F}\Big{(}t^{n+1/2},\mathbf{X}+\frac{\Delta t}{2}K_{1}\Big{)}$
(46b) $\displaystyle K_{3}=\mathbf{F}\Big{(}t^{n+1/2},\mathbf{X}+\frac{\Delta
t}{2}K_{2}\Big{)}$ (46c) $\displaystyle
K_{4}=\mathbf{F}(t^{n+1},\mathbf{X}+\Delta tK_{3}).$ (46d)

Since the coordinate argument of $K_{4}$ in (46d), namely
$\mathbf{X}^{n}+\Delta tK_{3}$, is supposed to be an estimate of the arrival
point at $t^{n+1}$, which coincides with the grid, we then replace (46d) by
$K_{4}=\mathbf{F}(t^{n+1},\mathbf{x}_{ij})$.
In its explicit form, the advection estimates take the following form,

$\displaystyle K_{1}=\mathbf{F}(t^{n},\mathbf{x}_{ij})$ (47) $\displaystyle
K_{2}=\mathbf{F}\Big{(}t^{n+1/2},\mathbf{x}_{ij}-\frac{\Delta
t}{2}K_{1}\Big{)}$ $\displaystyle
K_{3}=\mathbf{F}\Big{(}t^{n+1/2},\mathbf{x}_{ij}-\frac{\Delta
t}{2}K_{2}\Big{)}$ $\displaystyle
K_{4}=\mathbf{F}(t^{n+1},\mathbf{x}_{ij}-\Delta tK_{3}).$

Although it may seem more intuitive to have
$K_{1}=\mathbf{F}(t^{n+1},\mathbf{x}_{ij})$ and
$K_{4}=\mathbf{F}(t^{n},\mathbf{x}_{ij}-\Delta tK_{3})$ in the explicit scheme
(47), this leads to much poorer performance. The reason is most likely that
any initial inaccuracies in the estimate of
$\mathbf{F}(t^{n+1},\mathbf{x}_{ij})$ will affect all subsequent $K_{i}$.
Therefore, starting from the “exact” value
$K_{1}=\mathbf{F}(t^{n},\mathbf{x}_{ij})$ leads to an overall more accurate
scheme.

###### Remark 1

It should be noted that one can also construct a fully implicit Runge–Kutta
scheme of fourth order following the method by Hammer and Hollingsworth [56]
where

$\mathrm{RK}_{\mathbf{F}}^{(4)}(\Delta
t;\mathbf{X})=b_{1}\mathbf{F}(t^{n+1}-c_{1}\Delta t,\mathbf{x}_{ij}-\Delta
tK_{1})+b_{2}\mathbf{F}(t^{n+1}-c_{2}\Delta t,\mathbf{x}_{ij}-\Delta tK_{2}),$

and the fully implicit advection estimates take the following form,

$\displaystyle K_{1}=a_{11}\mathbf{F}(t^{n+1}-c_{1}\Delta
t,\mathbf{x}_{ij}-\Delta tK_{1})+a_{12}\mathbf{F}(t^{n+1}-c_{2}\Delta
t,\mathbf{x}_{ij}-\Delta tK_{2})$ $\displaystyle
K_{2}=a_{21}\mathbf{F}(t^{n+1}-c_{1}\Delta t,\mathbf{x}_{ij}-\Delta
tK_{1})+a_{22}\mathbf{F}(t^{n+1}-c_{2}\Delta t,\mathbf{x}_{ij}-\Delta
tK_{2}).$

Here, the coefficients $a_{ij}$ and $b_{i}$ are given by

$\displaystyle a_{11}=\frac{1}{4},\quad
a_{12}=\frac{1}{4}-\frac{\sqrt{3}}{6},\quad
a_{21}=\frac{1}{4}+\frac{\sqrt{3}}{6}\;,\quad a_{22}=\frac{1}{4},$
$\displaystyle b_{1}=b_{2}=\frac{1}{2},$ $\displaystyle
c_{1}=\frac{1}{2}-\frac{\sqrt{3}}{6},\quad
c_{2}=\frac{1}{2}+\frac{\sqrt{3}}{6}.$

We find that this fully implicit scheme offers a small performance increase
over the previous implicit scheme.

###  6.2 Adams–Moulton semi-Lagrangian schemes

Adams–Moulton schemes are a family of implicit schemes, which employ the
solutions at previous time steps to obtain a higher-order approximation for
the solution of ordinary differential equations. When applied to equation
(43), Adams–Moulton schemes require knowledge of the endpoint advection field
$\mathbf{F}^{n+1}$. To this end, a high-order reconstruction of
$\mathbf{F}^{n+1}$ can be obtained by using either an affine extrapolation in
time using the known advection field at previous times or by using the
recursice procedure of Section 3, which allows us to compute recursively the
time derivatives of the advection field, in order to construct the time-Taylor
series expansion of the advection field (36). As we will demonstrate, using
the recursive procedure of Section 3 generally produces more accurate results
than time extrapolations. Using the Adams–Moulton schemes of order $s$, that
we formally write as $\mathrm{AM}_{\mathbf{F}}^{(s)}(\Delta
t;t^{n},t^{n-1},\ldots,;\mathbf{X}(t^{n}))$, we obtain

$\displaystyle\mathbf{X}(t^{n+1})-\mathbf{X}(t^{n})$
$\displaystyle=\int_{t^{n}}^{t^{n+1}}\mathbf{F}(s,\mathbf{X}(s))\mathrm{d}s$
(48) $\displaystyle\simeq\Delta t\>\mathrm{AM}_{\mathbf{F}}^{(s)}\big{(}\Delta
t;t^{n+1},t^{n},t^{n-1},\ldots,;\mathbf{X}(t^{n})\big{)}.$

As before, we have to solve the following fixed-point problem to obtain the
departure points $\mathbf{X}(t^{n})$,

$\mathbf{X}\to G(\mathbf{X})\coloneqq\mathbf{x}_{ij}-\mathbf{X}-\Delta
t\>\mathrm{AM}_{\mathbf{F}}^{(s)}(\Delta
t;t^{n+1},t^{n},t^{n-1},\ldots;\mathbf{X}).$ (49)

Again, the problem $G(\mathbf{X})=0$ is solved by using the Picard iteration.

####  6.2.1 Second-order schemes (AM2)

Assume that for all times $t^{m}$, $m\leq n$, the approximate solutions of the
filling function $f^{m}$ and the corresponding advection field
$\mathbf{F}^{m}$ are known. The second-order Adams–Moulton scheme reads [46],

$\mathrm{AM}_{\mathbf{F}}^{(2)}(\Delta
t;\mathbf{X})=\frac{1}{2}\big{(}\mathbf{F}(t^{n+1},\mathbf{x}_{ij})+\mathbf{F}(t^{n},\mathbf{X})\big{)},$

where $\mathbf{F}(t^{n+1},\mathbf{x}_{ij})$ is obtained either via equation
(36), by using the known grid values
$\\{\partial_{t}^{(k)}\mathbf{F}_{ij}(t^{n})\\}_{k\geq 0}$, or through the
second-order extrapolation in time,

$\mathbf{F}(t^{n+1},\mathbf{x}_{ij})\simeq
2\mathbf{F}(t^{n},\mathbf{x}_{ij})-\mathbf{F}(t^{n-1},\mathbf{x}_{ij}).$

####  6.2.2 Fourth-order schemes (AM4)

Assume that for all times $t^{m}$, $m\leq n$, the approximate solutions of the
filling function $f^{m}$ and the corresponding advection field
$\mathbf{F}^{m}$ are known. The fourth-order Adams–Moulton scheme then reads
[46],

$\begin{dcases}\mathrm{AM}_{\mathbf{F}}^{(4)}(\mathbf{X})&=\frac{1}{24}\big{(}9\mathbf{F}(t^{n+1},\mathbf{x}_{ij})+19\mathbf{F}(t^{n},\mathbf{X})-5\mathbf{F}(t^{n-1},\mathbf{X}^{n-1})+\mathbf{F}(t^{n-2},\mathbf{X}^{n-2})\big{)}\\\
\mathbf{X}^{n-1}&=\mathbf{x}_{ij}+4(\mathbf{x}_{ij}-\mathbf{X})-2\Delta
t\big{(}2\mathbf{F}(t^{n},\mathbf{X})+\mathbf{F}(t^{n+1},\mathbf{x}_{ij})\big{)}\\\
\mathbf{X}^{n-2}&=\mathbf{x}_{ij}+27(\mathbf{x}_{ij}-\mathbf{X})-6\Delta
t\big{(}3\mathbf{F}(t^{n},\mathbf{X})+2\mathbf{F}(t^{n+1},\mathbf{x}_{ij})\big{)}.\end{dcases}$

The advection field $\mathbf{F}(t^{n+1},\mathbf{x}_{ij})$ is obtained either
via equation (36), by using the known grid values
$\\{\partial_{t}^{k}\mathbf{F}_{ij}(t^{n})\\}_{k\geq 0}$, or through the
fourth-order extrapolation in time,

$\mathbf{F}(t^{n+1},\mathbf{x}_{ij})\simeq
4\mathbf{F}(t^{n},\mathbf{x}_{ij})-6\mathbf{F}(t^{n-1},\mathbf{x}_{ij})+4\mathbf{F}(t^{n-2},\mathbf{x}_{ij})-\mathbf{F}(t^{n-3},\mathbf{x}_{ij}).$

###  6.3 Leap-frog schemes (LF2)

Unlike the Runge–Kutta schemes, leap-frog methods do not estimate the midpoint
advection field using a time-Taylor series expansion. Instead, they rely on
taking intermediate steps of the function $f$, from which the advection field
can be computed. Here, we present an implicit second-order leap-frog scheme
based on the semi-Lagrangian technique [46]. Assume that for all times
$t^{m}$, $m\leq n$, the approximate solutions of the filling function $f^{m}$
and the corresponding advection field $\mathbf{F}^{m}$ are known. First, we
estimate the midpoint solution $f^{n+1/2}$ by solving the fixed-point problem,

$\begin{dcases}\mathbf{X}^{n+1/2}-\mathbf{X}^{n}-\frac{\Delta
t}{2}\mathbf{F}(t^{n},\mathbf{X}^{n})=0\\\
\mathbf{X}^{n+1/2}=\mathbf{x}_{ij},\end{dcases}$

with the help of the Picard iteration and then interpolating the filling
function at $\mathbf{X}^{n}$ to the mesh points $\mathbf{x}_{ij}$, namely,

$f^{n+1/2}(\mathbf{x}_{i})=\Pi_{h}f^{n}(\mathbf{X}^{n}).$

Next, from $f^{n+1/2}$, the advection field $\mathbf{F}^{n+1/2}$ is
calculated, and then used to compute a new approximation of $\mathbf{X}^{n}$
by solving the fixed-point problem,

$\begin{dcases}\mathbf{X}^{n+1}-\mathbf{X}^{n}-\Delta
t\mathbf{F}\Big{(}t^{n+1/2},\frac{\mathbf{X}^{n+1}+\mathbf{X}^{n}}{2}\Big{)}=0\\\
\mathbf{X}^{n+1}=\mathbf{x}_{ij}.\end{dcases}$

Finally, the solution at time $t^{n+1}$ is given through interpolation, i.e.

$f^{n+1}(\mathbf{x}_{i})=\Pi_{h}f^{n}(\mathbf{X}^{n}).$

An explicit form of the second-order leap-frog scheme for solving the GHD
advection equation was proposed in [3]. This scheme requires knowledge of the
approximate solutions of the filling function both at time $t^{n}$ and
$t^{n+1/2}$. Given the latter, the advection field $\mathbf{F}^{n+1/2}$ is
calculated and the departure points $\mathbf{X}^{n}$ are obtained from the
following scheme,

$\begin{dcases}\mathbf{X}^{n}&=\mathbf{x}_{ij}-\Delta
t\mathbf{F}(t^{n+1/2},\mathbf{X}^{n+1/2})\\\
\mathbf{X}^{n+1/2}&=\mathbf{x}_{ij}-\frac{\Delta
t}{2}\mathbf{F}(t^{n+1/2},\mathbf{x}_{ij}).\end{dcases}$ (50)

Hence, the solution at time $t^{n+1}$ is given through interpolation by

$f^{n+1}(\mathbf{x}_{i})=\Pi_{h}f^{n}(\mathbf{X}^{n}).$ (51)

To obtain the solution at time $t^{n+3/2}$, required for the next iteration,
one must repeat the step described by equations (50) and (51) with the newly
obtained solution $f^{n+1}$ replacing $f^{n+1/2}$.

###  6.4 BSL-hybrid schemes for dissipative GHD

Given the departure points $\mathbf{X}^{n}$ found using BSL schemes for the
GHD advection equation (1), solving the extended GHD equation (21) for a
single time step requires integrating the nonlinear source term along the
characteristics, namely

$f^{n+1}(\mathbf{x}_{ij})=f^{n}(\mathbf{X}^{n})+\int_{t^{n}}^{t^{n+1}}\mathrm{d}\tau\>\mathcal{N}(\tau,\mathbf{X}(\tau),f(\tau)).$
(52)

To this end, various time discretization methods can be employed to create
BSL-hybrid schemes. In the following we will review a few such schemes.

####  6.4.1 The endpoint scheme

The simplest scheme consists of estimating the integral (52) by its endpoint,
giving

$f^{n+1}(\mathbf{x}_{ij})-f^{n}(\mathbf{X}^{n})=\Delta
t\mathcal{N}(t^{n},\mathbf{X}^{n}).$

####  6.4.2 The midpoint scheme

Alternatively, one can employ the source term at the midpoint of the time step
to obtain,

$f^{n+1}(\mathbf{x}_{ij})-f^{n}(\mathbf{X}^{n})=\frac{\Delta
t}{2}\big{(}\mathcal{N}(t^{n+1/2},\mathbf{x}_{ij})+\mathcal{N}(t^{n+1/2},\mathbf{X}^{n})\big{)}.$
(53)

Here, the midpoint along the trajectory is obtained by averaging the source
term at its start and endpoint. This procedure improves the stability, yet
retains the second-order accuracy, of the more standard scheme, which consists
in evaluating $\mathcal{N}(t^{n+1/2})$ at the middle of the trajectory
$(\mathbf{x}_{ij}+\mathbf{X}^{n})/2$, i.e. by substituting $\Delta
t\mathcal{N}(t^{n+1/2},(\mathbf{x}_{ij}+\mathbf{X}^{n})/2)$ to right-hand side
of (53) (see [59]). The midpoint in time is approximated by linear
extrapolation from the previous time step such as

$\mathcal{N}(t^{n+1/2},\mathbf{X})\simeq\frac{3}{2}\mathcal{N}(t^{n},\mathbf{X})-\frac{1}{2}\mathcal{N}(t^{n-1},\mathbf{X}).$

Higher-order estimates of the temporal midpoint can be obtained by including
additional previous time steps, however, such extrapolations may also be more
prone to noise [59].

####  6.4.3 The two-time-level semi-Lagrangian scheme (SETTLS)

The extrapolation of the source term employed in the midpoint scheme can
suffer from noise problems. An alternative scheme is the SETTLS scheme [44,
59, 85], given by

$f^{n+1}(\mathbf{x}_{ij})-f^{n}(\mathbf{X}^{n})=\frac{\Delta
t}{2}\big{(}\mathcal{N}(t^{n},\mathbf{x}_{ij})+\mathcal{N}(t^{n+1},\mathbf{X}^{n})\big{)}.$

Here, $\mathcal{N}(t^{n+1})$ is again approximated by linear extrapolation
from the previous time step such as

$\mathcal{N}(t^{n+1},\mathbf{X})\simeq
2\mathcal{N}(t^{n},\mathbf{X})-\mathcal{N}(t^{n-1},\mathbf{X}).$

####  6.4.4 The Crank–Nicolson scheme

Applying a trapezoidal quadrature rule to estimate the integral (52) yields
[97],

$f^{n+1}(\mathbf{x}_{ij})-f^{n}(\mathbf{X}^{n})=\frac{\Delta
t}{2}\big{(}\mathcal{N}(t^{n+1},\mathbf{x}_{ij})+\mathcal{N}(t^{n},\mathbf{X}^{n})\big{)}.$
(54)

Unlike the previous methods, the Crank–Nicolson method is implicit as the
source term $\mathcal{N}(t^{n+1},\mathbf{x}_{ij})$ depends on the unknown
$f^{n+1}(\mathbf{x}_{ij})$. The problem (54) can be solved by using the Picard
iteration. Note that repeatedly evaluating the source term is extremely
costly. Therefore, in the benchmarks below, the tolerance for the convergence
of the Picard iteration is set so that the scheme converges in few iterations,
typically two or three.

##  7 Benchmark of BSL and BSL-hybrid schemes for solving the dissipative GHD
equations

To benchmark the performance of the various schemes, we will examine the
evolution of three particular conserved quantities, namely the particle number
$N$, the energy $E$, and the entropy $S$ given by the expressions

$\displaystyle N(t)$
$\displaystyle=\frac{1}{2\pi}\int_{\mathbb{R}}\mathrm{d}z\int_{\mathbb{R}}\mathrm{d}\theta\>1^{\mathrm{dr}}(t,z,\theta)f(t,z,\theta),$
(55) $\displaystyle E(t)$
$\displaystyle=\frac{1}{2\pi}\int_{\mathbb{R}}\mathrm{d}z\int_{\mathbb{R}}\mathrm{d}\theta\>\big{(}\theta^{2}+V(z)\big{)}1^{\mathrm{dr}}(t,z,\theta)f(t,z,\theta),$
(56) $\displaystyle S(t)$
$\displaystyle=\frac{1}{2\pi}\int_{\mathbb{R}}\mathrm{d}z\int_{\mathbb{R}}\mathrm{d}\theta\>1^{\mathrm{dr}}(t,z,\theta)f^{2}(t,z,\theta).$
(57)

Note, the latter quantity corresponds to the entropy defined by (17) with
$\Psi(f)=f$ and is particularly sensitive to (numerical) diffusion. To
quantify the accuracy of each scheme, we compute the relative error of the
conserved quantities w.r.t. their initial value, namely

$\mathcal{E}_{Q}(t)=\frac{|Q(t)-Q(0)|}{Q(0)},$

where $Q(t)=\\{N(t),E(t),S(t)\\}$.

For the benchmark, we will solve the ODE (43) for a relatively difficult
problem, namely the Quantum Newton’s cradle setup [62, 68, 78]. The setup can
be realized experimentally using ultracold atomic gases and constitutes one of
the most remarkable demonstrations of integrability. In the Quantum Newton’s
cradle setup, the system is initialized in a (nearly) harmonic, external
potential $V(z)$. Next, two large, opposite momenta are imparted onto the
system, causing the initial (typically thermal) state to split into two atomic
clouds oscillating in the external confinement. Every half-period of the
cradle oscillation the clouds collide in the center of the potential. The
ensuing interactions, encoded in the effective velocity, deform the clouds.
Further, a small anharmonicity of the external potential can lead to the
formation of very fine structures in the filling function, thus demanding high
accuracy of the numerical solutions [26]. An example of this can be seen in
Fig. 3.

In GHD simulations of experimental systems, the initial filling function
$f_{0}$ is often given by a thermal state obtained via the thermodynamic Bethe
Ansatz [104]. However, for easier reproducibility, we instead let it be given
by the expression

$f_{0}(z,\theta)=0.9\left[\exp\left(-\frac{\left(\theta-\theta_{\mathrm{B}}\right)^{2}}{2\sigma}\right)+\exp\left(-\frac{\left(\theta+\theta_{\mathrm{B}}\right)^{2}}{2\sigma}\right)\right]\exp\left(-\frac{z^{2}}{2\sigma}\right),$
(58)

where $\theta_{\mathrm{{B}}}=2$ and $\sigma=1/\sqrt{2}$. Further, we employ an
external potential in the shape of a Gaussian, imitating the intensity profile
of a laser beam typically used to generate the potential in an experimental
setting

$V(z)=\frac{1}{8}\omega^{2}\eta^{2}\left[1-\exp\left(-\frac{2z^{2}}{\eta^{2}}\right)\right].$

For the benchmarks we set $\omega=2$ and $\eta=12$. The coupling constant is
set to $c=1$. We will be expressing time in units of the oscillation period
$\mathcal{T}=2\pi/\omega$, and the simulations will be carried out for a
duration of 10 periods. Whenever the time-Taylor expansion (36) is used to
estimate the advection fields at future times, only the first order time
derivatives of the advection fields are employed. All numerical schemes for
solving the GHD equations are implemented using the iFluid framework [80].

Figure 3:  Example of filling function $f$ calculated for the quantum Newton’s
cradle benchmark problem. The initial state (58) is propagated for 10
oscillation periods $\mathcal{T}$ of the cradle following BSL (hybrid)
methods. Top row: Purely advective GHD calculation performed using explicit
Runge–Kutta scheme of fourth order (RK4) with equally-spaced grids
$N_{\theta}=N_{z}=513$ and $N_{\mathrm{steps}}=2000$ time steps. Note all
solutions for $f$ obtained using higher-order methods appear visually very
similar. Middle and bottom rows: Dissipative GHD calculations performed using
BSL hybrid schemes (SETTLS + explicit RK4) with equally-spaced grids
$N_{\theta}=N_{z}=257$ and $N_{\mathrm{steps}}=3000$ time steps. The
dissipative GHD equation is solved for the diffusion operator (middle row) and
the Boltzmann-like collision integral (bottom row). Note the grids employed
for the computation extend further than shown in the figure. Further, all
solutions obtained for $f$ appear visually very similar.

###  7.1 BSL schemes for solving the advective GHD equation

For the benchmarks of the BSL schemes for solving the advective GHD equation,
we employ uniformly spaced position and rapidity grids discretized into
$N_{\theta}=N_{z}=513$ grid points. Such high resolution is necessary to
properly resolve the fine structures appearing in the filling function $f$
following dynamics (see Fig. 3). We test the accuracy of the various schemes
for different time step resolutions, discretizing the evolution duration into
$N_{\mathrm{steps}}=500$, 1000, and 2000 time steps.

However, before comparing their accuracy, it is worth discussing the numerical
cost of performing a single time step following the different schemes. Among
all operations needed for performing a BSL time step of the GHD equation, the
most numerically expensive one is, by far, solving the dressing equation (38).
Calculating the effective velocity requires solving the dressing equation for
both the rapidity derivative of the single-particle energy
$\partial_{\theta}\epsilon(\theta)$ and momentum $\partial_{\theta}p(\theta)$.
Further, if the time-Taylor series expansion (36) is used to compute the
advection fields at future times, additional dressing equations must be
solved, as calculating the (first order) time derivatives of the advection
fields requires dressing the source terms $S_{w}$ and $S_{b}$. Note, for the
Lieb-Liniger model whose single-particle momentum reads $p(\theta)=\theta$,
solving the dressing equation for the unit-cell $1^{\mathrm{dr}}$ is not
necessary, as it is already given by the dressed rapidity derivative of the
single-particle momentum.

When timing the different schemes, we find that the total cost of performing a
single BSL step is indeed completely dominated by solving the dressing
equation. In the implicit schemes, the advection fields (and their time
derivatives) are calculated only once at each propagation step and then
interpolated to the solution $\mathbf{X}^{n}$ at each iteration. In comparison
to computing the fields, performing the phase space interpolation is a
relatively cheap operation, resulting in a rather small difference in runtime
between implicit and explicit schemes. Restricting ourselves to first-order
time derivatives of the advection fields, we find the propagation step of
higher-order Runge–Kutta schemes (RK2 and RK4) taking roughly twice the time
to execute than the first-order scheme (RK1), as the latter does not employ
estimates of the advection fields at intermediate times. Meanwhile, the leap-
frog schemes perform an additional, intermediate step for every full time
step, thus having to evaluate the effective velocity twice. Hence, for a full
$\Delta t$ time step, the leap-frog schemes exhibit computation times similar
to the higher-order Runge–Kutta schemes. Lastly, the Adams–Moulton schemes
employing the time-Taylor series expansion show similar computational cost to
the higher-order Runge–Kutta schemes. However, when employing an extrapolation
of the advection fields using previous steps, the Adams–Moulton schemes obtain
a numerical efficiency compared to the RK1 schemes, since only a single
evaluation of the advection fields is needed.

Note that all these observations are specific to the problem of solving the
GHD equation. In most other problems, computing the advection fields does not
represent such a significant fraction of the computational load, whereby a
discernible difference in the runtime between implicit and explicit schemes
would be found. Furthermore, for problems involving higher-dimensional
systems, storing the full solution at multiple time-steps (needed for the
Adams–Moulton schemes) may introduce significant overhead.

Figure 4:  Benchmark of first-order BSL schemes for the quantum Newton’s
cradle problem. In the legend, the labels ”imp” and ”exp” indicate that the
methods are implicit or explicit, respectively. The calculation is performed
on equally-spaced grids with resolutions $N_{\theta}=N_{z}=513$ and for
various number of time steps $N_{\mathrm{steps}}$. Plotted is the relative
error of the particles number, total energy, and entropy.  Figure 5:
Benchmark of second-order BSL schemes for the quantum Newton’s cradle problem.
In the legend, the labels ”imp” and ”exp” indicate that the methods are
implicit or explicit, respectively. For the Adams–Moulton schemes (AM), the
label ”ext” denotes that the advection fields at the endpoint are obtained
using extrapolation from previous time steps. Meanwhile, ”T.e” indicates that
the time-Taylor expansion series has been employed. The calculation is
performed on equally-spaced grids with resolutions $N_{\theta}=N_{z}=513$ and
for various number of time steps $N_{\mathrm{steps}}$. Plotted is the relative
error of the particles number, total energy, and entropy.  Figure 6:
Benchmark of fourth-order BSL schemes for the quantum Newton’s cradle problem.
In the legend, the labels ”imp” and ”exp” indicate that the methods are
implicit or explicit, respectively. For the Adams–Moulton schemes (AM), the
label ”ext” denotes that the advection fields at the endpoint are obtained
using extrapolation from previous time steps. Meanwhile, ”T.e” indicates that
the time-Taylor expansion series has been employed. The calculation is
performed on equally-spaced grids with resolutions $N_{\theta}=N_{z}=513$ and
for various number of time steps $N_{\mathrm{steps}}$. Plotted is the relative
error of the particles number, total energy, and entropy.

In Figures 4, 5, and 6 we plot the benchmark results for the first order,
second order, and fourth order BSL schemes, respectively.

Starting with the first order schemes of Figure 4, we find that the explicit
form of the RK1 algorithm performs extremely poorly, which is, perhaps, to be
expected. In comparison, the implicit form yields much more accurate results,
although the relative change in the conserved quantities is still fairly
large. For the implicit scheme, we find the error of both the particle number
and the entropy to exhibit rapid oscillations at a frequency corresponding to
roughly half the period of the quantum Newton’s cradle. The oscillation
frequency coincides with the overlap of the two atomic clouds in the center of
the external trap, where interactions are much more significant than
throughout the rest of the evolution. Therefore, as the clouds overlap in
position, the numerical solution becomes more sensitive to errors and
additional precision is often needed.

Compared to the first order schemes, the second order schemes, whose
benchmarks are plotted in Figure 5, are much more accurate. Comparing the
implicit versus explicit schemes for the Runge–Kutta and leap-frog methods, we
generally find the implicit schemes to perform better. In particular the
implicit leap-frog scheme exhibits very low relative errors of all three
conserved quantities calculated. Furthermore, the magnitudes of its errors
remain fairly stable throughout the whole evolution duration, as opposed to
most of the other schemes, whose associated errors increase steadily. Lastly,
both of the two Adams–Moulton schemes exhibit worse performance than the
Runge–Kutta schemes of same order. Between the two of them, the scheme
employing the time-Taylor series expansion (36), as opposed to the linear
extrapolation from previous time steps, to obtain the advection fields at time
$t^{n+1}$ is more accurate. However, when solving the GHD equation, the latter
is much cheaper to perform numerically. Interestingly, all the implicit
schemes tested exhibit the same oscillating behaviour in the errors of the
benchmark quantities, while the explicit schemes appear much more stable in
this regard.

Finally, in Figure 6 we compare the different fourth-order schemes. In
general, all schemes exhibit very high accuracy. For increasing number of time
steps, we observe increasing accuracy and thus convergence of the schemes. At
fixed number of time steps, we observe that both of the Runge–Kutta schemes
perform better than the Adams–Moulton schemes. Contrary to the second-order
schemes, the explicit Runge–Kutta scheme now performs better than its implicit
counterpart, and extrapolating of the advection fields from previous time
steps yields better results for the the Adams–Moulton scheme than employing
the time-Taylor series expansion (36). The improvement of the explicit scheme
is most likely due to the trick that we discovered, where the initial
advection estimate of the Runge–Kutta method $K_{1}$ is computed using the
known state at time $t^{n}$. Meanwhile, the worse performance of the
Adams–Moulton scheme with time-Taylor series expanded advection fields might
be due to the limited phase-space resolution. Indeed, when computing the
successive time derivatives of the advection field $\mathbf{F}$ according to
the recursive procedure of Section 3, derivatives in the position and rapidity
variables of the filling function $f$ are required. Here, the phase-space
derivatives of $f$ are obtained through finite-difference methods. At high
enough time resolution, the limiting factor of the time-Taylor series
expansion scheme may become the error of said finite-difference derivatives
together with the limited phase-space resolution of the grid.

###  7.2 BSL-hybrid schemes for solving the dissipative GHD equations

For the benchmarks of the BSL-hybrid schemes for solving the dissipative GHD
equations, we employ uniformly spaced position and rapidity grids discretized
into $N_{\theta}=N_{z}=257$ grid points. Compared to the purely advective
case, coarser grids can typically be employed when accounting for dissipation,
as these effects tend to smooth out the filling function $f$ following
dynamics. We test two dissipative nonlinear source terms, namely the diffusion
operator (22) and the Boltzmann-type collision integral (27) governing
transverse state-changing collisions. The two source terms have very different
structures. The diffusion operator scales with the position gradient of the
filling function, thus leading to a smearing of its edges. Meanwhile, the
transverse excitations effectively destroy quasi-particles at high rapidities
and re-create them at low rapidities, while de-excitations from transverse
excited states lead to a general redistribution of quasi-particles across the
($z,\theta$) phase-space. Hence, although both source terms lead to eventual
thermalization of the system, the exact mechanisms, and thus dynamics towards
thermalization, differ. Examples of this can be seen in Fig. 3, where
solutions for the filling function $f$ for the two different source-terms are
plotted.

We test the accuracy of the various schemes for different time step
resolutions, discretizing the evolution duration into
$N_{\mathrm{steps}}=2000$, 2500, and 3000 time steps. To compute the departure
points of the characteristics, the explicit form of the Runge–Kutta scheme of
fourth order (RK4) is employed. In addition to computing the error of the
particle number, total energy, and entropy, we also calculate and plot the
minimum value of the filling function $f$.

Compared to the pure BSL schemes, each time step of the BSL-hybrid schemes is
much more numerically costly to perform due to the evaluation of the source
terms. Hence, the endpoint, midpoint and SETTLS schemes all exhibit very
similar runtimes, as the computation is completely dominated by the single
calculation of the source term needed. Conversely, following the Crank-
Nicolson scheme, the source term must be calculated anew at every iteration,
making this scheme extremely expensive compared to the others. In order to
keep the computational runtime at a reasonable limit, we set the limit of
convergence for the Crank-Nicolson scheme at a level attainable after just 2-3
iterations.

####  7.2.1 Diffusion

Figure 7:  Benchmark of BSL-hybrid schemes for solving the dissipative GHD
equation, here with the diffusion kernel as the nonlinear dissipative term.
The benchmarks are performed on quantum Newton’s cradle problem using the
explicit RK4 method. The calculation is performed on equally-spaced grids with
resolutions $N_{\theta}=N_{z}=257$ and for various number of time steps
$N_{\mathrm{steps}}$. Plotted is the relative error of the particle number,
total energy, and entropy, as well as the negative minimum value of the
filling function $f$.

Figure 7 shows the benchmark results of solving the quantum Newton’s cradle
problem accounting for diffusive effects. Comparing to the previous benchmarks
of the explicit RK4 method for solving the purely advective GHD equation, we
do not find any significant increase in the relative error of neither particle
number $N$ nor total energy $E$ when accounting for diffusion. Meanwhile, the
entropy increases significantly following diffusive dynamics, as expected.
Lastly, we do observe that the filling function $f$ at some points of the
phase-space assumes negative values. Recall that the filling function
represents the fraction of allowed local momentum states occupied. Hence, it
should only assume values between 0 and 1. However, the magnitude of the
negative values found following numerical propagation is very small. Comparing
the different hybrid schemes, not one single scheme consistently stands out in
terms of accuracy, although for diffusive dynamics the midpoint scheme
consistently appears slightly worse than the others.

####  7.2.2 Boltzmann-type collision integral

Figure 8:  Benchmark of BSL-hybrid schemes for solving the dissipative GHD
equation, here with the Boltzmann-type collision integral for transverse
state-changing collisions as the nonlinear dissipative term. The benchmarks
are performed on quantum Newton’s cradle problem using the explicit RK4 method
to obtain estimates of characteristics. The calculation is performed on
equally-spaced grids with resolutions $N_{\theta}=N_{z}=257$ and for various
number of time steps $N_{\mathrm{steps}}$. Plotted is the relative error of
the particle number, total energy, and entropy, as well as the negative
minimum value of the filling function $f$.

Unlike the case of diffusion (and most other instances of Boltzmann-type
collision integrals), accounting for transverse excitations through a
collision integral requires solving an additional set of differential
equations, namely the rate equations for the excitation probabilities
$\nu_{n}$ (28). For the endpoint scheme, eq. (28) is numerically solved as
follows

$\nu_{n}(t^{n+1})=\nu_{n}(t^{n})+\Delta
t\frac{1}{2}\beta_{n}\big{[}\Gamma_{\mathrm{h}}^{+}(t^{n})-\Gamma_{\mathrm{p}}^{+}(t^{n})\nu_{n}^{\beta_{n}}(t^{n})\big{]}\;,$
(59)

while the remaining schemes employ midpoint estimates, such that

$\nu_{n}(t^{n+1})=\nu_{n}(t^{n})+\Delta
t\frac{1}{2}\beta_{n}\big{[}\Gamma_{\mathrm{h}}^{+}(t^{n+1/2})-\Gamma_{\mathrm{p}}^{+}(t^{n+1/2})\>\nu_{n}^{\beta_{n}}(t^{n+1/2})\big{]},$
(60)

where the function values at time $t^{n+1/2}=t^{n}+\Delta t/2\,$ are obtained
by using the time extrapolation scheme

$\Gamma_{\mathrm{h}}^{+}(t^{n+1/2})=\frac{3}{2}\Gamma_{\mathrm{h}}^{+}(t^{n})-\frac{1}{2}\Gamma_{\mathrm{h}}^{+}(t^{n-1})\;.$
(61)

When accounting for atoms excited to transverse states beyond the ground
state, we must also consider the associated transverse potential energy of the
atoms can computing the total energy (56). An atom in the $n$’th transverse
state (where $n=0$ denotes the transverse ground state), has the transverse
potential energy $2nl_{\perp}^{-2}$ (in units $\hbar=2m=1$). Thus, we can
account for the transverse potential energy by adding the term
$\sum_{n}\nu_{n}2nl_{\perp}^{-2}$ to the external potential $V(z)$ of eq.
(56). For these benchmarks we set the parameter $l_{\perp}=1$.

Because the collision integral is originally constructed in the conservative
form, we here employ a slightly different BSL-hybrid scheme: First, the
advective GHD equation is solved using a BSL scheme to obtain the
characteristic origin. Next, the filling function $f$ is transformed into the
quasi-particle distribution $\rho_{\mathrm{p}}$, from which the collision
integral is computed and accounted for. Finally, the resulting quasi-particle
distribution is transformed back to the filling function, marking the end of a
single propagation step.

Figure 8 shows the benchmark results of solving the quantum Newton’s cradle
problem accounting for transverse state-changing collisions via the collision
integral. Comparing to the previous benchmarks of the explicit RK4 method for
solving the purely advective GHD equation, we do not find any significant
increase in the relative error of neither particle number $N$ when accounting
for transverse state-changing collisions. Meanwhile, for all the schemes we
find a very similar increase in the error of the total energy over time. As
already discussed, simplifications made in the construction of the collision
integral result in the energy no longer being explicitly conserved. The
similarity of the error across all the tested schemes indicates that the
source of the error is indeed due to the model, not the numerical solution.
Furthermore, we argued that the error associated with the simplifications
would be rather small in most cases. Indeed, the relative error of the total
energy seen in Figure 8 is only on the order of $0.1\%$ following 10
oscillation periods of the Newton’s cradle. Lastly, we again observe an
increase in the error of entropy, as expected, and the filling function $f$
assuming negative values at some points of the $(z,\theta)$ phase-space.
However, again the magnitude of the negative values is very small for all the
schemes.

##  8 Conclusion and outlook

We present dissipative Generalized HydroDynamic (GHD) equations as a kinetic
nonlinear PDE proven to be very useful to model accurately the dynamics
observed in quantum experiments for quasi one-dimensional systems of cold
atoms. Our main goal is to present, benchmark and discuss numerical methods
for solving the GHD equations with a diffusion term or a Boltzmann type
collission term. In particular, we design novel backward semi-Lagrangian (BSL)
schemes of high-order, namely implicit/explicit Runge–Kutta BSL schemes.
Despite a growing interest for the (dissipative) GHD equations, accompanied by
an abundant literature on the subject, it is the first time that these
equations are solved with such high-order numerical schemes. In order to prove
the accuracy of these new numerical methods, we perform some comparisons with
other schemes in the literature. Numerical simulations of the example of
”quantum Newton’s cradle” setup show that implicit/explicit Runge–Kutta BSL
schemes preserve a posteriori some integral invariants (number of quasi-
particles, energy, entropies) of the system in an accurate and satisfactory
manner. These observations prevail as well in the case of numerical methods
including the dissipative terms as without them. Finally, our numerical
simulations of the dissipative GHD equations are consistent with the
thermalization processes observed in experimental quasi one-dimensional
devices initialized with a quantum Newton’s cradle setup [78].

There are several future extensions of our work: First we can develop other
high-order numerical schemes for solving the dressing operation, such as the
continuous-Galerkin or discontinuous-Galerkin methods, but also some spectral
(Fourier-like) methods to obtain a faster dressing transformation. Second we
can devise other higher-order in time BSL-hybrid schemes for solving the
dissipative GHD equations by using the strategy of exponential integrators
[14, 32, 29, 58, 100, 30, 51, 73, 85]. We can also design higher-order schemes
for the advection (Vlasov) operator by considering discontinuous-Galerkin
schemes [31, 57] or, even better, semi-Lagrangian discontinuous-Galerkin
methods [91, 87, 10]. In addition, we can consider other high-order and fast
methods such as spectral methods (Fourier–Galerkin, discontinuous-Galerkin) or
real methods (discrete-velocity) for solving the Boltzmann-like collision and
diffusion operators [38, 81, 60, 48, 84]. Finally, it would be interesting to
prove a $H$-theorem for the Boltzmann-type collision integral introduced in
Section 2.2.2.

## Appendix

##  Appendix A Some properties of the GHD equations

In this appendix, we show some properties of the GHD equations written either
in conservative form through equations (8)-(10) or in advective form through
equations (1), (2) and (3).

###  A.1 Equation for the density of states $\rho_{\mathrm{s}}$

Here, we show that the density of states $\rho_{\mathrm{s}}$, defined by (11),
satisfies the same equation as $\rho_{\mathrm{p}}$, i.e.,

$\partial_{t}\rho_{\mathrm{s}}+\partial_{z}\big{(}v^{\mathrm{eff}}\,\rho_{\mathrm{s}}\big{)}+\partial_{\theta}\big{(}a^{\mathrm{eff}}\,\rho_{\mathrm{s}}\big{)}=0.$
(62)

Using definition (11) for $\rho_{\mathrm{s}}$ and definition (4) for the
dressing operation, the density of states $\rho_{\mathrm{s}}$ rewrites as

$\rho_{\mathrm{s}}(\theta)=\frac{1}{2\pi}+\int_{\mathbb{R}}\mathrm{d}\theta^{\prime}\,T(\theta-\theta^{\prime})\rho_{\mathrm{p}}(\theta^{\prime})=\frac{1}{2\pi}+T\ast\rho_{\mathrm{p}},$
(63)

where the symbol $\ast$ denotes the standard convolution in the rapidity
variable $\theta$. Using (63) equations (9)-(10) rewrites as

$\displaystyle\rho_{\mathrm{s}}v^{\mathrm{eff}}=\frac{\partial_{\theta}\epsilon}{2\pi}+T\ast(v^{\mathrm{eff}}\rho_{\mathrm{p}}),$
(64)
$\displaystyle\rho_{\mathrm{s}}a^{\mathrm{eff}}=-\frac{\partial_{z}\epsilon}{2\pi}+T\ast(a^{\mathrm{eff}}\rho_{\mathrm{p}}).$
(65)

Using (63)-(65) and the properties of the convolution operation, we obtain
from equation (8) for the density of quasi-particle $\rho_{\mathrm{p}}$,

$\displaystyle 0$ $\displaystyle=$ $\displaystyle
T\ast[\partial_{t}\rho_{\mathrm{p}}+\partial_{z}(v^{\mathrm{eff}}\rho_{\mathrm{p}})+\partial_{\theta}(a^{\mathrm{eff}}\rho_{\mathrm{p}})]$
$\displaystyle=$
$\displaystyle\partial_{t}(T\ast\rho_{\mathrm{p}})+\partial_{z}(T\ast[v^{\mathrm{eff}}\rho_{\mathrm{p}}])+\partial_{\theta}(T\ast[a^{\mathrm{eff}}\rho_{\mathrm{p}}])$
$\displaystyle=$
$\displaystyle\partial_{t}(\rho_{\mathrm{s}})+\partial_{z}(\rho_{\mathrm{s}}v^{\mathrm{eff}}-\partial_{\theta}\epsilon/(2\pi))+\partial_{\theta}(\rho_{\mathrm{s}}a^{\mathrm{eff}}+\partial_{z}\epsilon/(2\pi))$
$\displaystyle=$
$\displaystyle\partial_{t}(\rho_{\mathrm{s}})+\partial_{z}(\rho_{\mathrm{s}}v^{\mathrm{eff}})+\partial_{\theta}(\rho_{\mathrm{s}}a^{\mathrm{eff}}),$

which is equation (62).

###  A.2 Integral invariants

Here, we show that the GHD equations preserve in time the total number of
quasi-particles $N$ defined by (14), the total density of states $J$ defined
by (15), the total energy $E$ defined by (16), and the infinite number of
entropies $S$ defined by (16). We start with $N$ and $J$. Time invariance of
$N$ and $J$ is obvious from integration in the phase space $(z,\theta)$ of
respectively the conservation law (8) for $\rho_{\mathrm{p}}$ and the
conservation law (13) for $\rho_{\mathrm{s}}$. For the time invariance of $E$,
we first observe that for any functions $\theta\mapsto g\in
L^{2}({\mathbb{R}})$ and $\theta\mapsto h\in L^{2}({\mathbb{R}})$, we have

$\langle(1-\widehat{T}f)g,h\rangle=\langle g,(1-f\widehat{T})h\rangle,$ (66)

where

$\langle g,h\rangle=\int_{\mathbb{R}}d\theta\,g(\theta)h(\theta),$

is the standard scalar product in $L^{2}({\mathbb{R}})$. In other words, (66)
means that $(1-\widehat{T}f)^{T}=(1-f\widehat{T})$, where the operator
$\widehat{K}^{T}$ is the transpose or adjoint operator of the operator
$\widehat{K}$ induced by the scalar product $\langle\cdot,\cdot\rangle$.
Indeed relation (66) is a consequence of the following property of the
convolution operator, $\langle g\ast T,h\rangle=\langle g,T\ast h\rangle$,
with $T$ an even function. Using (66) and definition (7) of the dressing
operator, we obtain

$\langle fg^{\mathrm{dr}},h\rangle=\langle g,fh^{\mathrm{dr}}\rangle.$ (67)

Indeed, using (66) and definition (7), we obtain

$\displaystyle\langle fg^{\mathrm{dr}},h\rangle$ $\displaystyle=$
$\displaystyle\langle
f(1-\widehat{T}f)^{-1}g,h\rangle=\langle(1-\widehat{T}f)^{-1}g,fh\rangle=\langle
g,(1-\widehat{T}f)^{-T}fh\rangle=\langle g,(1-f\widehat{T})^{-1}fh\rangle$
$\displaystyle=$ $\displaystyle\langle
g,[f(1-\widehat{T}f)(1/f)]^{-1}fh\rangle=\langle
g,f(1-\widehat{T}f)^{-1}h\rangle=\langle g,fh^{\mathrm{dr}}\rangle.$

Now, using integration by parts, the relation
$\rho_{\mathrm{p}}=\rho_{\mathrm{s}}f$, the identity
$\rho_{\mathrm{s}}=1^{\mathrm{dr}}/(2\pi)$, the definitions (2)-(3), and the
property (67), we obtain from he conservation law (8) for $\rho_{\mathrm{p}}$,

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}E(t)$ $\displaystyle=$
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\Big{(}\int_{\mathbb{R}}\mathrm{d}z\int_{\mathbb{R}}\mathrm{d}\theta\,\rho_{\mathrm{p}}\epsilon\Big{)}$
$\displaystyle=$
$\displaystyle\int_{\mathbb{R}}\mathrm{d}z\int_{\mathbb{R}}\mathrm{d}\theta\,\partial_{t}\rho_{\mathrm{p}}\epsilon$
$\displaystyle=$
$\displaystyle-\int_{\mathbb{R}}\mathrm{d}z\int_{\mathbb{R}}\mathrm{d}\theta\,\\{\partial_{z}(\rho_{\mathrm{p}}v^{\mathrm{eff}})+\partial_{\theta}(\rho_{\mathrm{p}}a^{\mathrm{eff}})\\}\epsilon$
$\displaystyle=$
$\displaystyle\int_{\mathbb{R}}\mathrm{d}z\int_{\mathbb{R}}\mathrm{d}\theta\,\\{\rho_{\mathrm{p}}v^{\mathrm{eff}}\partial_{z}\epsilon+\rho_{\mathrm{p}}a^{\mathrm{eff}}\partial_{\theta}\epsilon\\}$
$\displaystyle=$
$\displaystyle\frac{1}{2\pi}\int_{\mathbb{R}}\mathrm{d}z\int_{\mathbb{R}}\mathrm{d}\theta\,f\\{1^{\mathrm{dr}}v^{\mathrm{eff}}\partial_{z}\epsilon+1^{\mathrm{dr}}a^{\mathrm{eff}}\partial_{\theta}\epsilon\\}$
$\displaystyle=$
$\displaystyle\frac{1}{2\pi}\int_{\mathbb{R}}\mathrm{d}z\int_{\mathbb{R}}\mathrm{d}\theta\,f\\{(\partial_{\theta}\epsilon)^{\mathrm{dr}}\partial_{z}\epsilon+(-\partial_{z}\epsilon)^{\mathrm{dr}}\partial_{\theta}\epsilon\\}$
$\displaystyle=$ $\displaystyle 0.$

For the time invariance of $S$, using an integration by parts, and equations
(1) and (8), we obtain

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}S(t)$ $\displaystyle=$
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\Big{(}\int_{\mathbb{R}}\mathrm{d}z\int_{\mathbb{R}}\mathrm{d}\theta\,\rho_{\mathrm{p}}\Psi(f)\Big{)}$
$\displaystyle=$
$\displaystyle\int_{\mathbb{R}}\mathrm{d}z\int_{\mathbb{R}}\mathrm{d}\theta\,\\{\Psi(f)\partial_{t}\rho_{\mathrm{p}}+\rho_{\mathrm{p}}\Psi^{\prime}(f)\partial_{t}f$
$\displaystyle=$
$\displaystyle-\int_{\mathbb{R}}\mathrm{d}z\int_{\mathbb{R}}\mathrm{d}\theta\,\\{\Psi(f)[\partial_{z}(v^{\mathrm{eff}}\,\rho_{\mathrm{p}})+\partial_{\theta}(a^{\mathrm{eff}}\,\rho_{\mathrm{p}})]+\rho_{\mathrm{p}}\Psi^{\prime}(f)[v^{\mathrm{eff}}\partial_{z}f+a^{\mathrm{eff}}\partial_{\theta}f]\\}$
$\displaystyle=$ $\displaystyle 0.$

##  Appendix B Integral invariants for the dissipative GHD equations with the
diffusion operator

In this appendix, we prove that the dissipative GHD equations, where the
right-hand side is given by the diffusion operator viz. equations (25)-(26),
preserves the total number of quasi-particles $N$ defined by (14), and the
total energy $E$ defined by (16). Since in Appendix A.2 we have already shown
the conservation of the total number of quasi-particle and of the total energy
for the transport part, i.e. for the left-hand side of equation (25), it only
remains to show that these quantities are preserved by the diffusion operator
(26). From the conservative form of the diffusion operator (26), an
integration of the latter in the phase space $(z,\theta)$ gives
straightforwardly the conservation of the total number quasi-particles $N$. We
now show the conservation of the total energy $E$ in the case where the
potential $W$, entering in the definition of the single-particle energy
$\epsilon$ given by (20), is independent of the rapidity variable $\theta$,
namely $W(z,\theta)=V(z)$. Taking the diffusion operator (26), multiplying it
by $\epsilon$, integrating the result in the phase space $(z,\theta)$ and
using an integration by parts in the position variable $z$, we obtain

$\int_{\mathbb{R}}\mathrm{d}z\int_{\mathbb{R}}\mathrm{d}\theta\,\epsilon\mathcal{D}_{\rho_{\mathrm{p}}}[\rho_{\mathrm{p}}]=\int_{\mathbb{R}}\mathrm{d}z\int_{\mathbb{R}}\mathrm{d}\theta\,(-\partial_{z}V)(1-f\widehat{T})^{-1}\rho_{\mathrm{s}}\widehat{D}\rho_{\mathrm{s}}^{-1}(1-f\widehat{T})^{-1}\rho_{\mathrm{p}}.$
(68)

Using transposition (66), definition (3), and since
$g^{\mathrm{dr}}=(1-\widehat{T}f)^{-1}g$, for any function $g$, we obtain from
(68),

$\displaystyle\int_{\mathbb{R}}\mathrm{d}z\int_{\mathbb{R}}\mathrm{d}\theta\,\epsilon\mathcal{D}_{\rho_{\mathrm{p}}}[\rho_{\mathrm{p}}]$
$\displaystyle=$
$\displaystyle\int_{\mathbb{R}}\mathrm{d}z\int_{\mathbb{R}}\mathrm{d}\theta\,\big{[}(1-\widehat{T}f)^{-1}(-\partial_{z}V)\big{]}\rho_{\mathrm{s}}\widehat{D}\rho_{\mathrm{s}}^{-1}(1-f\widehat{T})^{-1}\rho_{\mathrm{p}}$
(69) $\displaystyle=$
$\displaystyle\int_{\mathbb{R}}\mathrm{d}z\int_{\mathbb{R}}\mathrm{d}\theta\,(-\partial_{z}\epsilon)^{\mathrm{dr}}\rho_{\mathrm{s}}\widehat{D}\rho_{\mathrm{s}}^{-1}(1-f\widehat{T})^{-1}\rho_{\mathrm{p}}$
$\displaystyle=$
$\displaystyle\frac{1}{2\pi}\int_{\mathbb{R}}\mathrm{d}z\int_{\mathbb{R}}\mathrm{d}\theta\,a^{\mathrm{eff}}\rho_{\mathrm{s}}^{2}\widehat{D}\rho_{\mathrm{s}}^{-1}(1-f\widehat{T})^{-1}\rho_{\mathrm{p}}.$

Since the potential $V=V(z)$ depends only on the position variable $z$, we
have $a^{\mathrm{eff}}=-\partial_{z}V$. Using this and setting
$g:=\rho_{\mathrm{s}}^{-1}(1-f\widehat{T})^{-1}\rho_{\mathrm{p}}$, equation
(69) becomes

$\int_{\mathbb{R}}\mathrm{d}z\int_{\mathbb{R}}\mathrm{d}\theta\,\epsilon\mathcal{D}_{\rho_{\mathrm{p}}}[\rho_{\mathrm{p}}]=\frac{1}{2\pi}\int_{\mathbb{R}}\mathrm{d}z\,(-\partial_{z}V)\int_{\mathbb{R}}\mathrm{d}\theta\,\rho_{\mathrm{s}}^{2}\widehat{D}g.$
(70)

Now, inserting the definition (23) for the kernel $D$ into equation (70), we
obtain

$\displaystyle\int_{\mathbb{R}}\mathrm{d}z\int_{\mathbb{R}}\mathrm{d}\theta\,\epsilon\mathcal{D}_{\rho_{\mathrm{p}}}[\rho_{\mathrm{p}}]$
$\displaystyle=$
$\displaystyle\frac{1}{2\pi}\int_{\mathbb{R}}\mathrm{d}z\,(-\partial_{z}V)\int_{\mathbb{R}}\mathrm{d}\theta\,\int_{\mathbb{R}}\mathrm{d}\alpha\,g(\alpha)\Big{[}\delta(\theta-\alpha)\int_{\mathbb{R}}\mathrm{d}\gamma\,W(\gamma,\alpha)-W(\theta,\alpha)\Big{]}$
$\displaystyle=$
$\displaystyle\frac{1}{2\pi}\int_{\mathbb{R}}\mathrm{d}z\,(-\partial_{z}V)\bigg{(}\int_{\mathbb{R}}\mathrm{d}\gamma\int_{\mathbb{R}}\mathrm{d}\alpha\,g(\alpha)W(\gamma,\alpha)-\int_{\mathbb{R}}\mathrm{d}\theta\int_{\mathbb{R}}\mathrm{d}\alpha\,g(\alpha)W(\theta,\alpha)\bigg{)}=0,$

which ends the proof of the conservation of the total energy $E$.

##  Appendix C Construction of Boltzmann collision integral for transverse
state-changing collisions

In this appendix, we outline the construction of the Boltzmann-type collision
integral, originally originally featured in [78]. Given a perturbation to an
integrable Hamiltonian, which only respect a few of the conservation laws of
the unperturbed Hamiltonian, an associated collision integral can be formally
derived following perturbation theory [43]. The conservation laws of the
perturbation (typically number of particles, momentum, and energy) are
explicitly included in the derivation, often determining the rapidities
involved in the scattering process. However, for most experimentally relevant
processes, such formal derivation is infeasible and one must construct the
collision integral following a more phenomenological approach, as well shall
do here.

Assuming the majority of atoms remains in the transverse ground state, we can
restrict our treatment to the three lowest states of the transverse potential,
here denoted by the indices $n=0,1,2$. Within this effective three-level
system, two excitation (and de-excitation) events are possible by virtue of
parity: (1) Two atoms in the ground state collide and both are excited to the
first excited state, or (2) two atoms in the ground state collide and only one
is excited to the second excited state. Similar selection rules exist for de-
exciting collisions.

Consider a state-changing collision of two quasi-particles with rapidities
$\theta$ and $\theta^{\prime}$, and let $\theta_{\pm}$ and
$\theta_{\pm}^{\prime}$ be their rapidities following the collisions. Here,
the subscript ”$-$” denotes an exciting collision, while ”$+$” denotes a de-
exciting one. Following conservation of energy we must have

$\theta^{2}+(\theta^{\prime})^{2}=\theta_{\pm}^{2}+(\theta_{\pm}^{\prime})^{2}\mp
4l_{\perp}^{-2}+\Delta
E(\theta,\theta^{\prime},\theta_{\pm},\theta_{\pm}^{\prime}),$

The term $\mp 4l_{\perp}^{-2}$ denotes the change in transverse potential
energy following the state-changing collision. Meanwhile, the conservation of
momentum implies

$\theta+\theta^{\prime}=\theta_{\pm}+\theta_{\pm}^{\prime}+\Delta
P(\theta,\theta^{\prime},\theta_{\pm},\theta_{\pm}^{\prime}).$

The terms $\Delta E$ and $\Delta P$ describe a shift in the total energy and
momentum, respectively, resulting from a reconfiguration of all local
rapidities following the collision. This effect is captured by the so-called
backflow function, and is a result of the 1D system being intrinsically
strongly correlated [63]. Finding the post-collision rapidities in the
presence of backflow is rather difficult, as the backflow itself is a
functional of the current state of the system (characterized by
$\rho_{\mathrm{p}}$ or $f$). Therefore, to construct a numerically tractable
model, the effect of the backflow was ignored in [78]. Note that for most
experimentally relevant systems this is a reasonable approximation. Neglecting
the backflow, the post-collision rapidities becomes independent of the state
and read
$\theta_{\pm}=\frac{1}{2}(\theta+\theta^{\prime})+\frac{1}{2}(\theta-\theta^{\prime})\sqrt{1\pm
8/[(\theta-\theta^{\prime})l_{\perp}]^{2}}$ and
$\theta_{\pm}^{\prime}=\frac{1}{2}(\theta+\theta^{\prime})-\frac{1}{2}(\theta-\theta^{\prime})\sqrt{1\pm
8/[(\theta-\theta^{\prime})l_{\perp}]^{2}}$. Simplifying the scattering event
to that of just two atoms colliding in a waveguide, the resulting scattering
cross section reads

$P_{\updownarrow}(\theta_{1},\,\theta_{2})=\frac{4c^{2}\theta_{1}\theta_{2}}{\theta_{1}^{2}\theta_{2}^{2}+c^{2}(\theta_{1}+\theta_{2})^{2}}\;.$

Following these considerations, the Boltzmann collision integral can be
constructed. Within GHD, Boltzmann collision integrals are typically
constructed in the conservative form in order to explicitly assure the
conservation of particle number. For the transverse state-changing collisions,
each possible excitation (and de-excitation) event has its own corresponding
collision integral
$\mathcal{I}_{\mathrm{p},\mathrm{h}}^{\pm}[\rho_{\mathrm{p}}]$. Each of these
collision integrals describes the destruction of two quasi-particles (i.e. the
creation of two holes) at the incoming rapidities and creation of two quasi-
particles at the outgoing ones, thus conserving the total number of particles.
The expressions for these collisions integrals can be found in Section 2.2.2.
Although a number of simplifications were made in order to arrive at the
collision integral, it has been demonstrated to successfully describe
experimental observations [78].

## Acknowledgement

We acknowledge support from the Austrian Science Fund (FWF) via the grants SFB
F65 ”Taming Complexity in PDE systems” and, together with the German Research
Foundation (DFG), via the Research Unit FOR 2724 “Thermal machines in the
thermal world.” The work was further supported by the FQXI program on
“Informations as fuel” and the ESQ Discovery Grant “Emergence of physical
laws: from mathematical foundations to applications in many body physics” of
the Austrian Academy of Sciences (ÖAW), and by the Vienna Science and
Technology Fund (WWTF) project MA16-066 ”SEQUEX”.
N.B. acknowledges the hospitality of the Wolfgang Pauli Institute Vienna which
boosted this research.

## References

  * [1] P. H. Acioli.  Review of quantum Monte Carlo methods and their applications.  J. Mol. Struct.: THEOCHEM, 394(2):75–85, 1997.  Proceedings of the Eighth Brazilian Symposium of Theoretical Chemistry.
  * [2] V. Alba, B. Bertini, M. Fagotti, L. Piroli, and P. Ruggiero.  Generalized-hydrodynamic approach to inhomogeneous quenches: correlations, entanglement and quantum effects.  J. Stat. Mech. Theory Exp., 2021(11):114004, nov 2021.
  * [3] A. Bastianello, V. Alba, and J.-S. Caux.  Generalized Hydrodynamics with space-time inhomogeneous interactions.  Phys. Rev. Lett., 123:130602, Sep 2019.
  * [4] A. Bastianello, B. Bertini, B. Doyon, and R. Vasseur.  Introduction to the special issue on emergent hydrodynamics in integrable many-body systems.  J. Stat. Mech. Theory Exp., 2022(1):014001, jan 2022.
  * [5] A. Bastianello, A. De Luca, B. Doyon, and J. De Nardis.  Thermalization of a trapped one-dimensional Bose gas via diffusion.  Phys. Rev. Lett., 125:240604, Dec 2020.
  * [6] A. Bastianello, A. D. Luca, and R. Vasseur.  Hydrodynamics of weak integrability breaking.  J. Stat. Mech. Theory Exp., 2021(11):114003, nov 2021.
  * [7] B. Bertini, M. Collura, J. De Nardis, and M. Fagotti.  Transport in out-of-equilibrium XXZ chains: Exact profiles of charges and currents.  Phys. Rev. Lett., 117(20):207201, 2016.
  * [8] N. Besse.  Convergence of a semi-Lagrangian scheme for the one-dimensional Vlasov–Poisson system.  SIAM J. Numer. Anal., 42:350–382, 2004.
  * [9] N. Besse.  Convergence of a high-order semi-Lagrangian scheme with propagation of gradients for the Vlasov–Poisson system.  SIAM J. Numer. Anal., 46:639–670, 2008.
  * [10] N. Besse, E. Deriaz, and E. Madaule.  Adaptive multiresolution semi-Lagrangian discontinuous Galerkin methods for the Vlasov equations.  J. Comput. Phys., 332, 2017.
  * [11] N. Besse and M. Mehrenberger.  Convergence of classes of high-order semi-Lagrangian schemes for the Vlasov–Poisson system.  Math. Comp., 77:93–123, 2008.
  * [12] N. Besse and E. Sonnendrücker.  Semi-Lagrangian schemes for the Vlasov equation on an unstructured mesh of phase space.  J. Comput. Phys., 191:341–376, 2003.
  * [13] H. Bethe.  Zur Theorie der Metalle.  Zeitschrift für Physik, 71(3):205–226, Mar 1931.
  * [14] G. Beylkin, J. M. Keiser, and L. Vozovoi.  A new class of time discretization schemes for the solution of nonlinear PDEs.  J. Comput. Phys., 147:362–387, 1998.
  * [15] I. Bloch, J. Dalibard, and S. Nascimbène.  Quantum simulations with ultracold quantum gases.  Nature Physics, 8(4):267–276, Apr 2012.
  * [16] I. Bloch, J. Dalibard, and W. Zwerger.  Many-body physics with ultracold gases.  Rev. Mod. Phys., 80:885–964, Jul 2008.
  * [17] M. Borsi, B. Pozsgay, and L. Pristyák.  Current operators in integrable models: a review.  J. Stat. Mech. Theory Exp., 2021(9):094001, sep 2021.
  * [18] I. Bouchoule and J. Dubail.  Generalized Hydrodynamics in the one-dimensional Bose gas: theory and experiments.  J. Stat. Mech. Theory Exp., 2022(1):014003, jan 2022.
  * [19] V. B. Bulchandani, S. Gopalakrishnan, and E. Ilievski.  Superdiffusion in spin chains.  J. Stat. Mech. Theory Exp., 2021(8):084001, aug 2021.
  * [20] V. B. Bulchandani, R. Vasseur, C. Karrasch, and J. E. Moore.  Solvable hydrodynamics of quantum integrable systems.  Phys. Rev. Lett., 119:220604, Nov 2017.
  * [21] V. B. Bulchandani, R. Vasseur, C. Karrasch, and J. E. Moore.  Bethe-Boltzmann hydrodynamics and spin transport in the XXZ chain.  Phys. Rev. B, 97:045407, Jan 2018.
  * [22] B. Buča, K. Klobas, and T. Prosen.  Rule 54: exactly solvable model of nonequilibrium statistical mechanics.  J. Stat. Mech. Theory Exp., 2021(7):074001, jul 2021.
  * [23] P. Calabrese, F. H. L. Essler, and G. Mussardo.  Introduction to ‘quantum integrability in out of equilibrium systems’.  J. Stat. Mech. Theory Exp., 2016(6):064001, jun 2016.
  * [24] O. A. Castro-Alvaredo, B. Doyon, and T. Yoshimura.  Emergent hydrodynamics in integrable quantum systems out of equilibrium.  Phys. Rev. X, 6(4):041065, 2016.
  * [25] F. Cataldini, F. Møller, M. Tajik, J. a. Sabino, S.-C. Ji, I. Mazets, T. Schweigler, B. Rauer, and J. Schmiedmayer.  Emergent pauli blocking in a weakly interacting Bose gas.  Phys. Rev. X, 12:041032, Dec 2022.
  * [26] J.-S. Caux, B. Doyon, J. Dubail, R. Konik, and T. Yoshimura.  Hydrodynamics of the interacting Bose gas in the Quantum Newton Cradle setup.  SciPost Phys., 6:70, 2019.
  * [27] J.-S. Caux and F. H. L. Essler.  Time evolution of local observables after quenching to an integrable model.  Phys. Rev. Lett., 110:257203, Jun 2013.
  * [28] M. A. Cazalilla, R. Citro, T. Giamarchi, E. Orignac, and M. Rigol.  One dimensional bosons: From condensed matter systems to ultracold gases.  Rev. Mod. Phys., 83:1405–1466, Dec 2011.
  * [29] E. Celledoni and B. K. Kometa.  Semi-Lagrangian Runge–Kutta exponential integrators for convection dominated problems.  J. Sci. Comput., 41:139–164, 2009.
  * [30] E. Celledoni, B. K. Kometa, and O. Verdier.  High-order semi-Lagrangian methods for the incompressible Navier–Stokes equations.  J. Sci. Comput., 66:91–115, 2016.
  * [31] B. Cockburn and C.-W. Shu.  Runge–Kutta discontinuous Galerkin methods for convection-dominated problems.  J. Sci. Comput., 16, 2001.
  * [32] S. M. Cox and P. C. Matthews.  Exponential time differencing for stiff systems.  J. Comput. Phys., 176:430–455, 2011.
  * [33] A. C. Cubero, T. Yoshimura, and H. Spohn.  Form factors and Generalized Hydrodynamics for integrable systems.  J. Stat. Mech. Theory Exp., 2021(11):114002, nov 2021.
  * [34] A. J. Daley, C. Kollath, U. Schollwöck, and G. Vidal.  Time-dependent density-matrix renormalization-group using adaptive effective hilbert spaces.  J. Stat. Mech. Theory Exp., 2004(04):P04005, apr 2004.
  * [35] C. De Boor.  A practical guide to splines, volume 27 of Applied mathematical sciences.  Springer, 2001.
  * [36] J. De Nardis, D. Bernard, and B. Doyon.  Hydrodynamic diffusion in integrable systems.  Phys. Rev. Lett., 121:160603, Oct 2018.
  * [37] J. De Nardis, D. Bernard, and B. Doyon.  Diffusion in Generalized Hydrodynamics and quasiparticle scattering.  SciPost Phys., 6:1–72, 2019.
  * [38] G. Dimarco and L. Pareschi.  Numerical methods for kinetic equations.  Acta Numer., 23:369–520, 2014.
  * [39] B. Doyon.  Lecture notes on Generalised Hydrodynamics.  SciPost Phys. Lect. Notes, page 18, 2020.
  * [40] B. Doyon, J. Dubail, R. Konik, and T. Yoshimura.  Large-scale description of interacting one-dimensional Bose gases: Generalized Hydrodynamics supersedes conventional hydrodynamics.  Phys. Rev. Lett., 119:195301, Nov 2017.
  * [41] B. Doyon and T. Yoshimura.  A note on Generalized Hydrodynamics: inhomogeneous fields and other concepts.  SciPost Phys., 2:014, 2017.
  * [42] B. Doyon, T. Yoshimura, and J.-S. Caux.  Soliton gases and Generalized Hydrodynamics.  Phys. Rev. Lett., 120:045301, Jan 2018.
  * [43] J. Durnin, M. J. Bhaseen, and B. Doyon.  Nonequilibrium dynamics and weakly broken integrability.  Phys. Rev. Lett., 127:130601, Sep 2021.
  * [44] D. R. Durran.  Numerical methods for wave equations in geophysical fluid dynamics.  Texts in applied mathematics 32. Springer, New York, 1999.
  * [45] G. A. El.  Soliton gas in integrable dispersive hydrodynamics.  J. Stat. Mech. Theory Exp., 2021(11):114001, nov 2021.
  * [46] F. Filbet and C. Prouveur.  High order time discretization for backward semi-lagrangian methods.  J. Comput. Appl. Math., 303:171–188, 2016.
  * [47] W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal.  Quantum Monte Carlo simulations of solids.  Rev. Mod. Phys., 73:33–83, Jan 2001.
  * [48] I. M. Gamba, J. R. Hauck, C. D. Hauck, and J. Hu.  A fast spectral method for the Boltzmann collision operator with general collision kernels.  SIAM J. Sci. Comput., 39:B658–B676, 2017.
  * [49] M. Gaudin.  Un systeme a une dimension de fermions en interaction.  Physics Letters A, 24(1):55–56, 1967.
  * [50] M. Gaudin and J. Caux.  The Bethe Wavefunction.  Cambridge University Press, 2014.
  * [51] S. Gaudreault and J. A. Pudykiewicz.  An efficient exponential time integration method for the numerical solution of the shallow water equations on the sphere.  J. Comput. Phys., 322:827–848, 2016.
  * [52] S. Giorgini, L. P. Pitaevskii, and S. Stringari.  Theory of ultracold atomic fermi gases.  Rev. Mod. Phys., 80:1215–1274, Oct 2008.
  * [53] A. Görlitz, J. M. Vogels, A. E. Leanhardt, C. Raman, T. L. Gustavson, J. R. Abo-Shaeer, A. P. Chikkatur, S. Gupta, S. Inouye, T. Rosenband, and W. Ketterle.  Realization of Bose-Einstein condensates in lower dimensions.  Phys. Rev. Lett., 87(13):130402, Sep 2001.
  * [54] M. Greiner, I. Bloch, O. Mandel, T. W. Hänsch, and T. Esslinger.  Exploring phase coherence in a 2D lattice of Bose-Einstein condensates.  Phys. Rev. Lett., 87(16):160405, Oct 2001.
  * [55] C. Gross and I. Bloch.  Quantum simulations with ultracold atoms in optical lattices.  Science, 357(6355):995–1001, 2017.
  * [56] P. C. Hammer and J. W. Hollingsworth.  Trapezoidal methods of approximating solutions of differential equations.  Mathematical Tables and Other Aids to Computation, 9(51):92–96, 1955.
  * [57] E. Heath, I. M. Gamba, P. J. Morrison, and C. Michler.  A discontinuous Galerkin method for the Vlasov–Poisson system.  J. Comput. Phys., 231, 2012.
  * [58] M. Hochbruck and A. Ostermann.  Exponential integrators.  Acta Numer., 19:209–286, 2010.
  * [59] M. Hortal.  The development and testing of a new two-time-level semi-Lagrangian scheme (SETTLS) in the ECMWF forecast model.  Q. J. R. Meteorol. Soc., 128:1671–1687, 2002.
  * [60] S. Jaiswal, A. A. Alexeenko, and J. Hu.  A discontinuous Galerkin fast spectral method for the full Boltzmann equation with general collision kernels.  J. Comput. Phys., 378:178–208, 2019.
  * [61] T. Kinoshita, T. Wenger, and D. S. Weiss.  Observation of a one-dimensional Tonks-Girardeau gas.  Science, 305(5687):1125–1128, 2004.
  * [62] T. Kinoshita, T. Wenger, and D. S. Weiss.  A quantum Newton’s cradle.  Nature, 440(7086):900, 2006.
  * [63] V. E. Korepin, N. M. Bogoliubov, and A. G. Izergin.  Quantum Inverse Scattering Method and Correlation Functions.  Cambridge Monographs on Mathematical Physics. Cambridge University Press, 1993.
  * [64] P. Krüger, S. Hofferberth, I. E. Mazets, I. Lesanovsky, and J. Schmiedmayer.  Weakly interacting Bose gas in the one-dimensional limit.  Phys. Rev. Lett., 105(26):265302, Dec 2010.
  * [65] L. D. Landau.  Über Die Bewegung der Elektronen in Kristallgitter.  Phys. Z. Sowjetunion, 3:644–645, 1933.
  * [66] T. Langen, S. Erne, R. Geiger, B. Rauer, T. Schweigler, M. Kuhnert, W. Rohringer, I. E. Mazets, T. Gasenzer, and J. Schmiedmayer.  Experimental observation of a generalized Gibbs ensemble.  Science, 348(6231):207–211, 2015.
  * [67] P. A. Lee, N. Nagaosa, and X.-G. Wen.  Doping a mott insulator: Physics of high-temperature superconductivity.  Rev. Mod. Phys., 78:17–85, Jan 2006.
  * [68] 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.  Relaxation of Bosons in One Dimension and the Onset of Dimensional Crossover.  SciPost Phys., 9:58, 2020.
  * [69] E. Lieb and D. Mattis.  Mathematical Physics in One Dimension: Exactly Soluble Models of Interacting Particles.  Perspectives in physics. Elsevier Science, 2013.
  * [70] E. H. Lieb.  Exact analysis of an interacting Bose gas. II. the excitation spectrum.  Phys. Rev., 130:1616–1624, May 1963.
  * [71] E. H. Lieb and W. Liniger.  Exact analysis of an interacting Bose gas. I. the general solution and the ground state.  Phys. Rev., 130:1605–1616, May 1963.
  * [72] J. Lopez-Piqueres, B. Ware, S. Gopalakrishnan, and R. Vasseur.  Hydrodynamics of nonintegrable systems from a relaxation-time approximation.  Phys. Rev. B, 103:L060302, Feb 2021.
  * [73] V. T. Luan and A. Ostermann.  Parallel exponential Rosenbrock methods.  Comput. Math. Appl., 71:1137–1150, 2016.
  * [74] N. Malvania, Y. Zhang, Y. Le, J. Dubail, M. Rigol, and D. S. Weiss.  Generalized Hydrodynamics in strongly interacting 1d Bose gases.  Science, 373(6559):1129–1133, 2021.
  * [75] N. Mannella, W. L. Yang, X. J. Zhou, H. Zheng, J. F. Mitchell, J. Zaanen, T. P. Devereaux, N. Nagaosa, Z. Hussain, and Z.-X. Shen.  Nodal quasiparticle in pseudogapped colossal magnetoresistive manganites.  Nature, 438(7067):474–478, Nov 2005.
  * [76] MATLAB.  version 9.8.0 (R2020a).  The MathWorks Inc., Natick, Massachusetts, 2020.
  * [77] F. Møller, S. Erne, N. J. Mauser, J. Schmiedmayer, and I. E. Mazets.  Bridging effective field theories and Generalized Hydrodynamics.  arXiv:2205.15871, 2022.
  * [78] F. Møller, C. Li, I. Mazets, H.-P. Stimming, T. Zhou, Z. Zhu, X. Chen, and J. Schmiedmayer.  Extension of the Generalized Hydrodynamics to the dimensional crossover regime.  Phys. Rev. Lett., 126:090602, Mar 2021.
  * [79] F. S. Møller, G. Perfetto, B. Doyon, and J. Schmiedmayer.  Euler-scale dynamical correlations in integrable systems with fluid motion.  SciPost Phys. Core, 3:016, 2020.
  * [80] F. S. Møller and J. Schmiedmayer.  Introducing iFluid: a numerical framework for solving hydrodynamical equations in integrable models.  SciPost Phys., 8:41, 2020.
  * [81] C. Mouhot, L. Pareschi, and T. Rey.  Convolution decomposition and fast summation methods for discrete-velocity approximations of the Boltzmann equation.  ESAIM Math. Model. Numer. Anal., 47:1513–1531, 2013.
  * [82] J. D. Nardis, B. Doyon, M. Medenjak, and M. Panfil.  Correlation functions and transport coefficients in generalised hydrodynamics.  J. Stat. Mech. Theory Exp., 2022(1):014002, jan 2022.
  * [83] R. Orbach.  Linear antiferromagnetic chain with anisotropic coupling.  Phys. Rev., 112:309–316, Oct 1958.
  * [84] L. Pareschi and T. Rey.  Moment preserving Fourier–Galerkin spectral methods and application to the Boltzmann equation.  SIAM J. Num. Anal., 60, 2021.
  * [85] P. S. Peixoto and M. Schreiber.  Semi-Lagrangian exponential integrator with application to the rotating shallow water equations.  SIAM J. Sci. Comput., 41:B903–B928, 2019.
  * [86] L. Pollet.  Recent developments in quantum Monte Carlo simulations with applications for cold gases.  Rep. Progr. Phys., 75(9):094501, aug 2012.
  * [87] J.-M. Qiu and C.-W. Shu.  Positivity preserving semi-Lagrangian discontinuous Galerkin formulation: theoretical analysis and application to the Vlasov–Poisson system.  J. Comput. Phys., 230, 2011.
  * [88] M. Rigol.  Breakdown of thermalization in finite one-dimensional systems.  Phys. Rev. Lett., 103:100403, Sep 2009.
  * [89] M. Rigol, V. Dunjko, and M. Olshanii.  Thermalization and its mechanism for generic isolated quantum systems.  Nature, 452(7189):854–858, apr 2008.
  * [90] M. Rigol, V. Dunjko, V. Yurovsky, and M. Olshanii.  Relaxation in a completely integrable many-body quantum system: An ab initio study of the dynamics of the highly excited states of 1d lattice hard-core bosons.  Phys. Rev. Lett., 98:050405, Feb 2007.
  * [91] J. A. Rossmanith and D. C. Seal.  A positivity-preserving high-order semi-lagrangian discontinuous Galerkin scheme for the Vlasov–Poisson equations.  J. Comput. Phys., 230, 2011.
  * [92] P. Ruggiero, P. Calabrese, B. Doyon, and J. Dubail.  Quantum Generalized Hydrodynamics.  Phys. Rev. Lett., 124:140603, Apr 2020.
  * [93] M. Schemmer, I. Bouchoule, B. Doyon, and J. Dubail.  Generalized Hydrodynamics on an atom chip.  Phys. Rev. Lett., 122(9):090601, 2019.
  * [94] U. Schollwöck.  The density-matrix renormalization group.  Rev. Mod. Phys., 77:259–315, Apr 2005.
  * [95] U. Schollwöck.  The density-matrix renormalization group in the age of matrix product states.  Ann. Physics, 326(1):96–192, 2011.  January 2011 Special Issue.
  * [96] L. L. Schumaker.  Spline functions: basic theory.  Cambridge Mathematical Library. Cambridge University Press, 2007.
  * [97] M. Spiegelman and R. F. Katz.  A semi-Lagrangian Crank–Nicolason algorithm for the numerical solution of advection-diffusion problems.  Geochem. Geophys. Geosyst., 7:1–12, 2006.
  * [98] A. Staniforth and J. Côté.  Semi-Lagrangian integration schemes for atmospheric models–A review.  Month. Weather Rev., 119(9):2206 – 2223, 1991.
  * [99] M. Takahashi.  Thermodynamics of One-Dimensional Solvable Models.  Cambridge University Press, Cambridge, 1999.
  * [100] M. Tokman.  A new class of exponential propagation iterative methods of Runge–Kutta type (EPIRK).  J. Comput. Phys., 230:8762–8778, 2011.
  * [101] R. S. Watson, S. A. Simmons, and K. V. Kheruntsyan.  Benchmarks of Generalized Hydrodynamics for 1d Bose gases.  arXiv:2208.06614v1, 2022.
  * [102] D. Xiu and G. E. Karniadakis.  A semi-Lagrangian high-order method for Navier–Stokes equations.  J. Comput. Phys., 172(2):658–684, 2001.
  * [103] C. N. Yang.  Some exact results for the many-body problem in one dimension with repulsive delta-function interaction.  Phys. Rev. Lett., 19:1312–1315, Dec 1967.
  * [104] C. N. Yang and C. P. Yang.  Thermodynamics of a one-dimensional system of bosons with repulsive delta-function interaction.  J. Math. Phys., 10(7):1115–1122, 1969.