Percolation and cluster Monte Carlo dynamics for spin models

A general scheme for devising efficient cluster dynamics proposed in a previous letter [Phys.Rev.Lett. 72, 1541 (1994)] is extensively discussed. In particular the strong connection among equilibrium properties of clusters and dynamic properties as the correlation time for magnetization is emphasized. The general scheme is applied to a number of frustrated spin model and the results discussed.

A general scheme for devising efficient cluster dynamics proposed in a previous letter [Phys.Rev.Lett. 72, 1541] is extensively discussed. In particular the strong connection among equilibrium properties of clusters and dynamic properties as the correlation time for magnetization is emphasized. The general scheme is applied to a number of frustrated spin model and the results discussed.

I. INTRODUCTION
The cluster formalism introduced by Kasteleyn and Fortuin (KF) [1] and later developed by Coniglio and Klein (CK) [2] for the ferromagnetic Ising model has greatly enhanced the understanding of critical phenomena in terms of geometrical concept. Moreover, based on such a formalism Swendsen and Wang (SW) [3] have introduced a cluster dynamics which drastically reduces the critical slowing down in Monte Carlo (MC) simulation of ferromagnetic spin models. In the SW dynamics spins belonging to the same cluster are flipped in one step as opposed to the single spin dynamics, where spins are flipped one at time. The efficiency of the SW dynamics stems from the fact that in the KF-CK formalism the clusters represent correlated spins; therefore if one spin in a cluster is flipped all the other spins in that cluster will successively tend to flip coherently. Consequently by flipping in one move all the spins in the same cluster results in a much faster dynamics. The SW algorithm has been extended and applied efficiently to several unfrustrated spin models. [4][5][6][7][8] Unfortunately, the SW dynamics based on the direct extension of the KF-CK cluster formalism to frustrated spin systems does not show any reduction of the relaxation times. [9][10][11] The reason for that is due to the fact that now the KF-CK clusters do not represent any longer correlated spins. [12,13] Recently Kandel, Ben-Av and Domany have introduced a new type of MC cluster algorithm for the particular case of the fully frustrated Ising model [14] which is able to reduce the critical slowing down. Attempts to use their algorithm to other frustrated models has been satisfactory in few cases [15,16] and discouraging in others [17,18]. A different approach based on the definition of quasi-frozen clusters in spin glasses looks very promising, but the implementation of a related cluster dynamics has not been yet explored [19]. More recently, based on the approach of Kandel et al., we have proposed a general criterion and a systematic procedure to define clusters and related efficient cluster dynamics for frustrated spin models [20]. In particular, we have applied our general criterion to a class of fully frustrated Ising spin models on square lattices where the relative strength between the interactions can be varied. For any value of the relative strength, without invoking "ad hoc" algorithms for each case, the same general criterion generate a Monte Carlo dynamics which dramatically reduce the critical slowing down.
The aim of this paper is to illustrate in more details the criterion proposed in Ref. [20] and apply our method to a number of frustrated spin models. In section II, III we discuss the extension of cluster formalism to frustrated spin models. We stress that for the unfrustrated spin model the clusters percolate at a temperature T p that coincides with the critical temperature T c , while for the frustrated models T p is larger than T c . In section IV following the approach of Kandel et al. [14] we introduce a large variety of cluster definitions which contains as a particular case the KF-CK clusters. We then illustrate a procedure to approach systematically, in successive order of approximations, a cluster definition for which the percolation temperature T p becomes closer and closer to the critical temperature T c . In section V, VI and VII we check our procedure by comparing percolation quantities [21,22] and thermodynamic quantities on a variety of frustrated models using MC simulations. For a number of frustrated models without disorder we find to second order that the clusters percolate at a temperature T p numerically indistinguishable from the critical temperature T c . For disordered frustrated models like Spin Glass (SG) the convergence of T p towards the SG critical temperature T SG is very slow. In section VIII we implement a cluster MC dynamics based on the novel cluster definition. We show that the dynamics is very efficient with drastic reduction of the relaxation time for those frustrated systems, introduced in section V-VII, for which the cluster definition leads to T p ≃ T c . In the appendix we show that the cluster dynamics fulfills detailed balance and briefly discuss the ergodicity problem.

II. THE CLUSTER APPROACH IN FRUSTRATED SYSTEMS.
It is well known that he partition function of a ferromagnetic Ising model can be written in terms of the clusters of an equivalent percolation model [1,2]. A similar result can be also obtained for Ising systems where frustration is present [12]. The aim of this section is to recall and discuss those results which will be useful in the following.
Let us consider the Ising Hamiltonian where ǫ ij = ±1 is the sign of the quenched interaction and J ≥ 0 the interaction modulus. The interaction configuration {ǫ ij } can be a deterministic periodic structure or a disordered one. Hamiltonian (1) is said to contain frustration if at least one closed path L exists such that i,j ∈L ǫ ij = −1. Within the CK approach [2,12,13] we introduce "bonds" between nearest neighbor (nn) spins satisfying the interaction probability p = 1 − e −2βJ . The weight for each given configuration of spin {S i } and bond C is given by where C F (C A ) is the subset of C which covers ferromagnetic (antiferromagnetic) bonds, i.e. bonds with ǫ ij = 1 (ǫ ij = −1), and p = 1 − e −2βJ . The product takes into account the fact that bonds can link only spins which satisfy the interaction. The clusters are defined as maximal sets of spins connected by bonds. It can be shown [13] that in this case the partition function becomes where * C means that the sum is over the bond configurations which do not contain frustration. Furthermore, the presence of positive and negative interaction implies that the following relations can be proved and where γ ∞ i↑ (γ ∞ i↓ ) is 1 if the spin i is "up" ("down") and belongs to an infinite cluster, otherwise it is 0, and γ ij (γ ij ) is 1 if the spins i and j belong to the same cluster and are parallel (antiparallel), otherwise it is 0, and where ... is the usual thermodynamic average for a fixed configuration of interaction {Jǫ ij } and ... CK is the average over spin and bond configurations weighted with (2). The percolation quantities are instead given by where γ ∞ i = γ ∞ i↑ + γ ∞ i↓ and γ ij = γ ij + γ ij , and P i is the probability that the spin at site i belongs to the ∞ cluster and P ij is the probability that the spins at site i and j are in the same cluster.
It is clear that without frustration the relations valid in the ferromagnetic case are recovered with only trivial differences; in fact, in this case, eqs. (4,5) become (8) and Examples of models without frustration for which eqs. (8) and (9) are satisfied are the antiferromagnetic Ising model on a square lattice and an Ising model with interaction J i,j = Jσ i σ j where σ i = ±1 are quenched variables. It is easy to realize that for any fixed configuration of σ i , although the interactions may be positive and negative, all the loops are unfrustrated.
From (4) and (7) it follows that unlikely the ferromagnetic case the critical temperature does not coincide with the percolation temperature. In fact, defining the critical temperature T c as the temperature at which the Edwards-Anderson [24] order parameter q EA = 1 N S i 2 vanishes and the percolation temperature T p as the temperature at which the percolation probability ρ ∞ = 1 N i P i vanishes, from eq. (4) follows that T c < T p [12]. (In the definition of q EA the bar stands for the average over the configurations of interaction {Jǫ ij } and N is the number of spins. [25]) This result has been verified numerically for a number of different systems as 2d and 3d spin glass [9,10], fully frustrated model [10] and frustrated XY model [11].
Similarly from eq. (5) correlation and connectivity do not coincide any more. In fact two spins instead of being in the same cluster may be parallel in one configuration and antiparallel in another (see Fig. 1). Although these two configurations will both contribute to the connectivity they will interfere and strongly reduce correlations. Therefore for T > T c defining g ij = S i S j it follows from (5) and (7) |g ij | ≤ P ij .

III. GENERALIZATION OF THE CLUSTER APPROACH.
The main result of the previous section is that when frustration is present the KF-CK clusters do not represent anymore correlated spins. As a consequence the percolation temperature T p is higher than the thermodynamic critical temperature T c . We will generalize now the KF approach [1] in order to define new clusters which represent correlated spins even in the frustrated case. The new clusters will reduce to the usual KF-CK clusters in the unfrustrated case.
We will achieve this goal in two steps. First, in this section, following the approach introduced by Ben Av at al. [14,26] (BKD), we will consider a large class of clusters which contain as a particular case the KF-CK clusters. Second, in the following section, we will give a criterion to choose systematically the "right" clusters in successive order approximations in such a way that pair connectness and pair correlation functions tend to coincide. We will consider an approximation good enough when the percolation temperature T p has approached the critical temperature T c .
Let us consider a square lattice with ferromagnetic and antiferromagnetic interactions [27] and let us focus our attention on a single isolated plaquette. This can be either unfrustrated or frustrated. [28] Following Ref. [ [14]], it is possible to generalize the KF procedure by assigning globally to each bond configurations on a given plaquette or even on larger spin blocks taken as "elementary units" (see Fig. 2)] its own bond probability. Within each elementary unit or plaquette each bond configuration probability is independent on each other, and thus, it is not given, as in the CK formalism, by the product of bond probabilities of single pairs of spins.
We consider, to begin with a specific case, a checkerboard partition of the square lattice (see Fig.3) and we take one of the two sets of plaquettes (the shaded or the unshaded ones) as the set of "elementary units" on which we make our independent choices. We proceed further by generalizing the KF approach. Consequently we "dilute" the couplings on the plaquette, replacing them with a set of new interaction configurations which contain only J ′ → ∞ or J ′ = 0 interactions. [23] We consider the generic hamiltonian (1). In sec. V-VII we will consider specific examples. The possible interaction configurations are shown in Fig.4. We assign a bond weight w 0 to the interaction configuration c 0 of the first row (no bonds), we assign the same weight w 1 = w 2 ...= w 4 to all the configurations of bonds, c 1 ...c 4 in the second row, and so on. The symmetry of the plaquette allows us to choose the same weight for symmetric configurations, i.e. members of the same row in Fig.4.
The requirement of the equality between the partition function of the original model and that of the diluted model gives <i,j>∈plaq.
where the sum is over all the possible interaction configurations on the plaquette (M = 15 for the example in the Fig.4). With < i, j >∈ c α we mean the n.n. spins connected by J ′ interactions in c α and with < i, j >∈ plaq. the n.n. spins on the plaquette. The previous relation can be rewritten as and are the energy of the plaquette viewed as "elementary unit" in the new interaction configuration c α and in the original system. Of course such procedure can be repeated for every Hamiltonian that can be written as sum of elementary unit energies, i.e.

H({S
whose block l we dilute and obtain the stochastic Hamiltoniañ with C = l c α(l) . The equivalence between the original model and the diluted one is obtained imposing (cfr. eq.(11) ), where the w α(l) is weight with which the configuration c α(l) on the lth elementary unit occurs. Furthermore, the partition function can be written as Spins that are connected by infinite strength interaction are frozen while the others do not interact. Thus the dilution of the original Hamiltonian is also called a "freezing and deleting operation". Of course we have: e −βH l ({Si},cα) = δ {Si},c α(l) (18) where δ {Si},cα is 1 or 0 depending whether or not the spin configurations satisfy all the ∞ strength interactions in the interaction configuration c α (l) of the l-th plaquette. Two spins connected by an infinite strength interaction will be frozen in the configuration which satisfy the interaction. On the entire lattice we can define: where C F and C A are defined like in eq.(2). Let us observe that where n α (C) is the number of elementary units on which we have chosen the αth configuration (with α = 0, ...M ) in the given interaction configuration C (with α n α (C) = N u where N u is the number of elementary units). Therefore, from eqs. (17), (18) and (20), we obtain where * C stands for the sum over all the interaction configurations that are compatible with at the least one of the possible spin configurations (i.e. all the unfrustrated graphs). Eq. (21) is the generalization of eq.(3) which can be recovered considering a pair of n.n spins as elementary unit.
Equation (21) can be also obtained following the CK approach where the clusters are defined in the original system introducing fictitious bonds between spins satisfying the interaction. Given a spin configuration {S i }, the probability to realize a configuration of bonds c α(l) on each unit l is given by where w α satisfy eq.(16). Due to (16) and (18)these probabilities are normalized for any spin configuration α(l) P (c α(l) |{S i }) = 1. The Kronecker delta assures that the bonds are thrown only between spins satisfying the interaction. For the entire system the weight for a given configuration of spins {S i } and bond configuration C = l c α(l) is given by where (14) have been taken into account. Finally from (19) and (20) we have Summing over the spin and bond configurations we recover eq.(21). The advantage of this approach is to make clear that both spins and bonds can be defined in the original system where the clusters are defined as maximal sets of spins connected by bonds. To calculate the statistics of the CK clusters we have to generate equilibrium spin configurations first. Then, for each equilibrium configuration we can assign a bond configuration on each unit l with probability given by (24).
Summarizing, following the approach of BKD, we have defined a vast class of generalized percolation models equivalent to our original spin model eq.(1). When one reduces the elementary unit to a single spin pair one recovers KF-CK solution and that for larger units it is always possible to find a solution in the form of product of KF probability p = 1 − e −2βJ , or, generally, solutions that are factorization of probabilities for a subpartition of the elementary units.
The generalization discussed does not solve automatically our problem of recovering the identity between cluster connectivity and spin correlation function, eq.(9), when frustration is present. However, the great freedom given by eq.(16) still gives hope to find solutions for w α(l) , in such a way to achieve the equality (9) at least in an approximate way. Then it is crucial to find a criterion to select among the many different possibilities offered by eq.(16), those for which (9) holds. This is the second step needed to achieve our aim and will be discussed in the next section.

IV. CONDITIONS BETWEEN CORRELATION AND CONNECTIVITY.
For each percolation model defined in the previous section (satisfying eq.(16)), it is straightforward to generalize the relations (4) and (5) provided that now the average over the spin and bond configurations has to be computed with weights given in eq. (24). In particular relation (5) is also valid when we consider a sub-system made of a single unit l, namely for each i and j on the unit l. Here where the sum is over all possible spin and cluster configurations on the unit l and P (c α |{S i }) is given by (22). When the quantity to average is function only of the spin variables like s i s j l , due to eq. (22), eq.(26) simplifies Our aim is to find among the large class of solutions of eq.(16) a solution for w α(l) in such a way that the equality (9) is satisfied. In this way the percolation temperature T p would coincide with T c and the clusters so identified will be characterized by percolation critical exponents equal to thermodynamical critical exponents. Since this is not a trivial task, we seek solutions which fulfill approximately eq. (9) by imposing the condition on a sub-system made of a single unit, namely To find the solution of (16) and (28) may still be complicated due to the large number of unknown (2 n where n is the number of edges in the unit cell). Technically the problem can be simplified if we make use of eq. (25) which is always valid if eqs (16) are satisfied.
Eq. (28) together with eq. (25) implies for each pair ij in the elementary unit If (29) and (16) are satisfied, from (25) follows that also (28) is satisfied. As we will show in the next section in a specific example it is easier to impose (16) and (29) than to impose the equivalent conditions (16) and (28) .
It is clear that the larger the unit the better eq.(9) is satisfied. It also clear that there is no guarantee on how good the different approximations are and how fast eq.(9) is approached by increasing the unit size. However, this procedure, by increasing the size of the unit, allows a systematic way to improve the approximation.
Eqs. (16) and (28) introduce a set of independent conditions whose number depends on the size and the symmetry of the chosen elementary unit. In general we are not able to know if the conditions introduced by (16) and (28) have a solution and if it is unique. Of course only solutions such that w α(l) > 0 are acceptable and these conditions introduce further restriction.
In general, there are two possibilities. The first possibility is that there are no solutions which satisfy eqs. (28). In this case we can relax eq.(28) by imposing where r ij is the distance between spin i and j. Choosing k M in an appropriate way it is possible to reduce the number of conditions until a solution is found. If k M ≥ 2 we believe that the solution is rather reasonable since the conditions of Sec. II are satisfied on the unit. The second possibility is that the solution is not unique. In this case we expect only small differences among different solutions.
To show how this scheme works, in the following sections we will analyze a number of frustrated systems. In all cases we are always able to satisfy conditions (16) and (28).

V. A DECORATED ISING MODEL.
In order to check our approach we have considered a decorated Ising model with frustration. Starting with an Ising model on a square lattice we introduce between each pair of n.n. spins S i , S j on the square lattice two extra spins S k , S l (see Fig.5) modifying the interaction from to this generic set of four spins will be the elementary unit l. For simplicity in eqs. (31) and (32) and in the following of this section we omit the label l.
The partition function of this model is reducible to the Ising one via a "decimation" on spins S k , S l : . The critical temperature can therefore be calculated exactly T c = 2.24.... Monte Carlo estimation of the percolation temperature for the unmodified KF-CK clusters on such a system gives T p − T c ≃ 0.2 (here and in the following temperature are expressed in units of J/k B ). As expected the presence of frustration prevents the coincidence of percolation and thermodynamic properties for the unmodified KF-CK clusters.
We have to solve eqs. (16) and (29) for the unknown weights w α where α labels the bond configurations in the elementary unit (Fig. 5). The average in eq. (29) is over all spin and bond configurations with probability given by (26) where H l is given by (32).
The spin correlations can be easily calculated from (27) since they does not require the knowledge of the w α . We can immediately find (Fig.5) S k S l l = 0.
We can disregard bond configurations by inspection. For example the weight of a bond configuration which connects i and j through site l must be zero. In fact this bond configuration would correspond to S i and S j antiparallel resulting in γ ij l > 0 contrary to eq. (29). By imposing eqs. (29)-(35) we reduce the number of possible bond configurations to 12. Furthermore 3 of them have the same connectivity properties (i.e, they connect the same sites) as the configuration α = 3 (Fig. 6). Therefore they can be disregarded reducing the number of weights different from zero to 9. Eqs. (16) now read where u = e −2βJ . The structure of such an equation system still allows us to chose one unknown arbitrary; in particular we have imposed w 4 = 0 because it provides w i ≥ 0 ∀T .
Then the solution is In order to calculate the percolation temperature of the clusters defined by eq.(37), T p , we proceed in the following way. Given a spin configuration {S i } we assign to each plaquette l a bond configuration with the probability P (c α |{S i }) provided by eq. (22) and (37). Then we obtain clusters defined in the entire lattice as maximal sets of spins which are connected by bonds (CK-like cluster definition). It is, then, possible to measure percolation quantities.
We have estimated T p via a data collapsing of the probability P (T, L) of having a percolating cluster at temperature T in a system of size L (Fig.7). We have simulated the model with both standard spin-flip MC dynamics and the cluster MC dynamics which we will discuss in sec. VIII.A obtaining indistinguishable results within our numerical precision. In Fig. 7 we have reported the results obtained with the latter dynamics. Using the scaling "ansatz" that near the percolation point P (T, L) = f ((T − T p )L 1/νp ), where the functional shape of f is unknown, we have found that T p ≃ 2.25 and ν p ≃ 1 consistent with the critical temperature T c = 2.24 and with the critical exponent ν = 1 of a ferromagnetic Ising model.

VI. THE ASYMMETRIC FULLY FRUSTRATED MODEL (AFF).
Let us consider a less trivial model, the frustrated Ising model on a square lattice with periodic boundary conditions where each plaquette contains three equal ferromagnetic interactions J and one antiferromagnetic interaction −XJ (0 ≤ X ≤ 1). This model interpolates between the Fully Frustrated (FF) model (X = 1) and the diluted Ferromagnetic Ising model (X = 0) [29].
If we take w α according to the definition given in Fig.8 eqs. (16) specified to this case give The number of unknowns is larger than the number of equations. The percolation temperatures associated to these solutions will be, in general, higher than the thermodynamic critical temperature. In order to have these two temperature as close as possible and the cluster connectivity representing as better as possible the spin correlation, as discussed in sec. IV, we can impose either conditions (28) or more simply eq. (29). For pedagogical reasons we choose here eq. (28). These are three independent equations Eqs. (38) and (39) have the solution [29] and w 3 = w 7 = w 9 = 0. The equations (38) and (39) do not provide a unique solution. In fact a general solution can be found choosing w 3 as a free parameter. However, the further requirement w α ≥ 0 leads to w 3 = 0. The solution changes form for a temperature T * such that u * = e −2J/KB T * satisfies the equation (1 + 3u 2 )u X − u 3 − 3u = 0. This is due to the fact that at T = T * the correlation between spins linked by the interaction −XJ (let's call them S 1 and S 2 ) changes sign leading to different possible bond configurations. For example, the configuration associated to the weight w 5 which links antiferromagnetically S 1 and S 2 vanishes when S 1 S 2 > 0 [30]. It is interesting to specify the general solution for the diluted ferromagnetic Ising model (X = 0) and the FF model (X = 1).
In the first case (X = 0) it results u < u * for any temperature and we get the following non zero weights These weights reproduce the original KF-CK solution for a plaquette; in fact they are in the form of products of u and 1 − u which are the KF-CK weights for a single spin pair. Since for X = 0 the model reduces to a ferromagnetic Ising model where some interactions have been put equal to zero, the original KF-CK solutions reproduces the right clusters.
In the case X = 1 it results u > u * for any temperature and one obtains the following non zero weights which are in agreement with the cluster structure used in Ref. [31]. It is worth while to note that in this limit all the bond configurations which connect spins on opposite corners have weights equal to zero preventing the four spins of the plaquette to belong to the same cluster even at T = 0. In order to check if the percolation model we have defined has the expected properties, i.e. if T p = T c , we have studied the percolation and spin properties of the system with both spin-flip standard MC dynamics and the cluster MC dynamics which we will discuss in sec. VIII.B obtaining indistinguishable results within our numerical precision. In Figs. 9 and 10 we have reported the results obtained with the latter MC dynamics averaging at least over 6 · 10 4 MC sweeps after discarding the first 10 4 .
In Fig.9 we show the data collapsing for the mean cluster size defined as S = s n s s 2 where n s is the number of clusters of size s and the sum is over finite clusters for three values of X (X = .5, X = .75, X = 1). For comparison in Fig.10 we show, for the same X values, the data collapsing for χ = ( M 2 − |M | 2 ) where M is the magnetization. From those data we extract the percolation temperature, T p (X), the critical temperature T c (X) and the critical exponents γ p (X), ν p (X), γ(X), and ν(X) for any X value. Summarizing we find for the thermodynamic quantities. We note that T p (X) ≃ T c (X) within the estimated numerical precision. Furthermore the same critical exponents control the divergence of percolation and spin properties γ p (X) = γ(X) and ν p (X) = ν(X) for all X values but X = 1. In the X = 1 case we find ν(1) = ν p (1)  We have analyzed the cluster structure at T ≃ T p (X) ≃ T c (X) [20,35]. For X = 1 we have found that a typical configuration of critical clusters is a fractal made of a backbone and dangling bonds. The backbone is made of links and blobs as found in the ferromagnetic Ising model [34] where the fractal dimension of the entire cluster was found to be equal to D = 1 2 ( γ ν + d) which for dimension d = 2 gives D = 15/8 and the fractal dimension of the links or red bonds was found equal to D R = 13/24.
For the symmetric fully frustrated model X = 1 the structure changes drastically. In fact all the clusters are made of self-avoiding chains with fractal dimension given by the scaling relation D = 1 2 ( γp νp + 2). Using the numerical result γ p ≃ 2.0 and ν p ≃ 2.0 we find numerically D ≃ 2 in agreement with the result of Ref.BADK,Coddington,Kerler that for T = 0 predicts two percolating self-avoiding chains which fill up the entire system. It is interesting to note that if the approximation could be improved such that the condition (9) would be exactly satisfied we would obtain γ = γ p = 7/4 and ν = ν p = 1. In the plausible event that the exact clusters are still self-avoiding chains the fractal dimension would be D = 7/4 identical to the fractal dimension of a self-avoiding random walk at the θ point [37].

VII. ISING SPIN GLASS.
The most complex and interesting model in the class of spin systems described by the Hamiltonian (1) corresponds to the case in which ǫ ij is a random variable: this introduces quenched disorder together with frustration. We have studied the ±J Ising Spin Glass, in the case of a square lattice with the probability distribution p(ǫ ij ) = κδ ǫij ,−1 + (1 − κ)δ ǫij ,1 where κ is the concentration of antiferromagnetic interactions. The phase diagram of such system was described by Ozeki [38] and exhibits at low temperature a paramagnetic-ferromagnetic transition if the concentration of antiferromagnetic interactions κ is enough diluted (i.e. if 0 ≤ κ ≤ κ 0 , with κ 0 ∼ 0.1), otherwise there is a spin-glass transition at zero temperature (i.e. if κ 0 ≤ κ ≤ 0.5). The case κ = 1/2 correspond to the Edwards and Anderson model (EA). [24] As in the AFF case, we partition the lattice in square plaquettes of four spins. The system is then characterized by two kind of plaquettes: frustrated and unfrustrated. Frustrated plaquettes are those with one or three antiferromagnetic interactions [28]. We have analyzed such cases in section VI where we have obtained the probabilities P (c α |{S i }) for the bond configuration, c α , given a configuration of spin {S i } (see eq. (22) and (46)). We note that the weights given in eq.(46) for a frustrated plaquette containing three ferromagnetic and one antiferromagnetic interactions only depend on the interactions satisfied. Therefore the same weights can be used for the frustrated plaquette containing an odd number of antiferromagnetic interactions. We have used the unmodified KF-CK weights p = 1 − e −2Jβ for the unfrustrated ones (eq.(45)).
In order to compare the results obtained using the plaquette of four spins as a unit (we call this choice 2nd order approximation) with those obtained by using the single pair of spins as unit (1st order approximation) we have simulated the model with both standard spin-flip MC dynamics and the cluster MC dynamics which we will discuss in sec. VIII.C. The percolation quantities which we have measured, T p , γ p and ν p , do not depend, within our numerical accuracy, on the dynamics used.
In the case of EA model (κ = 0.5), we have found a percolation temperature T (2) p ≃ 1.40 higher than the critical temperature T c = 0, but lower than that obtained with the unmodified KF-CK clusters where T (1) p ≃ 1.80. [10] We have also estimated the percolation critical exponents ν p and γ p via a data collapsing (see Fig.11), obtaining the values ν p ≃ 1.33 and γ p ≃ 2.36, which are consistent with the random bond percolation values ν p = 4/3 and γ p = 43/18.
We have also studied the region of low κ where T c is finite and the transition is ferromagnetic. [38] In Fig.12 we show T (1) p (κ), T (2) p (κ) and T c (κ) for 0 ≤ κ ≤ 0.1 and κ = 0.5. The values have been obtained after a data collapsing for the mean cluster size S (eq. (19)) and susceptibility χ, respectively. It's clear from Fig.12 that the percolation temperature T  (2) p (0.5) ≃ 1.4, while T c abruptly goes to zero. [38] From this analysis it comes out that neither the unmodified KF-CK clusters nor our clusters are able to correctly represent spin correlations in Spin Glass systems. However since the percolation temperature T one might expect that a systematic improvement can be obtained if larger elementary units are used (Fig.2).

VIII. MONTE CARLO DYNAMICS ASSOCIATED TO PERCOLATION MODELS
It is possible to apply the general definition of cluster given above, to develop general MC cluster dynamics. The clusters are constructed assigning to each elementary unit one of the possible bond configurations according to the probability given in eq. (22). Then the usual SW algorithm can be applied to the clusters described above. Following Ref. [14] it can also be proven that detailed balance holds (see Appendix).
A. Decorated Ising model.
We partition the original square lattice in plaquettes as described in section V, and use the clusters defined there. With that cluster definition we have implemented the SW generalized cluster dynamics.
We have estimated the percolation temperature T p ≃ 2.25 and the percolation critical exponents γ p ≃ 1.77 and ν p ≃ 1. The values obtained are indistinguishable from those obtained by using spin-flip MC dynamics reported in sec. V. We have also estimated the corresponding thermodynamic quantities T c ≃ 2.24, γ ≃ 1.78 and ν ≃ 1.05 which are consistent with a ferromagnetic Ising critical point within our numerical accuracy.
In order to study the relaxation times we have computed the time dependent magnetization correlations where M (t) is the magnetization at time t. Using the new cluster dynamic we find a dramatic reduction of the slowing down which is present for the standard SW dynamics. The new dynamics has critical autocorrelation times of about 10 MCS (Monte Carlo Step per Spin), whose order of magnitude is comparable to those of the standard SW dynamics for a ferromagnet of the same size at criticality. On the contrary standard SW algorithm on Decorated Ising model shows very large correlation times at T c (see Fig.13). Our approach, then, reduces the critical slowing down in this system when compared to standard SW and local spin flip dynamics. B. AFF model.
We partition the original lattice in elementary units as described in section VI and use bond configurations and the associated probabilities introduced there to define clusters. With such cluster definition we have then implemented a SW generalized cluster dynamics. We have estimated percolation quantities T p (X), γ p (X) and ν p (X) for X = 0.5, 0.75, 1.0 obtaining values indistinguishable from those calculated by using spin flip MC dynamics. We have also computed the critical temperature T c (X) and critical exponents γ(X) and ν(X).
In order to study the relaxation times of our SW generalized cluster dynamics we have calculated the magnetization correlation function (53) versus Monte Carlo sweeps at the critical temperature T c (X). It shows a dramatic reduction for all the X values which we have studied (X = 0.5, 0.75, 1.0) with respect to SW unmodified CK-KF cluster and local MC dynamics. For a quantitative analysis of the critical dynamic exponent z(X) at T c (X), we have calculated the integrated autocorrelation time τ using the self-consistent procedure suggested in Ref. [39] with a window equal to 6. This method allows us to calculate τ for different system size L in a consistent way. We found a power law scaling of the form τ = kL z as shown in Fig.14. The estimated values of z(X) for X = 0.5, 0.75 are z(0.5) ≃ 0.30 and z(0.75) ≃ 0.46. The result definitely shows a strong systematic reduction of the critical dynamical exponent compared with those of standard MC dynamics. These results seem to indicate that the criterion to let T p (X) as near as possible to T c (X) (we have previously showed the coincidence of this two temperature for these models), allows to individuate efficient cluster dynamics. It is worth noting that even in the case X = 1 when T p = T c = 0, ν p = ν = 1 and γ p > γ the cluster dynamics exhibits a drastic reduction of the critical slowing down. Nevertheless, our analysis suggests that, for X = 1, it is possible improve further this result considering larger units as starting point of the proposed procedure.

C. ± J Ising SG model
The panorama is more variegated in the more complex case of Ising Spin Glass with varying ferromagnetic interactions concentration κ. This model exhibits in 2d a paramagnetic-ferromagnetic transition for 0 ≤ κ ≤ 0.1 and a spin-glass transition for 0.1 ≤ κ ≤ 0.5. Analogously to the other presented cases the cluster dynamic is realized by using the clusters defined in sec.VII. For frustrated plaquettes we have used the probabilities calculated in sec.VI for X = 1 (eqs. (46) and (22)), while for unfrustrated plaquettes we have used unmodified KF-CK clusters (eqs. (45) and (22)).
To check our simulation, we measured thermodynamic functions as energy E and specific heat C v . [40] They reproduce known data in literature up to a temperature, T f , under which our MC cluster dynamics freezes. We have also estimated T p , γ p and ν p obtaining good agreement with the values obtained by using spin flip MC dynamics.
We have already noted that in Ising Spin Glasses the percolation temperature of unmodified KF-CK clusters, T (1) p (κ), is higher than the percolation temperature, T (2) p (κ), of the clusters defined in section VII: we have where T c is the thermodynamical critical temperature. The equality holds only for κ = 0 and 1 (ferromagnetic and antiferromagnetic case respectively). It is possible to summarize the results in this way: as κ departs from the ferromagnetic Ising model κ = 0, relaxation times for temperature close to the critical temperature T c (κ) get longer, and in the region where ferromagnetic phase disappears (0.1 ≤ κ ≤ 0.5), they become extremely long, even if always shorter than those of both standard SW cluster dynamics (unmodified KF-CK clusters) and local spin flip dynamics.
Along the paramagnetic-ferromagnetic transition line (i.e. at T c (κ) with 0 ≤ κ ≤ 0.1) we have estimated the critical autocorrelation time τ (κ), defined as the time to reduce square magnetization correlation to 1/10 of its value at t = 0. These results are shown in Fig.15 for a square lattice of size L = 32.
In the region where the ferromagnetic phase disappears (0.1 ≤ κ ≤ 0.5) and the SG transition takes over at T SG = 0, our simulations get worse. We have studied for the case κ = 0.5 the following relaxation function as a function of time (Monte Carlo Step) for systems whose size is L = 80, 90, 100. Due to very long autocorrelation times we were able to perform simulation up to T f ≃ 0.8. Averages in eq.(55) were taken over 1 ÷ 4 · 10 4 MCS discarding the first 5 · 10 3 MCS. We observed that relaxation time of our cluster dynamics above T f is at least one order of magnitude lower than that of a standard spin flip dynamics (see Fig.16 and, for comparison, Ref. [41]).
In conclusion in the case of spin glass we see that to a lowering of the difference |T c − T p | corresponds a reduction of the relaxation times. However, there are indications that such a reduction exists only for T > T f . Our results suggest that taking a partition of the lattice made by larger "elementary" units the procedure we have discussed define clusters whose percolation temperature is closer to the critical temperature of the original spin model. The associated cluster dynamics is expected to be characterized by shorter autocorrelation times. Work is in progress in this direction.

IX. SUMMARY.
In this paper we have discussed a general scheme for devising efficient MC cluster dynamics for spin models. The scheme is based on three main steps. The first one consists in choosing a partition of the lattice into "elementary units". Then, using a method first introduced by Kandel, Ben-Av and Domany [14] which is based on independent choices on each elementary unit, it is possible to define a vast class of cluster models whose Free Energy is identical to the original spin model. Finally, among the many cluster models it is possible to choose the one which satisfies at the best the equality between cluster connectivity and spin correlation (cfr. eqs. (9) and (28)). This procedure defines clusters which can be used to implement a MC cluster dynamics. We have applied this method to a number of 2d frustrated spin models taking as elementary unit a single plaquette. We show that every time T c ≃ T p (i.e. the thermodynamic critical temperature of the spin model T c is equal to the percolation temperature for the equivalent percolation model) the associated MC cluster dynamics is characterized by very small autocorrelation times and a critical dynamic exponent z much smaller than the one obtained in local (Metropolis) MC dynamics. We have also shown that in the more complex case of a spin glass where disorder is added to frustration the percolation temperature T p results larger than spin glass temperature T SG and the percolation critical exponent are consistent with random bond percolation exponents. However we see that in this case the percolation temperature can be decreased up to T p ≃ 1.4 using a lattice partition based on a single plaquette. This result suggests that taking larger elementary units, like the ones in Fig.2d, T p can be further reduced. Then, we still find a lowering of the autocorrelation time for T > T f ≃ 0.8 but, this time, there are no indication of a lower z compared to standard Metropolis dynamics.
In conclusion our procedure allows for a systematic decrease of the autocorrelation times and, therefore, may serve as a general framework for the development of efficient MC dynamics in frustrated spin models.

APPENDIX
The aim of this appendix is to show that the MC dynamics defined in sec. VIII satisfies detailed balance. Following Ref. [14] we show that provided the mapping from H toH (see eq.(16)), a MC dynamic which verifies detailed balance principle forH verifies it also for H. In particular after executing a freezing and deleting operation on the original spin system, we can implement with the built clusters a cluster dynamic based for example on random flipping of independent clusters as in SW procedure. This is possible because such a dynamics certainly satisfies the detailed balance principle and is generally ergodic at finite temperature.
To prove that detailed balance is respected, let us make the following preliminary considerations. We can rewrite the relation (16) as This is the normalization condition for the conditioned probability to have the bond configuration c α(l) on the l-th elementary block, given the spin configuration {S i } on the system, i.e. the (56) expresses the normalization condition for the Since the choices on the elementary block are independent, the probability of the bond configuration C on the whole system, given the spin configuration {S i }, is the product P (C|{S i }) = l P (c α(l) |{S i }).
To obtain detailed balance principle, we must impose the following condition: may then be written as is the transition probability associated to the dynamic that we use on the dilute system with the HamiltonianH (for example this dynamic may be the simple random flipping of independent clusters). Generally let's suppose that theT C ({S ′ i } → {S i }) respects the detailed balance principle, i.e.
Therefore, from the (59) and the definition of P (C|{S i }), we have that, using the (60), becomes Now, multiplying and dividing for e −βH({Si}) , we obtain and by (59) and the definition of P (C|{S i }), we get This expression is the (58) and therefore the validity of the principle is demonstrated [14]. Summarizing, the main assumptions underling this proof consist in supposing that we can write the Hamiltonian as sum on elementary blocks H = l H l , the choices on the elementary blocks are independent from each other and we are using a dynamics for the dilute systemH that respects the detailed balance principle.
In particular a generalized cluster MC dynamics may be implemented with the following steps: individuate the clusters with "freezing and deleting" (i.e. to map H intoH); random flip of such clusters (this move certainly verifies detailed balance because clusters are not interacting inH), and iterate.
About the ergodicity we can say that the cluster dynamics here described is certainly ergodic for every finite temperature, because the probability to go from a given spin configuration to any other is always different from zero for every non zero temperature. In general ergodicity at zero temperature is difficult to prove: it must be checked specifically in each particular case. Nevertheless it is possible to guarantee ergodicity also in such extreme conditions, alternating cluster moves with a dynamics that certainly is ergodic at this temperature, without changing the qualitative features of the cluster dynamics [14].
In conclusion we have proven that adopting the proposed mapping of Hamiltonians, and in particular for our general definition of clusters, it is possible to develop MC dynamics which verify detailed balance principle.       6 The bond configurations of the decorated Ising model (sec. V) whose weights are different from zero (eq.57). We assign the same weight to all the elements belonging to the same group. In the figure each group is identified by a curl bracket. The conventional representation of interactions is the same as in the Fig.5.   (60)). In the figure we show only one element for each group; the other elements can be obtained trivially conserving the number of ferromagnetic and antiferromagnetic bonds. In the brackets it is shown the number of bond configurations belonging to the specific group. Full (dashed) lines are bonds between ferromagnetically (antiferromagnetically) interacting spins. Fig.9 The data collapsing for mean cluster size S of AFF models for systems with X = 0.5, 0.75, 1.0 and number of sizes L. The data have been obtained by using the cluster dynamics discussed in sec. VIII.B. Fig.10 The data collapsing for susceptibility for the AFF model, with the same values of X and L as in Fig.9. Fig.11 The data collapsing for mean cluster size S for EA Spin Glass model for systems with size L = 48, 64, 80, 100. Critical exponents and temperature are also reported. Fig.12 The critical temperature T c (κ) , the percolation temperature T