Solvable Hydrodynamics of Quantum Integrable Systems

TL;DR

This paper develops a hydrodynamic framework for quantum integrable systems using a kinetic equation, demonstrating accurate predictions for dynamics in models like XXZ and Lieb-Liniger, relevant to ultracold atom experiments.

Contribution

It introduces the Bethe-Boltzmann kinetic equation for quantum integrable systems and validates its predictions against numerical simulations and experimental models.

Findings

Hydrodynamics accurately predicts quantum integrable system dynamics.

Bethe-Boltzmann equation matches RDMG simulations.

Models replicate ultracold atom experiments.

Abstract

The conventional theory of hydrodynamics describes the evolution in time of chaotic many-particle systems from local to global equilibrium. In a quantum integrable system, local equilibrium is characterized by a local generalized Gibbs ensemble or equivalently a local distribution of pseudo-momenta. We study time evolution from local equilibria in such models by solving a certain kinetic equation, the "Bethe-Boltzmann" equation satisfied by the local pseudo-momentum density. Explicit comparison with density matrix renormalization group time evolution of a thermal expansion in the XXZ model shows that hydrodynamical predictions from smooth initial conditions can be remarkably accurate, even for small system sizes. Solutions are also obtained in the Lieb-Liniger model for free expansion into vacuum and collisions between clouds of particles, which model experiments on ultracold one-dimensional Bose gases.

Equations

$$H={\int }_0^{L}dx{\Psi }^{†}\left(-\frac{{\hslash }^2}{2m}{\nabla }^2-\mu \right)\Psi +c{\Psi }^{†}{\Psi }^{†}\Psi \Psi ,$$

$$Q_n(x,t)={\int }_{-\infty }^{\infty }dk{\rho }_k(x,t)q_n(k).$$

$${\partial }_t{Q}_n(x,t)+{\partial }_x{J}_n(x,t)=0.$$

$$J_n(x,t)={\int }_{-\infty }^{\infty }dk{\rho }_k(x,t)q_n(k)v_k[\rho (x,t)],$$

$${\partial }_t{\rho }_k(x,t)+{\partial }_x({\rho }_k(x,t)v_k[\rho ])=0.$$

$${\partial }_t{\vartheta }_k(x,t)+v_k[\hat{\vartheta }]{\partial }_x{\vartheta }_k(x,t)=0.$$

$${\vartheta }_k(x,t)={\vartheta }_k(x-v_k[\hat{\vartheta }(x,t)]dt,t-dt).$$

$$H=J\sum _j{S}_j^x{S}_{j+1}^x+S_j^y{S}_{j+1}^y+\Delta S_j^z{S}_{j+1}^z,$$

$$\beta (x)={\beta }_0-({\beta }_0-{\beta }_M){\mathrm{e}}^{-x^2/L^2},$$

$$Q^{\prime }(k)+{\int }_{-\infty }^{\infty }d{k}^{\prime }\mathcal{K}(k,k^{\prime })\vartheta (k^{\prime })Q^{\prime }(k^{\prime })=q^{\prime }(k).$$

$$\mathcal{K}(k,k^{\prime })=-\frac{c}{\pi }\frac{1}{c^2+(k-k^{\prime }{)}^2}.$$

$$\hat{K}[f](k)={\int }_{-\infty }^{\infty }d{k}^{\prime }\mathcal{K}(k,k^{\prime })f(k^{\prime }),\hat{\vartheta }[f](k)=\vartheta (k)f(k).$$

$$P^{\prime }(k)=(1+\hat{K}\hat{\vartheta }{)}^{-1}[1](k),E^{\prime }(k)=(1+\hat{K}\hat{\vartheta }{)}^{-1}[k^{\prime }](k),$$

$$v(k)=\frac{(\hat{1}+\hat{K}\hat{\vartheta }{)}^{-1}[k^{\prime }](k)}{(\hat{1}+\hat{K}\hat{\vartheta }{)}^{-1}[1](k)}.$$

$$\frac{\rho (x,t,k)}{\vartheta (x,t,k)}+{\int }_{-\infty }^{\infty }d{k}^{\prime }\mathcal{K}(k,k^{\prime })\rho (x,t,k^{\prime })=\frac{1}{2\pi },$$

$$v(x,t,k)=\frac{(\hat{1}+\hat{K}\hat{\vartheta }(x,t){)}^{-1}[k^{\prime }](k)}{(\hat{1}+\hat{K}\hat{\vartheta }(x,t){)}^{-1}[1](k)}.$$

$${\partial }_t\rho (x,t,k)+{\partial }_x(\rho (x,t,k)v(x,t,k))$$

$$\frac{2\pi \rho (x,t,k)}{1-2\pi \hat{K}[\rho (x,t,k^{\prime })](k)}$$

$$\frac{(\hat{1}+\hat{K}\hat{\vartheta }(x,t){)}^{-1}[k^{\prime }](k)}{(\hat{1}+\hat{K}\hat{\vartheta }(x,t){)}^{-1}[1](k)}$$

$${\partial }_t\vartheta (x,t,k)+v[\hat{\vartheta }(x,t)](k){\partial }_x\vartheta (x,t,k)$$

$$\frac{(\hat{1}+\hat{K}\hat{\vartheta }(x,t){)}^{-1}[k^{\prime }](k)}{(\hat{1}+\hat{K}\hat{\vartheta }(x,t){)}^{-1}[1](k)}$$

$$H=J\sum _{j=1}^{N-1}{S}_j^x{S}_{j+1}^x+S_j^y{S}_{j+1}^y+\Delta S_j^z{S}_{j+1}^z.$$

$${\lambda }_{\alpha ,a}^j={\lambda }_{\alpha }^j+i[(m_j+1-2a)\gamma +(1-v_j)\frac{\pi }{2}],a=1,2,\dots ,m_j,$$

$$a(\lambda ,n,v)$$

$$a_j(\lambda )=$$

$$T_{jk}(\lambda )=$$

$$+$$

$$\Delta Q_j^{\prime }(\lambda )+\sum _{k=1}^{N_t}{\int }_{-\infty }^{\infty }d{\lambda }^{\prime }{T}_{jk}(\lambda -{\lambda }^{\prime }){\vartheta }_k({\lambda }^{\prime }){\sigma }_k\Delta Q_k^{\prime }({\lambda }^{\prime })=Q_j^{\prime }({\lambda }^{\prime }).$$

$$P_j^{\prime }(\lambda )=2\pi (1+\hat{T}\hat{\vartheta }\hat{\sigma }{)}^{-1}[\vec{a}{]}_j(\lambda ),E_j^{\prime }(\lambda )=(1+\hat{T}\hat{\vartheta }\hat{\sigma }{)}^{-1}[-A{\vec{a}}^{\prime }{]}_j(\lambda ).$$

$$v_j(\lambda )=\frac{1}{2\pi }\frac{(\hat{\sigma }+\hat{T}\hat{\vec{\vartheta }}{)}^{-1}[-A{\vec{a}}^{\prime }{]}_j}{(\hat{\sigma }+\hat{T}\hat{\vec{\vartheta }}{)}^{-1}[\vec{a}]{)}_j}(\lambda ),$$

$${\partial }_t\vartheta (x,t,k)+v[\vartheta ](k){\partial }_x\vartheta (x,t,k)=0.$$

$$q_i=q(k_i),{\theta }_i(x,t)$$

$$\sum _j({\delta }_{ij}+A_{ij}(\theta ))Q_j^{\prime }(\theta )=q_j^{\prime }$$

Full text

Introduction.

Understanding the dynamics of interacting, many-body quantum systems far from equilibrium remains one of the most challenging problems in modern physics. In recent decades, this problem has taken on a new urgency thanks to rapid progress in the experimental construction of ultracold atomic systems. The tools available for strongly non-equilibrium dynamics with non-uniform initial conditions, even in integrable models whose equilibrium properties can be calculated exactly, have been restricted to low temperature and conformal invariance Sotiriadis and Cardy (2008); Bernard and Doyon (2012), to specific quantities Vasseur et al. (2015), or to long-time asymptotic behaviour Castro-Alvaredo et al. (2016); Bertini et al. (2016); De Luca et al. (2016); Bulchandani et al. (2017). Quantum integrable models include experimentally relevant examples like the Heisenberg antiferromagnet and the Lieb-Liniger gas in one dimension. They possess extensively many conserved quantities, which prevent them from thermalizing like generic ergodic systems and often result in ballistic transport properties.

The fact that integrable models can have unusual “generalized hydrodynamics” due to an infinite number of local conservation laws was first understood in the context of classical particle systems Boldrighini et al. (1983); Doyon and Spohn (2017). This stands in contrast to conventional hydrodynamics, which describes transport of only three conserved quantities, namely mass, momentum, and energy. The generalized hydrodynamics of quantum integrable models was developed recently in studies of the non-equilibrium steady state Sotiriadis and Cardy (2008); Bernard and Doyon (2012); Karrasch et al. (2013); Doyon (2015); De Luca et al. (2013); Bhaseen et al. (2015); Bernard and Doyon (2015); Vasseur et al. (2015); De Luca et al. (2014); Biella et al. (2016); Bernard and Doyon (2016); Vasseur and Moore (2016); Langmann et al. (2017) that is established at the junction between two infinite reservoirs Castro-Alvaredo et al. (2016); Bertini et al. (2016). An important insight is that making a local-density-type approximation for all local conserved charges implies a conservation law at the level of the local pseudo-momentum distribution. Thus in the context of integrable models, the hydrodynamic equations imply a fundamental “Bethe-Boltzmann equation”, which is an inversion of the logic familiar from conventional statistical mechanics.

The completeness of this equation for the two-reservoir steady state of the XXZ model, i.e., that the Bethe-Boltzmann equation correctly captures the physics of unusual quasilocal conservation laws Prosen (2011); Prosen and Ilievski (2013); Prosen (2014); Pereira et al. (2014); Ilievski et al. (2015a, b), was tested by comparing hydrodynamic predictions to known results for spin transport in the linear-response limit Ilievski and De Nardis (2017); Bulchandani et al. (2017). It was observed at the end of Bulchandani et al. (2017) that the ansatz for two reservoirs introduced in Ref. Bertini et al. (2016) was actually valid to first order in time for arbitrary smooth, locally equilibrated initial conditions. Hydrodynamics in the two-reservoir case is a function of only one variable (say $x/t$) because of the absence of any length scale in the initial condition, and consequently a first-order solution is sufficient; every other nontrivial initial condition yields dynamics that is a function of two variables, space and time. The present work builds on this earlier observation to develop converged solutions for finite-time hydrodynamical evolution from general smooth, locally equilibrated initial conditions.

This allows us to make novel and detailed physical predictions for finite-time dynamics in a wide range of physical systems. For example, we obtain the first practical hydrodynamic technique for the one-dimensional Bose gas Dunjko et al. (2001); Ohberg and Santos (2002); Pedri et al. (2003); Campbell et al. (2015); Bouchoule et al. (2016); De Rosi and Stringari (2016) that applies to arbitrary local GGE initial conditions and takes into account the higher conservation laws of the underlying quantum system. This allows us to obtain detailed profiles for the evolution of the Lieb-Liniger gas from collision type initial conditions, which could, in principle, be tested in the laboratory. At the same time, our approach allows for a hydrodynamic description of finite-time spin dynamics in the XXZ chain, in excellent agreement with results obtained from density matrix renormalization group (DMRG) techniques White (1992); Schollwoeck (2011); Karrasch et al. (2012); Barthel (2013); Kennes and Karrasch (2016).

Bethe-Boltzmann equation.

The Bethe-Boltzmann equation is a hydrodynamic description of quantum integrable systems Castro-Alvaredo et al. (2016); Bertini et al. (2016), which aims to capture non-equilibrium dynamics in such systems using the thermodynamic Bethe ansatz (TBA) Yang and Yang (1969); Gaudin (1971). The Bethe-Boltzmann equation takes its simplest form for the Lieb-Liniger interacting Bose gas, which was used in Castro-Alvaredo et al. (2016). We briefly summarize the physical assumptions leading to this equation below, following the presentation given in Bulchandani et al. (2017) for the XXZ chain. Thus consider a one-dimensional Bose gas with delta-function interactions, placed on a line of length $L$. This has Hamiltonian

[TABLE]

and the field operators satisfy canonical commutation relations $\left[\Psi^{\dagger}(x),\Psi(y)\right]=\delta(x-y)$. It is useful to set $\hbar=2m=1$. This system is called integrable for all values of the interaction strength $c$ because every $N$-body scattering process with $N>2$ factorizes as the product of two-body scattering processes. Integrability in this sense is reflected by the existence of infinitely many conserved charges. In a given macrostate of the Lieb-Liniger gas, with occupied density of states $\rho_{k}$ at pseudo-momentum $k$, these can be written as $Q_{n}=\int_{-\infty}^{\infty}dk\,\rho_{k}q_{n}(k)$,with $q_{n}(k)=k^{n}/n$ and $n=0,1,\ldots$. Let us now consider evolution of the Lieb-Liniger gas from locally equilibrated initial conditions; we assume that this is captured by a spatio-temporally varying pseudo-momentum distribution $\rho_{k}(x,t)$ 111This condition is equivalent to local quasi-stationarity as defined in Ref. Bertini and Fagotti (2016), modulo the subtleties discussed in Ref. Pozsgay et al. (2017). Thus the hydrodynamic approach is a priori inapplicable to initial states that do not admit such a description, such as Néel states.. This yields a spatio-temporal distribution of charge density, given by

[TABLE]

We note that in order for this expression to be defined (and indeed for the hydrodynamic approach as understood in Refs. Bertini et al. (2016); Castro-Alvaredo et al. (2016) to be consistent for the Lieb-Liniger gas), $\rho_{k}(x,t)$ must decrease more rapidly in $k$ than any power of $k$, for all $x$ and $t$ 222This property is preserved under time evolution from intial conditions with compact support in $x$-$k$ space (as are required for any numerical implementation) thanks to conservation of the mode occupancies, $\int dx\,\rho_{k}(x,t)$.. Motivated by conservation of $Q_{n}$ at the quantum mechanical level, let us postulate the local conservation law

[TABLE]

Surprisingly, the physically correct formula for $J_{n}(x,t)$ turns out to be given in terms of the quasiparticle velocity $v_{k}[\rho(x,t)]$ of collective excitations of the state with pseudo-momentum distribution $\rho_{k}(x,t)$, which is complicated but known from TBA Sup ; it has been found that Castro-Alvaredo et al. (2016); Bertini et al. (2016)

[TABLE]

in the hydrodynamic limit, which is connected to the validity of earlier conjectures for the Drude weight Zotos (1999). Substituting this expression into Eq. (3) and appealing to completeness of conserved charges in integrable models, one deduces a conservation law for the local pseudo-momentum distribution, given by

[TABLE]

We call this the Bethe-Boltzmann equation, as it has the structure of a dissipationless Boltzmann equation for the occupied pseudo-momentum density. Intuitively, the Bethe-Boltzmann equation has the meaning that “occupied quantum numbers are locally conserved”. We emphasize that for integrable systems, the generalized hydrodynamic equations (3) imply the fundamental Bethe-Boltzmann equation (5), in sharp contrast with the logic familiar from conventional statistical mechanics.

Finite time scheme.

In practice, it is useful to change variables to the local Fermi factor $\vartheta_{k}(x,t)$, defined as the ratio of occupied quantum numbers at pseudo-momentum $k$. This yields the advection form of the Bethe-Boltzmann equation Castro-Alvaredo et al. (2016); Bertini et al. (2016)

[TABLE]

Whereas this equation has so far only been used to analyze self-similar non-equilibrium steady states whose properties depend only on $x/t$, the purpose of this letter is to illustrate how it can be solved efficiently at finite time for arbitrary initial conditions. We propose a numerical solution to Eq. (6), based on a backwards implicit numerical scheme Bulchandani et al. (2017); Press (2007) which for time step $dt>0$, determines $\vartheta_{k}(x,t)$ from $\vartheta_{k}(x,t-dt)$ via the implicit equation

[TABLE]

This solves (6) up to order $\mathcal{O}(dt^{2})$. We emphasize that the velocity in the right-hand side of Eq. (7) depends non-linearly on the all of the Fermi factors at $(x,t)$, making Eq. (7) an implicit equation that can be solved by numerical iteration. Details of our implementation are deferred to the Supplemental Material Sup . Achieving convergence of numerical schemes for non-linear conservation laws in general, even in the low-dimensional setting, is known to be difficult Lev (2013). It is therefore remarkable that the above scheme, applied to an extremely high-dimensional system 333In practice, the exact dimension corresponds to the TBA discretization, typically around 4000 in our implementation., converges at all. Moreover, the scheme (7) is found to converge quickly as $dt\to 0$, so that one can obtain accurate results even for large time steps $dt$. From the solution of Eq. (6), one can readily compute physical quantities of interest (such as charge and current densities) using Eqs. (2) and (4).

We note that in general, non-linear systems of equations of conservation type (5) or advection type (6) are difficult to understand analytically, because of the possibility of shock formation from smooth initial conditions. From the viewpoint of mathematical rigour, the conservation form (5) is better defined, but existing analytical Lax (1957) methods for understanding conservation laws have little practical utility in the present high-dimensional limit. Ordinarily, one can make little analytical progress with non-linear advection equations. However, somewhat surprisingly, the advection form (6) lies in a special class of such systems which are possible to understand analytically. These are the “semi-Hamiltonian” or “rich” systems of hydrodynamic type Dubrovin and Novikov (1983); Tsarëv (1991); Serre (1989); Ferapontov (1991), and possess several interesting geometrical properties related to integrability (see Supplemental Material Sup ).

Hydrodynamics for the XXZ Spin Chain.

The Bethe-Boltzmann formalism can be extended to study non-equilibrium dynamics and transport in any integrable system or integrable quantum field theory. A particularly interesting example is provided by the spin-$1/2$ XXZ chain with Hamiltonian

[TABLE]

where predictions from hydrodynamics can be compared to density matrix renormalization group (DMRG) results White (1992); Schollwoeck (2011). Here, we set the coupling to $J=1$, and parameterize the anisotropy of the theory by $\Delta=\cos\gamma$. The Bethe-Boltzmann formalism for the gapless phase ($-1<\Delta<1$) of this model is discussed in detail elsewhere in the literature Bertini et al. (2016); Ilievski and De Nardis (2017); Bulchandani et al. (2017) (see Supplemental Material Sup ). For the purposes of comparison with DMRG, we restrict to $\Delta=\frac{1}{2}$ (other values of $\Delta$ were considered in previous works Bertini et al. (2016); Ilievski and De Nardis (2017); Bulchandani et al. (2017) for non-equilibrium steady-states). We also focus on non-equilibrium energy transport, in particular the evolution of local energy density, given by $n_{E}(x,t)=\sum_{j=1}^{N_{t}}\int d\lambda\,e_{j}(\lambda)\rho_{j}(x,t,\lambda)$ in the hydrodynamic limit (see Sup ).

To illustrate the range of validity of the method, we consider a strongly non-equilibrium example, namely the Gaussian initial temperature profile

[TABLE]

with $\beta_{0}>\beta_{M}$. Physically speaking, this corresponds to a perturbation $\beta_{M}^{-1}$ of a background temperature $\beta_{0}^{-1}$, localized over a typical length $\sim L$. We first illustrate the convergence of our numerical scheme (6) by taking such a Gaussian initial state and letting the time step $dt\to 0$. This is depicted in Fig. 1. As $dt$ is lowered, the numerical solution at long times (say, $t=20$) converges very quickly, and remarkably, even one-step or two-step schemes (e.g. $dt=20$ or $dt=10$) yield good approximations to the converged solution.

We now compare the predictions of the Bethe-Boltzmann equation against DMRG calculations White (1992); Schollwoeck (2011), with the initial condition (9) prepared using standard finite-temperature methods Karrasch et al. (2012); Barthel (2013); Kennes and Karrasch (2016). This is shown in Fig. 2. We find an excellent agreement between DMRG and hydrodynamic results with $dt=2.5$ for quite different initial temperature profiles — Gaussian (9) and $\tanh$ $\beta(x)=(\beta_{0}+\beta_{M})/2+(\beta_{0}-\beta_{M})/2\times\tanh(x/L)$ functions. Provided that the initial condition is smooth enough for the DMRG and the thermodynamic Bethe ansatz calculations to agree at $t=0$, subsequent agreement at later times is essentially perfect. In fact, at low temperatures where it is hard to obtain smooth initial conditions in DMRG, the main source of error comes from slight disagreements in the initial conditions between the two approaches.

Hydrodynamics for the interacting Bose gas.

An existing hydrodynamic description of the quasi-one-dimensional Bose gas, based on a local density approximation for the first three conserved charges of the Lieb-Liniger model, has proved effective for capturing non-equilibrium dynamics for such systems Kolomeisky and Straley (1992); Kolomeisky et al. (2000); Dunjko et al. (2001); Ohberg and Santos (2002); Pedri et al. (2003); Campbell et al. (2015); Bouchoule et al. (2016); De Rosi and Stringari (2016). In its domain of physical validity Kolomeisky et al. (2000); Kheruntsyan et al. (2003), the present approach improves the existing theory by allowing for local GGE initial conditions and respecting all higher conservation laws implied by integrability. The extension of this hydrodynamic approach to other important aspects of 1D Bose gas physics, including dynamics in external potentials Castin and Dum (1996); Kagan et al. (1997); Fang et al. (2014); Bouchoule et al. (2016); De Rosi and Stringari (2016); Atas et al. (2017) and correlation functions Gangardt and Shlyapnikov (2003); Kheruntsyan et al. (2003); Xu and Rigol (2017); Atas et al. (2017), is currently an active area of research  Doyon and Yoshimura (2016); Doyon and Spohn (2017).

Two non-equilibrium quenches which are of particular experimental interest are sudden expansions of Bose gases into vacuum and collisions of clouds of ultracold bosons Gobert et al. (2005); Kinoshita et al. (2006); Langer et al. (2009, 2011); Ronzheimer et al. (2013); Karrasch et al. (2014); Bolech et al. (2012); Langer et al. (2012); Vidmar et al. (2013); Steinigeweg et al. (2017); Karrasch et al. (2017). We find that our numerical solution to Eq. (6) converges for initial conditions modeling both of these scenarios. From the resulting evolution in $\theta(x,t,k)$, one can track the evolution of any local conserved charge of the model. In Fig. 3, we plot the time evolution of particle and energy densities (defined by $n=Q_{1}$ and $n_{E}=Q_{2}$ in Eq. (2) respectively) in a Lieb-Liniger model with interaction strength $c=1$ from free expansion and collision type initial conditions. The initial states are prepared at temperature $T=1$ using a chemical potential profile interpolating between $\mu=5$ inside a box and a large negative value $\mu=-50$ outside, with the edges of the box smoothed out using tanh functions. For the collision protocol, two clouds of bosons initially prepared as in the free expansion quench are given opposite momenta $k=\pm k_{0}$ with $k_{0}=2.5$.

General features.

We hope that the previous examples have established that the Bethe-Boltzmann equation is a valuable tool for specific computations. More generally, it is natural to ask which phenomena in integrable models are missed by this hydrodynamical approach and how it differs qualitatively from conventional hydrodynamics. Clearly there is a significant assumption that the initial condition is well described by a local generalized Gibbs ensemble (GGE) Rigol et al. (2008); Barthel and Schollwöck (2008). In fact, the “thermalization problem” for quantum integrable models, namely the question of determining the GGE to which a given quantum state converges in the long-time limit, is a difficult problem which remains unsolved in general. Nevertheless, when local equilibrium initial conditions can be imposed, hydrodynamics does seem to capture the leading behavior at long time and length scales, as illustrated in the examples above. There can also be other important subleading behaviors, beyond the approach to local equilibrium. An example is the behavior of the Lieb-Liniger model in the low-temperature limit, where it can be described by a conformal field theory or bosonization. Conformal invariance Sotiriadis and Cardy (2008) and other methods Langmann et al. (2017) both predict a Schwarzian derivative term in the time evolution from locally thermal initial conditions, which might be an example of subleading behavior beyond hydrodynamics.

An important difference between generalized hydrodynamics in integrable models and conventional hydrodynamics concerns reversibility. The collision term in the standard Boltzmann equation induces dissipation and an increase of entropy. The Bethe-Boltzmann equation is dissipationless, and in fact its time evolution is reversible. The action of microscopic time reversal on a particular state is to invert all pseudo-momenta $k$ while fixing space, so that the time-reversed pseudo-momenta are given by $\rho_{k}^{\prime}(x,t)=\rho_{-k}(x,t)$. One can show despite the complex form of the velocity $v_{k}$ in the Bethe-Boltzmann equations, it undergoes a simple sign change under this microscopic time reversal transformation 444This follows by the dressing equations for energy and momentum, together with the fact that TBA kernels are difference kernels in the rapidity parameterization..

So integrable models again present some surprises compared to ordinary kinetic theory: while the description by a local GGE is certainly a great reduction in complexity compared to an arbitrary quantum state and hence irreversible, the flow in the space of local GGEs described by the Bethe-Boltzmann equation is reversible. Presumably this means that truly diffusive behavior, as is believed to appear for example in the gapped phase of the XXZ model, lies beyond this equation; for linear-response spin transport in XXZ model, which involves both ballistic and diffusive components Sirker et al. (2009, 2011), hydrodynamics captures correctly the ballistic part Ilievski and De Nardis (2017); Bulchandani et al. (2017).

There are many possible mathematical questions regarding the existence and structure of solutions to the Bethe-Boltzmann equation Sup , but we hope that the above results demonstrate its practical utility for applications to physics. It can be used as a starting point for comparison for other methods for quantum dynamics, or for incorporating integrability-breaking terms or driving. The hydrodynamical theory of integrable models is one of many examples in recent years of how the old vine planted by Yang and Yang Yang and Yang (1969) continues to bear fruit.

Acknowledgements: The authors acknowledge support from NSF DMR1507141 and a Simons Investigatorship (V. B. B. and J. E. M.), the Emmy-Noether program of the Deutsche Forschungsgemeinschaft under KA 3360/2-1 (C. K.) and the Quantum Materials Program supported by the Director, Office of Science, Office of Basic Energy Sciences, Materials Sciences and Engineering Division, of the US DOE under Contract No. DE-AC02-05CH11231 (R. V.).

Appendix A Thermodynamic Bethe Ansatz and the Bethe-Boltzmann Equations

The Bethe-Boltzmann hydrodynamic approach is fundamentally based on the theory of quasiparticle excitations in thermodynamic Bethe ansatz Bertini et al. (2016); Castro-Alvaredo et al. (2016). In this Appendix, we summarize the thermodynamic Bethe ansatz description of quasiparticle excitations of the Lieb-Liniger model and gapless XXZ model, providing explicit formulas for the quasiparticle velocities required in the main text.

A.1 Lieb-Liniger model

In the Lieb-Liniger model, the group velocity of quasiparticle excitations about thermodynamic equilibrium is given in terms of the dressed energy and momentum by $v(k)=E^{\prime}(k)/P^{\prime}(k)$, with dressed charges $Q$ related to bare charges $q$ via the integral equation

[TABLE]

Here, $\vartheta$ denotes the Fermi factor of the underlying equilibrium state, and the Lieb-Liniger kernel is given in terms of the interaction strength $c$ by

[TABLE]

It is useful to introduce operators

[TABLE]

Then from the bare values $p(k)=k$ and $e(k)=k^{2}/2$ we can take the formal inverse of Eq. (10) (i.e. expand in a Neumann series) to yield the dressed values

[TABLE]

giving rise to the formula

[TABLE]

If we now consider a “local TBA” type approximation Castro-Alvaredo et al. (2016); Bertini et al. (2016), which amounts to the assumption that local energy density is well-defined Bulchandani et al. (2017), then we can postulate a local Bethe equation

[TABLE]

at every point; this may be taken as a definition of the local Fermi factor $\vartheta(x,t,k)$, which in turn yields the local quasiparticle velocity,

[TABLE]

In summary, the Bethe-Boltzmann equation is shorthand for the hierarchy of equations

[TABLE]

which together comprise a conservation law with self-consistently determined velocity. Using the local Bethe equations (14), we can instead change variables to the local Fermi factor $\vartheta(x,t,k)$, yielding the self-consistent system

[TABLE]

This form of the equation follows by various TBA identities Bertini et al. (2016); Bulchandani et al. (2017).

A.2 Gapless XXZ Model

Recall that the Hamiltonian for the spin-$1/2$ XXZ chain on $N$ sites, in zero external field, is given by

[TABLE]

We take periodic boundary conditions $S_{N}\equiv S_{N+1}$, set the coupling to $J=1$, and parameterize the anisotropy of the theory by $\Delta=\cos\gamma$. We also assume in the following that the model is in its gapless phase, i.e. $-1<\Delta<1$. The Bethe-Boltzmann equation for this system is discussed in detail in Bertini et al. (2016). Here, we summarize the main results. The only material difference compared to the Lieb-Liniger gas is that one must now account for the “strings” of bound states appearing in the thermodynamic limit. Let us therefore define a string of type $j$ to be a an ordered pair, $(n_{j},v_{j})$, where $n_{j}$ is the number of spin-flips comprising the string and $v_{j}$ is its parity. Suppose that there are $N_{t}$ string types in total so that $j\in\{1,2,\ldots,N_{t}\}$, and $M_{j}$ strings of type $j$. Let $M$ denote the total number of spin-flips in the system. Then by definition, $\sum_{j=1}^{N_{t}}M_{j}n_{j}=M$. Upon fixing a string type $j$, we denote the rapidities of a given string $\alpha\in\{1,2,\ldots,M_{j}\}$ within that type by

[TABLE]

where $m_{j}$ denotes the string length. It is useful to define functions

[TABLE]

and

[TABLE]

The quasiparticle velocities in each string are given in terms of dressed energy and momentum by the expression $v_{j}(\lambda)=E^{\prime}_{j}(\lambda)/P^{\prime}_{j}(\lambda)$. It can be shown Bertini et al. (2016) that the dressed charge for any given quasiparticle excitation is related to the bare charge via the $N_{t}$ coupled integral equations

[TABLE]

where the $\vartheta_{j}$ denote the Fermi factor for strings of type $j$ and $\sigma_{j}=\mathrm{sgn}(v_{j})$ (see Bertini et al. (2016) for details). Then from the bare values $p_{j}(\lambda)=2\pi a_{j}(\lambda)$, $e_{j}(k)=-Aa_{j}(\lambda)$ for the momentum and energy of string $j$, where $A=-2\pi J\sin\gamma/\gamma$, we obtain dressed momenta and energies

[TABLE]

One can deduce that

[TABLE]

and there are now $N_{t}$ Bethe-Boltzmann equations, one for each string, with associated quasimomentum densities $\rho_{j}(x,t,\lambda)$. Finally, we note that the XXZ chain considered in the main text, with $J=1$ and $\Delta=0.5$ possesses three strings in the thermodynamic limit, with parameters $(n_{1},v_{1})=(1,1)$, $(n_{2},v_{2})=(2,1)$ and $(n_{3},v_{3})=(1,-1)$.

Appendix B Details of Numerical Implementation

In this section, we describe the numerical scheme used in the main text to solve the Bethe-Boltzmann equation. For simplicity, we restrict our presentation to the Lieb-Liniger gas. Thus consider the advection form of the Bethe-Boltzmann equation

[TABLE]

As it stands, this is an infinite-dimensional system of coupled, non-linear equations. The first step in achieving a numerical solution is to discretize in the quasimomentum variable $k$. Let us therefore introduce a $k$-space cutoff $\Lambda$ and discretize the infinite-dimensional system (25) at $N$ k-space points $-\Lambda<k_{1}<k_{2}<\ldots<k_{N}=\Lambda$. It is also useful to set $k_{0}=-\Lambda$. Then functions map to vectors and kernels map to matrices; for example, we have discrete analogues

[TABLE]

of charges, Fermi factors and the Lieb-Liniger kernel respectively. Let us also introduce the matrix $A_{ij}(\theta)=K_{ij}\theta_{j}$. Then upon discretizing, derivatives of dressed charges $Q^{\prime}$ are related to bare ones $q^{\prime}$ via the matrix equation

[TABLE]

It is useful to introduce the dressing operator $U_{ij}(\theta)=(1+A(\theta))^{-1}_{ij}$ and dressed derivatives of energy and momenta $E^{\prime}_{i}(\theta)=\sum_{j}U_{ij}(\theta)k_{j}$, $P^{\prime}_{i}(\theta)=\sum_{j}U_{ij}(\theta)$. Then the discretized quasiparticle velocities can be written as $v_{i}(\theta)=E^{\prime}_{i}(\theta)/P^{\prime}_{i}(\theta)$. In terms of these functions, the discretized advection equation (25) reads

[TABLE]

We note that to derive the results of the main text, we used a Simpson’s rule method to compute the action of the dressing operator in evaluating the $v_{i}$. For large $N$, this matches the simpler discretization described here. Now consider the initial value problem for the equation (28), with initial condition $\theta_{i}(x,0)=\phi_{i}(x)$ at $t=0$. We propose the following backwards implicit scheme:

Choose a time step $dt>0$. 2. 2.

At step $n=0$, set $\theta^{0}_{i}(x)=\phi_{i}(x)$. 3. 3.

For steps $n\geq 1$, solve for $\theta^{n}_{i}(x)$ satisfying the implicit equation $\theta^{n}_{i}(x)=\theta^{n-1}_{i}(x-v_{i}[\theta^{n}]dt)$ by numerical iteration. 4. 4.

Identify $\theta_{i}(x,n\cdot dt)$ with $\theta_{i}^{n}(x)$.

This was shown to define a first-order numerical scheme for (28) in a previous work Bulchandani et al. (2017). In practice, some level of spatial discretization is also necessary, and to derive the results of the main text, we applied the scheme at time step $n$ to a cubic spline interpolation of the data at step $n-1$. This ensured accuracy of the past light cone traced out by the backwards scheme at each time step.

Appendix C Semi-Hamiltonian Structure of Bethe-Boltzmann Hydrodynamics

C.1 Definition and Properties

Rather surprisingly, the finite discretizations (28) of the Bethe-Boltzmann equation define a “semi-Hamiltonian” or “rich” system of non-linear equations. To see this, recall that a semi-Hamiltonian system of hydrodynamic type Dubrovin and Novikov (1983); Tsarëv (1991); Serre (1989); Ferapontov (1991) is a system of coupled non-linear PDEs

[TABLE]

whose characteristic velocities satisfy the properties

$\partial_{j}v_{j}=0$ for $j=1,2,\ldots,N$ 2. 2.

$\partial_{k}\left[\partial_{j}v_{i}/(v_{j}-v_{i})\right]=\partial_{j}\left[\partial_{k}v_{i}/(v_{k}-v_{i})\right]$ for all $i\neq j\neq k\neq i$,

where $\partial_{j}v_{i}\mathrel{\mathop{:}}=\partial v_{i}/\partial\theta_{j}$. If the system (29) is semi-Hamiltonian in this sense, then it possesses infinitely many flows commuting with the velocity field $\mathbf{v}$. This leads to a range of surprising geometrical and analytic properties Tsarëv (1991); Serre (1989). Remarkably, it can be shown that the system (28) is semi-Hamiltonian for any finite discretization length $N$. This is a non-trivial consequence of the structure of excitation dressing in TBA, and we briefly sketch its derivation. Let $\mathcal{S}\subset\mathbb{R}^{N}$ denote an open set on which the dressing operator $U$ exists. Defining the operator $\alpha_{ij}(\theta)=-\sum_{l}U_{il}(\theta)K_{lj}$, it can be shown that on $\mathcal{S}$,

[TABLE]

and similarly for $E^{\prime}$. From this, one may deduce that $\partial_{j}v_{j}=0$, together with the formula

[TABLE]

By an identical argument $\partial_{k}v_{i}/(v_{k}-v_{i})=\partial_{k}\log P^{\prime}_{i}$, and the second part of the semi-Hamiltonian property follows by commutativity of partial derivatives. In view of the theory of semi-Hamiltonian systems, we have actually shown slightly more. Recall that such systems have infinitely many flows $\mathbf{w}$ commuting with $\mathbf{v}$ and given by solutions to the system of linear equations

[TABLE]

with no summation. Since the result Eq. (31) above is independent of the bare energy, let us choose $N$ linearly independent vectors $\mathbf{a}^{(n)}\in\mathbb{R}^{N}$ such that $\mathbf{a}^{(0)}=\mathbf{p}^{\prime}$ and $\mathbf{a}^{(1)}=\mathbf{e}^{\prime}$. Then there are $N$ commuting flow fields $\mathbf{w}^{(n)}$ given by

[TABLE]

which satisfy (32) by our earlier reasoning. Moreover, these define linearly independent vector fields on $\mathcal{S}$ and each gives rise to a semi-Hamiltonian system in its own right. In the special case of Lieb-Liniger generalized hydrodynamics, a basis for these flows may be obtained from the first $N$ charges of the model, taking $a^{(n)}_{i}={q^{\prime}}^{(n)}_{i}=k_{i}^{n-1}$ for $n=0,1,\ldots,N-1$. This implies an intriguing geometrical relationship between the conserved charges of integrable models and the semi-Hamiltonian vector fields of their hydrodynamics, whose consequences in the limit $N\to\infty$ would be interesting to explore. It is also worth noting that semi-Hamiltonian structures arise naturally in the Bogoliubov or Whitham phase-averaging of classical integrable PDEs Dubrovin and Novikov (1983).

C.2 Solution by Quadrature

In general, a semi-Hamiltonian structure implies the possibility of a solution “by quadrature” Ferapontov (1991). Although this only works for a restricted class of initial conditions (those for which $\phi_{j}:\mathbb{R}\to\mathbb{R}$ has differentiable inverse), it is surprising to us that any exact solutions to the Bethe-Boltzmann equations exist at all. To motivate the construction, observe that we will have solved the discretized advection equation if we can find $\theta_{j}(x,t)$ and functions $f_{j}(\theta)$ solving the system

[TABLE]

of $2N$ equations. One can show that the semi-Hamiltonian property is equivalent to integrability of this system. Let us see how this works in practice. Assuming that the $f_{j}(y)$ are nowhere vanishing, one can apply the inverse dressing operator to find that

[TABLE]

Let us now consider the Cauchy problem for $t>0$, with initial state $\theta(x,0)=\phi(x)$ at $t=0$. Upon integrating the above relation along a path between $(0,0)$ and $(x,t)$, we obtain the implicit equation

[TABLE]

When the $f_{j}(y)$ are continuous, it follows by the chain rule and the implicit function theorem that this solves the system (34), which implies the Bethe-Boltzmann equation. If, in addition, the initial conditions $\phi_{j}$ each possess a differentiable inverse, one can set $t=0$ and differentiate in $x$ to fix the $N$ functional degrees of freedom $f_{i}(\theta)$.

C.3 Geometrical Interpretation

We now recall some geometrical properties of semi-Hamiltonian systems, due to Tsarëv Tsarëv (1991). The main result is that for a semi-Hamiltonian system, there exists a diagonal pseudo-Riemannian metric tensor $\sum_{i=1}^{N}g_{ii}(\theta)d\theta_{i}\otimes d\theta_{i}$ on the space of states, satisfying

[TABLE]

To derive this, one defines first defines coefficients of a symmetric connection according to

[TABLE]

The semi-Hamiltonian property then reads

[TABLE]

which is an integrability condition for the metric compatibility condition $\partial_{j}\log\sqrt{g_{ii}(\theta)}=\Gamma^{i}_{ij}(\theta)$. This corresponds to the vanishing of certain components of the Riemann curvature tensor. For Lieb-Liniger hydrodynamics, we can write Eq. (37) in terms of the total density of states $\rho^{t}$, as

[TABLE]

A solution is given by the metric tensor

[TABLE]

on $\mathcal{S}$. One can check that this is flat at the free Fermion point $c=\infty$, but not generically elsewhere. Thus interactions give rise to curvature on state space, which prevents Bethe-Boltzmann hydrodynamics from being Hamiltonian in the conventional sense Tsarëv (1991)555We were informed by B. Doyon before completion of our paper about forthcoming work on time-evolution quadrature that also points out the non-Hamiltonian property, but does not use the semi-Hamiltonian structure arising from integrability.. In the terminology of differential geometry, the components of the total density of states define the Lamé coefficients of the state-space metric Ibragimov (1995).

C.4 Continuum Limit

In fact, the continuum advection equation (17) appears to satisfy an infinite-dimensional analogue of the semi-Hamiltonian property with partial derivatives replaced by functional derivatives, namely:

$\frac{\delta v[\theta](k)}{\delta\theta(k)}=0$ for all $k\in\mathbb{R}$ 2. 2.

$\frac{\delta}{\delta\theta(k^{\prime\prime})}\left[\frac{1}{v(k^{\prime})-v(k)}\frac{\delta v[\theta](k)}{\delta\theta(k^{\prime})}\right]=\frac{\delta}{\delta\theta(k^{\prime})}\left[\frac{1}{v(k^{\prime\prime})-v(k)}\frac{\delta v[\theta](k)}{\delta\theta(k^{\prime\prime})}\right]$ for all $k\neq k^{\prime}\neq k^{\prime\prime}\neq k$.

To our knowledge, the geometrical theory of such systems has not yet been developed. The first property generalizes linear degeneracy in the sense of Lax Lax (1957), and when combined with monotonicity, implies existence and uniqueness for solutions to a certain self-consistent ansatz proposed for two-reservoir steady states Castro-Alvaredo et al. (2016); Bertini et al. (2016). The second property again appears to give rise to a solution by quadrature. To see this, one can simply take the continuum limit of (36), yielding

[TABLE]

These solutions to the Bethe-Boltzmann equation are parameterized by a functional degree of freedom $f(\theta,k)$, in two variables. Differentiating and using invertibility of the TBA kernel, one can check that this solves the system

[TABLE]

which implies the continuum Bethe-Boltzmann equation. We emphasize that this exact solution is, at present, mainly of formal interest. In particular, it is logically independent from the method used to derive the results of the main text, whose details are described in a previous work Bulchandani et al. (2017).