Ballistic transport in the one-dimensional Hubbard model: the hydrodynamic approach

TL;DR

This paper develops a hydrodynamic formalism for quantum integrable systems, specifically the one-dimensional Hubbard model, providing exact expressions for transport coefficients and revealing fundamental symmetries in non-equilibrium states.

Contribution

It introduces a universal hydrodynamic framework for multi-species quantum systems and derives the first exact analytic formulas for Drude weights in the Hubbard model.

Findings

Exact formulas for Drude weights in the Hubbard model

Hydrodynamic scale expressions for charge susceptibilities and correlators

Demonstration of Onsager relations in non-canonical states

Abstract

We outline a general formalism of hydrodynamics for quantum systems with multiple particle species which undergo completely elastic scattering. In the thermodynamic limit, the complete kinematic data of the problem consists of the particle content, the dispersion relations, and a universal dressing transformation which accounts for interparticle interactions. We consider quantum integrable models and we focus on the one-dimensional fermionic Hubbard model. By linearizing hydrodynamic equations, we provide exact closed-form expressions for Drude weights, generalized static charge susceptibilities and charge-current correlators valid on hydrodynamic scale, represented as integral kernels operating diagonally in the space of mode numbers of thermodynamic excitations. We find that, on hydrodynamic scales, Drude weights manifestly display Onsager reciprocal relations even for generic (i.e. non-canonical) equilibrium states, and establish a generalized detailed balance condition for a general quantum integrable model. We present the first exact analytic expressions for the general Drude weights in the Hubbard model, and explain how to reconcile different approaches for computing Drude weights from the previous literature.

Equations

$${\epsilon }_a^{\prime }={\Omega }_{ab}\star e_b^{\prime },p_a^{\prime }={\Omega }_{ab}\star k_b^{\prime }.$$

$${\left({\Omega }^{-1}\right)}_{ab}(u,w)={\delta }_{ab}\delta (u-w)+K_{ab}(u-w){\vartheta }_b(w){\sigma }_b.$$

$$q_{a,i}^{\mathrm{eff}}={\Omega }_{ab}\star q_{b,i}={\partial }_{{\mu }_i}\log \left({\vartheta }_a^{-1}-1\right).$$

$${\mathcal{D}}^{(i,j)}=\frac{\beta }{2}\underset{t\to \infty }{\lim }{\int }_{\tau =0}^t\mathrm{d}\tau ⟨{\hat{J}}_i(\tau ){\hat{j}}_j(0){⟩}_{\mathrm{c}},$$

$${\mathcal{A}}_{ij}=\sum _a\int \mathrm{d}u{q}_{a,i}^{\mathrm{eff}}(u)A_a(u)q_{a,j}^{\mathrm{eff}}(u),$$

$${\mathcal{D}}^{(i,j)}=\frac{\beta }{2}\underset{\delta {\mu }_j\to 0}{\lim }\frac{\partial }{\partial \delta {\mu }_j}\underset{t\to \infty }{\lim }\frac{J_i(t)}{t}.$$

$$D_a(u)={\rho }_a(u)(1-{\vartheta }_a(u)){\left(v_a^{\mathrm{dr}}(u)\right)}^2.$$

$${\mathcal{S}}_{\hat{q},a}(\kappa ,-\omega )=e^{-{\mathcal{F}}_a(k,\omega )}{\mathcal{S}}_{\hat{q},a}(\kappa ,\omega )+O({\kappa }^2)$$

$$\hat{H}=\sum _{x=1}^L{\hat{T}}_{x,x+1}+4\mathfrak{u}\sum _{x=1}^L{\hat{V}}_{x,x+1},$$

$$\hat{H}_{0}=-\sum_{x=1}^{L}\sum_{\sigma={\uparrow},{\downarrow}}\Big{(}\hat{c}^{\dagger}_{x,\sigma}\hat{c}_{x+1,\sigma}+\hat{c}^{\dagger}_{x+1,\sigma}\hat{c}_{x,\sigma}\Big{)}+4{\mathfrak{u}}\sum_{x=1}^{L}(\hat{n}_{x,{\uparrow}}-\tfrac{1}{2})(\hat{n}_{x,{\downarrow}}-\tfrac{1}{2}),$$

$$\hat{N}=\sum_{x=1}^{L}\Big{(}\hat{c}^{\dagger}_{x,{\uparrow}}\hat{c}_{x,{\uparrow}}+\hat{c}^{\dagger}_{x,{\downarrow}}\hat{c}_{x,{\downarrow}}\Big{)},\qquad\hat{S}^{z}=\frac{1}{2}\sum_{x=1}^{L}\Big{(}\hat{c}^{\dagger}_{x,{\uparrow}}\hat{c}_{x,{\uparrow}}-\hat{c}^{\dagger}_{x,{\downarrow}}\hat{c}_{x,{\downarrow}}\Big{)}.$$

$$e^{\mathrm{i}k(u_k)L}\prod _{j=1}^M\frac{u_k-w_j-\mathrm{i}\mathfrak{u}}{u_k-w_j+\mathrm{i}\mathfrak{u}}$$

$$\prod _{j=1}^N\frac{w_k-v_j-\mathrm{i}\mathfrak{u}}{w_k-u_j+\mathrm{i}\mathfrak{u}}\prod _{m=1}^M\frac{w_l-w_m+2\mathrm{i}\mathfrak{u}}{w_l-w_m-2\mathrm{i}\mathfrak{u}}$$

$$u_i=\frac{1}{2}\left(y_i+\frac{1}{y_i}\right).$$

$$x(u)=u+u\sqrt{1-\frac{1}{u^2}},$$

$$y_{+}(u)=x(u),y_{-}(u)=\frac{1}{x(u)}.$$

$$y_{\pm }(u)=x(u\pm \mathrm{i}0),$$

$$e^{\mathrm{i}k(u)}=\mathrm{i}y,$$

$$k_{\pm }(u)=-\mathrm{i}\log (\mathrm{i}x(u\pm \mathrm{i}0)).$$

$$k_{-}(u)$$

$$k_{+}(u)$$

$$k_{\pm }^{\prime }(u)=\pm \frac{1}{\mathrm{i}(u\pm \mathrm{i}0)\sqrt{1-1/(u\pm \mathrm{i}0{)}^2}}=\mp \frac{1}{\sqrt{1-u^2}},u\in (-1,1).$$

$$e_{\pm }(u)=-2(\cos (k_{\pm }(u))+\mathfrak{u})=\pm 2\sqrt{1-u^2}-2\mathfrak{u}.$$

$$u_{\pm j}=u\pm (M+2-2j)\mathrm{i}\mathfrak{u},w_j=u+(M+1-2j)\mathrm{i}\mathfrak{u},j=1,2,\dots ,M.$$

$$k_{M|uw}(u)$$

$$e_{M|uw}(u)$$

$$=2\sqrt{1-(u+M\mathrm{i}\mathfrak{u}{)}^2}+2\sqrt{1-(u-M\mathrm{i}\mathfrak{u}{)}^2}-4M\mathfrak{u},.$$

$$k_{M|uw}^{\prime }(u)=k_{+}^{\prime }(u+M\mathrm{i}\mathfrak{u})+k_{-}^{\prime }(u-M\mathrm{i}\mathfrak{u})=-\frac{1}{\sqrt{1-(u+M\mathrm{i}\mathfrak{u}{)}^2}}-\frac{1}{\sqrt{1-(u-M\mathrm{i}\mathfrak{u}{)}^2}}.$$

$$w_j=u+(M+1-2j)\mathrm{i}\mathfrak{u},j=1,2,\dots M.$$

$$N=N_{+}+N_{-}+\sum _{M=1}^{\infty }2M{N}_{M|uw},M=\sum _{M=1}^{\infty }M(N_{M|uw}+N_{M|w}).$$

$${\rho }_a\left(u_j^{(a)}\right)=\underset{L\to \infty }{\lim }\frac{1}{L\left(u_{j+1}^{(a)}-u_j^{(a)}\right)}.$$

$$(g\star h)(u,w)\equiv {\int }_{-\infty }^{\infty }\mathrm{d}zg(u,z)h(z,w),$$

Full text

In past few years, a lot of interest has been devoted to studying various paradigms of non-ergodic many-body physics, such as quantum quenches, equilibration to generalized Gibbs ensembles and phenomenon of pre-thermalization Polkovnikov et al. (2011); Calabrese et al. (2016); Gogolin and Eisert (2016). One of the prominent recent results is the formalism of generalized hydrodynamics developed in Bertini et al. (2016); Castro-Alvaredo et al. (2016), with a large number of subsequent studies investigating its various aspects and applications Doyon and Yoshimura (2017); B. Doyon, and H. Spohn, and T. Yoshimura (2017); B. Doyon and H. Spohn (2017); B. Doyon, T. Yoshimura, and J.-S. Caux (2017); B. Doyon, J. Dubail, R. Konik, and T. Yoshimura (2017); V. B. Bulchandani, R. Vasseur, C. Karrasch, and J. E. Moore (2017); Piroli, Lorenzo and Nardis, Jacopo De and Collura, Mario and Bertini, Bruno and Fagotti, Maurizio (2017); Alba, Vincenzo (2017), including the exact computation of Drude weights in the Heisenberg model XXZ spin-1/2 chain Ilievski and De Nardis (2017). In analogy to the conventional theory of hydrodynamics Spohn (1991), the authors of Doyon and Spohn (2017) just recently obtained a closed formula for Drude weights expressed in terms of local equilibrium state functions for the case of integrable Bose gas (Lieb–Liniger model) and conjectured that similar formulae may hold in quantum integrable models more generally. In this work, we go a step further and extend the formalism to integrable models which possess physical particles with internal degrees of freedom and are solvable by nested Bethe Ansatz. Nesting is referred to the situation when physical degrees of freedom are associated with a higher rank symmetry group, leading to eigenfunctions with a hierarchical structure of internal quantum numbers and elementary excitations of different flavours. While studies of such models has been traditionally focused on Gibbs equilibrium Essler et al. (1991, 1992); Essler and Korepin (1994); Fabian H. L. Essler and Holger Frahm and Frank Göhmann and Andreas Klümper and Vladimir E. Korepin (2005); Beisert et al. (2011); Frolov and Quinn (2012), they have also been recently studied in the nonequilibrium context Mestyan, Marton and Bertini, Bruno and Piroli, Lorenzo and Calabrese, Pasquale (2017); Robinson, Neil J. and Caux, Jean-Sébastien and Konik, M. Robert (2016).

The chief aspect in which interacting quantum integrable theories differ from widely studied noninteracting systems is the dressing of (quasi)particle excitations, i.e. a process in which bare properties of the particle-hole type of excitations renormalize in the presence of interactions with a non-trivial reference (vacuum) state. The task of classifying excitations has been traditionally restricted to ground states for some of the simplest Bethe Ansatz solvable models Korepin et al. (1993), and subsequently extended to some important examples of exactly solvable models of correlated electrons Bares et al. (1991); Essler and Korepin (1994); Fabian H. L. Essler and Holger Frahm and Frank Göhmann and Andreas Klümper and Vladimir E. Korepin (2005); Quinn and Frolov (2013). A comprehensive exposition of the dressing formalism for grand canonical ensembles in nested Bethe Ansatz models can be found in Quinn and Frolov (2013).

Dressing formalism.

Integrable theories exhibit a completely elastic (factorizable) scattering of particle-like excitations Zamolodchikov and Zamolodchikov (1979). Properties of such excitations represent the kinematic data of the theory. In particular, in Bethe Ansatz solvable models (see e.g. Korepin et al. (1993); Fabian H. L. Essler and Holger Frahm and Frank Göhmann and Andreas Klümper and Vladimir E. Korepin (2005)) thermodynamic excitations relative to a bare vacuum 111In fermionic interacting integrable models there exist distinct inequivalent possibilities of choosing a bare vacuum. Despite this results in different sets of excitations, various choices have no effect on the values of physical observables. can be inferred from the solutions to (nested) Bethe equations. The latter in a finite volume take the form $e^{{\rm i}p_{\alpha}(u^{(\alpha)}_{k})}\prod_{\beta}\prod_{j=1}^{N_{\beta}}S_{\alpha\beta}(u^{(\alpha)}_{k},u^{(\beta)}_{j})=1$, imposing single-valuedness of many-body eigenstates. Here the sets of quantum numbers $\{u^{(\alpha)}_{k}\}$ are called the Bethe roots and represent rapidity variables for distinct species (or flavours) of elementary excitations. The number and types of excitations depends on the model and can be inferred with aid of representation theory of the underlying quantized Lie (super)algebra. Elementary excitations typically form complexes which are interpreted as bound states. The emergent thermodynamic particle content, which can be inferred by e.g. analysing the $L\to\infty$ limit of Bethe equations, is generally different from elementary excitations and is labeled by a pair of mode numbers, a particle type index $a$ and a real rapidity variable $u$. The complete kinematic data are obtained from the bare momenta $k_{a}(u)$ and energies $e_{a}(u)$, and interparticle scattering phase shifts $\phi_{ab}(u,w)$. Once given these functions, no explicit operator representation of the Hamiltonian and its conservation laws is ever required. In this work we present the details of the entire formalism for the non-trivial case of the (fermionic) Hubbard model.

A distinguished feature of integrable systems is a macroscopic number of local conservation laws which can be formally expressed in terms of a discrete basis of local charges $Q_{i}=\sum_{x}\,q_{i}(x)$, with $x$ labelling lattice sites. The associated currents $J_{i}=\sum_{x}\,j_{i}(x)$ are defined with aid of the continuity equation, $\partial_{t}\hat{Q}_{i}+\partial_{x}\hat{J}_{i}=0$. The key concept of the hydrodynamic approach is the dressing of bare energies $e_{a}\mapsto\varepsilon_{a}$ and momenta $k_{a}\mapsto p_{a}$ of particle excitations, which can be presented in a compact form

[TABLE]

with convolution $(\Omega_{ab}\star f_{b})(u)=\sum_{b}\int{\rm d}w\,\Omega_{ab}(u,w)f_{b}(w)$. In interacting quantum integrable models solvable by (nested) Bethe Ansatz, the matrix convolution kernel $\Omega$ takes a universal form

[TABLE]

with kernels $K_{ab}(u,w)$ defined as derivatives of the scattering phase shifts $\phi_{ab}(u,w)=\phi_{ab}(u-w)$, $K_{ab}(u)=\tfrac{1}{2\pi{\rm i}}\partial_{u}\phi_{ab}(u)$, and $\sigma_{a}={\rm sign}(k^{\prime}_{a}(u))$. The (Fermi) filling functions $\vartheta_{a}(u)$ specify the fraction of occupied modes with rapidities inside a small interval around $u$.

Dispersion relations of excitations $\varepsilon_{a}(u)$ depend on a many-body vacuum which is uniquely specified by the rapidity distributions $\rho_{a}(u)$. In terms of (thermodynamic) particle excitations, the equilibrium averages of charge and current densities decompose as $q_{i}=\sum_{a}\int{\rm d}u\,q_{i,a}(u)\rho_{a}(u)$, $j_{i}=\sum_{a}\int{\rm d}u\,q_{i,a}(u)j_{a}(u)$, where $j_{a}(u)=\rho_{a}(u)v^{\rm dr}_{a}(u)$ are the current densities per mode Bertini et al. (2016); Castro-Alvaredo et al. (2016). The group velocities of propagating particles are thus state-dependent, $v^{\rm dr}_{a}(u)=\varepsilon^{\prime}_{a}(u)/p^{\prime}_{a}(u)$.

We furthermore introduce the effective charges as the bare charges renormalized under transformation $\Omega$, namely the effective value of a local charge density $q_{i}$ is obtained as

[TABLE]

Here parameters $\mu_{i}$ are the chemical potentials of a local (generalized) equilibrium ensemble parametrized as $\hat{\varrho}\simeq\exp{(-\sum_{i}\mu_{i}\hat{Q}_{i})}$ Vidmar and Rigol (2016); Essler and Fagotti (2016); Ilievski et al. (2015). It is important to emphasize that despite the derivatives of dressed energies satisfying $\varepsilon^{\prime}_{a}=\Omega_{ab}\star\,e^{\prime}_{b}=(e^{\prime}_{a})^{\rm eff}$, the effective charges are not the proper dressed charges associated with an excitation, and specifically $\varepsilon_{a}\neq e_{a}^{\rm eff}$. We moreover note that with aid of fusion identities among the scattering kernels, the transformation (2) can be decoupled to a quasi-local form in the mode space, cf. Supplemental Material SM (SM) for explicit form for Hubbard model.

Drude weights.

In this work, we shall mainly be concerned with general off-diagonal Drude weights

[TABLE]

which represent magnitudes of the singular parts of the zero-frequency generalized conductivities Kubo (1957); Mahan (2000), ${\rm Re}\,\sigma_{ij}(\omega)=2\pi\,\mathcal{D}^{(i,j)}\delta(\omega)+\sigma^{\rm reg}_{ij}(\omega)$. We use $\langle\cdot\rangle_{\rm c}$ to denote the connected part of the equilibrium expectation values. Although we shall restrict ourselves to grand canonical equilibria, our formalism applies (without modifications) to general local equilibrium states.

An exact representation for $\mathcal{D}^{(i,j)}$ can be given in terms of the static covariance matrix $\mathcal{C}$, $\mathcal{C}_{ij}=\langle Q_{i}q_{j}\rangle_{\rm c}$, with diagonal components $\chi_{i}=\mathcal{C}_{ii}$ representing (generalized) static susceptibilities, and charge-current correlators (overlaps) $\mathcal{O}$, $\mathcal{O}_{ij}=\langle Q_{i}j_{j}\rangle_{\rm c}$. Explicit expressions in terms of thermodynamic state functions can be found in SM . The time-averaged current–current correlator Eq. (4) can be projected onto the subspace formed by local conserved quantities which yields the well-known Mazur–Suzuki equality Mazur (1969); Suzuki (1971) and proves useful for bounding dynamical susceptibilities Zotos et al. (1997). In matrix notation the latter reads $\mathcal{D}^{(i,j)}=\tfrac{\beta}{2}\mathcal{O}_{ik}\,(\mathcal{C}^{-1})_{kl}\,\mathcal{O}_{lj}$ Ilievski and Prosen (2012).

A central result of our work is that on hydrodynamic scale, static charge-charge, charge-current correlations, and generic Drude weights, all assume a universal mode decomposition (writing formally $\mathcal{A}\in\{\mathcal{C},\mathcal{D},\mathcal{O}\}$)

[TABLE]

which has exactly the same form as in the case of a single-component interacting integrable Bose gas derived in a recent paper Doyon and Spohn (2017). Importantly, in the above formula the kernels $A_{a}(u)$ and effective charges $q^{\rm eff}_{a,i}$ are expressible in terms of properties of equilibrium states which can be efficiently computed within Thermodynamic Bethe Ansatz (TBA) method Yang and Yang (1969); Takahashi (1971); Gaudin (1971). It is noteworthy that Eq. (5) is written solely in the mode space, i.e. it acts (diagonally) on particle labels and rapidities, and that no explicit knowledge of a complete set of local charges is ever required in a computation. Indeed, thermodynamic expectation values of local charges are expressible as linear functionals of particles’ rapidity distributions (see e.g. Ilievski et al. (2015, 2016)) which are a natural extension of momentum distribution functions of free theories Ilievski et al. (2017).

Linearized hydrodynamics.

The hydrodynamic approach Bertini et al. (2016); Castro-Alvaredo et al. (2016) is based on the notion of local quasi-stationary states, characterized by the local continuity equation in the mode space $\partial_{t}\rho_{a}(u)+\partial_{x}j_{a}(u)=0$. In the simplest scenario, one can think of a quantum quench in which an inhomogeneous initial state is initialized as two homogeneous equilibrated macroscopic regions brought in contact at $t=0$, see  Spohn and Lebowitz (1977); Bernard and Doyon (2014, 2016). In such a scenario, an emergent nonequilibrium state remains confined to the light cone region determined by particles’ dressed velocities, leading eventually to a quasi-stationary state which depends on the ray coordinate $\zeta=x/t$ and is determined by the condition of vanishing convective derivative $(\partial_{t}+v^{\rm dr}_{a}(u)\partial_{x})\vartheta_{a}(u)=0$.

The setting proves particularly useful for studying nonequilibrium transport properties and, in particular, computation of Drude weights. The latter can be conveniently defined as asymptotic current rates in the limit of vanishing bias $\delta\mu_{j}$ (while keeping other chemical potentials fixed),

[TABLE]

The above prescription has been initially used in Vasseur et al. (2015) and employed in a recent numerical study Karrasch (2017a), while an analogous formula already appeared in an earlier work Ilievski and Prosen (2012). Equation (6) has been recently evaluated in Ilievski and De Nardis (2017); Bulchandani et al. (2017) using the hydrodynamic approach, transforming it first in the light cone coordinates, $\mathcal{D}^{(i,j)}=(\beta/2)\lim_{\delta\mu_{j}\to 0}\int{\rm d}\zeta\,\partial j_{i}(\zeta)/\partial\delta\mu_{j}$, and then computing quasi-stationary currents which are generated by joining together two nearly identical equilibrium states, i.e. imposing a small chemical potential drop at the origin $\mu^{\rm L}_{i}=\mu_{i}+\delta\mu_{i}/2$ and $\mu^{\rm R}_{i}=\mu_{i}-\delta\mu_{i}/2$. Here $\delta\mu_{i}$ has the role of a thermodynamic force, e.g. to study energy transport we identify $\mu_{e}=\beta$.

Just very recently in Doyon and Spohn (2017) the authors applied Eq. (6) to the Lieb–Liniger model and obtained closed-form expressions analogous to Eq. (5). Below we generalize this result to interacting quantum models which involve multiple species of excitations and internal degrees of freedom. It is quite remarkable however that the final outcome remains a bilinear functional operating diagonally in the mode-number space, while the effect of interparticle interactions gets absorbed into a universal renormalization of bare charges, see Eq. (3).

Eq. (6) indicates that Drude weights are expressible as the variation of the equilibrium expectation values of total current Ilievski and De Nardis (2017) with respect to thermodynamic forces $\delta\mu_{j}$, $\mathcal{D}^{(i,j)}=\tfrac{\beta}{2}(\partial J_{i}/\partial\,\delta\mu_{j})_{\delta\mu_{j}=0}=\tfrac{\beta}{2}\sum_{a}\iint{\rm d}\zeta\,{\rm d}u\,q_{a,i}(u;\zeta)(\partial j_{a}(u;\zeta)/\partial\,\delta\mu_{j})_{\delta\mu_{j}=0}$, being the susceptibility of a system to develop ballistic currents. On each ray $\zeta$, the averages of particle current densities are given by Bertini et al. (2016) $j_{a}(\zeta)=(\sigma_{a}\vartheta^{-1}_{a}(\zeta)\delta_{ab}+K_{ab})^{-1}\star e^{\prime}_{b}(\zeta)$, where rapidity dependence has been suppressed for brevity. Given the filling functions inside the light cone $\vartheta_{a}(u;\zeta)=\vartheta^{\rm L}_{a}(u)+\Theta(v^{\rm dr}_{a}(u)-\zeta)\left(\vartheta^{\rm R}_{a}(u)-\vartheta^{\rm L}_{a}(u)\right)$, with the left/right boundary conditions $\vartheta^{\rm L,R}_{a}$, and neglecting corrections of order $\mathcal{O}(\delta\mu^{2})$, one can integrate out the dependence on the light cone coordinates (see SM SM for details). This leads to the form of Eq. (5), with

[TABLE]

On detailed balance.

The symmetry under exchanging indices $i$ an $j$ in representation (5), $\mathcal{D}^{(i,j)}=\mathcal{D}^{(j,i)}$, indicates that the Onsager reciprocal relations Onsager (1931) remain valid for any stationary state, not only in thermal Gibbs equilibrium. This is indeed a general property of the hydrodynamic equation of motion Spohn (1991). Moreover we here show that in a general local equilibrium state of an integrable quantum model, there exist a generalized detailed balance condition on the hydrodynamic scale (i.e for small $\kappa$ and $\omega$), similarly as in the Lieb-Liniger model found recently in Nardis and Panfil (2016); De Nardis et al. (2017). More specifically, given a conserved quantity of the model $\hat{Q}=\sum_{x}\hat{q}_{x}$, the corresponding dynamical structure factor defined as $\mathcal{S}_{\hat{q}}(\kappa,\omega)=\sum_{x}\int{\rm d}te^{{\rm i}(\kappa\,x-\omega\,t)}\langle\hat{q}_{x}(t)\hat{q}_{0}(0)\rangle$ decomposes in terms of individual particle contributions, $\mathcal{S}_{\hat{q}}(\kappa,\omega)=\sum_{a}\mathcal{S}_{\hat{q},a}(\kappa,\omega)$. In the low-momentum limit $\kappa\to 0$, each term is determined by a single matrix element of a particle-hole excitation on a reference equilibrium state SM . Therefore, following the logic presented in Nardis and Panfil (2016), we derive the following generalized reversibility property

[TABLE]

with $\mathcal{F}_{a}(\kappa,\omega)=\kappa\tfrac{\partial}{\partial p_{a}(u)}\log\left(\vartheta_{a}^{-1}(u)-1\right)$, with $u$ fixed by the energy constraint $v_{a}^{\text{dr}}(u)\kappa=\omega$. In the case of thermal (canonical) equilibrium, given by $\vartheta_{a}=(1+\exp{(\beta\,\varepsilon_{a}+\sum_{i}\mu_{a,i}n_{a})})^{-1}$, we have $\mathcal{F}_{a}(\kappa,\omega)=\beta\omega$, which is the usual detailed balance relation.

Hubbard model.

The Hamiltonian of the 1D Hubbard model Gutzwiller (1963); Hubbard (1963) is given as

[TABLE]

where $\hat{T}_{x,x+1}=-\sum_{\sigma={\uparrow},{\downarrow}}\hat{c}^{\dagger}_{x,\sigma}\hat{c}_{x+1,\sigma}+\hat{c}^{\dagger}_{x+1,\sigma}\hat{c}_{x,\sigma}$ is electron hopping and $\hat{V}_{x,x+1}=(\hat{n}_{x,{\uparrow}}-\tfrac{1}{2})(\hat{n}_{x,{\downarrow}}-\tfrac{1}{2})$ is the Coulomb interaction. This model has received a lot of attention in the past decadesShastry (1988); Woynarovich (1989); Frahm and Korepin (1990); Ogata and Shiba (1990); Essler et al. (1991) as well as in the last years Prosen and Žnidarič (2012); Karrasch et al. (2014); Prosen (2014); Seabra et al. (2014); Neumayer et al. (2015); Reiner et al. (2016); Veness and Essler (2016); Veness et al. (2016); Tiegel et al. (2016); Karrasch et al. (2017); Lin et al. (2017); Karrasch (2017a, b); Bolens et al. (2017). We consider the repulsive case ${\mathfrak{u}}\geq 0$, featuring a ${\mathfrak{u}}$-dependent charge gap and gapless spin degrees of freedom.

The Hubbard model is diagonalized by means of nested Bethe Ansatz Essler et al. (1991); Fabian H. L. Essler and Holger Frahm and Frank Göhmann and Andreas Klümper and Vladimir E. Korepin (2005). Eigenstates in a finite system of length $L$ are characterized by quantum numbers which are solutions to Lieb–Wu equations Lieb and Wu (1968) (cf. SM SM ) The model involves two elementary degrees of freedom; the physical particles are momentum-carrying electrons, while spin degrees of freedom represent internal (non-dynamical) excitations described by auxiliary quantum numbers. In a thermodynamic system one finds various types of charge and/or spin-carrying bound states. Specifically, the thermodynamic particle content of the Hubbard model has been derived in Takahashi (1972) (see also Fabian H. L. Essler and Holger Frahm and Frank Göhmann and Andreas Klümper and Vladimir E. Korepin (2005); Quinn and Frolov (2013)) and comprises of (i) spin-up momentum-carrying electronic excitations which carry unit bare (electronic) charge (ii) spin-singlet electronic bound states and (iii) charge-neutral non-dynamical spin-carrying magnonic excitations. A detailed description of the particle content and other information, including explicit expressions for their bare momenta, energies, scattering phases and the dressing transformation, are reported in SM SM .

Numerical results.

We present temperature dependence of charge and spin, see Fig. 1, and thermal Drude weights, see Fig. 2, in grand canonical equilibrium $\hat{\varrho}_{GCE}\simeq\exp{(-\beta\,\hat{H}-\mu\,\hat{N}+B\,\hat{S}^{\rm z})}$, where $\hat{N}=\sum_{x=1}^{L}(\hat{c}^{\dagger}_{x,{\uparrow}}\hat{c}_{x,{\uparrow}}+\hat{c}^{\dagger}_{x,{\downarrow}}\hat{c}_{x,{\downarrow}})$ is total electron charge, and $\hat{S}^{z}=\tfrac{1}{2}\sum_{x=1}^{L}(\hat{c}^{\dagger}_{x,{\uparrow}}\hat{c}_{i,{\uparrow}}-\hat{c}^{\dagger}_{x,{\downarrow}}\hat{c}_{x,{\downarrow}})$, total magnetization. We compared our data with the recent DMRG computation presented in Karrasch (2017b, a). Most notably, at low temperatures appreciably below the charge gap we confirm the ‘Hubbard to Heisenberg crossover’ in the thermal Drude weight observed previously in Karrasch et al. (2016); Karrasch (2017a), see Fig. 2. In SM we also present an exact computation of the asymptotic charge and current profiles inside a light cone and make comparisons with the numerical results of Karrasch (2017b).

Conclusions.

We presented a general theoretical and computational framework to access the singular components (Drude weights) of generalized transport coefficients in quantum integrable models. We exemplified our approach by computing exact numerical values of (diagonal) charge, spin and thermal Drude weights in the one-dimensional fermionic Hubbard model in grand canonical equilibrium at finite temperatures and chemical potentials. Using the two-partition protocol, we additionally computed the quasi-stationary energy and charge density profile and the corresponding current  SM .

Our results finally permit to establish the equivalence of various approaches for computing the spin Drude weight employed in the previous literature: (i) using projections onto local conserved subspaces by virtue of Mazur–Suzuki equality Zotos et al. (1997); Prosen (2011); Prosen and Ilievski (2013), (ii) taking the linear-response limit of the asymptotic current rates Ilievski and De Nardis (2017); Bulchandani et al. (2017) and (iii) computing the energy-level curvatures Shastry and Sutherland (1990); Fujimoto and Kawakami (1998); Zotos (1999); Benz et al. (2005) under the twisted boundary conditions in accordance with Kohn formula Kohn (1964). The latter has been evaluated within the TBA framework in Fujimoto and Kawakami (1998); Zotos (1999), yielding a closed formula expressed in terms of filling functions, magnonic dispersion relations and $\mathcal{O}(1/L)$ corrections to the Bethe spectrum induced by the twist. Remarkably however, it is easy to see that the twist-dependence of the energy levels can be directly linked to the effective spin as given by Eq. (2). This in turn reconciles the results of Zotos (1999) with Eq. (7), representing the equilibrium analogue of definition (6) used previously in refs. Ilievski and De Nardis (2017); Bulchandani et al. (2017) (further details are given in SM SM , which also includes refs. Narozhny et al. (1998); Peres et al. (1999); Alvarez and Gros (2002); Heidrich-Meisner et al. (2003); Fujimoto and Kawakami (2003); Herbrych et al. (2011); Sirker et al. (2011); Žnidarič (2011); Gromov and Kazakov (2011); Karrasch et al. (2012); Cavaglià et al. (2015); De Luca et al. (2016)).

Finally, our results show that a generalized version of the detailed balance Nardis and Panfil (2016); De Nardis et al. (2017); Foini et al. (2017) is valid on hydrodynamic scales in any stationary state.

As a future task, it would be interesting to find an extension of the presented approach which would allow resolving the diffusive time-scale from the microscopic picture, see e.g. Medenjak et al. (2017a, b).

Authors contributions.

Both authors contributed equally to the theory. J. De Nardis performed the numerical computations.

Acknowledgements.

We are grateful to C. Karrasch for providing the tDMRG data for the Drude weights in Hubbard model and thank M. Van Caspel, M. Fagotti, E. Quinn and H. Spohn for valuable discussions and/or reading the manuscript. E.I. is supported by VENI grant number 680-47-454 by the Netherlands Organisation for Scientific Research (NWO). J.D.N. acknowledge support by LabEx ENS-ICFP:ANR-10-LABX-0010/ANR-10-IDEX-0001-02 PSL*.

Appendix A Thermodynamic Bethe Ansatz for Hubbard model

The Hamiltonian of 1D Hubbard model is of the form

[TABLE]

The model possesses two globally conserved charges associated with $U(1)$ symmetries, the total electron charge $\hat{N}$ and the total spin $\hat{S}^{z}$,

[TABLE]

which are sometimes included in the definition Hamiltonian, $\hat{H}=\hat{H}_{0}+\mu_{\rm c}\,\hat{N}+\mu_{\rm s}\,\hat{S}^{z}$.

Bethe equations for a finite system of length $L$ with periodic boundary conditions have been derived by Lieb and Wu Lieb and Wu (1968) and take the nested form,

[TABLE]

with $2M\leq N\leq L$. Bethe roots (rapidities) $u_{k}$ are related to electron (quasi)momenta, while $w_{k}$ are associated with their spin. The number of Bethe roots in Eqs. (13) in terms of the total charge and spin is $N$ and $(N-2M)/2$, respectively. Bethe roots are associated the bare charge $n_{u}=1$, $n_{w}=0$, and the bare spin $m_{u}=\tfrac{1}{2}$, $m_{w}=-1$. We note that parametrization of Lieb–Wu equations (13) in terms of $u$-roots is different from the conventional one given in terms of electron (quasi)momenta $p_{j}$ as in Takahashi (1972). While the two are simply related by $u_{j}=\sin{(k_{j})}$, rapidity parametrization proves more convenient since it renders all scattering amplitudes manifestly rational functions depending only on the difference of particles’ rapidities. A downside is that the momentum-dependent phase $e^{{\rm i}k(u)}$ as a function of momentum-carrying roots $u_{i}$ then becomes a double-valued function, meaning that each root $u_{i}$ gives two distinct values of momenta. It is thus convenient to introduce a new type of roots, referred to as the $y$-roots, by virtue of Zhukovsky transform

[TABLE]

The corresponding functional equation $\frac{1}{2}\left(x(u)+1/x(u)\right)=u$ has two solutions (branches), and presently we adopt

[TABLE]

with a square-root branch cut on the interval $\mathcal{I}=(-1,1)$. For any $u\in\mathbb{C}$ the two branches correspond to the values $y_{\pm}(u)$ given by

[TABLE]

When rapidity $u$ is taken from the branch cut, $u\in\mathcal{I}$, we adopt the following prescription

[TABLE]

i.e. we take the two values just above and just below the cut $\mathcal{I}$. Since we have

[TABLE]

the two branches of momenta are given by

[TABLE]

Thermodynamic solutions to Eqs. (13) – taking the limit $L\to\infty$ while keeping ratios $N/L$ and $M/L$ finite – can be inferred from the stability condition of the asymptotic solutions, and comprise of self-conjugate string-like patterns of regularly displaced complex-valued rapidities with equal real parts centred on the real axis. These are identified with the thermodynamic particle content of the model which in the Hubbard model and comprise of:

• •

$y$-particles, which are spin-up $m_{y}=\tfrac{1}{2}$ momentum-carrying electronic excitations which carry unit electron charge $n_{y}=1$. The $y$-particle excitations are split into two branches denoted by $y_{\pm}$ with the corresponding rapidities $u_{\pm j}$. The $y$-particles do not form bound states on their own. Their bare momenta are denoted by $k_{\pm}(u)\in\mathbb{R}$, and satisfy $k^{\prime}_{-}>0$ and $k^{\prime}_{+}<0$, with $k^{\prime}_{+}+k^{\prime}_{-}=0.$. The lower and upper momentum branches of the $y$-particle are

[TABLE]

with the corresponding derivatives

[TABLE]

The bare energies are $e_{\pm}(u)$ and read

[TABLE]

• •

$M|uw$-strings, which are bound states of $2M$ $u$-roots and $M$ $w$-roots, carrying charge $n_{M|uw}=2M$ and no spin $m_{M|uw}=0$. An $M|uw$-string is parametrized by $u\in\mathbb{R}$ and comprises of rapidities

[TABLE]

To find the corresponding $y$-roots we assign $y_{+j}=x(u_{j})$ and $y_{-j}=1/x(u_{-j})$. The corresponding momenta and energies are obtained by summing over all constituent $u$-roots (recall that $w$-roots carry no momentum)

[TABLE]

The derivative of their momenta satisfy $k^{\prime}_{M|uw}<0$, reading explicitly

[TABLE]

• •

$M|w$-strings, which are chargless ($n_{M|w}=0$) compounds made of $M$ $w$-roots, with spin $m_{M|w}=-M$. They are parametrized by $u\in\mathbb{R}$, and are of the form

[TABLE]

The total charge $N$ and total spin $M$ in a Bethe eigenstate in terms of numbers of string excitations $N_{a}$ ($a=\{\pm,M|uw,M|w\}$) is

[TABLE]

Thermodynamic limit.

In the thermodynamic limit, the solutions to Eqs. (13) become densely distributed on the rapidity axis and can be expressed in terms of particle densities $\rho_{a}$ which are defined as smooth densities of Bethe strings $u^{(a)}_{j}$,

[TABLE]

Given a set of string solutions $\{u^{(a)}_{j}\}$, the unoccupied solutions to Bethe equations $e^{{\rm i}p_{a}(u)}\prod_{b}\prod_{j=1}^{N_{b}}S_{ab}(u-u^{(b)}_{j})=1$ are understood as the holes. The hole densities in the thermodynamic limit are denoted by $\bar{\rho}_{a}(u)$, while the total densities of a state are denoted by $\rho^{t}_{a}(u)=\rho_{a}(u)+\bar{\rho}_{a}(u)$.

In the Hubbard model, the densities of $M|uw$- and $M|w$-strings, denoted by $\rho_{M|uw}(u)$ and $\rho_{M|w}(u)$, respectively, are supported on the whole real axis, $u\in\mathbb{R}$. On the other hand, the rapidity distributions of the special $y$-particles are split into two separate densities $\rho_{\pm}(u)$ which are compactly supported on the branch cut $\mathcal{I}$. To this end, we define two types of integral transformations. First, we introduce the standard convolution as

[TABLE]

taken with the convention that one drop $u$ and/or $w$ when $g$ and/or $h$ depend only on a single variable, and adopting the implicit summation convention for convolution expressions of the form $g_{ab}\star h_{bc}$ and $g_{ab}\star g_{b}$, namely summing and integrating over the domain of a $b$-string. Since the $y$-particles’ rapidity variable have a bounded integration domain, i.e. $z\in\mathcal{I}$, we introduced a restricted convolution operation $\,\hat{\star}\,$. The densities of $y$-particles satisfy the sum rule,

[TABLE]

Denoting $n_{M|uw}\equiv 1\star\rho_{M|uw}$, $n_{M|w}\equiv 1\star\rho_{M|w}$, and $n_{\pm}\equiv 1\,\hat{\star}\,\rho_{\pm}$, the electron charge and spin densities are expressed as

[TABLE]

Energy density of a macroscopic state is obtained by adding contributions of all energy-carrying particles,

[TABLE]

Takahashi’s equations for the densities in rapidity parametrization take the form

[TABLE]

The explicit form of integral kernels $K_{M}$ and $K_{MN}$ are given in section A.2. Equations (36) can be, using fusion identities among scattering kernels (cf. section A.2), decoupled in a quasi-local form

[TABLE]

where $\delta_{MN}$ is the Kronecker delta, $\delta$ the Dirac delta, and the $I$ is the adjacency (incidence) matrix for the model,

[TABLE]

Local statistical ensembles.

Thermodynamic Bethe Ansatz method is based on expressing the free energy density of a local statistical ensemble (a generalized Gibbs ensemble) as a set of coupled non-linear integral equations for the thermodynamic variables (e.g. Fermi filling functions of the thermodynamic excitations). Generalized Gibbs ensembles are conventionally expressed in the form

[TABLE]

for a suitable (discrete) basis of local conserved quantities $\hat{Q}_{i}$ and the corresponding chemical potentials $\mu_{i}$. By accounting for the fact that particles’ mode distributions $\rho_{a}(u)$ essentially contain the complete information about local correlations functions, it is convenient to consider as a starting point an analytic parametrization Ilievski et al. (2017)

[TABLE]

where $\hat{\varrho}_{a}(u)$ correspond formally to a continuous family of local conserved operators whose eigenvalues coincide with the particles’ rapidity distributions, and $\mu_{a}())u$ are the chemical potentials pertaining to individual modes. The partition sum $\mathcal{Z}_{\rm GGE}={\rm Tr}\,\hat{\varrho}_{\rm GGE}$ in the $L\to\infty$ limit is then evaluated with a saddle-point integration using Yang–Yang approach Yang and Yang (1969), where the entropy density per particle is the logarithm of the number states occupying an infinitesimal rapidity interval $[u,u+{\rm d}u]$ which (in models obeying the Fermi statistics) takes a universal form

[TABLE]

A solution to the variational problem $\delta f[\rho_{a}]=0$, with $f=\sum_{a}\int{\rm d}u\,(\mu_{a}(u)\rho_{a}(u)-\mathfrak{s}_{a}(u))$, yields canonical TBA equations

[TABLE]

where the TBA $Y$-functions are as usual defined as ratios of hole and particle densities for each thermodynamic excitation in the spectrum,

[TABLE]

The set of $Y$-functions is equivalent to the set of Fermi filling functions $\vartheta_{a}$, defined as $\vartheta_{a}(u)=\rho_{a}/\rho^{t}_{a}$. For later purposes we moreover introduce the filling functions of the holes, $\bar{\vartheta}_{a}(u)=1-\vartheta_{a}(u)$.

For instance, in canonical Gibbs equilibrium, $\hat{\varrho}_{\rm Gibbs}\sim\exp{(-\beta\hat{H}+\mu_{\rm c}\hat{N}+\mu_{\rm s}\hat{S}^{z})}$, the canonical source terms in terms of particles’ bare energies $e_{a}(u)$ and chemical potentials for the electronic charge and spin read

[TABLE]

Using the fusion identities (cf. A.2), the above set of equations can be brought to an equivalent quasi-local form, reading explicitly

[TABLE]

supplemented with the asymptotic conditions

[TABLE]

By making the identifications $Y_{M|uw}\equiv Y_{M+1,1}$ and $Y_{M|w}\equiv Y_{1,M+1}$, the Y-functions may be inscribed in the so-called $Y$-lattice (see e.g. Beisert et al. (2011); Cavaglià et al. (2015)), as shown on Fig. 3.

A.1 Dressing of excitations and effective charges

Excited states with respect to a reference macrosopic state (representing a many-body vacuum) are characterized in terms of the particle-hole type of excitations and a background of non-excited modes (quantum numbers) which experience a $\mathcal{O}(1/L)$ shift as a back-reaction to creating excitations. The difference between the rapidities of excited and reference states induced by $N^{b}_{\rm ex}$ particle-type of excitations of type $b$ can be expressed as

[TABLE]

while the hole-type excitations experience the same the shift of the opposite sign. The shift functions $F_{ab}(u,w)$, describing the back-flow of non-excited rapidities, satisfy a closed set of integral equations

[TABLE]

In the thermodynamic limit, Eq. (46) can be expressed as an integral equation which governs the dressing of bare quantities $q_{a}(u)$ (suppressing rapidity parameters)

[TABLE]

Differentiating this expression with respect to rapidity variable we find

[TABLE]

where the inverse of the dressing convolution kernel $\Omega$ explicitly reads

[TABLE]

Two special (but central) examples of the above transformation are the dressed energies $\varepsilon_{a}$ and dressed momenta $p_{a}$, providing dispersion relations of the particle-hole excitations with respect to a reference macrostate. The dressed velocities yield the group velocity of propagation and are given by

[TABLE]

Notice that particle velocities, which are given as a ratio of the derivatives of two dressed quantities, do not satisfy the universal dressing equation. Moreover, it is worthwhile stressing that $\varepsilon^{\prime}\neq(e^{\prime}_{a})^{\rm dr}$ since $[\partial_{u},\Omega_{ab}\star]\neq 0$. Hence, to avoid confusion, when applying the dressing transformation to the bare charge densities we shall speak of the effective charges, that is

[TABLE]

The effective charges for the charge density $q_{i}$ can be alternatively obtained from the $Y$-functions of an equilibrium state parametrized in the form of Eq. (38) as

[TABLE]

In the Hubbard model, the derivatives of the dress charges $q^{{\rm dr}\prime}_{a}$ are uniformly expressed as the solution to the following system of integral equations,

[TABLE]

For example, choosing $p^{\prime}_{a}=k^{{\rm dr}\prime}_{a}=2\pi\sigma_{a}\rho^{t}_{a}$ reduces Eq. (54) to Eq. (36), while Eqs. (44) are retrieved by plugging in $\varepsilon^{\prime}_{a}=e^{{\rm dr}\prime}_{a}$. Let us note that Eqs. (54) comply with the morphology of the $Y$-system lattice, see Fig. 3.

The electron charge and spin have an exceptional role since they do not depend on rapidities variables. Their effective values for the electron charge follow from the solution to

[TABLE]

along with the large-$M$ asymptotic condition $\lim_{M\to\infty}n^{\rm eff}_{M|uw}=2M$. Similarly, for the effective spin we have

[TABLE]

with the asymptotics $\lim_{M\to\infty}m^{\rm eff}_{M|w}=M$. Therefore, only $M|uw$-strings and $y$-particles yield non-vanishing effective charge, while $M|w$-stings and $y$-particle yield non-vanishing effective spin. On the other hand, all types of particles (inducing the auxiliary ones) have non-zero effective energies $e^{\rm eff}_{a}$ in general. The effective electron charges and spin can also be obtained from

[TABLE]

A.2 Scattering data and fusion identities

The elementary scattering amplitudes are

[TABLE]

The scattering amplitudes for the regular $M|uw$-strings and $M|w$-strings are given by

[TABLE]

whereas the scattering amplitudes between $y$-roots and $M|uw$-strings or $M|w$-strings are $S_{M}(u)$. The TBA integral kernels are defined as derivatives of the logarithmic scattering amplitudes,

[TABLE]

The kernels for $y$-particles are similarly given by

[TABLE]

for all types of particles $a$. Notice also that $K_{+a}=K_{-a}$.

Canonical TBA equations (41) can be cast in an equivalent quasi-local description by employing the following fusion identities among the scattering kernels,

[TABLE]

with $K_{0}\equiv 0$. In addition, the latter satisfy

[TABLE]

In the canonical (Gibbs) equilibrium, the action of $(K+1)^{-1}\star$ on the bare energies $e_{N|uw}$ yields

[TABLE]

Additionally, the terms which involve $k^{\prime}_{M|uw}$ can be simplified using

[TABLE]

Appendix B Static covariance matrix

To obtain the full static charge-current covariance matrix, i.e. overlap coefficients $\mathcal{O}_{ij}$, we consider

[TABLE]

To facilitate calculations, we introduce a compact vector notation for quantities depending solely on the mode labels $a$ and $u$, namely $\mathbf{q}=(q_{1}(u),q_{2}(u),\ldots)$, and similarly for other quantities. Thus, expressing particle currents as

[TABLE]

readily yields

[TABLE]

Employing $(\hat{\sigma}\hat{\vartheta}+\hat{K})^{-2}=\hat{\vartheta}^{2}\,\hat{\Omega}^{2}$, and noticing that $\hat{\vartheta}^{-1}$ is a diagonal operator in the mode space, we find

[TABLE]

The mode kernel $O_{a}(u)$ for the charge–current correlator $\mathcal{O}_{ij}=\sum_{a}\int{\rm d}u\,q^{\rm eff}_{a,i}(u)O_{a}(u)q^{\rm eff}_{a,j}(u)$ therefore takes the form

[TABLE]

Alternatively, the static covariances can also be derived from the second derivatives of a functional

[TABLE]

By setting $g=\{k,e\}$ and calculating the gradients one obtains the well-known mode resolutions

[TABLE]

Note that $f_{k}$ is the diagonal representation of the standard (generalized) free energy density $f$, which can be readily deduced from combining Bethe–Yang equations for the densities, Yang–Yang entropy and TBA equations for $\log Y_{a}$. The first derivatives of $\log Y_{a}$ with respect to chemical potentials give the effective charges,

[TABLE]

which readily implies

[TABLE]

From the identity

[TABLE]

the second derivatives of $f_{g}$ take the form

[TABLE]

In the second line we have used $\hat{1}-\hat{\sigma}\hat{\vartheta}(\hat{1}+\hat{K}\hat{\sigma}\hat{\vartheta})^{-1}\hat{K}=\hat{\Omega}$. After expressing the derivatives of the filling functions $\vartheta_{a}$ as

[TABLE]

we obtain

[TABLE]

in turn implying

[TABLE]

with the corresponding mode kernels

[TABLE]

Comment 1.

The Bare quantities which take constant values (i.e. do not depend on rapidities) require careful considerations. Considering spin of excitations as an example, the effective spin is determined via

[TABLE]

Since the bare spins $m_{a}$ of spin-carrying excitations are all non-zero, one would expect that the same holds automatically also for the corresponding effective values. While this is in general the case, in the limit of vanishing chemical potential $h\to 0$ the effective spin may tend to zero. The reason why this does not conflict with Eq. (84) is an infinite summation over the particle content. Thus, in the presence of infinitely many types of excitations, the dressing transformation is of infinite dimension and demands an appropriate regularization. To this end let us consider the Gibbs equilibrium state in $\beta\to 0$ limit where to the leading order in $h$ we have $m^{\rm eff(0)}_{a}(h)=\tfrac{2}{3}(a+1)^{2}h+\mathcal{O}(h^{3})$, i.e. yields a results which is well-behaved in the $h\to 0$ limit but diverges as $\sim a^{2}$ in the limit of large strings (large bare spin). The corresponding filling functions on the other hand converge asymptotically for large $a$ as $\sim 1/a^{2}$, and read $\vartheta^{(0)}_{a}=1/(a+1)^{2}$.

A handy way to regularize the canonical form of Eq. (84) is to use the fusion identities for the scattering kernels to transform it in the quasi-local form. Using the fact that $(K_{ab}+\delta_{ab}\delta)\star m_{b}=0$, Eq. (84) is readily transformed into

[TABLE]

The erased source term is substituted with an appropriate large-$a$ asymptotic condition. For finite $h$ we should impose $\lim_{a\to\infty}m^{\rm eff}_{a}=a$, which follows from the large-$a$ asymptotics of $\log Y_{a}$ in the canonical TBA equations. In the $\beta\to 0$ limit and finite $h$, the Y-functions $\log Y^{(\infty)}_{a}(h)$ take the form $1+\log Y^{(\infty)}_{a}(h)=\sinh^{2}{(h(a+1))}/\sinh^{2}{(h)}$. The $h\to 0$ limit is achieved by extrapolation.

Comment 2.

One often deal with a situation when the values of effective charges $\mathbf{q}^{\rm eff}_{j}$ exactly vanish, e.g. $m^{\rm eff}_{a}=0$ at zero chemical potential $h=0$. In such a case Eq. (78) has to be regularized by considering small $h$ and only taking the limit $h\to 0$ at the end of computation, after first performing the infinite mode summation and a non-compact integration in Eq. (81). In the opposite case, Eq. (78) would imply vanishing derivatives $\partial_{h}s^{\rm eff}_{a}$ and thus an incorrect result $\chi_{s}\equiv\mathcal{C}_{ss}=0$.

Appendix C Drude weights from linearized hydrodynamics

Drude weights can be defined as the variation of the equilibrium expectation values of the total currentIlievski and De Nardis (2017) with respect to thermodynamic forces $\delta\mu_{j}$,

[TABLE]

We imagine a bipartite initial state with a chemical potential drop $\delta\mu_{j}$ at the origin (while keeping other chemical potential fixed). The filling functions inside the light cone $\vartheta_{a}(u;\zeta)=\vartheta^{\rm L}_{a}(u)+\Theta(v^{\rm dr}_{a}(u)-\zeta)\left(\vartheta^{\rm R}_{a}(u)-\vartheta^{\rm L}_{a}(u)\right)$, with the left/right boundary conditions $\vartheta^{\rm L,R}_{a}$ which differ by amount $\mathcal{O}(\delta\mu_{j})$. The corrections due to the difference of particle velocities only enter in the sub-leading order and can be disregarded. On every ray $\zeta$, the particle current densities can be expressed as

[TABLE]

where $\hat{\Xi}=\hat{\Xi}(\zeta)$ is a diagonal operator which involves the jump discontinuity,

[TABLE]

Writing $\hat{A}=\hat{\sigma}\,\hat{\vartheta}\,\hat{\Omega}$, expanding the inverse $(\hat{A}^{-1}+\hat{\Xi})^{-1}=\hat{A}+\hat{\Xi}\hat{A}^{2}$, using the identities

[TABLE]

and finally integrating over the light cone region,

[TABLE]

yields

[TABLE]

In the last equality we have used $\varepsilon^{\prime}_{a}=2\pi\sigma_{a}\rho^{t}_{a}v^{\rm dr}_{a}$.

Appendix D On the detailed balance

It was shown in Nardis and Panfil (2016) that in the Lieb-Liniger model the dynamical density structure, in the low-momentum limit $\kappa\to 0$, is determined by a single particle-hole excitation,

[TABLE]

where $\hat{\rho}$ is the conserved density operator and $|\langle\vartheta|\hat{\rho}|\vartheta,u\to\tilde{u}\rangle|^{2}$ represents a matrix element of the single particle-hole excitation $u\to\tilde{u}$ with energy $\omega$ and small momentum $\kappa$. Recall that the matrix element is proportional to the available density of states for the particle-hole excitation, namely

[TABLE]

Therefore, expanding it to the first order in $\kappa$ around $\kappa=0$ we find

[TABLE]

implying a detailed balance condition in the order $\mathcal{O}(\kappa)$ of the form

[TABLE]

where the function $\mathcal{F}(\kappa,\omega)$ is given by

[TABLE]

The recent results of Doyon and Spohn (2017), lifting standard hydrodynamics results Spohn (1991) to the generalized hydrodynamic theory Bertini et al. (2016); Castro-Alvaredo et al. (2016), imply that the dynamical structure factor for any conserved local charge $\hat{q}$ is characterized by a single particle-hole contribution, with energy equal to $\kappa\,v^{\rm dr}(u)$. Models with multiple particle require an additional summation over all particle types, and formula (92) generalizes to

[TABLE]

where, using Eq. (5), the zero momentum limit of the matrix element is given by

[TABLE]

By repeating the logic of Nardis and Panfil (2016), we obtain the detailed balance expression for each particle type

[TABLE]

Appendix E Spin Drude weight in the anisotropic Heisenberg chain

We briefly revisit the exceptional case of ballistic spin transport in the anisotropic Heisenberg model at half filling, a phenomenon which has attracted considerable attention in the past, see e.g. Narozhny et al. (1998); Peres et al. (1999); Zotos (1999); Alvarez and Gros (2002); Heidrich-Meisner et al. (2003); Fujimoto and Kawakami (2003); Benz et al. (2005); Herbrych et al. (2011); Sirker et al. (2011); Žnidarič (2011); Karrasch et al. (2012); Karrasch (2017a). As explained and discussed in Ilievski and De Nardis (2017), the peculiar behaviour of the finite-temperature spin Drude weight at $\mu_{\rm s}=0$ – which is vanishing outside of the critical interval $|\Delta|<1$ where it exhibits a nowhere-continuous dependence on interaction anisotropy $\Delta$ – is directly related to the formation of an exceptional pair of excitations which are charged under a hidden non-unitary conservation law De Luca et al. (2016).

In Ilievski and De Nardis (2017), the anomaly has been explained on the basis of symmetry properties of thermodynamic states under the spin-reversal transformation, which rigorously confirmed the exact analytical high-temperature bound derived earlier in Prosen and Ilievski (2013). Given the that dressing of particle excitations is a property of the reference local equilibrium state, the aim of this section is to shortly revisit this interesting case from the point of view of the linearized hydrodynamics. As explained below, this allows to reconcile our results with the previously obtained analytical resulting derive from alternative approaches.

The finite-temperature spin Drude weights has been initially computed with aid of the Kohn formula, cf. Shastry and Sutherland (1990); Fujimoto and Kawakami (1998); Zotos (1999); Benz et al. (2005), expressing $\mathcal{D}^{(s)}$ as the thermal average of the energy-level curvatures with respect to applying twisted boundary conditions Kohn (1964) (or equivalently, piercing a ring with a magnetic flux $\phi$). By resolving the $\mathcal{O}(L^{-2})$ corrections of the spectrum using TBA approach, ref.Zotos (1999) finds an explicit expression for the spin Drude weight in the anisotropic Heisenberg model, which in our notation reads

[TABLE]

where functions $\gamma_{a}$ describe $\mathcal{O}(1/L)$ shifts of Bethe roots due to the twisted boundaries. The $\phi$-derivatives of $\gamma_{a}$ satisfy the following integral equations

[TABLE]

Identifying $m^{\rm eff}_{a}=2\pi\rho^{t}_{a}\partial_{\phi}\gamma_{a}$ now allows us to interpret Eq. (101) as the dressing of spin, which readily implies that $\mathcal{D}^{(s)}_{\rm Kohn}$ indeed agrees with

[TABLE]

In the Heisenberg spin chain, the thermodynamic particle are magnons and bound states thereof, carrying a finite amount of bare spin. Their effective spins with respect to a generic local equilibrium state are thus given by some non-trivial finite quantities. Nevertheless, when approaching half filling $h\to 0$ (at finite temperatures) their effective spin exactly vanishes in the gapped and isotropic regimes ($|\Delta|\geq 1$). In the gapless regime on the other hand, the number of distinct stable particles reduces to a finite set and consequently the effective spin cannot entirely vanish for all excitations in the $h\to 0$ limit. Remarkably, it turn our that the effective spin in a half-filled state exactly vanishes for all magnonic particles with an exception of a distinguished pair of particles, in Ilievski and De Nardis (2017) labelled by $a=\bullet,\circ$, which at $\Delta=\cos{(\pi m/\ell)}$ carry a finite (temperature-independent) effective spin $m^{\rm eff}_{\bullet,\circ}=\pm\ell/2$. It thus follows from the mode resolution (5) that only this special pair of magnonic bound states contribute to spin Drude weight, and since the latter have been shown to be the only excitations which transport non-unitary local conservation laws found in Prosen (2011); Prosen and Ilievski (2013), this automatically implies that the exact Mazur projection calculated in ref. Prosen and Ilievski (2013) is complete. The main conclusion of this section is thus that all three different definitions of spin Drude weights are equivalent.

Appendix F Hydrodynamic description of the Hubbard model

We consider a non-equilibrium protocol such that the initial state is a tensor product of two different macroscopic states joined at the origin ($x=0$). This is usually referred to as the bi-partite protocol. In the long-time limit $t\to\infty$, the system is described locally by a quasi-stationary state which depends on the ray direction $\zeta=x/t$. The properties of such states can be computed from the generalized hydrodynamic theory Bertini et al. (2016); Castro-Alvaredo et al. (2016) which is formally a kinetic theory for the thermodynamic degrees of freedom of a system. In the Hubbard model, thermodynamic excitations comprise of spin-up electronic excitations (called the $y$-particles), the $M|uw$-strings (spinless electonic bound state) and the $M|w$-strings (spin-carrying chargeless bound states). Their root densities in position $x=\zeta t$ satisfy the usual hydrodynamic continuity equation

[TABLE]

Solving the above set of equations automatically yields the correct distributions of the auxiliary particles $\rho_{M|w}$. The solution to Eqs. (103) is given as usual in terms of the filling functions $\vartheta_{a}$

[TABLE]

which completely determine the expectation values for local operators within the light cone emanating from the junction of the two initial equilibrium states. We wish to point out that the $M|w$-strings, being the auxiliary degrees of freedom, do not supply and dynamical information, but merely adapt their values to those of the physical ones, that is the $y$ particles and the $M|uw$ strings. This implies that one can either obtain functions $\vartheta_{M|w}(u;\zeta)$ with Eq. (105) or, alternatively, compute them via the TBA equations (suppressing dependence on $u$)

[TABLE]

with $Y_{y}(\zeta)=\vartheta_{y}^{-1}(\zeta)-1$ given by Eq. (104). These two schemes yield equivalent descriptions of the stationary state. To demonstrate the above procedure, we present the profiles of the energy density and its current,

[TABLE]

[TABLE]

and moreover of the charge density operator

[TABLE]

together with the charge current

[TABLE]

We consider a thermal state with chemical potentials set to $(\beta,\mu,B)=(1,0,0)$ on one side, and the bare vacuum state on the other (i.e. a thermal state with $\mu\to\infty$). Numerical results of our theoretical framework, obtained as the solution to Eqs. (104) perfectly reproduce the DMRG data at late times as shown on Fig. 4.