1 Introduction

The study of non-equilibrium phenomena in isolated quantum systems has received much attention in recent years [1,2,3,4]. A fruitful avenue into non-equilibrium dynamics has emerged from the study of integrable systems. These systems are characterised by a large number of conserved quantities and include several paradigmatic models, such as the Lieb-Liniger model describing a gas of one-dimensional bosons [5] and the Hubbard model describing electrons in a solid [6]. Integrability is known to fundamentally affect the dynamics of isolated quantum systems in contrast to that of generic non-integrable systems [2, 3, 7, 8]. The most prominent consequence of this is the preclusion of thermalisation, whereby a system does not relax to a thermal state. Rather, integrable systems will typically relax to nonthermal, equilibrium states, described by the so-called generalised Gibbs ensemble [9,10,11]. This has been demonstrated in various models and theories, such as conformal field theories [12, 13], the repulsive Lieb-Liniger gas [14, 15], the quantum Ising chain in a transverse field [16,17,18,19,20], and in the spin-1/2 anisotropic Heisenberg (XXZ) spin chain [21, 22].

In classical systems, thermalisation can be well understood in terms of dynamical chaos and the notions of ergodicity and mixing [23] in which a system explores its phase space uniformly and densely for almost all initial conditions. This enables a description of the system via classical statistical mechanics. However, the mechanisms of thermalisation for isolated quantum systems are not well-established. Understanding these mechanisms has been the subject of many theoretical and experimental investigations over the past two decades.

An indispensable protocol to induce non-equilibrium scenarios in isolated quantum systems is the so-called quantum quench [10, 12, 24] where an initial equilibrium state of a many-body Hamiltonian is suddenly made to evolve unitarily under another Hamiltonian. These quench scenarios constitute a vital outpost to address the lack of thermalisation in integrable quantum systems. However, as is the case for any quantum system, a direct theoretical approach with exact microscopic calculations of physically relevant observables involves the diagonalisation of the many-body Hamiltonian associated with the system. The dimension of such a many-body Hilbert space notoriously scales exponentially with the total number of particles. This often renders direct computational approaches intractable beyond more than a dozen particles. One can transcend this exponentially large space of quantum states by instead employing an aggregated, coarse-grained description of the system with a reduced number of degrees of freedom. One of the most prominent and successful such frameworks is hydrodynamics.

At the heart of any hydrodynamic theory is a set of conservation laws that govern a specific model’s dynamical degrees of freedom over vast space and time scales. These quantities typically include energy, momentum, and particle number density. Hydrodynamics has been particularly successful as an effective theory for the emergent dynamical behaviour of a variety of interacting many-body systems. For instance, Landau’s celebrated two-fluid model of superfluid helium [25], the description of electron currents in graphene [26], and in the magneto-hydrodynamic description of electrically conducting fluids [27].

Hydrodynamic theories are a valuable resource to address the mechanisms of thermalisation in isolated integrable quantum systems. Here, our interest lies in the rich many-body physics of one-dimensional integrable systems, which, besides being a mathematical curiosity, are also relevant to ultracold Bose and Fermi gases in highly anisotropic traps and to other low-dimensional condensed matter systems such as superconducting nanowires [28] and 1D polariton condensates [29]. In particular, this review focuses on the one-dimensional Bose gas. This is motivated by the high degree of experimental control over system parameters and dynamics, owing to unprecedented advancements over the past 20 years [30, 31].

By utilising the exact integrability of the one-dimensional Bose gas in the uniform limit, its equilibrium thermodynamic properties can be derived exactly using Yang-Yang thermodynamics [32]. One can then use these thermodynamic properties (more specifically, the thermodynamic equation of state for the pressure of the gas) in the equations of classical hydrodynamics that govern large-scale dynamics of the system. This approach has been dubbed “conventional hydrodynamics” (CHD) and provides an excellent description of the collective excitations of 1D bosons driven out of equilibrium in certain cases [33,34,35,36,37]. The equations of CHD take the form of Euler hydrodynamic equations, which express three conservation laws, namely of particle number (or the mass), momentum, and energy, and are valid over large space and time scales. This theory relies on the assumption of local thermal equilibrium: the system is divided into small cells, which are assumed to thermalise sufficiently fast with their environment. However, in isolated integrable and near-integrable quantum systems—wherein the 1D Bose gas is an example—this assumption of fast thermalisation (i.e. relaxation to the conventional Gibbs ensemble of statistical mechanics) is not justified; instead, the system is expected to relax to the generalized Gibbs ensemble [9,10,11], which respects the infinitely many conserved quantities of an integrable system, and not just three. Moreover, simulations of the Euler equations of conventional hydrodynamics are often plagued with the gradient catastrophe problem, preventing them from adequately describing scenarios involving, for example the formation of dispersive quantum shock waves [38,39,40].

In 2016, two papers sparked the development of a hydrodynamic theory specialised to integrable systems, known as the theory of generalised hydrodynamics (GHD) [41, 42]. In contrast to CHD, GHD does not rely on the assumption of local thermal equilibrium in the canonical Gibbs ensemble sense. Instead, systems are assumed to relax to equilibrium states, described by a generalised Gibbs ensemble (GGE)Footnote 1. GHD avoids the gradient catastrophe problem typical of scenarios involving, e.g. shock waves and, thus, can describe systems very far from equilibrium [38]. In particular, GHD has recently been shown to reproduce the most striking effects observed in the quantum Newton’s cradle experiment in a strongly interacting 1D Bose gas [38, 43], such as the undamped collisional oscillations and the lack of conventional thermalization even after hundreds of collisions. The applicability of GHD extends well beyond ultracold atomic physics. Indeed, GHD has proven to be a robust theory applicable in both classical and quantum models over large length and time scales [38, 41, 42, 44,45,46,47].

Since its original formulation, GHD has been extended to account for various experimentally relevant physical effects, such as weak integrability breaking in inhomogeneous systems [48], hydrodynamic diffusion [49, 50], and quantum fluctuations [51]. This review aims to overview the theory of generalised hydrodynamics and its recent developments, as well as its laboratory tests in ultracold atom experiments. In doing so, we restrict ourselves to the applications of GHD to describe the dynamics of the repulsive one-dimensional (1D) Bose gas described by the Lieb-Liniger model [5]. Our intention here is to give a brief yet reasonably comprehensive overview of this topical research area, intended for a broad, non-specialist audience. Upon doing so, we bring the attention of interested readers to the special issue [52] of the Journal of Statistical Physics, devoted to recent advances in GHD. The special issue includes (among many original articles) several comprehensive, in-depth review articles on GHD intended for specialists [53,54,55]; we also highlight a recent related review article by Guan et al. [56], which includes a very brief overview of GHD in the broader context of new trends in quantum integrability. We hope that our review will occupy the space between these two extremes, and that it will stimulate further interest in GHD in a broader physics community.

The organisation of this review is as follows. In Section 2, we introduce the Lieb-Liniger model and the equations of GHD at the Euler scale. We also describe the quantum Newton’s cradle experiment, its description via GHD, and the implications on the mechanisms of thermalisation. In Section 3, we outline several extensions of the original formulation of GHD. In Section 4, we review the experimental tests of GHD and the benchmarks against other theoretical approaches. Finally, Section 5 offers some perspectives on the future direction of GHD and some open problems regarding GHD.

2 Euler-scale generalised hydrodynamics

2.1 Theoretical considerations

The Lieb-Liniger Hamiltonian for N bosons in a uniform box of length L (with periodic boundary condition) interacting via the two-body contact interaction potential \(U(x,x') = g\delta (x-x')\) is given by [5]

\(\begin{aligned} H_{\text {LL}} = -\frac{\hbar ^2}{2m}\sum _{i=1}^{N} \frac{\partial ^2}{\partial x_i^2} + g\sum _{1\le i
(1)

where g quantifies the strength of interactions, assumed here to be repulsive (\(g>0\)). The Hamiltonian (1) is integrable and exactly solvable using Bethe ansatz, and as such, it admits an infinite number of conservation laws. In the second quantized form, it can be written as

\(\begin{aligned} \hat{H}_{\text {LL}}&=-\frac{\hbar ^{2}}{2m} \int dx\, \hat{\Psi }^{\dagger }(x) \frac{\partial ^{2}}{\partial x^{2}} \hat{\Psi }(x)\nonumber \\&\quad + \frac{g}{2} \int dx\, \hat{\Psi }^{\dagger }(x) \hat{\Psi }^{\dagger }(x) \hat{\Psi }(x) \hat{\Psi }(x). \end{aligned}\)
(2)

Such a one-dimensional system can be experimentally realised by confining an ultracold gas of bosons to, e.g. a highly elongated harmonic trap with transverse frequency \(\omega _{\perp }\) and axial frequency \(\omega \ll \omega _{\perp }\). When the transverse excitation energy is much larger than all other relevant energies of the problem, such as the average thermal energy and the chemical potential of the system, \(\hbar \omega _{\perp }\gg \max \{k_BT,\, \mu\)}, the transverse excitations are negligible, and the dynamics take place only along the longitudinal dimension while being frozen out in the transverse dimension. For sufficiently large systems, the boundary effects can be neglected, so that the above Hamiltonian can describe the properties of systems that are not necessarily periodic or even uniform, wherein the inhomogeneities due to the longitudinal trapping V(x) can be accounted for within the local density approximation (LDA) [57].

An important parameter encoding the strength of interactions between particles in a uniform 1D Bose gas is the dimensionless Lieb-Liniger parameter \(\gamma\) defined as

\(\begin{aligned} \gamma = \frac{mg}{\hbar ^2 n}, \end{aligned}\)
(3)

where \(n=N/L\) is the 1D (linear) density. When \(\gamma \ll 1\), the system is weakly interacting. Conversely, for \(\gamma \gg 1\), the interaction energy is large, and the system is strongly interacting. Note that one can equivalently enter the strongly interacting regime \(\gamma \gg 1\) by either increasing the interaction strength g or decreasing the density of the system n.

In addition to the dimensionless interaction strength \(\gamma\), one can also define a dimensionless temperature parameter, \(\mathcal {T}\), by scaling the temperature of the system T by the temperature of quantum degeneracy \(T_d = \hbar ^2 n^2/2mk_B\),

\(\begin{aligned} \mathcal {T} = \frac{T}{T_d} = \frac{2mk_BT}{\hbar ^2 n^2}. \end{aligned}\)
(4)

When \(T\sim T_d\), the thermal de Broglie wavelength of the particles is on the order of the mean interparticle separation. This represents the temperature regime below which quantum effects begin to dominate. Overall, finite temperature equilibrium properties of 1D Bose gas systems can be studied using the thermodynamics Bethe ansatz first derived and solved by Yang and Yang [32].

To model the dynamics of a 1D Bose gas in the hydrodynamic sense, we view the system as a continuum of mesoscopic fluid cells that are thermodynamically large but small enough to be considered spatially homogeneous. For equilibrium states, this assumption is equivalent to the LDA. Due to integrability, our system admits a set of conserved charges \(Q_i\) (\(i=1,2,3,...\)), such as energy, particle number, and momentum, each of which we assume can be written as an integral of a corresponding charge density, \(q_i(x,t)\), i.e. as

\(\begin{aligned} Q_i(t) = \int dx\, q_i(x,t). \end{aligned}\)
(5)

To each charge density, \(q_i(x,t)\), we have a corresponding current density \(j_i(x,t)\) satisfying the following continuity equation:

\(\begin{aligned} \partial _t q_i + \partial _x j_i = 0. \end{aligned}\)
(6)

The fundamental assumption of GHD is that after some relaxation time, an inhomogeneous nonstationary system approaches, within each fluid cell, states which have maximised entropy with respect to each of the conserved quantities. These maximal entropy states are described by generalised Gibbs ensembles (GGE) with the density matrix of the form:

\(\begin{aligned} \rho _{\text {GGE}}(x,t) \propto e^{-\sum _{i}\beta ^i(x,t)Q_i}, \end{aligned}\)
(7)

where \(\beta ^i(x,t)\) is the Lagrange multiplier corresponding to the charge \(Q_i\). The density matrix \(\rho _{\text {GGE}}\) is related to the average of the charge density \(q_i(x,t)\) via:

\(\begin{aligned} \mathfrak {q}_i(x,t) := \langle {q_i(x,t)}\rangle = \frac{\text {tr}\left( \rho _{\text {GGE}}\, q_i(x,t)\right) }{\text {tr}(\rho _{\text {GGE}})}. \end{aligned}\)
(8)

The set of average conserved densities \(\{\mathfrak {q}_i\}_{i=1}^{\infty }\) can be considered as a set of coordinates for the manifold of maximal entropy states. In principle, this gives a complete (coarse-grained) description of our system as the set of average conserved densities specifies a particular configuration of our system, described as a point in the manifold of maximal entropy states. An alternative set of coordinates can be obtained by utilising the quasi-particle formulation of the thermodynamic Bethe ansatz (TBA) [32, 58], in which the maximal entropy states are described via a phase space density of quasi-particles. As we will see, these coordinates offer a significant improvement over the set \(\{\mathfrak {q}_i\}\) as the density of quasi-particles satisfies a single partial differential equation, and once obtained, it is possible to construct each \(\mathfrak {q}_i\). We note that whilst the techniques of the thermodynamic Bethe ansatz enable a description of the equilibrium properties of a uniform system, the applicability of GHD extends beyond this to nonequilibrium scenarios in spatially inhomogeneous systems within the so-called local density approximation (LDA), which is implicitly assumed in any hydrodynamic theory.

The main quantity of GHD is the density of quasi-particles, also known as the rapidity distribution, denoted \(\rho _{p}(\theta ,x,t)\)Footnote 2, which carry quasi-momentum \(m\theta\) at the space-time point (x, t). This quantity is analogous to the equilibrium root density in the thermodynamic Bethe ansatz solution. However, in GHD, the rapidity distribution in our out-of-equilibrium problem evolves in time according to a classical Euler-like hydrodynamic equation. More specifically, the evolution of the rapidity distribution in GHD is governed by the set of integro-differential equations [41, 42, 59] as follows:

\(\begin{aligned} \partial _t\rho _p + \partial _x\left( v^{\text {eff}}(\theta )\rho _p\right) - \frac{1}{m}\left( \partial _x V\right) \partial _{\theta }\rho _p= 0, \end{aligned}\)
(9)

where V(x) is an external potential and \(v^{\text {eff}}(\theta )\) is the effective velocity, defined via the integral equation

\(\begin{aligned} v^{\text {eff}}(\theta ) = \theta + \int _{-\infty }^{\infty } d\theta '\, \Phi (\theta -\theta ')\left( v^{\text {eff}}(\theta ')- v^{\text {eff}}(\theta )\right) \rho _p(\theta '), \end{aligned}\)
(10)

where \(\Phi (\theta - \theta ')\) is the two-body scattering shift, which in the Lieb-Liniger model takes the following form:

\(\begin{aligned} \Phi (\theta ) = \frac{\hbar ^2}{m} \frac{2g}{g^2 + \hbar ^2\theta ^2}. \end{aligned}\)
(11)

Note we have suppressed the x and t dependence of \(\rho _p\) and \(v^{\text {eff}}\) above for simplicity. These are the main equations of GHD at the Euler scale, i.e. when quantities vary very slowly over space and time.

Note that the effective velocity is a functional of the rapidity distribution, and so the two equations are coupled. Physically, this effective velocity can be interpreted as the large-scale, coarse-grained velocity of a quasi-particle as it travels through the gas, taking into account the scattering shifts it accumulates at collisions with the other quasi-particles [60]. For the uniform Lieb-Liniger gas, the initial rapidity distribution \(\rho _p(\theta , x, 0)\) that is supplied to the GHD equations is often taken as the thermal equilibrium distribution obtained via the Yang-Yang thermodynamic Bethe ansatz [32]. For the nonuniform case, the initial rapidity distribution can be found locally for each x using the same thermodynamic Bethe ansatz but within the LDA [57].

We note that although GHD was first introduced for quantum field theories [41] and quantum chains [42], Eqs. (9) and (10) have previously been derived rigorously in the context of the classical hard rod gas [61, 62] and in soliton gases [63,64,65].

Given the rapidity distribution, \(\rho _p(\theta , x,t)\), at some fixed time, the expectation value of the conserved charge densities and their corresponding currents can then be computed via the following equations:

\(\begin{aligned} \mathfrak {q}_i(x,t)&= \int d\theta \, h_i(\theta )\, \rho _{p}(\theta ),\end{aligned}\)
(12)

\(\begin{aligned} \mathfrak {j}_i(x,t)&= \int d\theta \, h_i(\theta )\, v^{\text {eff}}(\theta )\, \rho _{p}(\theta ). \end{aligned}\)
(13)

where \(h_i(\theta )\) is the one-particle eigenvalue of the i-th conserved charge. For example, the average particle density can be obtained by setting \(h_i(\theta ) = 1\), whilst the average energy density is found by letting \(h_i(\theta ) = m\theta ^2/2\). Intuitively, Eqs. (12) and (13) can be understood as follows: the contribution of a quasi-particle of rapidity \(\theta\) to the i-th conserved charge is \(h_i(\theta )\); moreover, the number of quasi-particles with rapidity in the interval \([\theta ,\theta +d\theta )\) is given by \(\rho _p(\theta )d\theta\). Therefore, the contribution of all quasi-particles with rapidity in the range \([\theta , \theta +d\theta )\) to the i-th charge is given by \(h_{i}(\theta )\rho _{p}(\theta )d\theta\). Summing over all such rapidities yields Eq. (12). A similar interpretation holds for Eq. (13) for the average current density. We note also that the evolution equation for the rapidity distribution in Eq. (9) follows immediately by substituting Eqs. (12) and (13) into the respective continuity equations and using completeness of the set of functions \(\{h_i(\theta )\}\).

Although Eq. (13) was proposed in the original papers that introduced the theory of GHD [41, 42], its numerical verification [44, 46, 48], as well as rigorous derivations and proofs of \(v^{\text {eff}}(\theta )\) as the equation of state of GHD [66,67,68,69,70,71], was provided later. For a recent review of current operators of one-dimensional integrable models, we direct the reader to the review [72].

It is often more convenient to encode the thermodynamic properties of a system via the filling factor

\(\begin{aligned} \vartheta (\theta ,x,t) = \frac{\rho _{p}(\theta )}{\rho _p(\theta ) + \rho _h(\theta )} = \frac{\rho _p(\theta )}{\rho _s(\theta )}. \end{aligned}\)
(14)

where \(\rho _h\) is the density of "holes" and \(\rho _s(\theta )\) is the density of states which is related to the density of quasi-particles, \(\rho _p\), via the thermodynamic Bethe equation

\(\begin{aligned} 2\pi \rho _s(\theta ) = 1 + \int _{-\infty }^{\infty }d\theta '\, \Phi (\theta - \theta ')\rho _{p}(\theta '). \end{aligned}\)
(15)

One can then construct the average charge and current densities using

\(\begin{aligned} \mathfrak {q}_i(x,t)&= \int \frac{d\theta }{2\pi }\, \vartheta (\theta )\, h^{\text {dr}}_i(\theta ),\end{aligned}\)
(16)

\(\begin{aligned} \mathfrak {j}_i(x,t)&= \int \frac{d\theta }{2\pi }\, \vartheta (\theta )\, \theta \, h^{\text {dr}}_i(\theta ), \end{aligned}\)
(17)

where the dressing operation \(f\mapsto f^{\text {dr}}\) is defined via the following integral equation:

\(\begin{aligned} f^{\text {dr}}(\theta ) = f(\theta ) + \int \frac{d\theta '}{2\pi }\, \Phi (\theta ' - \theta )\vartheta (\theta ')f^{\text {dr}}(\theta '). \end{aligned}\)
(18)

In this formulation, Eq. (9) takes the following form:

\(\begin{aligned} \partial _t \vartheta + v^{\text {eff}} \partial _x \vartheta - \frac{1}{m}\left( \partial _x V\right) \partial _{\theta }\vartheta = 0, \end{aligned}\)
(19)

where the effective velocity is now outside the spatial derivative. Numerically, Eq. (19) is more convenient to work with than Eq. (9). We also advertise the open-source Matlab framework “iFluid” that numerically solves the equations of GHD [73].

2.2 The quantum Newton’s cradle experiment

Arguably, the lack of thermalisation in isolated quantum systems was best demonstrated in the hallmark quantum Newton’s cradle experiment of Kinoshita et al. [74]. There, clouds of strongly interacting rubidium atoms confined to a one-dimensional harmonic trap undergo repeated collisions without noticeably thermalising on observable time scales—corresponding to thousands of collisions. The lack of thermalisation can be attributed to the integrability of the underlying Lieb-Liniger model in the uniform limit and to weak integrability in the nonuniform (harmonically trapped) system owing to the applicability of the local density approximation [57].

In the decade that followed the Kinoshita experiment, a quantitatively accurate model of the quantum Newton’s cradle with experimentally relevant parameters remained elusive. GHD has emerged as an ideal tool to model this experiment. In 2018, Caux et al. used GHD to simulate the dynamics of a strongly interacting 1D Bose gas in the quantum Newton’s cradle experiment [43]. The evolution of the rapidity distribution is shown in Fig. 1 for the case of a harmonic potential and an anharmonic potential to mimic the trapping potential in the original quantum Newton’s cradle experiment.

Fig. 1
figure 1

Evolution of the rapidity distribution over the first ten oscillation cycles for a strongly interacting 1D Bose gas in the quantum Newton’s cradle scenario in a harmonic potential (top row) and with a small anharmonicity (middle row). Despite the presence of an integrability-breaking trapping potential and dephasing effects, neither system thermalises; the dephased states cannot be identified with a thermal state. The bottom row shows the corresponding density profiles, obtained by integrating the rapidity distribution \(\rho _p(\theta , x)\) over all rapidities \(\theta\). The blue curve corresponds to the harmonic potential, whilst the red curve is for the anharmonic trap. Adapted from [43]


Caux et al. found that even in the presence of a trapping potential that weakly breaks integrability, the 1D gas does not thermalise. Instead, the system relaxes to a generalised Gibbs ensemble. The preclusion of thermalisation here can be attributed to the existence of conserved quantities that are incompatible with convergence towards thermal equilibrium [43]. These quantities take the following form:

\(\begin{aligned} Q[f] = \int d\theta \, dx\, f(\vartheta (\theta , x, t)) \rho _p(\theta , x,t), \end{aligned}\)
(20)

where f is an arbitrary function and \(\rho _p(\theta , x, t)\) is continuous in \(\theta\) and x. The fact that these quantities are conserved under GHD evolution in a trap follows directly from Eqs. (9) and (19). We note that these quantities are only conserved at the Euler scale.

3 Beyond Euler-scale generalised hydrodynamics

Here, we present some of the extensions to the original formulation of GHD, namely, the incorporation of diffusive and quantum effects, the description of dimensional crossover in one-dimensional Bose gases, and the inclusion of space-time-dependent interactions.

3.1 Diffusive effects

Despite the excellent predictive power of standard GHD, it is often necessary to go beyond the lowest-order Euler scale. For instance, spin and charge transport in quantum chains have been shown to exhibit diffusion and other non-Eulerian behaviours not captured by Euler scale GHD [75,76,77].

As a first approximation, we assume that the state of the system and local observables at the point (x, t) can be described by the averages of conserved charges in a neighbourhood of this point. We are then permitted to take a gradient expansion of the local observables around (x, t). In particular, averages of the current densities, \(\mathfrak {j}_i(x,t)\) can be written in the following form:

\(\begin{aligned} \mathfrak {j}_i(x,t) = \mathcal {F}_{i}(x,t) + \sum _{j}\mathcal {F}_{ij}(x,t)\partial _x \mathfrak {q}_j(x,t) + \dots , \end{aligned}\)
(21)

where the space-time dependence of the functions \(\mathcal {F}_i\), \(\mathcal {F}_{ij}\), etc. is encoded through the charge densities \(\mathfrak {q}_i(x,t)\). Retaining only the first-order term \(\mathcal {F}_i\) corresponds to the so-called Euler-scale GHD. The next higher-order correction is the diffusive Navier-Stokes correction. This term introduces an arrow of time through the (irreversible) production of entropy. The corresponding Navier-Stokes diffusive GHD equation was first derived by De Nardis et al. [49] (see also [50, 78, 79]) and reads as follows:

\(\begin{aligned} \partial _t \rho _p + \partial _x\left( v^{\text {eff}}\rho _p\right) = \partial _x\left( \mathfrak {D}\partial _x\rho _{p}\right) + \frac{1}{m}\partial _x V\partial _\theta \rho _p, \end{aligned}\)
(22)

where \(\mathfrak {D}\) denotes the integral operator defined by the action,

\(\begin{aligned} \mathfrak {D}f(\theta ) = \int d\theta ' \, \mathcal {D}(\theta , \theta ') f(\theta '), \end{aligned}\)
(23)

and \(\mathcal {D}(\theta ,\theta ')\) is the diffusion kernel which satisfies the following relation:

\(\begin{aligned}{}[\mathcal {D}(\theta , \cdot )]^{\text {dr}}(\theta ')\rho _s(\theta ') = [\rho _s(\cdot )\tilde{\mathcal {D}}(\cdot , \theta ')]^{\text {dr}}(\theta ), \end{aligned}\)
(24)

with

\(\begin{aligned} \tilde{\mathcal {D}}(\theta , \theta ')&= \delta (\theta -\theta ')\left[ \int d\alpha \, \rho _p(\alpha )(1-\vartheta (\alpha ))\right. \nonumber \\&\quad \times \left. \left( \frac{\Phi ^{\text {dr}}(\alpha -\theta )}{\rho _s(\theta )}\right) ^2\left| v^{\text {eff}}(\alpha ) - v^{\text {eff}}(\theta ) \right| \right] \nonumber \\&\quad - \rho _p(\theta )(1-\vartheta (\theta ))\left( \frac{\Phi ^{\text {dr}}(\theta -\theta ')}{\rho _s(\theta )}\right) ^2\left| v^{\text {eff}}(\theta ) - v^{\text {eff}}(\theta ') \right| , \end{aligned}\)
(25)

which we note are all quantities from the Euler-scale GHD [49]. Physically, the diffusion kernel arises due to two-body scattering processes among quasi-particles [50] that are neglected in the Euler-scale hydrodynamics. Such two-body scattering processes, however, become important at smaller length scales, leading to the decay of current-current correlations and therefore to the presence of finite diffusion constants.

Fig. 2
figure 2

Simulation of diffusive GHD for a weakly interacting 1D Bose gas in the quantum Newton’s cradle setup initiated via a quench of a double-well trap to a single harmonic well. a Evolution of the density profile. b The density profile at three fixed times compared to its distribution in thermal equilibrium. In contrast to Euler-scale GHD, this system eventually thermalizes over sufficiently large time scales when diffusive effects are considered. Adapted from [78]


Using diffusive GHD, Bastianello et al. [78] studied the thermalisation of a 1D Bose gas in the quantum Newton’s cradle experiment; see Fig. 2. The authors found that diffusion was the leading mechanism that induced thermalisation towards a stationary state in the presence of an (integrability-breaking) trapping potential. The system initially relaxes to a pre-thermal state described by a GGE, and then at much longer, diffusive time scales (which scale with the length of the system L as \(\propto L^2\)), it eventually relaxes to the thermal state described by the standard Gibbs ensemble. These results demonstrate the fundamental role that diffusive effects play in the late-time dynamics of near-integrable systems.

Fig. 3
figure 3

Fermi contour of a system at zero temperature. The orange region \(\Gamma _t\) corresponds to a unit filling factor \(\vartheta (x,\theta ) = 1\). The dynamics of the system are reduced to that of the contour \(\partial \Gamma _t\). Locally, the region \(\Gamma _t\) is split into disjoint Fermi seas defined by the set of Fermi points \(\{\theta _1,\dots , \theta _{2n}\}\). Adapted from [51]


3.2 Quantum generalised hydrodynamics

In its original formulation, GHD neglects important quantum effects, such as quantum fluctuations and entanglement entropy. A fundamental assumption of GHD is that each fluid cell comprising the system is independent at any fixed time. This implies that equal-time correlations between fluid cells vanish. However, as with many quantum systems, these equal-time correlations are typically non-zero [80]; such effects occur beyond the Euler scale. One may resolve these shortcomings by re-quantising the theory of GHD. This was first achieved by Ruggiero et al. in 2020 in Ref. [51], in which they considered a system initially at zero temperature, for which quantum fluctuations are most significant. Since the entropy is initially zero and is conserved under the evolution of the standard GHD equations, the state remains at zero entropy at all times.

Zero-entropy states are characterised by the so-called Fermi contour \(\partial \Gamma _t\) enclosing the region \(\Gamma _t\) in \((x,\theta )\) phase space with unit filling

\(\begin{aligned} \vartheta (x, \theta , t) = \left\{ \begin{array}{ll} 1, &{} \text {if } (x,\theta )\in \Gamma _t,\\ 0, &{} \text {otherwise}. \end{array}\right. \end{aligned}\)
(26)

Locally, a zero-entropy state can be described by a split Fermi sea, consisting of a collection of disjoint regions in \((x,\theta )\) phase space for which the filling factor is given by

\(\begin{aligned} \vartheta (x,\theta , t) = \left\{ \begin{array}{ll} 1,&{} \text {if}\ \theta \in [\theta _1,\theta _2]\cup \dots \cup [\theta _{2p-1}, \theta _{2p}]\\ 0,&{}\text {otherwise}, \end{array}\right. \end{aligned}\)
(27)

where \(\theta _{n} \equiv \theta _n(x,t)\), \(n\in \{1,2,\dots , 2p\}\) denote the Fermi points at position x and time t, see Fig. 3.

For these zero-entropy states, it was shown by Doyon et al. that the equations of GHD, Eq. (19), reduce to a set of equations for the Fermi points [38]:

\(\begin{aligned} \partial _t \theta _n + v^{\text {eff}}\partial _x\theta _n = -\frac{1}{m}\partial _x V. \end{aligned}\)
(28)

One then considers small fluctuations about the Fermi contour. Locally, the Fermi points are modified via \(\theta _n \rightarrow \theta _n + \delta \theta _n\), where \(\delta \theta _n\) satisfies the following differential equation:

\(\begin{aligned} (\partial _t + v^{\text {eff}}\partial _x)\delta \theta _n(x,t) = 0, \end{aligned}\)
(29)

describing the propagation of linear sound waves [81].

To re-quantise this theory, one promotes the fluctuations about the Fermi contour to linear operators acting on a Hilbert space, i.e. \(\delta \theta _n \rightarrow \delta \hat{\theta }_n\). The key insight of Ruggiero et al. is that the problem of quantising sound waves above the classical GHD ground state can be recast as a problem of quantising incompressible regions in phase space, which is well-known in the literature on the quantum Hall effect [82]. The resulting re-quantised theory is a time-dependent, spatially inhomogeneous, multicomponent Luttinger liquid and has been dubbed “quantum generalised hydrodynamics”. It describes quantum fluctuations of non-equilibrium systems where conventional Luttinger liquid theory fails [51]. Additionally, unlike the Luttinger liquid theory, this quantum GHD is not restricted to low energies, and it also is no longer restricted to zero-entropy states [81, 83].

3.3 Dimensional crossover for the 1D Bose gas

One-dimensional atomic gases are often realised by using highly elongated “cigar-shaped” trapping potentials. It is assumed that due to the strong radial confinement, the transverse modes are effectively frozen out. However, when two atoms of sufficiently large momenta collide, the collisional energy may exceed the level spacing of the transverse confinement. This leads to a nonzero population of transverse excited states, which breaks integrability. To account for collisions of atoms in transverse excited states, Møller et al. introduced a phenomenological Boltzmann-type collision integral into the equations of GHD [84].

Møller et al. studied the dynamics of a quasi-1D Bose gas in the quantum Newton’s cradle scenario using this modified GHD. The evolution of the rapidity distributions under the equations of standard GHD and the extended GHD is shown in Fig. 4. By accounting for collisions with atoms in transverse excited states, the modified GHD predicts that the system thermalizes. This is in contrast to standard GHD which does not predict thermalization. The rate of thermalization of this extended GHD was also found to be consistent with previous experimental observations [85].

Fig. 4
figure 4

Evolution of the rapidity distribution for a 1D Bose gas containing 130 atoms at 94 nK in a quantum Newton’s cradle setup initiated with Bragg pulses. The top row shows the evolution during the first 100 oscillation periods according to the equations of standard GHD. The bottom row shows the evolution according to the extended model of Møller et al. that incorporates population in transversely excited states and hence leads to faster thermalization. The dashed lines mark the excitation threshold. The small fraction of transverse excited atoms has a strong influence on the dynamics of the system; the inclusion of the collision integral enables quasi-particles to redistribute across phase space and thermalize. The last two columns show the best fit for a thermal state at the temperature of the final evolved system. Adapted from [84]


A further test of GHD in the quasi-1D regime was provided by Cataldini et al. [86] which demonstrated that GHD can accurately describe the dynamics of a Bose gas whose chemical potential and thermal energy far exceed the conventional limits of one-dimensionality. In this experiment, a weakly interacting quasi-1D Bose gas, initially confined to a 1D box trap with a sinusoidal bottom, was suddenly quenched to a flat bottom trap. The evolution of the density profile was then measured and showed excellent agreement with the predictions of standard GHD. The observed agreement was attributed to the occupation of low rapidity states which, by virtue of the fermionic quasi-particle statistics, reduces the number of collision channels for transverse excitations, resulting in an emergent Pauli blocking of transverse excitations. The predictions of standard GHD and the extended GHD in this setup were almost identical. This is in contrast to the previously mentioned experiment of Møller et al. [84]. The difference here arises due to the initial double-peaked rapidity distribution in the quantum Newton’s cradle setup which permits virtual transitions to unoccupied transverse excited states, not blocked by the effective Pauli exclusion.

3.4 Space-time inhomogeneous interactions

Inhomogeneities are ubiquitous in experiments. The effects of a (spatially varying) trapping potential were first introduced into GHD by Doyon and Yoshimura in 2017 [59] in which the authors considered the addition of a generalised potential to the Hamiltonian:

\(\begin{aligned} H \rightarrow H + \sum _{k} \int dx\, V_k(x) q_k(x). \end{aligned}\)
(30)

This includes the effects of, e.g. the standard external trapping potential V(x) for \(k=0\), or perturbations by an inhomogeneous temperature field, by including the energy density \(q_1(x)\). However, there exists another relevant inhomogeneity not addressed in standard GHD, namely, that of space-time inhomogeneous interactions. Bastianello et al. were the first to incorporate these effects into GHD. They found that these inhomogeneities introduce an additional term into the standard GHD equations [48]:

\(\begin{aligned} \partial _t\rho _p + \partial _x\left( v^{\text {eff}}\rho _p\right) + \partial _{\theta }\left( \frac{f^{\text {dr}}\partial _t\alpha + \Lambda ^{\text {dr}}\partial _x\alpha }{(\partial _{\theta }p)^{\text {dr}}}\right) \nonumber \\ = \frac{1}{m}\left( \partial _x V\right) \partial _{\theta }\rho _p, \end{aligned}\)
(31)

where \(\alpha (x,t)\) is any parameter in the Hamiltonian that can be varied whilst still maintaining integrability (for instance, the interaction strength g in the Lieb-Liniger model), whereas the forces f and \(\Lambda\) are given by the following equations:

\(\begin{aligned} f(\theta )&= - \frac{\partial p(\theta )}{\partial \alpha } + \int d\theta '\, \frac{\partial \Phi (\theta -\theta ')}{\partial \alpha } \rho _p(\theta '),\end{aligned}\)
(32)

\(\begin{aligned} \Lambda (\theta )&= - \frac{\partial \varepsilon (\theta )}{\partial \alpha } + \int d\theta '\, \frac{\partial \Phi (\theta -\theta ')}{\partial \alpha } v^{\text {eff}}(\theta ')\rho _p(\theta '), \end{aligned}\)
(33)

where \(\varepsilon (\theta )\) and \(p(\theta )\) are the energy and momentum eigenvalues, respectively. We note that the above equations hold for general integrable models. For the Lieb-Liniger model, in particular, one has \(p(\theta ) = m\theta\) and \(\varepsilon (\theta ) = \frac{1}{2} m\theta ^2 - \mu\); hence, the first terms in Eqs. (32) and (33) vanish. However, such inhomogeneities are present in, for example the classical sinh-Gordon model. These results exhaust all possible inhomogeneities which can be considered on a purely hydrodynamical level [48].

4 Experimental tests and benchmarks of generalised hydrodynamics

To date, there have been four major experimental tests of GHD [84, 86,87,88]. These have all been within the context of 1D Bose gases confined to highly anisotropic trapping potentials as to realise the Lieb-Liniger model. In this section, we outline two of these experiments [87, 88] and their descriptions via GHD, while the other two experiments [84, 86] have already been mentioned in Section 3.3 in the context of dimensional crossover in quasi-1D Bose gases.

4.1 Tests of GHD in the weakly interacting regime

In 2019, Schemmer et al. experimentally demonstrated the validity of GHD for a system of one-dimensional bosons realised on an atom chip [87]. There, approximately \(6300\pm 200\) particles were confined to one dimension via magnetic trapping techniques developed by [89] and [90]. The particles were initially confined in a double-well potential at thermal equilibrium and were subsequently allowed to expand freely. The in situ density profile is shown at three times during the evolution in Fig. 5. One finds that GHD provides an excellent description of the evolving density profile, whilst conventional hydrodynamics incorrectly predicts the formation of two sharp density peaks.

Fig. 5
figure 5

Measurements of the in situ density profile for a system of \(N=6300\pm 200\) bosons, initially trapped in a 1D double-well potential and subsequently released to freely expand in 1D. These results (noisy lines) are compared against the predictions of conventional (dashed lines) and generalised hydrodynamics (solid lines). Conventional hydrodynamics incorrectly predicts two distinct density peaks at late times. By contrast, GHD shows excellent agreement with the experimental data. Adapted from [87]


Schemmer et al. also considered the dynamics of a system of 1D bosons induced via a quench from a double-well potential to a harmonic potential. They found that GHD reproduces most of the main features of the emergent dynamics; see Fig. 6. However, GHD predicts that the two density peaks persist much longer than the experimental data. Schemmer et al. attributed this discrepancy to atom losses in the experiment, such as three-body losses due to three-body recombination processes [91, 92], which break integrability. It was found that the total number of atoms in their system was reduced by 15% during the dynamics. It has been argued by Møller et al. that transverse excitations may have also influenced the observed dynamics [84]. By using the extension of GHD to the dimensional crossover regime, as outlined in Section 3.3, Møller et al. found that the density peaks were less pronounced than in standard GHD, in agreement with the experimental data.

Fig. 6
figure 6

Measurements of the in situ density profile for a system of \(N=3500\pm 140\) atoms, quenched from a 1D double-well potential to a harmonic potential [87]. These results are compared against the predictions of GHD shown as smooth lines on top of noisy experimental data. While the standard GHD reproduces most of the main features of the emergent dynamics, the minor disagreement with experimental data at later times can be rectified [84] by accounting for transverse excitations of the gas using the extension of GHD to dimensional crossover as outlined in Section 3.3. Adapted from [87]


4.2 Tests of GHD in the strongly interacting regime

In the strongly interacting regime—the opposite regime to that of the atom chip experiment of Schemmer et al.—Malvania et al. used a 2D lattice of independent harmonically trapped 1D Bose gases to test the predictions of GHD for the rapidity distributions following a strong confinement quench [88].

In this experiment, bundles of one-dimensional gases of rubidium-87 atoms were generated using optical trapping potentials. The average number of atoms per gas was on the order of 10. This ensures that the particle number density n is low, and hence, the system is well within the strongly interacting regime (\(\gamma \gg 1\)). The rapidity distribution was then measured by first switching off the axial trapping potential and letting the atoms freely expand purely in 1D until the rapidity distribution evolves into the momentum distributionFootnote 3 [93]. This is then followed by switching off the transverse confinement as well, hence allowing for the measurement of the resulting momentum distribution via the subsequent time-of-flight expansion [94]. The evolution of the space-integrated rapidity distribution [\(\rho _p(\theta ,t)=\int dx \rho _p(\theta ,x,t)\)] measured in this way is shown in Fig. 7. Note that the reported rapidity distribution is an average over all 1D tubes in the lattice. It was found that the dynamics of the gases were well-described by GHD. This is surprising as the low number of particles challenges the fundamental hydrodynamic assumption of GHD, which assumes that each fluid cell is thermodynamically large.

Fig. 7
figure 7

Comparison of experimental and GHD rapidity distributions of a strongly interacting 1D Bose gas following a strong confinement quench. The rapidity distributions are shown at different times, as indicated in each panel. The red curves show the experimental data, whilst the blue lines are the predictions of GHD, demonstrating good overall agreement between the theory and experiment. Adapted from [88]


To address the unreasonable predictive power of GHD at low particle numbers, Malvania et al. numerically simulated the experiment for a 1D Bose gas with a significantly stronger interaction strength, for which a comparison with exact results is possible. They found that the dynamics predicted by GHD were essentially identical to that obtained exactly, up to small ripples (Friedel oscillations) in the rapidity distribution, which arise at low particle numbers; see Fig. 8.

Fig. 8
figure 8

Comparison of the post-quench rapidity distributions from GHD (dashed lines) and from exact dynamics (solid lines) in the Tonks-Girardeau limit of infinite interaction strength for small atom numbers. The red, green, and blue lines correspond to 1D tubes with 5, 10, and 20 atoms, respectively. The black lines are averages of the other lines. Adapted from [88]


4.3 Benchmarks of generalised hydrodynamics

In this section, we compare the predictions of GHD against established theoretical approaches in the context of out-of-equilibrium dynamics of 1D Bose gases. Through these comparisons, we aim to elucidate the strengths and limitations of GHD.

An early comparison between GHD and CHD for the 1D Bose gas was provided by Doyon et al. in Ref. [38], in which they demonstrated that GHD surpasses the theory of conventional (or classical) Euler hydrodynamics in its ability to describe the dynamics of dissolution of an initial localised density bump on a uniform background. In Fig. 9, we show the results of such a comparison, where we see that at early times, the predictions of CHD and GHD are nearly identical. However, at longer times, CHD develops a shock that manifests itself as a derivative discontinuity or the so-called gradient catastrophe problem, common to classical Euler hydrodynamics that ignores dissipation and dispersion. This implies that Euler CHD is unable to describe the dynamics of the systems past this time. The standard GHD profile, by contrast, remains smooth and continues to describe the dynamics of dissolution past the classical shock time, despite the fact that it also ignores dissipation and dispersion effects. The reason for the success of GHD compared to CHD is that it respects the conservation of infinitely many local charges, instead of just three charges conserved in CHD which are the particle number, momentum, and energy. By doing so, GHD avoids the development of the unphysical gradient catastrophe problem.

Fig. 9
figure 9

Free evolution of a localized density bump on a uniform background in a finite temperature 1D Bose gas. The predictions of CHD (black, dashed) are compared with those of GHD (blue). The position (x) and time (t) variables are in dimensionless units corresponding to \(\hbar =m=1\), and the dimensionless temperature and interaction strength are chosen to be \(T=1\) and \(c=mg/\hbar ^2=2\). At large times, CHD develops two shocks (derivative discontinuities) on the two counter-propagating wave fronts and is unable to describe subsequent dynamics. GHD remains smooth at all times and continues to evolve past the onset of the classical shock. Adapted from [38]


Related comparisons of GHD with CHD were done also in Refs. [87, 95] (see also Fig. 5) which demonstrated similar behaviours and arrived at the same conclusions. Other benchmarks of GHD for 1D Bose gases included comparisons with full quantum simulations at zero temperature using the numerical ABACUS algorithm [38], the mean-field Gross-Pitaevskii (or nonlinear Schrödinger) equation in the weakly interacting regime [40, 47, 87, 96], exact numerics in the Tonks-Girardeau regime of infinitely strong interaction strength [40, 47, 88] (see also Fig. 8), and the matrix product state (MPS) or time-dependent density matrix renormalisation group (tDMRG) methods at intermediate interactions [47, 51, 81].

The initial benchmarks of GHD were in dynamical scenarios where generalised hydrodynamics was expected to be a valid theory. In all these scenarios, GHD has indeed demonstrated very good agreement with the alternative approaches. On the other hand, in scenarios involving small length-scale phenomena (which are not captured by GHD), it was conjectured that GHD would nevertheless adequately describe coarse-grained averages of the more accurate microscopic theories [38, 84, 87]. A specific model of such coarse-grained averaging that mimics finite imaging resolution in ultracold atom experiments was proposed and analysed by Watson et al. [47]. Watson et al. have also scrutinised the performance of GHD by benchmarking it against an array of alternative theoretical methods (including those applicable for finite temperature initial states) in some of the most challenging dynamical scenarios, such as the propagation of dispersive quantum shock waves [39, 40] and collisional dynamics in various quantum Newton’s cradle setups [74, 87].

Examples of simulations from Ref. [47] of dispersive quantum shock waves initiated from dissolution of a localised density bump are shown in Fig. 10 a–d for a weakly interacting regime of the 1D Bose gas, \(\gamma _{\textrm{bg}}\ll 1\), and in Fig. 10e in the Tonks-Girardeau regime of infinitely strong interactions (\(\gamma _{\textrm{bg}}\rightarrow \infty\)). Here, \(\gamma _{\textrm{bg}}=mg/\hbar ^2\rho _{\textrm{bg}}\) is the dimensionless interaction strength at background density \(\rho _{\textrm{bg}}\). In (a) and (b), the predictions of GHD are compared with the mean-field Gross-Pitaevskii equation (GPE) and the truncated Wigner approach (TWA) which takes into account the effect of quantum fluctuations on top of the mean-field description. In (c), the prediction of GHD at time \(\tau =0.0007\) is compared with coarse-grained averages of the GPE and TWA results and shows good agreement. Such coarse-graining that involves convolution averaging  [47] mimics finite imaging resolution in ultracold atom experiments and may explain the success of GHD when compared to experiments, even though its predictions may depart from other theoretical approaches that are valid at short wavelengths.

Fig. 10
figure 10

Dynamics of dispersive quantum shock waves forming in the evolution of an initial localized density bump in an otherwise uniform 1D Bose gas. Shown are snapshots of the evolving density profile at different dimensionless times \(\tau =t\hbar /mL^2\), where L is the length of the system, for different total atom numbers N and dimensionless interaction strength \(\gamma _{\textrm{bg}}\) in the background. Due to the reflectional symmetry about the origin, we only show the densities for \(x > 0\). See text for further discussion. Adapted from [47]


In Fig. 10d, the GHD results at \(\tau =0.0007\) are compared with the predictions of stochastic projected Gross-Pitaevskii equation (SPGPE) for a finite-temperature initial system, for two different initial dimensionless temperatures \(\overline{\mathcal {T}}\) (where \(\overline{\mathcal {T}}\!=\!T/T_d\), with \(T_d\!=\!\hbar ^2\rho _{\textrm{bg}}^2/2mk_B\) being the temperature of quantum degeneracy), again showing excellent agreement. The performance of GHD generally improves with temperature as short-wavelength interference fringes present in an equivalent GPE (\(T=0\)) simulation get washed out by thermal fluctuations so that the microscopic dynamics better conform the scope of the GHD. Finally, in Fig. 10e, we compare a snapshot of the density profile at \(\tau =0.00004\) (evolved from the same initial density bump) in the Tonks-Girardeau regime calculated using GHD and exact diagonalization (ED) of an equivalent free-fermion problem to which a Tonks-Girardeau gas can be mapped. Apart from missing density oscillations on short (microscopic)-length scales, GHD agrees very well with the ED result in this strongly interacting regime as well.

For intermediate interaction strengths, Watson et al. [47] have also compared the GHD results for the evolution of an initial bump with the results of the numerically exact MPS method for \(N=50\) particles at zero temperature. They found that the performance of GHD generally improves with stronger interactions as short-wavelength interference phenomena become increasingly suppressed due to the reduced phase coherence in the system at stronger interactions.

Watson et al. have also shown that GHD fails to capture interference phenomena at short-length scales, even in coarse-grained convolution averaging sense, for very low temperatures and very weak interaction strengths, where the local density approximation, intrinsic to any hydrodynamic approach, also fails.

Similar conclusions have been obtained from simulations of an initial density dip (see Fig. 11), known to shed a train of gray solitons in the weakly interacting regime in the mean-field GPE approximation. Here, the situation is similar to that of Fig. 10b and c, while GHD fails to capture individual solitons, whose characteristic width (on the order of the microscopic healing length) lies beyond the intended range of applicability of any hydrodynamic theory; it adequately captures the coarse-grained average density over the soliton train.

Fig. 11
figure 11

Same as in Fig. 10, except for the evolution of an initial density dip. a shows the initial density profile (\(\tau =0\)), and a snapshot of the time-evolved profiled at an early time (\(\tau =0.0005\)), evaluated using GPE, TWA, and GHD. b is a snapshot of the density profile at a later time (\(\tau =0.002\)). A train of three grey solitons can be clearly identified in the mean-field GPE results; their visibility diminishes once quantum fluctuations are taken into account through the TWA approach. While GHD fails to capture individual soliton profiles or short-wavelength structures, it agrees well with the coarse-grained averages of the GPE and TWA results, shown in c and d at \(\tau =0.0005\) and \(\tau =0.002\), respectively. e is a comparison of GHD and ED results in the Tonks-Girardeau regime, showing excellent agreement. Finally, f shows a comparison of GHD with exact MPS calculations for \(N=42\) in the nearly ideal Bose gas regime, where LDA is not applicable, hence explaining the disagreement of GHD with MPS results. Adapted from [47]


In the benchmarks of GHD in the quantum Newton’s cradle setup, performed for a double-well to single-well trap quench of a weakly interacting quasicondensate, Watson et al. [47] observed excellent agreement with the SPGPE results in both the transient dynamics and the final relaxed state, as well as in the overall characteristic thermalization time scale. However, in the comparison for the quantum Newton’s cradle initiated by Bragg pulses, they observed disagreement with SPGPE in the characteristic thermalization time scale [97], with the GHD generally predicting slower thermalization. This discrepancy, however, can be attributed to qualitatively very different ways that the Bragg pulses are implemented in GHD and the SPGPE (see Ref. [47] for further details).

5 Conclusion and perspectives

Generalised hydrodynamics has recently emerged as a broadly applicable hydrodynamic theory for modelling the quantum many-body dynamics of integrable and near-integrable systems on a large scale. Since its inception, GHD has been extended to account for various experimentally relevant effects such as diffusion, dimensional crossover in the 1D Bose gas, inhomogeneous interactions, and quantum effects. Experimentally, tests of GHD are still in their infancy. However, the first few experimental investigations show excellent agreement between the observed results and GHD predictions.

The unreasonable effectiveness of GHD at low particle numbers in [88] warrants further investigation, as do tunnelling effects between 1D tubes. Tunnelling between tubes breaks integrability and is therefore expected to lead to complete thermalization over sufficiently large time scales [98, 99]. Future experimental tests of GHD would benefit by exploring more challenging scenarios that push the assumptions of GHD, such as by straying further from integrability via long-range dipolar interactions. Experimental verification of diffusive GHD and hence characteristic thermalization rates [47, 78, 97] in different regimes of the 1D Bose gas would also be an important achievement.

In the context of verifying the required conditions for GHD, we also mention the recent work by Le et al. [100], which probed the rapid onset of hydrodynamics —referred to as hydrodynamisation—in an array of 1D Bose gases in a strongly interacting regime. Hydrodynamisation precedes local prethermalization, and it was conjectured that GHD can be applied immediately after hydrodynamisation, i.e. during local prethermalization, even though the local GGE is not yet established. If confirmed, this conjecture would imply a further relaxation of the conditions required for the applicability of GHD, which in turn would explain the success of GHD under a broader range of experimental conditions.

While we have reviewed GHD in the context of a repulsive 1D Bose gas, GHD is also applicable to the attractive 1D Bose gas [101, 102] and to other integrable models, such as the Hubbard model [103,104,105], the XXZ chain [42, 55, 106], and the Yang-Gaudin model, describing the spin-1/2 Fermi gas [107,108,109]. (For a recent review on advances in the study of transport in spin chains, sparked by the theory of GHD, we direct the reader to Ref. [55]). As such, it would be beneficial to see experimental verification of GHD in these and other condensed matter systems. Finally, it would be interesting to use GHD to model the non-equilibrium dynamics of more general integrable systems such as the multicomponent Bose gas [110] and Bose-Fermi mixtures [111].