Structure Factors of the Kagome-Lattice Heisenberg antiferromagnets at finite temperatures
We compute the real-space spin correlations and frquency and wave-vector resolved dynamic structure factors for the nearest-neighbor Kagome-Lattice Heisenberg Antiferromagnet (KLHM) at finite temperatures using Numerical Linked Cluster Expansion (NLCE) method. A triangle-based 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 fluctuation-dissipation relation is used to to reconstruct the frequency dependence. We find that some features of the low temperature KLHM structure-factors 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 Brillouin-Zone 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.
The spin-half nearest-neighbor Kagome Lattice Heisenberg Model (KLHM) is one of the best studied models of quantum magnetism misguich (); balents (). Recent computational studies have established a quantum spin-liquid ground state for the model, although the full nature of the quantum spin-liquid phase and the existence of a spin-gap remains under debate dmrg (); ran07 (); gapless (); recent (). While the breakthrough DMRG studies suggested a gapped spin-liquid with a robust spin-gap of order or larger than a tenth of the exchange constant dmrg (), several recent studies suggest a Dirac spin-liquid 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 imai-lee (); klhm-e (). The lack of structural distortions and magnetic isolation of copper based Kagome planes by intervening non-magnetic zinc planes, makes these materials well suited to the exploration of the rich quantum spin-liquid physics in these systems. Recent neutron-scattering neutron12 () and NMR Fuscience () measurements on large single-crystal materials find no evidence for magnetic order down to temperatures many orders of magnitude below the exchange energy scale. While the as obtained neutron-scattering spectra clearly shows gapless excitations, a recent analysis of the spectra taking into account the anti-site 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 real-space 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 nmr-theory (). The static structure factors, we calculate, should be very accurate down to the lowest temperatures studied. The prominent features of the wave-vector 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 short-time dynamics gelfand ().
We find that, in the static structure factors as well as the low-energy structure factors, the intensity is mostly spread near the extended Brillouin-Zone boundary once the temperature is below the exchange energy sclae . However, we find that the intensity peaks at the K-point 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:
where the sum runs over all nearest-neighbor bonds of the Kagome lattice. The () represent spin-half operators associated with the spin at site .
The Kagome Lattice consists of 3-sublattices, which we can label by . A site on the Kagome lattice has location:
where and are the lattice translational vectors of the underlying triangular Bravais-lattice, 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 cross-section 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:
Fourier transforming in time gives
where is the time Fourier transform of , which by translational symmetry only depends on vector distance given by , . The neutron scattering cross-section is the sum over all matrix elements of the matrix.
The equal-time correlation function is obtained by summing over all frequencies, and leads to the expression
In this work, we calculate real-space spin-spin 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
Note that given an intensive property such as spin-spin 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 linked-cluster 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 sub-clusters. It is defined as
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. nlc-et, 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 non-zero. The range of real-space 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 real-space 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 Kagome-lattice. 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.
Coordinate-dependent embeddings of the graphs in the lattice are needed for calculating distance dependent spin-spin 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 vector-distance the pair belongs to in the embedding, is determined.
Using an exact diagonalization program, spin-spin correlations are determined for every pair of spins of a topological graph.
Using the assignment of vector-distances to each pair, the spin-spin 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 spin-spin correlation functions for the infinite lattice. This process gives us a sum over all symmetry related spin-spin correlations, per lattice-site. Knowing the number of equivalent vectors in each case, the spin-spin 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 on-site spin operators and the Hamiltonian. This ensures that linked-cluster 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,
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 wave-vector space . To obtain the frequency dependence for any , we use the Gaussian approximation. For this, we first introduce the spectral density defined by
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 fluctuation-dissipation relations.
We begin with the correlations in real space. There are 29 relevant vector distances, and Table-1 gives representative vectors for the first four of them. Note in particular that the vector distances with labels 0 and 1 correspond to on-site and nearest neighbor vector distances respectively.
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 real-space correlations between spins on the Kagome lattice in Figure 3. In these plots we illustrate the equal-time, spin-spin 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 zeroth-moment yields the full q-dependence for the equal-time 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 drop-off 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 on-site 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 high-frequency 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 nmr-theory (). 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 equal-time 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
Thus, they are easy to obtain from the real-space 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 powder-average 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 real-space 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 fluctuation-dissipation relations.
The development of short-range 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 nmr-theory ().
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, large-S 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 spin-half Heisenberg model does not have long-range order in the ground state at all. Nevertheless, there may be competition for short-range 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 spin-liquid to an algebraic spin-liquid 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 wave-vector resolved measurements to be done at higher temperatures.
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 DMR-1306048.
- (1) G. Misguich and C. Lhuillier, arXiv:cond-mat/0310405, Frustrated Spin Systems Edited by H. T. Diep World Scientific, Singapore (2005).
- (2) L. Balents, Nature 464, 199-208 (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); J-W. Mei et al, Phys. Rev. B 95, 235107 (2017); Y-C. 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) T-H. 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, 289-293 (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, cond-mat/arXiv:0901.1065 (unpublished).
- (26) X. Chen, S. Ran, T. Liu, C. Peng, Y. Huang, and G. Su, cond-mat/arXiv:1711.01001.