Hubbard-to-Heisenberg crossover (and efficient computation) of Drude weights at low temperatures

We illustrate how finite-temperature charge and thermal Drude weights of one-dimensional systems can be obtained from the relaxation of initial states featuring global (left–right) gradients in the chemical potential or temperature. The approach is tested for spinless interacting fermions as well as for the Fermi-Hubbard model, and the behavior in the vicinity of special points (such as half filling or isotropic chains) is discussed. We present technical details on how to implement the calculation in practice using the density matrix renormalization group and show that the non-equilibrium dynamics is often less demanding to simulate numerically and features simpler finite-time transients than the corresponding linear response current correlators; thus, new parameter regimes can become accessible. As an application, we determine the thermal Drude weight of the Hubbard model for temperatures T which are an order of magnitude smaller than those reached in the equilibrium approach. This allows us to demonstrate that at low T and half filling, thermal transport is successively governed by spin excitations and described quantitatively by the Bethe ansatz Drude weight of the Heisenberg chain.


Introduction
Computing correlation effects on static or dynamic transport properties at finite temperature such as charge or thermal conductivities [1,2] s w p d w s w = + n n n ( ) ( ) ( ) ( ) D 2 1 reg generally poses a daunting task for theorists. Even in linear response and at low energies (temperatures), the DC conductivity of gapless systems is usually not governed by the Luttinger liquid fixed point alone but influenced by the existence of conserved quantities [3]. In order to connect to actual experimental transport measurements on (quasi) 1D systems such as carbon nanotubes or strongly anisotropic 3D materials, it is thus essential to directly study microscopic Hamiltonians such as Heisenberg spin chains, spin ladders, or Hubbard models. This is a hard task even for seemingly 'simple', Bethe-ansatz solvable systems (such as the Heisenberg spin chain), because-similarly to correlation functions-transport coefficients are determined by couplings between all excitations. Whether or not a physical system exhibits dissipationless transport is signaled by the Drude weight D ν in equation (1). For ¹ n D 0, an initially excited current does not fully decay but will survive to infinite time. If for a given model the current operator I is conserved by the Hamiltonian, = [ ] I H , 0, transport is dissipationless at any temperature T. If I does not commute with H but has a finite overlap with (quasi-) local conserved quantities, dissipation processes are restricted and the Drude weight is non-zero; this can be shown strictly using the Mazur inequality [4][5][6]. While the question of dissipationless transport is mainly investigated for closed systems within linear response, it can also be studied in non-equilibrium setups [7][8][9][10] or for open quantum systems [11][12][13][14].
A more complex (and more experimentally relevant) system is the 1d Fermi-Hubbard model, which is again integrable via the Bethe ansatz and can thus in principle support dissipationless transport at finite temperature; however, neither the charge nor the energy current are fully conserved, . Most prior studies concentrated on the charge Drude weight D c [39][40][41][42][43][44][45][46], which is finite away from half filling by virtue of the Mazur inequality [6]. Directly at half filling, most works point towards [42,43,45,46], but this issue is not fully resolved yet. The thermal Drude weight of the Hubbard model has attracted far less attention: while the Mazur inequality [6] can again be used to show that > ( ) D T 0 th for arbitrary fillings, quantitative numbers for D th were only computed recently for large-to-intermediate temperatures [47]. It is one goal of this work to obtain the thermal Drude weight via the DMRG for temperatures which are an order of magnitude smaller and to demonstrate that one successively recovers the exact form of the Heisenberg chain's thermal Drude weight at low T and half filling. We can reach such small temperatures by extracting D th using a novel numerical protocol (see the next paragraph) which differs from the standard one employed in [47].
The 'standard route' to compute the Drude weight n ( ) D T numerically is provided by the linear response expression Equation (3) was discussed briefly in [48], and its validity was tested explicitly for the XXZ spin chain. The aim of the present paper is to expand on the ideas of [48], to study the practical relevance of equation (3) in more detail, and, as an application, to extract D th of the Fermi-Hubbard model at low T. After briefly introducing our methodology in sections 2 and 3, we extensively compare the real-time dynamics of equations (2) and (3) in section 4. One particular focus is on charge and thermal transport in the Hubbard model (which was not considered in [48]). We discuss the behavior in the vicinity of special points such as half filling or 'isotropic chains'. In section 5, practical aspects are presented on how to implement the calculation of equation (3) numerically. In particular, we document that equation (3) is often less demanding to simulate and features less complex finite-time transients than equation (2). As an application, we exploit this simplicity to determine the thermal Drude weight of the Hubbard model for temperatures which are an order of magnitude lower than those reached in the linear response calculation of [47], allowing us to access the regime where thermal transport is successively governed by spin excitations and described quantitatively by the exact ( ) D T th of the Heisenberg chain (section 6).

Model
The first model we consider describes spinless, interacting lattice fermions, whose Hamiltonian is given by The one-dimensional Fermi-Hubbard model is governed by

Density matrix renormalization group
In order to calculate the real time evolution of the one-dimensional quantum-mechanical systems introduced in equations (4) and (6), we employ the time-dependent [50][51][52][53][54] DMRG method [55][56][57] implemented using matrix product states [58][59][60][61]. Finite temperatures [62][63][64][65][66][67] are incorporated via a purification Yñ | T of the thermal density matrix r~e H T . The state Yñ | T can be obtained from the (known) Yñ ¥ | via an evolution e H T 2 in the inverse temperature [57]. Bothe H T 2 as well as the real time evolution operatore Ht i are factorized by a fourth order Trotter-Suzuki decomposition. We keep the discarded weight during each individual 'bond update' below a threshold value  . This leads to an exponential increase of the bond dimension χ during the real time evolution. In order to access time scales as large as possible, we employ the finite-temperature disentangler introduced in [68], which exploits the fact that purification is not unique to slow down the growth of χ. Our calculations are performed using a system size of the order of~( ) L O 100 sites. By comparing to other values of L, we have ensured that L is large enough for the results to be effectively in the thermodynamic limit [69].

Motivation of the non-equilibrium expression for the linear Drude weight
For reasons of completeness, we briefly motivate the origin of equation (3) for the charge case-more details can be found in [48]. Linear response theory predicts that local currents ( ) i x c are related to gradients of an applied potential m ( ) . The spatially integrated current flowing in a large but finite system (see below for comments on this issue) is then given by where dm is the total potential difference, and we have exploited that finite times serve as an infrared cutoff and regularize the δ-function via The ellipsis in equation (8) denotes a contribution from the regular part of the conductivity which can be neglected in the asymptotic limit of large times. Hence, equation (8) suggests that the total non-equilibrium current flowing in the presence of an initial potential gradient should increase linearly for large times and that the prefactor is determined by the Drude weight. If the regular contribution to the conductivity vanishes and transport is purely ballistic, the finite-time transients should vanish and linear behavior should manifest even for small t. We will explicitly verify this picture for the XXZ chain as well as for the Hubbard model.
One might wonder why for a fully ballistic system whose total current I commutes with the Hamiltonian, á ñ I is not constant but increases linearly with time. This confusion can be resolved by recapitulating the meaning of boundary conditions. In equation (8), we have implicitly assumed that our system has open boundaries, that the potential gradients occur in its center, and that the system size L is large enough so that at a time t the perturbations spreading out from the center have not yet reached the boundaries (so that the system is practically in the thermodynamic limit). Put differently, the global current I is effectively determined by the integral over a finite region, , whose size a fulfills   vt a L, with v being the Lieb-Robinson velocity. Importantly, 0 generally holds only for systems with periodic (not open) boundary conditions; the standard example is the energy current in the XXZ chain. If in our setup the left and right ends of the open system are connected, this creates a second, identical potential gradient, and the total current flowing in its vicinity is up to a minus sign identical to the one flowing in the center. Hence, the total current for a system with periodic boundary conditions is indeed constant. Note that this intuitive argument can be confirmed explicitly using the DMRG.
As a side remark, we note that if that Drude weight vanishes and transport is purely diffusive, the present setup contains information about the diffusion constant [10] via the second moment of the spatial profile of the local currents in equation (7). We leave a more thorough study for future work.

Numerical details
The linear response expression for the Drude weight, which is given by equation (2), can be simulated directly using the DMRG. It is advantageous to 'exploit time translation invariance' [67], c,th eq 1,2 and to carry out two independent calculations for n ( ) I t 2 as well asn ( ) I t 2 . Combined with the finitetemperature disentangler [68], this allows one to reach time scales which are roughly four times as large as the ones accessible by a 'standard' DMRG approach [64]. In principle, a similar 'trick' can be implemented when calculating the out-of-equilibrium expression in equation (3) [70]; however, this is not necessary for our purposes.
In order to compute the Drude weight via equation (3), we initialize the system in a state featuring a gradient in the temperature or chemical potential. The latter case is straightforward: we simply prepare the system using a thermal density matrix for the Hubbard model) on sites  l L 2 and > l L 2, respectively. Furthermore, one can distinguish the cases in which the central bond connecting sites L 2 and + L 2 1is cut (h L 2 set to zero) or not cut inH . The time evolution is then calculated using the original Hamiltonian (the potential is switched off).
The simplest way to compute the thermal Drude weight via equation (3) is to prepare the system in a state where r L R , are thermal density matrices of separated left and right systems (sites £ / l L 2 and > / l L 2), respectively. Their temperatures are chosen as In this setup, the bond between the sites L 2 and + L 2 1is naturally cut. This can be circumvented (which will turn out to be advantageous numerically; see below) by preparing the system using a density matrix

Comparison with linear response
In this section, we explicitly verify the validity of equation (3) for spinless fermions as well as for the Hubbard model and show that linear response Drude weights can indeed be obtained from the evolution of an out-ofequilibrium initial state featuring global gradients dm or db in the chemical potential or temperature, respectively. We discuss the finite-time dynamics of equations (2) and (3) and demonstrate that they exhibit different decay rates as one approaches special points of vanishing Drude weights.
In figure 1, we show the time evolution of the total charge and energy currents á at any Δ, I c is not fully conserved, but the charge Drude weight is finite for D < | | 1 at any > T 0 [27] and zero if D > | | 1. Hence, one expects that the non-equilibrium charge current grows linearly only for large times (since there is a non-vanishing regular contribution to the conductivity in equation (1)) but that á ñĨ t th for all t (see the discussion at the end of section 3 to resolve a potential confusion related to the choice of boundary conditions). This is indeed the case for all parameters that we studied (and illustrated explicitly in figure 1; this figure will be discussed in more detail in section 5).  T . In all cases, the long-time asymptotes of the linear response correlator and the out-of-equilibrium current agree. This confirms the validity of equation (3). At infinite temperature, an exact analytic solution for the charge Drude weight was constructed by Prosen [27,30] and is shown as a reference in figure 2(c); see also [19,20]. While the linear response correlators converge to the exact asymptote on a time scale which seems to be roughly independent of Δ, the currents á ñ m I t c decay more slowly (towards zero for D > 1 or towards the Prosen bounds for D < 1) as one approaches the isotropic point D = 1 from either side. This is interesting since it is still debated whether or not ( ) D T c is finite at D = 1 [28,29,31,38], and the out-of-equilibrium setup discussed in this paper might provide a new route to investigate this issue. If one fits the data at large times to an exponential function, one observes that the corresponding decay rate increases as one approaches D = 1. However, at the same time the quality of the fit  (4)). At time t=0, the chain is prepared in a state featuring a small, sharp gradient dm or db in the chemical potential or temperature between the left and right halves; we explicitly compare the cases in which the bond connecting the two halves is cut (or not cut) in the preparation of this state (see section 3 for details). Inset: evolution of the bond dimension during the different simulations.  (2) and (3), respectively. At infinite temperature, the exact solution constructed by Prosen for D < 1 [27,30] is shown as a reference; for finite T, we include Zotos' Bethe ansatz result from [19]. The behavior in the vicinity of the isotropic point D = 1 is discussed in the main text. worsens since the time scale accessible by the DMRG decreases (compare the curves at D = 0.707 and D = 0.901 in figure 2(c)). Analyzing this more quantitatively is not straightforward and left for future work.

Charge case
Next, we study charge transport in the Fermi-Hubbard model. While the Drude weight is finite away from half filling by virtue of the Mazur inequality [6], most works point towards = ( ) D T 0 c directly at half filling [42,43,45,46], but this issue is not finally resolved. Figure 3 shows the linear response correlators as well as the non-equilibrium currents for an on-site interaction of strength for three values of the filling. Both quantities converge to the same asymptotic value, which again validates equation (3). Moreover, we observe that the currents á ñ m I t c follow a simple exponential decay at large times, and sufficiently away from half filling one can more reliably determine D c by fitting to this form (see the dotted lines in figure 3). Interestingly, it seems that while the non-equilibrium currents decay more slowly as one approaches half filling, the linear response correlators do not exhibit a similar qualitative change but level off on comparable time scales. This scenario is analogous to what happens in the vicinity of D = 1 for spinless fermions and might again be used to gain further insights about the-still not fully resolved-issue of the charge Drude weight at half filling (future work).

Thermal case
Next, we turn to the thermal Drude weight. For spinless fermions, the energy current operator commutes with the Hamiltonian-transport is always purely ballistic. This is no longer the case in the Hubbard model, but the Mazur inequality can be used to prove that > ( ) D T 0 th for arbitrary fillings [6]. Hence, no subtleties occur at special points (in contrast to the charge case), and the asymptotic behavior of á ñ ( ) I t I th th eq and á ñ ( ) I t t T th can be determined straightforwardly. This is illustrated for two sets of parameters in figures 4(a), (b). We therefore do not present real-time data in more detail but directly discuss results for the Drude weight obtained via equations (2) and (3), respectively. T th th eq, 2 . In both models and for all temperatures and interactions, the Drude weight extracted using equation (3) agrees with the linear response prediction. This again confirms the validity of the non-equilibrium approach.

Final thoughts
If the integrability of the model at hand is broken, charge and thermal Drude weights become zero, and the nonequilibrium currents á ñ m ( ) I t t T c,th , decay to zero. We have verified this explicitly for charge and thermal transport in the Hubbard model in presence of an additional nearest-neighbor interaction V; representative results are presented in figure 4(c).
To summarize, we have shown that charge and thermal Drude weights can be obtained either from the linear response correlators using equation (2) or from out-of-equilibrium currents via equation (3). While both expressions yield the same asymptotic value, the finite-time transients do not necessarily agree. This becomes particularly obvious as one approaches special points of potentially vanishing Drude weights. Pragmatically, the non-equilibrium currents often exhibit a simpler (e.g., non-oscillatory) transient behavior (see, e.g., figure 4(a)), which renders it simpler to extract the Drude weight away from those special points.

Computational details
In this section, we present data for different initial states and illustrate how small the gradients in the chemical potential or temperature need to be chosen in practice in order to recover the linear response prediction. Moreover, we compare the numerical effort necessary to simulate equations (2) and (3), respectively. Figure 1 shows the charge and energy currents á ñ m I c and á ñ I T th for spinless fermions and two different initial states. The bond connecting the left and right regions (between which the initial gradients dm and db occur) is cut in one of them by setting = h 0 L 2 but left unchanged in the other (see section 3 for details on how the state is actually prepared). The currents feature the same asymptotic behavior in both cases, and even the finite-time transients (which appear in the charge case) are small. However, the numerical effort is drastically reduced if the bond is not cut in the preparation of the state (see the inset to figure 1), which one can understand intuitively from the fact that by setting = h 0 L 2 , one chooses an initial state which is further away from the stationary one. Hence, it is numerically advantageous to not 'cut the bond' in the preparation of the initial state, and all data in this work was obtained using this setup.  In figures 2(a) and 5(a), we explicitly show non-equilibrium data for spinless fermions calculated for different strength dm and db of the initial potential and temperature gradients. One can see that dm db, 0 .0 1is small enough to reproduce the linear response result with an accuracy that is beyond the resolution of the figure; deviations only occur for db = 0.1 in figure 5(a). All other data in this work was obtained using dm db~-, 0.01 0.1, and we checked (in representative cases) that decreasing the gradients even further does not influence the results.
It is instructive to recall that the order of limits dm  0 and  ¥ t in equation (3) is defined on an operational level: one first prepares a gradient dm and then time-evolves until the DMRG breaks down (at a finite time scale). This procedure is repeated with successively decreasing dm (starting from fairly large dm) until the results (on the accessible time scales) no longer change. This is illustrated in figure 2 for dm Î { } 0.1 , 0.003, 0.001 . Applying the real-and imaginary time evolution operatorse Ht i ande H T to a matrix product state involves singular value decompositions which lead to an increase of the bond dimension. The key approximation of the DMRG is to truncate this state by discarding the singular values below of a given threshold. The allowed discarded weight is the central parameter which controls the accuracy of the method.
In practice, we choose some representative sets of physical parameters and carry out calculations using different values of the discarded weight  ò during the real time evolution. An example is shown in figure 6(a), which displays the data of figure 2(a)  : we start from a large  1 and then lower this value succesively until the physical quantity at hand is computed with the desired accuracy. Note that (i) the bond dimension grows faster for smaller  1 , and hence the accessible time scales are reduced, and (ii) the linear response and non-equilibrium calculations are generally performed using a different  1 chosen such that the corresponding curves á ñ n n ( ) I t I eq and á ñ n m ( ) I t T , eventually reach the desired accuracy. In this work, the desired accuracy is set by the scale of each plot: in the case of figure 6(a), no deviation between the data calculated for   = 10 1 and   = 100 1 can be observed (on the scale of figure 6(a)); hence, the former value is a reasonable choice.
In figure 6(b), we illustrate how the bond dimension χ grows if the smallest value  100 1 is chosen as the discarded weight. We compare c ( ) t for the simulations of (i) the linear response expression á ñ ( ) I t I  1, n o t r a n s . i n v . see [70] for details). If translation invariance is exploited, the bond dimension χ at a time t is identical to c ( ) t 2 of the standard, single-time approach (1-time and 2-time curves in figure 6(b)); it still grows significantly faster than in the non-equilibrium approach. If translation invariance is not exploited in the linear response simulation, the growth of the bond dimension is comparable to the one of the non-equilibrium calculation. However, the former simulation is much more demanding, especially at low temperatures (we postpone arguments to the next paragraph). Hence, one can conclude that for this set of parameters the non-equilibrium calculation is the least computationally challenging one and can therefore be performed up to larger times. From a purely pragmatic standpoint, one should note that in order to obtain á ñ n m ( ) I t T , , one simply needs to time-evolve a state which is determined by the purification of the initial, nonequilibrium density matrix. In contrast, the linear response approach in its two-time version requires the calculation of a correlation function á -ñ n n ( ) ( ) I t I t 2 2 eq , which is more difficult to implement numerically.