# Spin hydrodynamics in the S = 1/2 anisotropic Heisenberg chain

## Abstract

We study the finite-temperature dynamical spin susceptibility of the one-dimensional (generalized) anisotropic Heisenberg model within the hydrodynamic regime of small wave vectors and frequencies. Numerical results are analyzed using the memory function formalism with the central quantity being the spin-current decay rate . It is shown that in a generic nonintegrable model the decay rate is finite in the hydrodynamic limit, consistent with normal spin diffusion modes. On the other hand, in the gapless integrable model within the XY regime of anisotropy the behavior is anomalous with vanishing , in agreement with dissipationless uniform transport. Furthermore, in the integrable system the finite-temperature dynamical conductivity reveals besides the dissipationless component a regular part with vanishing .

###### pacs:

05.60.Gg, 71.27.+a, 75.10.Pq## I Introduction

The anisotropic Heisenberg model (AHM) on a one-dimensional (1D) chain is one of the prominent models describing the physics of strongly interacting fermions on a lattice. The model is also well realized in several novel materials, e.g., in the 1D Mott-Hubbard insulator the relevant low-energy degrees of freedom remain spin excitations, and the closest realizations correspond to the isotropic case . The integrability of the model via the Bethe Ansatz has so far led to a number of exact results regarding the spin dynamics at ,(1) as well as some thermodynamic results at . (2) On the other hand, it has been also realized that just the integrability itself, of this particular model or more generally, can be the origin of the anomalous behavior of transport quantities and low-energy dynamics. (3); (4); (5); (6)

A specific consequence of model integrability on a chain of sites is the existence of a macroscopic number of conserved local quantities and related operators commuting with the Hamiltonian and with each other . For the 1D AHM a nontrivial example is representing the energy current and leading directly to its non-decaying behavior(7) and within the linear response to infinite thermal conductivity at any , being also the explanation for the very large spin contribution to heat conductivity in spin-chain materials. (8)

We concentrate in this paper rather on spin dynamics and response, as evidenced by the dynamical spin susceptibility , with the emphasis on the hydrodynamic regime of small . The behavior is also expected to be anomalous since it has been shown that spin conductivity has a ballistic (dissipationless) contribution – spin stiffness – again at any within the gapless XY regime . (6) Although the corresponding spin current is not a conserved quantity, the ballistic component has been well established via several connections: (a) the spin current is closely related to many-body (MB) level dynamics with respect to flux (9) and their independent-like character in an integrable MB model(3); (10) in contrast to level repulsion in a generic nonintegrable system, (b) via the limited decay of correlation functions due the overlap with conserved quantities (7), , at least for a partly polarized system with nonzero magnetization, , (c) the overlap with a steady state solution(11) which – unlike other local conservation laws – is not orthogonal to the spin current in the sector. While the behavior is evidently different(12) but as well highly nontrivial in the gapped Ising-like regime where results seem to favor , (13); (10); (11) there seems to be most controversy on the strict vicinity of the isotropic point . (6); (14); (10) There is much less known about the spin dynamics extended to finite wave vectors as well as frequencies , which has been considered so far in finite systems using exact diagonalization(15); (16) and quantum Monte Carlo(17); (18) (QMC). While for a generic nonintegrable spin system the hydrodynamic regime should reveal diffusive behavior for , in the gapless integrable AHM a singular approach at is expected. In the latter case one of the open questions is the possible coexistence of a ballistic and diffusive processes at . (19); (20) Clearly, the central challenge is the design of a proper phenomenological description of the hydrodynamic response of the gapless integrable as well as weakly perturbed nonintegrable spin system.

In this paper we analyze spin hydrodynamics in the integrable and nonintegrable AHM by the numerical calculation of dynamical spin correlations on finite chains. Numerical data are used as the input to the memory-function representation where the crucial result is the spin-current decay rate at . The latter quantity should be constant within the diffusion regime, as is shown to be for the nonintegrable case. For the integrable gapless AHM within the XY regime, however, our results indicate vanishing as well as . This is well consistent with the dissipationless transport at but puts also restrictions on possible proposed theoretical scenarios.

The paper is organized as follows: In Sec. II we present the model, spin correlations and the memory-function formalism for the dynamical spin susceptibility. Section III is devoted to the presentation of numerical methods. We focus also on the evaluation of the decay rate from dynamical spin correlation functions. Results are presented in Sec. IV and Sec. V. First we investigate transport properties in a nonintegrable case with a next-nearest neighbor interaction. Next we repeat our investigation for high- and low- behavior of dependent in the case of the integrable AHM. We also extend our analysis for the case of finite magnetization . Finally, we focus in Sec. VI on the low- behavior of the uniform dynamical conductivity.

## Ii Memory-function analysis

In this paper we study the anisotropic Heisenberg (XXZ) model on a chain with sites and periodic boundary conditions (PBC),

(1) |

where () are spin operators at site and represents the anisotropy. We allow also for a next-nearest neighbor -interaction with breaking the integrability of the model. It should be reminded that the Hamiltonian (1) can be mapped on a -- model of interacting spinless fermions with hopping and inter-site interactions , . In the fermionic representation, spin transport corresponds to charge transport. We further on use as well as as the unit of energy.

Our aim is to investigate the dynamical spin susceptibility given by

(2) |

with

(3) |

where due to PBC and denotes the thermodynamic average at temperature . Note that the dynamical spin structure factor is related to the imaginary part of Eq. (2) by .

While the results for are obtained numerically for finite systems, we use for further analysis the memory-function (MF) formalism(21); (22); (23) with the central quantity being the relaxation function

(4) |

where is the static spin susceptibility, is the Liouville operator, , and within the MF formalism the scalar product between operators at and corresponding is defined by

(5) |

Applying Mori-Zwanzig reduction(21); (22) we get

(6) |

Here, is the projected dynamical spin conductivity,

(7) |

where the spin current is defined by the conservation law (written at small momentum ) and the orthogonal projection reads . The spin diffusion constant (if it exists) is equal to .

On the second level we can represent

(8) |

where , such that

(9) |

The memory function , and in particular its imaginary part representing the spin-current decay rate, is the central quantity of our analysis. In the generic case of normal transport it should be well behaved in the hydrodynamical regime with a finite value . This is clearly not the case for the dissipationless (integrable) case where one can expect .

On the other hand, and are both finite in the limit . Moreover, both are even independent at high (), as can be easily shown via the high- expansion,

(10) |

where is the magnetization and () is number of up (down) spins.

It should be pointed out that, due to the projection in Eq. (7), is not equal to the “standard” spin conductivity ,

(11) |

Note that the latter quantity is directly related to the relaxation function and thus differs essentially from for , in particular in the regime . On the other hand, for strictly (in a finite system with PBC) both quantities should be equal, , but in this case one has to calculate numerically and corresponding separately [not from ] since is a conserved quantity.

## Iii Numerical methods

In this paper we present results obtained using different numerical methods:

(a) For high studies we use the Microcanonical Lanczos Method(26); (24) (MCLM). The choice is motivated by the fact that in the thermodynamic limit the microcanonical ensemble should yield the same results as the canonical one. For a large enough system and temperature , dynamical autocorrelations can be evaluated with respect to a single wave function characterized by the energy uncertainty

(12) |

where the parameter determines the temperature for which is a relevant representative. Such can be generated via a first Lanczos procedure using instead of . The dynamical correlation functions are then calculated using the standard Lanczos method, where the modified wave function is the starting point for the second Lanczos iteration. Reachable finite-size systems, for magnetization and for , have very high density of states for high , hence statistical fluctuations are effectively smoothed out. This is in contrast to low- properties, dominated by a small number of low-lying MB states, therefore a different numerical method should be applied.

(b) For low we choose the Finite-Temperature Lanczos Method(25); (26) (FTLM) which is reliable for where the finite size temperature is typically in the AHM at available chain lengths . To reduce statistical fluctuations at low , one can introduce additional sampling over initial random vectors and increase the number of Lanczos steps. Due to the memory requirement and the CPU time, reachable system sizes are for and for .

Within the framework of linear response theory we numerically calculate the imaginary part of the dynamical susceptibility and for the spin dynamical conductivity . The MF can be calculated directly from Eq. (9),

(13) |

However, the denominator of this equation is highly influenced by the Kramers-Kronig transformation, e.g., by a small imaginary process added in the numerical realization of the Kramers-Kronig relation,

(14) |

To resolve this problem it is better to perform the Hilbert transform on or and use one of the following expressions for the MF,

(15) |

which comes from simple evaluation of Eq. (8) and Eq. (9) respectively. Finally, we can calculate using one of the sum rules

(16) |

where the integration goes over the whole, numerically obtained MB spectrum, i.e., . It is also worth mentioning that and are weakly dependent in the hydrodynamic limit, as evident for , cf. Eq. (10).

As a final remark of this section we show in Fig. 1 result for the decay rate , evaluated for different system sizes . We choose here and case, but this behavior is generic for “conducting” regime. It is evident that limit of decay rate revel only systematic size () dependence, without sizeable statistical error or deviations.

## Iv Nonintegrable case

First, let us address the question of spin transport in the generic “normal” case, i.e., nonintegrable case as introduced by . In Fig. 2 we present characteristic numerical results within the XY regime at high , as obtained via the MCLM method for a chain of sites and different allowed smallest .

Shown are the dynamical relaxation function , and the extracted projected dynamical conductivity , Eq. (8), as well as the decay rate , Eq. (15). Results obtained here for a substantial integrability breaking term can be easily interpreted with normal diffusion behavior. We note that the projected conductivity reveals a nearly independent Lorentzian form, which is then reflected in , being effectively constant in both and . In particular, it is evident that the result obtained directly from is essentially the same as extracted from . Deviations for are understandable since for our restricted chain lengths the latter ’s are actually not very small.

As a final numerical result of this section we present in Fig. 3 the finite -scaling of the d.c. rate for three different values of the anisotropy , , in the case of non-integrability of the model (1). As already evident from Fig. 2, the variation with is modest and is finite in the limit . Furthermore, the value agrees very well with a perturbation theory at on the basis of Ref. (27), on which we comment in more detail for the remainder of this section.

To this end, it is convenient to turn to the time domain for the moment and to consider an integro-differential equation of the form

(17) |

describing the decay of the dynamical conductivity in time and involving a time-convolution with a memory kernel . Such an equation, and particularly the memory kernel, can be obtained by an application of the Nakajima-Zwanzig (NZ) projection operator technique. (28) For the application of NZ an unperturbed Hamiltonian has to be chosen, where a natural choice is given by the XX model, i.e., . Then the interaction may play the role of a small perturbation and NZ allows to obtain the memory kernel as a lowest-order truncation of a systematic series expansion (in even powers of the perturbation strength involving and ). The somewhat lengthy and subtle calculation of the lowest-order truncation has already been undertaken and is given in Eq. (30) of Ref. (27) by . While the final expression for still needs to be evaluated numerically, the evaluation can be done for several thousands sites, e.g., as chosen here. Thus, numerically integrating the above Eq. (17) by the use of for eventually leads to a lowest-order prediction in the thermodynamic limit for , and after Fourier transforming, for as well. The obvious agreement in Fig. 2 also includes small deviations from a strict Lorentzian, which is a non-Markovian effect since the memory kernel cannot be considered as a mere function on the characteristic scale set by the decay time of . Notably, directly yields a lowest-order prediction on the frequency-dependent decay rate as well, via the relation

(18) |

with an as convincing agreement in Figs. 2 and 3. It is further worth mentioning that the success of NZ relies not only on the limit of small perturbations but also crucially on the non-integrability, allowing to assume vanishing higher-order contributions (equivalent to neglecting the coupling to other observables in the NZ equation). In the case of integrability, or in the limit of strong perturbations, the incorporation of higher order contributions is indispensable and also other variants of projection operator techniques may become convenient (27), not containing a memory kernel at all.

## V Integrable model

In the following we analyze the spin dynamics within the integrable AHM (XXZ) model as realized for . We concentrate here on the gapless (“conducting”) regime, i.e. on the XY regime at zero magnetization , but extending also to for . Note that a different behavior can emerge in the “insulating” (spin-gap) phase at and (even at ).

### v.1 XX model

In the particular case the AHM (i.e. the XX model) can be mapped on the tight-binding 1D model of noninteracting spinless fermions. In this case the calculation of the dynamical susceptibility, Eq. (2), reduced to the well known Lindhard formula and for the relaxation function can be evaluated analytically for arbitrary ,

(19) |

Using the above expression and relation (13) one can express directly the decay rate , and in particular

(20) |

Since the prefactor at any is a constant at , we note a characteristic linear variation and its (nonanalytic) vanishing for . As we see in the following, the XX model can serve as the guideline and explanation for the observed behavior for more general .

### v.2 XXZ model

Next let us focus on the case of a more general AHM with (XXZ model) which corresponds, via the Jordan-Wigner transformation, to a model of interacting spinless fermions and analytical results are not available for and . We are therefore restricted to the numerical calculation of . As described in Sec. III, two approaches are used: (a) the MCLM for high enough and (b) the FTLM for lower .

In Fig. 4 we first present the high- results for and . Again we show first , then the extracted , and finally . The difference to a nonintegrable case is quite evident. Both as well as are strongly dependent, in particular at low , most pronounced for . Here, exhibits a dissipationless component reflected in the vanishing .

To stress the correspondence of the results to the XX case we present in Figs. 4 and 5 results for and in a direct comparison with the analytical expressions Eq. (19) with the same parameters, same , as well as and , respectively. Similarities and differences to the XX case are quite visible: (a) singularities and cutoffs in at persist at as maxima diminishing with increasing , (b) the qualitative behavior of is quite similar for the and cases for lower , (c) for decay rate vanishes above the threshold, i.e., , for on the other hand becomes independent but nonzero in a broad range.

Next we focus on the finite -scaling of . In Fig. 6 we present for the zero-magnetization as a function of for various values of the anisotropy () and for two temperatures and . As we mention before, the and results are calculated from different (numerically obtained) quantities, dynamical conductivity and spin susceptibility , respectively. It is indicative that within the XY regime is linear function of . Also our results show vanishing for . This type of behavior is well consistent with dissipationless (ballistic) transport for zero wave vector . It is, however, also visible that the isotropic case does not follow the simple behavior (similar deviations can be observed also for and low which could be also due to finite-size or FTLM numerical uncertainty). As discussed later the conclusion could be that remains finite (29) or even more that the scaling is different (30).

We also repeat the same analysis for the case of finite magnetization . Figure 7 depicts the -scaling for . Since we are away from the singular isotropic point and for all , (7) here the case does not show any deviations from a general rule.

Finally let us focus on the temperature dependence of . As mentioned, the FTLM is reliable for while presented data for are only estimates. Figure 8 depicts results for , characteristic results for the XY regime. Several conclusions can be drawn directly from the obtained results: (a) For the decay rate is effectively independent, (b) In the regime but where FTLM is reliable, it appears that for all but the -dependence reveals the vanishing decay rate for as expected for low where the model becomes the one of free quasiparticles within the Luttinger 1D liquid.

In the inset of Fig. 8 we show a comparison with the bosonization result in Ref. (29) for and ,

(21) |

where the running coupling constant is determined by the equation

(22) |

We point out that our result, for the smallest possible at , is smaller in the entire regime. It should be also noted is clearly an upper bound to the desired .

## Vi Uniform dynamical spin conductivity

It is by now well accepted that the transport in an integrable 1D model can be dissipationless at any in the “conducting” regime. (7); (6) For the uniform spin current, given explicitly as

(23) |

one can express the conductivity for [] at as

(24) |

where the regular part can be expressed in terms of eigenstates and eigenenergies ,

(25) |

while the dissipationless component with the Drude weight (spin stiffness) is related to matrix elements between degenerated states,

(26) |

where are corresponding Boltzmann factors and is the partition function.

While there are several analytical as well as numerical results supporting in the XY regime of the integrable AHM, (4); (13); (6); (10); (11) there is still controversy on the behavior of the regular part in the same regime. The upper part of Fig. 9 shows two possible scenarios for differing essentially in the behavior at . Note that both scenarios differ also in the integrated normalized spectra

(27) |

also presented in Fig. 9, which are much more reliable (monotonically increasing function) when numerically dealing with finite-system results. Some applications of the MF formalism combined with a coupling to conserved quantities and perturbative scattering of nonconserved quantities (similar to the NZ approach in this paper) indicate on . (29); (31) We will in the following present both an analytical argument as well as numerical results leading to the conclusion that this is not the case.

As the basis of the reasoning we follow Ref. (32) for the integrable 1D model, characterized by a macroscopic number of local conserved quantities . Let us consider a perturbed Hamiltonian by a fictitious magnetic flux, (1); (9) which modifies the exchange (hopping) term in Eq. (1) by the Peierls phase factor . It can be directly verified that the perturbed Hamiltonian is still characterized by the same conserved quantities, i.e., . In particular, we employ here only representing within the integrable AHM the conserved energy current at . The Taylor expansion in of and leads to an exact relation,

(28) |

where and . Evaluating the matrix element of the above equation between eigenstates and , we find

(29) |

where are eigenvalues of and . It seems to be plausible (32); (4) that do not have the same crossing points as , such that the denominator in Eq. (29) remains finite at . Taking this into account as well as Eq. (25) we see that as and finally

(30) |

This scenario is depicted in Fig. 9. This behavior of at low frequency is also clearly visible in the normalized integrated dynamical conductivity , see lower part of Fig. 9. Such a behavior has been already observed in Ref (4); (15); (10); (33).

To strengthen our arguments we numerically investigate and the normalized integrated , obtained with exact (full) diagonalization where -peaks are binned in windows . In Fig. 10 we present characteristic results on and for anisotropies and different magnetizations . As clearly visible, the regular part starts with zero, consistent with Eq. (30). However, due to the small system sizes reachable by ED, it is hard to differentiate the predicted dependence in Eq. (30) from some more general power law of with .

Let us comment on the above results and the form of the memory function , as deduced from the relation Eq. (8). If one would insert a general form

(31) |

into Eq. (8) and assume a nonsingular behavior of , it is evident that the dissipationless component and are only consistent for and . In particular, for reproducing the form (30) we need . Such a result clearly puts strong restrictions to proper analytical approaches to the dynamical response of integrable systems such as the AHM in the “conducting” regime. Still some caution is welcome in the interpretation of our result since both analytical argument as well as numerical results should be firm against finite-size scaling. In particular, the relation (29) requires the uncorrelated eigenvalues of and beyond the finite size (some our preliminary results confirm this conjecture) while certain numerical results also reveal low-frequency finite-size anomalies (10) although not for presented examples.

## Vii Conclusions

Our results confirmed the essential difference of the finite- hydrodynamic behavior between integrable and generic nonintegrable quantum many-body systems. For the example of the dynamical spin susceptibility and spin relaxation function in the 1D generalized AHM we showed that a term breaking integrability, in our case non-zero next-neighbor repulsion , induces a “normal” diffusive behavior in the hydrodynamics, i.e. small regime. This is well reflected in the memory-function analysis where the spin-current decay rate as the central quantity is effectively constant, i.e. in a range of small . Results are in addition well captured within the perturbation-theory approach starting from the XX noninteracting-fermion model.

The integrable model with we investigated in the “conducting” (gapless) regime, i.e. for , and , in this way avoiding the anomalous diffusion as established in the “Mott-insulating” state at .(34) From the existence of dissipationless spin transport and finite stiffness it is evident that the current decay rate should vanish in the limit , i.e. . A nontrivial result of our analysis is, however, that in this respect the universality of the noninteracting XX model is followed throughout the whole gapless regime revealing at . While the similarity to the XX model occurs in the low- regime, this is not the case at larger . In the latter regime for the XX model while induces but only weakly -dependent. In fact, results in this regime again resemble the perturbation-treatment analysis, which at least within the lowest order clearly fails to capture the low-frequency dynamics.

Nontrivial are also conclusions for the uniform case where the results for are extracted from the uniform conductivity . The only consistent possibility with the dissipationless component and vanishing regular part is with , e.g. for the analytical argument is correct. This again indicates that in an integrable system the current scattering is ineffective at low at any in spite of even strong fermion interactions at .

Although we analyzed only spin-density hydrodynamics, one can easily speculate on the behavior of other transport properties, in particular the energy density in the context of heat diffusion. Since is a conserved quantity, the dynamical thermal conductivity is particularly simple with only a dissipationless part and the corresponding thermal-current memory function at any . The extension to finite- energy-density response is not straightforward, but from the analogy to spin hydrodynamics and the noninteracting case one can firmly predict that the decay rate should vanish again as . In a similar way one can possibly speculate on the hydrodynamics of other 1D integrable models such as the 1D Hubbard model.

Experimentally, the most relevant model is the isotropic 1D Heisenberg model realized in several novel materials (8). But in the absence of external magnetic field case it is also the most controversial one. Our results in Fig. 6 both for high as well as for lower are not conclusive due to the restricted system sizes. Nevertheless, together with the assumption (anomalous diffusion) one could expect the dependence with . Since the effective (momentum-dependent) diffusion coefficient , one could conclude that the scaling with the system size should follow . Indeed at high temperatures such a scaling and divergence of has been recently found in nonequilibrium bath scenarios with . (30)

###### Acknowledgements.

This research was supported by the RTN-LOTHERM project and the Slovenian Agency grant No. P1-0044.### References

- B. S. Shastry and B. Sutherland, Phys. Rev. Lett. 65, 243 (1990).
- M. Takahashi and M. Suzuki, Prog. Theor. Phys. 48, 2187 (1972).
- H. Castella, X. Zotos, and P. Prelovšek, Phys. Rev. Lett. 74, 972 (1995).
- X. Zotos and P. Prelovšek, Phys. Rev. B 53, 983 (1996).
- for a review, see X. Zotos and P. Prelovšek, Strong Interactions in Low Dimensions, eds. D. Baeriswyl and L. Degiorgi (Kluwer Academic Publishers), p. 347-382 (2004).
- for a review, see F. Heidrich-Meisner, A. Honecker, and W. Brenig, Eur. Phys. J. Special Topics 151, 135 (2007).
- X. Zotos, F. Naef, and P. Prelovšek, Phys. Rev. B 55, 11029 (1997).
- C. Hess, Eur. Phys. J. Special Topics 151, 73 (2007).
- W. Kohn, Phys. Rev. A 133, 171 (1964).
- J. Herbrych, P. Prelovšek, X. Zotos, Phys. Rev. B 84, 155125 (2011).
- T. Prosen, Phys. Rev. Lett. 106, 217206 (2011).
- A. J. A. James, W. D. Goetze, and F. H. L. Essler, Phys. Rev. B 79, 214408 (2009).
- X. Zotos, Phys. Rev. Lett. 82, 1764 (1999).
- C. Karrasch, J. H. Bardarson, and J. E. Moore, Phys. Rev. Lett. 108, 227206 (2012).
- S. Mukerjee, and B. S. Shastry, Phys. Rev. B 77, 245131 (2008).
- F. Naef and X. Zotos, J. Phys: Condens. Matter 10, L183 (1998).
- S. Grossjohann and W. Brenig, Phys. Rev. B 81, 012404 (2010).
- O. A. Starykh, A. W. Sandvik, and R. R. P. Singh, Phys. Rev. B 55, 14953 (1997).
- M. Takigawa, N. Motoyama, H. Eisaki, and S. Uchida , Phys. Rev. Lett. 76, 4612 (1996).
- K. R. Thurber, A. W. Hunt, T. Imai, and F. C. Chou, Phys. Rev. Lett. 87, 247202 (2001).
- H. Mori, Prog. Theor. Phys. 33, 423 (1965).
- D. Forster, Hydrodynamic Fluctuations, Broken Symmetry, And Correlation Functions, edited by D. Pines (Westview Press, New York, 1995), p. 95-120.
- P. Jung and A. Rosch, Phys. Rev. B 75, 245104 (2007).
- M. W. Long, P. Prelovšek, S. El Shawish, J. Karadamoglou, and X. Zotos, Phys. Rev. B 68, 235106 (2003).
- J. Jaklič and P. Prelovšek, Phys. Rev. B 49, 5065 (1994).
- for a recent review see P. Prelovšek, and J. Bonča, arXiv:1111.5931 (2011).
- R. Steinigeweg, Phys. Rev. E 84, 011136 (2011).
- H.-P. Breuer and F. Petruccione, The Theory of Open Quantum Systems (Oxford University Press, New York, 2007).
- J. Sirker, R. G. Pereira, and I. Affleck, Phys. Rev. Lett. 103, 216602 (2009); J. Sirker, R. G. Pereira, and I. Affleck, Phys. Rev. B 83, 035115 (2011).
- M. Žnidarič, Phys. Rev. Lett. 106, 220601 (2011).
- A. Rosch and N. Andrei, Phys. Rev. Lett. 85, 1092 (2000).
- E. A. Yuzbashyan, B. L. Altshuler, and B. S. Shastry, J. Phys. A: Math. Gen. 35, 7525 (2002).
- M. Rigol and B. S. Shastry, Phys. Rev. B 77, 161101(R) (2008).
- R. Steinigeweg, J. Herbrych, P. Prelovšek, and M. Mierzejewski, Phys. Rev. B 85, 214409 (2012).