A Cluster Method for the Ashkin–Teller Model

Shai Wiseman^{∗} and Eytan Domany^{†}

Department of Electronics, Weizmann Institute of science, Rehovot 76100 Israel (October 15, 1993)

Abstract

A cluster Monte Carlo algorithm for the Ashkin-Teller (AT) model is constructed according to the guidelines of a general scheme for such algorithms. Its dynamical behaviour is tested for the square lattice AT model. We perform simulations on the line of critical points along which the exponents vary continuously, and find that critical slowing down is significantly reduced. We find continuous variation of the dynamical exponent z along the line, following the variation of the ratio α/ν, in a manner which satisfies the Li-Sokal bound z_{cluster} ≥ α/ν, that was so far proved only for Potts models.

75.40M, 75.10H, 05.50

Typeset using REVTEX

I. INTRODUCTION

The Ashkin-Teller (AT) model^{1} has been studied extensively, ever since its introduction, by a variety of methods. In two dimensions, in particular, much is known about the phase diagram and critical behavior of the model . Nevertheless there are problems that were not addressed extensively before; the critical behavior of random AT models is one such issue about which not much is known, hope for analytic treatment is slim, and therefore one expects numerical simulations to be the main tool of investigation. With this aim in mind, we set out to develop an efficient Monte Carlo (MC) cluster algorithm for the AT model.

A convenient representation of the generalized AT model (gAT) is in terms of two Ising spin variables, σ_{i} and τ_{i}, placed on every site of a lattice. Denoting by < ij > a pair of nearest neighbor sites, the Hamiltonian is given by

H = − ∑_{<ij>}[K_{σ}σ_{i}σ_{j} + K_{τ}τ_{i}τ_{j} + Lσ_{i}σ_{j}τ_{i}τ_{j}]

Here K_{σ} (K_{τ}) are the strengths of the interactions between neighboring σ (τ) spins, and L is a four-spin coupling. The phase diagram of the ferromagnetic general AT model is known in two dimensions from duality transformations and renormalization group studies^{2,3}. The three dimensional model has been studied as well^{4}. In this paper we are concerned with the Z(4) subspace of the general model, in which K_{σ} = K_{τ} = K. The phase diagram in this subspace is reviewed in Sec. 2. The critical properties of the model in this subspace are of special interest, since it has a line of critical points, along which the exponents vary continuously, and have been determined analytically^{5}, interpolating between Ising and fourstate Potts exponents. For instance, the value of the ratio α/ν varies from 0 at the Ising (L = 0) critical point to α/ν = 1 at the four-state Potts point K = L. This exponent is of special interest to us since it has been proved that for Potts models it serves as lower bound to the dynamic exponent z of cluster algorithms^{6}. Furthermore, the bound seems to be reached in the case of the Ising and Potts models^{6}.

Cluster algorithms^{7} as introduced by Swendsen and Wang^{8}, and extended by Wolff^{9}, are reviewed briefly in Sec. 3. These algorithms give rise to dynamics whose relaxation time, τ_{SW}, is significantly shorter than that of standard, single-spin flip MC methods. This is most important at a critical point, where the relaxation time of a finite system grows with its linear size L according to

τ ∼ L^{z} .

Cluster algorithms have a significantly lower dynamic exponent z_{SW} than that of standard MC. Hence if one is interested in performing extensive simulations of a model such as AT, it is well worth to spend time on developing an appropriate cluster algorithm.

Creating an efficient algorithm can be a challenge. Naive application of the original SW scheme does not work (as explained in Sec. 3 and demonstrated in Sec. 4). We set out to construct an efficient cluster algorithm using the guidelines and methods that were presented by Kandel and Domany^{10} . This general scheme is guaranteed to yield an algorithm that satisfies detailed balance, and once the important excitations of the model have been identified and incorporated, we are guaranteed to get a good algorithm. This general formalism is briefly reviewed, and the resulting algorithm is presented in Sec. 3. Interestingly, the cluster algorithm found this way is identical to what would have been obtained had we used Wolff’s Ising embedding method^{9,11}, as shown in the Appendix. Numerical results are given in Sec. 4; in particular, efficiency of our algorithm is tested by measuring the dynamic exponent along the critical line, which is indeed significantly lower than that of standard MC.

An interesting question we set out to resolve concerns comparison of z_{SW} with α/ν along the AT critical line. We found that the Li-Sokal bound^{6}

z_{cluster} ≥ \frac{α}{ν}

is satisfied, variation of the dynamic exponent follows that of α/ν, and within our numerical accuracy and limitations due to finite size effects, our results indicate that the two are equal.

II. PHASE DIAGRAM OF THE ASHKIN TELLER MODEL

We now review the phase diagram of the square lattice AT model^{4} and some of its critical properties. Figure 1 gives the phase diagram of the model plotted as a function of the parameters X and Z where

Z = exp^{−4K} ; X = exp^{−2(K+L)} .  (1)

At times we use the notation X⃗ = (X, Z) to denote a point in phase space. Full lines of phase transitions separate three phases:

(a) Paramagnetic, labeled as ”P” . The couplings are sufficiently weak so that the system is in a paramagnetic phase in which neither σ nor τ (nor στ) are ordered.

(b) Ferromagnetic phase, labeled as ”F” . The couplings are sufficiently strong so that σ and τ independently order in a ferromagnetic fashion so that ⟨σ⟩ = ±⟨τ⟩. In this phase ⟨στ⟩ is also different from zero and has the same sign as ⟨σ⟩⟨τ⟩.

(c) A phase labeled ⟨στ⟩, in which στ is ordered ferromagneticaly but ⟨σ⟩ = ⟨τ⟩ = 0. This phase arises only for L > K.

On the dashed line Z = X^{2} we have L = 0 ; and thus it is obviously a subspace of two decoupled Ising models, having an Ising transition at the point X_{0}. The dashed line Z = X has L = K, in which case the AT model becomes the 4-state Potts model. The point X_{4} is a 4-state Potts multicritical point.

A marginal operator generates a continuous variation of critical exponents along the line X_{0}X_{4}, isomorphic to the known critical line of the eight-vertex model^{12}. Through an exact duality type transformation this line is mapped onto the critical line of a staggered 8-vertex model and through a relation with the Coulomb gas its critical exponents are known exactly^{5}:

y_{t} = 2 − (2/g_{R})  (2)

where

g_{R} = \frac{8}{π} sin^{−1}[\frac{1}{2} coth(2K)] ,  (3)

and \frac{γ}{ν} = \frac{7}{4} all along the line. Lastly the lines X_{4}B and X_{4}C flow under renormalization to Ising type fixed points.

The exact location of the transition line X_{0}X_{4} can be found through the duality transformation of the AT model^{3}. It is given by the self-dual line Z = 1 − 2X.

III. CLUSTER METHOD FOR THE ASHKIN TELLER MODEL

A. SW cluster method

Cluster algorithms have proved to be a useful method of reducing critical slowing down in MC simulations. For completeness we review here the pioneering cluster algorithm of Swendsen and Wang (SW)^{8} for the Ising model^{13} with the Hamiltonian

H = −J ∑_{<k,j>} σ_{k}σ_{j}  σ = ±1 .  (4)

The SW procedure stochastically identifies clusters of aligned spins, and then flips whole clusters simultaneously. Starting from a given configuration u, SW go over all the bonds, and either ”freeze” or ”delete” them. A bond connecting two neighboring sites k and j, is deleted with probability P_{d}, and frozen with probability P_{f} = 1 − P_{d}, where:

P_{d} = e^{−β(Jσ_{j}σ_{k}+J)} .  (5)

Having gone over all the bonds, all spins which have a path of frozen bonds connecting them, are identified as being in the same cluster. Now the new configuration is generated by flipping every cluster with probability 1/2 . Note that according to (5), only spins of the same sign can be frozen in the same cluster. SW identify correctly the elementary large scale excitations as clusters of aligned spins. This correct identification is essential to the success of a cluster method for any more general model.

In sec. III D we shall present a naive implementation of the original SW scheme to the AT model, and explain why it is not expected to work well. This necessitates search for a different clustering rule; the one we found is based on a general scheme which we now describe.

B. General scheme for cluster methods

A unifying view of all Cluster Algorithms has been given by Kandel and Domany (KD)^{10}, which we now review. The general scheme consists of two steps: given a spin configuration u, the first step consists of stochastically generating a new Hamiltonian H̃. The second step consists of simulating the model with the new Hamiltonian, thereby bringing it to a new configuration u′. To carry out the first step, KD write the Hamiltonian as

H = ∑_{l} V_{l}.  (6)

For example V_{l} can be the energy of a single bond in the nearest neighbor Ising Hamiltonian. Then to each V_{l} they assign one of N possible integers i, in a stochastic manner that depends on the starting configuration u. That is, the probability to assign i to l is written as P_{i}^{l} = P_{i}^{l}(u). This probability is normalized, i.e.

∑_{i} P_{i}^{l}(u) = 1  (7)

for any term l and configuration u. Then they construct a new Hamiltonian:

H̃_{{i}} = ∑_{l} Ṽ_{i(l)}^{l},  (8)

where (for any spin configuration ũ):

Ṽ_{i}^{l}(ũ) = V_{l}(ũ) − \frac{1}{β} ln[P_{i}^{l}(ũ)] + C_{i}^{l}.  (9)

The free parameters C_{i}^{l} are configuration independent.

The second step consists of simulation of the model with any procedure whose transition probability T̃_{{i}}(u → u′) satisfies the detailed balance condition with respect to the new Hamiltonian, i.e.

e^{−βH̃_{{i}}(u)}T̃_{{i}}(u → u′) = e^{−βH̃_{{i}}(u′)}T̃_{{i}}(u′ → u).  (10)

After completing the second step, a new configuration u′ is arrived at, the original Hamiltonian is restored and the process is repeated. Equations ( 9) and ( 10) ensure that the whole procedure satisfies the detailed balance condition, with respect to the original Hamiltonian (for the proof see^{10} ), but ergodicity needs to be proved for each application separately.

We give now two types of modifications to the Hamiltonian. Consider a term V_{l}(u) that can take M distinct energy values E_{i}, i = 1,…,M. The first is the deletion operation, used by SW, which eliminates the interaction V_{l}(ũ) that gets replaced by

Ṽ_{d}^{l}(ũ) = 0,  (11)

for any configuration ũ. To get this, we must have ( see eq. ( 9)):

P_{d}^{l}(u) = e^{β[V_{l}(u)+C_{d}^{l}]}.  (12)

The constant C_{d}^{l} must be chosen so that P_{d}^{l}(u) ≤ 1 for any u.

The second modification is a ‘generalized’ freezing operation which we will later use in our scheme for the Ashkin-Teller model. Its probability is

P_{g}^{l}(u) = \{ \begin{matrix}e^{β[V_{l}(u)+C_{g}^{l}]} & if V_{l}(u) = E_{m} and 1 ≤ m ≤ μ \\ 0 & if V_{l}(u) = E_{m} and μ < m ≤ M \end{matrix} ,  (13)

and the modified interaction according to (9) is

Ṽ_{g}^{l}(ũ) = \{ \begin{matrix}0 & if V_{l}(ũ) = E_{m} and 1 ≤ m ≤ μ \\ ∞ & if V_{l}(ũ) = E_{m} and μ < m ≤ M \end{matrix} .  (14)

Again the constant C_{g}^{l} must be chosen so that P_{g}^{l}(u) ≤ 1 for any u. This type of modification allows free movements between some of the states u of V_{l}(u) which didn’t have the same energy in the original Hamiltonian, and we will use it for the Ashkin Teller model in section III C. Perhaps it can be useful in general, in cases where the V_{l}s can have more than two possible energies. This type of operation has been used for example in^{14}.

A particular case of this operation is obtained if μ = 1, and we set

C_{1}^{l} = \frac{1}{β} ln(p_{1}) − E_{1}.  (15)

Substitution into (13) yields

P_{1}^{l}(u) = \{ \begin{matrix}p_{1} & if V_{l}(u) = E_{1} \\ 0 & otherwise \end{matrix} ,  (16)

and from (14), the modified interaction takes the form :

Ṽ_{1}^{l}(ũ) = \{ \begin{matrix}0 & if V_{l}(ũ) = E_{1} \\ ∞ & otherwise \end{matrix} .  (17)

This operation assigns infinite energy to any configuration ũ for which V_{l}(ũ) ≠ E_{1}. That is, in the ensuing simulation the interaction V_{l} is frozen at energy E_{1}. Freezing of V_{l} is assigned with probability p_{1}, and only when V_{l}(u) = E_{1}. It is easy to see that the SW freezing is precisely of this form.

C. Cluster Method for the Ashkin Teller model

We describe now the cluster algorithm we devised for the gAT model and the considerations that lead us to it. It consists of a freeze-delete scheme which generates non-interacting clusters of aligned σ spins and of aligned τ spins. We will phrase it in terms of the general scheme described in sec. III B.

The first decision one needs to make when coming to design a freeze-delete scheme is the choice of the basic interaction term V_{l}. Our choice is to associate all interactions that reside on an edge < jk > of the lattice to one V_{l} :

V_{l} = −[K_{σ}σ_{k}σ_{j} + K_{τ}τ_{k}τ_{j} + Lσ_{k}τ_{k}σ_{j}τ_{j}].  (18)

Since the model is invariant under any permutation of K_{σ}, K_{τ} and L ( to make this symmetry explicit, define a new Ising spin s_{j} = σ_{j}τ_{j} along with the constraint s_{j}σ_{j}τ_{j} = 1 ), it is possible to choose L ≤ K_{σ}, K_{τ} (the reason for this choice will become clear later). The interaction V_{l} depends on four independent Ising spins that can have 16 states. Every V_{l} can take one of four possible energy values:

\begin{matrix}E_{1} = −K_{σ} − K_{τ} − L \\ E_{2} = K_{σ} − K_{τ} + L \\ E_{3} = −K_{σ} + K_{τ} + L \\ E_{4} = K_{σ} + K_{τ} − L \end{matrix} .  (19)

Every energy is four-fold degenerate; we denote by u_{i} all (four) states for which V_{l} = E_{i}. Let u_{i} represent the state of the spins σ_{j}, σ_{k}, τ_{j}, τ_{k} of a particular pair of n.n. sites < j, k >. Four representative possible states are depicted in figure 2. u_{1} is a ground state in which both the σ and the τ bonds are satisfied. u_{2} and u_{3} are excited states in which one Ising bond is broken and the other is satisfied. Our choice L ≤ K_{σ}, K_{τ} makes u_{4} the highest energy state^{15}, in which both Ising bonds are broken.

The general philosophy of our freeze-delete scheme is as follows. We wish to build clusters of σ spins and of τ spins since we know that these clusters are the basic excitations of the model. Clusters of s = στ spins are not important because L, the bond between them, is the weakest. In order to build clusters, one needs to freeze parallel σ spins to each other and parallel τ spins to each other, and delete the bonds between antiparallel spins. For example in the state u_{2} we wish to freeze the bond between the τ spins and delete the bond between the σ spins. This consideration leads us to include in our scheme two operations. The operation which we will identify as i=2 freezes the bond between τ_{j} and τ_{k} and deletes the bond between σ_{j} and σ_{k} (see fig. 3 ). According to the discussion above, we want to perform this operation with some probability p_{2} ≠ 0 for the state u_{2}, and with probability 0 for the states u_{3} and u_{4}. This operation allows one to move from u_{2} to u_{1}, so in order to maintain detailed balance, we must perform it with some non vanishing probability q_{2} ≠ 0 on u_{1} too. All this is achieved by assigning, in the modified interaction term Ṽ_{2}^{l}, infinite energy to u_{3}, u_{4} and zero energy to u_{1}, u_{2}. This modification is precisely of the form of eq.

(14) so it is of the generalized freezing type described at the end of section III B and hence the probabilities must be (see eq. (13))^{16}:

P_{2}^{l}(u) = \{ \begin{matrix}e^{V_{l}(u)+C_{2}} & if V_{l}(u) = E_{2} or V_{l}(u) = E_{1} \\ 0 & if V_{l}(u) = E_{3} or V_{l}(u) = E_{4} \end{matrix} .  (20)

The choice of C_{2} will set the values of p_{2} and q_{2} :

p_{2} = e^{V_{l}(u_{2})+C_{2}} q_{2} = e^{V_{l}(u_{1})+C_{2}} .  (21)

The operation which we identify as i=3 freezes the bond between σ_{j} and σ_{k} and deletes the bond between τ_{j} and τ_{k}. It follows the same logic as the operation i=2 and its probability is:

P_{3}^{l}(u) = \{ \begin{matrix}e^{V_{l}(u)+C_{3}} & if V_{l}(u) = E_{3} or V_{l}(u) = E_{1} \\ 0 & if V_{l}(u) = E_{2} or V_{l}(u) = E_{4} \end{matrix}  (22)

Again we define the probabilities:

p_{3} = e^{V_{l}(u_{3})+C_{3}} q_{3} = e^{V_{l}(u_{1})+C_{3}} .  (23)

The third operation we wish to perform is a freezing operation of the form defined in section III B equations (16-17). In our case it freezes the bond between σ_{j} and σ_{k} and the bond between τ_{j} and τ_{k}, so it also builds the τ and σ clusters. We’ll denote it by i=1 and its probability is:

P_{1}^{l}(u) = \{ \begin{matrix}p_{1} & if V_{l}(u) = E_{1} \\ 0 & otherwise \end{matrix} .  (24)

The fourth and last operation is the deletion operation defined in section III B equations (11,12) with the probability:

P_{d}^{l}(u) = e^{V_{l}(u)+C_{d}}.  (25)

These four operations define completely our freeze-delete scheme, summarized by fig.

3. In the first column, in each row, one or two representative configurations of u_{i} appear,depicting the state of σ_{j}, τ_{j}, σ_{k} and τ_{k}. In the first row the four modified interactions Ṽ_{i}^{l} are depicted. An upper double line for example, denotes a frozen bond between σ_{j} and σ_{k}, a blank denotes a deleted bond. We have yet to determine the constants C_{2}, C_{3}, p_{1}, C_{d}. In the case of u_{4}, we want to delete the bond between the τ spins and the bond between the σ spins with probability 1, as we never want an unsatisfied bond to be frozen, so we choose C_{d} = −E_{4}, to get P_{d}(u_{4}) = 1. Having chosen C_{d}, all deletion probabilities are set. The rest of the constants are determined by the normalization condition (7), i.e.

\begin{matrix}p_{2} = 1 − P_{d}(u_{2}) \\ p_{3} = 1 − P_{d}(u_{3}) \end{matrix} .  (26)

Consequently, equations (21, 23) determine the constants C_{2}, C_{3} and the probabilities q_{2}, q_{3}. Finally, p_{1} is determined again by the normalization condition

p_{1} = 1 − P_{d}(u_{1}) − q_{2} − q_{3}.  (27)

For completeness we’ll list all the probabilities of our scheme which follow from equations (20–27):

\begin{matrix}P_{d}(u_{i}) = e^{E_{i}−E_{4}} & ∀i \\ p_{i} = 1 − e^{E_{i}−E_{4}} & i = 2, 3 \\ q_{i} = p_{i}e^{E_{1}−E_{i}} & i = 2, 3 \\ p_{1} = 1 − P_{d}(u_{1}) − q_{2} − q_{3} \end{matrix} .  (28)

Checking that all the probabilities of the scheme fulfill the condition 0 < P_{i}(u) < 1 for all i and u is trivial for all probabilities except for p_{1} which is a bit more laborious, but still straightforward.

We’ve described how we generate a new Hamiltonian of non-interacting clusters of σ spins and of τ spins. Now we have to choose some legitimate MC procedure to simulate this Hamiltonian. We chose to do it in a similar fashion to the single cluster algorithm of Wolff^{9}, but the SW version is equally applicable. To be more explicit, we choose a site j of the lattice at random and a random spin, either σ_{j} or τ_{j} ( The choice of σ or τ was done with probability 1/2 which is sensible in the case K_{σ} = K_{τ}, but in general any probability is acceptable, and an optimal choice can be made to minimize the autocorrelation time τ). This spin will belong to the cluster we will flip, and the cluster will contain only τ (σ) spins if we initially chose τ_{j} (σ_{j}). We perform the freeze-delete operations only on bonds belonging to the surface of the cluster. For example suppose at some stage τ_{j} was joined to the cluster, and suppose that the term V_{j,k} was modified according to i=3, then τ_{k} will not be joined to the cluster and we need not perform freeze-delete operations on the other three bonds connecting the site k. If it is modified according to i=1 or i=2 then τ_{k} is joined to the cluster. Since our cluster algorithm fits into the general scheme, we do not need to prove detailed balance, while ergodicity is ensured by the deletion operation.

We’ve explained the reasoning behind our algorithm. We believe that the clusters we build of parallel σ spins and of parallel τ spins are the basic excitations of the model. We therefore believe that it will be efficient. It is encouraging to notice that in the decoupled Ising subspace and the 4-state Potts subspace our freeze/delete scheme is identical with SW’s freeze/delete operations for these models. In the appendix we show that our algorithm is equivalent to the idea of embedding into the AT model an Ising model and simulating it by Wolff’s single cluster procedure. In sec. IV we present numerical evidence for our algorithm’s efficiency.

D. ‘Naive’ SW option

We find it illuminating to compare our algorithm to a cluster algorithm which one could regard as the naive generalization of the SW method to the AT model. Such an algorithm would define V_{l} in the same manner as our scheme does ( see eq. 18). For each u_{i}, the bonds between the two neighboring sites get either deleted ( with our P_{d}(u_{i}) ) or frozen with P_{f}(u_{i}) = 1 − P_{d}(u_{i}). Fig. 4 can clarify how this scheme fits into the general one of sec. III B and how it compares with our scheme. Since it fits into the general scheme we do not need to prove detailed balance, while the deletion operation ensures its ergodicity except for T = 0.

This scheme also generates clusters of σ spins and clusters of τ spins, but with the ‘naive’ scheme, antiparallel Ising spins can be found in the same cluster, so this scheme does not identify the elementary excitations of the model. Moreover the σ clusters’ structure is forced to match the τ clusters ( and vice versa). This is totally unphysical since in practice it is energetically favorable for the two cluster structures \underline{not} to match ( for L < K_{σ}, K_{τ} ). At T = 0 this scheme could freeze the whole lattice into a single cluster even when it is not in the ground state. This is in contrast to the ability of our scheme to relax excitations on any scale even at T = 0 ( in other words, At T = 0 our algorithm will bring any initial configuration to the ground-state in a short time). At a finite temperature the ‘naive’ SW scheme will produce clusters that are too large. In sec. IV A we list results of simulations using the single cluster (1C) version of this ‘naive’ SW scheme, which show that it is indeed much less efficient than our algorithm.

IV. SIMULATIONS OF THE AT CRITICAL LINE

In addition to checking the efficiency of our method at the AT critical line, we also wanted to check a prediction made by Li and Sokal^{6}, who have proved a rigorous lower bound

z_{SW} ≥ α/ν  (29)

for the dynamical critical exponent of the SW algorithm for the ferromagnetic q-state Potts model. More precisely, the Li-Sokal bound is τ_{exp/int} ≥ const × C_{h}, from which eq. (29) follows (for definitions of τ_{exp/int} and z_{exp/int} see^{17}). This bound relates dynamics to the static properties of a model and is thus of great importance. Their proof is for both z_{int} and z_{exp}. We wanted to check whether this bound holds all along the AT critical line, which connects the decoupled Ising critical point with \frac{α}{ν} = 0 at one end to the 4-state Potts critical point with \frac{α}{ν} = 1 at the other end. At both of these points our algorithm is identical to Wolff’s 1C version of the SW method for Potts models. The Li-Sokal bound was proved for the SW dynamics, but perhaps it is valid for Wolff’s 1C dynamics as well (at least for d = 2).

We are unaware of any rigorous proof for that, but our results for z for q = 2 and 4 seem to indicate z_{1C} = z_{SW} and are in accordance with previous results for the Ising critical point (see^{8} and^{18}). Besides estimating the dynamical exponents, we’ve estimated, using finite size scaling, the critical exponents \frac{α}{ν} and \frac{γ}{ν} along the AT critical line. As stated in sec. II the AT critical line is given by

Z = 1 − 2X.  (30)

Our measurements were done at the decoupled Ising critical point X⃗_{0} and at the 4-state Potts critical point X⃗_{4} (see figure 1 and exact definition in sec. II. Three additional measurements were carried out at three intermediate equidistant points in the X − Z plane, on the AT critical line. So a total of 5 measurements were done at the points

X⃗_{i} = X⃗_{0} + \frac{i}{4}(X⃗_{4} − X⃗_{0}) i = 0 … 4.  (31)

All the points X⃗_{i} are marked in fig. 1 . We simulated lattices with periodic boundary conditions of up to size 128 × 128, and up to 5 × 10^{5} clusters were flipped for each lattice size.

We calculated the energy^{19} per site:

E = ⟨E⟩ = −\frac{1}{L^{2}} ⟨ ∑_{<k,j>}[Kσ_{k}σ_{j} + Kτ_{k}τ_{j} + Lσ_{k}τ_{k}σ_{j}τ_{j}]⟩  (32)

where L is the linear lattice size, and the angular brackets denote the usual thermal MC average. The specific heat per site follows from the energy fluctuations

C = \frac{C_{h}}{k_{B}} = L^{2}(⟨E^{2}⟩ − ⟨E⟩^{2}).  (33)

We calculated the magnetic susceptibility of the σ spins defined as^{19}

χ = L^{2}⟨M^{2}⟩  M = \frac{1}{L^{2}} ∑_{k} σ_{k}  (34)

and also measured c, the size of the cluster flipped at each step, and calculated ⟨c⟩. There is a connection between the size of the clusters and the susceptibility^{20} which is common for algorithms which build non-interacting clusters of spins (if all spins in a cluster have the same value). For 1C algorithms it has the simple form^{9}

χ = ⟨c⟩.  (35)

To get the dynamic properties we calculated the time dependent autocorrelation functions ϕ_{E}(t) and ϕ_{χ}(t) defined as

ϕ_{A}(t) = \frac{⟨A(0)A(t)⟩ − ⟨A⟩^{2}}{⟨A^{2}⟩ − ⟨A⟩^{2}} ,  (36)

where A stands for E or χ. A typical plot of ϕ_{χ}(t), measured from the X_{2} model at L = 128 is presented in fig. 6.

A. MC results and discussion

1. Susceptibility and Specific heat

According to finite-size scaling theory^{21} one expects:

χ ∼ L^{\frac{γ}{ν}} C ∼ L^{\frac{α}{ν}}  (37)

for large enough L. Fitting our measurements of χ to eq. (37) for lattice sizes L ≥ 16 fits the exact universal value \frac{γ}{ν} = \frac{7}{4} within errors (see table I) as expected. Our measurements confirm the equality (35), where ⟨χ⟩ and ⟨c⟩ have an error of the same magnitude. Plots of log C vs. log L for the five models can be seen in fig. 5. From linear fits to the log − log plots, estimates for \frac{α}{ν} are obtained which do not agree with the exact known values. For comparison with the exact values see table I and fig. 13. The differences may be due to finite size effects and corrections to scaling . An exception to this mismatch is the decoupled Ising point for which the value of \frac{α}{ν} = 0 fits nicely according to the semi-log plot in fig. 8 (for completeness we also quote in table I an estimate for \frac{α}{ν} at X_{0} from the log − log fit of fig. 5). The slope of the C curves at X_{1} and X_{2} does tend to a lower value with increasing lattice size. The 4-state Potts model is known to have a logarithmic correction to scaling: C ∼ L/ log^{3/2} L^{6}; in any case, our result \frac{α}{ν} = .747(3) agrees with previous MC results obtained from lattices of sizes up to L = 256^{6}.

2. τ_{exp} and τ_{int}

In order to check the efficiency of our algorithm and to see whether the Li-Sokal bound holds, we measured τ_{exp} and τ_{int} for both the energy E and the susceptibility χ. This was done by the following procedure. ϕ_{E/χ}(t) was plotted on a semi-log plot, the unit of time being a single cluster flip. τ_{exp} was then consistently estimated by the slope of ln ϕ(t) in a window from τ_{exp} to 3τ_{exp}. A typical example for ϕ(t) and the linear fit to extract τ_{exp} can be seen in figure 6. To calculate τ_{int}, we integrated numerically ϕ(t) in the interval from t = 0 up to t ≈ 1.5τ_{exp} and estimated the tail (t > 1.5τ_{exp} ) of ϕ(t) by the exponential fit. The error of both τ’s was estimated from repeated experiments when the statistics were large enough and otherwise from subjective estimates from the fluctuations in τ_{exp} in the window τ_{exp} − 3τ_{exp}. Each MC step, or cluster flip, involves on the average the flipping of ⟨c⟩ spins. Since the natural unit of time is a sweep over all spins of the lattice, we multiplied τ by ⟨c⟩L^{−2}.

In figures 7 through 12 we plot measured values of τ_{exp} and τ_{int} for E and χ as a function of the linear lattice size L. In figures 7–10 we display the measurements of τ_{χ} at different models X_{i} separately, while in figures 11 and 12 we display τ_{exp,E} and τ_{int,E} for the five models jointly. Fitting our results to the forms

τ_{exp} ∼ L^{z_{exp}} ; τ_{int} ∼ L^{z_{int}} ,  (38)

for large enough L, yielded the dynamical critical exponents listed in table I. Figures 7,8 are both from the decoupled Ising model point X_{0}. In table I we list values of z obtained from the power law fits of fig. 7. None the less, comparing these fits to the log L fits of fig. 8, we feel that the latter ones are better and that the true value is z = 0 for all four autocorrelation times τ. At X_{0} ( and all along the Z = X^{2} or L = 0 line ) our algorithm is identical with the Wolff 1C algorithm for the Ising model. Our results for z_{int} from the fits in fig. 7 are consistent with those of Wolff^{18}. A log L behaviour has been measured for τ_{exp,M} and τ_{int,M} for SW dynamics in^{22} .

At the point X_{4} , the 4-state Potts critical point (and all along the X = Z line), our algorithm is again identical with the Wolff 1C algorithm for the 4-state Potts model. Our result for z_{int,E}, obtained from the fit in figure 12 is consistent with the result of ref.^{6} for SW dynamics. This result can be added to the accumulating evidence indicating: z_{SW} = z_{W olff} for d = 2.

The continuous variation of z along the AT critical line can be explicitly seen in figures 11 and 12 where we plot τ_{exp,E} and τ_{int,E} for the 5 points X_{i}. In fig. 12 we also plot results from simulations at X_{2} using the Metropolis method and the ‘naive’ SW method described in sec. III D. For the Metropolis method, we measured a value of z_{int,E} = 1.64(8) , which should be compared with z_{int,E} = .542(8) using our method. The ‘naive’ SW method almost froze the whole lattice, into a single cluster, the size of which was almost independent of the lattice size L. For example, the average cluster size for L = 32 was \frac{⟨c⟩}{L^{2}} = .868(2) , as compared with \frac{⟨c⟩}{L^{2}} = .481(1) using our algorithm. We found determination of τ for the ‘naive’ SW method of lattice sizes L > 32 so time consuming that it was impractical. The advantage of our method is clearly demonstrated. The three methods yielded the same results for the static observables (within errors).

The main question we wish to answer is whether the Li-Sokal bound is fulfilled and whether it is sharp. In fig. 13 we compare, for the 5 points X_{i} i = 0 … 4, the exact values of \frac{α}{ν} and our estimated values of \frac{α}{ν}, z_{int,χ} and z_{exp,χ}. We see that the rise of our estimated values of z follows that of \frac{α}{ν}, from the decoupled Ising point X_{0} to the 4-state Potts point X_{4}. Except for the point X_{4} the Li-Sokal bound is fulfilled with respect to (\frac{α}{ν})_{exact}. The anomalously low values for z at X_{4} are probably caused by the multiplicative logarithmic correction for C described in the discussion of the specific heat results. This explanation has been suggested by Li and Sokal^{6} to explain the low value of z_{SW} for this model. From table I we see that z ≥ (\frac{α}{ν})_{estimate} with only one exception at the point X_{1} (see fig. 9). The fact that z_{int,χ} < (\frac{α}{ν})_{estimate} could be due to a difference in the finite size corrections or in the corrections to scaling. τ_{int,χ} is probably less influenced by one or both of these two factors, in comparison with the large difference between (\frac{α}{ν})_{estimate} and (\frac{α}{ν})_{exact}. At the point X_{3}, where the difference between (\frac{α}{ν})_{estimate} and (\frac{α}{ν})_{exact} is the smallest, no anomalies in the Li-Sokal bound occur (see fig. 10).

We conclude that the Li-Sokal bound is fulfilled in a moderately sharp manner. The smallest z (which is z_{int,χ}) fulfills (\frac{α}{ν})_{estimate} ≤ z_{int,χ} ≤ (\frac{α}{ν})_{estimate} + 0.1 (not including the anomalies of the models X_{4} and X_{1} discussed above). Note that we are comparing estimated values of z with \underline{estimated} values of \frac{α}{ν} . This is the correct comparison to make, assuming that τ and C have similar finite size corrections and similar corrections to scaling.

V. SUMMARY AND DISCUSSION

The correct identification of the basic excitations of the model lead, along with the guidelines of a general scheme for cluster algorithms, to the construction of a cluster algorithm for the AT model. The algorithm was shown to be identical with the one obtained by embedding Ising spins into the AT model. Our algorithm is suitable for spatially varying coupling constants under the restriction L^{k,j} ≤ K_{σ}^{k,j}, K_{τ}^{k,j} everywhere on the lattice. We are currently carrying out intensive simulations of a random-bond version of the AT model.

The dynamical behaviour of the cluster algorithm was examined on the AT critical line. Critical slowing down of the algorithm ( 0 ≤ z < 1. ) was found to be significantly smaller than that of the standard Metropolis method. A continuous variation of the dynamical exponent z along the AT critical line was seen, along with continuous variation of the static exponent \frac{α}{ν}. The Li-Sokal bound z ≥ \frac{α}{ν}, that was proved only for q-state Potts models with SW dynamics, is satisfied for the AT model with single cluster dynamics. The bound is moderately sharp. Another new result is z_{W olff} = z_{SW} for the 4-state Potts model, which can be added to the accumulating results indicating that in two dimensions z_{W olff} = z_{SW} .

ACKNOWLEDGMENTS

We would like to thank Robert Swendsen and Ulli Wolff for helpful discussions. This research has been supported by the US-Israel Bi-national Science Foundation (BSF).

APPENDIX: ISING EMBEDDING

Our algorithm can also be seen from a totally different point of view; as an example of an embedding algorithm^{9,11}. The main idea is to embed into the AT model an Ising model of space dependent couplings J_{jk} and simulate it using the SW or Wolff procedure for the Ising model. To be explicit, consider the Hamiltonian (I), and take the τ variables as fixed, so we can write

H = H_{1} + H_{2} = − ∑_{<k,j>}(K_{σ} + Lτ_{k}τ_{j})σ_{k}σ_{j} − ∑_{<k,j>} K_{τ}τ_{k}τ_{j}.  (A1)

H_{2} represented by the second sum is a constant, and remembering that we chose L < K_{σ}, K_{τ}, H_{1} is a ferromagnetic Ising model in the σ variables with couplings J_{jk} = K_{σ} + Lτ_{k}τ_{j}. Simulating this Ising model with any procedure that will maintain detailed balance with respect to H_{1} and will not change the value of H_{2}, will also maintain detailed balance with respect to H. So we can use for example the SW or Wolff procedure for the Ising model, explained in sec. III B. This by itself will maintain detailed balance but will not be ergodic since the τ variables will not be updated. Obviously to update the τ variables the same process should be repeated , holding the σ variables fixed and simulating an Ising Hamiltonian with the τ variables. To summarize, a possible procedure would go as follows: Choose at random whether to embed into the AT Hamiltonian an Ising Hamiltonian in the σ spins or in the τ spins. Pick a random site in the lattice, grow a cluster of σ or τ spins using the Wolff (1C) procedure using the Ising Hamiltonian and flip it. As we will now show, this process ( we’ll call it for shortness IE – Ising embedding) is exactly equivalent to the AT algorithm (ATA) we suggested in section III C .

A first argument which is really sufficient goes as follows. Suppose we perform the IE freeze-delete scheme for the σ spins for the whole lattice, viewing the τ spins as constant. We could say we’ve generated a new Hamiltonian of noninteracting clusters of σ spins and which assigns infinite energy to any configuration which differs from the current configuration of τ spins. Now suppose we perform the ATA freeze-delete scheme for the whole lattice. Then we get a new Hamiltonian of noninteracting clusters of σ spins and of τ spins, but if at this stage we decide to flip only σ spins, then in practice we’ve assigned infinite energy to any configuration which differs from the current configuration of τ spins. In practice we flip only one cluster of the σ spins, so what we do is identical to the Wolff (1C) version of the IE. Since both processes maintain detailed balance using the same new Hamiltonian they must do so using the same probabilities (the general scheme in sec. III B shows a one to one correspondence between probabilities and the new Hamiltonian.), which completes our argument that the two methods are actually identical.

One can, of course, check that the probabilities of the two procedures are the same. For example, denote the probability to delete a bond between σ spins in ATA as P_{d,σ}^{AT A}(u). for u_{1} for example, From fig. 3:

P_{d,σ}^{AT A}(u_{1}) = P_{d}(u_{1}) + q_{2} = e^{E_{1}−E_{2}} = e^{−2(K_{σ}+L)}.  (A2)

Now with the Ising Embedding algorithm, according to (5):

P_{d,σ}^{IE}(u) = e^{−(K_{σ}+Lτ_{k}τ_{j})(σ_{k}σ_{j}+1)},  (A3)

so P_{d,σ}^{IE}(u_{1}) = P_{d,σ}^{AT A}(u_{1}). This check can be carried out for all u_{i}.

REFERENCES

∗ email address: shai@elect1.weizmann.ac.il † email address: fedomany@WEIZMANN.weizmann.ac.il 1 J. Ashkin and E. Teller, Phys. Rev.64, 178 (1943). 2 F.Y. Wu and K.Y. Lin, J.Phys. C7, L181 (1974). 3 E. Domany and E.K. Riedel, Phys. Rev. B19, 5817 (1979). 4 R.V. Ditzian, J.R. Banavar, G.S. Grest and L.P. Kadanoff, Phys. Rev. B22, 2542, (1980). 5 B. Nienhuis , in Phase Transitions and Critical Phenomena, edited by C. Domb and J.L. Lebowitz(Academic, New York, 1987), vol. 11. 6 X.J. Li and A.D. Sokal, Phys. Rev. Lett.63, 827 (1989). 7 A.D. Sokal, Nuc. Phys. B (Proc. Suppl.),20,55, (1991) 8 R.H. Swendsen and J.S. Wang, Phys. Rev. Lett.58, 86 (1987). 9 U. Wolff, Phys. Rev. Lett.62, 361 (1989). 10 D. Kandel and E. Domany, Phys. Rev. B43, 8539 (1991). 11 S. Caracciolo, R.G. Edwards, A. Pelissetto, and A.D. Sokal, Nuc. Phys. B (Proc.

Suppl.)20, 72 (1991). 12 R.J. Baxter, Ann. Phys. (N.Y.)70, 193 (1972). 13 Actually the algorithm is for the q-state Potts model. 14 H. G. Evertz, G. Lana, and M. Marcu, Phys. Rev. Lett.70, 875 (1993). 15 If L = K_{σ} or L = K_{τ} then E_{4} = E_{3} or E_{4} = E_{2}. We do not discuss here the behaviour of

our algorithm in such a case. 16 From now on we’ll absorb β into the coupling constants of the Hamiltonian. We’ll also

omit the l dependence of all constants of the algorithm since we are considering, for the time being, constant K_{σ}, K_{τ} and L. Generalization to spatially varying couplings is trivial as long as K_{σ}, K_{τ} ≥ L everywhere.

17 A.D. Sokal Nuc. Phys. B (Proc. Suppl.),20:55, (1991). 18 U. Wolff, Phys.Lett. B,228:379, (1989). 19 A factor of k_{B}T has been omitted, but is irrelevant since all measurements are at a single

temperature. 20 U. Wolff Nucl. Phys. B,300:501,(1988). 21 M.N. Barber, in Phase Transitions and Critical Phenomena, edited by C. Domb and J.L.

Lebowitz, vol. 8. 22 D.W. Heermann and A.N. Burkitt, Physica A162, 210 (1990).

FIGURES

FIG. 1. Phase diagram of the Ashkin Teller (Z(4) ) model.

FIG. 2. States of spins at a pair of nearest-neighbor sites j, k.Each state u_{i} represents one out of four states with the same energy E_{i}.

FIG. 3. Freeze/delete scheme for the AT model. The four states of two neighboring sites are denoted in the left column; Four freeze/delete operations are denoted in the top row, and the respective probabilities are denoted in the table. A frozen bond is denoted by a double thick line, while a deleted bond is denoted by a blank.

FIG. 4. ‘Naive’ SW freeze/delete scheme for the AT model. The four states of two neighboring sites are denoted in the left column; Four freeze/delete operations are denoted in the top row, and the respective probabilities are denoted in the table. A frozen bond is denoted by two thick lines, while a deleted bond is denoted by a blank.

FIG. 5. A log-log plot of specific heat vs. linear lattice size at the five critical models X_{0} … X_{4}. The solid lines are the linear fits. the critical exponents \frac{α}{ν} are given in table 1.

FIG. 6. The ln of the time auto-correlation function of the susceptibility ϕ_{χ}(t), measured for the X_{2} model at L=128. The unit of time is a single MC step or a single cluster flip. The vertical lines are the error bars. The exponential fit is also shown.

FIG. 7. Log-log plots of auto-correlation times τ_{χ} and the specific heat C as a function of linear lattice size L, at the decoupled Ising point X_{0}. The values of z listed in table 1 are the slopes of the fits.

FIG. 8. Semi-log plots of auto-correlation times τ_{χ} and the specific heat C as a function of linear lattice size L, at the decoupled Ising point X_{0}. This fit seems to be better than the log-log fit, yielding values of z = \frac{α}{ν} = 0.

FIG. 9. Log-log plots of auto-correlation times τ_{χ} and the specific heat C as a function of linear lattice size L, at the point X_{1}. The values of z listed in table 1 are the slopes of the fits.

FIG. 10. Log-log plots of auto-correlation times τ_{χ} and the specific heat C as a function of linear lattice size L, at the point X_{3}. The values of z listed in table 1 are the slopes of the fits.

FIG. 11. A log-log plot of τ_{exp,E} vs. linear lattice size at the five critical models X_{0} … X_{4}. The solid lines are the linear fits. the critical exponents z_{exp,E} are given in table 1.

FIG. 12. Log-log plots of τ_{int,E} as a function of linear lattice size L, at all 5 points X_{i} i = 0 … 4. The continuous variation of the slope z along the AT critical line is easily seen. Metropolis results and ‘naive’ SW results for X_{2} are also plotted.

FIG. 13. Comparison of z and \frac{α}{ν} at the 5 points X_{i} i = 0 … 4. The line denotes the exact value of \frac{α}{ν}, while our estimated values for the five models are denoted by full circles. values of z_{int,χ}/z_{exp,χ} are denoted by empty circles/ crosses.

TABLES

TABLE I. Results from the AT critical line. The errors in parentheses are only the statistical errors of the fit in our fitting interval. They do not include systematic errors stemming from finite size effects and corrections to scaling.

\frac{γ}{ν}  \frac{α}{ν}  z_{int,E}  z_{exp,E}  z_{int,χ}  z_{exp,χ}  \frac{α}{ν}_{exact}
X_{0}  1.751(4)  .23(1) / 0  .26(3)  .23(2)  .13(1)  .267(4)  0
X_{1}  1.751(1)  .38(2)  .396(5)  .40(3)  .273(3)  .47(2)  .2096
X_{2}  1.752(3)  .542(8)  .61(1)  .68(3)  .532(5)  .62(1)  .4182
X_{3}  1.763(6)  .655(13)  .719(2)  .74(2)  .717(8)  .72(3)  .6383
X_{4}  1.752(4)  .747(3)  .92(1)  .99(5)  .99(3)  .94(3)  1