# Transport properties of the one-dimensional Hubbard model at finite temperature

###### Abstract

We study finite-temperature transport properties of the one-dimensional Hubbard model using the density matrix renormalization group. Our aim is two-fold: First, we compute both the charge and the spin current correlation function of the integrable model at half filling. The former decays rapidly, implying that the corresponding Drude weight is either zero or very small. Second, we calculate the optical charge conductivity in presence of small integrability-breaking next-nearest neighbor interactions (the extended Hubbard model). The DC conductivity is finite and diverges as the temperature is decreased below the gap. Our results thus suggest that the half-filled, gapped Hubbard model is a normal charge conductor at finite temperatures. As a testbed for our numerics, we compute for the integrable XXZ spin chain in its gapped phase.

###### pacs:

71.10.Fd, 71.10.Pm, 71.27.+a## I Introduction

The physics of one-dimensional (1D) systems is strongly influenced by electronic correlations. E.g., the low-energy behavior of a large class of 1D models is not described by a Fermi liquid but exhibits bosonic excitations. Thermodynamic properties or equilibrium correlation functions of this so-called Luttinger liquid can be obtained elegantly using field theory. Transport properties, however, are usually not governed by the low-energy Luttinger liquid fixed point but by an interplay between dangerously irrelevant operators scattering the currents and conserved quantities protecting them.andrei (); xxz_sirker (); xxz_prosen () 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 study generic microscopic models.

The Hubbard modelhubbook () plays a fundamental role in the physics of correlated electrons as it provides the most transparent realization of Mott physics: it describes spinful electrons moving on a lattice with on-site interaction. From an experimental perspective, the Hubbard model has been employed to describe a wide variety of crystals (both insulating and conducting); in one dimension, it is also used as a starting point for polymers. Despite the Hubbard model’s apparent simplicity and the existence of an exact solution for its thermodynamics in 1D, little is known about its transport properties. Above one dimension even the phase diagram of the Hubbard model with doping is a matter of debate, including the regime of greatest interest where charged excitations are gapped while the neutral sector includes gapless spin excitations.

From a theoretical perspective, it is less demanding (both computationally and analytically) to study a spinless fermion model with nearest-neighbor interactions, which has therefore been investigated extensively despite the fact that it is of less relevance for experimental setups. Even though the model (which can be mapped to a XXZ spin chain via a Jordan-Wigner transformation) allows for a Bethe ansatz solution which is ‘simpler’ than that of the Hubbard model, extracting transport coefficients, which are determined by couplings between all excitations, remains a formidable task. Over the last decades a significant number of worksbethe0 (); xxz_betheT1 (); xxz_betheT2 (); xxz_qmcgros (); xxz_edmillis (); xxz_qmcsorella (); xxz_fieldtheory (); xxz_fabian (); xxz_rosch (); xxz_edherbrych (); xxz_sirker () investigated the equilibrium transport properties of the XXZ chain, but it was only recently proven rigorouslyxxz_prosen (); xxz_prosen2 () that the half-filled system can support dissipationless currents at finite temperature ; this corresponds to a finite Drude contribution to the conductivity

(1) |

Numerical values for the Drude weight can be obtained, e.g., via the density matrix renormalization group.drudepaper (); drudepaper2 ()

It is the main goal of this paper to investigate finite-temperature linear-response transport properties of the Hubbard model, about which comparably little is known (we will give a more detailed overview of previous works in Sec. III). In particular, we employ the density matrix renormalization groupwhite1 (); dmrgrev (); dmrgrev2 () to compute real time charge- and spin current correlation functions whose Fourier transform determines . Our focus is two-fold: First, we demonstrate that the charge current correlators at half filling decay rapidly, suggesting that the charge Drude weight is either zero or very small at any finite temperature ( is known to vanish at ). Second, we study the optical conductivity in presence of small next-nearest neighbor interactions which break integrabilityintegrab1 (); integrab2 () – and thus eliminate a potentially small, finite Drude weight – but do not trigger a phase transition.extendedhub () As the temperature is lowered below the charge gap, successively develops a sharp increase at the optical absorption threshold. More importantly, the DC conductivity is finite and diverges approximately as (our data is insufficient to rule out a weakly -dependent prefactor). This analysis suggests that the half-filled, gapped Hubbard model is a normal charge conductor at finite temperatures. As a testbed for our numerics, we also study the integrable XXZ chain in its gapped phase.diff_zotos (); diff_prosenznidaric (); diff_steinigeweggemmer (); diff_steinigeweg1 (); diff_znidaric (); diff_znidaric2 (); diff_steinigeweg2 (); finiteTquench (); diff_typicality () The DC conductivity again grows monotonously for all temperatures from to far below the gap.

## Ii Model and Method

Models — The prime interest of this work is the extended Hubbard model governed by

(2) |

with , and being fermionic annihilation operators. We solely focus on the case of half filling and zero magnetic field (the reason for this will be outlined below). At , the model is symmetric under if charge and spin degrees of freedom are interchanged, and we thus stick to repulsive interactions only. The Hubbard model is integrable via Bethe ansatz for . A charge gap opens for while the spin sector remains gapless.hubbook ()

The second model we study is the XXZ chain defined by

(3) |

where are spin- operators. The model is integrable via Bethe ansatz; a gap opens for . At and , the Hubbard model can be mapped to an isotropic, antiferromagnetic XXZ chain ().

Transport coefficients — Both the Drude weight and the regular part of the conductivity of Eq. (1) can be obtained from the current correlation function,

(4) |

and

(5) |

where a potentially finite Drude weight has been subtracted in . Only finite times can be reached in the DMRG and it is thus imperative to estimate the associated error of . In the DC limit where ‘finite-time’ effects are largest, only contributes to Eq. (5), and one can therefore estimate the error by comparing with obtained from ,

(6) |

which follows from a Kramers-Kronig relation.

The current operator is defined via a continuity equation. The local charge and spin current of the Hubbard model read

(7) |

and for the XXZ chain one finds

(8) |

DMRG — We compute the current correlation function

(9) |

using the time-dependenttdmrg1 (); tdmrg2 (); tdmrg3 (); tdmrg4 (); tdmrg5 (); tdmrg6 () density matrix renormalization groupwhite1 (); dmrgrev (); dmrgrev2 () in a matrix product statemps1 (); mps2 (); mps3 (); mps4 () implementation. Finite temperaturesdmrgT (); barthel (); verstraete (); vidalop (); tmrg1 (); metts (); trick2a () are incorporated via purification of the thermal density matrix. The real- and imaginary time evolution operators in Eq. (9) are factorized by a fourth order Trotter decomposition with a step size of . The discarded weight during each individual ‘bond update’ is kept beneath a threshold value , which is the key parameter controlling the accuracy of the simulation. The system size, however, can easily be chosen large enough for the results to be effectively in the thermodynamic limit. We carry out most of our calculations for or and exemplary compare against other values. The dependence of the numerical data on and is illustrated in the insets to Figs. 1 and 2(b).

The bond dimension increases exponentially fast during the real time evolutions, and the simulation is stopped once numerical resources are exhausted. We pursue several strategies in order to access time scales as large as possible. First, we employ the finite-temperature disentangler introduced in Ref. drudepaper, , which uses the fact that purification is not unique to slow down the growth of . Second, we ‘exploit time translation invariance’,trick2a () rewrite , and carry out two independent calculations for as well as . This allows to access time scales roughly twice as large. Third, for the XXZ chain we recast with , and similarly for the charge and spin currents of the Hubbard model. We exploit symmetries (e.g., both charge and spin conservation) during the time evolution and stop the calculation once the bond dimension reaches for the local and for the global (, in some exemplary cases).

## Iii Drude weight

We start this section with a brief summary of what is know about the Drude weight of the Hubbard model. For (free fermions), both the spin and charge currents are conserved by , and thus trivially . At but , the model is integrable via Bethe ansatz,hub0 () and away from half filling it follows from the Mazur inequality that the Drude weights are finite at any temperature.integrab2 () At half filling, however, all known local conserved quantitieshub1 () have zero overlap with the current operators for symmetry reasons. At one can use Bethe ansatzbethe0 () to show that the charge Drude weight vanishes while the spin Drude weight is finite (recall that the charge sector is gapped while the spin sector is gapless). Further studies of ground state transport properties can be found in Refs. hubT0a, ; hubT0b, ; hubT0c, ; hubT0d, ; hubT0e, ; hubT0f, ; hubT0g, ; hubT0h, ; hubT0i, ; hubT0j, ; hubT0k, ; hubT0l, ; hubT0m, ; hubT0n, .

The most interesting situation is thus the half-filled Hubbard model at finite temperature, where there is still controversy about whether or not the Drude weight is finite. In Refs. integrab1, ; integrab2, it was conjectured that for an integrable model if and only if , which holds true for the XXZ chain whose finite- Drude weight is non-zero in the gapless phase but vanishes in the gapped phase. For the Hubbard model, Bethe ansatz resultsbetheThub () suggest that despite the charge gap. However, this calculation is very involved: two different waysxxz_betheT1 (); xxz_betheT2 () to approximately solve the Bethe ansatz equations for the XXZ chain yield results which disagree, and one might expect the situation for the Hubbard model to be similarly subtle. While quantum Monte Carlo (QMC) numericsqmchub () and an analysis of low-energy excitationsexhub () support a finite Drude weight, exact diagonalization (ED) data,edhub () large- analytics,largeUhub () and symmetry argumentssymhub () seems to favor . Non-equilibrium DMRG calculations (at only) suggest a zero Drude weight but a finite diffusion constant .prosenhub () The situation for the Hubbard model is thus completely analogous to the gapless XXZ chain where – after decades of dispute – it was only recently shown rigorously that the Drude weight at half filling and is nonzero at finite .xxz_prosen (); xxz_prosen2 (); xxz_prosenD1 ()

The controversial status of the Drude weight of the half-filled Hubbard model requires further work in this direction. The DMRG allows to obtain results directly on the real time axis (in contrast to QMC), and the thermodynamic limit can be accessed easily (in contrast to ED). Its drawback is that only finite time scales can be reached. We show DMRG results for the current correlation functions of the half-filled, integrable Hubbard model at (where spin and charge are symmetric) in Fig. 1. Finite temperatures are shown in Fig. 2. The charge current correlation functions fall off rapidly at any , suggesting that the charge Drude weight vanishes at finite temperatures in agreement with Refs. edhub, ; largeUhub, ; symhub, but in contrast to the Bethe ansatzbetheThub () and QMCqmchub () predictions. However, the times reached in our numerics are too small to unambiguously rule out a small but finite . It is thus instructive to establish an upper bound to as a reference for future works. Such a bound can be defined from the value of at the largest time in case that it falls off monotonously or from the value at the last maximum in case oscillates (this is illustrated by the dashed lines in Fig. 1). The upper bound for the charge Drude weight at finite and is shown in the inset to Fig. 2(a). It is significantly smaller than the values estimated by QMC.qmchub ()

In contrast, the spin current correlation functions decay on an increasingly larger time scales as the temperature is decreased. This is consistent with the spin Drude weight being nonzero at .bethe0 () At infinite temperature, however, spin and charge degrees of freedom are symmetric, . If we assume that is zero at all , this leaves two possible scenarios for the temperature dependence of : either is nonzero at any finite but vanishes in the limit of faster than , or . Our numerical data is insufficient to answer this question conclusively. For reasons of completeness, an upper bound to is shown in the inset to Fig. 2(a).

## Iv Optical Conductivity

XXZ chain — In this section we investigate the regular part of the optical conductivity. It is instructive to first study the simpler case of the XXZ chain in its gapped phase , particularly in order to illustrate how to assess the error when is computed from real time DMRG data. Fig. 3 again shows that the current correlation decays to zero for and that thus the Drude weight vanishes. The DC conductivity is determined by the integral of via Eq. (5). For anisotropies where the spectral gap is of order one, we can simulate up to times where this integral can be computed without having to resort to any extrapolation algorithms for all temperatures from to (we explain below how to estimate the error). This is illustrated in Fig. 4(a) where and . Interestingly, diverges even as is decreased below the gap [see the inset to Fig. 4(b)]. The temperature-dependence is approximately , but we cannot rule out an additional prefactor that varies weakly with . Such a divergence is consistent with a semiclassical analysis by Damle and Sachdevgaparg (); gaparg2 () who show that below the gap an exponentially small quasiparticle density is compensated by an exponentially long life time; for a model with spin- symmetry; this yields . However, to the best of our knowledge this picture was never confirmed conclusively for a microscopic model and all temperatures ranging from to . Previous studies of the diffusion constant, which is related to via an Einstein relation, can be found in Refs. diff_prosenznidaric, ; diff_steinigeweggemmer, ; diff_znidaric, ; diff_znidaric2, ; diff_steinigeweg2, ; finiteTquench, for as well as in Refs. diff_steinigeweg1, ; diff_znidaric2, for finite but large (compared to the gap) temperatures. Finally, we show the full frequency-dependent AC conductivity in Fig. 4(b). When the temperature is decreased, one observes two distinct features: a narrowing Lorentzian peak around as well as a sharp increase of at .

Only finite times can be reached within the DMRG and it is thus essential to establish a controlled way to assess the associated error of . This can be achieved in various ways. First, one can compare the conductivities obtained via Eqs. (5) and (6), which only coincide if the time integral extends to infinity, as exemplified by solid and dashed lines in Fig. 4(b). As expected, the finite-time error systematically becomes larger at small frequencies. Moreover, the missing contribution from larger times can be estimated by carrying out calculations with bigger discarded weights (which trades accuracy of the real-time data for reaching longer time scales) or by employing extrapolation schemes such as linear prediction (see Refs. barthel, ; scaling, for details). It seems reasonable to combine both strageties to define an error bar as twice as difference between obtained with and without extrapolation or twice the difference between Eqs. (5) and (6) (with extrapolation), whatever is larger. An example is shown in the right inset to Fig. 4(b). Finally, it is instructive is to verify the sum rule

(10) |

where the kinetic energy can be obtained easily from infinite-system DMRG. The left inset to Fig. 4(b) illustrates that Eq. (10) is fulfilled with great accuracy.

Hubbard model — We now turn to the Hubbard model. In order to rule out any subtleties due to a small but potentially finite Drude weight, we switch on a small next-nearest neighbor interaction , which is far from the phase transition occurring at .extendedhub () Fig. 5(a) shows for (where the spectral gap is )hubbook () that the charge current correlation functions decay to zero, and that their integral can be obtained for all temperatures from to . At smaller , the gap becomes exponentially small, and at low one can no longer reach time scales large enough to compute without using extrapolation schemes.

The optical charge conductivity of the almost-integrable Hubbard model is shown in Fig. 5(b). As the temperature is decreased below the gap, a sharpening Lorentzian peak develops around – the DC conductivity grows monotonously. Moreover, increases sharply for frequencies above the optical absorption threshold (which is twice as large as the charge gaphubbook ()) and features a long tail that decays on a scale set by the width of the Hubbard bands ( in our units). One successively approaches the zero temperature form of calculated in Ref. hubT0f, (see also Refs. hubT0b, ; hubT0d, ; hubT0e, ; hubT0g, ; hubT0h, ; hubT0i, ; hubT0j, ; hubT0k, ; hubT0l, for more results on at ; finite but small temperatures were considered in Ref. sigmahub, ; an exact diagonalization study of small systems at can be found in Ref. edhub, ; a lower bound on the infinite-temperature diffusion constant was established in Ref. prosenhub, ).

## V Summary and Outlook

We studied finite-temperature linear response transport properties of the one-dimensional fermionic Hubbard model at half filling. Using real time DMRG numerics, we showed that the charge Drude weight is either zero or small; we established upper bounds. The optical charge conductivity was investigated in presence of small next-nearest neighbor interactions. Its DC part is finite and increases even as the temperature is lowered below the gap (the same holds for the integrable XXZ chain in its gapped phase). Our analysis thus suggests that the half-filled, gapped Hubbard model is a normal charge conductor at finite temperatures.

Aside from finding limits on the Drude weight and estimating transport properties of this long-standing model for comparison to experiment, we believe that the present work can serve as a starting point for considering the effects of additional perturbations, including those that break integrability. Our approach using time-dependent correlation functions can also be extended to compute nonequilibrium physics beyond linear response, as previously done for the XXZ model. Such nonequilibrium calculations could be compared to optical pump-probe experiments on correlated materials. Finally, generalizing the calculations here to energy currents would allow calculations of the thermopower and the electronic contribution to thermal conductivity, which are particularly interesting in the context of proposals to use conducting polymers as low-cost, flexible thermoelectric materials.

Acknowledgments — We are grateful to Fabian Heidrich-Meisner for useful comments. We acknowledge support by the Nanostructured Thermoelectrics program of LBNL (CK), by the Forschergruppe 723 of the DFG (DMK), and by AFOSR MURI as well as the Simons Foundation (JEM).

## References

- (1) A. Rosch and N. Andrei, Phys. Rev. Lett. 85, 1092 (2000).
- (2) J. Sirker, R. G. Pereira, and I. Affleck, Phys. Rev. Lett. 103, 216602 (2009).
- (3) T. Prosen, Phys. Rev. Lett. 106, 217206 (2011).
- (4) F. H. L. Essler, H. Frahm, F. Göhmann, A. Klümper, and V. E. Korepin, The One-Dimensional Hubbard Model (Cambridge University Press, Cambridge, 2005).
- (5) B. S. Shastry and B. Sutherland, Phys. Rev. Lett. 65, 243 (1990)
- (6) X. Zotos, Phys. Rev. Lett. 82, 1764 (1999).
- (7) J. Benz, T. Fukui, A. Klümper, and C. Scheeren, J. Phys. Soc. Jpn. Suppl. 74, 181 (2005).
- (8) B. N. Narozhny, A. J. Millis, and N. Andrei, Phys. Rev. B 58, R2921 (1998).
- (9) J. V. Alvarez and C. Gros, Phys. Rev. Lett. 88, 077203 (2002).
- (10) S. Fujimoto and N. Kawakami, Phys. Rev. Lett. 90, 197202 (2003).
- (11) F. Heidrich-Meisner, A. Honecker, D. C. Cabra, and W. Brenig, Phys. Rev. B 68, 134436 (2003).
- (12) D. Heidarian and S. Sorella, Phys. Rev. B 75, 241104(R) (2007).
- (13) P. Jung and A. Rosch, Phys. Rev. B 76, 245108 (2007).
- (14) J. Herbrych, P. Prelovšek, and X. Zotos, Phys. Rev. B 84, 155125 (2011).
- (15) T. Prosen and E. Ilievski, Phys. Rev. Lett. 111, 057203 (2013).
- (16) C. Karrasch, J. H. Bardarson, and J. E. Moore, Phys. Rev. Lett. 108, 227206 (2012).
- (17) C. Karrasch, J. Hauschild, S. Langer, and F. Heidrich-Meisner, Phys. Rev. B 87, 245128 (2013).
- (18) S. R. White, Phys. Rev. Lett. 69, 2863 (1992).
- (19) U. Schollwöck, Ann. Phys. 326, 96 (2011).
- (20) U. Schollwöck, Rev. Mod. Phys. 77, 259 (2005).
- (21) H. Castella, X. Zotos, and P. Prelovšek, Phys. Rev. Lett. 74, 972 (1995).
- (22) X. Zotos, F. Naef, and P. Prelovšek, Phys. Rev. B 55, 11029 (1997).
- (23) S. Glocke, A. Klümper, and J. Sirker, Phys. Rev. B 76, 155121 (2007).
- (24) X. Zotos and P. Prelovšek, Phys. Rev. B 53, 983 (1996).
- (25) T. Prosen and M. Žnidarič, J. Stat. Mech.: Theory Exp. (2009) P02035.
- (26) R. Steinigeweg and J. Gemmer, Phys. Rev. B 80, 184402 (2009).
- (27) R. Steinigeweg and W. Brenig, Phys. Rev. Lett. 107, 250602 (2011).
- (28) M. Žnidarič, Phys. Rev. Lett. 106, 220601 (2011).
- (29) S. Jesenko and M. Žnidarič, Phys. Rev. B 84, 174438 (2011).
- (30) R. Steinigeweg, J. Herbrych, P. Prelovšek, and M. Mierzejewski, Phys. Rev. B 85, 214409 (2012).
- (31) C. Karrasch, J. E. Moore, and F. Heidrich-Meisner, Phys. Rev. B 89, 075139 (2014).
- (32) R. Steinigeweg, J. Gemmer, and W. Brenig, Phys. Rev. Lett. 112, 120601 (2014).
- (33) G. Vidal, Phys. Rev. Lett. 93, 040502 (2004).
- (34) S. R. White and A. E. Feiguin, Phys. Rev. Lett. 93, 076401 (2004).
- (35) A. Daley, C. Kollath, U. Schollwöck, and G. Vidal, J. Stat. Mech. (2004) P04005.
- (36) P. Schmitteckert, Phys. Rev. B 70, 121302(R) (2004).
- (37) G. Vidal, Phys. Rev. Lett. 98, 070201 (2007).
- (38) M. C. Bañuls, M. B. Hastings, F. Verstraete, and J. I. Cirac, Phys. Rev. Lett. 102, 240603 (2009).
- (39) M. Fannes, B. Nachtergaele, and R. F. Werner, J. Phys. A: Math. Gen. 24, L185 (1991).
- (40) S. Östlund and S. Rommer, Phys. Rev. Lett. 75, 3537 (1995).
- (41) F. Verstraete, J. I. Cirac, and V. Murg, Adv. Phys. 57, 143 (2008).
- (42) F. Verstraete and J. I. Cirac, Phys. Rev. B 73, 094423 (2006).
- (43) A. E . Feiguin and S. White, Phys. Rev. B 72, 220401(R) (2005).
- (44) T. Barthel, U. Schollwöck, and S. R. White, Phys. Rev. B 79, 245101 (2009).
- (45) F. Verstraete, J. J. García-Ripoll, and J. I. Cirac, Phys. Rev. Lett. 93, 207204 (2004).
- (46) M. Zwolak and G. Vidal, Phys. Rev. Lett. 93, 207205 (2004).
- (47) J. Sirker and A. Klümper, Phys. Rev. B 71, 241101(R) (2005).
- (48) S. R. White, Phys. Rev. Lett. 102, 190601 (2009).
- (49) T. Barthel, U. Schollwöck, and S. Sachdev, arXiv:1212.3570.
- (50) E. H. Lieb and F. Y. Wu, Phys. Rev. Lett. 20, 1445 (1968).
- (51) B. S. Shastry, Phys. Rev. Lett. 56, 1529 (1986).
- (52) C. A. Stafford, A. J. Millis, and B. S. Shastry, Phys. Rev. B 43, 13660 (1990).
- (53) R. M. Fye, M. J. Martins, D. J. Scalapino, J. Wagner, and W. Hanke, Phys. Rev. B 45, 7311 (1992).
- (54) C. A. Stafford and A. J. Millis, Phys. Rev. B 48, 1409 (1993).
- (55) F. B. Gallagher and S. Mazumdar, Phys. Rev. B 56, 15025 (1997).
- (56) J. M. P. Carmelo, N. M. R. Peres, and P. D. Sacramento, Phys. Rev. Lett. 84, 4673 (2000).
- (57) E. Jeckelmann, F. Gebhard, and F. H. L. Essler, Phys. Rev. Lett. 85, 3910 (2000).
- (58) Y. Mizuno, K. Tsutsui, T. Tohyama, and S. Maekawa, Phys. Rev. B 62, 4769(R) (2000).
- (59) F. H. L. Essler, F. Gebhard, and E. Jeckelmann, Phys. Rev. B 64, 125119 (2001).
- (60) S. S. Kancharla and C. J. Bolech, Phys. Rev. B 64, 085119 (2001).
- (61) D. N. Aristov, V. Cheianov, and A. Luther, Phys. Rev. B 66, 073105 (2002).
- (62) E. Jeckelmann, Phys. Rev. B 67, 075106 (2003).
- (63) H. Matsueda, T. Tohyama, and S. Maekawa, Phys. Rev. B 70, 033102 (2004).
- (64) N. Maeshima and K. Yonemitsu, J. Phys. Soc. Japan 74, 2671 (2005).
- (65) T. Shirakawa and E. Jeckelmann, Phys. Rev. B 79, 195121 (2009).
- (66) S. Fujimoto and N. Kawakami, J. Phys. A: Math. Gen. 31 (1998) 465.
- (67) S. Kirchner, H. G. Evertz, and W. Hanke, Phys. Rev. B 59, 1825 (1999).
- (68) S-J. Gu, N. M. Peres, and J. M. P. Carmelo, J. Phys.: Condens. Matter 19, 506203 (2007).
- (69) P. Prelovšek, S. El Shawish, X. Zotos, and M. Long, Phys. Rev. B 70, 205129 (2004).
- (70) N. M. R. Peres, R. G. Dias, P. D. Sacramento, and J. M. P. Carmelo, Phys. Rev. B 61, 5169 (2000).
- (71) J. M. P. Carmelo, S-J. Gu, and P. D. Sacramento, Ann. Phys. 339 (2013), 484.
- (72) T. Prosen and M. Žnidarič, Phys. Rev. B 86, 125118 (2012).
- (73) J. M. P. Carmelo, T. Prosen, and D. K. Campbell, arXiv:1407.0732.
- (74) S. Sachdev and K. Damle, Phys. Rev. Lett. 78, 943 (1997).
- (75) K. Damle and S. Sachdev, Phys. Rev. B 57, 8307 (1998).
- (76) Y. Huang, C. Karrasch, and J. E. Moore, Phys. Rev. B 88, 115126 (2013).
- (77) S. Sota and T. Tohyama, Phys. Rev. B 82, 195130 (2010).