Recent progress on cluster and meron algorithms for strongly correlated systems
NNoname manuscript No. (will be inserted by the editor)
Recent progress on cluster and meron algorithms forstrongly correlated systems
Debasish Banerjee
Received: date / Accepted: date
Abstract
Ab-initio studies of strongly interacting bosonic and fermionic sys-tems is greatly facilitated by efficient Monte Carlo algorithms. This articleemphasizes this requirement, and outlines the ideas behind the constructionof the cluster algorithms and illustrates them via specific examples. Numericalstudies of fermionic systems at finite densities and at real-times are sometimeshindered by the infamous sign problem, which is also discussed. The construc-tion of meron cluster algorithms, which can solve certain sign problems arediscussed. Cluster algorithms which can simulate certain pure Abelian gaugetheories (realized as quantum link models) are also discussed.
Keywords
Cluster Algorithms · Fermionic and Spin models · Lattice GaugeTheories
The ab-initio studies of strongly correlated systems occurring in Nature, whetherin particle physics or in condensed matter physics, is an extremely challengingtopic. Analytical solutions are difficult to find in most cases, and perturba-tion theory fails for strong couplings. Certain weak coupling methods (suchas the epsilon-expansion [1,2]) have been successful in addressing the exis-tence of fixed points in renormalization group flow, as well as compute criticalexponents at phase transitions by systematically improving over the meanfield estimates. Large-N methods have provided another analytic handle onsome interesting quantum field theories (QFTs) [3]. In the case of conformalfield theories (CFTs) there has been exciting developments through the use ofAdS/CFT [4], and the more recent conformal bootstrap [5] and the large charge
Saha Institute of Nuclear Physics, HBNI, 1/AF Bidhannagar, Kolkata 700064, India.E-mail: [email protected] a r X i v : . [ h e p - l a t ] J a n Debasish Banerjee expansion [6]. Tensor Networks [7] have significantly contributed to comput-ing paradigms in lower dimensional systems, by opening up the possibility tosimulate a large range of strongly interacting systems, even in real-time.However, in the overwhelmingly large majority of cases, the Markov ChainMonte Carlo (MCMC) methods, starting from the initial proposal of Metropo-lis et al [8] have provided an unbiased ab-initio method to numerically samplemultidimensional integrals and evaluate the expectation values of physical op-erators. Let us denote the degrees of freedom in a (classical or quantum) systemby { q i } , i = 1 , · · · , V , and the partition function for the system as Z = (cid:90) D q e − S ( { q i } ) = Tr exp( − βH ) , (1)where S ( { q i } ) is an action functional corresponding to the Hamiltonian H atan inverse temperature β . The quantity e − S ( { q i } ) = W ( { q i } ) is the Boltzmannweight and is typically positive for the chosen computational basis of { q i } . Ifthis is not the case, we encounter the sign problem and importance samplingfails. We will assume the positivity of the Boltzmann weight for the moment,and will come back to the exceptions later. The expectation value of a physicaloperator O (for example an order parameter or a correlation function) is (cid:104)O(cid:105) = 1 Z (cid:90) D q O e − S ( { q i } ) (2)Starting from an initial probability distribution, a Markov Chain is a se-ries of probability distributions Π k ( { q i } ), which steadily converge to the fixedpoint distribution Π (cid:63) ( { q i } ), which is the equilibrium distribution W ( { q i } ).We will assume that the reader is familiar with the basics of Markov Chains,and local Monte Carlo algorithms and point to excellent textbooks which coverthis topic in depth [9,10]. While it is not always difficult to construct a MarkovChain with the necessary properties, the difficulty lies ensuring that the rateof convergence is insensitive to V. Further, even if a Markov chain were to con-verge to the equilibrium distribution, one needs to sample this distribution toobtain uncorrelated samples on which expectation values and correlation func-tions of relevant operators can be measured. This is exactly where improvedalgorithms, such as the cluster algorithm, assert their importance. To makethis notion quantitative, we first introduce the concept of autocorrelationtime .Operationally, after equilibrium is reached, (local) Monte Carlo algorithmgenerates a set of configurations of ( { q i } ) according to the Boltzmann weight,on which physical operators are measured. Let us denote the distribution atthe l -th Monte-Carlo time as ( { q i } ) l . Typically the distributions ( { q i } ) l and( { q i } ) l +1 are highly correlated since only a small change is involved in one stepto the next. Therefore, measurement of the physical operator on these twosubsequent configurations, O l ( { q i } ) and O l +1 ( { q i } ) are also not independent.Numerical calculations necessarily deal with finite data, let us denote thisby N conf . The sample mean of the dataset is O = N conf (cid:80) l O l . According to the ecent progress on cluster and meron algorithms for strongly correlated systems 3 central limit theorem, the sample mean O has a Gaussian distribution aboutthe exact expectation value (cid:104)O(cid:105) , such that (cid:104)O(cid:105) = O ± σ O . With uncorrelatedmeasurements, the best unbiased estimate of the error from the finite data-setis σ O = (cid:34) N conf ( N conf − (cid:88) k ( O k − O ) (cid:35) . (3)Note that we have used O as the best estimate of the true mean, (cid:104)O(cid:105) . Forlarge N conf , σ O = (cid:34) N conf (cid:88) m ( O m − O ) 1 N conf (cid:88) n ( O n − O ) (cid:35) = 1 N N conf (cid:88) m,n Γ O ( m − n ) . (4)The autocorrelation function, Γ O ( m ), is defined through the average over theequilibrium distribution Γ O ( m ) = (cid:104)O m (cid:48) O m + m (cid:48) (cid:105) − (cid:104)O(cid:105) = Γ O ( − m ) m (cid:29) ∼ C exp( − m/τ exp ) (5)and decays exponentially at asymptotic times, with the decay time τ exp corre-sponding to the slowest mode in the Monte Carlo dynamics. The integratedautocorrelation time , τ int , defined as2 τ int ≈ M (cid:88) m = − M Γ O ( m ) = Γ O (0) + 2 M (cid:88) m =1 Γ O ( m ) (6)and the cut M is due to finiteness of the dataset, and we have ignored cor-rections of order τ exp /N conf , for large N conf . Note that Γ O (0) is the samplevariance. The expression for the error on the mean is then σ O = 2 τ int Γ O (0) N conf = Γ O (0) N eff (7)implying that there are only N eff independent configurations in the Markovchain. Unlike the equilibrium values, τ int depends on the algorithm. Near acritical point, for example, autocorrelation times diverge with the physicalcorrelation length ξ as τ ∼ cξ z , and z ∼ z is calledthe dynamical critical exponent). Such studies are therefore numerically chal-lenging especially near the critical point, and the phenomena is called criti-cal slowing down . Cluster algorithms can guarantee extremely small τ int , andsometimes even z ∼
0, thereby practically eliminating this problem.
Using the formulation of Fortuin and Kasteyln [11,12], the first cluster al-gorithm for simulating classical spin (the Ising and the Potts) models wereconstructed by Swendsen and Wang [13], and then extended by Niedermayer
Debasish Banerjee [14] and Wolff [15] for continuous spins systems. We refer to the review [16] formore details. Cluster algorithms for continuous systems, such as hard spheres,have also been developed. We do not discuss them here, but refer to [17] for amore detailed exposition. We instead focus on cluster algorithms for quantumsystems, starting with quantum spins.Cluster algorithms for quantum spin systems, in the world line formu-lation of quantum spins, were first developed in [18,19], and extended to thecontinuous time version in [20]. We note that the Stochastic Series Expansion(SSE) of the quantum spin Hamiltonian as developed in [21] leads to a loopMonte Carlo updating method [22], similar to the cluster algorithm. Therehas been significant developments in generalizing this algorithm for a varietyof models, larger quantum spins, and we refer the reader to [23] for detailsand a more complete set of references. We note that the worm algorithm isanother powerful idea which has led to the development of efficient algorithmfor a large class of models, and in many cases is as powerful as the clusteralgorithms. We refer to the interested reader to [35] for details.Let us illustrate the construction of a cluster algorithm for the Heisenbergmodel, which not only has a very rich physical pedagogy behind it, but is alsouseful to understand the magnetic properties in certain electronic systems[25]. This method is independent of spatial dimensions, but we will consider(2 + 1) − d, anti-ferromagnetic version J > L , the partition function Z atan inverse temperature β is: Z = Tr exp( − βH ); H = J (cid:88) x, ˆ i =1 , (cid:126)S x · (cid:126)S x +ˆ i . (8)Using the Suzuki-Trotter formula [26], we separate the Hamiltonian into 4parts (2 d parts in d − spatial dimensions), H = H + H + H + H , such thatthe spins contained in each part mutually commute. H = J (cid:88) x =(2 m,n ) (cid:126)S x · (cid:126)S x +ˆ1 ; H = J (cid:88) x =(2 m +1 ,n ) (cid:126)S x · (cid:126)S x +ˆ1 ; H = J (cid:88) x =( m, n ) (cid:126)S x · (cid:126)S x +ˆ2 ; H = J (cid:88) x =( m, n +1) (cid:126)S x · (cid:126)S x +ˆ2 . (9)We construct Z by inserting intermediate time-slices, such that β = 4 N (cid:15) ,and (cid:15) is the temporal lattice spacing. Using I = | n k (cid:105)(cid:104) n k | , we expand Z asa product over the matrix elements of the transfer matrix, expressed in the S z = ± basis, which is chosen as the computational basis: Z = Tr exp( − βH ) = lim N →∞ (cid:15) → (cid:88) n (cid:104) n | (cid:104) e − (cid:15) ( H + H + H + H ) (cid:105) N | n (cid:105) = lim N →∞ (cid:15) → (cid:88) { n } n N = n N − (cid:89) k =0 (cid:104) n k | e − (cid:15)H | n k +1 (cid:105)(cid:104) n k +1 | e − (cid:15)H | n k +2 (cid:105)(cid:104) n k +2 | e − (cid:15)H | n k +3 (cid:105)(cid:104) n k +3 | e − (cid:15)H | n k +4 (cid:105) (10) ecent progress on cluster and meron algorithms for strongly correlated systems 5 We have used the Trotter formula in the second line, and the matrix elementsconnect only the specific bonds. Finally, we can express the matrix elementsvia the local action Z = (cid:89) x,t (cid:88) s ( x,t )= ± e − S [ s ,s ,s ,s ] (11)where the term S [ s , s , s , s ] connects two neighboring spins, s , s at a time-slice, t , with their forward-in-time partners s , s . Note that only a set ofbonds are active at a given time-slice, and all others are passive, and this4-spin interaction traces an active plaquette . A schematic figure for theTrotter decomposition in (1 + 1) − d is shown in Fig 1 (left), where the shadedplaquettes are the active ones, and indicate the spin pairs which are interactingat that given time. A representative s , s , s , s are also shown in the samefigure. For the anti-ferromagnetic Heisenberg model, the explicit values of thetransfer matrix (in the basis where the states of the two spins at sites x and y are | ↑↑(cid:105) , | ↑↓(cid:105) , | ↓↑(cid:105) , | ↓↓(cid:105) ) are:e − S [ s ,s ,s ,s ] = (cid:104) s s | e − (cid:15)J (cid:126)S x (cid:126)S y | s s (cid:105) = e − (cid:15)J (1 + e (cid:15)J ) (1 − e (cid:15)J ) 00 (1 − e (cid:15)J ) (1 + e (cid:15)J ) 00 0 0 1 (12)Note that there are only two off-diagonal elements, both of which are nega-tive. For a bipartite lattice, the negative sign can be eliminated by doing aunitary transformation on every alternative spin, such that the off-diagonalelements become (e (cid:15)J − sign problem . In the chosen basis the probability weightsare not positive definite, and Monte Carlo simulation cannot be performed.Next, to construct a cluster algorithm we have to expand the configurationspace by including bond variables, together with the quantum spins. Thisimplies choosing the bond variables [ b ] such that the following equation issatisfied for the different possibilities of s , s , s , s ,e − S [ s ,s ,s ,s ] = (cid:88) [ b = A,B ] e − S [ s ,s ,s ,s ; b ] . (13)In this particular example, as shown in Fig 1, the linear equations can be sat-isfied by the use of two bond breakups, A and B . The clusters are constructedas follows: (a) start from an initial site, (b) follow the chosen breakups asexplained in the caption of Fig 1 (right) until the initial point is reached. Inthis case, the resulting cluster is the same as a loop, and hence the relation tothe loop algorithm. Flipping the cluster implies the operation | ↑ k (cid:105) ↔ | ↓ k (cid:105) ,which preserves the Boltzmann weight of the configuration, and yet changesthe configuration globally, thereby decorrelating them very fast. This is thereason for the efficiency of the cluster algorithm — the clusters correspond to Debasish Banerjee xt ↑ ↓ ↑ ↓ ↑↑ ↓ ↑ ↓ ↑↑ ↑ ↓ ↓ ↑↑ ↑ ↓ ↓ ↑↑ ↓ ↑ ↓ ↑ s s s s A B ↑ ↑↑ ↑ ↓ ↓↓ ↓ = 1 = A ↑ ↓↓ ↑ ↓ ↑↑ ↓ = (e (cid:15)J −
1) = B ↑ ↓↑ ↓ ↓ ↑↓ ↑ = (e (cid:15)J + 1) = A + B Fig. 1 ( left ) An example of a world-line configuration of the anti-ferromagnetic Heisenbergmodel. There is one transition and one anti-transition which spatially displaces a spin-up (down) and brings it back. ( right ) Possible bond-types for this model: (a) the A-typebreakup connects identical spins forward in time, (b) the B-type breakup connects opposite(spin-up with spin-down) spins sideways. We can figure out the probability of the weightsby computing which breakup can be applied for which interaction plaquette. For thosewhere the spin-states do not, and cannot change (top), only A can be applied, whereaswhere the spin state flips (middle) only B can be applied. The last case is where the spinstate can change but did not. Here, the A-breakup can be applied with the probability W A = A/ ( A + B ) = 2 / (exp( (cid:15)J ) + 1). correlated degrees of freedom, which get efficiently updated by flips. In fact,representing the configurations in terms of clusters allows one to constructimproved estimators for quantities such as the magnetization, susceptibilityand the correlation function.Two variants of this cluster algorithm are well-known: multi-cluster algo-rithm, which proceeds to construct all the possible clusters on a given config-uration of spins and bonds, and then flips each cluster with a 50% probability,and the single-cluster , or the Wolff algorithm which constructs a single clus-ter at always flips it. It is also possible to construct this cluster algorithm incontinuous time.It must be emphasized that for the cluster algorithm to be successful, cer-tain special configurations called reference configurations need to exist. Inthe case of the Heisenberg anti-ferromagnet the reference configuration(s) arethe ones where the spins are staggered on the bi-partite lattice. Any configu-ration can then be decomposed into a certain set of clusters, and the clustersare accordingly flipped to bring the configuration into (one of) the referenceconfigurations. This trivially proves ergodicity. In many cases (examples are ecent progress on cluster and meron algorithms for strongly correlated systems 7 frustrated magnets, spin glasses, gauge theories), while it is conceivable tobuild clusters, absence of reference configurations causes the clusters to growtoo big, almost filling up the whole volume. In such a case, the cluster algo-rithm does not perform any better than a local update algorithm. Attempting to simulate fermions brings us head on with the fermion signproblem . This is most directly seen in the world-line representation of thefermions in the occupation number basis. In dimensions d ≥
2, it is easy toconstruct world-lines of identical fermions which twist around each other asthey evolve in Euclidean time. Configurations which have fermions exchangingpositions an odd number of times have an overall negative sign compared withthose where fermions exchange positions an even number of times (everythingelse being identical), due to the Pauli exclusion principle. A particularly clearstatement of this fermion sign problem and what it entails to solve the problemis given in [27].As in the case of quantum spins, to simulate fermionic systems we constructthe partition function as usual. In the Lagrangian formulation, one uses Grass-mann variables, which is typically integrated out to leave behind a determinantinvolving auxiliary fields (if the theory has four-fermi interactions, or Yukawacouplings), or even gauge fields (as in quantum electrodynamics, or quantumchromodynamics). Updating procedures either involve updating the determi-nant directly [28,29], or recast the determinant in terms of bosonic fields,which are then updated [30]. To construct cluster algorithms for fermions, wefirst obtain a bosonic representation using the Jordan-Wigner transformation.If the resulting bosonic system can be efficiently updated using a cluster algo-rithm, then one can identify types of interactions for which the fermion signproblem can be solved [31]. When such an approach succeeds, it is known as the meron algorithm first introduced in [32], and will be discussed below. Othernovel approaches, such as the fermion bag [33,34] or the diagrammaticMonte Carlo [35,36] are getting increasingly popular in simulating fermionicsystems. The former method expands the fermionic action and groups regionswhere the fermions can propagate freely against regions where they are boundinto monomers and dimers. The Monte Carlo procedure then updates theseregions, which are denoted as fermion bags. The latter approach directly sam-ples the Feynman diagrams associated with the interactions, and in certainregimes can be related to the bag approach.Before discussing the construction of the meron algorithm, we elaboratethe nature of the fermionic sign problem, which will also serve to understandthe solution. First, we set up the fermionic path integral in the occupationnumber basis { n k } : Z f = Tr exp( − βH f ) = (cid:88) { n k } p ( { n k } ); (cid:104) A (cid:105) = 1 Z f (cid:88) { n k } A ( { n k } ) p ( { n k } ) , (14) Debasish Banerjee where p ( { n k } ) is the probability for a given configuration, and A is an operatorin the occupation number basis. To qualify the severity of this problem, one canconsider the sign of the configuration as a part of the observable A through thedecomposition p ( { n k } ) = sign( { n k } ) | p ( { n k } ) | , the sign being ± (cid:104) A (cid:105) f = (cid:80) { n k } A ( { n k } ) p ( { n k } ) (cid:80) { n k } p ( { n k } )= (cid:80) { n k } A ( { n k } )sign( { n k } ) | p ( { n k } ) | (cid:80) { n k } | p ( { n k } ) | × (cid:80) { n k } | p ( { n k } ) | (cid:80) { n k } sign( { n k } ) | p ( { n k } ) | = (cid:104) A · sign (cid:105) b (cid:104) sign (cid:105) b (15)The subscript b indicates that the averaging is done over a bosonic systemwhose weights, | p ( { n k }| , are positive definite by construction. The quantity (cid:104) sign (cid:105) b measures the severity of the sign problem. It can shown that (cid:104) sign (cid:105) b =exp( − βV ∆f ), where ∆f = f f − f b is the difference in free energy densitybetween the fermion and bosonic ensembles. At low temperatures and largevolumes, this quantity is exponentially small:( σ sign ) b (cid:104) sign (cid:105) b = (cid:2) (cid:104) sign (cid:105) b − (cid:104) sign (cid:105) (cid:3) √ N (cid:104) sign (cid:105) b ≈ exp( βV ∆f ) √ N , (16)where N is the number of uncorrelated data, and we have used (cid:104) sign (cid:105) b =1 , (cid:104) sign (cid:105) b ≈
0. Clearly, N ∼ exp(2 βV ∆f ) to determine the sign. Thus even if (cid:104) A (cid:105) f ∼
1, it is measured as the ratio of two exponentially small signals, andrequires exponential effort to extract from statistical noise.One could aim to cancel the negative signs by pairing the configurationswhich carry an − = sign, sincesign = 0 , σ sign ) b (cid:104) sign (cid:105) b = (cid:2) (cid:104) sign (cid:105) b − (cid:104) sign (cid:105) (cid:3) √ N (cid:104) sign (cid:105) b ≈ exp( βV ∆f / √ N , (17)Note that we have achieved an exponential gain in statistics, but one stillneeds N ∼ exp( βV ∆f ) for any meaningful result. The meron algorithm in-corporates the two steps together: it identifies the configurations which can beanalytically cancelled, and never generates them. Thus, one always generatesconfigurations which contribute non-trivially to Z f . The algorithm, however,only works for a restricted class of interactions. ecent progress on cluster and meron algorithms for strongly correlated systems 9 Let us illustrate this for the t − V model in (2+1) − d, extensively used in thestudy of certain electronic systems (though the algorithm can be constructedin any dimensions). The Hamiltonian isH f = − t2 (cid:88) x, ˆ i =1 , ( c † x c x +ˆ i + c † x +ˆ i c x ) + V (cid:88) x, ˆ i =1 , ( n x −
12 )( n x +ˆ i −
12 ) , (18)with V ≥ t >
0. The fermion creation (annhiliation) operators at site x, c x ( c † x ) satisfy the following anti-commutation relations: { c x , c y } = { c † x , c † y } =0; { c † x , c y } = δ xy . The corresponding bosonic system is obtained by using theJordan-Wigner transformation. An ordering of the lattice points is defined: l = x + y · L (in d = 2), and the fermionic operators are expressed as follows: c † x = σ · σ · · · σ l − σ + l ; c x = σ · σ · · · σ l − σ − l ; n x = c † x c x = 12 ( σ l +1) . (19)This preserves the anti-commutation relations. We use the (fermion) occupa-tion number basis ( n x | (cid:105) = 0 , n x | (cid:105) = | (cid:105) , with 0 denoting an unoccupied and1 an occupied site) to construct Z f by splitting the Hamiltonian into differentparts and invoking the Suzuki-Trotter formula. With some algebra to keeptrack of the negative signs, the transfer matrix between the adjacent time-slices can be expressed through the following 4 × x and y denote | (cid:105) , | (cid:105) , | (cid:105) , | (cid:105) ):e − S [ n ,n ,n ,n ] = (cid:104) n n | e (cid:15) t2 ( c † x c y + c † y c x ) − (cid:15)V ( n x − )( n y − ) | n n (cid:105) = e (cid:15)V e − (cid:15)V/ (cid:15) t / Σ sinh( (cid:15) t /
2) 00 Σ sinh( (cid:15) t /
2) cosh( (cid:15) t /
2) 00 0 0 e − (cid:15)V/ (20)The structure of the transfer matrix is very similar to the one of the quantumspins, the main difference is the string operator Σ = σ l +1 · σ l +2 · · · σ m − be-tween the sites x = l and x + ˆ i = m , where l < m in the lattice ordering. Thisallows the rewriting of Z f in the occupation number basis [ { n x = 0 , } ] as Z f = (cid:88) { n } sign[ { n } ]e − S [ { n } ] . (21)Some example configurations which contribute to Z f are shown in Fig. 2.Without the sign factor, the weights of the resulting bosonic system are iden-tical to that of the XXZ-Hamiltonian: H = (cid:88) x, ˆ i =1 , (cid:104) t( S x · S x +ˆ i + S x · S x +ˆ i ) + V S x · S x +ˆ i (cid:105) , (22)which reduces to the anti-ferromagnet for t = V = J . The Z f is decomposedin terms of bonds, in addition to spins: Z f = (cid:88) [ n,b ] sign[ n, b ] exp([ − S [ n, b ]]) (23) C1 C2C3 C4C5 C1C1 C10 1 2 3 4 xt f f f f f f f f Fig. 2
Example configurations which contributes to Z f . ( left ) The two fermions f and f are static during their entire Euclidean time-evolution, and trace out vertical worldlines.For t = J , only A (vertical) and B (horizontal) bonds are required. The figure shows one bond configuration. The oriented clusters are constructed by starting from a single pointand following the bonds – upwards (downwards) for a filled (empty) site for A − breakups,sideways for a filled (empty) site touching B − breakups. This choice of breakups gives riseto 5 different clusters, as marked in the figure. ( middle ) On flipping clusters marked 1, 3and 5 we obtain a completely different world-line configuration, but one where the fermionsdo not exchange positions. ( right ) Flipping cluster 5 gives rise to a configuration wherethe fermions interchange positions, and hence due to the anticommutation relations, thisconfiguration has an overall negative sign compared to the one on the left. The cluster 5 isa meron . Interestingly, in the limit t = V the only breakups that are needed are the A and the B breakup, with same probabilities as derived in Fig. 1 (right). Fig 2shows a typical fermion configuration, and a particular set of breakups, andhow clusters are identified. Cluster flips are operations { n k ↔ − n k | n k ∈ C i } ,where C i is the i-th cluster. Cluster flips give rise to new worldlines significantlydifferent from their parent ones, but carry the same weights.The crucial consequence of these cluster rules is that the sign of the wholeconfiguration factorizes into a product of the signs associated with each indi-vidual cluster. An example is already seen in Fig 2 (right). The clusters whichcan change the sign of the configuration are called merons . It is possible toreach a reference configuration by flipping clusters appropriately. In thisexample, the configuration in Fig 2 (left) is the reference configuration, whichhas a positive sign (sign C i = 1) for all clusters i. When a cluster is flipped, itscontribution to sign is sign C i = 1( −
1) if the quantity N w + N h / N w is the temporal winding number of the cluster, while N h is the numberof spatial hops [32]. The meron concept allows exact pairing of even and odd ecent progress on cluster and meron algorithms for strongly correlated systems 11 signs: Z f = (cid:88) [ n,b ] sign[ b ] exp( − S [ n, b ]) = (cid:88) [ b ] , zero − meron N C exp( − S [ b ]);sign[ b ] = 12 N C (cid:88) clusterflips sign[ n, b ] . (24)sign[ b ] = 0 if at least one of the clusters is a meron. Further, the existenceof the reference configuration with a positive Boltzmann weight, guaranteesthat the unpaired configurations all carry positive weight, and can be reachedby flipping clusters from an initial configuration which has zero-merons. TheMonte Carlo sampling proceeds by generating configurations which only liein the zero-meron sector. If a proposed flip gives rise to meron cluster, it isrejected. Thus, one is able to completely solve the sign problem with the meronconcept.Let us emphasize that for generic interactions the meron concept does nothold, and whether a cluster is a meron depends on the orientation of otherclusters. Therefore, much more (exponential) computational effort is neededto identify the merons. However, there is a large class of interactions whichcan be solved with the meron concept, and we refer the reader to [31] for areview. Cluster algorithms for gauge theories run into problems since the relevantvariables that need to be updated get frustrated, and very often the clustersoccupy the entire volume. The resulting algorithm is not much better thana local update algorithm. One direction of research tried to identify the rel-evant degrees of freedom, which may or may not be identical to the degreesof freedom in which the action is constructed. For example, for a φ theory,the φ field can be separated into an Ising ( Z ) variable and a modulus vari-able. The Z degrees of freedom were updated using a cluster algorithm whilethe modulus of the φ field was updated with a local update algorithm [37],which resulted in an efficient algorithm with a low ( z <
1) dynamical expo-nent. Similar methods have been applied to the SU (2) lattice gauge theory[38,39] targeting the Polyakov loop as the embedded degree of freedom. Whilesome success on the theories with one and two time-slices was reported, thesemethods do not work well for larger lattices.A completely different formulation of gauge theories, the quantum linkmodel formulation [40,41,42], which realizes continuous gauge symmetrieswith finite dimensional Hilbert spaces and generalizes the Wilson construc-tion of lattice gauge theories, has recently been more amenable to simulationwith cluster algorithms. As an illustrative example, we will describe the U (1)quantum link model (QLM). The degrees of freedom of this Abelian latticegauge theory are quantum links, U x, ˆ i and electric fluxes, E x, ˆ i , both of which U (cid:3) U x +ˆ i, ˆ j U x +ˆ j, ˆ i U x, ˆ j U x, ˆ i U (cid:3) = U x, ˆ i U x +ˆ i, ˆ j U † x +ˆ j, ˆ i U † x, ˆ j H J = − JH λ = λH = 0 Fig. 3 ( left ) Layout of the plaquette, and the four spin interaction. ( right ) The action ofthe operators in the Hamiltonian on the possible plaquette states. The U (cid:3) ( U † (cid:3) ) term flips aclockwise (anti-clockwise) oriented plaquette to an anti-clockwise (clockwise) oriented one,and annhiliates the remaining 14 possible states. This is the J − term of the Hamiltonian,and is an off-diagonal operator. The λ − term counts the total number of flippable (both inthe clockwise and anti-clockwise orientation) plaquettes. are operators defined on the bonds connecting neighboring sites, x and x + ˆ i and ˆ i = 1 , d = 2. These operators satisfy canonical commutation relations:[ E x, ˆ i , U y, ˆ j ] = U x, ˆ i δ xy δ ij ; [ E x, ˆ i , U † y, ˆ j ] = − U † x, ˆ i δ xy δ ij ; [ U x, ˆ i , U † y, ˆ j ] = 2 E x, ˆ i δ xy δ ij , (25)The operators E, U can be chosen to be the components of a spin- S object: U = S + , U † = S − , E = S . The parameter S can be thought of regulatingthe local Hilbert space dimension. In the limit S → ∞ [46], these modelsreproduce the Hamiltonian formulation of Wilson gauge theories [47]. TheHamiltonian we will be interested in only contains operators which commutewith the Gauss’ Law: H = − J (cid:88) (cid:3) (cid:16) U (cid:3) + U † (cid:3) (cid:17) + λ (cid:88) (cid:3) (cid:16) U (cid:3) + U † (cid:3) (cid:17) ; G x = (cid:88) i ( E x,x + i − E x − i,x )(26)In addition, we could have included any function of the electric fluxes, and inparticular, the kinetic energy of the gauge fields (cid:80) x,i E x,i . When the quantumlink operators are considered in the spin- representation, this model hasseveral interesting physical applications. In (2+1) − d, the physics of this modelis relevant for the understanding the behavior of low-temperature frustratedmagnets [43], and spin-liquid phases [44,45]. A close cousin of this model isthe quantum dimer model (QDM), especially well-known in condensed matterphysics as a toy model to describe certain aspects of superconductivity [48,49,50]. Therefore, we continue with the gauge links in the spin- representation,and hence ignore the kinetic energy term of the gauge links, which is a constant.The local link Hilbert space is 2-dimensional, and in the electric flux basis they ecent progress on cluster and meron algorithms for strongly correlated systems 13 Fig. 4 ( left ) The Gauss’ Law condition for the QLM, in d = 2, this allows 6 possible stateson each vertex, while ( right ) the Gauss’ Law condition for the QDM in d = 2 allows for 4states for each vertex. The red circle represents Q = +1, while the blue circle the Q = − are denoted by upward (downward) arrows for E = ( − ) for vertical links,and right (left) arrow for E = ( − ) for horizontal links. The action of theHamiltonian on the electric flux states is shown in Fig 3.Specifying the Gauss’ Law further specifies which basis states can con-tribute to the partition function, and which not. For the QLM, the Gauss’Law is chosen as G x | ψ (cid:105) = 0, for every site x , for every eigenstate | ψ (cid:105) of theHamiltonian. The QDM chooses a different set of states, which are specifiedby G x | ψ (cid:105) = ( − x + x | ψ (cid:105) . Physically, this implies that the vacuum selected bythe QLM is charge neutral at each vertex, while the QDM chooses a vacuumwith staggered positive and negative unit charges. The link states allowed atthe vertex for the QLM is shown in Fig 4 (left) and for the QDM in Fig 4(right). The Gauss’ Law generates the gauge transformations, which can begenerically denoted as V = (cid:81) x exp( iθ x G x ), where θ x ∈ (0 , π ] is the parameterof the transformation. The Hamiltonian is invariant under this transformation: (cid:101) H = V † · H · V = H , since [ G x , H ] = 0, for all x. Physically, G x labels dif-ferent super-selection sectors of the theory which do not mix under unitarydynamics.The partition function of this model is defined by projecting the eigenstatesof the Hamiltonian into a specified Gauss’ Law sector, Z = Tr[exp( − βH ) P G ].The projection operator P G = (cid:81) x G x , therefore constrains the Hilbert spacefurther. To construct the cluster algorithm, we note that we can divide theHamiltonian into two parts H = H A + H B , such that the all terms in eachpart mutually commute. Thus, the Trotterization divides the lattice into evenand odd sub-lattices, such that at a given time-step only a single sub-lattice ++ −− − ++ − −− ++ + −− + Fig. 5
Mapping of an electric flux configuration (shown with arrows on the links) to aheight configuration (shown with + and − variables at the centers of the plaquettes). The+ and − are placeholders for 0 , ± ) on A (B) sub-lattices. Every time a flux pointingright or upwards (corresponding to E = ) is crossed, the height variable is changed, whileit remains unchanged if a left or downward pointing flux (corresponding to E = − ) iscrossed. The configuration shown has all plaquettes flippable, and the corresponding heightvariables are in their reference configuration . needs to be updated. For a fixed β = (cid:15)N , we have a two-step transfer matrix Z = Tr (cid:2) ( T A T B ) N P G (cid:3) . (27)The single plaquette transfer matrix can be expressed as T (cid:3) = 1 + ( U (cid:3) + U † (cid:3) )e − (cid:15)Jλ sinh( (cid:15)J ) + ( U (cid:3) + U † (cid:3) ) (cid:2) e − (cid:15)Jλ cosh( (cid:15)J ) − (cid:3) (28)The resulting transfer matrix is 16-dimensional (i.e. a 16 ×
16 matrix), and itis easy to read off the Boltzmann weights from this equation. For plaquetteswhich are not flippable, only the diagonal element is unity, and all other off-diagonal elements are zero. For the two flippable plaquette configurations,the diagonal contribution has the weight e − (cid:15)λ cosh( (cid:15)J ) while the off-diagonalelements (which indicate the plaquette flips) carry the weight e − (cid:15)λ sinh( (cid:15)J ).The idea behind the construction of the cluster algorithm is to dualize thegauge theory in (2 + 1) − d, which gives rise to a Z (2) quantum height modelin (2 + 1) − d. The dualization is an exact rewriting of the partition functionin terms of new degrees of freedom which are Z (2) degrees of freedom locatedat dual sites. As shown in Fig 5, every flux configuration can be mapped to aheight configuration. A configuration of quantum height variables is assignedthe values h A ˜ x = 0 , h B ˜ x = ± . located at the dual sites ˜ x = ( x + , x + ),and is associated with a flux configuration E x, ˆ i = (cid:104) h X ˜ x − h X (cid:48) ˜ x +ˆ i − ˆ1 − ˆ2 (cid:105) mod2 = ±
12 ;
X, X (cid:48) ∈ {
A, B } . (29) ecent progress on cluster and meron algorithms for strongly correlated systems 15 m m m m m m m , m , m , m height variableslives at timeslice t m ( m ) lives at timeslice t-1 (t+1)Bonds for top-right panel A join m and m B keep m and m separateBonds for bottom-right panel A join m , m , m , m B keep m , m , m , m separate 011 0+ + e − (cid:15)Jλ cosh( (cid:15)J ) A + B
011 0+ − e − (cid:15)Jλ sinh( (cid:15)J ) B
011 0+ + e − (cid:15)Jλ cosh( (cid:15)J ) A + B
010 1+ + B Fig. 6 ( left top ) Layout of the height variables across the timeslices t-1, t and t+1. ( leftbottom ) Two kinds of bonds are enough to solve the linear equations: the A-breakup joinsheights, while the B-breakup keeps the heights separate. Different scenarios need to be con-sidered — out-of-plane (temporal) bonds, and in-plane (spatial) bonds. ( right top ) Clusterrules for temporal bonds: with m , m , m , m in a reference configuration and m = m ,both A- and B-breakups can contribute, while if m (cid:54) = m then only the B-breakup con-tributes. To satisfy detailed balance, we must have P A = A/ ( A + B ) = e − (cid:15)Jλ / cosh( (cid:15)J )and P B = tanh( (cid:15)J ). In the case when m = m and m , m , m , m are not in a referenceconfiguration, then we have to bind m and m with probability 1, otherwise a forbiddenconfiguration will be generated. ( right bottom ) Spatial breakups: if m (cid:54) = m then connect m , m , m , m with probability 1, to prevent the generation of a forbidden configuration.If m = m then two different cases can occur: if m , m , m , m are in a reference config-uration, then the Boltzmann weight is e − (cid:15)Jλ cosh( (cid:15)J ) and we can apply either the A, or theB-breakup. Otherwise, the weight is 1 and only B breakup is applied ( m and m are notbound together). Solving the equations, we get P A = 1 − e − (cid:15)Jλ / cosh( (cid:15)J ) and P B = 1 − P A . The modulo function acts by adding or subtracting 2 until the result is inthe desired range ( − , Q = ± Z (2) nature of theheight variables.Before deriving the cluster rules, it is useful to understand the Boltzmannweights of the configurations in terms of the height variables, as shown Fig6. In terms of the height variables, the transfer matrix gives rise to a 6-height interaction between the heights m , · · · , m . Of these, m , m , m , m live on a timeslice t, m lives on the timeslice t-1 and m on t+1. Depend-ing on the values of m , m , m , m , the height m can undergo a transi-tion to m , or not (corresponding to a plaquette flip). In this language if( m , m , m , m ) = (1 , , ,
0) or (0 , , ,
1) the plaquette is flippable, andform the reference configurations. The weights of these configurations andthe derivation of the cluster rules is illustrated in Fig 6 and the figure caption.With these cluster rules implemented in terms of the height variables, thesuper-selection sectors which have charges Q x = ± Q x = ±
2, since they are compatiblewith the height representation. Therefore, in addition to the cluster rules forthe Hamiltonian, the ones for the Gauss’ Law need to additionally imple-mented. To implement the Gauss’ Law note that we need to consider fourheight variables meeting around a site: the top-left and the bottom-right be-long to a one sub-lattice (say A), and top-right and bottom left belong tothe other sub-lattice (B). It can be easily verified that if the site contains acharge ±
2, then both the height variables in both sub-lattices are locally outof their reference configuration. Thus, this situation should always be avoided.This can be implemented as follows: when updating the A sub-lattice at agiven site, it is checked if the corresponding sites of the B sub-lattice are outof the reference configuration. If this is the case then we bind the two heightvariables together, so that they are always flipped together, and a forbiddenconfiguration is not generated. If the relevant height variables in the B sub-lattice are in a reference configuration then no additional bonds are put onthe A sub-lattice heights. An identical procedure is followed while updatingthe other sub-lattice. It is sufficient to check this on a single timeslice, sincethe Gauss’ Law commutes with the Hamiltonian.We emphasize that this cluster algorithm for the gauge theory is somewhatdifferent from the usual type of cluster algorithms since clusters are separatelybuilt on each sub-lattice depending on the values of frozen height variablesin the other sub-lattice. In the next step, the frozen sub-lattice is updated.This algorithm was first implemented in [51] to study the physics of the U(1)QLM and it uncovered new phases of confined gauge theories which breaktranslation and charge conjugation symmetry [52].It is interesting to note that, even though this algorithm is very efficient onthe QLM, extending this to the QDM is not as useful. Since the Hamiltonian isthe same, the same cluster rules can be applied. The problem is however withthe different Gauss’ Law, which now frustrates the reference configuration. Dueto the absence of the reference configuration in the chosen basis of the quantumheight model, the clusters grow quite large to occupy about 80 −
90% of thespace-time volume, which makes the cluster algorithm not much better thana local update Metropolis algorithm. However, a combination of both thesealgorithms was used to study the physics of the QDM [53,54] and reliablyextract the phase diagram at zero temperature. For the QDM certain othercomputational methods have been reported [55,56], including certain noveleffective theory approaches [57]. ecent progress on cluster and meron algorithms for strongly correlated systems 17
Cluster algorithms are extremely efficient tools which are very useful to speedup the Monte Carlo sampling of strongly interacting systems wherever theyare applicable. These algorithms act by placing bonds between degrees of free-dom which are correlated, and build up patches in configuration space whichcan be flipped . The flipping process gives rise to a new configuration whichis very different from the parent one, but has the same Boltzmann weight.Thus, the correlation in between the subsequent Monte Carlo configurationsare removed, since the clusters themselves are physical degrees of freedom. Inthe case of fermions, configurations often come with negative signs associatedwith fermions interchanging their respective positions. A novel cluster algo-rithm — the meron algorithm, is able to analytically cancel clusters which areresponsible for the negative worldlines, and only operate in the Hilbert spacewhich contribute non-trivially to the partition function. Cluster algorithms forgauge theories are even more challenging than other systems since one has toidentify proper degrees of freedom which can be used to build clusters thatdo not occupy the whole space-time volume. In the case for certain quantumlink models, this has been achieved by dualizing the original Hamiltonian, andrewriting it exactly in terms of dual quantum height variables. The resultingcluster algorithm is every efficient. Further, for gauge theories, the Gauss’ Lawneeds to be implemented, which may or may not be possible with the clusterrules that allow the sampling the full Hilbert space. For example, in the caseof QLM where the vacuum is charge neutral, this is not a problem, but theQDM which has a staggered background charge this does not work so well, andthe clusters become large. This can be understood as the lack of a referenceconfiguration for the QDM, which is an essential ingredient for the success ofa cluster algorithm.Demanding the existence of reference configurations, it is possible to con-struct models which are guaranteed to have efficient cluster algorithm. SuchHamiltonians are often called designer Hamiltonians [58], and can be studiedin any space-time dimensions. Since the exact form of the microscopic interac-tions is unimportant to study physics associated with breaking of symmetriesacross thermal or quantum phase transitions, thanks to universality, cluster al-gorithms can be used to study a wide range of physical phenomena in naturallyoccurring strongly correlated systems.In the past few years, there has been renewed interest with cluster andmeron algorithms. In the context of classical spin models, the meron algorithmis being used to study the topological charge and the vorticity properties ofnon-linear sigma models [59]. For lattice fermions, the meron algorithm hasbeen used to study a Hamiltonian spin-half lattice fermions that displays sym-metry enhancement for certain coupling regimes [60]. In the context of gaugetheories, the dualization technique has been extended to a non-Abelian SU (2)quantum link model on the honeycomb lattice and an efficient cluster algo-rithm has been constructed in terms of the dual height variables [61]. This clus-ter algorithm operates on a four sub-lattice structure, where the clusters are built using the heights on one sub-lattice, while the heights on the other threesub-lattices decide the nature of bonds. This study identified more crystallineconfined phases in non-Abelian theories, as well as showed fractionalization ofa Z (2) flux.Another very novel and interesting development with cluster algorithms isthe simulation of strongly correlated systems in real-time and far out ofequilibrium . The first such study considered interesting initial states whichare in contact with a thermal reservoir, and are completely driven via dissipa-tion. The authors of [62] devised a cluster algorithm to simulate the Lindbladevolution in an anti-ferromagnetic initial state. More extensive studies revealedwhat kinds of measurements give rise to sign problems, and which measure-ments did not [63]. Transport phenomena in strongly correlated spin- drivenvia Lindblad dynamics was studied using cluster algorithms in [64]. Acknowledgements
This article is written in the memory of the late Pushan Majumdar.Pushan taught me the intricacies of the multi-level algorithm [65,66], a highly efficient MonteCarlo method which greatly accelerates the simulation of pure Wilson-type gauge theories.Besides, I would like to acknowledge Pushan for the many discussions on physics, his adviceon parallel and GPU programming that we have had, and for the many cups of coffee hewould prepare after lunch.I would like to thank Shailesh Chandrasekharan and Uwe-Jens Wiese for teaching meso much about cluster and meron algorithms. I would also like to acknowledge SaumenDatta, F.-J. Jiang, Kieran Holland, Emilie Huffman, Ferenc Niedermayer, Stefan Schaefer,Arnab Sen, Rainer Sommer, Urs Wenger, and Ulli Wolff for various discussion on improvedalgorithms.
References
1. J. Cardy,
Conformal Field Theory and Statistical Mechanics , Les Houches SummerSchool on Exact Methods in Low-Dimensional Statistical Physics and Quantum Com-puting (2008).2. A. Pelissetto, E. Vicari,
Critical Phenomena and renormalization group theory , Phys.Rept., 368, 549 (2002).3. M. Moshe. J. Zinn-Justin,
Quantum Field Theory in the Large-N limit: A Review , Phys.Rept., 385, 69 (2003).4. J. Maldacena,
The Large-N Limit of Superconformal Field Theories and Supergravity ,Intl. J. of Theo. Phys., 38, 1113 (1999).5. R. Rattazzi, V. Rychkov, E. Tonni, A. Vichi,
Bounding scalar operator dimensions in4D CFT , JHEP 12, 31 (2008).6. S. Hellerman, D. Orlando, S. Reffert, M. Watanabe,
On the CFT Operator Spectrum atLarge Global Charge , JHEP 17, 12 (2015).7. U. Schollwoeck,
The density-matrix renormalization group in the age of matrix productstates , Ann. Phys. 326, 96 (2011).8. N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller, E. Teller,
Equation of StateCalculations by Fast Computing Machines , J. Chem. Phys. 21, 1087 (1953).9. D. Landau, K. Binder,
A Guide to Monte Carlo Simulations in Statistical Physics , Cam-bridge University Press (2014).10. M. Newman, G. Barkema,
Monte Carlo Methods in Statistical Physics , Clarendon Press(1999).11. P. W. Kasteleyn, C. M. Fortuin,
Phase transitions in lattice systems with random localproperties , J. Phys. Soc. Jpn. Suppl., 26, 11 (1969).12. C. M. Fortuin, P. W. Kasteleyn,
On the random-cluster model. I. Introduction andrelation to other models , Physica 57, 536 (1972).ecent progress on cluster and meron algorithms for strongly correlated systems 1913. R. H. Swendsen, J. -S. Wang,
Non-universal critical dynamics in Monte Carlo simula-tions , Phys. Rev. Lett. 58(2), 86 (1987).14. F. Niedermayer,
General Cluster Updating Method for Monte Carlo Simulations , Phys.Rev. Lett. 61, 2026 (1988).15. U. Wolff,
Collective Monte Carlo updating for spin systems , Phys. Rev. Lett.62(4), 361(1989).16. F. Niedermayer,
Cluster Algorithms , Advances in Computer Simulations, Lecture Notesin Physics, Vol 501, 36 (1997).17. W. Krauth,
Cluster Monte Carlo Algorithms , New Optimization Algorithms in Physics,(2003).18. H. G. Evertz, G. Lana, M. Marcu,
Cluster algorithm for vertex models , Phys. Rev. Lett.70, 875 (1993).19. H. G. Evertz, M. Marcu,
The Loop-Cluster Algorithm for the Case of the 6 VertexModel , Nucl. Phys. Proc. Suppl. 30, 277 (1993).20. B. B. Beard, U.-J. Wiese,
Simulations of Discrete Quantum Systems in ContinuousEuclidean Time , Phys. Rev. Lett. 77, 5130 (1996).21. A. W. Sandvik, J. Kurkij¨arvi,
Quantum Monte Carlo simulation method for spin sys-tems , Phys. Rev. B 43, 5950 (1991).22. O. F. Sylju˚asen, A. W. Sandvik,
Quantum Monte Carlo with Directed Loops , Phys.Rev. E 66, 046701 (2002).23. H. G. Evertz,
The loop algorithm , Advances in Physics 52 (2003).24. N. V. Prokof’ev, B. V. Svistunov,
Worm Algorithms for Classical Statistical Models ,Phys. Rev. Lett. 87, 160601 (2001).25. W. Heisenberg,
On the theory of ferromagnetism , Zeit. f¨ur Phys., 49 (9), 619 (1928).26. M. Suzuki,
Relationship between d-dimensional quantum spin system and (d+1)-dimensional Ising systems , Prog. Theor. Phys., 56, 1454 (1976).27. M. Troyer, U.-J. Wiese,
Computational Complexity and Fundamental Limitations toFermionic Quantum Monte Carlo Simulations , Phys. Rev. Lett. 94, 170201 (2005).28. R. Blankenbecler, D. J. Scalapino, R. L. Sugar,
Monte Carlo calculations of coupledboson-fermion systems. I , Phys. Rev. D. 24 (8), 2278 (1981).29. F. Assaad, H. G. Evertz,
World-line and determinantal quantum Monte Carlo Methodsfor spins, phonons and electrons , Lecture Notes in Physics, Vol. 739, pp. 277–356, Springer,Berlin Heidelberg (2008).30. S. Duane, A. D. Kennedy, B. J. Pendleton, D. Roweth,
Hybrid Monte Carlo , Phys. Lett.B 195 (2), 216 (1987).31. S. Chandrasekharan, J. Cox, J. C. Osborn, U.-J. Wiese,
Meron cluster approach tosystems of strongly correlated electrons , Nucl. Phys. B 673, 405 (2003).32. S. Chandrasekharan, U.-J. Wiese,
Meron cluster solution of fermion sign problems ,Phys. Rev. Lett. 83, 3116 (1999).33. S. Chandrasekharan,
The Fermion bag approach to lattice field theories , Phys. Rev. D82, 025007 (2010).34. E. Huffman, S. Chandrasekharan,
Fermion bag approach to Hamiltonian lattice fieldtheories in continuous time , Phys. Rev. D 96, 114502 (2017).35. M. Boninsegni, N. V. Prokof’ev, B. V. Svistunov,
Worm algorithm and diagrammaticMonte Carlo: a new approach to continuous-space path integral Monte Carlo simulations ,Phys. Rev. E 74, 036701 (2006).36. K. van Houcke, E. Kozik, N. V. Prokof’ev, B. V. Svistunov,
Diagrammatic Monte Carlo ,Phys. Procedia 6, 95 (2010).37. R. Brower, P. Tamayo,
Embedded dynamics for φ theory , Phys. Rev. Lett. 62, 1087(1989).38. R. Ben-Av, H. G. Evertz, M. Marcu, S. Solomon, Critical acceleration of finite temper-ature SU(2) gauge simulations , Phys. Rev. D 44, 2963 (1991).39. W. Kerler, T. Metz,
Ising embedding for cluster algorithms in finite temperature SU(2)gauge theory , Phys. Rev. D 50, 7082 (1994).40. D. Horn,
Finite Matrix Models with Continuous Local Gauge Invariance , Phys. Lett.B 100, 149 (1981).41. P. Orland, D. Rohrlich,
Lattice gauge magnets: Local isospin from spin , Nucl. Phys. B338, 647 (1990).0 Debasish Banerjee42. S. Chandrasekharan, U.-J. Wiese,
Quantum Link Models: A Discrete Approach toGauge Theories , Nucl. Phys. B 492, 455 (1997).43. N. Shannon, G. Misguich, K. Penc,
Cyclic exchange, isolated states and spinon de-confinement in an XXZ Heisenberg model on the checkerboard lattice , Phys. Rev. B 69,220403 (2004).44. M. Hermele, M. Fisher, L. Balents,
Pyrochlore Photons: The U (1) Spin-Liquid in aS=1/2 Three Dimensional Frustrated Magnet , Phys. Rev. B 69 064404 (2004).45. L. Balents,
Spin liquids in frustrated magnets , Nature 464, 199 (2010).46. B. Schlittgen, U.-J. Wiese,
Low-energy effective theories of quantum spin and quantumlink models , Phys. Rev. D 63, 085007 (2001).47. J. B. Kogut, L. Susskind,
Hamiltonian Formulation of Wilson’s Lattice Gauge Theories ,Phys. Rev. D 11, 395 (1975).48. D. S. Rokhsar, S. A. Kivelson,
Superconductivity and the Quantum Hard-Core DimerGas , Phys. Rev. Lett. 61, 2376 (1988).49. S. Sachdev,
Spin-Peierls ground states of the quantum dimer model: A finite-size study ,Phys. Rev. B 40, 5204 (1989).50. R. Moessner, K. S. Raman,
Quantum Dimer Models , Lecture Notes, Trieste 2007(arXiv:0809.3051).51. D. Banerjee, F.-J. Jiang, P. Widmer, U.-J. Wiese,
The (2 + 1)-d U(1) quantum linkmodel masquerading as deconfined criticality , J. Stat. Mech. 1312, P12010 (2013).52. D. Banerjee, P. Widmer, F.-J. Jiang, U.-J. Wiese,
Crystalline Confinement , PoS LAT-TICE2013, 333 (2014).53. D. Banerjee, M. B¨ogli, C. P. Hofmann, F.-J. Jiang, P. Widmer, U.-J. Wiese,
Interfaces,Strings, and a Soft Mode in the Square Lattice Quantum Dimer Model , Phys. Rev. B 90,24, 245143 (2014).54. D. Banerjee, M. B¨ogli, C. P. Hofmann, F.-J. Jiang, P. Widmer, U.-J. Wiese,
Finite-Volume Energy Spectrum, Fractionalized Strings, and Low-Energy Effective Field Theoryfor the Quantum Dimer Model on the Square Lattice , Phys. Rev. B 94, 11, 115120 (2016).55. T. Oakes, S. Powell, C. Castelnovo, A. Lamacraft, J. P. Garrahan,
Phases of quantumdimers from ensembles of classical stochastic trajectories , Phys. Rev. B 98, 064302 (2018).56. Z. Yan, Y. Wu, C. Liu, O. F. Sylju˚asen, J. Lou, Y. Chen,
Sweeping cluster algorithm forquantum spin systems with strong geometric restrictions , Phys. Rev. B 99, 165135 (2019).57. J. Herzog-Arbeitman, S. Mantilla, I. Sodemann,
Solving the quantum dimer and sixvertex models one electric field line at a time , Phys. Rev. B 99, 245108 (2019).58. R. K. Kaul, R. G. Melko, A. W. Sandvik,
Bridging lattice-scale physics and continuumfield theory with quantum Monte Carlo simulations , Annu. Rev. Con. Mat. Phys. 4, 179(2013).59. W. Bietenholz, J. C. Pinto Barros, S. Caspar, M. Hornung, U.-J. Wiese,
Meron- andSemi-Vortex-Clusters as Physical Carriers of Topological Charge and Vorticity , PoS LAT-TICE2019 288, (2019).60. H. Liu, S. Chandrasekharan, R. K. Kaul,
Hamiltonian models of lattice fermions solvableby the meron-cluster algorithm , arXiv:2011.13208.61. D. Banerjee, F.-J. Jiang, T. Z. Olesen, P. Orland, U.-J. Wiese,
From the SU(2) quantumlink model on the honeycomb lattice to the quantum dimer model on the kagom´e lattice:Phase transition and fractionalized flux strings , Phys. Rev. B 97, 20, 205108 (2018).62. D. Banerjee, F.-J. Jiang, M. Kon, U.-J. Wiese,
Real-Time Simulation of Large OpenQuantum Spin Systems driven by Measurements , Phys. Rev. B 90, 241104 (2014).63. F. Hebenstreit, D. Banerjee, M. Hornung, F.-J. Jiang, F. Schranz, U.-J. Wiese,
Real-time dynamics of open quantum spin systems driven by dissipative processes , Phys. Rev.B 92, 035116 (2015).64. D. Banerjee, F. Hebenstreit, F.-J. Jiang, U.-J. Wiese,
Real-time simulation of non-equilibrium transport of magnetization in large open quantum spin systems driven bydissipation , Phys. Rev. B 92, 121104 (2015).65. M. L¨uscher, Peter Weisz,
Locality and exponential error reduction in numerical latticegauge theory , JHEP 0109, 010 (2001).66. P. Majumdar,