A geometric viewpoint on generalized hydrodynamics

TL;DR

This paper presents a geometric reformulation of generalized hydrodynamics (GHD), revealing a new perspective that simplifies solving the equations and enhances understanding of many-body integrable systems.

Contribution

It recasts GHD equations into a geometric framework with state-dependent metrics, providing a novel solution method for initial value problems in integrable systems.

Findings

Derived a geometric interpretation of GHD equations.

Developed a new integral equation-based solution algorithm.

Demonstrated the efficiency of the geometric approach for solving GHD.

Abstract

Generalized hydrodynamics (GHD) is a large-scale theory for the dynamics of many-body integrable systems. It consists of an infinite set of conservation laws for quasi-particles traveling with effective ("dressed") velocities that depend on the local state. We show that these equations can be recast into a geometric dynamical problem. They are conservation equations with state-independent quasi-particle velocities, in a space equipped with a family of metrics, parametrized by the quasi-particles' type and speed, that depend on the local state. In the classical hard rod or soliton gas picture, these metrics measure the free length of space as perceived by quasi-particles, in the quantum picture, they weigh space with the density of states available to them. Using this geometric construction, we find a general solution to the initial value problem of GHD, in terms of a set of integral equations where time appears explicitly. These integral equations are solvable by iteration and provide an extremely efficient solution algorithm for GHD.

Equations

$$2\pi {\rho }_{\mathrm{s}}(x;\theta )=p^{\prime }(\theta )+\int \mathrm{d}\alpha \phi (\theta -\alpha ){\rho }_{\mathrm{p}}(x;\alpha ).$$

$$n(x;\theta )=\frac{{\rho }_{\mathrm{p}}(x;\theta )}{{\rho }_{\mathrm{s}}(x;\theta )}.$$

$$2\pi {\rho }_{\mathrm{p}}(x;\theta )$$

$$2\pi {\rho }_{\mathrm{s}}(x;\theta )$$

$$h_{[n]}^{\mathrm{dr}}(\theta )=h(\theta )+\int \frac{\mathrm{d}\alpha }{2\pi }\phi (\theta ,\alpha )n(\alpha )h_{[n]}^{\mathrm{dr}}(\alpha ).$$

$$v_{[n]}^{\mathrm{eff}}(\theta )=\frac{(E^{\prime }{)}_{[n]}^{\mathrm{dr}}(\theta )}{(p^{\prime }{)}_{[n]}^{\mathrm{dr}}(\theta )}.$$

$$\partial_{t}\rho_{\rm p}(x,t;\theta)+\partial_{x}\big{(}v^{\rm eff}_{[n(x,t)]}(\theta)\rho_{\rm p}(x,t;\theta)\big{)}=0.$$

$${\partial }_{t}n(x,t;\theta )+v_{[n(x,t)]}^{\mathrm{eff}}(\theta ){\partial }_{x}n(x,t;\theta )=0.$$

$$\partial_{t}\rho_{\rm s}(x,t;\theta)+\partial_{x}\big{(}v^{\rm eff}_{[n(x,t)]}(\theta)\rho_{\rm s}(x,t;\theta)\big{)}=0.$$

$${\rho }_{\mathrm{p}}(x,t;\theta )={\rho }_{\mathrm{p}}^{-}(\theta )\text{mbox}forallx\le x_0,t\in [0,T],\theta \in \mathbb{R}.$$

$$\mathrm{d}V={\rho }_{\mathrm{s}}(x;\theta )\mathrm{d}x2\pi \mathrm{d}\theta .$$

$$\mathrm{d}q={\mathcal{K}}_{[n(x)]}(\theta )\mathrm{d}x.$$

$${\mathcal{K}}_{[n(x)]}(\theta )=\frac{2\pi }{p^{\prime }(\theta )}{\rho }_{\mathrm{s}}(x;\theta )=\frac{(p^{\prime }{)}_{[n(x)]}^{\mathrm{dr}}(\theta )}{p^{\prime }(\theta )}.$$

$$q=q(x;\theta )=\frac{2\pi }{p^{\prime }(\theta )}{\int }_{x_0}^x\mathrm{d}y{\rho }_{\mathrm{s}}(y;\theta ).$$

$$n(x;\theta )\mathrm{d}V={\rho }_{\mathrm{p}}(x;\theta )\mathrm{d}x2\pi \mathrm{d}\theta$$

$$n(x;\theta )=\frac{{\rho }_{\mathrm{p}}(x;\theta )}{{\rho }_{\mathrm{s}}(x;\theta )}.$$

$$q(x,t;\theta )=\frac{2\pi }{p^{\prime }(\theta )}{\int }_{x_0}^x\mathrm{d}y{\rho }_{\mathrm{s}}(y,t;\theta ).$$

$$v^{-}(\theta )=v_{[n^{-}]}^{\mathrm{eff}}(\theta ){\mathcal{K}}_{[n^{-}]}(\theta )=\frac{(E^{\prime }{)}_{[n^{-}]}^{\mathrm{dr}}(\theta )}{p^{\prime }(\theta )}.$$

$${\partial }_t\hat{n}(q,t;\theta )+v^{-}(\theta ){\partial }_q\hat{n}(q,t;\theta )=0$$

$$\hat{n}(q(x,t;\theta ),t;\theta )=n(x,t;\theta ).$$

$${\partial }_{t}n(x,t;\theta )={\partial }_t\hat{n}(q,t;\theta )+{\partial }_{t}q(x,t;\theta ){\partial }_q\hat{n}(q,t;\theta ).$$

$${\partial }_{t}n(x,t;\theta )+\frac{p^{\prime }(\theta )({\hat{v}}^{-}(\theta )-{\partial }_{t}q(x,t;\theta ))}{2\pi {\rho }_{\mathrm{s}}(x,t;\theta )}{\partial }_{x}n(x,t;\theta )=0.$$

$${\partial }_{t}q(x,t;\theta )$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{d}V$$

$$\displaystyle=\bigl{[}\partial_{t}\mathcal{K}_{[n(x,t)]}(\theta)+v^{\rm eff}_{[n(x,t)]}(\theta)\partial_{x}\mathcal{K}(\theta)\bigr{]}{\rm d}x_{t}\wedge{\rm d}p_{t}+\mathcal{K}_{[n(x,t)]}(\theta)\partial_{x}v^{\rm eff}_{[n(x,t)]}{\rm d}x_{t}\wedge{\rm d}p_{t},$$

$$\displaystyle=\bigl{[}\partial_{t}\mathcal{K}_{[n(x,t)]}(\theta)+\partial_{x}(v^{\rm eff}_{[n(x,t)]}(\theta)\mathcal{K}_{[n(x,t)]}(\theta))\bigr{]}{\rm d}x_{t}\wedge{\rm d}p_{t},$$

$$\left(\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{d}{x}_t\right)\wedge \mathrm{d}{p}_t=\mathrm{d}{\dot{x}}_t\wedge \mathrm{d}{p}_t=({\partial }_x{v}_{[n(x,t)]}^{\mathrm{eff}}(\theta )\mathrm{d}{x}_t+{\partial }_p{v}_{[n(x,t)]}^{\mathrm{eff}}(\theta )\mathrm{d}{p}_t)\wedge \mathrm{d}{p}_t={\partial }_x{v}_{[n(x,t)]}^{\mathrm{eff}}\mathrm{d}{x}_t\wedge \mathrm{d}{p}_t$$

$$n(x,t;\theta )=n^0(u(x,t;\theta );\theta )$$

$${\int }_{x_0}^{u(x,t;\theta )}\mathrm{d}y{\mathcal{K}}_{[n^0(y)]}(\theta )+v^{-}(\theta )t={\int }_{x_0}^x\mathrm{d}y{\mathcal{K}}_{[n(y,t)]}(\theta ).$$

$$\hat{n}(q,t;\theta )={\hat{n}}^0(q-v^{-}(\theta )t;\theta )$$

$$q(u(x,t;\theta ),0;\theta )=q(x,t;\theta )-v^{-}(\theta )t.$$

$${\mathcal{K}}_{[n^0(u(x,t;\theta ))]}(\theta ){\partial }_{t}u(x,t;\theta )+{\mathcal{K}}_{[n(x,t)]}(\theta )v_{[n(x,t)]}^{\mathrm{eff}}(\theta ){\partial }_{x}u(x,t;\theta )=0$$

Full text

A geometric viewpoint on generalized hydrodynamics

Benjamin Doyon∗, Herbert Spohn*†* and Takato Yoshimura∗

∗* Department of Mathematics, King’s College London, Strand, London WC2R 2LS, U.K.

† Physik Department and Zentrum Mathematik, Technische Universität München, Boltzmannstrasse 3, 85748 Garching, Germany*

Generalized hydrodynamics (GHD) is a large-scale theory for the dynamics of many-body integrable systems. It consists of an infinite set of conservation laws for quasi-particles traveling with effective (“dressed”) velocities that depend on the local state. We show that these equations can be recast into a geometric dynamical problem. They are conservation equations with state-independent quasi-particle velocities, in a space equipped with a family of metrics, parametrized by the quasi-particles’ type and speed, that depend on the local state. In the classical hard rod or soliton gas picture, these metrics measure the free length of space as perceived by quasi-particles; in the quantum picture, they weigh space with the density of states available to them. Using this geometric construction, we find a general solution to the initial value problem of GHD, in terms of a set of integral equations where time appears explicitly. These integral equations are solvable by iteration and provide an extremely efficient solution algorithm for GHD.

1 Introduction

Generalized hydrodynamics (GHD) [1, 2] is a hydrodynamic theory where the notion of local equilibration to a Galilean (or relativistic) boost of a Gibbs state, is replaced by that of local relaxation to generalized Gibbs ensembles [3, 4, 5]. It is expected to emerge in appropriate hydrodynamic limits in both quantum and classical integrable many-body systems, including field theory and spin chains, and describes time evolution in inhomogeneous backgrounds or from inhomogeneous states. It is in the context of steady states arising from domain-wall initial conditions (see, for instance, a review [6]) that it was originally introduced, where it solved the long-standing problem of obtaining full density and current profiles in interacting integrable quantum systems [1, 2]. It was generalized to include inhomogeneous force fields [7]. It is seen as emerging in integrable classical systems such as the hard rod fluid [8, 9] and soliton gases [10, 11, 12, 13, 14], and a classical molecular dynamics solver has been developed for the general form of GHD [14]. It was applied to study spin transport [15, 16, 17], transport in the hard rod fluid [18], and quantum dynamics of density profiles such as propagating waves in interacting Bose gases [19, 20]. Most of these studies concentrate on the emerging hydrodynamics at the Euler scale, at which GHD was originally formulated, but see [7, 22, 18] for discussions of viscosity effects. Here we do not discuss such effects.

The goal of this paper is twofold. First, we provide a geometric interpretation of GHD. We show that it is a theory for a gas of freely (inertially) propagating particles, but within a space whose metric depends both on the type and velocity of the particle, and on local distribution of particles in the gas. In the hard rod fluid, the metric has a clear interpretation: it measures the free space available between the rods, a notion that was used in [18] in a derivation of the exact solution to the domain wall initial problem. The observation that this generalizes to soliton gases, still described by GHD, then suggests the metric construction proposed here. It is worth noting some similarity, in spirit, to Einstein’s theory of general relativity, where currents are conserved in a metric that is determined by the matter content. Second, we use this geometric construction in order to provide a system of integral equations that solve the initial-value problem of GHD in full generality. The integral equations involve the initial condition and the time parameter in an explicit fashion, essentially integrating out the time direction. They can be solved by iteration. We confirm their validity by providing comparisons with direct solutions of the GHD partial differential equations. It is surprising that a general solution to a hydrodynamic equation can be obtained, and this might connect with the integrability of the GHD equations themselves, as found in soliton gases [10, 11, 12, 13, 21].

The paper is organized as follows. In Section 2 we review some of the main features of GHD. In Section 3 we develop the geometric interpretation of GHD in generality. In Section 4, we derive the integral equations that solve the initial value problem. Finally, in Section 5 we conclude, and in Appendix A we briefly explain the specialization to the hard rod problem.

2 Overview of generalized hydrodynamics

The most powerful formulation of GHD to date [1, 2] uses the physical notion of quasi-particles (although other formulations will doubtless come to the fore in the future). Integrable models solvable by Bethe ansatz, and integrable classical gases such as the hard rod model [8, 9] and soliton gases [14], can be seen as models for interacting quasi-particles. A quasi-particle is specified by a spectral parameter $\theta\in{\mathbb{R}}$, parametrizing its momentum $p(\theta)$ (it can be taken as the velocity in the Galilean case, or the rapidity in the relativistic case; in general it is just a parametrization of the momentum). The scattering kernel $S(\theta,\alpha)$ characterizes the interaction amongst quasi-particles, and it is customary to define the differential scattering phase $\varphi(\theta,\alpha)=-{\rm i}\,{\rm d}\log S(\theta,\alpha)/{\rm d}\theta$, which we assume to be symmetric. In this paper we assume for simplicity that there is a single particle species, but all equations are directly generalizable to models with many species such as the Heisenberg chain, simply by viewing the rapidity $\theta$ as a multi-index ${\bm{\theta}}=(\theta,j)$.

In the quasi-particle formulation of GHD, the local fluid state is described by a function specifying the distribution of quasi-particles. This function can be taken as the density $\rho_{\rm p}(x;\theta)$: the infinitesimal $\rho_{\rm p}(x;\theta)\,{\rm d}x\,{\rm d}\theta$ is the number of quasi-particles in the phase space volume element $[x,x+{\rm d}x]\times[p(\theta),p(\theta)+p^{\prime}(\theta){\rm d}\theta]$ (here and below the prime symbol (′) denotes a derivative with respect to the spectral parameter).

It is conventional to define the state density $\rho_{\rm s}(x;\theta)$ via

[TABLE]

In terms of these densities, one defines the occupation function

[TABLE]

The occupation function can be taken, instead of the quasi-particle density, as a characterization of the local fluid state. One may go from occupation function to quasi-particle and state densities as follows:

[TABLE]

where the dressing operation depends on the function $n(x):\theta\mapsto n(x;\theta)$ (seen as an $x$-dependent function of the spectral parameter), and is defined, for any function $n$ (of the spectral parameter), by the solution to the following linear integral equation:

[TABLE]

As in any hydrodynamic theory, a model of generalized hydrodynamics also necessitates the equations of state: relations between the conserved currents and the conserved densities that emerge due to the interactions and dynamics of the constituents. This is determined by the energy function $E(\theta)$. More precisely, it was found in [1, 2] that, in the quasi-particle picture of GHD, for each $\theta$, the function $\rho_{\rm p}(x;\theta)$ is a conserved density with associated current $v^{\rm eff}_{[n(x)]}(\theta)\rho_{\rm p}(x;\theta)$, where the effective velocity is given by [1, 2, 23]

[TABLE]

Equation (5) can be seen as the equation of state.

The Euler-type hydrodynamic equations of GHD, in the quasi-particle language, are therefore

[TABLE]

These hold at the Euler scale, and thus omit viscosity terms; see [18] for explicit viscosity terms in the hard rod fluid, and [7] for an analysis of viscosity in integrable quantum models. Surprisingly, it turns out that the occupation function provides the normal modes of GHD, and thus is convectively conserved: Eq. (6) is equivalent to [1, 2]

[TABLE]

Note that these imply the continuity equation for the state density,

[TABLE]

The generalization to include force terms has also been obtained, see [7]. In quantum field theory models, Eqs (6) and (7) are seen [1] to emerge from the conservation of local densities $\partial_{t}\langle q_{i}(x,t)\rangle+\partial_{x}\langle j_{i}(x,t)\rangle=0$, in the hydrodynamic approximation where averages are evaluated within entropy-maximized local states. A different argument [2] based on a kinetic theory leads to the same equations in quantum chains.

The above formulation is valid in complete generality. However, there are conditions on $\varphi(\theta,\alpha)$, $p(\theta)$, $E(\theta)$ and $\rho_{\rm p}(\theta)$ for it to provide physically sensible results. We do not know exactly what these conditions are, but the derivation presented below makes sense if the momentum derivative and the state density are strictly positive $p^{\prime}(\theta)>0$, $\rho_{\rm s}(\theta)>0$.

3 Geometry of GHD

Eqs (5) and (7) can be viewed as arising from a gas of colliding classical particles [14], generalizing the hard rod problem studied a long time ago [9]. Let us recall the main features of these models.

In the hard rod fluid, segments, all of fixed length $a$, move on the line, freely (inertially) except for collisions at which they exchange their velocities. By definition, a quasi-particle is a tracer of a given velocity: it follows the trajectory of a given velocity. The quasi-particle therefore travels like a free particle, except at collisions. Two quasi-particles, with velocities $v,w$, collide if their actual positions, $x_{w},x_{v}$ and velocities satisfy either $x_{v}-x_{w}=a$ and $v<w$, or $x_{w}-x_{v}=a$ and $w<v$. Because of the collision, in the first case $x_{v}$ jumps to $x_{v}-a$ and $x_{w}$ to $x_{w}+a$, while in the second case $x_{v}$ jumps to $x_{v}+a$ and $x_{w}$ to $x_{w}-a$. The size of the jump is independent of $v,w$.

As a natural generalization, the jump size may be assumed to depend on the velocities of both collision partners. In other words, in the rules above, $-a$ is replaced by the general kernel $\varphi(v,w)$ without any particular sign restrictions. That is, $\varphi(v,w)$ is seen as an effective rod length, as perceived by the quasi-particles at velocities $v$ and $w$ when they collide. In order to take into account correctly cases where many particles interact within short time intervals, care must be taken in implementing the jumps. However, the above rules provide the sufficient features for our analysis, and for connecting with soliton scattering. See [14] for more details. Heuristically, $\varphi(v,w)<0$ corresponds to a repulsive and $\varphi(v,w)>0$ to an attractive interaction. In our context, velocity is replaced by the spectral parameter. Thus in a collision two quasi-particles at spectral parameters $\theta$ and $\alpha$ jump by the distance $-\varphi(\theta,\alpha)$ with a sign governed by the rules from above.

There is a natural geometric re-interpretation of the Euler equations for the hard rod or soliton gas. The idea is that the quasi-particles travel as free particles, which do not seem to interact, if one “shrinks” the effective rods’ lengths to zero. Indeed the motion of a test quasi-particle in the free space available in-between the (effective) rods in the gas in which it moves, is that of a free point-like particle at the group velocity $v^{\rm gr}(\theta)=E^{\prime}(\theta)/p^{\prime}(\theta)$. This free space is, however, state-dependent, as it depends on the density of rods in the fluid, and in the general case it also depends on the spectral parameters (the type and speed) of the test quasi-particle itself. Therefore, we may put a family of metrics on the one-dimensional space in which quasi-particles propagate, parametrized by the spectral parameters, whose associated length is the available length as perceived by the test quasi-particles. Each such metric is determined by the local state, and in this metric, the test quasi-particle propagates freely. The relation that fixes the metric as a function of the local state may be interpreted as an “Einstein equation”, and the free propagation is a conservation equation within the metric. The full system is therefore, in spirit, somewhat similar to Einstein’s theory of general relativity.

3.1 Metric and continuity equation

In this subsection, we construct a dynamics on quasi-particles as per the above geometric arguments, and show that it gives rise to GHD.

Let us assume that there is some point $x_{0}$ such that the densities are independent of time $t$ for all $x\leq x_{0},\;t\in[0,T]$. They are therefore also independent of $x$ in this region by the continuity equation:

[TABLE]

Then the derivation below will be valid for all times up to $T$. The choice of the left side is purely conventional. We will use the symbol $n^{-}(\theta)$ for the associated occupation function, and $\rho_{\rm s}^{-}(\theta)$ for the state density.

This assumption is immediately satisfied if all particles are initially located in a bounded region of phase space, as well as for domain-wall initial conditions [6, 1, 2], with baths that extend both on the left and the right (in both cases the assumption holds for every finite $T$ by choosing $|x_{0}|$ large enough).

Let us consider the phase space $(x;\theta)\in{\mathbb{R}}\times{\mathbb{R}}$. We put on this space the volume element ${\rm d}V$ equal to the “available”, or “free”, volume: the space available within ${\rm d}x$ as seen by a test quasi-particle at rapidity $\theta$, times $2\pi{\rm d}\theta$:

[TABLE]

Convenient coordinates are the $(q,p)$ coordinates, where $p$ is the momentum and the volume element is ${\rm d}V={\rm d}q\,{\rm d}p$. This volume element therefore gives a relation between ${\rm d}q$ and ${\rm d}x$, inducing a one-dimensional, $\theta$- and state-dependent metric on $x$-space,

[TABLE]

One takes ${\rm d}q^{2}=g(x;\theta)\,{\rm d}x^{2}$ as the infinitesimal square length at constant $\theta$ with metric $g(x;\theta)=\mathcal{K}_{[n(x)]}(\theta)^{2}$, and

[TABLE]

The quantity ${\rm d}q$ may be interpreted as the infinitesimal number of states available within the interval $[x,x+{\rm d}x]$ per unit momentum, as seen by a test quasi-particle at spectral parameter $\theta$. In this point of view, there is a one-parameter family of metrics, in bijection with the one-parameter family of conserved quasi-particle numbers $n(\theta)$. Each metric relates to the one-dimensional space on which the quasi-particle at spectral parameter $\theta$ moves. Clearly, these metrics depend on the local state. The coordinate $q$ is defined by integrating from $x_{0}$:

[TABLE]

Consider the free density of particles $n(x;\theta)$, the density per unit free volume $V$. The total number of particles within ${\rm d}V$ is

[TABLE]

and therefore

[TABLE]

We thus see that the occupation function (7) is equal to the particle density per unit free volume in the geometric language. We define $\hat{n}(q;\theta)$ by $\hat{n}(q(x;\theta);\theta)=n(x;\theta)$.

We now put a dynamics on this system, evolving in time $t$. All quantities acquire a $t$ dependence, including the relation between the $q$ coordinate and the $x$ coordinate,

[TABLE]

The quasi-particles are deemed to propagate ballistically in the free space. From the above construction, it seems natural that the velocity of propagation depend on the asymptotic particle density $\rho_{\rm p}^{-}(\theta)$. We choose the following definition of the propagation velocity:

[TABLE]

This has the interpretation as the asymptotic effective velocity as expressed with respect to $q$-space, thus multiplied by the square-root of the metric. The fluid variable $\hat{n}(q,t;\theta)$ therefore satisfies the equation for free particle motions at the velocity $v^{-}(\theta)$ as function of the position in the free space $q$ and of time $t$:

[TABLE]

where $\hat{n}(q,t;\theta)$ is defined by

[TABLE]

We now show that (18) is equivalent to (7).

Proof.

The proof we present goes in one direction of the equivalence, the other direction being immediate. From (19), we have

[TABLE]

Therefore, using (18), as well as (11) and (12) in order to relate $\partial_{q}$ to $\partial_{x}$, we obtain

[TABLE]

We evaluate $\partial_{t}q(x,t;\theta)$ by using (8):

[TABLE]

from which (7) follows.

If the asymptotic particle density is zero, then $v^{-}(\theta)=v^{\rm gr}(\theta)$ is the group velocity: in this case, the particles in $q$ space propagate as free particles at the group velocity. The derivation above is unaffected by the change $\theta\mapsto(\theta,j$) where $j$ is an internal index for particle types; the only difference, in the derivation, is that integrals over spectral parameters are augmented by sums over particle types. Therefore the result stays the same in models with many particle species. The set of “Einstein-Euler” equations (16) and (18) can be interpreted as an infinite set of Euler-type conservation laws, one for each quasi-particle spectral parameter, in one-dimensional spaces whose metrics, one for each spectral parameter, are determined by the local state.

3.2 GHD and invariance of volume form

The choice of the metric (12) is also natural if we consider GHD as a non-Hamiltonian system [24], which is a classical system whose determinant of the Jacobian associated with the coordinate transformation $(x_{0},p_{0})\mapsto(x_{t},p_{t})$ is not unity: the standard volume element ${\rm d}x\,{\rm d}p$ is not preserved in the course of time-evolution. Volume preservation is recovered under an appropriate choice of metric. Here, we show that preservation of the volume element ${\rm d}q\wedge{\rm d}p=\mathcal{K}_{[n(x,t)]}(\theta)\,{\rm d}x\wedge{\rm d}p$ along the path $(\dot{x}_{t},\dot{p}_{t})=(v^{\rm eff}_{[n(x,t)]}(\theta),0)$ of a test particle within GHD, amounts to the continuity equation (8) for $\rho_{\rm s}(x,t;\theta)$

Let us denote the standard volume form as a differential 2-form ${\rm d}x_{t}\wedge{\rm d}p_{t}$. Then,

[TABLE]

where in the second line we used

[TABLE]

and ${\rm d}x_{t}\wedge\left(\frac{{\rm d}}{{\rm d}t}{\rm d}p_{t}\right)=0$. Thus $\frac{{\rm d}}{{\rm d}t}{\rm d}V_{t}=0$ implies (8).

Likewise, the invariance of the density element $n(x,t;\theta){\rm d}V_{t}$ under the time-evolution leads to the continuity equation (6). This is a natural generalization of the Liouville theorem to the non-Hamiltonian systems. It is also a simple matter to extend this derivation to the case where a system is subject to an external force.

4 Integral equations solving the initial value problem

We have shown that the set of equations (16), (18) and (19) is equivalent to (7). Therefore, we may look for solving these simpler-looking equations. This turns out to be possible, providing a scheme that is a very efficient numerical method for solving the GHD initial value problem for a time evolution without inhomogeneous external force field.

4.1 Solution to the initial value problem

Let $n^{0}(x,\theta)$ be the initial occupation function, with associated state density $\rho_{\rm s}^{0}(x,\theta)$. Recall the metric (12). Clearly, if there is no interaction, $\varphi(\theta,\alpha)=0$, then $\mathcal{K}_{[n]}(\theta)=1$. We show that the solution to the initial value problem for (7) is given by

[TABLE]

where $u(x,t;\theta)$ is determined implicitly by the integral equation

[TABLE]

These simply indicate that the particles move at the velocity $v^{-}(\theta)$ in the free space. Note that $u(x,t;\theta)$ is monotonically increasing as a function of $x$ because $\mathcal{K}_{[n]}(\theta)$ is positive. Equations (25) and (26) are derived as follows. Solving (18) is simple:

[TABLE]

where $\hat{n}^{0}(q(u,0;\theta),\theta)=n^{0}(u;\theta)$ is the initial condition in the $q$-space. Therefore $n(x,t;\theta)=\hat{n}^{0}(q(x,t;\theta)-v^{-}(\theta)t;\theta)$. Hence we must define $u=u(x,t;\theta)$ as solving

[TABLE]

This immediately yields (26).

The explicit time parameter appears in (26). Given the initial fluid states at every point $y$ we may construct $\mathcal{K}_{[n_{0}(y)]}$. With this, we then must solve simultaneously the system of integral equations (25) and (26) for the unknowns $n(x,t;\theta)$ and $u(x,t;\theta)$ as functions of $(x;\theta)$ for $t$ fixed. Here (25) involves $n(x,t;\theta)$ and $u(x,t;\theta)$, and (26) involves $u(x,t;\theta)$, and involves $n(y,t;\alpha)$ for all $y\in[x_{0},x]$ and all $\alpha\in{\mathbb{R}}$. In the case without interaction, we immediately find the expected relation $u(x,t;\theta)=x-v^{-}(\theta)t$ with $v^{-}(\theta)=v^{\rm gr}(\theta)$. In general, $u(x,t;\theta)$ satisfies the differential equation

[TABLE]

and, using (25), this can be seen as a generalization to GHD of the solution of a single-component fluid by the method of characteristics. In particular, in a uniform, stationary state where $n(x,t;\theta)=n(\theta)$ is independent of $x$ and $t$, we have $u(x,t;\theta)=x-v^{\rm eff}_{[n]}(\theta)\,t$, representing the propagation of a particle at rapidity $\theta$ within a state characterized by $n$.

In the case of the domain-wall initial condition, this general solution specializes to that provided in [1, 2, 18]. Indeed, assume that $n^{0}(x;\theta)=n^{-}(\theta)$ for $x<0$, and $n^{0}(x;\theta)=n^{+}(\theta)$ for $x>0$. Since $u(x,t;\theta)$ is monotonic with $x$, we must find $x_{*}=x_{*}(t;\theta)$ solving $u(x_{*},t;\theta)=0$; we have $n(x,t;\theta)=n^{-}(\theta)\chi(x<x_{*}(t;\theta))+n^{+}(\theta)\chi(x>x_{*}(t;\theta))$ where $\chi$ is the indicator function. By (28), differentiating with respect to $t$ and taking into account (16), this implies

[TABLE]

From the result (22) this gives

[TABLE]

With scale-invariant initial condition, it is natural to assume that the solution to the scale-invariant equations (7) is self-similar (that is, depends on the ray $x/t$ only). If we set $x_{*}(t;\theta)=t\xi_{*}(\theta)$, then the solution presented here is self-similar, and $n(x_{*}(t;\theta),t)=n(\xi_{*}(\theta),1)$. In this case we find

[TABLE]

which is indeed the solution given in [1, 2, 18].

We do not know yet how to address the uniqueness of the solution to (25), (26), and in particular how to show that self-similarity must hold in the domain-wall problem. However, we note that, in the hard rod problem, the idea of “shrinking” the rods to points and then freely evolving them in time, which directly corresponds to the integral equations (25), (26), was used in [9] in order to show uniqueness properties. Thus the above may give in the future some insight into this problem.

Finally, again, it is clear that all the above results hold under $\theta\mapsto(\theta,j)$ when many particle types are involved.

4.2 Numerical method

The set of equations (25), (26) are of interest from a theoretical perspective. But at the same time they furnish an efficient algorithm for numerically solving the initial value problem of GHD.

In general, Equations (25), (26) can be solved by the following iteration scheme. First, one sets, as an initial condition to the iteration scheme, $n(x,t;\theta)=n_{0}(x;\theta)$. Note that with the knowledge of $n(x,t;\theta)$ and $n_{0}(x;\theta)$, both integrands in (26) can be evaluated. One then solves (26) for $u(x,t;\theta)$, and constructs a new iteration of $n(x,t;\theta)$ using (25). The process is then repeated until $n(x,t;\theta)$ converges. Alternatively, one may set, as an initial condition to the iteration scheme, $n(x,t;\theta)$ to the right-hand side of (25) with $u(x,t;\theta)=x-v^{-}(\theta)t$. The latter only depends on the asymptotic velocity $v^{-}(\theta)$ (17), whose evaluation only requires the knowledge of the initial state.

In order to explicitly check that solving the integral equations by iteration coincides with directly solving the GHD equation (7), we focus on the “bump-release problem” studied for instance in [20]. In this problem, a Lieb-Liniger (LL) gas, with $\varphi(\theta,\alpha)=2c/((\theta-\alpha)^{2}+c^{2})$, is initially in the ground state or a finite-temperature state of a potential with an inverted Gaussian centered at the origin, where a density accumulation occurs. The potential is then suddenly flattened, so that the density accumulation is released into two oppositely propagating waves. At zero temperature, GHD results have been explicitly checked in [20] to agree with numerical quantum evolution, and it has been shown that GHD equations (7) simplify to conventional hydrodynamics for finite times, providing a simple alternative for the time evolution. This allows for compelling comparison with the current iterative method. We also study the bump release problem at finite temperature where conventional hydrodynamics fails, comparing with the molecular dynamics developed in [14] and thus giving further compelling evidence.

Quite surprisingly, we observe that a few iterations, of the order of 2 to 4, are enough to accurately reproduce the space-time profile of the density $\rho(x,t)=\int{\rm d}\theta\,\rho_{\rm p}(x,t;\theta)$.

In Fig. 1, we depict the time-evolution of the LL gas with a coupling constant $c=0.5$ (with mass set to 1) upon a release from the ground state of an energy potential of the form $V(x)=-0.1e^{-x^{2}/8}-\mu_{\infty}$ (coupling to the density of particles). This is a perturbation of a background chemical potential $\mu_{\infty}$ corresponding to the background density $\rho_{\infty}=0.512$. The initial fluid state is obtained by a local density approximation of the ground state. With this choice of parameters, the interaction is strong enough to render the dynamics nontrivial. The solid lines represent curves from a direct solution of conventional hydrodynamics using finite-element approximations, and dots are from the second iteration of (25) and (26). Agreement is excellent. We show results only for times before the formation of sharp structures, where conventional hydrodynamics fails (these are “dissolving shock,” and the solution beyond this time was first obtained in [20]). Naturally, the steep wave front that develops would necessitate a smaller finite element $\Delta x$ in order to improve precision.

A similar observation is made at finite temperature. In Fig. 2, we depict the time-evolution of the LL gas with $c=2$, and initial state obtained by a local density approximation of the potential $V(x)=-5\,e^{-(x/50)^{2}}-1$ at temperatures $T=1$ (a) and $T=2$ (b). We used 4 iterations, and enough finite-element divisions for high precision results. It required about two minutes of computer time on a standard laptop for every of the four curves in each graph. Using the iterative method with $n_{0}(x;\theta)$ set to the $t=0$ initial value, we did not see any degradation in precision as a function of time $t$. Comparison is made against a simulation using molecular dynamics (the flea gas algorithm) [14], with excellent agreement. There are only small differences between the cases $T=1$ and $T=2$, which seem to be well captured by the iterative method. The iterative method is much more precise than molecular dynamics, for less computer time.

Other algorithms for solving GHD are (1) the approximate scheme used in [19], essentially based on an approximation of the present algorithm that is valid for small times; (2) the direct solution to the finite hydrodynamics to which GHD reduces in the zero-entropy subspace [20]; and (3) the molecular dynamics of [14]. The latter two are the only ones valid also for evolution within external, inhomogeneous force fields. However, without such force fields, the numerical algorithm for the GHD initial value problem developed here appears to be very efficient, providing accurate solutions for the full distribution $n(x,t;\theta)$, at any given time $t$, within at most few minutes of laptop computer time.

5 Conclusion

We have shown that the equations of GHD (those valid at the Euler scale), in the quasi-particle formulation, can be recast into hydrodynamic conservation equations for a gas of particles freely propagating at velocities that are determined by the homogeneous and stationary state asymptotically far in space. This gas propagates in a space whose metric is proportional to the local state density. In the classical interpretation of GHD, it measures the effective space available for quasi-particles to travel between collisions. In the quantum interpretation, it measures distances by weighing with the number of Bethe roots available in local states, as if the actual space in which a particle propagates were proportional to the “Bethe space” 111We thank Bruno Bertini for suggesting this interpretation.. The structure is somewhat similar to that of Einstein’s equations for general relativity, where a metric is dynamical, determined by the particle content at that space-time point, and particles satisfy conservation equations within this metric. In the present case, because of the infinity of conservation laws, there are infinitely many conserved quasi-particle numbers, parametrized by the spectral parameter $\theta$, and there is a one-parameter family of metrics, as perceived by test quasi-particles at spectral parameters $\theta$.

Importantly, we have shown that this new viewpoint allows for an exact solution to the initial value problem of GHD. The exact solution is expressed as a system of nonlinear integral equations, Equations (25) and (26), which are readily solvable on standard laptop computer.

This geometric construction seems to connect naturally with the recent work [25] showing how inhomogeneous states and evolutions in free fermion models can be recast into field theory on curved space. It would be interesting to analyze the uniqueness properties of solutions to Equations (25), (26). Using these equations, the approach to stationary states could also be analyzed. It would also be interesting to generalize the above geometric construction to the presence of force fields. Finally, there is an elegant symplectic structure underlying the above geometric construction, with symplectic form ${\rm d}q\wedge{\rm d}p$, which we plan to develop in the future.

Acknowledgments. We thank B. Bertini for discussions and J. Dubail for insightful comments on the manuscript. BD and TY thank City University New York (workshop “Dynamics and hydrodynamics of certain quantum matter”, March 2017) for hospitaliy. TY is grateful for the support from the Takenaka Scholarship Foundation.

Note added: As the first version of this paper was in its last stage of preparation, the work [19] appeared which, amongst other things, develops aspects of the semi-Hamiltonian and integrability structures of GHD, see also [21]. A natural geometry appears within this context, with, surprisingly, a metric of similar form to that introduced here, and a certain formal solution by quadrature. It would be very interesting to understand the relation between these two viewpoints.

Appendix A The hard rod fluid

The geometric interpretation developed in the main text can be specialized to the hard rod fluid; it is instructive to see this explicitly. The system is Galilean with a single quasi-particle specie, and we choose unit mass. Let us use the more standard notation $v=\theta$ for the velocity, and choose $p(v)=v$ and $\varphi(v,w)=-a$, where $a$ is the length of the rods.

The particle density is $\rho(x)=\int{\rm d}v\rho_{\rm p}(x;v)$. The density of state simplifies and is independent of $v$, with $2\pi\rho_{\rm s}(x)=1-a\rho(x)$. We denote the average velocity by

[TABLE]

The effective velocity was shown in [18] to be $v^{\rm eff}_{[n]}(v)=(v-a\rho u)/(1-a\rho)$. Similarly to what we did in the main text, we assume that there is some point $x_{0}$ such that the densities on its left are all independent of space-time $t$ up to some time $T$.

The geometry is very clear in this case, as the coordinate $q$ corresponds to shrinking the rod length $a$ to zero. The length element ${\rm d}q$ is related to the original coordinate’s infinitesimal ${\rm d}x$ as

[TABLE]

with $0\leq a\rho(x)\leq 1$. Contrary to the general case, it does not depend on the velocity. The density of particle with respect to the new metric is

[TABLE]

which was referred to as the free density in [18] (the specialization of (2) to the hard rod case). In the $q$ coordinate, taking into account the asymptotic bath with density $\rho^{-}$ and average velocity $u^{-}$, particles move at velocities $v^{-}(v)=v-a\rho^{-}u^{-}$. We define

[TABLE]

and $\hat{n}(q,t;v)$ satisfying $\hat{n}(q(x,t),t;v)=n(x,t;v)$. As in the main text, one can show that the trivial evolution

[TABLE]

is equivalent to $\partial_{t}n(x,t;v)+v^{\rm eff}_{[n(x,t)]}(v)\partial_{x}n(x,t;v)=0$.

The exact solution to the initial value problem simplifies thanks to the lack of $v$ dependence of the state density. Let $n^{0}(x;v)$ be the initial free density, with associated total density $\rho^{0}(x)$. The solution is

[TABLE]

where $y=y(x,t;v)$ is determined implicitly by the equation

[TABLE]