Nested Cluster Algorithm for Frustrated Quantum Antiferromagnets
aa r X i v : . [ c ond - m a t . s t r- e l ] M a r Nested Cluster Algorithm for Frustrated Quantum Antiferromagnets
M. Nyfeler, F.-J. Jiang, F. K¨ampfer, and U.-J. Wiese
Institute for Theoretical Physics, Bern University, Sidlerstrasse 5, 3012 Bern, Switzerland
Simulations of frustrated quantum antiferromagnets suffer from a severe sign problem. We solvethe ergodicity problem of the loop-cluster algorithm in a natural way and apply a powerful strategyto address the sign problem. For the spin Heisenberg antiferromagnet on a kagom´e and on afrustrated square lattice, a nested cluster algorithm eliminates the sign problem for large systems.The method is applicable to general lattice geometries but limited to moderate temperatures.
PACS numbers: 75.10Jm, 75.40Mg, 75.50.Ee
Monte Carlo calculations are a powerful tool for firstprinciples investigations of strongly coupled quantumsystems. Early simulations of quantum spin systemswere performed with local Metropolis-type algorithms[1]. They suffered from critical slowing down and werethus limited to rather high temperatures. Cluster al-gorithms perform non-local updates and are capable ofsubstantially reducing critical slowing down. Such al-gorithms were first developed by Swendsen and Wangfor discrete classical Ising and Potts spins [2] and thengeneralized by Wolff [3] to classical spins with a con-tinuous O ( N ) symmetry. Improved estimators whichaverage over an exponentially large number of configu-rations at polynomial cost are an additional benefit ofcluster algorithms. The first cluster algorithm for thespin quantum Heisenberg model was developed in [4].While that algorithm works efficiently only for quantumspin chains, the loop-cluster algorithm [5] is efficient alsoin higher dimensions, and was first applied to the 2-dspin Heisenberg antiferromagnet on a square lattice in[6]. The continuous-time variant of the algorithm elim-inates the Suzuki-Trotter time-discretization error andcan reach very low temperatures [7]. This method hasalso been used to simulate systems on very large lattices[8] and with very long correlation lengths [9]. An elegantand powerful related method based on stochastic seriesexpansion is available as well [10].Unfortunately, in many cases of physical interest,including frustrated quantum spin systems, quantumMonte Carlo calculations suffer from a very severe signproblem. Using an improved estimator, the sign problemof the 2-d classical O (3) model at vacuum angle θ = π hasbeen addressed with a variant of the Wolff cluster algo-rithm [11]. In that case, some clusters are half-instantonsalso known as merons. Flipping a meron-cluster leadsto a sign-change of the Boltzmann weight and hence toan exact cancellation between two configurations. As aconsequence, only configurations without meron-clusterscontribute to the partition function. Restricting the sim-ulation to those configurations eliminates the sign prob-lem, since all configurations in the zero-meron sector havea positive sign. The meron concept has been gener-alized to fermionic systems [12] and the meron-cluster algorithm has been used to solve a number of very se-vere fermion sign problems [13, 14, 15]. Unfortunately,the meron-cluster algorithm is not generally applicable.In fact, as shown in [16], some sign problems are NP-complete. Hence, a hypothetical method that can solveany sign problem would solve all NP-complete problemsin polynomial time. This would imply the equality of thecomplexity classes NP=P. Since it is generally believedthat NP =P, it is expected that a universally applicablemethod that solves all sign problems cannot exist. In thispaper we construct a nested cluster algorithm which, forthe first time, is capable of eliminating severe sign prob-lems of frustrated antiferromagnets at least at moderatetemperatures.Let us consider the antiferromagnetic spin quantumHeisenberg model with the Hamiltonian H = X x,y ∈ Λ J xy ~S x · ~S y . (1)Here ~S x is a quantum spin operator located at the site x of a lattice Λ, and J xy > x and y . Although our method can be applieddirectly in the Euclidean time continuum, in order toease its implementation we work in discrete time. De-pending on the lattice geometry, the Hamiltonian H = H + H + ... + H M is expressed as a sum of M terms H i which leads to a Suzuki-Trotter decomposition of thepartition function Z = Tr exp( − βH )= lim ε → Tr [exp( − εH ) exp( − εH ) ... exp( − εH M )] N . (2)Here the inverse temperature β = 1 /T = N ε representsthe extent of a periodic Euclidean time interval, whichis divided into N discrete time steps of size ε . Each H i is a sum of mutually commuting pair interactions h xy = J xy ~S x · ~S y on a set of disconnected bonds. Insert-ing complete sets of spin states s x,t = ± = ↑ , ↓ betweenthe factors exp( − εH i ) in eq.(2), the partition function isexpressed as a path integral over all spin configurations[ s ] [6] Z = X [ s ] Sign[ s ] exp( − S [ s ]) . (3)The weight of a configuration is a product of con-tributions from individual space-time plaquettes cor-responding to the two-spin transfer matrix elements h s x,t s y,t | exp( − εh xy ) | s x,t +1 s y,t +1 i . In the basis | ↑↑i , | ↑↓i , | ↓↑i , | ↓↓i , up to an irrelevant overall factor, thetwo-spin transfer matrix takes the formexp( − εh xy ) = A A + B − B − B A + B
00 0 0 A , (4)with A = exp( − εJ xy /
2) and B = sinh( εJ xy / s ] = ± − S [ s ]) represents a positiveBoltzmann weight which can be interpreted as a proba-bility and thus can be used for importance-sampling in aMonte Carlo simulation.When one samples the system using the positive weightexp( − S [ s ]), one must include Sign[ s ] in the measuredobservables O [ s ] and expectation values are given by h O i = 1 Z X [ s ] O [ s ] Sign[ s ] exp( − S [ s ]) = h O Sign i + h Sign i + . (5)Here the index + refers to expectation values in the simu-lated ensemble with positive Boltzmann weights and par-tition function Z + = P [ s ] exp( − S [ s ]) such that h Sign i + = 1 Z + X [ s ] Sign[ s ] exp( − S [ s ]) = ZZ + ∼ exp( − ∆ f βV ) . (6)Here V is the spatial volume and ∆ f is the differencebetween the free energy densities of the original ensem-ble with the weight Sign[ s ] exp( − S [ s ]) and the simulatedensemble with the positive weight exp( − S [ s ]). The ex-pectation value of the sign is exponentially small in thespace-time volume βV . Since it is obtained as a MonteCarlo average of contributions Sign[ s ] = ±
1, one needs anexponentially large statistics in order to accurately mea-sure h Sign i + . This is impossible in practice and gives riseto a very severe sign problem.How can one increase the statistics by an exponentialfactor without investing more than a polynomial numeri-cal effort? The meron-cluster algorithm [11, 12] achievesthis by constructing an improved estimator for the sign.Like the meron-cluster algorithm, the method presentedhere is based on the loop-cluster algorithm [5] which deco-rates a spin configuration with bonds connecting spins to closed loop-clusters. The four spins on a space-time pla-quette are connected in pairs. In fact, A and B in eq.(4)represent weights of two possible bond configurations ona space-time plaquette. The weight A corresponds tobonds connecting the spins s x,t and s y,t with their time-like neighbors s x,t +1 and s y,t +1 , while B corresponds tospace-like bonds connecting s x,t with s y,t and s x,t +1 with s y,t +1 . Sites connected by bonds form a closed orientedloop-cluster. Up to an overall spin-flip of the entire clus-ter, the spin configuration on a cluster is determined bythe cluster geometry. Time-like bonds connect parallelspins, while space-like bonds connect anti-parallel spins.Integrating out the spins, the partition function can beexpressed as a sum over bond configurations [ b ] Z = X [ b ] Sign[ b ] A n A B n B N C . (7)Here n A is the number of time-like and n B is the numberof space-like plaquette break-ups, while N C is the num-ber of loop-clusters. The factor 2 N C arises because eachcluster has two possible spin orientations. The partitionfunction can be sampled by a Metropolis update of theplaquette break-ups. Remarkably, while the original clus-ter algorithm which operates on spins and bonds neverchanges the sign and is thus not ergodic [17], the algo-rithm which operates only on bonds (after the spins havebeen integrated out) is ergodic and still avoids unnaturalfreezing. Interestingly, Sign[ s ] remains invariant undercluster flips, i.e. all clusters are non-merons. However, inthis case the meron-cluster algorithm does not solve thesign problem because almost half of the configurationsin the zero-meron sector have a negative sign [17]. Sinceit does not change under spin flips, Sign[ s ] = Sign[ b ] isuniquely determined by the bond configuration. It isimportant to note that the sign can be expressed as aproduct of cluster signs Sign[ b ] = Q C Sign C . Depend-ing on the orientation of a cluster, each space-like break-up contributes a factor ± i to the two clusters traversingthe corresponding space-time plaquette. By construc-tion, each cluster traverses an even number of space-likebreak-ups, and hence Sign C = ± C defines the set of lattice sites Λ C contained in C . Theinner Monte Carlo algorithm generates clusters with dif-ferent orientations that visit all sites of Λ C in differentorders, thus contributing different values of Sign C . Inthis process, break-ups that lead to the decompositionof Λ C into separate clusters must be rejected. The innerMonte Carlo algorithm estimates an average h Sign C i i foreach set of sites Λ C . Since the different sets are indepen-dent, the improved estimator of the sign is given by h Sign i i = Y Λ C h Sign C i i . (8)Remarkably, the nesting of an outer and an inner clusteralgorithm achieves exponential error reduction at poly-nomial cost. A similar strategy was very successfullyapplied to the measurement of exponentially suppressedWilson loops in lattice gauge theory [18] as well as toquantum impurity models [19]. Correlation functions andsusceptibilities can also be measured with improved es-timators. Let us consider the staggered magnetizationoperator ~M s = P x z x ~S x . Here z x is a stagger factor de-pending on the sub-lattice to which the site x belongs.The corresponding staggered susceptibility χ s = h M s Sign i + βV h Sign i + = hh M s Sign i i i + βV hh Sign i i i + . (9)is obtained from an improved estimator which is given interms of M s = P C M s C with M s C = P ( x,t ) ∈C z x s x,t as h M s Sign i i = X Λ C h M s C Sign C i i Y Λ C′ =Λ C h Sign C ′ i i . (10)In which cases will the nested cluster algorithm elim-inate or at least substantially reduce the sign problem?Since some sign problems are NP-hard, it is expected thatany method will fail at least in those cases. The nestedcluster algorithm fails to solve the sign problem when acluster fills almost the entire volume, because then theinner Monte Carlo algorithm becomes inefficient. Sincelarge clusters necessarily arise in the presence of largecorrelation lengths, the nested cluster algorithm does notwork efficiently in low-temperature ordered phases.Even in the absence of long-range order, cluster al-gorithms may become inefficient if the clusters grow tounphysically large sizes beyond the physical correlationlength. This potential problem is prevented when there isa reference configuration that limits cluster growth [15].For the antiferromagnet on the square lattice the refer-ence configuration is given by the classical N´eel state, i.e.all spins in a loop-cluster are in a staggered pattern. Thecluster-size squared is then tied to the staggered suscep-tibility which protects the clusters from growing to un-physically large sizes. Also for frustrated systems it isnatural to consider a classical ground state as a refer-ence configuration. When one quantizes the spins alonga local quantization axis in the direction of the spin ori-entation in the classical ground state, an interesting al-gorithm with open string-clusters emerges. The spins ineach cluster are in the reference configuration and hence FIG. 1: kagom´e lattice (left) and frustrated square (oranisotropic triangular) lattice (right) consisting of three sub-lattices
A, B, C . these clusters are protected from becoming unphysicallylarge. However, the meron-concept does not apply to theopen string-clusters, i.e. when these clusters are flipped,they are not independent but affect each other in theireffect on the sign. Remarkably, one can still integrateout the spins analytically. This glues the open string-clusters together to the closed loop-clusters of the algo-rithm discussed before. While typical closed loop-clustersare hence larger than the correlation length correspond-ing to the classical order, they still represent physicalcorrelated regions. In fact, they grow up to the lengthscale at which the signs, which are a manifestation ofquantum entanglement, decorrelate.Even if the typical cluster-size is moderate, the innerMonte Carlo algorithm may not lead to an efficient can-cellation of signs. For example, there are cases in whichthe improved estimator h Sign i i is not positive. Still, ifsuch cases are rare, the sign problem is substantially re-duced. In order to optimize the performance of the al-gorithm, the numerical effort invested in the inner andouter Monte Carlo procedures must be properly balancedagainst each other. It pays off to invest a larger numberof inner Monte Carlo sweeps on the larger sets Λ C . Inany case, the efficiency of the nested cluster algorithmmust be investigated on a case by case basis.We now consider the Heisenberg antiferromagnet withuniform nearest-neighbor coupling J xy = J on the lat-tices illustrated in figure 1. The frustrated square latticehas an additional coupling J ′ along the diagonals. Wehave simulated large kagom´e lattices with up to V ≈ βJ ≈
1. Figure 2shows the probability distribution of the improved es-timator h Sign i i . Although sometimes it is negative, itstill leads to an accurate determination of the averagesign. We consider M s with z x = 1 , − , A, B, C , respectively, which may signal coplanar spin or-der. As shown in figure 3, with increasing volume V both h Sign i + and h M s Sign i + decrease dramatically over nu-merous orders of magnitude, but are still accurately ac-counted for by the nested cluster algorithm. For example,with V = 882 spins h Sign i + = 2 . × − . A bruteforce approach would require an astronomical statistics of -1.0 × -9 -5.0 × -10 × -10 × -9 < Sign > i β J = 1V = 576 spins
FIG. 2:
Probability distribution of h Sign i i for the kagom´e lat-tice with V = 576 spins and βJ = 1 . V -18 -12 -6
FIG. 3:
Volume-dependence of h Sign i + and h M s Sign i + (rescaled by − ) for the kagom´e lattice with βJ = 1 . about 10 sweeps in order to achieve a similar precision.Figure 4 shows the coplanar staggered susceptibility χ s compared to the collinear N´eel susceptibility χ N . On thesquare lattice, frustration reduces the N´eel order, while(at least for J ′ = J/
4) the coplanar order is as weak ason the kagom´e lattice (and practically indistinguishablefrom it in figure 4).To conclude, in contrast to other Monte Carlo meth-ods, the nested cluster algorithm is capable of eliminatingvery severe sign problems for large systems, at least atmoderate temperatures. This is useful, for example, fordetermining the couplings of frustrated magnets by com-parison with experimental finite temperature data. Aswe have demonstrated, although the nested cluster algo-rithm cannot reach very low temperatures, by studyingappropriate susceptibilities one may still obtain valuableinsights concerning possible types of order. Applicationsto frustrated antiferromagnets on various lattice geome-tries are currently in progress. U.-J. W. likes to thank S. Chandrasekharan for a long-lasting fruitful and very inspiring collaboration on thesign problem. We also have benefited from interestingdiscussions with M. Troyer. This work was supported bythe Schweizerischer Nationalfonds. [1] J. D. Reger and A. P. Young, Phys. Rev. B37 (1988)5978.[2] R. H. Swendsen and J.-S. Wang, Phys. Rev. Lett. 58(1987) 86.[3] U. Wolff, Phys. Rev. Lett. 62 (1989) 361.[4] U.-J. Wiese and H.-P. Ying, Phys. Lett. A168 (1992) 143.[5] H. G. Evertz, G. Lana, and M. Marcu, Phys. Rev. Lett.70 (1993) 875.[6] U.-J. Wiese and H.-P. Ying, Z. Phys. B93 (1994) 147.[7] B. B. Beard and U.-J. Wiese, Phys. Rev. Lett. 77 (1996)5130.[8] J.-K. Kim and M. Troyer, Phys. Rev. Lett. 80 (1998)2705.[9] B. B. Beard, R. J. Birgeneau, M. Greven, and U.-J.Wiese, Phys. Rev. Lett. 80 (1998) 1742.[10] A. W. Sandvik, Phys. Rev. B56 (1997) 11678.[11] W. Bietenholz, A. Pochinsky, and U.-J. Wiese, Phys.Rev. Lett. 75 (1995) 4524.[12] S. Chandrasekharan and U.-J. Wiese, Phys. Rev. Lett.83 (1999) 3116.[13] S. Chandrasekharan, J. Cox, K. Holland, and U.-J.Wiese, Nucl. Phys. B576 (2000) 481.[14] S. Chandrasekharan and J. C. Osborn, Phys. Rev. B66(2002) 045113.[15] S. Chandrasekharan, J. Cox, J. C. Osborn, and U.-J.Wiese, Nucl. Phys. B673 [FS] (2003) 405.[16] M. Troyer and U.-J. Wiese, Phys. Rev. Lett. 94 (2005)170201.[17] P. Henelius and A. W. Sandvik, Phys. Rev. B62 (2000) β V χ N unfrustrated square (J’= 0) χ N frustrated square (J’= J/4) χ S frustrated square (J’= J/4) χ S kagome FIG. 4:
Coplanar staggered susceptibility χ s and collinearN´eel susceptibility χ N as functions of the space-time volume βV for the kagom´e as well as the frustrated ( J ′ = J/ ) andunfrustrated ( J ′ = 0 ) square lattice at fixed space/time aspectratio √ V /βJ = 20 ..