Structure Factors of the KagomeLattice Heisenberg antiferromagnets at finite temperatures
Abstract
We compute the realspace spin correlations and frquency and wavevector resolved dynamic structure factors for the nearestneighbor KagomeLattice Heisenberg Antiferromagnet (KLHM) at finite temperatures using Numerical Linked Cluster Expansion (NLCE) method. A trianglebased NLCE is used to calculate frequency moments of the dynamic structure factors in the thermodynamic limit, which show excellent convergence for . A Gaussian approximation and the fluctuationdissipation relation is used to to reconstruct the frequency dependence. We find that some features of the low temperature KLHM structurefactors begin to set in at temperatures of order . Our results are in very good agreement with powder diffraction measurements reported earlier on the Herbertsmithite materials ZnCu(OH)Cl . However, the calculated properties differ from the low temperature () experimental measurements in one important regard. In line with the experimental observations, the spectral weight has a diffuse nature, which is predominantly spread along the extended BrillouinZone boundary. However, the maximum intensity is found in our calculations to be at the point of the extended Brillouin Zone in contrast to the low temperature experiments, where it is at the point. We suggest that experiments should be done at various temperatures to look for such a crossover of the maximum from the point to the point. In the absence of such a crossover, the Herbertsmithite materials must differ from KLHM in a significant manner.
I Introduction
The spinhalf nearestneighbor Kagome Lattice Heisenberg Model (KLHM) is one of the best studied models of quantum magnetism misguich (); balents (). Recent computational studies have established a quantum spinliquid ground state for the model, although the full nature of the quantum spinliquid phase and the existence of a spingap remains under debate dmrg (); ran07 (); gapless (); recent (). While the breakthrough DMRG studies suggested a gapped spinliquid with a robust spingap of order or larger than a tenth of the exchange constant dmrg (), several recent studies suggest a Dirac spinliquid with gapless excitations recent ().
On the experimental front, the Herbertsmithite materials ZnCu(OH)Cl have been celebrated as possibly a nearly ideal realization of a Kagome antiferromagnet imailee (); klhme (). The lack of structural distortions and magnetic isolation of copper based Kagome planes by intervening nonmagnetic zinc planes, makes these materials well suited to the exploration of the rich quantum spinliquid physics in these systems. Recent neutronscattering neutron12 () and NMR Fuscience () measurements on large singlecrystal materials find no evidence for magnetic order down to temperatures many orders of magnitude below the exchange energy scale. While the as obtained neutronscattering spectra clearly shows gapless excitations, a recent analysis of the spectra taking into account the antisite copper impurities in the zinc planes, assigns the low energy scattering entirely to these impurities han16 (). The authors conclude a gap of order for the KLHM, in agreement with the eariler NMR study and somewhat below the DMRG calculations of Yan et al dmrg ().
In this work, we aim to calculate the realspace spin correlations and strcture factors at intermediate temperatures, where the gap issue is not relevant. Our main goal is to benchmark the structure factor for the KLHM, in a temperature regime where they can be accurately calculated in the thermodynamic limit, and thus be directly compared to the experiments. Despite the many theoretical studies, this quantity has not been calculated before apart from a high temperature expansion study at selected wavevectors elstner (), and is important for addressing the question of how good the Heisenberg Model is for these materials nmrtheory (). The static structure factors, we calculate, should be very accurate down to the lowest temperatures studied. The prominent features of the wavevector dependence of the structure factor begin to develop at relatively high temperatures of order . The frequency dependence is obtained through the Gaussian approximation, which should be a good approximation for the shorttime dynamics gelfand ().
We find that, in the static structure factors as well as the lowenergy structure factors, the intensity is mostly spread near the extended BrillouinZone boundary once the temperature is below the exchange energy sclae . However, we find that the intensity peaks at the Kpoint in the extended Brillouin Zone. This is in contrast to the low temperature experimental observation at , where the peak is at the M point neutron12 (); sachdev (); hao (). We note that the finite size calculation of Shimokawa and Kawamura kawamura () also found a crossover of the maximum in the structure factor from the K point to the M point at . Thus, our results are fully consistent with their studies. Our work suggests that measuring the temperature dependence of the structure factors as a function of temperature can help clarify how good the KLHM is for these materials and determine an important crossover energy scale.
Ii Model and Methods
We consider the Heisenberg model with Hamiltonian:
(1) 
where the sum runs over all nearestneighbor bonds of the Kagome lattice. The () represent spinhalf operators associated with the spin at site .
The Kagome Lattice consists of 3sublattices, which we can label by . A site on the Kagome lattice has location:
(2) 
where and are the lattice translational vectors of the underlying triangular Bravaislattice, and are integers, and for are the basis vectors in a unit cell. One explicit representation (taking nearest neighbor distance of unity) is:
with basis vectors,
Neutron scattering measures the scattering crosssection resolved by momentum transfer and energy transfer (we set ). Let us begin with the correlations in time instead of energy transfer . For momentum transfer , the dynamic structure factor is a matrix:
(3)  
Fourier transforming in time gives
(4) 
where is the time Fourier transform of , which by translational symmetry only depends on vector distance given by , . The neutron scattering crosssection is the sum over all matrix elements of the matrix.
The equaltime correlation function is obtained by summing over all frequencies, and leads to the expression
(5)  
In this work, we calculate realspace spinspin correlation functions , as well as static and dynamic structure factors and using the Numerical Linked Cluster method.
Iii Numerical Linked cluster method
The essence of the Numerical Linked Cluster (NLC) method is to express an extenive property for a large lattice with sites as
(6) 
Note that given an intensive property such as spinspin correlation function, one can always construct an extensive property by defining . Here, the sum over runs over all distinct linked clusters of the lattice . is called the lattice constant of the cluster c, and is the number of embeddings of the linkedcluster in the lattice per site. The quantity is called the weight of the cluster and is determined entirely by a calculation of the property on the finite cluster and all its subclusters. It is defined as
(7) 
where the sum over is over all proper subclusters of the cluster . In a high temperature expansion, the property is expanded in powers of inverse temperature. In NLC, one carries out a calculation at a given temperature by exact diagonalization of the finite system.
For the Kagome lattice, it is useful to consider clusters made up of complete triangles only. It was found in Ref. nlcet, that whereas an NLC based on bond or site based graphs starts to breakdown as soon as the high temperature expansion diverges, the triangle based NLC converges down to much lower temperatures. We define the order of the calculation by the clusters with largest number of triangles included in the calculation. We carry out complete calculations for up to triangles or th order.
We first calculate the spin correlations in real space and then Fourier transform to get the correlations in momentum space. The order of the NLC calculation limits the largest vector distance for which the correlations can be nonzero. The range of realspace correlations studied in th order NLC is depicted in Fig. 1. The figure depicts all sites that are within triangles of a given site. We will find that these distances are large enough that, at least at temperatures of interest in this work, the correlations become very small well before the largest distances accessed in the study.
In order to calculate the realspace correlations for all , and all relevant vector distance , we group them into symmetry distinct sets. The correlations will be identical for two vector distances that are related by a symmetry of the Kagomelattice. Up to th order, there are 29 distinct vector distances. The following steps are needed to carry out the calculations:

First we prepare a list of topological graphs, their subgraphs and their lattice constants up to some order . An th order graph, with triangles will have sites and bonds. The topology of the graph is fully specified by the connectivity or adjacency matrix of the graph and this information is sufficient for calculating distance independent properties.

Coordinatedependent embeddings of the graphs in the lattice are needed for calculating distance dependent spinspin correlation functions. For each topological graph, all possibe lattice embeddings are determined up to symmetries of the lattice, together with their symmetry related count.

A list is prepared of all relevant vector distances divided into 29 distinct sets (for order ). That is, every vector between pair of spins in all graphs must be in one and only one of the vector distances in the set.

For all embeddings, an identification for every pair of sites with one of the 29 elements is made. That is, the vectordistance the pair belongs to in the embedding, is determined.

Using an exact diagonalization program, spinspin correlations are determined for every pair of spins of a topological graph.

Using the assignment of vectordistances to each pair, the spinspin correlation sum (and frequency moments) for all the distinct vector distances are calculated for each graph. These define the extensive properties for which weights can now be obtained by subgraph subtraction.

Once weights have been determined, summing over all topological graphs gives us the spinspin correlation functions for the infinite lattice. This process gives us a sum over all symmetry related spinspin correlations, per latticesite. Knowing the number of equivalent vectors in each case, the spinspin correlation between pairs of spins follows.

Fourier Transforming the results gives us the wavevector dependence.
The frequency moments of , can be written as thermal expectation values of commutation relations of onsite spin operators and the Hamiltonian. This ensures that linkedcluster expansion exists. In our NLC calculations, we do not use the commutation relations. We have the exact eigenstates of the graph. Then, for spins at site and , the frequency moments can be calculated from the expression,
(8) 
Where are the eigenvectors and eigenvalues of the Hamiltonian, is inverse temperature, and is the partition function. Due to the spin rotational symmetry of the Heisenberg model, it suffices to calculate the zz correlation functions. We note that the zeroth moment is the equal time correlation function. All results are Fourier transformed to wavevector space . To obtain the frequency dependence for any , we use the Gaussian approximation. For this, we first introduce the spectral density defined by
(9) 
which is an even function in , and also shares all of its even moments with . Since this function is even, we assume that it is a gaussian with zero mean. Thus it is determined from its zeroth and second moment. After NLC and fourier transformation is performed for the zeroth and second moments, we construct in the gaussian approximation and use equation (9) to determine . The benefit of going through the spectral density is that this function is even in , and so a gaussian with mean zero preserves the fluctuationdissipation relations.
Iv Results
We begin with the correlations in real space. There are 29 relevant vector distances, and Table1 gives representative vectors for the first four of them. Note in particular that the vector distances with labels 0 and 1 correspond to onsite and nearest neighbor vector distances respectively.
Label  

0  0  0 
1  1  1 
2  0  2 
3  4  0 
The zeroth and second moment for the first four vector distances are shown in Figure 2. Here we show the result for 5th, 6th and 7th order. It is clear that convergence is excellent, with hardly much difference between 6th and 7th orders down to a temperature of . Thus, for the remainder of this paper we only show results from calculations for our highest order that is 7th order (that is up to 7 triangle graphs), and do not show temperatures below where differences between th and th order begin to arise.
Next we show pictorially the realspace correlations between spins on the Kagome lattice in Figure 3. In these plots we illustrate the equaltime, spinspin correlation between a given site with the site labelled by a green star. On each site we draw a colored circle, whose area quantitatively illustrates the magnitude of the correlation. The color of the circle specifies the sign of the correlation, with red signifying a negative correlation and blue positive. We see that even down to a temperature of , significant correlations do not extend beyond a few lattice constants.
Since these correlations die off so rapidly, it means terminating the Fourier transform over space at finite distances, as we have from only considering up to 7 triangles, is a very good approximation. The Fourier transform of the zerothmoment yields the full qdependence for the equaltime correlations as shown in Figure 4. We see the characteristic development of dark nearly circular patches at the centers of the extended Brillouin Zone (BZ) at relatively high temperatures of order , and the intensity starts to concentrate on the boundaries of the extended BZ. Focusing on the bright zone boundary regions, for all temperatures in this study, the maximum spectral weight is found at the K point, with a decrease in magnitude as we move towards the point, and a rapid dropoff as we move from the K point to the origin .
We calculate the dynamic structure factors using the Gaussian approximation. To assess the limitations of this approximation, we show in Figure 5 a comparison between higher moments obtained from NLC with those obtained by the gaussian approximation. We show this comparison for , as well as for the and points, as defined in Figure 4. We find that for the onsite calculation the gaussian approximation reproduces the higher moments very well and gives a good approximation over the temperature range. However, for the and points, we find that the deviations from gaussianity changes sign at a temperature below . At high temperature the skew is towards lower frequencies, where as it develops a highfrequency asymmetry at lower temperatures.
For the dynamic structure factor in the gaussian approximation, the intensity accumulates most prominently at the K point, and the line towards the point from K holds the most spectral weight. The frequency dependence for the K and M points for several temperature values is shown in Figure 6. The intensity peaks around . On general grounds, one expects the spectral weights to decrease rapidly above nmrtheory (). That rapid decrease is ensured by the Gaussian approximation. For the same q values, we show the temperature dependence of for several values of in Figure 7. The intensity grows monotonically with decreasing temperature for the temperatures shown, and the overall behavior is very similar to the equaltime correlation functions.
We also compare directly with powder experiments by integrating over all points at equal . Here it is important to perform a 3D powder average as appropraite for the experiments of de Vries et al. Vries () The angular averages reduce to
(10) 
Thus, they are easy to obtain from the realspace correlations. Assuming the lattice spacing between neighboring copper ions in Herbertsmithite is as stated by de Vries et al Vries (), we show vs in Figure 8. We find peaks at and , and a trough near , in very good agreement with the experiments Vries () which found a peak at . We also show plots of powderaverage for various temperatures in Figure 9. We show a range of values of , and find peaks developing at and and this intensity diminishes in all direction in the plane from there.
V Discussions and Conclusions
In this paper, we have calculated the realspace spin correlations and wavevector and frequency dependent structure factors of the Kagome Lattice Heisenberg Model (KLHM) at finite temperatures using the Numerical Linked Cluster (NLC) method. These calculations should be very accurate for the static structure factors, in the thermodynamic limit, for . The frequency dependence is obtained using the Gaussian approximation maintaining the fluctuationdissipation relations.
The development of shortrange antiferromagnetic order sets in at temperatures of order and leads to drak patches at the extended Brillouin Zone centers and enhanced spectral weights along the boundaries of the extended Brillouin Zone. Our results for powder diffraction are in very good agreement with the Neutron spectra on the Herbertsmithite materials obtained earlier by de Vries et al Vries (). Earlier NMR relaxation rates, which are a sum over wavevectors, calculated using the NLC, were also found to be in good agreement with experiments nmrtheory ().
Recently full wavevector resolved neutron spectra were measured on single crystals of Herbertsmithites at low temperatures () neutron12 (). We are not aware of any such measurements at higher temperatures. Comparison of our calculated spectra at much higher temperatures than the experiments show agreement with the broad features where spectral weight begins to get concentrated at the boundaries of the extended Brillouin Zone. However, our results differ from the experiments in one rather striking regard. We find the intensity maximum at the point on the corners of the extended Brillouin Zone. In contrast, the low temperature experiments show a peak at the points at the middle of the boundary of the extended Brillouin Zone. We suggest that these measurements should be done as a function of temperature. In the absence of a crossover of a maximum from the point to the point as a function of temperature kawamura (); hao (), the Herbertsmithite materials must differ from the KLHM in some important regard. However, if such a crossover is found, its temperature will provide an important crossover energy scale for the material.
From a theoretical point of view, a peak in the structure factor at the point is consistent with order in the pattern, which is favored in classical, largeS and many computational approacheshao (); harris (); rutenberg (); chernyshev (); richter (); oitmaa () , where as a peak at the point is consistent with order at . As the computational studies show, these two classical patterns are very close in energy but the spinhalf Heisenberg model does not have longrange order in the ground state at all. Nevertheless, there may be competition for shortrange order reflected in the dominance of different q points sachdev (). The higher temperature behavior favors the semiclassical and perturbative approaches. But, the ground state studies of the largest clusters studied show otherwise kawamura (); lauchli (). A recent theoretical work discusses a crossover from a conventional spinliquid to an algebraic spinliquid at a temperature above xichen (). It would be interesting to have experimental input to the competition and crossover between these ordering tendencies, which would require the wavevector resolved measurements to be done at higher temperatures.
Vi Acknowledgements
We would like to thank Jeff Rau for a careful reading of the manuscript. This work is supported in part by the National Science Foundation grant number DMR1306048.
References
 (1) G. Misguich and C. Lhuillier, arXiv:condmat/0310405, Frustrated Spin Systems Edited by H. T. Diep World Scientific, Singapore (2005).
 (2) L. Balents, Nature 464, 199208 (2010).
 (3) S. Yan, D. A. Huse and S. R. White, Science 332, 1173 (2011); S. Depenbrock, I. P. McCulloch and U. Schollwock, Phys. Rev. Lett. 109, 067201 (2012); H. C. Jiang, Z. H. Wang and L. Balents, Nat. Phys. 8, 902 (2012);
 (4) Y. Ran, M. Hermele, P. A. Lee, and X. G. Wen, Phys. Rev. Lett. 98, 117205 (2007).
 (5) Y. Iqbal, D. Poilblanc and F. Becca, Phys. Rev. B 91, 020402(R) (2015).
 (6) H. J. Liao et al, Phys. Rev. Lett. 118, 137202 (2017); JW. Mei et al, Phys. Rev. B 95, 235107 (2017); YC. He et al, Phys. Rev. X 7, 031020 (2017).
 (7) T. Imai and Y. Lee, Physics Today, 69, 30 (2016).
 (8) J. S. Helton et al, Phys. Rev. Lett. 98, 107204 (2007); P. Mendels et al Phys. Rev. Lett. 98, 077204 (2007); A. Zorko et al Phys. Rev. Lett. 101, 026405 (2008); T. Han, S. Chu, Y. S. Lee, Phys. Rev. Lett. 108, 157202 (2012).
 (9) T. H. Han et al, Nature 492, 406 (2012).
 (10) M. Fu et al., Science 350, 655 (2015)
 (11) TH. Han et al, Phys. Rev. B 94, 060409(R) (2016).
 (12) N. Elstner and A. P. Young Phys. Rev. B 50, 6871 (1994).
 (13) N. E. Sherman, T. Imai, and R. R. P. Singh Phys. Rev. B 94, 140415(R) (2016).
 (14) M. P. Gelfand and R. R. P. Singh, Phys. Rev. B 47, 14413 (1993); R. R. P. Singh and M. P. Gelfand, Phys. Rev. B 42, 996 (1990).
 (15) M. Punk, D. Chowdhury, and S. Sachdev, Nature Physics 10, 289293 (2014).
 (16) Z. Hao and O. Tchernyshyov, Phys. Rev. B 81, 214445 (2010).
 (17) M. Rigol, T. Bryant, R. R. P. Singh, Phys. Rev. E 75, 061118 (2007); ibid. Phys. Rev. Lett. 97, 187202 (2006).
 (18) T. Shimokawa and H. Kawamura, J. Phys. Soc. Jpn. 85, 113702 (2016)
 (19) M. A. de Vries et al, Phys. Rev. Lett. 103, 237201 (2009).
 (20) A. B. Harris, C. Kallin and A. J. Berlinsky, Phys. Rev. B 45, 2899 (1992).
 (21) D. A. Huse and A. D. Rutenberg, Phys. Rev. B 45, 7536 (1992).
 (22) A. L. Chernyshev, M. E. Zhitomirsky Phys. Rev. Lett. 113, 237202 (2014); Phys. Rev. B 92, 144415 (2015).
 (23) O. Goetze and J. Richter, Phys. Rev. B 91, 104402 (2015).
 (24) J. Oitmaa and R. R. P. Singh, Phys Rev B 93, 014424 (2016).
 (25) A. Laeuchli and C. Lhuillier, condmat/arXiv:0901.1065 (unpublished).
 (26) X. Chen, S. Ran, T. Liu, C. Peng, Y. Huang, and G. Su, condmat/arXiv:1711.01001.