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 behaviour 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 which are an order of magnitude smaller than those reached in the equilibrium approach. This allows us to demonstrate that at low and half filling, thermal transport is successively governed by spin excitations and described quantitatively by the Bethe ansatz Drude weight of the Heisenberg chain.
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 . 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 in Eq. (1). For , an initially excited current does not fully decay but will survive to infinite time. If for a given model the current operator is conserved by the Hamiltonian, , transport is dissipationless at any temperature . If does not commute with 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].
One prototypical low-dimensional model is given by spinless, interacting fermions (equivalently, a Heisenberg XXZ spin chain); it can be diagonalized exactly using the Bethe ansatz [15, 16] and possesses an infinite number of local conserved quantities. The energy current operator commutes with the Hamiltonian, and the corresponding Drude weight was computed analytically [17, 18]. The charge Drude weight has attracted considerable attention [19, 20, 21, 22, 23, 24, 25, 26], in particular at half filling where the current operator has no overlap with any of the standard local conserved quantities. Considerable progress has been made within the recent years [3, 10, 27, 28, 29, 30, 31, 32, 33, 34, 35]. In particular, Prosen constructed quasi-local conserved quantities [27, 30, 36] to show analytically that is finite throughout the gapless phase (excluding the isotropic point); quantitative numbers can be obtained, e.g., using the real-time density matrix renormalization group (DMRG) [29, 37] or dynamical typicality . Whether or not the Drude weight is finite for an isotropic chain is still debated [28, 29, 31, 38].
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 [39, 40, 41, 42, 43, 44, 45, 46], which is finite away from half filling by virtue of the Mazur inequality . 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  can again be used to show that for arbitrary fillings, quantitative numbers for were only computed recently for large-to-intermediate temperatures . 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 and half filling. We can reach such small temperatures by extracting using a novel numerical protocol (see the next paragraph) which differs from the standard one employed in Ref. .
The ‘standard route’ to compute the Drude weight numerically is provided by the linear response expression
where the real-time current correlation function can be obtained directly using the DMRG. It was recently demonstrated  that can alternatively be calculated from the non-equilibrium current flowing in the presence of a small chemical potential or temperature gradient via 
Eq. (3) was discussed briefly in Ref. , and its validity was tested explicitly for the XXZ spin chain. The aim of the present paper is to expand on the ideas of Ref. , to study the practical relevance of Eq. (3) in more detail, and, as an application, to extract of the Fermi-Hubbard model at low .
After briefly introducing our methodology in Secs. 2 and 3, we extensively compare the real-time dynamics of Eqs. (2) and (3) in Sec. 4. One particular focus is on charge and thermal transport in the Hubbard model (which was not considered in Ref. ). We discuss the behaviour in the vicinity of special points such as half filling or ‘isotropic chains’. In Sec. 5, practical aspects are presented on how to implement the calculation of Eq. (3) numerically. In particular, we document that Eq. (3) is often less demanding to simulate and features less complex finite-time transients than Eq. (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 Ref. , allowing us to access the regime where thermal transport is successively governed by spin excitations and described quantitatively by the exact of the Heisenberg chain (Sec. 6).
2 Model and Method
The first model we consider describes spinless, interacting lattice fermions, whose Hamiltonian is given by
with being a fermionic annihilation operator acting on site , , and denoting the hopping amplitude and nearest-neighbor interaction strength, respectively. Eq. (4) can be mapped to an XXZ spin chain with an exchange coupling and anisotropy via a Jordan-Wigner transformation. The charge and energy current of this model take the standard form
The one-dimensional Fermi-Hubbard model is governed by
where annihilates a fermion with spin on site , and . The interaction strength and the hopping matrix element are denoted by and , respectively. The charge and energy current are given by
2.2 Density Matrix Renormalization Group
In order to calculate the real time evolution of the one-dimensional quantum-mechanical systems introduced in Eqs. (4) and (6), we employ the time-dependent [50, 51, 52, 53, 54] density matrix renormalization group 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 of the thermal density matrix . The state can be obtained from the (known) via an evolution in the inverse temperature . Both as well as the real time evolution operator 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 Ref. , 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 sites. By comparing to other values of , we have ensured that is large enough for the results to be effectively in the thermodynamic limit .
3 Computation of the Drude weight
3.1 Motivation of the non-equilibrium expression for the linear Drude weight
For reasons of completeness, we briefly motivate the origin of Eq. (3) for the charge case – more details can be found in Ref. . Linear response theory predicts that local currents are related to gradients of an applied potential via . The spatially integrated current flowing in a large but finite system (see below for comments on this issue) is then given by
where 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 Eq. (8) denotes a contribution from the regular part of the conductivity which can be neglected in the asymptotic limit of large times.
Hence, Eq. (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 behaviour should manifest even for small . 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 commutes with the Hamiltonian, is not constant but increases linearly with time. This confusion can be resolved by recapitulating the meaning of boundary conditions. In Eq. (8), we have implicitly assumed that our system has open boundaries, that the potential gradients occur in its center, and that the system size is large enough so that at a time 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 is effectively determined by the integral over a finite region, , whose size fulfills , with being the Lieb-Robinson velocity. Importantly, 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.
3.2 Numerical details
and to carry out two independent calculations for as well as . Combined with the finite-temperature disentangler , this allows one to reach time scales which are roughly four times as large as the ones accessible by a ‘standard’ DMRG approach . In principle, a similar ‘trick’ can be implemented when calculating the out-of-equilibrium expression in Eq. (3) ; however, this is not necessary for our purposes.
In order to compute the Drude weight via Eq. (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
where is the Hamiltonian of Eq. (4) or Eq. (6) complemented by a term for the XXZ chain [and similarly for the Hubbard model] on sites and , respectively. Furthermore, one can distinguish the cases in which the central bond connecting sites and is cut ( set to zero) or not cut in . 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 Eq. (3) is to prepare the system in a state
where are thermal density matrices of separated left and right systems (sites and ), respectively. Their temperatures are chosen as
In this setup, the bond between the sites and is naturally cut. This can be circumvented (which will turn out to be advantageous numerically; see below) by preparing the system using a density matrix
4 Comparison with linear response
In this section, we explicitly verify the validity of Eq. (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-of-equilibrium initial state featuring global gradients or in the chemical potential or temperature, respectively. We discuss the finite-time dynamics of Eqs. (2) and (3) and demonstrate that they exhibit different decay rates as one approaches special points of vanishing Drude weights.
In Fig. 1, we show the time evolution of the total charge and energy currents and for spinless interacting fermions. While at any , is not fully conserved, but the charge Drude weight is finite for at any  and zero if . 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 Eq. (1)] but that for all (see the discussion at the end of Sec. 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 Fig. 1; this Figure will be discussed in more detail in Sec. 5).
4.1 Charge case
We now explicitly compare the finite-time behaviour of the linear response charge current correlation function and the non-equilibrium current induced by an initial gradient in the chemical potential. The large- asymptote of both quantities determines the Drude weight via Eqs. (2) and (3), respectively.
Fig. 2 shows data for spinless fermions at half filling and (a) at , (b) at , and (c) at . In all cases, the long-time asymptotes of the linear response correlator and the out-of-equilibrium current agree. This confirms the validity of Eq. (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 Fig. 2(c); see also Refs. [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 decay more slowly (towards zero for or towards the Prosen bounds for ) as one approaches the isotropic point from either side. This is interesting since it is still debated whether or not is finite at [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 . However, at the same time the quality of the fit worsens since the time scale accessible by the DMRG decreases [compare the curves at and in Fig. 2(c)]. Analysing this more quantitatively is not straightforward and left for future work.
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 , most works point towards directly at half filling [42, 43, 45, 46], but this issue is not finally resolved. Fig. 3 shows the linear response correlators as well as the non-equilibrium currents for an on-site interaction of strength at a temperature of for three values of the filling. Both quantities converge to the same asymptotic value, which again validates Eq. (3). Moreover, we observe that the currents follow a simple exponential decay at large times, and sufficiently away from half filling one can more reliably determine by fitting to this form (see the dotted lines in Fig. 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 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).
4.2 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 for arbitrary fillings . Hence, no subtleties occur at special points (in contrast to the charge case), and the asymptotic behaviour of and can be determined straightforwardly. This is illustrated for two sets of parameters in Fig. 4(a,b). We therefore do not present real-time data in more detail but directly discuss results for the Drude weight obtained via Eqs. (2) and (3), respectively.
Fig. 5 shows linear response and out-of-equilibrium data for as a function of the temperature for (a) spinless fermions with , and (b) the Hubbard model with , both at half filling. In (a), we plot the exact Bethe ansatz solution  for comparison; note that in (b), the point can be solved analytically, and we show the exact linear response result instead of the DMRG data. The high- asymptote (dashed lines) displays . In both models and for all temperatures and interactions, the Drude weight extracted using Eq. (3) agrees with the linear response prediction. This again confirms the validity of the non-equilibrium approach.
4.3 Final thoughts
If the integrability of the model at hand is broken, charge and thermal Drude weights become zero, and the non-equilibrium currents 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 ; representative results are presented in Fig. 4(c).
To summarize, we have shown that charge and thermal Drude weights can be obtained either from the linear response correlators using Eq. (2) or from out-of-equilibrium currents via Eq. (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 behaviour [see, e.g., Fig. 4(a)], which renders it simpler to extract the Drude weight away from those special points.
5 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 Eqs. (2) and (3), respectively.
Fig. 1 shows the charge and energy currents and for spinless fermions and two different initial states. The bond connecting the left and right regions (between which the initial gradients and occur) is cut in one of them by setting but left unchanged in the other (see Sec. 3 for details on how the state is actually prepared). The currents feature the same asymptotic behaviour 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 Fig. 1), which one can understand intuitively from the fact that by setting , 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 Figs. 2(a) and 5(a), we explicitly show non-equilibrium data for spinless fermions calculated for different strength and of the initial potential and temperature gradients. One can see that is small enough to reproduce the linear response result with an accuracy that is beyond the resolution of the Figure; deviations only occur for in Fig. 5(a). All other data in this work was obtained using , 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 and in Eq. (3) is defined on an operational level: One first prepares a gradient and then time-evolves until the DMRG breaks down (at a finite time scale). This procedure is repeated with successively decreasing (starting from fairly large ) until the results (on the accessible time scales) no longer change. This is illustrated in Fig. 2 for .
Applying the real- and imaginary time evolution operators and 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 Fig. 6(a), which displays the data of Fig. 2(a) for three different, successively decreasing : We start from a large 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 , and hence the accessible time scales are reduced, and (ii) the linear response and non-equilibrium calculations are generally performed using a different chosen such that the corresponding curves and eventually reach the desired accuracy. In this work, the desired accuracy is set by the scale of each plot: In the case of Fig. 6(a), no deviation between the data calculated for and can be observed [on the scale of Fig. 6(a)]; hence, the former value is a reasonable choice.
In Fig. 6(b), we illustrate how the bond dimension grows if the smallest value is chosen as the discarded weight. We compare for the simulations of (i) the linear response expression calculated in the standard way from a single time evolution, (ii) the same, but writing and carrying out two individual time evolutions for , and (iii) the non-equilibrium approach . The fastest growth occurs in (i). Using (ii), one can access a time scale which is (roughly) twice as large at the same computational cost. More precisely, translation invariance (in space) can only be exploited in one of the ; their calculations thus exhibit different (and are also performed using different individual discarded weights ; see Ref.  for details). If translation invariance is exploited, the bond dimension at a time is identical to of the standard, single-time approach [1-time and 2-time curves in Fig. 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 , one simply needs to time-evolve a state which is determined by the purification of the initial, non-equilibrium density matrix. In contrast, the linear response approach in its two-time version requires the calculation of a correlation function , which is more difficult to implement numerically. Extracting Drude weights via Eq. (3) thus seems to be a viable alternative to the standard linear response route.
We conclude this Section with a few more technical remarks; additional details can again be found in Ref. . If one does or does not exploit translation invariance in the linear response approach, the calculation of amounts to time-evolving locally or globally quenched states or , respectively. The non-equilibrium setup always corresponds to time-evolving a locally quenched state. In local quenches, perturbations spread with a finite Lieb-Robinson velocity, and the bond dimension does not increase significantly outside of this ‘light cone’. This is one reason why the linear response calculation is more demanding if one cannot use translation invariance. Moreover, one needs to perform the time evolution of , which requires the application of a global operator to the state . This increases the bond dimension instantaneously by a factor which is determined by the matrix product operator representation of ; all additional symmetries (such as spin-flip symmetry) should hence be exploited in this simulation and not in . Since finite temperatures are reached via an evolution in starting from , grows with decreasing . In practice, for the Hubbard model at moderately low , can reach values of , and applying the global energy current operator to and subsequently computing its time evolution becomes no longer feasible.
6 Hubbard model: thermal Drude weight at low
In this Section, we revisit the realm of the thermal Drude weight of the Fermi-Hubbard model. As mentioned above, the Mazur inequality  stipulates that is finite at any filling, and quantitative values were recently obtained  using the linear response expression of Eq. (2). In Fig. 5(b), we explicitly demonstrated that the Drude weight extracted from the time evolution of an initial, small temperature gradient via Eq. (3) coincides with the linear response prediction.
We can now exploit the computational simplicity of the non-equilibrium approach as well as the fact its finite-time transients have a simpler form (see Fig. 4) to determine the thermal Drude weight for temperatures which are an order of magnitude lower than those reached in Ref. . The results are shown in Fig. 7 at half filling and for two values of the on-site interaction. At high temperatures, increases quadratically, becomes maximal when reaches the charge gap, and subsequently decreases since charge excitations are frozen out. One expects that at low , transport is governed by spin excitations whose dynamics are described by a Heisenberg spin chain [or equivalently, Eq. (4)] with an exchange coupling of strength . In other words, one expects to recover a second peak in at low whose form quantitatively follows the exact (Bethe ansatz or DMRG) Drude weight of the Heisenberg chain [the curve at in Fig. 5(a) with units rescaled to ]. This is indeed the case. To the best of our knowledge, this two-peak structure constitutes the first quantitative observation of a full Hubbard-to-Heisenberg crossover for a transport quantity within the one-dimensional Fermi-Hubbard model at finite temperature.
In this paper, we have investigated how the linear response charge and thermal Drude weights of integrable one-dimensional systems can be computed from the relaxation of initial states featuring small gradients in the chemical potential or temperature. Using density matrix renormalization group numerics for spinless fermions as well as for the Hubbard model, we extensively compared the real-time dynamics of the currents flowing in this non-equilibrium setup with the linear response correlators . While both quantities determine in the limit , the finite-time behaviour differs. Only seems to exhibit diverging decay rates in the vicinity of points where it is still debated whether or not the Drude weight vanishes; we explicitly demonstrated this for charge transport in an XXZ spin chain near by comparing with Prosen’s exact solution. Away from such special points, the non-equilibrium currents often exhibit simpler (e.g., non-oscillatory) transients and are less demanding to compute numerically. We exploited this to extract the thermal Drude weight of the Hubbard model for temperatures which are an order of magnitude lower than those reached in the linear response approach. At half filling and sufficiently large on-site interactions, features a two-peak structure and at low is quantitatively described by the exact Bethe ansatz Drude weight of the Heisenberg spin chain.
It would be interesting to generalize our approach in order to efficiently extract transport properties beyond the Drude weight (such as the low-frequency behaviour of the regular part of the conductivity). Moreover, the vicinity of special points (e.g., isotropic XXZ chains) certainly deserves further attention.
I thank F. Heidrich-Meisner and D. M. Kennes for discussions and acknowledge support by the DFG through the Emmy Noether program (KA 3360/2-1).
-  Mahan G D 1990 Many-Particle Physics (New York London: Plenum Press)
-  Luttinger J M 1964 135 A1505
-  Sirker J, Pereira R G and Affleck I 2009 Phys. Rev. Lett. 103 216602
-  Suzuki M 1971 Physica 51 277
-  Mazur P 1969 Physica (Amsterdam) 43 533
-  Zotos X, Naef F and Prelovšek P 1997 55 11029
-  Langer S, Heidrich-Meisner F, Gemmer J, McCulloch I and Schollwöck U 2009 Phys. Rev. B 79 214409
-  Jesenko S and Žnidarič M 2011 84 174438
-  Langer S, Heyl M, McCulloch I P and Heidrich-Meisner F 2011 Phys. Rev. B 84 205115
-  Karrasch C, Moore J E and Heidrich-Meisner F 2014 Phys. Rev. B 89(7) 075139
-  Prosen T and Žnidarič M 2009 J. Stat. Mech 2009 P02035
-  Žnidarič M 2011 Phys. Rev. Lett. 106 220601
-  Žnidarič M 2013 Phys. Rev. B 88(20) 205135
-  Mendoza-Arenas J J, Al-Assam S, Clark S R and Jaksch D 2013 J. Stat. Mech. P07007
-  Essler F, Frahm H, Göhmann F, Klümper A and V E Korepin V E 2005 The one-dimensional Hubbard model (Cambride University Press)
-  Giamarchi T 2004 Quantum Physics in One Dimension (Oxford: Clarendon Press)
-  Klümper A and Sakai K 2002 J. Phys. A 35 2173
-  Sakai K and Klümper A 2003 J. Phys. A 36 11617
-  Zotos X 1999 Phys. Rev. Lett. 82 1764
-  Benz J, Fukui T, Klümper A and Scheeren C 2005 J. Phys. Soc. Jpn. Suppl. 74 181
-  Jung P and Rosch A 2007 Phys. Rev. B 76 245108
-  Narozhny B N, Millis A J and Andrei N 1998 Phys. Rev. B 58 R2921–R2924
-  Heidrich-Meisner F, Honecker A, Cabra D C and Brenig W 2003 Phys. Rev. B 68 134436
-  Mukerjee S and Shastry B S 2008 Phys. Rev. B 77 245131 (pages 7)
-  Alvarez J V and Gros C 2002 88 077203
-  Heidarian D and Sorella S 2007 Phys. Rev. B 75 241104(R)
-  Prosen T 2011 Phys. Rev. Lett. 106 217206
-  Herbrych J, Prelovšek P and Zotos X 2011 Phys. Rev. B 84 155125
-  Karrasch C, Hauschild J, Langer S and Heidrich-Meisner F 2013 Phys. Rev. B 87 245128
-  Prosen T and Ilievski E 2013 Phys. Rev. Lett. 111(5) 057203
-  Steinigeweg R, Gemmer J and Brenig W 2014 Phys. Rev. Lett. 112 120601
-  Prelovšek P, ElShawish S, Zotos X and Long M W 2004 Phys. Rev. B 70 205129
-  Steinigeweg R and Gemmer J 2009 Phys. Rev. B 80 184402
-  Sirker J, Pereira R G and Affleck I 2011 Phys. Rev. B 83 035115
-  Steinigeweg R and Brenig W 2012 Phys. Rev. Lett. 107 250602
-  Mierzejewski M, Prelovšek P and Prosen T 2015 Phys. Rev. Lett. 114 140601
-  Karrasch C, Bardarson J and Moore J E 2012 Phys. Rev. Lett. 108 227206
-  Carmelo J M P, Prosen T and Campbell D K 2015 Phys. Rev. B 92(16) 165133
-  Fujimoto S and Kawakami N 1998 Journal of Physics A: Mathematical and General 31 465
-  Kirchner S, Evertz H G and Hanke W 1999 Phys. Rev. B 59 1825
-  Peres N M R, Dias R G, Sacramento P D and Carmelo J M P 2000 Phys. Rev. B 61 5169–5183
-  Carmelo J, Gu S J and Sacramento P 2013 Annals of Physics 339 484
-  Prosen T and Žnidarič M 2012 Phys. Rev. B 86(12) 125118
-  Prosen T 2014 Phys. Rev. Lett. 112 030603
-  Karrasch C, Kennes D M and Moore J E 2014 Phys. Rev. B 90 155104
-  Jin F, Steinigeweg R, Heidrich-Meisner F, Michielsen K and De Raedt H 2015 Phys. Rev. B 92 205103
-  Karrasch C, Kennes D M and Heidrich-Meisner F 2016 Phys. Rev. Lett. 117(11) 116401
-  Vasseur R, Karrasch C and Moore J E 2015 Phys. Rev. Lett. 115(26) 267201
-  Ilievski E and Prosen T 2013 Communications in Mathematical Physics 318 809–830
-  Vidal G 2004 Phys. Rev. Lett. 93 040502
-  White S R and Feiguin A E 2004 Phys. Rev. Lett. 93 076401
-  Daley A, Kollath C, Schollwöck U and Vidal G 2004 J. Stat. Mech.: Theory Exp. 2004 P04005
-  Schmitteckert P 2004 Phys. Rev. B. 70 121302(R)
-  Vidal G 2007 Phys. Rev. Lett. 98 070201
-  White S R 1992 Phys. Rev. Lett. 69 2863
-  Schollwöck U 2005 Rev. Mod. Phys. 77 259
-  Schollwöck U 2011 Ann. Phys. (NY) 326 96
-  Fannes M, Nachtergaele B and Werner R F 1991 J. Phys. A: Math. Gen. 24 L185
-  Östlund S and Rommer S 1995 Phys. Rev. Lett. 75 3537–3540
-  Verstraete F and Cirac J I 2006 Phys. Rev. B 73 094423
-  Verstraete F, Cirac J I and Murg V 2008 Adv. Phys. 57 143
-  Verstraete F, Garcia-Ripoll J J and Cirac J I 2004 Phys. Rev. Lett. 93 207204
-  White S 2009 Phys. Rev. Lett. 102 190601
-  Barthel T, Schollwöck U and White S R 2009 Phys. Rev. B 79 245101
-  Zwolak M and Vidal G 2004 Phys. Rev. Lett. 93 207205
-  Sirker J and Klümper A 2005 Phys. Rev. B 71 241101(R)
-  Barthel T 2013 New J. Phys. 15 073010
-  Karrasch C, Bardarson J and Moore J E 2012 108 227206
-  Karrasch C, Hauschild J, Langer S and Heidrich-Meisner F 2013 Phys. Rev. B 87 245128
-  Kennes D and Karrasch C 2016 Computer Physics Communications 200 37 – 43 ISSN 0010-4655