Bethe-Boltzmann Hydrodynamics and Spin Transport in the XXZ Chain

TL;DR

This paper develops a numerical hydrodynamic approach to study spin and heat transport in the XXZ quantum spin chain, providing new insights into non-equilibrium steady states and spin Drude weight calculations.

Contribution

It introduces a numerical scheme for hydrodynamic equations in integrable systems and applies it to accurately compute spin Drude weight in the XXZ chain.

Findings

Hydrodynamic methods agree with DMRG results for spin Drude weight.

The approach describes non-equilibrium steady states with ballistic transport.

Numerical scheme enables finite-time analysis of integrable quantum systems.

Abstract

Quantum integrable systems, such as the interacting Bose gas in one dimension and the XXZ quantum spin chain, have an extensive number of local conserved quantities that endow them with exotic thermalization and transport properties. We discuss recently introduced hydrodynamic approaches for such integrable systems from the viewpoint of kinetic theory and extend the previous works by proposing a numerical scheme to solve the hydrodynamic equations for finite times and arbitrary locally equilibrated initial conditions. We then discuss how such methods can be applied to describe non-equilibrium steady states involving ballistic heat and spin currents. In particular, we show that the spin Drude weight in the XXZ chain, previously accessible only by rigorous techniques of limited scope or controversial thermodynamic Bethe ansatz arguments, may be evaluated from hydrodynamics in very good agreement with density-matrix renormalization group calculations.

Equations

$$\hat{H}=-\sum _{j=1}^N{\partial }_j^2+\sum _{j<k}2c\delta (x_j-x_k).$$

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

$${\partial }_t\rho (x,t,k)+{\partial }_{x}j(x,t,k)=0,$$

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

$${\varphi }_t=6\varphi {\varphi }_x-{\varphi }_{xxx},$$

$${\partial }_t\rho (x,t,\lambda )+{\partial }_x(\rho (x,t,\lambda )v[\rho ](x,t,\lambda ))=0,$$

$$v(\lambda )=v_0(\lambda )+{\int }_{-\infty }^{\infty }d\mu \Delta x(\lambda ,\mu )\rho (\mu )[v(\lambda )-v(\mu )],$$

$$v(\lambda )=v_0(\lambda )+\frac{1}{p^{\prime }(\lambda )}{\int }_{-\infty }^{\infty }d\mu {\phi }^{\prime }(\lambda -\mu )\rho (\mu )[v(\lambda )-v(\mu )],$$

$${\phi }^{\prime }(\lambda -\mu )=i\frac{d}{d\lambda }\ln S(\lambda -\mu ).$$

$$\Delta x(\lambda ,\mu )=\frac{1}{p^{\prime }(\lambda )}{\phi }^{\prime }(\lambda -\mu ).$$

$$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\rho +{\partial }_x(\rho v[\rho ])=0.$$

$$\rho [\hat{\vartheta }(x,t)](k)$$

$$(\rho v)[\hat{\vartheta }(x,t)](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)}$$

$${\partial }_t\vartheta +v[\vartheta ]{\partial }_x\vartheta =0.$$

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

$$\vartheta (x,0,k)$$

$$\varphi (x,k)={\theta }_L(k)H(-x)+{\theta }_R(k)H(x),$$

$$\vartheta (x,t,k)=\varphi (x-v[\hat{\vartheta }](k)t,k),$$

$$\vartheta (x,t,k)$$

$$+{\theta }_R(k)H(x-v[\hat{\vartheta }(x,t)](k)t).$$

$$v[\hat{\vartheta }(x,t)](k)t-x=0,$$

$$v[\hat{\vartheta }(x,t)](k^{\ast }(x,t))=x/t.$$

$$\vartheta (x,t,k)={\theta }_L(k)H(k-k^{\ast }(x,t))+{\theta }_R(k)H(k^{\ast }(x,t)-k).$$

$$\vartheta (\zeta ,k)=$$

$$v[\hat{\vartheta }(\zeta )](k^{\ast }(\zeta ))$$

Full text

I Introduction

The study of many-body quantum systems far away from equilibrium conditions poses a considerable challenge for theory, even for quantum integrable systems whose equilibrium properties may be computed exactly. Such systems, which include the Heisenberg antiferromagnet and the Lieb-Liniger gas in one dimension, possess an extensive number of conserved quantities, which prevent them from thermalizing like generic ergodic systems and lead to dissipation-less transport properties. Under unitary evolution, the local properties of these systems are believed to tend to a generalized Gibbs ensemble (GGE) Rigol et al. (2008); Barthel and Schollwöck (2008) at long times, containing in principle all the independent conserved quantities, not just particle number and energy as in the standard Gibbs ensemble. The rest of the system acts as an “unusual bath”.

However, this convergence can be rather subtle, as the GGEs constructed using the standard conserved quantities of the spin-$1/2$ XXZ chain Pozsgay (2013); Fagotti and Essler (2013); Fagotti et al. (2014) were shown Wouters et al. (2014); Pozsgay et al. (2014) to fail to reproduce the correct steady states obtained either numerically or using the so-called quench action method Caux and Essler (2013). This paradox can only be resolved Ilievski et al. (2015a) by taking into account “hidden” (non-standard) quasi-local conserved quantities Prosen (2011); Prosen and Ilievski (2013); Prosen (2014); Pereira et al. (2014); Ilievski et al. (2015b) in the GGE. More recently, the nature of these new conserved quantities and their relation to the pseudo-momentum distributions used in equilibrium thermodynamic Bethe ansatz calculations was clarified Ilievski et al. (2016, 2017); De Luca et al. (2017).

Nevertheless, when combined with the infinite number of Lagrange multipliers which must be fixed from the initial state, this complication makes the GGE approach quite cumbersome for the study of interesting non-equilibrium dynamics in integrable models, particularly those arising from spatially non-uniform states. Moreover, generalized Gibbs ensembles are by definition restricted to describe steady-states, making a general approach able to deal with non-equilibrium finite-time dynamics desirable. Very recently, a more practical hydrodynamic approach based on a semi-classical quasi-particle picture was introduced Castro-Alvaredo et al. (2016); Bertini et al. (2016), and conjectured to yield exact results for the long-time scaling limit reached from the “two-reservoir quench,” an initial condition of two semi-infinite reservoirs connected at the origin. One way to understand these approaches is that they reflect an equivalence between the hydrodynamical and Boltzmann-type descriptions of integrable models.

In standard statistical mechanics, the Boltzmann equation involves the full one-particle distribution function, i.e., a function of momentum at each point in space and time. Standard hydrodynamics contains considerably less information, as only three quantities survive at each point in space and time: the local particle density, momentum density, and energy density. For integrable models, there is a more fundamental relationship between the distribution function $\rho(x,t,k)$ over pseudo-momentum $k$ (the analogue of ordinary momentum for integrable models) and the full set of conserved quantities. This is most easily seen in the Lieb-Liniger model, where the conserved quantities are just moments of the pseudo-momentum distribution Korepin et al. (1997). It is much less clear that it holds for the XXZ model, whose conserved quantities have a rather subtle structure (see Refs. Ilievski et al., 2016, 2017; De Luca et al., 2017 for recent developments), but we show in this paper that a hydrodynamical description is successful even for observables that are sensitive to the newly discovered quasi-local charges. This result is consistent with recent work demonstrating the equivalence of the GGE and thermodynamic Bethe ansatz (TBA) pictures of equilibrium states in the XXZ model De Luca et al. (2017), although some subtleties of interpretation still remainPozsgay et al. (2017).

The first derivations of the recent hydrodynamic formalism were based upon the observation 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. The latter conservation law, previously obtained in Refs. Castro-Alvaredo et al., 2016; Bertini et al., 2016, can alternatively be viewed as a dissipationless Boltzmann equation for the pseudo-momentum distribution, where the effect of collisions is to dress the bare velocity of particle-type excitations. It is precisely this kinetic equation that was solved in Refs. Castro-Alvaredo et al., 2016; Bertini et al., 2016, rather than an infinite set of hydrodynamic equations, and in recognition of its crucial role within the hydrodynamic formalism, we call this equation the “Bethe-Boltzmann equation”. Our terminology seeks to emphasize that this equation can be understood phenomenologically as a Boltzmann equation whose collision term captures the Wigner time-delay Wigner (1955) due to Bethe ansatz phase shifts (this is one way to interpret the velocity dressing equation, Eq. (8)).

From this perspective, the Bethe-Boltzmann equation emerges as the quantum analogue of an established kinetic theory of classical soliton gasesPercus (1969); Zakharov (1971); Boldrighini et al. (1983); El (2003); El and Kamchatnov (2005). Our viewpoint is also consistent with a recently obtained “molecular dynamics” realization of this equationDoyon et al. (2017). This viewpoint is useful because it provides some additional physical justification for applying this equation at finite- time and length- scales, which, although admittedly less rigorous than the discussion of quasi-stationary states in Refs. Castro-Alvaredo et al., 2016; Bertini et al., 2016, is closer in spirit to the semi-classical theory of quasiparticle transport that underpins more conventional solid-state physicsAshcroft et al. (2016).

This relatively simple picture of the Bethe-Boltzmann equation raises the question of whether predictions obtained from the hydrodynamic formalism are in fact correct. Fortunately, there exist a plethora of predictions for non-equilibrium evolution of quantum integrable models from a range of initial conditions, against which hydrodynamics can be compared. Among these, the two-reservoir quench mentioned above has had an enduring popularity; this set-up consists of two half-infinite systems prepared in thermal equilibrium with different temperatures and chemical potentials (or more generally, two half-infinite homogeneous reservoirs), joined together at time $t=0$ and allowed to evolve unitarily for $t>0$ according to the Schrödinger equation Sotiriadis and Cardy (2008); Bernard and Doyon (2012); Karrasch et al. (2013a); 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). In light of generalized hydrodynamics, this gives rise to the intriguing possibility of deriving aspects of non-equilibrium quantum transport from an essentially classical equation. Non-equilibrium thermal transport in the XXZ spin chain from hydrodynamics was studied in detail in Ref. Bertini et al., 2016 and was thoroughly compared to matrix-product state numerics White (1992); Schollwoeck (2011).

The case of spin transport is somewhat more complex theoretically, even at the level of linear response, as the presence of ballistic spin currents requires the existence of quasi-local conserved quantities Prosen (2011) going beyond the standard local conserved quantities of the XXZ spin chain. For conventional (linear response) transport, the ballistic component of the spin current is characterized by the spin Drude weight, which measures the degree of divergence of the zero-frequency spin conductivity. The Bethe ansatz calculation of the spin Drude weight Zotos (1999); Benz et al. (2005) has attracted a lot of attention in the past, and remains controversial. In fact, before the discovery of the quasi-local conserved quantities mentioned above, it was even debated whether this Drude weight was non-zero at finite temperature (see e.g. Refs. Heidrich-Meisner et al., 2003; Sirker et al., 2009, 2011; Karrasch et al., 2012, 2013b and references therein). In this paper, we devise a method to compute the spin Drude weight from the Bethe-Boltzmann equation, and show that the value of the resulting Drude weight is compatible with known exact results and with density-matrix renormalization group (DMRG) calculations. This is remarkable because it indicates that the non-equilibrium steady-state predicted by hydrodynamics takes all quasi-local conserved charges into account (see also Ref.  De Luca et al., 2017), even though such charges can effectively be ignored at the level of the Bethe-Boltzmann equation, which could have been deduced on phenomenological grounds long before the discovery of quasi-local conservation laws.

We proceed as follows. In Section II, we summarize the background from Refs. Castro-Alvaredo et al., 2016; Bertini et al., 2016 needed to formulate our main results on the spin Drude weight. We also discuss how the Bethe-Boltzmann equation can be viewed as a phenomenological extension of older kinetic theories for classical soliton gases to quantum integrable models (see also Doyon et al., 2017), which provides additional physical motivation for applications of the hydrodynamic approach at finite time- and length- scales. We mostly defer the detailed calculation of this finite-time hydrodynamics in physically relevant examples for a subsequent publication Bulchandani et al. (2017), except for one example, a thermal expansion in the XXZ chain, in Figure 4 below. This extends the analysis in previous worksCastro-Alvaredo et al. (2016); Bertini et al. (2016) by providing the first numerical checks of the hydrodynamic approach beyond scale-invariant steady-states. In Section III we present our main results, on the evaluation of the spin Drude weight $D_{s}(T)$.

II The Bethe-Boltzmann Equation

The reader familiar with Refs. 18 and 19 can skip Section II without loss of continuity. Here, we summarize the background from Refs. 18 and 19 that is necessary for formulating our main results on the spin Drude weight, in Section III. We also discuss the intimate connection between the Bethe-Boltzmann equation and various existing kinetic equations describing the dynamics of classical soliton gasesPercus (1969); Zakharov (1971); Boldrighini et al. (1983); El (2003); El and Kamchatnov (2005).

II.1 Motivation

To provide some intuition, we first sketch the Bethe-Boltzmann formalism for the one-dimensional Bose gas with delta-function interactions. For $N$ bosons on a line, with interaction strength $c$, the Hamiltonian may be written as

[TABLE]

This model, also called the Lieb-Liniger gas, is the simplest non-trivial integrable model, and amongst the entire class of such models has the merit of being the most relevant to experimental physics. The physics of quasiparticle excitations in this model, as obtained from TBA, is summarized in Appendix B. Let us now consider a Lieb-Liniger gas on a line of length $L$, consisting of $N$ particles. In the thermodynamic limit as $N,L\to\infty$, we assume that the system may be characterized by a local density of occupied pseudo-momenta, $\rho(x,t,k)$, giving rise to a locally varying particle density

[TABLE]

Physically speaking, this amounts to coarse-graining our line into cells of length $l\ll L$, such that on each cell, the gas lies in a macrostate fixed by the average particle and hole densities over that cell. We additionally postulate that “occupied quantum numbers are locally conserved”. This implies a local conservation law of the form

[TABLE]

for some current $j(x,t,k)$ to be determined. To obtain a specific form for $j$, observe that a physically natural velocity scale for the transport of quantum numbers is given by the quasiparticle velocity $v[\rho](x,t,k)$ in each cell, which generally depends on the local occupation number $\{\rho(x,t,k):k\in\mathbb{R}\}$. This yields

[TABLE]

which we shall henceforth refer to as the Bethe-Boltzmann equation. This equation has the structure of a conservation law for the local pseudo-momentum distribution $\rho$, which is how it was first introduced in Refs. Castro-Alvaredo et al., 2016; Bertini et al., 2016, where it was derived from an infinite set of “generalized hydrodynamic” equations. However, Eq. (4) can be obtained directly by the above phenomenological arguments, and from this viewpoint, defines a Boltzmann-type equation for the dissipationless transport of quasiparticles. By varying the number of local densities appropriately, an equation of this type may be formulated for any quantum integrable model. We now discuss this formalism in more detail, and in particular, show how it fits into an existing kinetic theory for solitons of classically integrable PDEs.

II.2 Relation with the Kinetic Theory of Solitons of Classically Integrable PDEs

It has long been known that there exists a class of non-linear partial differentiable equations, called “integrable”, which admit stably propagating “solitonic” solutions. For example, the KdV equation

[TABLE]

is of this type. In 1971, shortly after the discovery of the classical inverse scattering method for solving such equations, ZakharovZakharov (1971) considered the kinetic theory of their solitons and proposed an integro-differential equation for a gas of KdV solitons in the dilute limit, of the form

[TABLE]

where $\rho(x,t,\lambda)$ denotes the local density of solitons with spectral parameter $\lambda$ and $v$ their effective velocity, given as a certain functional linear in $\rho$. Several years later, an extension of this equation beyond the dilute limit was obtained from an infinite-genus limit of the Whitham equationsEl (2003). It was found that the equation (6) held in general, provided that Zakharov’s explicit formula for $v$ was replaced by an implicit integral equation

[TABLE]

with $v_{0}(\lambda)$ the velocity of a single soliton with spectral parameter $\lambda$ and $\Delta x(\lambda,\mu)$ the asymptotic position shift after a collision of two KdV solitons with spectral parameters $\lambda$ and $\mu$. A similar formula was obtained rigorously in the hydrodynamic limit of a classical hard-rod gasBoldrighini et al. (1983). Despite the mathematical complexities of deriving Eq. (7), whether for classical solitons or hard rods, its physical interpretation in terms of two-body phase shifts is straightforward. This interpretation led El and Kamchatnov to propose an equation of the form (7) for arbitrary soliton gases with two-body elastic scattering, arising from a given classically integrable PDEEl and Kamchatnov (2005). More recently, the system of equations (6), (7) was obtained in the context of Euler-scale hydrodynamics for quantum integrable systems Castro-Alvaredo et al. (2016); Bertini et al. (2016). In particular, one can show that in the quantum setting, the dressed velocities satisfy the integral equationCastro-Alvaredo et al. (2016)

[TABLE]

where $p(\lambda)$ denotes the bare quasiparticle momentum, $e(\lambda)$ the bare energy and $v_{0}(\lambda)=e^{\prime}(\lambda)/p^{\prime}(\lambda)$ the bare group velocity. Also, the differential scattering phase is given in terms of the two-particle S-matrix by

[TABLE]

Upon comparing equations (7) and (8), it is clear that the kinetic theory for quantum solitons with differential scattering phase $\varphi^{\prime}$ may be expressed as a kinetic theory for classical solitons, upon making the identification

[TABLE]

This equivalence was also obtained in Ref. Doyon et al., 2017, independently of the existing theory in Ref. El and Kamchatnov, 2005. Thus, although the original derivation of the equations (6), (8) in the quantum setting made use of an infinite system of Euler equations Castro-Alvaredo et al. (2016); Bertini et al. (2016), these equations are intimately related with a kinetic theory framework for solitons of classically integrable PDEs that is many decades older.

It is worth noting that at present, both the hydrodynamic and kinetic theory viewpoints for quantum integrable models are approximate to the same degree; on the hydrodynamic side, the approximation lies in the “Euler-scale” assumption that diffusive terms are negligible while on the kinetic theory side, the approximation lies in neglecting all higher-order collision terms which could lead to entropy generation. Although it remains to be seen which viewpoint is better suited to incorporating higher-order effects, it appears that all of the proposed extensions of the hydrodynamic formalism for quantum integrable models that take into account new physics, such as external potentials Doyon and Yoshimura (2017) and collision terms Doyon and Spohn (2017) (and indeed the recent analogy with classical solitonsDoyon et al. (2017)), fall naturally into a Boltzmann paradigm rather than a hydrodynamic one.

A natural question to ask is whether there exist gases of solitons arising from classically integrable PDEs whose kinetic theory is captured by the Bethe-Boltzmann equation for some quantum integrable system. For example, there is a well-known mappingIshikawa and Takayama (1980) between hole-type excitations of the Lieb-Liniger gas and dark solitons of the defocusing non-linear Schrödinger equation in the weak-coupling limit, $c\to 0^{+}$. We have found that the position shift resulting from a collision of two dark solitonsFaddeev et al. (2007) coincides with an expression for the inverse of the Lieb-Liniger kernel obtained by Gaudin in the weak-coupling limitGaudin (1971).

II.3 The Bethe-Boltzmann Equation

II.3.1 Formulation

We now return to the precise formulation of the Bethe-Boltzmann equation (4) for the Lieb-Liniger model. Although the integral equation Eq. (7) is helpful for developing the analogy with kinetic theory, it is easier in practice to use the following explicit formula for $v$. Recall that in TBA, the group velocity of dressed excitations on a given equilibrium state with Fermi factors $\{\theta(k):k\in\mathbb{R}\}$ for each pseudo-momentum $k$, is given by

[TABLE]

The manipulations required to derive this equation are summarized in Appendix B. In order to be able to use this result in the Bethe-Boltzmann formalism, we must impose an additional assumption on the local pseudo-momentum density $\rho(x,t,k)$. In particular, we need to assume that each microstate corresponding to the set of local occupation numbers $\{\rho(x,t,k):k\in\mathbb{R}\}$ defines an eigenstate of the Lieb-Liniger Hamiltonian. Thus we demand that the local Bethe equation

[TABLE]

holds 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]

To summarize, the Bethe-Boltzmann equation is shorthand for the hierarchy of equations

[TABLE]

which together comprise a conservation law with self-consistently determined velocity. We can write this schematically as

[TABLE]

This turns out not to be the most useful form of the Bethe-Boltzmann equation. As in the previous works Castro-Alvaredo et al. (2016); Bertini et al. (2016), we find it more convenient to view $\vartheta$ as the fundamental degree of freedom rather than $\rho$. This is possible because the local Bethe equation (12) allows one to express $\rho$ and $\rho v$ as functionals of the local Fermi factor $\vartheta$, namely

[TABLE]

Upon substituting these expressions into the Bethe-Boltzmann equation, a surprising simplification occurs, and a conservation law for $\rho$ is replaced by a simpler advection equation for $\vartheta$. One finds that

[TABLE]

or schematically,

[TABLE]

This is an example of a quasilinear advection equation. In general, the time evolution of such equations rapidly leads to shock formation, even for smooth initial conditions, and without additional assumptions the initial value problem is ill-posed. For two-reservoir initial conditions, this issue was circumvented Bertini et al. (2016); Castro-Alvaredo et al. (2016) by imposing a certain self-consistent ansatz, which for TBA-based hydrodynamics picks out a unique weak solution $\vartheta(x,t,k)$.

Before the present work, the equation (18) had only been verified numerically for long-time steady states of scale-invariant quenches, a regime in which the validity of hydrodynamics can be argued using Bethe-ansatz techniquesCastro-Alvaredo et al. (2016); Bertini et al. (2016). Although the possibility of applications to finite-time dynamics had been raised in previous worksCastro-Alvaredo et al. (2016); Doyon and Yoshimura (2017), there had been no means to actually solve the hydrodynamic equations in the finite-time regime.

In order to address this shortcoming, we have developed a numerical scheme to solve Eq. (18) at finite times for evolution from arbitrary smooth, locally equilibrated initial conditions, whose details are given in Appendix C. The simplest implementation of this scheme yields a self-consistent ansatz that generalizes the two-reservoir case. We were able to obtain physically reasonable predictions for a range of models and initial conditions, which are mostly deferred to a companion paper Bulchandani et al. (2017) except for one example, a thermal expansion in the XXZ model. As discussed in Ref. Castro-Alvaredo et al., 2016, we expect on general physical grounds that Eq. (18) should be accurate at time-scales which are long compared to the time-scale of local thermalization and at length scales which are long compared to the length scale over which the system is locally in thermal equilibrium, even though rigorous justification for applications beyond the scaling limit is lacking at present.

II.4 The Two-Reservoir Quench

Let us now outline the ansatz for non-equilibrium steady states of the two-reservoir quench already presented in Refs. Castro-Alvaredo et al., 2016; Bertini et al., 2016. Thus consider solving the initial value problem

[TABLE]

for $t>0$, with discontinuous initial conditions

[TABLE]

where $\theta_{L}$ and $\theta_{R}$ denote the Fermi factors for initial equilibrium states with temperatures and chemical potentials $\{T_{L},\mu_{L}\}$ and $\{T_{R},\mu_{R}\}$ respectively, as given by the Yang-Yang equations Yang and Yang (1969), and $H$ denotes the Heaviside step function. By analogy with the solution by characteristics for the Burgers equation, one can write down the ansatz

[TABLE]

which yields111This only works for two reservoirs due to (i) scale invariance and (ii) the fact that TBA-based hydrodynamics is linearly degenerate, in the sense of LaxLax (1957). Otherwise, existence and uniqueness of solutions to this ansatz are not guaranteed, even if one assumes monotonicity of $v[\vartheta](k)$ in $k$.

[TABLE]

In the special case that $v[\hat{\vartheta}(x,t)](k)$ is monotonic in $k$, we can write this as a step-function of the wavenumber $k$, as was done in Ref. Castro-Alvaredo et al., 2016. To see this, suppose that for fixed $x,t$, the function $v[\hat{\vartheta}(x,t)](k)$ is monotonic in $k$. Then the equation

[TABLE]

has a unique solution, $k^{*}(x,t)$, such that

[TABLE]

Thus for example, if $v[\hat{\vartheta}(x,t)](k)$ increases with $k$, we can write

[TABLE]

As they stand, equations (22) and (25) both appear intractable. However, when solving for long time steady-states, we can exploit the crucial property of scale-invariance. In particular, at long times we may suppose that $\vartheta(x,t,k)$ depends on position and time via their ratio $\zeta=x/t$ alone. In (22), this yields

[TABLE]

which is essentially equation (16) of Ref. Bertini et al., 2016. Assuming that $v$ increases with $k$, this may be recast as the pair of self-consistent equations

[TABLE]

which is essentially equation (35) of Ref. Castro-Alvaredo et al., 2016. This form is particularly amenable to iterative solution.

II.5 Hydrodynamic Charges and Currents

In order to make contact with direct numerical simulations of the two-temperature quench, we must develop hydrodynamic expressions for local charges and currents of the model. Therefore suppose that $\mathbf{Q}$ is a conserved charge operator of the model, with single-particle eigenvalue $q(k)$. Then the total charge carried by a Bethe wavefunction with limiting density of states $\rho(k)$ is given by

[TABLE]

in the thermodynamic limit. Whilst this is a standard result, the surprising claim of Ref. Bertini et al., 2016 (proven in the case of Lorentz or Galilean-invariant theories in Ref. Castro-Alvaredo et al., 2016) is that the local currents associated with such charges may also be written in terms of the local density of states, via the formula

[TABLE]

where $v(k)$ denotes the quasi-particle velocity. In the hydrodynamic approximation, this allows us to write down expressions for the local charge and current densities associated with the operator $\mathbb{Q}$, namely

[TABLE]

Given these definitions, the Bethe-Boltzmann equation immediately implies local conservation of charge, in the form

[TABLE]

The equation (32), together with a “completeness property” of local conserved charges, is used to derive the Bethe-Boltzmann equation by both Ref. Castro-Alvaredo et al., 2016 and Ref. Bertini et al., 2016. However, attempting to derive the Bethe-Boltzmann equation in this way at finite length- and time- scales ends up being equivalent to the phenomenological derivation we provide in Section II.1, as one must a priori make a local-density-type approximation for $\rho$, whose rigorous justification for non-stationary states does not yet exist.

III Linear Response and Drude Weights in the spin-$1/2$ XXZ Chain

Now that we are equipped with the Bethe-Boltzmann or generalized hydrodynamics framework, we illustrate how it can be applied to study energy and spin transport in the spin-$1/2$ XXZ spin chain, and compare our results to density-matrix renormalization group calculations. Even near equilibrium (linear response), we argue that this framework allows one to compute transport quantities that were previously inaccessible, including the spin Drude weight at arbitrary temperature.

III.1 Bethe-Boltzmann Formalism for the XXZ Chain

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

[TABLE]

Here, 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 assume in the following that $-1<\Delta<1$; the behavior outside this regime is mentioned briefly in the final Discussion. The Bethe-Boltzmann equation for the Hamiltonian (33) is discussed in detail in Ref. Bertini et al., 2016. The derivation proceeds almost exactly as for the Lieb-Liniger gas, except that one must now account for the “strings” of bound states appearing in the thermodynamic limit. We 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, and assume that there are $N_{t}$ string types in total. One can show that their dressed velocities are given by the formula

[TABLE]

whose motivation is sketched in Appendix B. Since there are now multiple branches of quasiparticle excitations, we must postulate multiple Bethe-Boltzmann equations; in abridged form, these read

[TABLE]

with advection formulation

[TABLE]

where the $v_{j}$ are given by  (34). For steady-states of the two-reservoir quench, the latter may be solved formally to yield the $2N_{t}$ coupled equations

[TABLE]

where we again assumed that the $v_{j}$ were monotonic in $\lambda$.

III.2 Linear Response

We now turn to conventional (linear response) quantum transport in integrable systems. Linear response transport coefficients are given by the Kubo formula, which relates the conductivity at zero frequency to the integral over equilibrium dynamical correlation functions describing the return to equilibrium of a spontaneous fluctuation. In an integrable system that does not thermalize in a conventional way, energy or spin (charge) currents may not be able to relax if they have a non-zero overlap with conserved quantities, leading to a divergent zero-frequency conductivity and dissipation-less transport Zotos et al. (1997). The degree of divergence of the DC conductivity may be characterized by considering the conductivity at finite frequency and defining the Drude weight $D(T)$ such that

[TABLE]

with $D(T)$ given by the long-time behavior of the equilibrium dynamical correlation function

[TABLE]

For energy transport in the XXZ spin chain, this is especially simple as the energy current $J_{E}$ is a conserved quantity ($[H,J_{E}]=0$), so that $\sigma(\omega)=\pi D_{E}\sigma(\omega)$ with $D_{E}=\beta\langle J_{E}^{2}\rangle/L$. In this case, the Drude weight can be computed using standard Bethe Ansatz techniques Klümper and Sakai (2002); Zotos (2016). In general, the Drude weight for a current is bounded from below by a sum over conserved quantities that have a nonzero projection on to the current, via the Mazur inequality Mazur (1969); Zotos et al. (1997)

[TABLE]

where $Q_{\alpha}$ are independent local conserved quantities Suzuki (1971). In the XXZ spin chain at zero magnetic field, the conventional (strictly local) conserved quantities give zero contribution to the Drude weight by symmetry Sirker et al. (2009), but a new set of conserved quantities Prosen (2011) (see also Prosen and Ilievski (2013); Prosen (2014); Pereira et al. (2014)) that are given by sums of quasi-local operators (local up to exponential tails) do contribute. At least at high temperatures and some values of anisotropy $\Delta$, these new integrals of motion appear to saturate the numerical value of the Drude weight Karrasch et al. (2012) obtained from time-dependent density-matrix renormalization group simulations. Two different thermodynamic Bethe ansatz expressions for the spin Drude weight have been proposed Zotos (1999); Benz et al. (2005), yielding contradictory results. These Bethe ansatz results are controversial and they were argued to violate exact results at high temperatures (see Ref. Sirker et al., 2011 and references therein). Going beyond linear response, the description of non-equilibrium spin transport in the XXZ model is described by generalized Gibbs ensembles that include non-standard quasi-local conserved quantities De Luca et al. (2017).

III.3 Spin Drude weights from hydrodynamics

Given the long history of linear-response spin transport in the XXZ spin chain, it seems intriguing that the hydrodynamic approach introduced above could describe non-equilibrium spin transport exactly. In particular, it is interesting to note that the hydrodynamic approach discussed above could have been discovered, in principle, shortly after the development of the TBA formalism Yang and Yang (1969) — years before the subtleties related to quasi-local charges in the XXZ chain were uncovered by Prosen Prosen (2011). Nevertheless, we will argue below that the Bethe-Boltzmann equation (4) can be used to compute both energy and spin Drude weights in agreement with exact low and high temperature results, and with density-matrix renormalization group (DMRG) calculations.

We start by considering energy transport for a two-temperature quench with a left reservoir at temperature $T_{L}$ and a right reservoir at temperature $T_{R}$, joined together at time $t=0$. For a small temperature difference $T_{L}=T+\Delta T/2$ and $T_{R}=T-\Delta T/2$ ($\Delta T\ll T$), it is natural to expect that the energy current in the steady state should be described by linear-response theory. In fact, because the energy current is itself a conserved quantity, one can show that the spatial integral of the energy current at long times is determined by the equilibrium Drude weight even far from equilibrium Vasseur et al. (2015)

[TABLE]

for arbitrary $T_{L}$ and $T_{R}$. In particular, for a small temperature gradient $\Delta T\ll T$, the Drude weight immediately follows from the value of the energy current in the steady state $D_{E}=\lim_{\Delta T\to 0}\frac{1}{\Delta T}\int d\xi j_{E}(\xi=x/t)$. Interestingly, equation (41) is exactly satisfied by the Bethe-Boltzmann hydrodynamic approach: in a way similar to the Lieb-Liniger discussion above (Sec. II.5), one can show Bertini et al. (2016) that the expectation value of the local energy current $\langle j_{E}(x,t)\rangle$ from the hydrodynamic framework coincides with the local conserved quantity $\langle q_{3}(x,t)\rangle$, as it should for the XXZ spin chain where the energy current is a conserved quantity. This is a non-trivial check on the hydrodynamic approach. Integrating spatially the conservation law $\partial_{t}q_{3}+\partial_{x}G=0$, with $G$ the current associated with the conserved charge $q_{3}$, thus yields $\int dx\,j_{E}/t=G(T_{L})-G(T_{R})$, where the “state function” $G$ can be determined for a small temperature gradient from linear response Vasseur et al. (2015) . One finds $G(T)=\int^{T}dTD_{E}(T)$, from which Eq. (41) follows.

The hydrodynamic approach is therefore consistent with the linear-response energy Drude weight. This is perhaps not especially surprising, given that the energy Drude weight has a very simple equilibrium expression that is accessible from quantum transfer matrix or thermodynamic Bethe ansatz Klümper and Sakai (2002); Zotos (2016). The case of spin transport between two reservoirs prepared at different magnetic fields $h_{L}$, $h_{R}$ and uniform temperature $\beta^{-1}$ is much more interesting. In this case, there is no simple relation like (41) relating non-equilibrium transport and equilibrium Drude weights, but for small fields $h_{L}=-h_{R}=\frac{h}{2}\ll J=1$, the spin Drude weight can be expressed as (see Refs. Vasseur et al., 2015; Karrasch, 2017)

[TABLE]

which is a spatial integral over the steady-state spin current. This formula allows us to extract the spin Drude weight by solving (37) iteratively, and using the hydrodynamic formula

[TABLE]

for the steady-state spin current in a two-reservoir quench with $T_{L}=T_{R}$ and $h_{R}=-h_{L}=\frac{h}{2}\ll 1$.

Results are shown in Fig. 1 for different values of the XXZ anisotropy parameter $\Delta$. The small field gradient was taken to be either $h=10^{-3}$ or $h=10^{-4}$ with no significant difference in the value of the Drude weight, and we carefully checked for $\Delta=\frac{1}{2}$ that all the numerical errors associated with the numerical discretization of the hydrodynamic equations lead to relative errors on the Drude weight below $10^{-4}$. Our results are in good agreement with exact asymptotic results at low Shastry and Sutherland (1990) and high Prosen and Ilievski (2013) temperature – note that the high temperature results are strictly speaking a lower bound that is believed to be saturated for the values of $\Delta$ that we consider Karrasch et al. (2012); Prosen and Ilievski (2013). The speed of convergence of the spin Drude weight to the $T=0$ result decreases as $\Delta$ approaches $\Delta=1$, consistent with earlier results Zotos (1999).

We further checked these results against density-matrix renormalization group (DMRG) calculations. To this end, we used a finite-temperature version Karrasch et al. (2012); Barthel (2013); Kennes and Karrasch (2016) of the real-time DMRG White (1992); Schollwoeck (2011) to compute the linear response current-correlation function $\langle J(t)J(0)\rangle_{\beta}$ that appears in Eq. (39). The key parameter governing the accuracy of this method is the so-called discarded weight, which we chose such that the error of $\langle J(t)J(0)\rangle_{\beta}$ was smaller than the linewidth. The system size was taken large enough for all results to be effectively in the thermodynamic limit (a typical value is $L\sim 200$).

The results of the hydrodynamic approach turn out to match the DMRG predictions to within numerical accuracy (see Fig. 2 and Appendix A), except at low temperatures where extracting the Drude weight from DMRG calculations would require accessing longer time scales. Thus it seems that the non-equilibrium steady-state predicted by hydrodynamics is cognizant of the quasi-local conserved charges discovered in Refs. Prosen, 2011; Prosen and Ilievski, 2013; Prosen, 2014; Pereira et al., 2014 (see also De Luca et al., 2017).

It is also instructive to discuss the Drude weight for negative values of $\Delta$, for which bosonizationSirker et al. (2011) predicts a small decay rate at low temperatures and finite-time data need to be interpreted with caution. In this regime, hydrodynamics can be viewed as a check on results obtained from the DMRG. In Fig. 3, we show that the DMRG curves saturate to the hydrodynamics prediction even at low $T$. This shows that almost all of the spectral weight is concentrated in the Drude peak and that the DMRG data for the Drude weight Karrasch et al. (2013b) is not misleading.

IV Discussion and a strongly non-equilibrium example

Having treated linear response in the XXZ chain in the previous section, we now consider the following far-from-equilibrium example of the implementation described in Section II. Thus suppose that a system is prepared with a Gaussian temperature profile

[TABLE]

where $L$ characterizes the spread of the initial localized temperature inhomogeneity. We compared the results of the one-step implicit hydrodynamic solution (70) to DMRG simulations starting from the same initial profile (44) for various values of $\beta_{0}$, $\beta_{M}$ and $\Delta$ and found an excellent agreement (Fig. 4).

The agreement is perhaps surprising given the way that Eq. (70) is motivated in Appendix C, as an implementation of a finite-difference scheme with only one time step. We do not yet understand this level of agreement. It could arise because there are a number of limits or simple cases for which the finite-difference scheme is actually exact, which constrains the behavior even in other cases such as the Gaussian initial profile. It is also possible that the scheme is accurate to higher order. We will discuss more examples and mathematical aspects of this approach in a future publication Bulchandani et al. (2017).

A number of natural questions arise for future work. The first concerns the limitations of the hydrodynamical approach. For example, the Bethe-Boltzmann equation is effectively classical in the sense that nearly the same mathematical structure would arise in the dynamics of an integrable classical particle gas. (For a recent study of the two-reservoir quench in classical systems and references to earlier classical work, see Ref. Mendl and Spohn, 2016.) Where did quantum-mechanical effects, reflecting the wave nature of the particles, go? Obviously the quantum-mechanical interactions between particles determine the phase shifts that underlie the Bethe ansatz, but it is apparently true that long-time dynamics in many situations is effectively classical. This includes situations such as the two-reservoir quench that are sufficient to determine the Drude weights. It should be possible to discern quantum effects in short-time or short-length behavior, which could be viewed as transients before the local GGE assumption of generalized hydrodynamics becomes justified. We should also note that there is no reason to expect the striking agreement in Fig. 4 for observables such as spin density, which has two transport channels (ballistic and diffusive) even in the gapless regime Sirker et al. (2011). This is in contrast to energy density, whose transport is purely ballistic across all regimes, and is related to the fact that the energy current is a conserved charge of the XXZ Hamiltonian whereas the spin current is not.

Aside from numerical studies, a complementary approach that might be useful for understanding the scope of hydrodynamical techniques is based on exact solutions for time evolution in certain limits, such as the Luttinger-liquid two-reservoir quench studied in recent work Langmann et al. (2017). This should be comparable to the XXZ model studied in this paper in the low-temperature regime. At intermediate times, features are seen in the time evolution of densities that are compared to those in numerical calculations Karrasch et al. (2013a) and may reflect finite-time corrections to the hydrodynamic description. For the case of energy transport, similar terms, also involving the Schwarzian derivative of the initial temperature distribution, appeared previously in a calculation based on conformal invariance Sotiriadis and Cardy (2008).

Another obvious question concerns the mathematical existence and physical validity of the hydrodynamical solutions in Section II.3, beyond the special case of the two-reservoir quench. The two-reservoir quench is quite special for a number of reasons: for example, it has no intrinsic time or length scale, which means that the scaling limit is effectively a complete description of its universal properties. We have found that for at least some cases of practical importance, computations based on the implicit method described here yield physically plausible results, even after the time at which shocks from different initial discontinuities coincide. However, this is very far from a mathematical demonstration of existence of solutions, which may be unreasonable to expect given that no such proof exists even for the venerable equations of standard hydrodynamics.

A deeper physical understanding of the behavior of the Bethe-Boltzmann equation in various important contexts is a more feasible goal. Even for the two-reservoir quench, important questions remain. We have limited ourselves in this paper to the regime $-1<\Delta<1$, when the dynamics have a ballistic component (e.g., the energy and spin Drude weights are nonzero). A very recent numerical study Ljubotina et al. (2017) finds diffusive behavior for $|\Delta|>1$ (see also Refs. Langer et al., 2009; Steinigeweg and Brenig, 2011; Žnidarič, 2011; Karrasch et al., 2014) and superdiffusive behavior at the Heisenberg points $|\Delta|=1$, and it would be desirable yet difficult to capture this behavior using the hydrodynamic formalism. While it is no doubt challenging to capture the entire diversity of dynamical behavior in integrable models within a single formalism, the Bethe-Boltzmann equation, or equivalently generalized hydrodynamics, is at an exciting stage of its development, with important results for some long-standing problems and tantalizing hints for others.

Acknowledgments. We thank Fabian Essler for several valuable discussions. We also wish to thank J. De Nardis, B. Doyon, M. Fagotti, and E. Ilievski for comments on the manuscript. J.E.M. and V.B.B. are grateful to the Institut d’Études Scientifiques de Cargèse and the organizers of the conference “Exact Methods in Low-Dimensional Statistical Physics 2017” for their hospitality. We acknowledge support from the Department of Energy through the LDRD program of LBNL (R.V.), from NSF DMR-1507141 and a Simons Investigatorship (J.E.M.), and center support from the Moore Foundation’s EPiQS initiative. C.K. is supported by the DFG via the Emmy-Noether program under KA 3360/2-1.

Note added. While this work was being completed, we became aware of a related paper Ilievski and De Nardis (2017), which also showed that the spin Drude weight could be obtained from hydrodynamics.

Appendix A Spin Drude weight from DMRG

In the main text, we compared DMRG calculations of the spin Drude weight to the hydrodynamic approach for $\Delta=\frac{1}{2}$. We also performed detailed comparisons for $\Delta=\frac{1}{\sqrt{2}}$ (left) and $\Delta=\cos\frac{\pi}{5}$ with good agreement, even though the convergence of the DMRG calculations to the asymptotic values becomes slower as $\Delta$ gets closer to the isotropic value $\Delta=1$.

Appendix B Summary of Thermodynamic Bethe Ansatz

B.1 Lieb-Liniger Gas

In this section, we briefly recall the main features of thermodynamic Bethe ansatz for the Lieb-Liniger gas. Recall that the Lieb-Liniger Hamiltonian for $N$ bosons on a line, with interaction strength $c$, may be written as

[TABLE]

Away from the collision planes $\{x_{j}=x_{k}:j\neq k\}$, the eigenfunctions of this Hamiltonian are superpositions of plane wavefunctions. For example, in the fundamental chamber $\mathcal{D}_{0}=\{x_{1}<x_{2}<\ldots<x_{N}\}$ we may write

[TABLE]

where the sum is over all permutations on $N$ letters. The model is called integrable because all $N$-body scattering processes factorize into two-body processes. This property turns out to be sufficiently restrictive to solve the model. In particular, upon imposing periodic boundary conditions

[TABLE]

one may use certain combinatorial relations between the amplitudes $A(P)$ to deduce the Bethe equations

[TABLE]

where $\varphi(k)$ denotes the two-body phase shift. The Bethe quantum numbers $I_{a}$ are integers for fermions or an odd number of bosons, and half-integers for an even number of bosons. Generally one finds that only certain allowed Bethe quantum numbers are occupied in a given eigenstate of the Hamiltonian. Passing to the thermodynamic limit $N,L\to\infty$, with $N/L$ fixed, this allows one to define densities of particles, holes and vacancies via

[TABLE]

respectively. By definition, we have

[TABLE]

and it is useful to define the non-equilibrium Fermi factor, given by

[TABLE]

One can then show Takahashi (1999) that the continuum version of the Bethe equations  (48) reads

[TABLE]

where the kernel $\mathcal{K}$ is given by

[TABLE]

Integral equations of this type occur sufficiently frequently in the non-equilibrium theory that we define operators $\hat{K}$ and $\hat{\vartheta}$ which act on functions via

[TABLE]

For example, in operator notation, equation  (52) reads

[TABLE]

Now consider particle-type excitations on such states. These arise when one introduces an additional quantum number in the Bethe equations (48) causing a shift in wavenumber $k_{j}\to k_{j}^{\prime}$ across all $k$-space, which reflects the collective nature of the underlying excitation. The physics of such excitations is captured by the pair of integral equations,

[TABLE]

where $F(k^{\prime}|k)\mathrel{\mathop{:}}=L\rho^{t}(k)\Delta k$ denotes the backflow function222Our convention for the backflow function follows Ref. Bertini et al., 2016, and in that respect is different from the convention $F(k^{\prime}|k)=L\rho(k)\Delta k$ used, for example, in Ref. Takahashi, 1999. for a particle-type excitation with momentum $k$, and $Q$ and $\Delta Q$ denote the bare and dressed values respectively for any given conserved charge of the model. Rather conveniently, the equations (57) together imply the closed integral equation

[TABLE]

for the dressed charge. In the context of equilibrium TBA, this equation is used to justify the interpretation of $\vartheta$ as a Fermi factor Korepin et al. (1997); Takahashi (1999). In the present context, it allows us to determine the dressed charges carried by a quasiparticle excitation directly from the bare charges. Now recall that the group velocity of a quasiparticle excitation is given by

[TABLE]

where $\Delta E$ and $\Delta P$ denote its dressed energy and momentum respectively. From the bare values $P(k)=k,\,E(k)=k^{2}/2$, we can take the formal inverse333This is shorthand for the infinite (Neumann) series obtained by solving  (58) iteratively. of the kernel appearing in  (58) to yield the dressed values

[TABLE]

so that the quasiparticle velocity is given by

[TABLE]

B.2 The Gapless XXZ Chain

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 that the model is in its gapless phase, i.e. $-1<\Delta<1$. The thermodynamic limit of this phase is believed to be characterized by a finite number of “strings” of bound states, consisting of families of rapidities in the complex plane possessing the same real part Takahashi (1999). We 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. 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]

The TBA equations then read

[TABLE]

where $j\in\{1,2,\ldots,N_{t}\}$; for definitions of the various terms see the book of Takahashi Takahashi (1999). We define quantities

[TABLE]

and the Fermi factors

[TABLE]

for strings of type $j$. 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]

These imply the formula

[TABLE]

for the velocity of quasiparticle excitations within each string type.

Appendix C Finite-time Numerical Scheme

Consider the infinite-dimensional initial value problem

[TABLE]

Motivated by analogy with the solution by characteristics for the Burgers equation, we propose the compact approximate ansatz

[TABLE]

Here the velocity functional $v[\hat{\vartheta}]$ is evaluated using all values of $\vartheta$ at $(x,t)$. This ansatz is presumably not exact in general but arises naturally as the one-step version of a backwards implicit scheme to solve the advection equation Eq. (69). 444We thank B. Doyon for helpful discussions and for sharing unpublished work on a different approach to this problem. We solve it numerically for an example in the closing discussion section of this paper (see Fig. 4 above).

First, in order to provide some intuition, we illustrate the meaning of this ansatz for the finite-dimensional version of Eq. (69) (as arises in practice when discretizing on a computer). Thus consider the initial-value problem

[TABLE]

for $t>0$, with $n=1,2,\ldots,N$. We will have found a solution of this equation (ignoring the possibility of discontinuities in any derivative) if we can find $\theta_{n}$ satisfying

[TABLE]

for every $x_{0},t_{0}$ and small $dt$, because expanding the left-hand side in a Taylor series yields the advection equation. Since the error is of order $dt^{2}$, this could be the basis of a convergent finite-difference scheme: given $\theta_{n}$ at time $t_{0}$, one steps forward in time repeatedly by $dt$ using Eq. (72), and the global error in advancing by time $t$ is small as $dt\rightarrow 0$. However, it is well-known from the theory of simpler differential equations that this type of Euler method can be quite unstable. A safer option is to use an implicit or backwards method Press (2007): for a time step $dt$, we solve for the current values of $\theta_{n}$ by looking back at earlier times.

[TABLE]

where the notation means that the velocities are evaluated at ${x,t}$. This is just the infinitesimal form of Eq. (70); in other words, Eq. (70) is obtained by covering the entire desired integration range of time in a single step. Another advantage in practice of the backwards scheme is that, unlike for the forward scheme, it is easy to get all the $\theta_{n}$ at a particular $(x,t)$ point without interpolation if the initial condition is known everywhere.

We have found that one can solve the ansatz (73) by numerical iteration to obtain detailed predictions for finite-time evolution, from a range of initial conditions Bulchandani et al. (2017). This is of particular experimental relevance for the Lieb-Liniger model, which provides a physically realistic description of quasi one-dimensional Bose gases. We report one such calculation for energy expansion in the XXZ chain at the end of this paper using the simpler approximate ansatz (70) (corresponding to a single time step with $dt=t$ in eq. (73)), after introducing that model and discussing its linear response, and discuss there its (perhaps unexpected) agreement with a DMRG solution. A detailed discussion of the range of validity of the predictions obtained using this method is given in Sections II and IV of the main text.