# Simulating generalised fluids via interacting wave packets evolution

Andrew Urilyon, Leonardo Biagetti, Jitendra Kethepalli, Jacopo de Nardis

► To cite this version:

Andrew Urilyon, Leonardo Biagetti, Jitendra Kethepalli, Jacopo de Nardis. Simulating generalised fluids via interacting wave packets evolution. Physical Review B, 2026, 113, pp.014314. 10.1103/b587-8yyt. hal-05475848

HAL Id: hal-05475848

https://hal.science/hal-05475848v1

Submitted on 25 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Simulating generalised fluids via interacting wave packets evolution

Andrew Uriloni ◦,¹ Leonardo Biagetti ◦,¹ Jitendra Kethepalli ◦,¹ and Jacopo De Nardis ◦,¹

¹Laboratoire de Physique Théorique et Modélisation, CNRS UMR 8089,
CY Cergy Paris Université, 95302 Cergy-Pontoise Cedex, France

One-dimensional integrable and quasi-integrable systems display, on macroscopic scales, a uni-
versal form of transport known as *Generalized Hydrodynamics* (GHD). In its standard Euler-scale
formulation, GHD mirrors the equations of a two-dimensional compressible fluid but ignores fluctu-
ations and becomes numerically unwieldy as soon as integrability-breaking perturbations are intro-
duced. We show that GHD can be efficiently simulated as a gas of *semiclassical wave packets*—a
natural generalisation of hard-rod particles—whose trajectories are efficiently mapped onto those
of point particles. This representation (i) provides a transparent route to incorporate integrability-
breaking terms, and (ii) automatically embeds the exact fluctuating-hydrodynamics extension of
GHD. The resulting framework enables fast, large-scale simulations of quasi-integrable systems even
in the presence of complicated integrability-breaking perturbations. It also manifest the pivotal role
of two-point correlations in systems confined by external potentials: we demonstrate that situations
where local one-point observables appear thermalised can nevertheless sustain long-lived, far-from-
equilibrium *long-range* correlations for arbitrarily long times, signaling that, differently from what
previously stated, true thermalisation is *not* reached at diffusive time-scales.

I. Introduction

Generalized hydrodynamics (GHD) [1, 2] has emerged
as a powerful framework for studying large-scale trans-
port and relaxation in one-dimensional (1D) systems
that are either integrable or near-integrable [3–32],
with extensive experimental applications to quasi-one-
dimensional cold atomic and condensed matter setups
[3, 21, 33–41]. While GHD accurately captures the evo-
lution of average densities and currents, it does not fully
describe fluctuations and can become numerically cum-
bersome when including small integrability-breaking per-
turbations. Such perturbations, ubiquitous in realistic
cold-atom or condensed-matter setups, demand complex
analytical and computational corrections to the contin-
uum description.

In this paper, building on Refs. [42–44], and in par-
ticular on [19, 45], we adopt the alternative picture in
which the generalized fluid is viewed as a gas of semi-
classical point-like wave packets whose trajectories map
onto those of noninteracting (“bare”) particles; see Fig. 1.
This particle-based representation recovers the standard
hydrodynamic equations under averaging over initial re-
alisations, but it also goes much beyond by incorpo-
rating correlations and offering a transparent way to
describe integrability-breaking effects. Unlike conven-
tional GHD approaches, which become intractable when
adding noise and higher-order interactions, the soliton-
gas formalism naturally encodes fluctuating hydrody-
namics, since interactions appear as kinetic corrections
in the bare–interacting mapping. Consequently, both av-
erage properties and fluctuations can be treated on an
equal footing.

The bare–interacting soliton-gas mapping thus pro-
vides a versatile framework for exploring a broad class
of phenomena in 1D fluids. It smoothly interpolates
between strictly integrable models and quasi-integrable
regimes dominated by realistic perturbations. We

<figure><img src="image_2.png" /><figcaption>Figure 1. Example trajectories of the WPG gas, with coordinates $x$, as obtained by mapping from the bare coordinates $X$. By choosing an appropriate scattering shift $\mathbf{a}$ the interacting coordinates are realized for different models: pictured are the repulsive hard rod gas $\mathbf{a}(\theta, \theta') = -1$ (upper left), attractive hard rod gas with $\mathbf{a}(\theta, \theta') = 1$ (center), and Lieb-Liniger WPG with $\mathbf{a}(\theta, \theta') = 2c/((\theta - \theta')^2 + c^2)$ and $c = 1$ (upper right), all quenched from a homogeneous initial state to a cosine potential with the interacting dynamics determined from the bare coordinates (bottom). Notice how the attractive hard rod gas has similar scattering dynamics as the Lieb-Liniger gas and therefore can be used as a simpler theoretical toy model.</figcaption></figure>

demonstrate that our approach quantitatively repro-
duces known GHD results in integrable limits and
extends readily to systems with broken integrability
(overcoming the issues coming from solving complicated
partial differential equations [46]). Moreover, we show
that when integrable systems are subjected to external
trapping potentials, the rapidity (the quasiparticle momentum) distribution appears thermalized on diffusive time-scales (as demonstrated in previous works [47, 48]); yet, we also show the surprising and remarkable fact that long-range correlations generated during intermediate times remain finite even at very large times. We argue that intermediate-time ballistic transport builds up these long-range correlations, as first noted in [49–53], (and which are also the source of different anomalous fluctuations in integrable systems, see, for example, [51, 54–57]) and that integrability-breaking terms need additional time-scales in order to fully dissipate them. In cases such as harmonic or quartic traps, these correlations can be strong enough to prevent true thermalization also in the momentum distribution. Our findings thus prompt a reassessment of thermalization in the presence of integrability-breaking. Nevertheless, the relevance of our results is beyond integrable systems, since thermalization in 1D systems is a problem shared by different fields such as plasma physics and astrophysics [58–67].

This paper is organized as follows: Sec. II provides an introduction on Generalized Hydrodynamics for purely integrable systems; in Sec. III we show GHD equations in the presence of spatially inhomogeneous perturbations; in Sec. IV we introduce the wave packet gas formulation of GHD and present the algorithm to implement it in classical and quantum integrability breaking protocols; in Sec. V we show the numerical analysis for the evolution of WPG in inhomogeneous external traps; Sec. VI explains how WPG algorithm can be implemented for more generic two-body integrability breaking potentials; Sec. VII contains our conclusions with an outlook on potential applications of WPG.

## II. Generalised Hydrodynamics and effective viscosities

GHD is known to be the correct hydrodynamic theory in the context of 1D quantum integrable systems, where it describes the non-equilibrium evolution of density fields of conserved charges, provided that the inhomogeneity is sufficiently smooth. There is a particularly natural way to understand GHD by beginning with the continuity equation for the aforementioned local density fields and associated local current

$$ \partial_t q_i(x, t) + \partial_x j_i(x, t) = 0. \quad (1) $$

These local density fields are defined within fluid cells, assumed to have a hydrodynamic length scale $\ell$, that have relaxed such that around their space-time positions, scaled by this hydrodynamic length, these cells are relaxed. Note that $\ell$ sets the length scale, however, by appropriate rescaling of $x, t$ this scale will be exactly canceled and the hydrodynamic equation of the corresponding scale follows. Within each of these fluid cells, a probability distribution is approximately given by a

slowly varying generalized Gibbs ensemble (GGE)

$$ \rho \sim \exp\left[-\ell \int dx \, \beta^i(x, t) q_i(x)\right], \quad (2) $$

where from now on we shall use the repeated-index convention. From this fluid-cell assumption, a hydrodynamic expansion can be performed, which, provided the generalized temperature is sufficiently smooth, determines the current as a functional of the averages of all the conserved charge densities $\langle q_k(x, t) \rangle$,

$$ \langle j_i(x, t) \rangle = j_i[\bar{q}](x, t) + \frac{\delta j_i}{\delta \beta^j} \Big|_{[\bar{q}]} \partial_x \langle q_j(x, t) \rangle + \dots \quad (3) $$

At zeroth order, $O(\ell^0)$, this yields Euler GHD, and at $O(\ell^{-1})$ a Navier–Stokes-type correction associated with diffusive GHD can be found [5, 68]. Importantly, Euler GHD in the presence of external trapping potentials lacks any mechanism for thermalisation (which is expected as the underlying integrability of the system is broken) and requires the presence of the Navier–Stokes terms before any thermalisation can be observed [47, 48, 69]. Namely, it is the combination of internal viscosities and external breaking of integrability (via external forces, for example) that induces thermalization.

Although diffusive GHD successfully captures certain integrability-breaking relaxation dynamics, the Navier–Stokes derivation from the hydrodynamic expansion relies on the assumption that no long-range correlations can develop between fluid cells. This is valid at very short times, within the linear-response regime; however, at intermediate times long-range correlations are known to develop [49, 50]. These transient long-range correlations can be derived using the cumulant expansion, as described in Ref. [53, 70]. There, it was shown that typically integrable systems fail to locally equilibrate, invalidating the applicability of the hydrodynamic states (2). Therefore the correct equilibrium states are those that retain non-trivial information about two-point functions,

$$ \rho_{LR} = Z^{-1} \exp \left[ -\ell \int dx \, \beta_{(1)}^i(x) q_i(x) - \ell \int dx \, dx' \beta_{(2)}^{ij}(x, x') E_{ij}(x, x') \right], \quad (4) $$

where $E_{ij}(x, x') \sim \ell^{-1}$ are the regularised (namely, with the divergent part in $x = x'$ subtracted out) two-point correlations of the hydrodynamic charges; namely, the dynamical correlation functions are composed of a local $\delta$-peaked component, whose amplitude is given by the local susceptibility matrix $C_{ij}$, and a regular long-range part which we denote as $E_{ij}$, namely

$$ S_{ij}(x, y, t) = \delta(x - y) C_{ij}(x, t) + E_{ij}(x, y, t). \quad (5) $$

Using the local-equilibrium states (4) one can expand the current to include $1/\ell$ corrections and derive the hydro-
dynamic equation

$$ \partial_t q_i + A_i^k \partial_x q_k = -\partial_x \left( \frac{1}{2\ell} \frac{\delta^2 j_i}{\delta q_r \delta q_k} \Big|_{\bar{q}} E_{rk}^{\text{sym}}(x) \right) \\ \stackrel{t \ll 1}{\approx} \frac{1}{2\ell} \partial_x (\mathcal{D}_i^k \partial_x q_k), \quad (6) $$

with the point-splitting defined as $E_{ij}^{\text{sym}}(x) = (E_{ij}(x^+, x^-) + E_{ij}(x^-, x^+))/2$. At intermediate times the regular correlation must be evolved simultaneously, according to

$$ \partial_t E_{ij} + \partial_x (A_i^k E_{kj}) + \partial_y (E_{ik} A_j^k) \\ = -\delta(x-y) \langle j_h^- , q_i , q_j \rangle C^{hk} \partial_x q_k. \quad (7) $$

Evolution at the diffusive scale from an initial locally equilibrated GGE state thus requires the long-range correlation functions to be evolved simultaneously, making the dynamics substantially more complicated than those of Euler GHD.

It is often convenient to recast the hydrodynamic description in terms of quasiparticles instead of local charges. A charge density $q_i(x)$ can be written in terms of the quasiparticle density $\rho(\theta, x)$ as

$$ q_i(x) = \int d\theta h_i(\theta) \rho(\theta, x), \quad (8) $$

where the form factors $\{h_i(\theta)\}$ form a complete basis, so that Eq. (8) is invertible. Quasiparticles may be regarded as traveling wave packets (see the following sections) whose group velocity $v_\theta^{\text{bare}}$ is a continuous function of the rapidity $\theta$. Their dynamics involve only elastic two-body scattering events, characterised by the phase shift $\varphi_{\theta,\theta'}$. The bare velocity $v_\theta^{\text{bare}}$, together with the bare momentum $k_\theta$ and the scattering phase, fully specifies the microscopic model at hand.

In the quasiparticle representation, the naive Navier-Stokes equation reads

$$ \partial_t \rho_\theta + \partial_x (v_\theta^{\text{eff}} \rho_\theta) = \frac{1}{2} \partial_x [\mathfrak{D} \partial_x \rho]_\theta, \quad (9) $$

where the effective velocity and the diffusion constant are given in terms of the local quasiparticle density,

$$ v_{\theta}^{\text{eff}} = v_{\theta}^{\text{bare}} + \int d\mu \mathfrak{a}_{\theta\mu} \rho_{\mu} (v_{\mu}^{\text{eff}} - v_{\theta}^{\text{eff}}), \quad (10) $$

$$ \mathfrak{D}_{\alpha\gamma} = R_{\alpha\theta}^{-t} \frac{\delta_{\theta\mu} \left[ \int d\lambda W_{\lambda\theta} \right] - W_{\theta\mu}}{\rho_{T,\theta} \rho_{T,\mu}} R_{\mu\gamma}^{t}, \quad (11) $$

$$ W_{\theta\mu} = \rho_{\theta} f_{\theta} [\varphi_{\theta\mu}^{\text{dr}}]^2 |v_{\theta}^{\text{eff}} - v_{\mu}^{\text{eff}}|, \quad (12) $$

with $R_{\alpha\theta} \equiv \delta_{\alpha\theta} - \varphi_{\alpha\theta} n_\theta$ and the rescaled scattering shift $\mathfrak{a}_{\theta\mu} = 2\pi\varphi_{\theta\mu}/k'_\theta$. Here $\rho_{T,\theta} = k'_\theta/(2\pi) + \int d\mu \varphi_{\theta\mu} \rho_\mu$ is the total density, $g_\theta^{\text{dr}} = g_\theta + \int d\mu \varphi_{\theta\mu} n_\mu g_\mu^{\text{dr}}$ denotes dressing, $n_\theta = \rho_\theta/\rho_{T,\theta}$ is the occupation function, and $f_\theta = \{1, 1 - n_\theta, 1 + n_\theta\}$ is the statistical factor (classical, Fermi-Dirac or Bose-Einstein, respectively). As specified above, the Navier-Stokes equation is only an approximation, valid

at short times, whereas the correct diffusive dynamics are described by the two coupled equations (6) and (7), which, in quasiparticles, read

$$ \partial_t \rho_\theta + \partial_x (v_\theta^{\text{eff}} \rho_\theta) = \frac{1}{2} \partial_x (H_\theta^{\alpha\beta} E_{\alpha\beta}^{\text{sym}}(x)), \quad (13) $$

$$ \begin{aligned} \partial_t E_{\theta\theta'}(x, x') + \partial_x (A_\theta^\alpha E_{\alpha\theta'}) + \partial_{x'} (A_{\theta'}^\alpha E_{\theta\alpha}) \\ = \delta(x - x') M_{\theta\theta'}^\alpha \partial_x \rho_\alpha, \end{aligned} \quad (14) $$

with repeated indices integrated. The tensor operators are

$$ A_{\theta}^{\alpha} = R_{\theta\mu}^{-t} v_{\mu}^{\text{eff}} R_{\mu\alpha}^{t}, \quad (15) $$

$$ M_{\theta\theta'}^{\alpha} = [R_{\theta\mu}^{-t} R_{\theta'\gamma}^{-t} \varphi_{\mu\gamma}^{\text{dr}} \frac{\rho_{\gamma} f_{\gamma}}{\rho_{T,\mu}} (v_{\mu}^{\text{eff}} - v_{\gamma}^{\text{eff}}) R_{\mu\alpha}^{t}]_{(\theta,\theta')}, \quad (16) $$

$$ H_{\theta}^{\alpha\beta} = [R_{\theta\gamma}^{-t} (v_{\mu}^{\text{eff}} - v_{\gamma}^{\text{eff}}) \frac{\varphi_{\gamma\mu}^{\text{dr}}}{\rho_{T,\gamma}} R_{\mu\alpha}^{t} R_{\gamma\beta}^{t}]_{(\alpha,\beta)}, \quad (17) $$

where $[\bullet]_{(\theta,\theta')}$ denotes the sum over permutations of the indices $\theta, \theta'$. The inclusion of two-point correlations is crucial for accurately capturing hydrodynamic evolution up to diffusive time scales. As demonstrated in Ref. [53, 70], the Navier-Stokes equation applies only in the linear-response regime—i.e. for hydrodynamic evolution of initial states weakly perturbed from equilibrium—or at very short times, when the long-range part $E_{\theta\theta'}(x, x')$ can be neglected. However, in generic systems that support ballistic modes—provided the mode velocity depends non-linearly on the local state—hydrodynamic evolution generates long-range correlations that do not decay in time. As argued below, the same remains true in the presence of weak integrability breaking, such as that induced by a trapping potential, making the inclusion of long-range correlations necessary in those cases as well.

### III. Trap-induced integrability breaking and thermalisation

As shown in Refs. [47, 48, 69, 71], starting from a generic integrable Hamiltonian $H_0$ one may introduce a spatially inhomogeneous perturbation of the form

$$ H = H_0 + \int dx V(x) q_{i_0}(x), \quad (18) $$

where $q_{i_0}(x)$ is an arbitrary conserved charge density and $V(x)$ is a smooth external potential.

Under this perturbation quasiparticles not only drift in real space but also acquire a flow in rapidity (momentum) space. Their current can be written as a two-dimensional vector in this space [71],

$$ \mathbf{J} = (v_{\theta}^{\text{eff}} \rho_{\theta}, a_{\theta}^{\text{eff}} \rho_{\theta}), \quad (19) $$

where $a^{\text{eff}}$ is defined analogously as the effective velocity, bare velocity replaced by the bare acceleration $a^{\text{bare}} = -h_{i_0} \partial_x V(x)$. Because transport now proceeds
in both real and rapidity space, the Euler–Navier–Stokes equation generalizes to

$$ \partial_t \rho + \nabla \cdot \mathbf{J} = \frac{1}{2} \nabla \cdot (\mathbf{D} \nabla \rho), \quad \nabla \equiv (\partial_x, \partial_\theta), \quad (20) $$

where $\mathbf{D}$ denotes the rapidity-resolved diffusion tensor.

Equation (20) possesses the thermal local-density-approximation Gibbs ensemble as a fixed point of evolution, namely the thermal state modulated at each point by the external potential,

$$ \rho_{\text{th}}(x) = Z^{-1} \exp[-\beta(H_0 + V(x)q_{i_0}(x) - \mu N)]. \quad (21) $$

Although a formal proof is still missing, it is widely believed to be the unique steady state, apart from exceptional choices of $V(x)$. Accordingly, integrable systems subjected to external potentials are expected to thermalise on diffusive time scales—a prediction borne out numerically except for a few notable deviations in trapped hard-rod gases [72, 73].

However, recent developments have shown that the Navier–Stokes equation (20) is not correct on pre-thermalisation time scales [48]. We therefore have to formulate the problem including two-point correlation functions. Below we derive coupled evolution equations for the density and its connected two-point function. They coincide with those in the previous section, but include the acceleration terms induced by the external force field

$$ \begin{aligned} \partial_t \rho_\theta + \partial_x (v_\theta^{\text{eff}} \rho_\theta) + \partial_\theta (a_\theta^{\text{eff}} \rho_\theta) = & \\ = & \frac{1}{2} \partial_x (H_\theta^{\alpha\beta} E_{\alpha\beta}^{\text{sym}}) + \frac{1}{2} \partial_\theta (H_\theta^{f,\alpha\beta} E_{\alpha\beta}^{\text{sym}}) \end{aligned} \quad (22) $$

$$ \begin{aligned} \partial_t E_{\theta\theta'}(x, x') + \partial_x(A_\theta^\alpha E_{\alpha\theta'}) + \partial_{x'}(A_{\theta'}^\alpha E_{\theta\alpha}) + & \\ + \partial_\theta(A_\theta^{f,\gamma} E_{\gamma\theta'}) + \partial_{\theta'}(A_{\theta'}^{f,\gamma} E_{\theta\gamma}) = & \\ = & \delta(x - x')(M_{\theta\theta'}^\alpha \partial_x \rho_\alpha + M_{\theta\theta'}^{f,\alpha} \partial_\alpha \rho_\alpha). \end{aligned} \quad (23) $$

The matrix definitions are

$$ A_{\theta}^{f, \alpha} = R_{\theta \mu}^{-t} a_{\mu}^{\text{eff}} R_{\mu \alpha}^{t}, \quad (24) $$

$$ M_{\theta\theta'}^{f,\alpha} = [R_{\theta\mu}^{-t} R_{\theta'\gamma}^{t} \varphi_{\mu\gamma}^{\text{dr}} \frac{\rho_{\gamma} f_{\gamma}}{\rho_{T,\mu}} (a_{\mu}^{\text{eff}} - a_{\gamma}^{\text{eff}}) R_{\mu\alpha}^{t}]_{(\theta,\theta')}, \quad (25) $$

$$ H_{\theta}^{f, \alpha\beta} = [R_{\theta\gamma}^{-t} (a_{\mu}^{\text{eff}} - a_{\gamma}^{\text{eff}}) \frac{\varphi_{\gamma\mu}^{\text{dr}}}{\rho_{T,\gamma}} R_{\mu\alpha}^{t} R_{\gamma\beta}^{t}]_{(\alpha,\beta)}. \quad (26) $$

In Appendix C we explicitly derive Eq. (22) and (23) and show that, for uncorrelated initial states $E(t = 0) = 0$, in the small time limit $t \to 0^+$ they reproduce Eq. (20).

Two remarks are necessary. First, a faithful simulation of the system requires evolving the fully inhomogeneous correlation function, which is numerically challenging. Second, one should ask about the fate of thermalisation in this system of coupled equations. We cannot provide a definitive answer here and leave it as an open question. However, since there is no term that dissipates correlations, the system cannot thermalise at order $O(\ell^2)$, contrary to what the Navier–Stokes equation would suggest. In fact, full thermalisation of the two-point functions appears only at dispersive order in the hydrodynamic expansion, i.e. on time scales $\sim O(\ell^3)$. More precisely, the

relaxation of the two-point function involves three-point functions [52] and can be captured by a diffusive contribution proportional to

$$ \partial_t E_{\theta\theta'} \propto H_{\theta}^{\alpha\beta} E_{\alpha\beta\theta'}^{\text{sym}}(x^+, x^-, x') + \dots, \quad (27) $$

where $E_{\alpha\beta\theta'}^{\text{sym}}$ is the three-point function, symmetrised with respect to the first two spatial variables. In this hierarchical structure, the three-point function $E_{\alpha\beta\theta'}^{\text{sym}}$ is expected to evolve through an equation coupled to the first two layers of the hierarchy. The resulting cascade of correlations is clearly too complex to numerically simulate, the following section provides a numerical method that circumvents the problem by simulating the effective wave-packet gas, which provides the correct discretisation and contains by definition the correct correlations’ dynamics.

## IV. Wave packet gas description of GHD

As we have seen in the previous section, a complete solution of the Generalized hydrodynamics (GHD) including diffusive effects that is capable of describing a system’s thermalisation would require solving two coupled equations, one of which describes two-point functions. From a numerical standpoint, the two-point function is quite a complex object. Therefore, especially to go beyond the realm of strictly integrable systems and to be able to efficiently incorporate integrability-breaking perturbations, it is desirable to introduce a particle model whose averages over realizations coincide with those of integrable models in the hydrodynamic limit.

Following recent works, [19, 42–45] we introduce the following gas of wave packets (WPG) for this purpose. First, we observe that for a simple model such as the repulsive hard rods, the wave packets coincide with the rods themselves. For a generic interacting integrable model, the WPG is simply described by a gas of interacting solitonic-like particles that scatter via the model-dependent scattering shift. Although another effective particle gas was already introduced in [74], known as the “flea gas”, the latter missed the correct hydrodynamic fluctuations. Our method instead generalizes similar approaches already employed for evolving the repulsive hard rod gas to more general settings and models, and, as we shall show, guarantees the correct viscosities, thermalisation times under integrability breaking, and correlations.

### IV.1. The classical hard rods gas

To elucidate our method, we consider a classical hard rod gas, consisting of either attractive or repulsive particles. In the repulsive case this is the most familiar form of a hard rod gas [75–80], whose extended particles have a finite length $a > 0$ that are unable to overlap and instead
bounce off of one another exchanging momenta whenever their edges come into contact. In the attractive case the particles all have regions around their center such that all neighbors within a distance $|a|$ away are pulled into the center, as a consequence when a collision occurs these attractive particles will stick together for a period of time based on the relative speed of the colliding particles. The position of the i-th particle, $x_i$, is transformed to a corresponding bare particle position, $X_i$, by use of the transformation

$$X_i = x_i \mp \frac{a}{2} \sum_k \mathrm{sign}(x_i - x_k). \quad (28)$$

with the sign corresponding to repulsive (-) and attractive (+) case, respectively, and where the bare particles coordinates are given by free ballistic propagation

$$X_i(t) = X_i(0) + \theta_i t. \quad (29)$$

For the repulsive case this transformation is easily inverted, since the bare particle positions {$X_i$} and WPG {$x_i$} coordinates are identically ordered, which implies $\mathrm{sign}(x_i - x_k) = \mathrm{sign}(X_i - X_k)$. Hence we can write

$$x_i = X_i + \frac{a}{2} \sum_k \mathrm{sign}(X_i - X_k) \quad (\text{repulsive HR}). \quad (30)$$

For attractive hard rods, a naive inversion will result in cases with bare coordinates $X_i < X_k$ and interacting coordinates $x_i > x_k$, hence differently ordered, whenever two particles are too close together. Thus, additional care is needed here to keep the mapping consistent. As a resolution, whenever wave packets sit too close to one another they merge into a single position in an example of clustering. This behavior is also observed for more general WPGs, where the trajectories stick together before separating at a later time.

While the dynamics given by (28) is deterministic, one then needs to sample over an ensemble of initial states in order to describe a thermodynamic state. However the advantage of working with (28) is crucial: one does not need to know or use the invariant or thermal distribution of interacting rods, but one can simply initialise the bare particles variables $X_i, \theta_i$, whose statistics is known to be simple Poissonian, and then perform the mapping to interacting coordinates. We shall show below the algorithm to efficiently sample generic classical hard rods particles.

To make use of the transformation for attractive hard rods the transformation must be inverted, a gradient transform can be implemented; however, for the attractive hard rods a series of inequalities can also be used to determine the interacting positions from the bare coordinates. The key to this method is the introduction of a naive transformation, where $\mathrm{sign}(x_i - x_j) \to \mathrm{sign}(X_i - X_j)$, recalled to be the correct inversion for repulsive hard rods, although generally incorrect in the attractive case. In addition to this naive transformation, the ordering of the bare positions $X_j < X_k$ implies $x_j \le x_k$ when $x_j = x_k$ these particles are said to be within

a cluster. With this in mind it is possible to identify the positions of clusters by finding where this weak ordering fails after a naive transformation.

Concretely, the bare coordinates $X_{I_1} < \dots < X_{I_m}$ belonging to a single cluster are ordered as $I_1 < \dots < I_m$. Surprisingly, but usefully, the average position of clustered particles in both the true and naive transformations are identical. As a consequence the naive transformation implies the following condition on the first, $x_{I_1}^{\text{naive}}$, and last particle positions, $x_{I_m}^{\text{naive}}$, within a cluster of m-particles relative to the true cluster position $x_I$

$$x_{I_m}^{\text{naive}} < x_I < x_{I_1}^{\text{naive}}. \quad (31)$$

This condition can be used to group the particles into clusters by hand, which effectively inverts the transformation for the attractive hard rod gas. With this approach numerical solutions could be reproduced more quickly and stably for the attractive hard rod gas.

## IV.2. Generic integrable models and quantum statistics

Generic integrable systems can be viewed as deformations of the hard-rod gas. In particular, they are characterised by effective quasi particles whose large-scale generalised hydrodynamics coincide with that of hard rods under the substitution

$$a \to a_{ij} = \frac{2\pi \varphi(\theta_i - \theta_j)}{k'(\theta_i)}. \quad (32)$$

Some relevant examples, especially its experimental applications, are provided by the Lieb-Liniger model, which describe interacting (quantum) bosonic atoms [81–84] with Hamiltonian given by

$$H_{\text{Lieb-Liniger}} = \sum_j \frac{p_j^2}{2m} + 2c \sum_{i<j} \delta(x_i - x_j), \quad (33)$$

and associated scattering shift $a_{ij} = 2c/((\theta_i - \theta_j)^2 + c^2)$, or the celebrated Toda model [85, 86], for its numerous studies in classical integrability-to-chaos transitions [87, 88], with Hamiltonian

$$H_{\text{Toda}} = \sum_j \frac{p_j^2}{2m} + \sum_{j=1}^{N-1} e^{-(x_{j+1}-x_j)}, \quad (34)$$

and associated scattering shift $a_{ij} = -2\log(\theta_i - \theta_j)$. We also set the mass as $m = 1$.

The hard-rod mapping of Eq. (28) then generalises to generic models, relating interacting and bare particle coordinates through

$$X_i = x_i + \frac{1}{2} \sum_{k=1}^{N} a_{ik} \mathrm{sign}(x_i - x_k), \quad \dot{X}_i = v^{\text{bare}}(\theta_i), \quad (35)$$
where the bare velocity parameterization $v_{\theta}^{\text{bare}}$ can account for different types of models (in generic Galilean invariant models one typically has $v_{\theta}^{\text{bare}} = \theta$). This equation was first introduced in [19, 45]. We shall show in Appendix A the argument showing that the dynamics induced by (35) leads indeed to GHD type of hydrodynamics (notice that a similar argument is also provided in [44, 45]).

A crucial difference compared to the hard rods mapping is that it cannot be inverted explicitly. Therefore, in order to invert the transformation a convex surface defined by the phase space position of each particle, $(\vec{X}, \vec{\theta})$, is introduced such that the interacting positions are those positions that minimize the surface

$$ S_{\vec{X},\vec{\theta}}(\vec{x}) = \sum_{i} \frac{(x_i - X_i)^2}{2} + \frac{1}{2} \sum_{i,j} a_{ij} |x_i - x_j|, \\ \vec{x}(t) = \arg \min S_{\vec{X}(t),\vec{\theta}}(\vec{x}). \qquad (36) $$

Numerically, a gradient descent method is capable of carrying out this minimization, thus identifying a set of interacting coordinates from a set of bare coordinates. By inverting this formula an initial fluctuating state of bare particles is mapped to an initial fluctuating state of a general quantum integrable system. Averages over these initial states reproduce the correct expectation values and correlation functions.

A less complete particle model of GHD has been realized in the literature prior to both this work and Ref. [45]. In Ref. [74] the authors introduced the so-called flea gas algorithm, where hard rod-like particles scatter when they find themselves at distance $a_{ij}$. However, such an algorithm does not properly treat multi-scattering processes and, even if it gives the correct average effective velocities [74, 89, 90], it fails to give the correct diffusive terms as both initial fluctuations and correct dynamics are not properly described.

The WPG dynamics of Eq. (35) obviously reproduces the correct effective velocity, defined as $t^{-1}\langle x_i \rangle = v^{\text{eff}}(\theta_i)$ at large times, and since the WPG contains accurate initial fluctuations it also reproduces additional hydrodynamical effects such as the diffusive spreading of quasiparticles, namely the quantity $\delta x_i = t^{-1}\langle (x_i - \langle x_i \rangle)^2 \rangle$ at large times, see Fig. 2. The latter is known to be related to the diffusion constant given by [6, 68]

$$ \delta x_{\theta} = \int d\alpha \rho_{\alpha} f_{\alpha} |v_{\theta}^{\text{eff}} - v_{\alpha}^{\text{eff}}| \left( \frac{\varphi_{\theta,\alpha}^{\text{dr}}}{\rho_{T,\theta}} \right)^2, \qquad (37) $$

which clearly depends, through the factor $f_{\alpha}$, on the statistics of the initial fluctuations, which can be chosen to be either quantum or classical. In Appendix B, we sketch the particle sampling procedure used for both classical and Fermi-Dirac statistics and demonstrate that the correct initial fluctuations are present.

<figure><img src="image_3.png" /><figcaption>Figure 2. Evolution of a homogeneous WPG sampled with Fermi-Dirac statistics using the scattering kernel of the Lieb-Liniger model with $\beta = c = \mu = 1$ and with $dX = 40$. A tracing particle is placed in the center of the system with $\theta = -2, -1, 0, 1, 2$ and the displacement $\Delta x = x(t) - x(0)$ of this particle from its initial position is considered (left) the average displacement $\langle x \rangle$ approaches the corresponding effective velocity from GHD (right) the second moment of displacement approaches the diffusive constant predicted from Eq. (37).</figcaption></figure>

## V. Generalised fluids (I): integrability breaking via external potentials

Now that we have described how to properly initialize an initial state, we can consider the dynamics, especially in the presence of integrability breaking terms. When the dynamics is given entirely by integrable Hamiltonian, the corresponding bare particles positions $X_i(t)$ evolve freely. However, in the presence of an external potential as given by Eq. (18), these trajectories are now modified. The logic is simple: consider the interacting Hamiltonian and map it to bare particles using the bare-to-interacting mapping of Eq. (35). While this mapping makes the interacting Hamiltonian $H_0$ trivial, it gives a new effective potential for the bare particles

$$ H_0 + \int dx V(x) q_{i0}(x) \rightarrow \sum_i \left[ \frac{(v^{\text{bare}}(\theta_i))^2}{2} + V(X_i(x)) \right]. \qquad (38) $$

As the potential is computed at position $X_i(x)$ we are required to know the position of the interacting particles $x_i$. The result is a system of interacting particles with an interaction mediated by the external potential.

As an example, let us consider the case of a repulsive hard rods gas in the presence of a quadratic external potential $V(x) = \omega^2 x^2/2$. In this case, the Hamiltonian evolution can be directly mapped to an interacting Hamiltonian for the bare particles as

$$ H_{\text{HR(a)}}[X, \theta] = \sum_{i} \left[ \frac{\theta_{i}^{2}}{2} + \frac{\omega^{2} X_{i}^{2}}{2} + \frac{a\omega^{2}}{2} \sum_{j \neq i} |X_{i} - X_{j}| \right]. \qquad (39) $$

The latter corresponds to the Hamiltonian of a self-gravitating gas [63, 64, 91–95] confined in a harmonic
trap, whose equation of motions are written as

$$ \dot{X}_i = \theta_i, \quad \dot{\theta}_i = -\omega^2 X_i - \frac{a\omega^2}{2} \sum_{i \neq j} \operatorname{sgn}(X_i - X_j). \quad (40) $$

In an analogous way, in generic interacting models the equations of motion can be written via the mapping (35) but with the bare particles now satisfying the following equations of motions, whose trajectories are shown in Fig. 1:

$$ \dot{X}_i = v^{\text{bare}}(\theta_i), \quad \dot{\theta}_i = -(\partial_x V(x))\Big|_{x=x_i}. \quad (41) $$

In particular, Eq. (41) shows the main complication of this approach: bare particles have a trivial kinetic spreading but they feel a complicated, state dependent, non-local force that is a functional of the interacting particles coordinates. Hence, it requires transforming between the bare and interacting coordinate systems for each time step and represents the main computational bottleneck of the numerical algorithm. This phenomenon resembles the intrinsic many-body complexity of interacting integrable systems.

It is important now to stress that, although the single realization of microscopic time evolution is fully deterministic, the average over the initial state ensemble introduces the statistical complexity through the initial state fluctuations, that are transported through the dynamics according to Eq. (23). Hence, this method is not only capable of reproducing the Euler evolution of the gas, but in principle any higher order correction to the hydrodynamics. On the other hand, since the n-th hydrodynamic corrections produces effects suppressed as $1/\ell^{n-1}$, the expected number of initial states over which perform the average is expected to grow as $\sim \ell^{2(n-1)}$ in order to be able to observe the effect of n-th term in the hydrodynamic expansion. Thus, to show that this method reproduces the correct diffusive order hydrodynamics, we focus on two toy models of WPG: attractive and repulsive hard rods, for which the mapping can be more easily inverted. More precisely, in Fig. 3, we compare diffusive hydrodynamics obtained from WPG simulations with Navier-Stokes (20) for the evolution in the periodic potential defined in Eq. (42). As expected [48], for both attractive and repulsive HR the WPG simulations match almost perfectly with the diffusive equation at small times $t/\ell < 1$, i.e. when Eq. (20) is expected to be predictive. This must be considered as a non trivial benchmark that confirms the validity of the numerical method.

## V.1. Simulating trapped integrable systems: thermalisation and long-range correlations

We shall now use the methods introduced in the previous section to study the dynamics of integrable systems in

<figure><img src="image_4.png" /><figcaption>Figure 3. The two columns show respectively the dynamics of classical WPG with $a_{ij} = -1$ (left, repulsive) and $a_{ij} = +1$ (right, attractive). The initial state is a Bragg pulse, given by $\rho = \sum_{s=\pm} \exp\{-(\theta - \theta_s)^2/(2T_0)\}$ with $\theta_{\pm} = \pm 1$. The system is quenched into a cosine potential $V(x) = V_0 \cos(2\pi x/\ell)$ with $\ell = 10$ and $V_0 = 1.0$. The top plots show the dynamics of kurtosis of the rapidity distribution as function of rescaled time $t/\ell$ for $T_0 = 1.0$ and $T_0 = 0.1$, where solid lines represent the WPG numerics and dashed lines the solution of Navier-Stokes GHD (20). The insets show the dynamics of rescaled rapidity distribution for $T_0 = 0.01$, where different colors represent different times, according to the colorbars shown in the plots below. The solid black line shows a featureless Gaussian distribution, while the dashed line represents the initial distribution. The second and third lines show respectively the dynamics of density and energy correlations, defined in (44). Different colors represent different increasing times from light to dark, according to the colorbars.</figcaption></figure>

the presence of external traps. We consider here, for convenience, the case of attractive and repulsive hard rods, with scattering shift $a_{ij} = \pm 1$ and with classical statics. The attractive hard-rod model can be considered as a toy model for the Lieb-Liniger gas at high temperature and therefore its dynamics is directly relevant for experimental realisation. We start by considering particles in a cosine potential

$$ V(x) = V_0 \cos(2\pi x/\ell). \quad (42) $$

This potential was also considered in previous works [47, 48], where it was concluded that the system thermalises on, diffusive, time scales proportional to $(\ell/V_0)^2$. Moreover, it was remarked in Ref. [48] that diffusive hydrodynamics with Navier-Stokes (20) can capture only the
<figure><img src="image_5.png" /><figcaption>Figure 4. The four columns show respectively the dynamics of repulsive hard rods evolved in harmonic potential (I), repulsive hard rods evolved in quartic potential (II), attractive hard rods evolved in harmonic potential (III) and attractive hard rods evolved in quartic potential (IV). In all the case we set hard rods length $|a| = 1$ and we use $\omega_0 = \sqrt{50}$ for the repulsive case and for the attractive quartic case, while $\omega_0 = \sqrt{5}$ for the attractive harmonic. The initial state is homogeneous in the interval $[-\ell/2, \ell/2]$ with an average number of particles $N = \ell/2 = 50$ and is prepared in a Bragg pulse, given by rapidity distribution $\rho = \sum_{s=\pm} \exp\{-(\theta - \theta_s)^2/(2T_0)\}$, with $\theta = \pm 1$ and $T_0 = 0.1$. The first line shows the WPG dynamics of kurtosis of rapidity distribution against rescaled time $t/\ell$. The insets show the dynamics of rescaled rapidity distribution, where different colors represent different time data according to the colorbars shown in the plots below. The solid black line shows a featureless Gaussian distribution, while the dashed line represent the initial distribution. The second and third lines respectively show the dynamics of space and energy correlations defined in (44), computed at $y = 0$. Different colors represent different increasing times from light to dark, according to the colorbars. Dashed and solid lines are respectively the initial and final time data.</figcaption></figure>

initial Euler-like dynamics and the final thermalisation dynamics, while there is an intermediate time window in which the dynamics are strongly discontinuous (such a region was denoted a turbulent phase) and where hydrodynamics and exact simulation do not fully agree. Such an effect is enhanced by larger values of $\ell$ which prolong the turbulent phase in time.

To benchmark the WPG method a cosine potential is simulated with both the WPG and Navier-Stokes GHD. In this cosine potential case a full simulation with Navier-Stokes GHD is possible, as the state remains smooth during time evolution and no strong singularities in the distribution $\rho_\theta$ are created. We can see that the prediction for the point function agree during the dynamics. In both the repulsive and attractive regimes the distribution of rapidities $\rho(\theta, x)$ converges to the expected thermal distribution (which is also the fixed point of the Navier-Stokes GHD equation Eq. (20), see [48]), namely

a Gaussian (thermal) function,

$$ \rho(\theta, x) \underset{t \sim \ell^2}{\rightarrow} e^{-\beta(\theta^2 + V(x) - \mu)}, \quad (43) $$

up to very small deviations. The thermalization predicted by Navier-Stokes GHD can be further interrogated with the WPG method, where higher-point correlations are accessible. For a thermal state these higher-point correlations are expected to become local as the thermal state is approached at the diffusive scale, however an analysis of the two-point correlations suggests a different picture. For simplicity, we focus on two different types of connected correlations, the density-density $\langle N(x)N(y) \rangle^c$ and energy-energy $\langle e(x)e(y) \rangle^c$, namely

$$ E_{0,0} = \int d\theta d\theta' E_{\theta\theta'}(x, y) = \langle N(x)N(y) \rangle_{\text{reg}}^c, \quad (44) $$

$$ E_{2,2} = \int d\theta d\theta' \theta^2 \theta'^2 E_{\theta\theta'}(x, y) = \langle e(x)e(y) \rangle_{\text{reg}}^c, $$
where for regularised correlations $\langle \cdot \rangle_{\text{reg}}$ we mean that they can be obtained by binning the microscopic system with some mesoscopic space interval, (in this paper we take quite a large one to improve statistics, $dx = 3\ell \sim 30$), and we exclude the singular contribution at $x = y$. We observe that as the system evolves these correlations develop a finite long-range part (namely, particles correlate over multiple units of $\ell$) and, even when the one-point function appears thermal, such correlations remain finite, signaling that the state is not truly thermal but is instead a long-range state of the type (4). Unless stated otherwise, true thermalization would imply a correlation length of the order of the thermal correlation length, which is $O(1)$, whereas here we observe correlation lengths of order $O(10\ell)$, which therefore signals that the system has not reached a proper thermal equilibrium even at times of order $\ell^2$. As specified in the previous section, one should therefore expect an even larger time scale of order $\ell^3$ in order, at least in principle, to see two-point functions vanish. Notice that such an effect is clearly due to inter-particle interactions, as in the limit of vanishing interaction one would have only independent trajectories inside the potential.

We then move to the case of harmonic and quartic traps, namely with $V(x) = \omega^2 x^2/2$ and $V(x) = \omega^4 x^4/4$. To ensure extensive energy, the trap frequency is scaled as $\omega = \omega_0/\ell$ and $\ell = N$, which corresponds to weak-integrability breaking. We should first stress that a full simulation with Navier-Stokes GHD in these cases is not possible, as the dynamics become increasingly "turbulent" as function of time, generating strong gradients that are difficult to resolve with finite difference methods. We first study the repulsive ($a = -1$) and attractive hard rod ($a = +1$) gases prepared in an initial Bragg state, as shown in Figs. 4, and released into the respective traps. For the harmonic trap, the one-point function struggles to thermalise, a phenomenon that was well discussed in previous works [72, 73]. Here we show that, while the distribution of rapidities fails to approach a Gaussian, even the two-point correlations remain finite and large at times of order $\ell^2$. Meanwhile, in the quartic case, there are some indications that at times larger than $\ell^2$ the one-point function thermalise. This seems consistent with the results of Ref. [73] for the quartic potential, where due to strong integrability breaking ($\omega \sim \omega_0$) the one-point function approaches a Gaussian at $t \sim O(N)$. Yet, even in this case correlations functions do not give signs of decaying at these diffusive times, signaling that true thermalisation is not really reached.

To amplify and elucidate the role of these correlations, we implement a second protocol (Fig. 5): starting with a Gaussian rapidity distribution and uniform spatial density, we remove the particles from the center of the system, creating a Newton’s cradle setup, but which can be more efficient realizable experimentally. Strikingly, both harmonic and quartic traps now exhibit a non-Gaussian

rapidity distribution and long-range two-point correlations that fail to vanish even at time-scales that are one order of magnitude more than diffusive ones.

Comparisons of the pure trap to the chop protocol make it clear that the two-point correlations are clearly modified between the two settings and that they are strong enough to be detectable in current state-of-the-art experimental platforms. Particularly notable is the very large energy-energy correlator, which clearly supports the idea that there are preferred initial states for investigating these long-range correlators. Therefore, testing different initial states could be used to identify promising experimental protocols that could in turn be realized in cold atom experiments to investigate the emergence of these long-range correlations.

## VI. Generalised fluids (II): integrability breaking via two-body potential

The approach introduced in the previous sections allows us to tackle more complex integrability-breaking potentials. In this section, we consider integrable systems perturbed by two-body interactions, a scenario motivated by recent cold-atom experiments. For example, in [96], a gas of ultracold bosonic dysprosium atoms is tightly confined in one-dimensional tubes formed by a two-dimensional optical lattice. These atoms exhibit dipole-dipole interactions (DDI) roughly 100 times stronger than those of rubidium. Whereas ultracold bosons usually interact locally and are well described by the integrable Lieb-Liniger model, the DDI introduces a long-range, two-body, integrability-breaking interaction that drives late-time thermalization of the gas, namely the system is described by the Hamiltonian

$$H = H_{\text{Lieb-Liniger}} + \sum_{i} V(x_i) + \sum_{i<j} U(x_i - x_j). \quad (45)$$

As shown in [52], the hydrodynamic evolution of such a quasi-integrable gas is governed by a hierarchy of nonlinear integro-differential equations, analogous to the BBGKY hierarchy. Under the assumptions of weak interactions and a large hydrodynamic scale $\ell$, this hierarchy can be closed, but its numerical solution remains computationally prohibitive. We now illustrate how the numerical method developed here can efficiently simulate the hydrodynamics of these systems.

Consider the long-range potential

$$U(x - x') = \frac{U_0}{\ell} \left[ 1 + \left( \frac{|x - x'|}{\ell} \right)^3 \right]^{-1}, \quad (46)$$

written in interacting coordinates and evolved via the bare-particle mapping. Our wave-packet gas (WPG) algorithm applies directly, requiring only the coordinate transformations between bare and interacting frames. In the presence of an external trap $V(x)$, the equations of
<figure><img src="image_6.png" /><figcaption>Figure 5. Same as in Fig. 4 with the only difference in the initial state preparation: we initialize the system in a thermal state in the interval [-ℓ, ℓ] with temperature T₀ = 0.01 and N = ℓ = 100, i.e. homogeneous in space and with Gaussian rapidity distribution. Then, at time t = 0, the particles in [-ℓ/2, ℓ/2] are removed and we evolve the two chopped clouds.</figcaption></figure>

motion become

$$ \dot{X}_i = v^{\text{bare}}(\theta_i), \quad (47) $$

$$ \dot{\theta}_i = -\partial_x V(x)|_{x=x_i} - \sum_{k \neq i} \partial_x U(x)|_{x=x_i-x_k}. \quad (48) $$

In bare coordinates, the two-body integrability-breaking potential acquires dependence on the interacting-particle positions, yielding a 1/ℓ expansion of many-body interactions in which the n-body term scales as 1/ℓⁿ⁻².

As an illustrative example, Fig. 6 shows the dynamics of an attractive WPG with $a_{ij} = c$, interacting via the potential (46). This toy model approximates the high-temperature limit of a DDI-interacting Lieb–Liniger gas. Starting from a Bragg-pulse initial state, we track the density and connected correlations. The kurtosis of the rapidity distribution relaxes rapidly on a time scale set by the integrability-breaking strength. After this fast relaxation, the density reaches a quasi-steady plateau and then drifts slowly toward the thermal distribution. Simultaneously, intermediate-time correlations build up and relax to a correlation length of order ℓ. Since the true thermal state exhibits correlations of the same range, Fig. 6 confirms that the system indeed fully thermalizes in this case.

## VII. Conclusion

In this work, we have shown that the semiclassical wave-packet formulation of GHD [19, 42–45] provides a transparent and efficient algorithm for numerically evolving both classical and quantum integrable systems, including all orders of fluctuating hydrodynamics. Importantly, our method can incorporate arbitrary integrability-breaking terms. On one hand, it overcomes the challenge of solving the full hydrodynamic hierarchy for a quasi-integrable model, yielding all-order hydrodynamics efficiently; on the other, it embeds generic perturbations—many-body interactions, particle loss, or impurity scattering—simply by adding the corresponding forces to the bare-particle dynamics.

Motivated by recent cold-atom experiments, we applied the WPG algorithm to two paradigmatic scenarios: an external periodic potential (Eq. (42)) and a long-range dipolar two-body interaction. In the periodic trap, at times $O((ℓ/V_0)^2)$ the density relaxes to a nearly Gaussian form, consistent with previous findings [48], however, correlations reveal persistent far from thermal long-range terms that are built-up at intermediate time-scales and that fail to dissipate at large times. According to Eq. (27), two-point functions are expected to relax on the much longer scale $ℓ^3$. On the other hand, for polyno-
<figure><img src="image_7.png" /><figcaption>Figure 6. Dynamics of classical attractive WPG with $a_{ij} = 1$ and long-range potential (46) with $\ell = 20$. The initial state is a Bragg split state with $\theta = \pm 1$ and $T_0 = 0.01$. The top left panel shows an example of microscopic trajectories for both interacting and bare coordinates. Bottom left figure shows the dynamics of kurtosis of rapidity distribution. The inset shows the rapidity distribution at different times, compared with the expected Gaussian thermal shape (black dashed line). The two figures on the right show the dynamics of density and energy correlations defined in (44). Different colors are associated with different increasing times from light to dark, as explained by the colorbar. The dashed and solid lines are respectively the initial and final time data.</figcaption></figure>

mial traps, in agreement with Refs. [72, 73], the rapidity distribution remains non-Gaussian and long-range correlations persist for even longer times, indicating that true thermalization is still yet to be reached. The final reason for such a failure of dissipation of correlations is left for future works, even if the bare-to-interacting mapping shown here seems to point to the fact that the integrable anyway the weakly mixing nature or trajectories of the associated bare particles may play an important role. Moreover, varying initial states strongly affect the generation of long-range correlations, suggesting that tailored initial conditions could enable experimental protocols to probe correlation dynamics in cold-atom setups.

Our findings raise new questions about the role of long-range many-body correlations in the relaxation of quasi-integrable systems, particularly regarding the absence of thermalization in harmonic confinement. They also open exciting directions for studying nearly integrable systems with diverse integrability-breaking mechanisms, potentially modeling realistic features of cold-atom experiments.

## Acknowledgments

A.U., L.B, J.K and J.D.N. are funded by the ERC Starting Grant 101042293 (HEPIQ) and the ANR-22-CPJ1-0021-01. This work was granted access to

the HPC resources of IDRIS under the allocation AD010613967R2. We thank Alvise Bastianello, Frederik Møller, Jorg Schmiedmayer, Benjamin Doyon, Romain Vasseur, Sarang Gopalakrishnan, Takato Yoshimura, Friedrich Hübner and Herbert Spohn for valuable feedbacks.

## Data availability

The data supporting the findings of this study are not publicly available due to the costs associated with its hosting and preparation. The data are available from the corresponding author upon reasonable request.

## A. Dynamics induced by bare-to-interacting mapping and GHD

We present a microscopic derivation of the Euler GHD equation for the evolution of quasiparticle density

$$ \rho(x, \theta) = \sum_i \delta(x - x_i)\delta(\theta - \theta_i) \quad (A1) $$

Here, the quasiparticle position $x_i$ evolves according to the mapping of Eq. (35), where the bare particles' positions $X_i$ evolve with the bare velocities $\dot{X}_i = v^{\text{bare}}(\theta_i)$. We note that the motion of the quasiparticles is given by the dressed velocity $\dot{x}_i = v^{\text{eff}}_{[\rho(x_i, \cdot)]}(\theta_i)$.

Taking the time derivative of Eq. 35 we get

$$ v^{\text{bare}}(\theta_i) = \dot{X}_i = \dot{x}_i + \sum_{j \neq i} a_{ij} \delta(x_i - x_j)(\dot{x}_i - \dot{x}_j), \quad (A2) $$

where $\dot{x}_i$ is a function of all the particles' coordinates defined by Eq. (A2). Taking the time derivative of the density of particles, we get

$$ \begin{aligned} \partial_t \rho(x, \theta) &= \sum_i \partial_t (\delta(x - x_i)\delta(\theta - \theta_i)) && (A3) \\ &= -\sum_i \partial_x (\delta(x - x_i)\delta(\theta - \theta_i)) v^{\text{eff}}_{[\rho(x_i, \cdot)]}(\theta_i) \\ &= -\partial_x \left( v^{\text{eff}}_{[\rho(x, \cdot)]}(\theta) \sum_i (\delta(x - x_i)\delta(\theta - \theta_i)) \right), \end{aligned} $$

where we used the following delta function property $\partial_x \delta(x - y)f(y) = \partial_x(\delta(x - y)f(x))$. We therefore obtain

$$ \partial_t \rho(x, \theta) + \partial_x \left( v^{\text{eff}}_{[\rho(x, \cdot)]}(\theta) \rho(x, \theta) \right) = 0 \quad (A4) $$

Also, we have introduced the effective velocity

$$ \begin{aligned} v^{\text{eff}}_{[\rho(x_i, \cdot)]}(\theta_i) &= \dot{x}(x_i, \theta_i) \\ &= v^{\text{bare}}(\theta_i) \\ &+ \sum_{j \neq i} a_{ij} \delta(x_i - x_j)(\dot{x}(x_i, \theta_i) - \dot{x}(x_j, \theta_j)). \end{aligned} \quad (A5) $$
By replacing the sum with a coarsegrained integral i.e., $\sum_{i,j} \rightarrow \int dx' \int d\theta' \rho(x', \theta')..$ we get

$$v_{[\rho(x_i, \cdot)]}^{\text{eff}}(\theta_i) = v^{\text{bare}}(\theta_i) - \int d\theta' \mathfrak{a}(\theta_i - \theta') \rho(x_i, \theta')(v_{[\rho(x_i, \cdot)]}^{\text{eff}}(\theta_i) - v_{[\rho(x_i, \cdot)]}^{\text{eff}}(\theta')). \quad (A6)$$

Equations (A4) and (A6) are valid even at the microscopic level, namely, they are exact for any single trajectory of particle coordinates.

## B. Sampling initial state for the WPG

Below we sketch our procedures for sampling particle position and rapidity such that the initial state contain the correct fluctuations, which can be checked against known initial correlation functions. These procedures are carried out for homogeneous and inhomogeneous states, separately considering both classical and quantum (Fermi-Dirac) systems.

### 1. Sampling a classical homogeneous stationary state:

We here describe how to sample the initial particle positions and rapidities of an initial stationary homogeneous state, described by a given density of rapidities $\rho_{\theta}$

1. Identify the target density for a particle species, $\rho = \int d\theta \rho(\theta)$, along with a fixed length $L$ for the interacting model, such as a hard rod gas (HR). From this given target density, the bare particle density $\rho^{\text{bare}}$ is determined as

$$\rho_{\text{bare}}(\theta) = \frac{\rho(\theta)}{2\pi\rho_T(\theta)} = \frac{n(\theta)}{2\pi} \quad (B1)$$

with $\rho_T(\theta)$ being the total density of states.

2. Initialize the bare positions $X_i$: draw a distance $r_i$ from an exponential distribution weighted by the integrated bare density,

$$P^{(s)}(r) \sim e^{-r} \int d\theta \rho_{\text{bare}}(\theta), \quad (B2)$$

a procedure known as a Poisson point process. A bare particle is then placed a distance $r_i$ away from the previously placed particle, or the origin. The Poisson point process generates an initially uncorrelated state in the bare coordinates.

3. Assign a momentum to the newly placed bare particle by sampling the distribution of momenta $\rho(\theta)$.

4. Transform from bare coordinates to interacting ones by inverting Eq. (35). Break the loop if the transformed position of the new particle is outside the system size $L$; if not, return to step 2.

<figure><img src="image_8.png" /><figcaption>Figure 7. The figure shows the density profile $\rho(x)$ (top row) and the susceptibility (bottom row) i.e. $\langle \delta\rho(x', \theta_0) \delta\rho(x, \theta_0) \rangle_c = C_{\theta_0, \theta_0}(x)\delta(x - x')$ for WPG gas with the scattering kernel of the Toda model at equilibrium with $\beta = 1.0$. We confine the WPG in two different traps: harmonic (left column) and quartic (right column) with $\ell = 16$ and $\omega_0 = 1/16$ for $N = 16$ particles. These profiles (star symbols) are obtained using the sampling procedure for classical inhomogeneous systems described in the text, and the dashed lines represent the MC simulation of the Toda model. For computing correlation, we set $dx = 5$ and average over $10^5$ samples.</figcaption></figure>

### 2. Sampling a classical inhomogeneous equilibrium state:

We here describe how to correctly sample the initial particle positions and rapidities of an initial stationary inhomogeneous state, as for example the stationary state of a classical system inside an external potential $V(x)$

1. To obtain thermal density profile at a temperature $\beta$ when the system is confined in a trap $V(x)$. We use Monte-Carlo (MC) simulations, with an energy function

$$H(\{\theta_i, X_i\}) = \sum_i \frac{\theta_i^2}{2} + \sum_i V(x_i), \quad (B3)$$

where $x_i$ is a function of $\{\theta_j, X_j\}$ and obtained by inverting Eq. (35).

2. In the MC procedure, we randomly choose a bare particle and perturb its rapidity and position $(\theta_j, X_j) \rightarrow (\theta'_j, X'_j)$. Under this perturbation, we compute the energy difference $\delta H = H' - H$ where
$H$ and $H'$ are the energies of the old and new configurations, respectively.

3. Sample a random variable $r$ from uniform distribution between $r \in [0, 1]$, if $r < \exp(-\beta\delta H)$ accept the configuration otherwise reject.

4. Save the interacting particle configuration $\{\theta_i, x_i\}$ after $N$ such steps.

5. Bin the particles in the resulting interacting phase space with discretisation $d\theta dx = 2\pi$.

This protocol is verified in Fig. (7) which shows the density profile and two-point correlations for a Toda chain at thermal equilibrium inside a confining potential. We compare the result obtained from WPG and from exact Monte Carlo (MC) simulations, obtaining excellent agreement.

### 3. Sampling a quantum (fermionic) homogenous stationary state:

We here introduce a procedure analogous to the classical homogenous case but done to incorporate fermionic quantum statistics.

1. Introduce a trapping potential so that outside of a fixed region there will be no particles.

2. Here, the initial input will be the occupation function, $n(\theta, X) = n(\theta)$, determined from TBA and determines the bare particle density as

$$ \rho^{\text{bare}}(\theta) = \frac{n(\theta)}{2\pi}. \qquad (B4) $$

Chop the bare phase space into a set of bins subject to the quantization condition $d\theta dX = 2\pi$.

3. Sample a random number $p_{ij}$, whenever $p_{ij} < n(\theta_i, X_j)$ a particle is inserted into the corresponding phase space bin.

4. Transform from bare coordinates to interacting ones by inverting Eq. (28). Re-bin the particles in the resulting interacting phase space according to the same discretization $d\theta dX = d\theta dx = 2\pi$.

This protocol is verified in Fig. (2) which shows how this initialization of the state not only guarantees that the effective velocities are correct but also that its fluctuations, giving the diffusion spreadings, are.

### 4. Sampling a quantum (fermionic) inhomogeneous stationary state:

Similarly to the classical case, we consider a stationary state of a system with fermionic statistics in an external potential $V(x)$.

1. Compute the occupation function $n(\theta, x)$ using the TBA with the local-density-approximation (LDA) chemical potential $\mu \to \mu - V(x)$. The trap $V(x)$ is such that it naturally confines the particles.

2. Discretize the phase space into a set of bins with the quantization condition $d\theta dX = 2\pi$.

3. Compute the bare particle density at any bin located at $(\theta_i, X_j)$ using

$$ \rho^{\text{bare}}(\theta_i, X_j) = n(\theta_i, x) \qquad (B5) $$

where $x$ is a function of $(\theta_i, X_j)$ and is obtained by numerically solving coarsegrained Eq. (35) which is given by

$$ x = X_j - \frac{1}{2} \iint dx'd\theta' \rho(\theta', x') a(\theta_i - \theta') \text{sign}(x - x'). \qquad (B6) $$

Here $\rho(\theta, x)$ is the phase space density. Note that we use interpolated $n(\theta, x)$ in Eq. (B5).

4. For each $i$ and $j$, sample a random number $p_{ij}$ from uniform distribution between $p_{i,j} \in [0, 1]$, whenever $p_{ij} < \rho^{\text{bare}}(\theta_i, X_j)$, a particle is inserted into the corresponding phase space bin.

5. Transform from bare coordinates to interacting ones by inverting Eq. (35). Re-bin the particles in the resulting interacting phase space according to the same discretization $d\theta dX = d\theta dx = 2\pi$.

This protocol is verified in Fig. (8) which shows the density profile and two-point correlations for the Lieb-Liniger model at thermal equilibrium in the presence of an external potential. We compare the result obtained from WPG and from local density approximation (LDA), obtaining excellent agreement.

### C. Diffusive hydrodynamics from long-range correlations in inhomogeneous Hamiltonians

In this section we explicitly derive the Eq. (22) and (23) introduced in the main text. Then, we show that, in the small time limit $t \to 0^+$, they reproduce the same diffusive equation shown in [69], therefore confirming the general statement that Navier-Stokes equation only applies in the limit of short times.
<figure><img src="image_9.png" /><figcaption>Figure 8. The figure shows the density profile $\rho(x)$ (top row) and the density-density susceptibility (bottom row) i.e. $\langle\delta\rho(x', \theta_0)\delta\rho(x, \theta_0)\rangle_c = C_{\theta_0, \theta_0}(x)\delta(x - x')$ for WPG gas with the scattering kernel of the Lieb-Liniger model with $c = 1$ at equilibrium with $\beta = 1.0$ and $\mu = 1.0$. We confine the WPG in two different traps: harmonic (left column) and quartic (right column) with $\ell = 100$ and $\omega_0 = 1.0$. These profiles (star symbols) are obtained using the sampling procedure for quantum inhomogeneous systems described in the text, with $dX = 5$ and the average is over $10^5$ samples. The local density approximation (LDA) prediction is obtained by considering a gas at equilibrium with $\mu \to \mu - V(x)$.</figcaption></figure>

1. Derivation of Eq. (22) and (23)

We firstly consider the equation for the quasiparticles density $\rho$. In the same spirit as [70], we aim to write the diffusive correction to the Euler dynamics using the hydrodynamic expansion of currents. In the case of inhomogeneous forcing fields, we must consider the double projection of both currents proportional to the effective velocity and rapidity

$$ \partial_t \rho_\theta + \partial_x (v_\theta^{\text{eff}} \rho_\theta) + \partial_\theta (a_\theta^{\text{eff}} \rho_\theta) = \frac{1}{2} \partial_x (H_\theta^{\alpha\beta} E_{\alpha\beta}^{\text{sym}}) + \frac{1}{2} \partial_\theta (H_\theta^{f,\alpha\beta} E_{\alpha\beta}^{\text{sym}}). \quad (C1) $$

In particular we have defined (see Appendix's section 3 for a derivation)

$$ H_{\theta}^{\alpha\beta} = \frac{\delta^2 v_{\theta}^{\text{eff}}}{\delta\rho_{\alpha}\delta\rho_{\beta}} \quad (C2) $$

$$ = R_{\theta\gamma}^{-t} \Delta v_{\mu,\gamma}^{\text{eff}} \frac{\varphi_{\gamma\mu}^{\text{dr}}}{\rho_{T,\gamma}} (R_{\mu\alpha}^{t} R_{\gamma\beta}^{t} + R_{\mu\beta}^{t} R_{\gamma\alpha}^{t}), $$

$$ H_{\theta}^{f,\alpha\beta} = \frac{\delta^2 a_{\theta}^{\text{eff}}}{\delta\rho_{\alpha}\delta\rho_{\beta}} \quad (C3) $$

$$ = R_{\theta\gamma}^{-t} \Delta a_{\mu,\gamma}^{\text{eff}} \frac{\varphi_{\gamma\mu}^{\text{dr}}}{\rho_{T,\gamma}} (R_{\mu\alpha}^{t} R_{\gamma\beta}^{t} + R_{\mu\beta}^{t} R_{\gamma\alpha}^{t}), $$

with the definitions $\Delta v_{\mu,\gamma}^{\text{eff}} = v_{\mu}^{\text{eff}} - v_{\gamma}^{\text{eff}}$ and $\Delta a_{\mu,\gamma}^{\text{eff}} = a_{\mu}^{\text{eff}} - a_{\gamma}^{\text{eff}}$ and where we considered the point splitting regularization $E_{\alpha,\beta}^{\text{sym}} = (E_{\alpha,\beta}(x^+, x^-) + E_{\alpha,\beta}(x^-, x^+))/2$. Also, in the whole section we use Einstein index notation, assuming integration over repeated indexes and more precisely considering unintegrated any indexes appearing on both l.h.s and r.h.s.

The evolution equation for the connected two-point function reads

$$ \partial_t \langle \rho_\theta(x) \rho_{\theta'}(x') \rangle^c + \left[ \partial_x (A_\theta^\gamma \langle \rho_\gamma(x) \rho_{\theta'}(x') \rangle^c) + \partial_\theta (A_\theta^{f,\gamma} \langle \rho_\gamma(x) \rho_{\theta'}(x') \rangle^c) + (x, \theta \leftrightarrow x', \theta') \right] = 0, \quad (C4) $$

together with the matrices

$$ A_{\theta}^{\theta'} = R_{\theta,\gamma}^{-t} v_{\gamma}^{\text{eff}} R_{\gamma,\theta'}^{t}, \quad A_{\theta}^{f,\theta'} = R_{\theta,\gamma}^{-t} a_{\gamma}^{\text{eff}} R_{\gamma,\theta'}^{t}. \quad (C5) $$

In order to derive an equation for the regular part of the connected correlator we must use the decomposition

$$ \langle \rho_{\theta}(x) \rho_{\theta'}(x') \rangle^c = \delta(x - x') C_{\theta,\theta'}(x) + E_{\theta,\theta'}(x, x'), \quad (C6) $$

being $C_{\theta,\theta'} = R_{\theta,\gamma}^{-t} f_{\gamma} \rho_{\gamma} R_{\gamma,\theta'}^{-1}$, the system equilibrium susceptibility matrix. Hence, we can write

$$ \begin{aligned} & \partial_t E_{\theta,\theta'}(x,x') + \left[ \partial_x (A_{\theta}^{\gamma} E_{\gamma,\theta'}(x,x')) + \partial_\theta (A_{\theta}^{f,\gamma} E_{\gamma,\theta'}(x,x')) + (x,\theta \leftrightarrow x',\theta') \right] \\ & = \delta(x-x') \left( \partial_t C_{\theta,\theta'} + \partial_x B_{\theta,\theta'} + (\partial_\theta + \partial_{\theta'}) B_{\theta,\theta'}^f \right), \end{aligned} \quad (C7) $$

with the matrices

$$ B_{\theta,\theta'} = A_{\theta}^{\gamma} C_{\gamma,\theta'} = R_{\theta,\gamma}^{-t} v_{\gamma}^{\text{eff}} f_{\gamma} \rho_{\gamma} R_{\gamma,\theta'}^{-1}, \quad (C8) $$

$$ B_{\theta,\theta'}^{f} = A_{\theta}^{f,\gamma} C_{\gamma,\theta'} = R_{\theta,\gamma}^{-t} a_{\gamma}^{\text{eff}} f_{\gamma} \rho_{\gamma} R_{\gamma,\theta'}^{-1}. \quad (C9) $$

We now proceed by computing explicitly the term proportional to the Dirac delta function in terms of derivatives of the quasiparticles density

$$ \begin{aligned} \partial_t C_{\theta,\theta'} = & \frac{\delta C_{\theta,\theta'}}{\delta \rho_{\gamma}} \partial_t \rho_{\gamma} = \\ & - \frac{\delta (R_{\theta,\gamma}^{-t} f_{\gamma} \rho_{\gamma} R_{\gamma,\theta'}^{-1})}{\delta \rho_{\zeta}} [\partial_x (v_{\zeta}^{\text{eff}} \rho_{\zeta}) + \partial_{\zeta} (a_{\zeta}^{\text{eff}} \rho_{\zeta})], \end{aligned} \quad (C10) $$
where we use the leading order of (C1). In addition, we also compute the following derivatives (see Appendix's section 4 for a derivation)

$$ \partial_x B_{\theta, \theta'} = R_{\theta, \gamma}^{-t} f_{\gamma} \partial_x (v_{\gamma}^{\text{eff}} \rho_{\gamma}) R_{\gamma, \theta'}^{-t} + \quad (\text{C11}) $$

$$ + \frac{\delta(R_{\theta, \gamma}^{-t} f_{\gamma} R_{\gamma, \theta'}^{-1})}{\delta \rho_{\zeta}} v_{\gamma}^{\text{eff}} \rho_{\gamma} \partial_x \rho_{\zeta}, $$

$$ (\partial_{\theta} + \partial_{\theta'}) B_{\gamma, \theta'}^{\dagger} = R_{\theta, \gamma}^{-t} f_{\gamma} \partial_{\gamma} (a_{\gamma}^{\text{eff}} \rho_{\gamma}) R_{\gamma, \theta'}^{-t} \quad (\text{C12}) $$

$$ + \frac{\delta(R_{\theta, \gamma}^{-t} f_{\gamma} R_{\gamma, \theta'}^{-1})}{\delta \rho_{\zeta}} a_{\gamma}^{\text{eff}} \rho_{\gamma} \partial_{\zeta} \rho_{\zeta}. $$

Summing the results (C10), (C11) and (C12) we find

$$ \begin{aligned} \partial_t C_{\theta, \theta'} + \partial_x B_{\theta, \theta'} + (\partial_{\theta} + \partial_{\theta'}) B_{\theta, \theta'}^{\dagger} = & \\ = \frac{\delta(R_{\theta, \gamma}^{-t} f_{\gamma} R_{\gamma, \theta'}^{-1})}{\delta \rho_{\zeta}} \rho_{\gamma} R_{\zeta, \xi}^{-t} \Delta v_{\gamma, \xi}^{\text{eff}} R_{\xi, \alpha}^{t} \partial_x \rho_{\alpha} \quad (\text{C13}) \\ & + \frac{\delta(R_{\theta, \gamma}^{-t} f_{\gamma} R_{\gamma, \theta'}^{-1})}{\delta \rho_{\zeta}} \rho_{\gamma} R_{\zeta, \xi}^{-t} \Delta a_{\gamma, \xi}^{\text{eff}} R_{\xi, \alpha}^{t} \partial_{\alpha} \rho_{\alpha}, \end{aligned} $$

where we used the relation

$$ \begin{aligned} \partial_x (v_{\theta}^{\text{eff}} \rho_{\theta}) + \partial_{\theta} (a_{\theta}^{\text{eff}} \rho_{\theta}) = & \\ = R_{\theta, \gamma}^{-t} v_{\gamma}^{\text{eff}} R_{\gamma, \zeta}^{t} \partial_x \rho_{\zeta} + R_{\theta, \gamma}^{-t} a_{\gamma}^{\text{eff}} R_{\gamma, \zeta}^{t} \partial_{\zeta} \rho_{\zeta}. \quad (\text{C14}) \end{aligned} $$

Using the result (C13) into Eq. (C7), we finally find the explicit form for the r.h.s

$$ \begin{aligned} \partial_t E_{\theta\theta'}(x, x') + \partial_x (A_{\theta}^{\alpha} E_{\alpha\theta'}) + \partial_{x'} (A_{\theta'}^{\alpha} E_{\theta\alpha}) + & \\ + \partial_{\theta} (A_{\theta}^{f, \gamma} E_{\gamma\theta'}) + \partial_{\theta'} (A_{\theta'}^{f, \gamma} E_{\theta\gamma}) = & \\ = \delta(x - x') (M_{\theta\theta'}^{\alpha} \partial_x \rho_{\alpha} + M_{\theta\theta'}^{f, \alpha} \partial_{\alpha} \rho_{\alpha}) & \quad (\text{C15}) \end{aligned} $$

together with the definition

$$ M_{\theta\theta'}^{\alpha} = (R_{\theta\mu}^{-t} R_{\theta'\gamma}^{-t} + \theta \leftrightarrow \theta') \varphi_{\mu\gamma}^{\text{dr}} \frac{\rho_{\gamma} f_{\gamma}}{\rho_{T,\mu}} \Delta v_{\mu,\gamma}^{\text{eff}} R_{\mu\alpha}^{t}, \quad (\text{C16}) $$

$$ M_{\theta\theta'}^{f, \alpha} = (R_{\theta\mu}^{-t} R_{\theta'\gamma}^{-t} + \theta \leftrightarrow \theta') \varphi_{\mu\gamma}^{\text{dr}} \frac{\rho_{\gamma} f_{\gamma}}{\rho_{T,\mu}} \Delta a_{\mu,\gamma}^{\text{eff}} R_{\mu\alpha}^{t}. \quad (\text{C17}) $$

More precisely, we made use of the following variational relations

$$ \frac{\delta \rho_{\beta}}{\delta n_{\alpha}} = R_{\beta, \alpha}^{-t} \frac{n_{\alpha}}{\rho_{\alpha}}, \quad \frac{\delta R_{\alpha, \beta}^{-1}}{\delta n_{\zeta}} = -T_{\alpha, \zeta}^{\text{dr}} R_{\zeta, \beta}^{-1}. \quad (\text{C18}) $$

## 2. Limit of small times and comparison with Navier-Stokes eq. (20)

In this section we derive the diffusive equation for density of quasiparticles emerging at initial time and we show that this is equivalent with the result of [69]. More precisely, we consider uncorrelated local equilibrium initial states, i.e. with $E_{\theta, \theta'}(x, x'; t = 0) = 0$. Firstly, we solve Eq. (C15) in normal modes $E^{(n)}$, reading

$$ \begin{aligned} \partial_t (R^{-t} E_{\theta\theta'}^{(n)} R^{-1}) + & \left[ \partial_x (R^{-t} v_{\theta}^{\text{eff}} E_{\theta\theta'}^{(n)} R^{-1}) + \right. \\ & \left. + \partial_{\theta} (R^{-t} a_{\theta}^{\text{eff}} E_{\theta\theta'}^{(n)} R^{-1}) + (x, \theta \leftrightarrow x', \theta') \right] = \\ & = \delta(x - x') (M_{\theta\theta'}^{\alpha} \partial_x \rho_{\alpha} + M_{\theta\theta'}^{f, \alpha} \partial_{\alpha} \rho_{\alpha}), \quad (\text{C19}) \end{aligned} $$

where in this symbolic notation the rotation matrices are acting from left and right on the inner object. The latter equation can be solved in the small time limit through the following ansatz

$$ \begin{aligned} E_{\theta, \theta'}^{(n)}(x, x', t \to 0^+) = & \frac{1}{2} K_{\theta, \theta'}^{(I)} \text{sgn}(x - x') + \\ & + \frac{1}{2} K_{\theta, \theta'}^{(II)} \text{sgn}(x - x' + \Delta v_{\theta, \theta'}^{\text{eff}} t). \quad (\text{C20}) \end{aligned} $$

where $K^{(I)}$ and $K^{(II)}$ are two time evolving functions. This solution represent the small ballistic evolution of initial fluctuation through the characteristics, that for $t \to 0^+$ are not perturbed by the force (the acceleration is suppressed as $t^2$).

In particular, the initial time condition imposes $K^{(II)} = -K^{(I)}$. For small times, this ansatz exactly solves the r.h.s of (C19). In particular, we find

$$ K_{\theta, \theta'}^{(I)}(x, x) = \frac{1}{\Delta v_{\theta, \theta'}^{\text{eff}}} \left[ \varphi_{\theta\theta'}^{\text{dr}} \frac{\rho_{\theta} f_{\theta'}}{\rho_{T,\theta'}} \Delta v_{\theta, \theta'}^{\text{eff}} R_{\theta\alpha}^{t} \partial_x \rho_{\alpha} + (\text{C21}) \right. \\ \left. + \varphi_{\theta\theta'}^{\text{dr}} \frac{\rho_{\theta} f_{\theta'}}{\rho_{T,\theta'}} \Delta a_{\theta, \theta'}^{\text{eff}} R_{\theta\alpha}^{t} \partial_{\alpha} \rho_{\alpha} + (\theta \leftrightarrow \theta') \right]. $$

Hence, we can now write the initial time regularized correlation as

$$ \begin{aligned} R^t E_{\theta, \theta'}^{\text{sym}}(x^+, x^-; t \to 0^+) R = & \\ = \frac{1}{2} \frac{\text{sgn}(\Delta v_{\theta, \theta'}^{\text{eff}})}{\Delta v_{\theta, \theta'}^{\text{eff}}} \left[ \varphi_{\theta\theta'}^{\text{dr}} \frac{\rho_{\theta'} f_{\theta'}}{\rho_{T,\theta'}} (\Delta v_{\theta, \theta'}^{\text{eff}} R_{\theta\alpha}^{t} \partial_x \rho_{\alpha} + \right. \\ & \left. + \Delta a_{\theta, \theta'}^{\text{eff}} R_{\theta\alpha}^{t} \partial_{\alpha} \rho_{\alpha}) + (\theta \leftrightarrow \theta') \right]. \quad (\text{C22}) \end{aligned} $$

Plugging this solution into Eq. (C1) we find four additive diffusive terms involving derivatives in space and rapidity. Firstly we consider the two terms coming from the multiplication with tensor $H$

$$ \begin{aligned} H_{\theta}^{\alpha\beta} E_{\alpha, \beta}^{\text{sym}}(t \to 0^+) = & \\ = R_{\theta\gamma}^{-t} \frac{\varphi_{\gamma\mu}^{\text{dr}}}{\rho_{T,\gamma}} |\Delta v_{\mu,\gamma}^{\text{eff}}| (\rho_{\mu} f_{\mu} \frac{\varphi_{\gamma\mu}^{\text{dr}}}{\rho_{T,\mu}} R_{\gamma\alpha}^{t} - \mu \leftrightarrow \gamma) \partial_x \rho_{\alpha} & \\ + R_{\theta\gamma}^{-t} \frac{\varphi_{\gamma\mu}^{\text{dr}}}{\rho_{T,\gamma}} \frac{\Delta a_{\mu,\gamma}^{\text{eff}}}{\text{sgn}(\Delta v_{\mu,\gamma}^{\text{eff}})} (\rho_{\mu} f_{\mu} \frac{\varphi_{\gamma\mu}^{\text{dr}}}{\rho_{T,\mu}} R_{\gamma\alpha}^{t} - \mu \leftrightarrow \gamma) \partial_{\alpha} \rho_{\alpha} & \\ = -\mathfrak{D}_{\theta}^{\gamma} \partial_x \rho_{\gamma} - \mathfrak{D}_{f,\theta}^{\gamma} \partial_{\gamma} \rho_{\gamma}, & \end{aligned} \quad (\text{C23}) $$

where we used the definition of inverse rotation matrix $R^{-t}R^t = 1$ and where the factor 1/2 is compensated by the symmetry of $H$. Similarly we now consider the two terms coming from the multiplication of $E^{\text{sym}}(t = 0^+)$ with the tensor $H^\dagger$

$$ \begin{aligned} H_{\theta}^{f, \alpha \beta} E_{\alpha, \beta}^{\text{sym}}(t \to 0^+) = & \\ = R_{\theta\gamma}^{-t} \frac{\varphi_{\gamma\mu}^{\text{dr}}}{\rho_{T,\gamma}} \frac{\Delta a_{\mu,\gamma}^{\text{eff}}}{\text{sgn}(\Delta v_{\mu,\gamma}^{\text{eff}})} (\rho_{\mu} f_{\mu} \frac{\varphi_{\gamma\mu}^{\text{dr}}}{\rho_{T,\mu}} R_{\gamma\alpha}^{t} - \mu \leftrightarrow \gamma) \partial_x \rho_{\alpha} & \\ + R_{\theta\gamma}^{-t} \frac{\varphi_{\gamma\mu}^{\text{dr}}}{\rho_{T,\gamma}} \frac{(\Delta a_{\mu,\gamma}^{\text{eff}})^2}{|\Delta v_{\mu,\gamma}^{\text{eff}}|} (\rho_{\mu} f_{\mu} \frac{\varphi_{\gamma\mu}^{\text{dr}}}{\rho_{T,\mu}} R_{\gamma\alpha}^{t} - \mu \leftrightarrow \gamma) \partial_{\alpha} \rho_{\alpha} & \\ = -\mathfrak{D}_{f,\theta}^{\gamma} \partial_x \rho_{\gamma} - \mathfrak{D}_{f^2,\theta}^{\gamma} \partial_{\gamma} \rho_{\gamma}. & \end{aligned} \quad (\text{C24}) $$
Finally, using the results (C23) and (C24) into (C1) we find the diffusion equation emerging at $t \to 0^+$ for an integrable systems in inhomogeneous fields quenched from an uncorrelated local equilibrium state

$$ \partial_t \rho_\theta + \partial_x (v_\theta^{\text{eff}} \rho_\theta) + \partial_\theta (a_\theta^{\text{eff}} \rho_\theta) = \frac{1}{2} \partial_x \mathcal{D}_\theta^\gamma \partial_x \rho_\gamma + \\ + \frac{1}{2} \partial_x \mathcal{D}_{f,\theta}^\gamma \partial_\gamma \rho_\gamma + \frac{1}{2} \partial_\theta \mathcal{D}_{f,\theta}^\gamma \partial_x \rho_\gamma + \frac{1}{2} \partial_\theta \mathcal{D}_{f^2,\theta}^\gamma \partial_\gamma \rho_\gamma, \quad (\text{C25}) $$

being in perfect agreement with the result presented in [69].

### 3. Derivation of Eq. (C2) and (C3)

In this section we show the steps to derive the tensors shown in Eq. (C2) and (C3). Since the derivation in the two cases is completely analogous, we will only show how to compute $H$, meanwhile the steps for $H^\dagger$ will straightforwardly follow the same logic. Hence, we aim to compute the following variation

$$ H_{\theta}^{\theta' \theta''} = \frac{\delta A_{\theta}^{\theta'}}{\delta \rho_{\theta''}} = \frac{\delta A_{\theta}^{\theta'}}{\delta n_{\gamma}} \frac{\delta n_{\gamma}}{\delta \rho_{\theta''}}. \quad (\text{C26}) $$

Thus, we consider

$$ \frac{\delta A_{\theta}^{\theta'}}{\delta n_{\theta''}} = \left[ \frac{\delta R_{\theta, \gamma}^{-t}}{\delta n_{\theta''}} v_{\gamma}^{\text{eff}} R_{\gamma, \theta'}^{t} + R_{\theta, \gamma}^{-t} v_{\gamma}^{\text{eff}} \frac{\delta R_{\gamma, \theta'}}{\delta n_{\theta''}} \right] + \\ + \left[ R_{\theta, \gamma}^{-t} \frac{\delta v_{\gamma}^{\text{eff}}}{\delta n_{\theta''}} R_{\gamma, \theta'}^{t} \right] \equiv W_{\theta}^{\theta' \theta''} + Z_{\theta}^{\theta' \theta''}, \quad (\text{C27}) $$

where the tensors $W$ and $Z$ are respectively defined by the arguments inside the first and second square bracket. We now compute separately the two contributions using the useful relations [7, 52, 97]

$$ \frac{\delta n_{\theta}}{\delta \rho_{\theta'}} = \frac{R_{\theta, \theta'}^{t}}{\rho_{T, \theta}}, \quad \frac{\delta R_{\theta, \theta'}^{-1}}{\delta n_{\gamma}} = -\varphi_{\theta, \gamma}^{\text{dr}} R_{\gamma, \theta'}^{-1}, \quad (\text{C28}) $$

$$ \frac{\delta v_{\theta}^{\text{eff}}}{\delta n_{\theta'}} = - \frac{\rho_{T, \theta'}}{\rho_{\theta}} \Delta v_{\theta, \theta'} R_{\theta, \theta'}^{-t}. \quad (\text{C29}) $$

They explicitly read

$$ Z_{\theta}^{\theta' \theta''} = R_{\theta, \gamma}^{-t} \frac{\rho_{T, \theta''}}{\rho_{\gamma}} \Delta v_{\theta'' \gamma} R_{\gamma \theta''}^{-t} R_{\gamma \theta'}^{t}, \quad (\text{C30}) $$

$$ W_{\theta}^{\theta' \theta''} = R_{\theta\theta''}^{-t} \varphi_{\theta''\gamma}^{\text{dr}} \Delta v_{\gamma\theta''}^{\text{eff}} R_{\gamma\theta'}^{t} \quad (\text{C31}) $$

$$ H_{\theta}^{\theta' \theta''} = (W_{\theta}^{\theta' \gamma} + Z_{\theta}^{\theta' \gamma}) \frac{R_{\gamma, \theta''}^{t}}{\rho_{T, \theta''}} = \quad (\text{C32}) \\ = R_{\theta, \gamma}^{-t} \Delta v_{\mu, \gamma}^{\text{eff}} \frac{\varphi_{\gamma \mu}^{\text{dr}}}{\rho_{T, \gamma}} (R_{\mu \theta'}^{t} R_{\gamma \theta''}^{t} + R_{\mu \theta''}^{t} R_{\gamma \theta'}^{t}), $$

being the result already shown in (C2).

### 4. Derivation of Eq. (C11) and (C12)

In this section we derive Eq. (C11) and (C12). In particular, Eq. (C11) is a trivial consequence of the fact that the only independent space varying functions in the theory are $\rho(x)$ and $V(x)$, which enters in the effective velocity through the derivative of the energy density. Hence, the space derivative of any object can be simply written as variations with respect to these two functions. On the other hand, the argument to show (C12) is more involved, since in principle the rapidity is also entering non-trivially in the rotation matrix through the scattering shift $\varphi$. In this section we specialize to Galilean invariant systems where the scattering shift depends only on the difference of rapidities $\varphi_{\theta\theta'} = \varphi(\theta - \theta')$ and $\rho_{T,\theta} = 1_\theta^{\text{dr}}$. We now derive the some useful relation for the rapidity derivatives. Firstly we compute the derivative of the filling function $n_\theta$

$$ \begin{aligned} \partial_{\theta} n_{\theta} = & \frac{\partial_{\theta} \rho_{\theta}}{1^{\text{dr}}} - \frac{n_{\theta}}{1^{\text{dr}}} \varphi_{\theta, \gamma}^{\text{dr}} 1_{\gamma}^{\text{dr}} \partial_{\gamma} n_{\gamma}, \\ & (\delta_{\theta, \gamma} + n_{\theta} \varphi_{\theta, \gamma}^{\text{dr}}) 1_{\gamma}^{\text{dr}} \partial_{\gamma} n_{\gamma} = \partial_{\theta} \rho_{\theta}, \\ \partial_{\theta} n_{\theta} = & \frac{R_{\theta, \gamma}^{t}}{1^{\text{dr}}} \partial_{\gamma} \rho_{\gamma}, \end{aligned} \quad (\text{C33}) $$

where in the first line we used $R_{\theta\gamma}^{-1} = \delta_{\theta,\gamma} + \varphi_{\theta\gamma}^{\text{dr}} n_\gamma$ and that $1_\theta^{\text{dr}} = R_{\theta\gamma}^{-1} 1_\gamma$. We now can compute the derivative of a generic dressed function $f_\theta^{\text{dr}}$, where $f_\theta$ is a state independent function of rapidity

$$ \partial_{\theta} f_{\theta}^{\text{dr}} = R_{\theta, \gamma}^{-1} \partial_{\gamma} f_{\gamma} + \varphi_{\theta, \gamma}^{\text{dr}} \frac{f_{\gamma}^{\text{dr}}}{1_{\gamma}^{\text{dr}}} R_{\gamma, \alpha}^{t} \partial_{\alpha} \rho_{\alpha}, \quad (\text{C34}) $$

in which we mainly used the property $\partial_\theta \varphi_{\theta\theta'} = -\partial_{\theta'} \varphi_{\theta\theta'}$ and integration by parts. In particular, this crucially permits to compute the derivatives of the dressing function

$$ \partial_{\theta} R_{\theta\theta'}^{-t} = R_{\theta\theta'}^{-t} \partial_{\theta'} + R_{\theta\gamma}^{-t} \frac{T_{\gamma\theta'}^{\text{dr}}}{1_{\gamma}^{\text{dr}}} R_{\gamma\alpha}^{t} \partial_{\alpha} \rho_{\alpha}, \quad (\text{C35}) $$

$$ \partial_{\theta} R_{\theta\theta'}^{-1} = -\partial_{\theta'} R_{\theta,\theta'}^{-1} + T_{\theta,\gamma}^{\text{dr}} \frac{R_{\gamma,\alpha}^{t}}{1_{\gamma}^{\text{dr}}} \partial_{\alpha} \rho_{\alpha} R_{\gamma,\theta'}^{-1}. \quad (\text{C36}) $$

We finally can compute the derivative of the matrix $B^\dagger$ defined in Eq. (C9)

$$ \begin{aligned} (\partial_{\theta} + \partial_{\theta'}) B^{\dagger} = R_{\theta, \gamma}^{-t} f_{\gamma} \partial_{\gamma} (a_{\gamma}^{\text{eff}} \rho_{\gamma}) R_{\gamma, \theta'}^{-t} + & \\ + (R_{\theta\beta}^{t} R_{\theta'\gamma}^{-t} + R_{\theta\gamma}^{-t} R_{\theta'\beta}^{-t}) a_{\gamma}^{\text{eff}} \rho_{\gamma} f_{\gamma} \frac{\varphi_{\gamma\beta}^{\text{dr}}}{1_{\beta}^{\text{dr}}} R_{\beta\alpha}^{t} \partial_{\alpha} \rho_{\alpha} & \\ + R_{\theta, \gamma}^{-t} a_{\gamma}^{\text{eff}} \rho_{\gamma} f_{\gamma} \partial_{\theta'} (R_{\gamma, \theta'}^{-1} - R_{\gamma, \theta' \theta''}^{-1}), & \end{aligned} \quad (\text{C37}) $$

where the last line vanishes. Hence, we observe that the last equation is equivalent to (C12).
[1] O. A. Castro-Alvaredo, B. Doyon, and T. Yoshimura, Emergent hydrodynamics in integrable quantum systems out of equilibrium, Phys. Rev. X **6**, 041065 (2016).

[2] 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**, 207201 (2016).

[3] M. Schemmer, I. Bouchoule, B. Doyon, and J. Dubail, Generalized hydrodynamics on an atom chip, Phys. Rev. Lett. **122**, 090601 (2019).

[4] 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 (2017).

[5] J. De Nardis, D. Bernard, and B. Doyon, Hydrodynamic diffusion in integrable systems, Phys. Rev. Lett. **121**, 160603 (2018).

[6] J. D. Nardis, D. Bernard, and B. Doyon, Diffusion in generalized hydrodynamics and quasiparticle scattering, SciPost Phys. **6**, 49 (2019).

[7] B. Doyon, Lecture Notes On Generalised Hydrodynamics, SciPost Phys. Lect. Notes, **18** (2020).

[8] B. Doyon and J. Durnin, Free energy fluxes and the kubo–martin–schwinger relation, Journal of Statistical Mechanics: Theory and Experiment **2021**, 043206 (2021).

[9] B. Doyon and H. Spohn, Drude Weight for the Lieb-Liniger Bose Gas, SciPost Phys. **3**, 039 (2017).

[10] B. Doyon, Diffusion and superdiffusion from hydrodynamic projections, J. Stat. Phys. **186**, 25 (2022).

[11] J. Durnin, M. J. Bhaseen, and B. Doyon, Nonequilibrium dynamics and weakly broken integrability, Phys. Rev. Lett. **127**, 130601 (2021).

[12] B. Doyon, H. Spohn, and T. Yoshimura, A geometric viewpoint on generalized hydrodynamics, Nucl. Phys. B **926**, 570 (2018).

[13] 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**, 014001 (2022).

[14] 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**, 16 (2020).

[15] A. Bastianello, B. Doyon, G. Watts, and T. Yoshimura, Generalized hydrodynamics of classical integrable field theory: the sinh-Gordon model, SciPost Phys. **4**, 45 (2018).

[16] B. Bertini, F. Heidrich-Meisner, C. Karrasch, T. Prosen, R. Steinigeweg, and M. Žnidarič, Finite-temperature transport in one-dimensional quantum lattice models, Rev. Mod. Phys. **93**, 025003 (2021).

[17] A. Bastianello and A. De Luca, Integrability-protected adiabatic reversibility in quantum spin chains, Phys. Rev. Lett. **122**, 240606 (2019).

[18] V. B. Bulchandani, S. Gopalakrishnan, and E. Ilievski, Superdiffusion in spin chains, J. Stat. Mech. Theory Exp. **2021**, 084001 (2021).

[19] B. Doyon and F. Hübner, Ab initio derivation of generalised hydrodynamics from a gas of interacting wave packets, (2023), arXiv:2307.09307.

[20] F. Møller, P. Schüttelkopf, J. Schmiedmayer, and S. Erne, Whitham approach to generalized hydrodynamics, Phys. Rev. Res. **6**, 013328 (2024).

[21] N. Malvania, Y. Zhang, Y. Le, J. Dubail, M. Rigol, and D. S. Weiss, Generalized hydrodynamics in strongly interacting 1d bose gases, Science **373**, 1129 (2021).

[22] F. Møller, F. Cataldini, and J. Schmiedmayer, Identifying diffusive length scales in one-dimensional Bose gases, SciPost Phys. Core **7**, 025 (2024).

[23] P. Ruggiero, P. Calabrese, B. Doyon, and J. Dubail, Quantum generalized hydrodynamics, Phys. Rev. Lett. **124**, 140603 (2020).

[24] F. Møller, N. Besse, I. E. Mazets, H.-P. Stimming, and N. J. Mauser, The dissipative Generalized Hydrodynamic equations and their numerical solution, Journal of Computational Physics **493**, 112431 (2023).

[25] 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 (2018).

[26] F. Møller, S. Erne, N. J. Mauser, J. Schmiedmayer, and I. E. Mazets, Bridging effective field theories and generalized hydrodynamics, (2022), arXiv:2205.15871.

[27] A. Bastianello, V. Alba, and J.-S. Caux, Generalized hydrodynamics with space-time inhomogeneous interactions, Phys. Rev. Lett. **123**, 130602 (2019).

[28] F. Møller, B. C. Nagy, M. Kormos, and G. Takács, Anomalous charge transport in the sine–Gordon model, (2024), arXiv:2411.11473.

[29] I. Bouchoule, B. Doyon, and J. Dubail, The effect of atom losses on the distribution of rapidities in the one-dimensional Bose gas, SciPost Phys. **9**, 44 (2020).

[30] B. Doyon, S. Gopalakrishnan, F. Møller, J. Schmiedmayer, and R. Vasseur, Generalized Hydrodynamics: A Perspective, Phys. Rev. X **15**, 010501 (2025).

[31] T. Bonnemain, B. Doyon, and G. El, Generalized hydrodynamics of the kdv soliton gas, Journal of Physics A: Mathematical and Theoretical **55**, 374004 (2022).

[32] A. Bastianello, Ž. Krajnik, and E. Ilievski, Landau–lifschitz magnets: Exact thermodynamics and transport, Phys. Rev. Lett. **133**, 107102 (2024).

[33] L. Dubois, G. Théméze, F. Nogrette, J. Dubail, and I. Bouchoule, Probing the local rapidity distribution of a one-dimensional bose gas, Phys. Rev. Lett. **133**, 113402 (2024).

[34] 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 (2021).

[35] W. Kao, K.-Y. Li, K.-Y. Lin, S. Gopalakrishnan, and B. L. Lev, Topological pumping of a 1d dipolar gas into strongly correlated prethermal states, Science **371**, 296–300 (2021).

[36] Y. Le, Y. Zhang, S. Gopalakrishnan, M. Rigol, and D. S. Weiss, Observation of hydrodynamization and local prethermalization in 1d bose gases, Nature **618**, 494–499 (2023).

[37] 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 (2022).
[38] P. Schüttelkopf, M. Tajik, N. Bazhan, F. Cataldini, S.-C. Ji, J. Schmiedmayer, and F. Møller, Characterising transport in a quantum gas by measuring drude weights, (2024), arXiv:2406.17569.

[39] D. Wei, A. Rubio-Abadal, B. Ye, F. Machado, J. Kemp, K. Srakaew, S. Hollerith, J. Rui, S. Gopalakrishnan, N. Y. Yao, I. Bloch, and J. Zeiher, Quantum gas microscopy of kardar-parisi-zhang superdiffusion, Science 376, 716-720 (2022).

[40] M. Horvath, A. Bastianello, S. Dhar, R. Koch, Y. Guo, J.-S. Caux, M. Landini, and H.-C. Nägerl, Observing bethe strings in an attractive bose gas far from equilibrium, (2025), arXiv:2505.10550.

[41] L. Dubois, G. Thémèze, J. Dubail, and I. Bouchoule, Experimental investigation of a bipartite quench in a 1d bose gas, (2025), arXiv:2505.05839.

[42] B. Doyon, J. Durnin, and T. Yoshimura, The Space of Integrable Systems from Generalised $t\bar{t}$-Deformations, SciPost Phys. 13, 072 (2022).

[43] T. Bonnemain, V. Caudrelier, and B. Doyon, Hamiltonian formulation and aspects of integrability of generalised hydrodynamics, Annales Henri Poincaré 10.1007/s00023-025-01546-2 (2025).

[44] A. Aggarwal, Asymptotic scattering relation for the toda lattice, (2025), arXiv:2503.08018.

[45] B. Doyon, F. Hübner, and T. Yoshimura, New classical integrable systems from generalized $t\bar{t}$-deformations, Phys. Rev. Lett. 132, 251602 (2024).

[46] F. Møller, N. Besse, I. Mazets, H. Stimming, and N. Mauser, The dissipative generalized hydrodynamic equations and their numerical solution, Journal of Computational Physics 493, 112431 (2023).

[47] 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 (2020).

[48] L. Biagetti, G. Cecile, and J. De Nardis, Three-stage thermalization of a quasi-integrable system, Phys. Rev. Res. 6, 023083 (2024).

[49] B. Doyon, G. Perfetto, T. Sasamoto, and T. Yoshimura, Emergence of hydrodynamic spatial long-range correlations in nonequilibrium many-body systems, Phys. Rev. Lett. 131, 027101 (2023).

[50] B. Doyon, G. Perfetto, T. Sasamoto, and T. Yoshimura, Ballistic macroscopic fluctuation theory, SciPost Phys. 15, 136 (2023).

[51] T. Yoshimura and Ž. Krajnik, Anomalous current fluctuations from euler hydrodynamics, Phys. Rev. E 111, 024141 (2025).

[52] L. Biagetti, M. Lebek, M. Panfil, and J. De Nardis, Generalised bbgky hierarchy for near-integrable dynamics, To appear (2025).

[53] F. Hübner, L. Biagetti, J. D. Nardis, and B. Doyon, Diffusive hydrodynamics of hard rods from microscopics (2025), arXiv:2503.07794.

[54] S. Gopalakrishnan, A. Morningstar, R. Vasseur, and V. Khemani, Distinct universality classes of diffusive transport from full counting statistics, Phys. Rev. B 109, 024417 (2024).

[55] S. Gopalakrishnan, E. McCulloch, and R. Vasseur, Non-gaussian diffusive fluctuations in dirac fluids, Proceedings of the National Academy of Sciences 121, 10.1073/pnas.2403327121 (2024).

[56] Ž. Krajnik, E. Ilievski, and T. Prosen, Absence of normal fluctuations in an integrable magnet, Phys. Rev. Lett. 128, 090604 (2022).

[57] E. McCulloch, R. Vasseur, and S. Gopalakrishnan, Ballistic modes as a source of anomalous charge noise, Phys. Rev. E 111, 015410 (2025).

[58] O. C. Eldridge and M. Feix, Numerical experiments with a plasma model, The Physics of Fluids 6, 398-406 (1963).

[59] J. M. Dawson, Thermal relaxation in a one-species, one-dimensional plasma, The Physics of Fluids 7, 419-425 (1964).

[60] J. L. Rouet and M. R. Feix, Relaxation for one-dimensional plasma: Test particles versus global distribution behavior, Physics of Fluids B: Plasma Physics 3, 1830-1834 (1991).

[61] J.-B. Fouvry, B. Bar-Or, and P.-H. Chavanis, Kinetic theory of one-dimensional homogeneous long-range interacting systems sourced by $j\mathbf{math}_j j\mathbf{mrow}_j j\mathbf{mn}_j l_j j\mathbf{msup}_j j\mathbf{mi}_j n_j j\mathbf{mn}_j 2j j\mathbf{msup}_j j\mathbf{mrow}_j j\mathbf{math}_j$ effects, Physical Review E 100, 052142 (2019).

[62] J.-B. Fouvry, P.-H. Chavanis, and C. Pichon, Kinetic theory of one-dimensional homogeneous long-range interacting systems with an arbitrary potential of interaction, Physical Review E 102, 052110 (2020).

[63] G. B. Rybicki, Exact statistical mechanics of a one-dimensional self-gravitating system, Astrophysics and Space Science 14, 56-72 (1971).

[64] J. Reidl, Charles J. and B. N. Miller, Gravity in one dimension - quasi-equilibrium, The Astrophysical Journal 332, 619 (1988).

[65] L. Milanović, H. Posch, and W. Thirring, Statistical mechanics and computer simulation of systems with attractive positive power-law potentials, Physical Review E 57, 2763-2775 (1998).

[66] T. Tsuchiya, N. Gouda, and T. Konishi, Relaxation processes in one-dimensional self-gravitating many-body systems, Physical Review E 53, 2210-2216 (1996).

[67] B. N. Miller, K. R. Yawn, and V. P. Youngkins, Recent progress in studies of one-dimensional gravitating systems, in Dynamical Systems, Plasmas and Gravitation (Springer Berlin Heidelberg, 1996) p. 164-175.

[68] S. Gopalakrishnan, D. A. Huse, V. Khemani, and R. Vasseur, Hydrodynamics of operator spreading and quasiparticle diffusion in interacting integrable systems, Phys. Rev. B 98, 220303 (2018).

[69] J. Durnin, A. De Luca, J. De Nardis, and B. Doyon, Diffusive hydrodynamics of inhomogenous hamiltonians, Journal of Physics A: Mathematical and Theoretical 54, 494001 (2021).

[70] F. Hübner, L. Biagetti, J. De Nardis, and B. Doyon, Diffusive hydrodynamics from long-range correlations, Phys. Rev. Lett. 134, 187101 (2025).

[71] B. Doyon and T. Yoshimura, A note on generalized hydrodynamics: inhomogeneous fields and other concepts, SciPost Phys. 2, 014 (2017).

[72] X. Cao, V. B. Bulchandani, and J. E. Moore, Incomplete thermalization from trap-induced integrability breaking: Lessons from classical hard rods, Phys. Rev. Lett. 120, 164101 (2018).

[73] D. Bagchi, J. Kethepalli, V. B. Bulchandani, A. Dhar, D. A. Huse, M. Kulkarni, and A. Kundu, Unusual ergodic and chaotic properties of trapped hard rods, Phys. Rev. E 108, 064130 (2023).

[74] B. Doyon, T. Yoshimura, and J.-S. Caux, Soliton gases and generalized hydrodynamics, Phys. Rev. Lett. 120,
045301 (2018).

[75] J. L. Lebowitz, J. K. Percus, and J. Sykes, Time evolution of the total distribution function of a one-dimensional system of hard rods, Phys. Rev. **171**, 224 (1968).

[76] J. Percus, Exact solution of kinetics of a model classical fluid, *The Physics of Fluids* **12**(8), 1560–1563 (1969).

[77] C. Boldrighini, R. L. Dobrushin, and Y. M. Sukhov, One-dimensional hard rod caricature of hydrodynamics, J. Stat. Phys. **31**, 577 (1983).

[78] H. Spohn, *Large Scale Dynamics of Interacting Particles* (Springer Berlin Heidelberg, 1991).

[79] C. Boldrighini and Y. M. Sukhov, One-dimensional hard rod caricature of hydrodynamics: Navier-Stokes correction for locally-equilibrium initial states, Commun. Math. Phys. **189**, 577 (1997).

[80] B. Pozsgay, T. Gombor, and A. Hutsalyuk, Integrable hard-rod deformation of the heisenberg spin chains, Phys. Rev. E **104**, 064124 (2021).

[81] 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 (1963).

[82] M. A. Cazalilla, Bosonizing one-dimensional cold atomic gases, *Journal of Physics B: Atomic, Molecular and Optical Physics* **37**, S1 (2004).

[83] 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 (2011).

[84] I. Bouchoule and J. Dubail, Generalized hydrodynamics in the one-dimensional bose gas: theory and experiments, *Journal of Statistical Mechanics: Theory and Experiment* **2022**, 014003 (2022).

[85] J. Ford, S. D. Stoddard, and J. S. Turner, On the integrability of the toda lattice, Progress of Theoretical Physics **50**, 1547 (1973).

[86] M. Toda, Vibration of a chain with nonlinear interaction, Journal of the Physical Society of Japan **22**, 431 (1967).

[87] S. Lepri, R. Livi, and A. Politi, Too close to integrable: Crossover from normal to anomalous heat diffusion, Phys. Rev. Lett. **125**, 040604 (2020).

[88] T. Goldfriend and J. Kurchan, Equilibration of quasi-integrable systems, Phys. Rev. E **99**, 022146 (2019).

[89] M. Mestyán and V. Alba, Molecular dynamics simulation of entanglement spreading in generalized hydrodynamics, SciPost Phys. **8**, 055 (2020).

[90] 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**, 070 (2019).

[91] F. Hohl and J. W. Campbell, Statistical mechanics of a collisionless self-gravitating system, *The Astronomical Journal* **73**, 611 (1968).

[92] H. L. Wright, B. N. Miller, and W. E. Stein, The relaxation time of a one-dimensional self-gravitating system, Astrophysics and Space Science **84**, 421–429 (1982).

[93] M. Joyce and T. Worrakitpoonpon, Relaxation to thermal equilibrium in the self-gravitating sheet model, *Journal of Statistical Mechanics: Theory and Experiment* **2010**, P10012 (2010).

[94] T. N. Teles, Y. Levin, and R. Pakter, Statistical mechanics of 1d self-gravitating systems: the core—halo distribution, *Monthly Notices of the Royal Astronomical Society: Letters* **417**, L21–L25 (2011).

[95] M. Roule, J.-B. Fouvry, C. Pichon, and P.-H. Chavanis, Long-term relaxation of one-dimensional self-gravitating systems, Phys. Rev. E **106**, 044118 (2022).

[96] Y. Tang, W. Kao, K.-Y. Li, S. Seo, K. Mallayya, M. Rigol, S. Gopalakrishnan, and B. L. Lev, Thermalization near integrability in a dipolar quantum newton's cradle, Phys. Rev. X **8**, 021030 (2018).

[97] M. Medenjak, J. D. Nardis, and T. Yoshimura, Diffusion from convection, SciPost Phys. **9**, 075 (2020).