Soliton gases and generalized hydrodynamics

TL;DR

This paper reveals a quantum-classical equivalence in generalized hydrodynamics, showing that quantum integrable systems' hydrodynamics can be described by classical gases with soliton-like scattering, enabling efficient numerical simulations.

Contribution

It introduces a classical gas model with velocity-dependent jumps that exactly reproduces the hydrodynamics of quantum integrable systems, providing a new computational approach.

Findings

Quantum GHD equations are equivalent to classical soliton-like gases.

Classical gas model allows efficient numerical simulations of quantum dynamics.

Applicable to cold-atom experiments and inhomogeneous quantum systems.

Abstract

Dynamical equations in generalized hydrodynamics (GHD), a hydrodynamic theory for integrable quantum systems at the Euler scale, take a rather simple form, even though an infinite number of conserved charges are taken into account. We show a remarkable quantum-classical equivalence: we demonstrate the equivalence between the equations of GHD, and the Euler-scale hydrodynamic equations of a new family of classical gases which generalize the gas of hard rods. In this family, the "quasi-particles", upon colliding, jump forward or backward by a distance that depends on their velocities, generalizing the jump forward by the rods' length of the fixed-velocity tracer upon elastic collision of two hard rods. Such velocity-dependent position shifts are characteristics of classical soliton scattering. The emerging hydrodynamics of a quantum integrable model is therefore that of the classical gas of its solitons. This provides a "molecular dynamics" for GHD which is numerically efficient and flexible. This is directly applicable, for instance, to the study of inhomogeneous dynamics in integrable quantum chains and in the Lieb-Liniger model realized in cold-atom experiments.

Equations

$${\partial }_t{\rho }_{\mathrm{p}}(\mathbf{\theta })+{\partial }_x(v^{\mathrm{eff}}(\mathbf{\theta }){\rho }_{\mathrm{p}}(\mathbf{\theta }))$$

$${\partial }_{t}n(\mathbf{\theta })+v^{\mathrm{eff}}(\mathbf{\theta }){\partial }_{x}n(\mathbf{\theta })$$

$$v^{\mathrm{eff}}(\mathbf{\theta })=v^{\mathrm{gr}}(\mathbf{\theta })+\int \mathrm{d}\mathbf{\alpha }\frac{\phi (\mathbf{\theta },\mathbf{\alpha })}{p^{\prime }(\mathbf{\theta })}{\rho }_{\mathrm{p}}(\mathbf{\alpha })(v^{\mathrm{eff}}(\mathbf{\alpha })-v^{\mathrm{eff}}(\mathbf{\theta }))$$

$$\mathtt{q}=\int \mathrm{d}\mathbf{\theta }{\rho }_{\mathrm{p}}(\mathbf{\theta })h(\mathbf{\theta }),\mathtt{j}=\int \mathrm{d}\mathbf{\theta }{v}^{\mathrm{eff}}(\mathbf{\theta }){\rho }_{\mathrm{p}}(\mathbf{\theta })h(\mathbf{\theta })$$

$$\phi (\theta ,\alpha )=2c/({\theta }^2+c^2)\text{mbox}(Lieb-Liniger)$$

$${\partial }_t{\rho }_{\mathrm{p}}(\mathbf{\theta })+{\partial }_x(v^{\mathrm{eff}}(\mathbf{\theta }){\rho }_{\mathrm{p}}(\mathbf{\theta }))+{\partial }_{\theta }(a^{\mathrm{eff}}(\mathbf{\theta }){\rho }_{\mathrm{p}}(\mathbf{\theta }))$$

$${\partial }_{t}n(\mathbf{\theta })+v^{\mathrm{eff}}(\mathbf{\theta }){\partial }_{x}n(\mathbf{\theta })+a^{\mathrm{eff}}(\mathbf{\theta }){\partial }_{\theta }n(\mathbf{\theta })$$

$${\partial }_t{\rho }_{\mathrm{cl}}(v)+{\partial }_x(v_{\mathrm{cl}}^{\mathrm{eff}}(v){\rho }_{\mathrm{cl}}(v))=0,$$

$$v_{\mathrm{cl}}^{\mathrm{eff}}(v)=v+d\int \mathrm{d}w{\rho }_{\mathrm{cl}}(w)(v_{\mathrm{cl}}^{\mathrm{eff}}(v)-v_{\mathrm{cl}}^{\mathrm{eff}}(w)).$$

$$\Delta x=\Delta t{v}_{\mathrm{cl}}^{\mathrm{eff}}(v).$$

$$\int \mathrm{d}wd(v,w){\rho }_{\mathrm{cl}}(w)\Delta t(v_{\mathrm{cl}}^{\mathrm{eff}}(v)-v_{\mathrm{cl}}^{\mathrm{eff}}(w)).$$

$$v_{\mathrm{cl}}^{\mathrm{eff}}(v)=v+\int \mathrm{d}wd(v,w){\rho }_{\mathrm{cl}}(w)(v_{\mathrm{cl}}^{\mathrm{eff}}(v)-v_{\mathrm{cl}}^{\mathrm{eff}}(w)).$$

$${\rho }_{\mathrm{cl}}(v)\mathrm{d}v={\rho }_{\mathrm{p}}(\theta )\mathrm{d}\theta ,v=v^{\mathrm{gr}}(\theta ),v_{\mathrm{cl}}^{\mathrm{eff}}(v)=v^{\mathrm{eff}}(\theta )$$

$$d(v,w)=-\phi (\theta ,\alpha )/p^{\prime }(\theta ).$$

$${\rho }_{\mathrm{free}}(\mathbf{v})=1-\int \mathrm{d}\mathbf{w}d(\mathbf{v},\mathbf{w}){\rho }_{\mathrm{cl}}(\mathbf{w}).$$

$${\rho }_{\mathrm{free}}(\mathbf{v})=2\pi {\rho }_{\mathrm{s}}(\mathbf{\theta })/p^{\prime }(\mathbf{\theta }).$$

$$\frac{{\rho }_{\mathrm{cl}}(\mathbf{v})\mathrm{d}v}{{\rho }_{\mathrm{free}}(\mathbf{v})}=n(\mathbf{\theta })\frac{p^{\prime }(\mathbf{\theta })\mathrm{d}\theta }{2\pi }.$$

$$v_{\mathrm{cl}}^{\mathrm{eff}}(\mathbf{v})=\frac{v^{\mathrm{gr}}(\mathbf{v})-\int \mathrm{d}\mathbf{w}d(\mathbf{v},\mathbf{w}){\rho }_{\mathrm{cl}}(\mathbf{w})v_{\mathrm{cl}}^{\mathrm{eff}}(\mathbf{w})}{1-\int \mathrm{d}\mathbf{w}d(\mathbf{v},\mathbf{w}){\rho }_{\mathrm{cl}}(\mathbf{w})}.$$

$$H=-\frac{1}{2m}\sum _{j=1}^N\frac{{\partial }^2}{\partial x_j^2}+c\sum _{j_1<j_2}\delta (x_{j_1}-x_{j_2}).$$

$$E(\lambda ,j)=j{\lambda }^2-\frac{{\bar{c}}^2}{12}j(j^2-1),p(\lambda ,j)=j\lambda .$$

$$v^{\mathrm{gr}}(\lambda ,j)=2\lambda$$

$$E(\mathbf{v})=\frac{j{v}^2}{2}-\frac{{\bar{c}}^2}{12}j(j^2-1),p(\mathbf{v})=\frac{jv}{2}.$$

$$\phi ((\lambda ,j),({\lambda }^{\prime },j^{\prime }))$$

$$a_j(\lambda )=\frac{j\bar{c}}{{\lambda }^2+(j\bar{c}/2{)}^2}.$$

$$v^{\mathrm{eff}}(\theta )=v^{\mathrm{gr}}(\theta )+\int \mathrm{d}\alpha \tilde{\phi }(\theta ,\alpha ){\rho }_{\mathrm{p}}(\alpha )(v^{\mathrm{eff}}(\alpha )-v^{\mathrm{eff}}(\theta )).$$

$$(\theta -\alpha ){\tilde{\phi }}^{\prime }(\theta ,\alpha )\le 0$$

$$(v^{\mathrm{eff}}{)}^{\prime }(\theta )=(v^{\mathrm{gr}}{)}^{\prime }(\theta )+$$

$$\left(1+\int \mathrm{d}\alpha \tilde{\phi }(\theta ,\alpha ){\rho }_{\mathrm{p}}(\alpha )\right)(v^{\mathrm{eff}}{)}^{\prime }(\theta )$$

$$1+\int \mathrm{d}\alpha \tilde{\phi }(\theta ,\alpha ){\rho }_{\mathrm{p}}(\alpha )>0$$

Full text

Introduction. It is widely believed and acknowledged that the late-time and large-scale dynamics of interacting systems, weather quantum or not, is well described by hydrodynamics. The applicability of hydrodynamics encompasses a large number of many-body systems, from classical gases and interacting quantum field theories spohn ; jeon where few hydrodynamic variables are necessary, to more exotic systems such as the classical hard-rod model bds , where velocities are preserved in each collision. Recently, the realm of hydrodynamics was extended to integrable quantum models by accounting for the infinity of nontrivial conservation laws they admit ghd ; bertini1 . In this context, on large (Eulerian) scales one assumes that the system reaches, locally at each points, a generalized Gibbs ensemble (GGE). The theory describing this was dubbed generalized hydrodynamics (GHD). This theory is applicable to many integrable models, including integrable quantum chains and field theory. In particular, it is applicable to the Lieb-Liniger model, and can therefore be used to describe the inhomogeneous dynamics in quasi-one-dimensional cold atom setups review_giamarchi such as that of the celebrated quantum Newton cradle qnc .

In this paper, we first observe that a special case of the GHD equations reproduces exactly the equations for a gas of hard rods on the line, whose dynamics is free except for elastic collisions where velocities are exchanged. The fluid equations for such a gas were mathematically derived by Boldrighini, Dobrushin, Sukhovin in 1983 bds , and the domain wall problem studied recently ds paralleling the recent solution in the quantum context from GHD ghd . We then show that a certain modification of the hard-rod dynamics leads exactly to the general form of GHD. In the modified problem, colliding rods are replaced by point-like “quasi-particles” which upon colliding jump, backward or forward, by a distance that depends on their velocities. These quasi-particles generalize the tracer of a given velocity in the hard rod problem, which jump forward by the fixed rods’ length at rods collisions. This thus gives an explicit quantum-classical correspondence at the hydrodynamic level. The velocity-dependent spatial shift of trajectories is an effect that is known to occur in classical soliton scattering. GHD can therefore be interpreted as hydrodynamic equations for a classical soliton gas, and thus should be applicable to classical field theory problems. The spatial shift exactly agrees with that of wave-packet constructions of solitons in quantum integrable models, studied and numerically observed in vlijm . This means that, from the viewpoint of Eulerian hydrodynamics, all quantum effects can be accounted for by considering the two-body classical-like scattering of soliton wave packets. This provides a nontrivial connection between the equations of state, or the hydrodynamics, of a quantum gas, and the classical soliton picture that underlies it. The soliton gas problem is extremely easy to implement on the computer, and thus this gives rise to a classical “molecular dynamics” that is numerically efficient and flexible enough to be able to account for external forces and, eventually, other effects such as integrability breaking and viscosity. The classical gas picture also gives a geometric framework allowing for exact solutions to GHD (developed in DSY ).

Generalized hydrodynamics. In the context of many-body quantum physics, hydrodynamics is a theory for the dynamics of quantum averages of local observables in non-homogeneous, non-stationary states. The principle of hydrodynamics is that when profiles of local averages become smooth enough in space-time, at every point, local averages may be evaluated, to a good approximation, within an entropy-maximized homogeneous state. Local averages are thus related to each other in the same way they are in entropy-maximized states: the equations of state hold. Entropy maximization occurs with respect to the set of conservation laws afforded by the dynamics, and the associated Lagrange parameters, which determine the local state, depend on space-time. Eulerian hydrodynamics (that is, hydrodynamics neglecting viscosity effects) is obtained simply by imposing the conservation laws for local averages of densities and currents. This thus transfers quantum dynamics into a classical dynamics for the Lagrange parameters, or for any other hydrodynamic variables (parametrizations of the state).

These concepts may be applied to one-dimensional integrable models, where infinitely many local and quasi-local conserved charges are considered in entropy maximization. In this context, local averages are to be evaluated in generalized Gibbs ensembles (GGEs) (for reviews see GGEreviews ), and the hydrodynamic theory based on GGEs is referred to as generalized hydrodynamics (GHD) ghd . The most convenient hydrodynamic variables of GHD are defined using the quasi-particle description of Bethe ansatz integrable models. Consider the set of quasi-particle states $|\theta_{1}\ldots\theta_{n}\rangle_{a_{1},\ldots,a_{n}}$. We parametrize each quasi-particle by its “rapidity” $\theta$ and specie $a$. The rapidity can be chosen as any parametrization of the two fundamental characteristics of the model, the energy $E(\bm{\theta})$ and the momentum $p(\bm{\theta})$, which form the group velocity $v^{\rm gr}(\bm{\theta})=E^{\prime}(\bm{\theta})/p^{\prime}(\bm{\theta})$ (here and below we use boldface greek letters to represent a pair of spectral parameter and particle type, $\bm{\theta}=(\theta,a)$, and the prime symbol (′) denotes rapidity derivatives ${\rm d}/{\rm d}\theta$). For instance, in relativistic models, one takes $v^{\rm gr}(\bm{\theta})=\tanh\theta$ and $p(\bm{\theta})=m_{a}\sinh(\theta)$, while in Galilean models $v^{\rm gr}(\bm{\theta})=\theta$ and $p(\bm{\theta})=m_{a}\theta$. In integrable models, the interaction can be fully characterized by the two-particle scattering amplitude $S(\bm{\theta},\bm{\alpha})$, and we define the differential scattering phase, $\varphi(\bm{\theta},\bm{\alpha})=-{\rm i}\,{\rm d}S(\bm{\theta},\bm{\alpha})/{\rm d}\theta$. We assume $\varphi(\bm{\theta},\bm{\alpha})$ to be symmetric.

In the thermodynamic limit, one represents a GGE state by its quasi-particle density $\rho_{\rm p}(\bm{\theta})$; the number of quasi-particles of type $a$ in the phase space element $[x,x+{\rm d}x]\times[p(\bm{\theta}),p(\bm{\theta})+p^{\prime}(\bm{\theta}){\rm d}\theta)]$ is $\rho_{\rm p}(\bm{\theta}){\rm d}\theta{\rm d}x$. This quantity fully determines the GGE state, and any local average can be evaluated in terms of $\rho_{\rm p}(\bm{\theta})$. Within this picture, one also introduces the density of states $\rho_{\rm s}(\bm{\theta})$ defined by $2\pi\rho_{\rm s}(\bm{\theta})=p^{\prime}(\bm{\theta})+\int{\rm d}\bm{\alpha}\,\varphi(\bm{\theta},\bm{\alpha})\rho_{\rm p}(\bm{\alpha})$, and the occupation function characterizing the proportion of available states that are occupied by quasi-particles, $n(\bm{\theta})=\rho_{\rm p}(\bm{\theta})/\rho_{\rm s}(\bm{\theta})$.

It was shown in ghd ; bertini1 that the infinity of hydrodynamic conservation laws of GHD lead to the following equivalent continuity and convective equations:

[TABLE]

where the effective velocity $v^{{\rm eff}}(\bm{\theta})$ solves the equation

[TABLE]

(here and below $\int{\rm d}\bm{\theta}=\sum_{a}\int_{\mathbb{R}}{\rm d}\theta$, and we suppress the explicit $x,t$ dependence for lightness of notation). The effective velocity $v^{{\rm eff}}(\bm{\theta})$ bonnes ; ghd ; bertini1 serves as a velocity of the quasi-particle $\bm{\theta}$ in the fluid cell at $(x,t)$ (microscopically, it is the velocity of the elementary excitation carrying the state $\bm{\theta}$). The quasi-particle spectral density $\rho_{\rm p}(\bm{\theta})$ is a conserved fluid density, while $n(\bm{\theta})$ is a convected fluid variable. In particular, one may evaluate averages of any local conserved density $q$ and its current $j$ via ghd ; bertini1

[TABLE]

where $h(\bm{\theta})$ is the quantity of charge carried by the quasi-particle $\bm{\theta}$ (the one-particle eigenvalue of the associated conserved charge).

For instance, in the Lieb-Liniger (LL) model LLmodel , a Galilean invariant model for interacting Bose gases experimentally realizable review_giamarchi , there is a single particle specie in the repulsive regime, and the differential scattering phase takes the form

[TABLE]

where $c$ is the coupling strength (see the Supplementary Material (SM) for a description of the LL model in the attractive regime). The observable for the number of particles in the gas corresponds, in the above formulae, to the choice $h(\bm{\theta})=1$. In the LL model, Eqs. (1)-(4) were derived in ghd , and in the XXZ quantum spin chains, they were first obtained in bertini1 .

GHD was also considered for evolution within inhomogeneous fields, such as force or temperature fields dy . We may represent such a field by simply taking the energy to be space-dependent $E(\bm{\theta})=E(x;\bm{\theta})$. Then the following equivalent equations hold:

[TABLE]

where $v^{\rm eff}(\bm{\theta})$ is still given by (3), but with space-dependent group velocity. The effective acceleration is the time derivative $a^{\rm eff}(\bm{\theta})={\rm d}v^{\rm eff}(\bm{\theta})/{\rm d}t|_{\rho_{\rm p}}$ under the time change of the group velocity due to the force $F(\bm{\theta})=-\partial_{x}E(\bm{\theta})$ at constant fluid state $\rho_{\rm p}$, obtained from differentiating (3) with ${\rm d}v^{\rm gr}(\bm{\theta})/{\rm d}t=F(\bm{\theta})/p^{\prime}(\bm{\theta})$. It simplifies to the ordinary accleration $a^{\rm eff}(\bm{\theta})=F(\bm{\theta})/p^{\prime}(\bm{\theta})$ when the latter quantity is independent of $\bm{\theta}$. In the single-specie Galilean case with a force potential $V$, such as in the repulsive LL model, we have $E(\theta)=m\theta^{2}/2+V$ and $a^{\rm eff}(\theta)=-\partial_{x}V/m$. See dy for details.

Classical gases. Let us first recall the classical hard rod model spohn ; bds . Rods (non-intersecting one-dimensional segments) of a fixed length $d$ move inertially at various velocities $v$ on the infinite line, except for elastic collisions at which they exchange their velocities. The emergence of hydrodynamic equations on large scales in this model for a large class of initial conditions was rigorously demonstrated bds . Let $\rho_{\rm cl}(v)$ be the density of rods with velocity $v$ (that is, $\rho_{\rm cl}(v)\,{\rm d}x\,{\rm d}v$ is the number of rod centers lying within the phase space element $[x,x+{\rm d}x]\times[v,v+{\rm d}v]$). The hydrodynamic equations read

[TABLE]

where $v^{\rm eff}_{\rm cl}(v)$ satisfies bds

[TABLE]

These are exactly the same equations as (3) and (1), in the Galilean case ($v=\theta$), without inhomogeneous fields, with a single unit-mass particle specie, and with the identifications $d=-\varphi(v,w)$, $\rho_{\rm p}(v)=\rho_{\rm cl}(v)$, and $v^{\rm eff}(v)=v^{\rm eff}_{\rm cl}(v)$. In the quantum model, this corresponds to a purely exponential scattering phase, $S(\theta,\alpha)=e^{-{\rm i}d\,(\theta-\alpha)}$. This is the non-relativistic limit of a scattering phase which was considered in GL in a different context.

The above observation suggests that if we allow the rods to collide more “softly”, in such a way that $d$ becomes velocity dependent, the hydrodynamics of the emerging gas might be identical to that of GHD.

Rod lengths are distances at which two rod centers exchange their velocities. One could attempt to directly extend this to velocity dependent rod lengths: neighboring rods with velocities $v$ (left) and $w$ (right), with $v>w$, exchange their velocities when their centers are at distance $d(v,w)$ (as if rods were elastically contractible instead of hard). However this causes difficulties with respect to many-body scattering: after the instantaneous exchange of velocities, other rods in the immediate neighborhood might be nearer than the velocity-dependent minimal distance. It also causes difficulties for negative lengths. Instead, let us modify the microscopic dynamics to one that is more easily generalizable, keeping the same hydrodynamics.

Consider the trajectory of a velocity-tracer, following the center of the rod carrying velocity $v$. Its trajectory is that of a free particle, except for jumps by a distance $d$ when rod collisions occur. The tracer is a quasi-particle of velocity $v$. In this language, quasi-particles are point-like objects, and two quasi-particles at positions $x_{1}<x_{2}$ “collide” if their positions satisfy $x_{2}-x_{1}=d$ and their velocities $v_{1}>v_{2}$, at which instant $x_{1}\mapsto x_{1}+d$ and $x_{2}\mapsto x_{2}-d$, and quasi-particle trajectories cross. The important property of this microscopic dynamics is that every crossing of two quasi-particles’ trajectories comes with such a shift by $d$, in microscopic time (here, instantaneously). As is clear from the derivation below, any microscopic dynamics with this property will lead to the same hydrodynamics.

We may therefore modify this classical gas by proclaiming collisions to occur at $x_{2}=x_{1}$ instead. At a collision, the involved quasi-particles instantaneously jump, like fleas, by a distance $d$ as above. The jump is “forward”: the quasi-particle on the left (right) jumps towards the right (left). This classical “flea gas” is generalized easily: a collision between quasi-particles of velocities $v$ and $w$ occasion a jump by $d(v,w)$, forward if positive, backward if negative. A jump should be seen as an infinitely-fast displacement, during which more collisions can occur, occasioning new jumps in a chain reaction that re-organizes the quasi-particles’ positions in the local neighborhood. See the SM for the precise flea-gas algorithm and other microscopic procedures.

We now argue that this reproduces GHD. We are looking for the effective velocity $v^{\rm eff}_{\rm cl}(v)$ of a test quasi-particle, defined through the actual distance that it travels in a macroscopic time $\Delta t$,

[TABLE]

The gas is characterized by the density $\rho_{\rm cl}(v)$, and by standard arguments, the continuity equation holds, $\partial_{t}\rho_{\rm cl}(v)+\partial_{x}(v^{\rm eff}_{\rm cl}(v)\rho_{\rm cl}(v))=0$. The quantity $\Delta x$ results from the linear displacements at velocity $v$, given by $\Delta t\,v$, along with the accumulation of jumps the quasi-particle undergoes as it travels through the gas. The oriented distance jumped due to hitting a quasi-particle that has velocity $w$ is ${\rm sign}(v-w)d(v,w)$. The average number of quasi-particles of velocity between $w$ and $w+{\rm d}w$ that has been crossed, is the total number ${\rm d}w\,\rho_{\rm cl}(w)\Delta x$ present within the length $\Delta x$, times the probability $\Delta t/\Delta x\times|v^{\rm eff}_{\rm cl}(v)-v^{\rm eff}_{\rm cl}(w)|$ that the test particle crosses such a quasi-particles in time $\Delta t$. Assuming that the effective velocity is monotonic with $v$ (see the SM), the total jumped distance is therefore

[TABLE]

Equating (10) with the sum of $\Delta t\,v$ and (11) we find

[TABLE]

It is clear that the GHD equations (1), in the case of a single specie, reproduce the hydrodynamics of the above classical gas under the following identification:

[TABLE]

along with

[TABLE]

The above derivation is readily generalizable to many species, with, in (12), the velocity parameters $v,w$ replaced by doublets $\bm{v}=(v,a),\,\bm{w}=(w,b)$, and the driving velocity value $v$ replaced by the appropriate value $v^{\rm gr}(\bm{v})$ for specie $a$. We recover (3) by reparametrization. It is also clear that (6) holds in the single-specie case if an external force field is applied, as the acceleration times $\rho_{\rm cl}(\bm{\theta})$ directly gives the current of quasi-particles in velocity space. We leave the general case of (6) for a future work.

The quantum-classical dictionary. The above shows that GHD, derived in quantum integrable models using quantum integrability, also describes the hydrodynamic limit of classical models of integrable dynamics: the hard rod gas and the flea gas generalization. Eqs. (13) and (14) give a quantum-classical dictionary. Further elements of the dictionary are as follows. Consider the free space fraction (in the hard rod gas, the fraction of a unit length where there is no rod at all; in the general case, that available omitting the distances jumped if forward, or adding them if backward):

[TABLE]

We recognize this quantity as, up to a factor, the density of state in the quantum problem,

[TABLE]

For the occupation function $n(\bm{\theta})$, we then have

[TABLE]

The left-hand side is the classical number of quasi-particles per unit length of free space. Note that for negative jumped distances, the “fraction” of a unit length that is free is greater than one: the effect of quasi-particle scattering is to increase the space available.

The classical picture provides a clear interpretation of the effective velocity. Let us write it in the form

[TABLE]

Suppose $d(v,w)<0$. In this case, the gas slows down a test quasi-particle, as it is affected by backwards jumps at collisions. The denominator accounts for this friction effect. If the gas also carries a drift, then the test particle is carried by the gas itself for the same reason. The second term in the numerator accounts for this. These two competing effects were noticed in ghd when studying the steady state emerging from by a domain-wall initial condition in the sinh-Gordon and Lieb-Liniger models. In the case $d(v,w)>0$, the opposite occurs.

Finally, it is worth noting that in the classical problem, the rod length $d(v,w)$ can be chosen at will, and in particular the full dynamics is invariant under simultaneous scaling of space, time and rod lengths. In the quantum case, there is no such freedom: the differential scattering phase $\varphi(\bm{\theta},\bm{\alpha})$ is fixed by the scattering of quasi-particles (the a fundamental length scale is determined by $\hbar$).

Soliton gases. The above intriguing quantum-classical correspondence might be explained in terms of soliton gases. In classical soliton scattering, two solitons retain, asymptotically, their form and their speeds, the only change being in shifts of their trajectories. These shifts are velocity dependent, and thus, up to microscopic differences that do not affect the large-scale hydrodynamics, the flea-gas dynamics is equivalent to classical soliton scattering. Repulsive (attractive) soliton scattering correspond to positive (negative) position shifts, in agreement with the intuition from the flea gas. It is tempting to propose that GHD be the hydrodynamics arising in models of classical integrable field theory possessing soliton modes.

GHD, however, was derived in quantum integrable systems. In this context, quasi-particles are often interpreted as the quantization of classical solitons. This interpretation was recently numerically observed by forming wave packets of quasi-particle excitations in the Heisenberg quantum chain vlijm . It was confirmed that the trajectory shifts are given by the differential scattering phase of the quantum model, in agreement with (14). Therefore, we conclude that at the Euler scale, quantum gases can be seen as the gas of their classical soliton wave packets. This is in agreement with the picture according to which in quantum models, any multi-particle scattering process may be seen as a sequence of well separated two-body scattering processes (at the basis of the (generalized) thermodynamic Bethe ansatz GTBA ).

Numerical checks. We have numerically simulated the classical gas corresponding to the LL model (5) with $c=1$, $m=1$, without force field. We have verified the form of the effective velocity by evaluating explicitly, in a homogeneous stationary gas, the total displacement of a test quasi-particle divided by the time spent, and comparing with the result of solving numerically the integral equation (18). See Fig. 1a; the agreement is excellent. This is well within the strongly interacting regime, with LL coupling parameter $\gamma=mc\rho^{-1}\approx 1.1$ ($\rho=\int{\rm d}v\,\rho_{\rm cl}(v)$ is the particle density). We have also implemented a domain wall initial condition in the LL model, and checked that its dynamics reproduces the self-similar solution derived in ghd . See Fig. 1b, as well as Fig. 3 in the SM.

Conclusion. We have developed a classical gas dynamics that reproduces, at the Euler scale, the equations of GHD for arbitrary differential scattering phase. This gives an efficient way of simulating full space-time dependent profiles solving GHD including with force fields, complementing exact solutions found in BVKM ; DSY . It is applicable to experimentally relevant models such as the LL model from any state and in any trap potential, or quantum spin chains. It also shows that the GHD equations describe classical soliton gases, and it provides a clear interpretation of friction and drag effects. In the present formulation, scattering must be diagonalized before it is classically implemented; however perhaps non-diagonal scattering could be directly implemented in the classical gas. It would be interesting to analyze the large-deviation theory of these new classical gases, especially their relation with that of quantum problems, as well as classical gases reproducing critical quantum models (see bdreview ).

Our observations may have far-reaching implications. Soliton gases may be seen as wide generalizations of the semiclassical picture proposed in SY , and may lead to efficient ways of evaluating correlations in certain regimes. Viscosity or other higher-derivative effects in the quantum problems will have many sources, including the finite scattering length taken into account by the classical gas, but also soliton wave packet spreading. By appropriately modifying the classical algorithm, it might be possible to phenomenologically account for such corrections to GHD, as well as for integrability-breaking processes, which would otherwise be extremely difficult to numerically implement.

Acknowledgment. We thank Robert Konik and Herbert Spohn for useful discussions. Hospitality is acknowledged as follows: all authors thank SISSA and ICTP, Trieste (visits and workshop “Entanglement and non-equilibrium physics of pure and disordered systems”, July 2016), and BD and TY thank City University New York (workshop “Dynamics and hydrodynamics of certain quantum matter”, March 2017). TY is grateful for the support from the Takenaka Scholarship Foundation, and thanks Tokyo Institute of Technology for hospitality.

I Flea gas algorithm

The gas is represented by a chain of cells, each cell representing a quasi-particle $A$ and containing all necessary information pertaining to it (its velocity, its type, its position $A.{\tt pos}$ on the line), ordered from left to right. The procedure distance$(A,B)$ returns the oriented distance of the jump for a collision of $A$ against $B$ (it depends on the velocities and types of $A$ and $B$); it was denoted $d(\cdot,\cdot)$ in the text, and it is positive for forward jumps, negative for backward jumps. In order not to perform a given jump, associated to a give crossing, twice, we need to mark the pairs of quasi-particle at their first collision; marked pairs, if they collide another time, just past through each other. This is of course essential if jumps are backwards, and in general has the effect that under the re-organization of quasi-particles’ positions in a neighborhood of a collision, the exact distance $d(\cdot,\cdot)$ has been jumped by every quasi-particle affected by a collision. A picture for possible collisions between quasi-particles with velocities $v$ and $w$ are depicted in the Fig.2, and a precise algorithm for the flea gas is as follows.

procedure Evolve:

Displace all until next collision;

Collide left-particle against right-particle;

repeat until evolution time has elapsed;

end.

procedure Collide $A$ against $B$:

if marked $(A,B)$:

Exchange $A$, $B$ in the chain

else:

Mark $(A,B)$;

Jump $A$ by distance$(A,B)$;

Jump $B$ by $-$distance$(A,B)$;

end.

procedure Jump $A$ by $D$:

if $D<0$ then side is left;

else side is right;

repeat:

$B$ := quasi-particle to the side of $A$;

if $B$ exists and $|A.{\tt pos}-B.{\tt pos}|<|D|$:

$D:=D-B.{\tt pos}+A.{\tt pos}$;

$A$.pos := $B$.pos;

if side is left: Collide $B$ against $A$;

else: Collide $A$ against $B$;

else:

$A$.pos := $A$.pos + $D$;

break;

end.

Above, side refers to the side within the chain. Note that, because of the recursive process (the chain reaction of collisions), the event Jump $A$ by distance$(A,B)$ may also affect the position of $B$. This is fine, as long as $B$ then jumps by the correct distance afterwards. Of course, the choice of making $A$ jump before $B$ in the Collide procedure is arbitrary. Different orderings lead to different microscopic re-organizations of quasi-particles’ positions at collisions, but to the same large-scale Euler hydrodynamics.

An acceleration due to an external field is implicitly implemented with the Displace procedure. There, the evolution of each quasi-particle changes, in general, both its position and its velocity, as per the usual physical laws for particles within external force fields $F(\bm{\theta})$.

Other algorithms are possible. In particular, it is not necessary to perform jumps “instantaneously” at collisions; as long as the distance $d(\cdot,\cdot)$ is added to the trajectory in a microscopic time, the large-scale effect is the same. For instance, one may associate to each quasi-particle a “ghost velocity change” $\Delta v$, used for time evolution in order to add an appropriate distance, but not for the calculation of the jump distance (which uses the original spectral parameter and quasi-particle type). We choose an overall time $t_{\rm micro}>0$. A quasi-particle either is in a “ghost” state ($\Delta v>0$), or not ($\Delta v=0$). We note the time $t_{0}$ of the start of a ghost state. Upon a collision of $A$ with $B$, we do $\Delta v:=\Delta v+\Delta t/{\tt distance}(A,B)$, so that within a time $\Delta t$, the correct distance is added. When entering a ghost state (at a collision), we set $\Delta t=t_{\rm micro}$, while for further collisions during the ghost state, we choose $\Delta t=t_{\rm micro}-t+t_{0}$, where $t$ is the time. The particle reverts to its normal state when $t-t_{0}=t_{\rm micro}$.

II The Lieb-Liniger model in the attractive regime

In first-quantized form, the Lieb-Liniger (LL) model is described by the Hamiltonian

[TABLE]

With repulsive interaction, $c>0$, the spectrum of Bethe excitations is composed of a single particles, and the differential scattering phase given by (5). The case of attractive interaction $c=-2\bar{c}<0$ is however not simply obtained by replacing $c$ by $-2\bar{c}$ in (5). Instead, the spectrum of Bethe excitations is more complicated: it is composed of an infinity of Bethe strings, one for every positive length $j=1,2,3,\ldots$. The attractive LL model is therefore classically implemented as a gas with infinitely-many species, representing the infinitely-many string lengths, all interacting with each other via a string-length- and velocity-dependent shift.

A thermodynamic Bethe ansatz analysis gives the following (we take $m=1/2$ for simplicity). The energy and momentum functions, for a string of length $a$ with “pseudo-momentum” $\lambda$, are given by

[TABLE]

The velocity is therefore

[TABLE]

which allows us to parametrize energy and momentum functions in terms of velocities,

[TABLE]

The differential scattering phase, expressed in the pseudo-momentum coordinates $\lambda$, is the following function:

[TABLE]

where

[TABLE]

The above can be directly applied to the formulae presented in the main text, by taking pseudo-momenta $\lambda$ as “rapidities” (this is a good definition as the differential scattering phase depends on differences of pseudo-momenta), and identifying gas species $a$ with string lengths $j$.

III Monotonicity of the effective velocity

Here we provide a demonstration that, under certain conditions on the differential scattering phase, the effective velocity is monotonic with $\theta$ (this is a fact that is used in the derivation in the main text).

Consider the case of a single particle in the spectrum:

[TABLE]

where $\tilde{\varphi}(\theta,\alpha)=\frac{\varphi(\theta,\alpha)}{p^{\prime}(\theta)}$. Assume first that $\tilde{\varphi}(\theta,\alpha)$ is positive (and also that $\rho_{\rm p}(\theta)$ is positive) and that its derivative has the following sign (with ${}^{\prime}={\rm d}/{\rm d}\theta$):

[TABLE]

Assume also that $\theta$ parametrizes the group velocity in a monotonic fashion, $(v^{\rm gr})^{\prime}(\theta)>0$. Taking a derivative, we find

[TABLE]

and thus

[TABLE]

Thanks to (26), the right-hand side of the latter equation is positive if $(v^{\rm eff})^{\prime}(\theta)\geq 0$ for all $\theta$, and thus this latter condition is consistent with this integral equation. That is, if we assume that the solution is obtainable by recursion, starting with $(v^{\rm eff})^{\prime}(\theta)=(v^{\rm gr})^{\prime}(\theta)$, then the solution is monotonic. With decaying asymptotic for the particle density $\rho_{\rm p}(\alpha)$ it is easy to see that $v^{\rm eff}(\theta)\sim v^{\rm gr}(\theta)$ as $\theta\to\pm\infty$, and thus $v^{\rm eff}(\theta)$ covers full range of velocities.

All conditions, including (26), are satisfied in the Lieb-Liniger model and the sinh-Gordon model.

The condition on the sign of $\tilde{\varphi}(\theta,\alpha)$ may be relaxed, as long as

[TABLE]

for all $\theta$.

IV Additional numerical checks

In addition to the density, we have also verified numerically the accuracy of the current simulated by the classical flea gas as compared to the exact self-similar solution, in the Lieb-Liniger model. See Fig. 3. Again, the agreement is excellent, confirming the adequacy of the classical gas for representing GHD.