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

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 ¹¹1In 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^{i{p}_{\alpha }\left(u_k^{\left(\alpha \right)}\right)}{\prod }_{\beta }{\prod }_{j=1}^{N_{\beta }}{S}_{\alpha \beta }(u_k^{\left(\alpha \right)},u_j^{\left(\beta \right)})=1$, imposing single-valuedness of many-body eigenstates. Here the sets of quantum numbers $\left\{u_k^{\left(\alpha \right)}\right\}$ 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\left(u\right)$ and energies $e_a\left(u\right)$, and interparticle scattering phase shifts ${\varphi }_{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\left(x\right)$, with $x$ labelling lattice sites. The associated currents $J_i={\sum }_x{j}_i\left(x\right)$ 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 {\epsilon }_a$ and momenta $k_a\mapsto p_a$ of particle excitations, which can be presented in a compact form

with convolution $({\Omega }_{ab}\star f_b)\left(u\right)={\sum }_b\int dw{\Omega }_{ab}(u,w)f_b\left(w\right)$. In interacting quantum integrable models solvable by (nested) Bethe Ansatz, the matrix convolution kernel $\Omega$ takes a universal form

with kernels $K_{ab}(u,w)$ defined as derivatives of the scattering phase shifts ${\varphi }_{ab}(u,w)={\varphi }_{ab}(u-w)$, $K_{ab}\left(u\right)=\frac{1}{2\pi i}{\partial }_u{\varphi }_{ab}\left(u\right)$, and ${\sigma }_a=sign\left(k_a^{\prime }\right(u\left)\right)$. The (Fermi) filling functions ${\vartheta }_a\left(u\right)$ specify the fraction of occupied modes with rapidities inside a small interval around $u$.

Dispersion relations of excitations ${\epsilon }_a\left(u\right)$ depend on a many-body vacuum which is uniquely specified by the rapidity distributions ${\rho }_a\left(u\right)$. In terms of (thermodynamic) particle excitations, the equilibrium averages of charge and current densities decompose as $q_i={\sum }_a\int du{q}_{i,a}\left(u\right){\rho }_a\left(u\right)$, $j_i={\sum }_a\int du{q}_{i,a}\left(u\right)j_a\left(u\right)$, where $j_a\left(u\right)={\rho }_a\left(u\right)v_a^{dr}\left(u\right)$ 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_a^{dr}\left(u\right)={\epsilon }_a^{\prime }\left(u\right)/p_a^{\prime }\left(u\right)$.

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

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 ${\epsilon }_a^{\prime }={\Omega }_{ab}\star e_b^{\prime }=(e_a^{\prime }{)}^{eff}$, the effective charges are not the proper dressed charges associated with an excitation, and specifically ${\epsilon }_a\ne e_a^{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.

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

which represent magnitudes of the singular parts of the zero-frequency generalized conductivities Kubo (1957); Mahan (2000), $Re{\sigma }_{ij}\left(\omega \right)=2\pi D^{(i,j)}\delta \left(\omega \right)+{\sigma }_{ij}^{reg}\left(\omega \right)$. We use $⟨\cdot {⟩}_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 $D^{(i,j)}$ can be given in terms of the static covariance matrix $C$, $C_{ij}=⟨Q_i{q}_j{⟩}_c$, with diagonal components ${\chi }_i=C_{ii}$ representing (generalized) static susceptibilities, and charge-current correlators (overlaps) $O$, $O_{ij}=⟨Q_i{j}_j{⟩}_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 $D^{(i,j)}=\frac{\beta }{2}{O}_{ik}(C^{-1}{)}_{kl}{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 $A\in \{C,D,O\}$)

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\left(u\right)$ and effective charges $q_{a,i}^{eff}$ 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).

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\left(u\right)+{\partial }_x{j}_a\left(u\right)=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_a^{dr}(u\left){\partial }_x\right){\vartheta }_a\left(u\right)=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),

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, $D^{(i,j)}=\left(\beta /2\right){lim}_{\delta {\mu }_j\to 0}\int d\zeta \partial j_i\left(\zeta \right)/\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 }_i^L={\mu }_i+\delta {\mu }_i/2$ and ${\mu }_i^R={\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$, $D^{(i,j)}=\frac{\beta }{2}(\partial J_i/\partial \delta {\mu }_j{)}_{\delta {\mu }_j=0}=\frac{\beta }{2}{\sum }_a\iint d\zeta du{q}_{a,i}(u;\zeta \left)\right(\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\left(\zeta \right)=\left({\sigma }_a{\vartheta }_a^{-1}\right(\zeta ){\delta }_{ab}+K_{ab}{)}^{-1}\star e_b^{\prime }(\zeta )$, where rapidity dependence has been suppressed for brevity. Given the filling functions inside the light cone ${\vartheta }_a(u;\zeta )={\vartheta }_a^L\left(u\right)+\Theta \left(v_a^{dr}\right(u)-\zeta )\left({\vartheta }_a^R\left(u\right)-{\vartheta }_a^L\left(u\right)\right)$, with the left/right boundary conditions ${\vartheta }_a^{L,R}$, and neglecting corrections of order $O\left(\delta {\mu }^2\right)$, 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

The symmetry under exchanging indices $i$ an $j$ in representation (5), $D^{(i,j)}=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 $S_{\hat{q}}(\kappa ,\omega )={\sum }_x\int dt{e}^{i(\kappa x-\omega t)}⟨{\hat{q}}_x\left(t\right){\hat{q}}_0\left(0\right)⟩$ decomposes in terms of individual particle contributions, $S_{\hat{q}}(\kappa ,\omega )={\sum }_a{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

with $F_a(\kappa ,\omega )=\kappa \frac{\partial }{\partial p_a\left(u\right)}\log \left({\vartheta }_a^{-1}\left(u\right)-1\right)$, with $u$ fixed by the energy constraint $v_a^{\text{dr}}\left(u\right)\kappa =\omega$. In the case of thermal (canonical) equilibrium, given by ${\vartheta }_a=(1+\exp (\beta {\epsilon }_a+{\sum }_i{\mu }_{a,i}{n}_a){)}^{-1}$, we have $F_a(\kappa ,\omega )=\beta \omega$, which is the usual detailed balance relation.

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

where ${\hat{T}}_{x,x+1}=-{\sum }_{\sigma =\uparrow ,\downarrow }{\hat{c}}_{x,\sigma }^{†}{\hat{c}}_{x+1,\sigma }+{\hat{c}}_{x+1,\sigma }^{†}{\hat{c}}_{x,\sigma }$ is electron hopping and ${\hat{V}}_{x,x+1}=({\hat{n}}_{x,\uparrow }-\frac{1}{2})({\hat{n}}_{x,\downarrow }-\frac{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 $u\ge 0$, featuring a $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 ().

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}}^z)$, where $\hat{N}={\sum }_{x=1}^L({\hat{c}}_{x,\uparrow }^{†}{\hat{c}}_{x,\uparrow }+{\hat{c}}_{x,\downarrow }^{†}{\hat{c}}_{x,\downarrow })$ is total electron charge, and ${\hat{S}}^z=\frac{1}{2}{\sum }_{x=1}^L({\hat{c}}_{x,\uparrow }^{†}{\hat{c}}_{i,\uparrow }-{\hat{c}}_{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).

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 $O\left(1/L\right)$ 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).

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

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*.