# Dynamic and static properties of the invaded cluster algorithm

## I Introduction

Monte Carlo (MC) methods used to simulate classical spin systems, such as Potts models, fall primarily into two broad classes: local-update algorithms and cluster algorithms. Algorithms with local update rules, such as the Metropolis algorithm, provide an efficient means of simulating these spin systems in non-critical regions. Near a second-order phase transition, however, where long-range correlations are present, relaxation times increase rapidly with system size. This phenomenon, known as critical slowing down, may be characterized by a dynamic exponent $`z`$ according to $`\tau L^z`$, where $`\tau `$ is the autocorrelation time measured at criticality (roughly, the time necessary to generate a statistically independent configuration) and $`L`$ is the system size. Local-update algorithms typically have values of $`z`$ slightly greater than 2 and, therefore, are impractical for simulating large systems near a critical point.

Cluster algorithms, on the other hand, such as the Swendsen-Wang (SW)  algorithm, employ non-local update moves, flipping clusters of spins of linear extent comparable to the correlation length. This technique significantly reduces critical slowing down, and thus makes cluster algorithms preferable for simulating spin systems near a critical phase transition. Recent numerical estimates of the SW dynamic exponent for two-dimensional ferromagnetic $`q`$-state Potts models are $`z^{\mathrm{SW}}.25`$  for the Ising ($`q=2`$) case and $`z^{\mathrm{SW}}.52`$  for $`q=3`$.

Although there currently exists no theoretical means by which the dynamic exponent of a SW-type algorithm may be calculated, there is a rigorous lower bound. The Li-Sokal  bound, as it has come to be known, states that $`z^{\mathrm{SW}}\alpha /\nu `$ for $`q`$-state Potts models, where $`\alpha `$ and $`\nu `$ are the usual static critical exponents for the specific heat and correlation length, respectively. We note that the numerical values given above are consistent with this bound, since, in two dimensions, $`\alpha /\nu =0`$ $`(\mathrm{log})`$ for $`q=2`$, and $`\alpha /\nu =2/5`$ for $`q=3`$. For the remainder of this paper we will continue to focus on the Ising and 3-state Potts models in two-dimensions, since these are the two most carefully studied cases.

We now turn to the invaded cluster (IC)  algorithm, a recent approach based on invasion percolation, for simulating equilibrium critical points. This algorithm has the unique property that it “self-organizes” to the critical point. Therefore, no a priori knowledge of the critical temperature is required; instead, $`T_\mathrm{c}`$ is an output of the algorithm. In addition, due to an intrinsic negative-feedback mechanism, the IC algorithm equilibrates very quickly in the sense that thermodynamic quantities are measured to be near their equilibrium values within a few MC steps (after starting, say, from a completely ordered state).

Initial studies  seemed to indicate that the IC algorithm suffers no critical slowing down for the Ising model. For both $`d=2`$ and $`d=3`$, the integrated autocorrelation time $`\tau _\epsilon `$ was observed to decrease with $`L`$, while $`\tau _m`$ remained constant (within error bars), where $`\epsilon `$ is the energy per spin, and $`m`$ is the fraction of spins in the largest cluster. (We omit the usual “int” subscript on $`\tau `$, since we deal almost exclusively with integrated as opposed to exponential autocorrelation times.) The decrease of $`\tau _\epsilon `$ with $`L`$ was also observed for $`3`$\- and $`4`$-state Potts models in two dimensions, but in these cases critical slowing was evident in the behavior of $`\tau _m`$. In Ref.  the dynamic exponents were estimated to be $`z_m.28`$ for $`q=3`$ and $`z_m.63`$ for $`q=4`$, each of which is less than the Li-Sokal bound on $`z^{\mathrm{SW}}`$ for its respective value of $`q`$.

In this paper, we repeat these studies for $`q=2`$ and $`q=3`$ in two dimensions using larger lattice sizes and an improved method  of estimating $`\tau `$. We also investigate the $`L`$-dependence of the “specific-heat-like” quantity $`c(L)L^d\mathrm{var}(\epsilon )`$. (In the canonical ensemble $`c`$ is the specific heat, but the same is not true for the IC ensemble.) In addition we measure $`\tau _\epsilon (l,L)`$, the integrated autocorrelation time for the energy per spin $`\epsilon (l,L)`$ measured on an $`l\times l`$ sub-lattice of the whole $`L\times L`$ lattice, as well as $`c(l,L)l^d\mathrm{var}(\epsilon (l,L))`$. While the motivations for these experiments are discussed in more detail in Sec. III, the central idea is that we expect the negative feedback to diminish for length scales $`l<L`$, leading to sub-system behavior that differs from that of the whole system. We argue in Sec. III that, for $`l`$ sufficiently less than $`L`$, we should observe a crossover to SW behavior for observables measured on a sub-lattice of size $`l`$. Upon investigating the crossover region, we find that $`c(l,L)`$ has a simple scaling form and give an estimate of the length scale at which $`c(l,L)`$ crosses over to its SW analog, namely, the specific heat for an $`l\times l`$ system.

The results are less conclusive in the case of the dynamic variable $`\tau _\epsilon (l,L)`$, but it appears that crossover to SW behavior does occur and that the crossover length for $`\tau _\epsilon (l,L)`$ differs from that for $`c(l,L)`$ and is likely the same for $`q=2`$ and $`q=3`$. In addition, the crossover phenomenon leads to a maximum in $`\tau _\epsilon (l,L)`$ as $`l`$ is varied for a given $`L`$. We argue in Sec. III that the dynamic exponent $`z^{\mathrm{IC}}`$ of the IC algorithm is appropriately defined by $`\tau _{\epsilon ,\mathrm{max}}L^{z^{\mathrm{IC}}}`$. For the 3-state Potts model we give a numerical estimate for $`z^{\mathrm{IC}}`$ that is smaller than a recent estimate  of $`z^{\mathrm{SW}}`$ and also less than the Li-Sokal bound on $`z^{\mathrm{SW}}`$. For the Ising model the results are less conclusive, but it appears that $`z^{\mathrm{IC}}<.21`$ and possibly that $`\tau _{\epsilon ,\mathrm{max}}\mathrm{log}L`$ so that $`z^{\mathrm{IC}}=0`$—similar to the state of affairs for the SW algorithm .

The remainder of this paper is organized as follows. In Sec. II we provide some background on the invaded cluster algorithm and discuss results of previous IC simulations. In Sec. III we discuss the crossover phenomenon in greater detail and propose a scaling form to relate observables measured on sub-lattices using the IC algorithm to corresponding quantities for the SW algorithm. We describe our IC simulations of the Ising and 3-state Potts models in Sec. IV and discuss the results in Sec. V. Sec. VI contains our conclusions.

## II Invaded cluster algorithm for critical Potts models

In order to understand how the IC algorithm works, we first review the SW algorithm for Potts models. The ferromagnetic $`q`$-state Potts hamiltonian is

$$H=\underset{i,j}{}\delta _{\sigma _i,\sigma _j},$$
(1)   
where $`\sigma _i\{0,1,\mathrm{},q1\}`$ and the sum is over nearest-neighbor spin pairs. (Note that the Ising model is just the special case $`q=2`$.)

Given an initial spin configuration, the SW algorithm proceeds as follows: First, satisfied bonds are occupied with probability $`p=1e^\beta `$, with $`\beta =1/T`$, where a bond joining two spins $`i`$ and $`j`$ is defined to be satisfied if and only if $`\sigma _i=\sigma _j`$. Unsatisfied bonds are never occupied. Next, clusters of spins connected by occupied bonds are identified, and each cluster is “flipped”, i.e., independently and uniformly assigned a new random spin value from $`\{0,1,\mathrm{},q1\}`$. (Note that a cluster can consist of a single spin.) Finally, after statistics have been collected, occupied bonds are erased and the whole process is repeated. It can easily be shown that the SW algorithm satisfies detailed balance for the canonical ensemble.

The IC algorithm uses invasion percolation to generate the spin clusters to be flipped. Given an initial spin configuration, the first step is to assign a random order to the bonds of the lattice. The bonds are then examined, and satisfied bonds occupied one at a time in this order. If a bond joining two clusters is occupied, they are combined into one. Cluster growth continues until some stopping condition is fulfilled. In this paper, we consider the topological spanning condition, which dictates that growth be stopped as soon as some cluster winds around the system in one of the $`d`$ directions. As soon as spanning is detected, clusters (including the spanning cluster) are flipped exactly as in the SW algorithm, statistics are collected, bonds erased, and the process repeated.

To understand why the IC algorithm self-organizes to the critical point, we define $`f`$ to be the ratio of the number of occupied bonds to the number of satisfied bonds when some cluster first spans the system. It has been argued  that, as the system size $`L`$ approaches infinity, the distribution of $`f`$ approaches a delta function at $`p_\mathrm{c}1e^{\beta _\mathrm{c}}`$, where $`\beta _\mathrm{c}`$ is the inverse critical temperature. Though not a rigorous proof, the argument proceeds as follows. First, we note that $`p_\mathrm{c}`$ is the threshold for percolation on the satisfied bonds of a critical spin configuration . Thus, given a spin configuration that is typical of the critical point, the fraction $`f`$ of satisfied bonds that must be occupied to achieve spanning is close to $`p_\mathrm{c}`$. Second, we observe that each iteration of the IC algorithm is identical to an iteration of the SW algorithm with $`p=f`$. Therefore, performing an iteration of the IC algorithm on a critical spin configuration is equivalent to performing an iteration of the SW algorithm with $`pp_\mathrm{c}`$, and thus the system will remain near the critical point.

If, instead, the system is started in the low-temperature phase, the number of satisfied bonds will be larger than is typical of $`T_\mathrm{c}`$. Therefore, a smaller fraction $`f`$ will need to be occupied to achieve spanning. In this case, an IC iteration is equivalent to a SW iteration with $`p<p_\mathrm{c}`$, i.e., $`T>T_\mathrm{c}`$, and therefore the system is pushed toward $`T_\mathrm{c}`$ from below. Similarly, if the system is started in the high-temperature phase, it is pushed toward $`T_\mathrm{c}`$ from above. Thus, in summary, a system in a non-critical state is pushed toward criticality, while a system in a critical state remains near criticality with $`T`$ fluctuating about $`T_\mathrm{c}`$.

Because of this negative-feedback mechanism, the IC algorithm self-organizes to the critical point with no a priori knowledge of $`T_\mathrm{c}`$. Instead, $`T_\mathrm{c}`$ is obtained as an output of the algorithm, via the relation $`f=1e^{1/T_\mathrm{c}}`$. For example, results of IC simulations for the $`2d`$ Ising model yield an estimate of $`T_\mathrm{c}1.1355`$ when extrapolated to $`L=\mathrm{}`$, as compared to the exact result for an infinite system, $`T_\mathrm{c}=1.1346\mathrm{}`$.

Since every configuration generated by the IC algorithm must contain a spanning cluster, it is clear that the algorithm does not sample the canonical ensemble for a system of finite volume. We refer to the stationary distribution sampled by the IC algorithm as the IC ensemble. If we assume that as $`L\mathrm{}`$ the distribution of $`f`$ approaches a delta function at $`p_\mathrm{c}`$ and that the volume fraction of the spanning cluster goes to zero, local observables such as internal energy and magnetization will approach their infinite-volume critical values in this limit. Simulation results  support this hypothesis. For example, results of IC simulations for the $`2d`$ Ising model yield an estimate for the energy per spin of $`\epsilon _\mathrm{c}1.706`$  when extrapolated to $`L=\mathrm{}`$, as compared to the exact result for an infinite system, $`\epsilon _\mathrm{c}=1.7071\mathrm{}`$.

Some finite-volume fluctuations in the IC ensemble, however, are very different from those in the canonical ensemble. For example, in the canonical ensemble the quantity $`c(L)L^d\mathrm{var}(\epsilon )`$ is the specific heat which diverges as $`L^{\alpha /\nu }`$ at the critical point—a logarithmic divergence for the Ising model in two dimensions. In the IC ensemble, however, $`c(L)`$ is observed  to increase roughly linearly with $`L`$ for the $`2d`$ Ising model. These differences can be traced to fluctuations in the effective temperature (measured by $`f`$) in the equilibrium state. In the next section, we will examine more closely the roles played by temperature fluctuations and the negative-feedback mechanism in determining the properties of the IC algorithm.

## III Crossover to Swendsen-Wang behavior

As described in the previous section, the negative-feedback mechanism, inherent in the IC algorithm, drives the system to criticality by effectively adjusting the temperature after each iteration. As previously mentioned, this mechanism leads to differences in the $`L`$-dependence of several dynamic and static quantities from that observed for the SW algorithm. Now, however, we consider an $`l\times l`$ sub-lattice within the $`L\times L`$ lattice. Since the energy of a sub-system of size $`lL`$ is weakly correlated with that of the whole system, the negative feedback mechanism is less effective for the sub-system. A “warm” (relative to $`T_\mathrm{c}`$) sub-system in a “cool” system will be further warmed by the next IC iteration. As a result, the energy autocorrelation time for the sub-system may be longer than for the whole system.

The observation that an iteration of the IC algorithm is equivalent to an iteration of the SW algorithm with $`p=f`$ provides further insight into IC dynamics. In particular, if the distribution of $`f`$ approaches a delta function as $`L\mathrm{}`$, then any finite sub-system of an infinite system will behave exactly as it would under SW dynamics. In short, for $`l`$ sufficiently smaller than $`L`$, the sub-system doesn’t “know” it is being updated by the IC and not the SW algorithm. Thus we expect that, in the limit $`L\mathrm{}`$, all static and dynamic quantities measured on a sub-system of finite size $`l`$ will approach the values measured for the SW algorithm for a sub-system of the same size. It also follows that, for fixed $`L`$, there is a crossover from SW to IC behavior at intermediate values of $`l`$. For example, the integrated autocorrelation time $`\tau _\epsilon (l,L)`$ for the energy per spin in a sub-system of size $`l`$ should initially increase with $`l`$ as $`l^{z^{\mathrm{SW}}}`$ for $`lL`$, reach a maximum, and then decrease as $`l`$ is increased further into the range where the negative-feedback mechanism becomes significant.

We can estimate the length scale at which the crossover occurs as follows. In the canonical ensemble, the temperature uncertainty of the critical region scales with system size $`L`$ as $`\delta TL^{1/\nu }`$. In the IC ensemble, however, temperature fluctuations are governed by the negative-feedback mechanism as described above. In Ref.  it was found that the standard deviation of $`f`$ scales as $`\sigma (f)L^b`$ with $`b.46(.30)`$ for $`q=2(3)`$ as compared with $`1/\nu =1(6/5)`$ respectively. Since the SW algorithm samples from the canonical ensemble and since we expect SW behavior for sub-systems of size $`lL`$, crossover between the IC and SW regimes should occur when the temperature uncertainties from the two sources are comparable. Thus, for a given $`L`$, we expect crossover at a sub-system size $`l_\mathrm{c}`$ given by $`l_\mathrm{c}L^{b\nu }`$.

Therefore, in light of the previous arguments, we hypothesize that the crossover from IC to SW behavior may be described by the scaling relationship

$$A^{\mathrm{IC}}(l,L)=A_{\mathrm{}}^{\mathrm{SW}}(l)F_A(l/L^y),$$
(2)   
where $`A_{\mathrm{}}^{\mathrm{SW}}(l)`$ is any observable measured for the SW algorithm on an $`l\times l`$ sub-lattice immersed in an infinite system, $`A^{\mathrm{IC}}(l,L)`$ is the same observable measured for the IC algorithm on an $`l\times l`$ sub-lattice of an $`L\times L`$ lattice, $`F_A`$ is a scaling function with the property that $`F_A(x)1`$ as $`x0`$, and $`y=b\nu `$. A word of caution concerning boundary conditions is in order here. We define $`A_0^{\mathrm{SW}}(L)`$ to be an observable measured for the SW algorithm on an $`L\times L`$ lattice with periodic boundary conditions. Although we expect $`A_{\mathrm{}}^{\mathrm{SW}}(l)/A_0^{\mathrm{SW}}(l)`$ to approach a constant for $`l\mathrm{}`$, the constant will in general not be exactly 1 due to the different boundary conditions.

In this paper, we also seek to define a meaningful dynamic exponent $`z^{\mathrm{IC}}`$ for the IC algorithm that may be compared with exponents for other algorithms as well as with the Li-Sokal bound. It is not obvious how to do this since the energy autocorrelation time $`\tau _\epsilon `$ for the entire system was observed  to decrease with $`L`$. However, this is not the whole story since we have argued above that correlations between successive IC configurations should decay more slowly on length scales $`l<L`$. Since we would like $`z^{\mathrm{IC}}`$ to describe the $`L`$-dependence of the slowest mode, we suggest that it is the maximum value of $`\tau _\epsilon (l,L)`$ attained for a given $`L`$ that is relevant. Thus we define $`z^{\mathrm{IC}}`$ according to

$$\tau _{\epsilon ,\mathrm{max}}L^{z^{\mathrm{IC}}}.$$
(3)   
In the next section, we describe simulations designed to measure $`z^{\mathrm{IC}}`$ and to test the scaling hypothesis (Eq. (2)) for the static variable $`c(l,L)`$ as well as for the dynamic variable $`\tau _\epsilon (l,L)`$.

## IV Description of simulations

We used the invaded cluster algorithm with the topological spanning rule to simulate the two-dimensional Ising and 3-state Potts models at their critical points for systems ranging in size from $`L=32`$ to $`L=1024`$. Starting from a completely ordered state, we performed a number of relaxation steps to allow the system to reach equilibrium at $`T_\mathrm{c}`$ and then collected data for four observables: the energy per spin $`\epsilon `$, the ratio $`f`$ of occupied to satisfied bonds, the fraction $`m`$ of spins in the largest cluster, and the susceptibility $`\chi `$, given by the sum of the squared cluster sizes divided by the total number of spins. In addition to measuring the mean value and variance for each observable, we also measured the autocorrelation function and used this to calculate integrated autocorrelation times. For a given observable $`A`$, the (normalized) autocorrelation function at a given time step $`t`$ can be calculated from a sequence of $`n`$ MC measurements $`\{A(j),j=1,\mathrm{},n\}`$ according to

$$\mathrm{\Gamma }_A(t)\frac{_{j=1}^{nt}(A(j)A)(A(j+t)A)}{_{j=1}^n(A(j)A)^2},$$
(4)   
where $`A`$ is the mean value of $`A`$.

The integrated autocorrelation time for the observable $`A`$ is defined by

$$\tau _A\frac{1}{2}+\underset{t=1}{\overset{\mathrm{}}{}}\mathrm{\Gamma }_A(t).$$
(5)   
Obviously, in practice, the sum must be truncated at some reasonable value of $`t`$. Following the recommendation of Ref. , we define

$$\tau _A(t_A^{})\frac{1}{2}+\underset{t=1}{\overset{t_A^{}}{}}\mathrm{\Gamma }_A(t)$$
(6)   
and choose the cutoff $`t_A^{}`$ to be the smallest integer such that $`t_A^{}\kappa \tau _A(t_A^{})`$, where $`\kappa `$ is a constant whose value will be discussed shortly. If the autocorrelation function has the scaling form $`\mathrm{\Gamma }_A(t)=G(t/\tau _{\mathrm{exp}})`$, where $`\tau _{\mathrm{exp}}`$ is the exponential autocorrelation time, then choosing the cutoff in this manner will insure that $`\tau _A(t_A^{})`$ is proportional to $`\tau _A`$. Thus estimates of $`z_A`$ will not be biased by truncating the sum at $`t=t_A^{}`$.

One also would like the values of $`\tau _A(t_A^{})`$ to approximate $`\tau _A`$ as accurately and precisely as possible, and here there is a tradeoff between excluding noise and including as much of the signal as possible. In Ref.  it is shown that if $`\mathrm{\Gamma }_A(t)`$ is roughly a single exponential, then choosing a value of $`\kappa `$ in the range 4–6 would achieve the optimal compromise for $`n/\tau `$ in the range $`10^4`$$`10^6`$ that we used in our simulations. However, although $`\mathrm{\Gamma }_A(t)`$ is well approximated by a single exponential in the case of the SW algorithm, this is not true for the IC algorithm . For this reason, and since we are willing to accept slightly larger statistical uncertainties in order to reduce systematic errors, we used $`\kappa =10`$ in all our calculations.

In light of the discussion in Sec. III, we also collected data for several sub-system sizes for each $`L`$. Here we concentrated on the energy per spin, measuring $`\epsilon (l,L)`$, $`\mathrm{var}(\epsilon (l,L))`$, and $`\tau _\epsilon (l,L)`$ for sub-system sizes ranging from $`l=1`$ to $`l=L/2`$ for each $`L`$. The sub-systems are squares all sharing a single corner of the lattice.

Error bars on all quantities were calculated using the blocking method. Each run was partitioned into $`k`$ contiguous blocks of $`n`$ MC steps each, and the individual blocks treated as independent runs. Although this is an approximation, it will be a good one provided that $`n`$ is large compared to the system’s longest relaxation time. As an example of the blocking method, we obtain the value for, say, $`\tau _m`$ by first calculating $`\tau _m^{(i)}`$ for each block $`i`$. We then calculate the mean

$$\overline{\tau _m}\frac{1}{k}\underset{i=1}{\overset{k}{}}\tau _m^{(i)}$$
(7)   
and its standard error

$$\sigma (\overline{\tau _m})\sqrt{\frac{_{i=1}^k(\tau _m^{(i)}\overline{\tau _m})^2}{k(k1)}}$$
(8)   
and report the result $`\tau _m=\overline{\tau _m}\pm \sigma (\overline{\tau _m})`$.

For both $`q=2`$ and $`q=3`$, one long run was initially performed for each lattice size $`L`$, and the blocking method implemented as just described. (Three independent runs were performed for the case $`q=2`$, $`L=1024`$.) The number of blocks used ranged from 100 for $`L=32`$ down to 10 for $`L=1024`$, and the number $`n`$ of MC steps per block ranged from $`5\times 10^3`$ to $`1\times 10^5`$. In each case, $`n`$ was greater than the longest observed autocorrelation time for the given system by at least a factor of $`10^3`$ ($`10^4`$ for the smaller lattices), thereby making the assumption of independent blocks, used in calculating the error bars, a reasonably good one. The number of equilibration steps performed at the beginning of each run also exceeded the longest observed $`\tau `$ by a factor of $`10^3`$ in all cases.

In these initial runs, we collected data for sub-system sizes $`l\{1,2,4,\mathrm{},L/2\}`$ as well as for the whole system ($`l=L`$). In order to estimate $`z^{\mathrm{IC}}`$ as defined in Eq. (3), we sought to obtain an accurate value for $`\tau _{\epsilon ,\mathrm{max}}(L)`$ for each $`L`$. Therefore, once we had learned, from the initial runs, the approximate sub-system size $`l_{\mathrm{max}}`$ at which $`\tau _\epsilon (l,L)`$ attains a maximum, we then performed between one and three additional independent runs for each system and collected data for evenly spaced values of $`l`$ near our rough estimate of $`l_{\mathrm{max}}`$. The entire experiment required about 5 months of CPU time on a single processor of a dual-processor 266MHz Pentium II Linux workstation. We used the machine-supplied random number generator coupled with a shuffling procedure as described in Ref. .

## V Discussion of results

The first quantity we consider is the static variable $`c(L)L^d\mathrm{var}(\epsilon )`$ (see Table I). In the canonical ensemble, $`c(L)`$ is the specific heat which diverges, at the critical point, as $`\mathrm{log}L`$ for $`q=2`$ and as $`L^{\alpha /\nu }`$ with $`\alpha /\nu =2/5`$ for $`q=3`$. We see from Fig. 1, however, that the situation is quite different for the IC ensemble, as first observed in Ref. . We assume that the asymptotic behavior is given by a power law $`c(L)L^w`$ and fit a line to a plot of $`\mathrm{log}_{10}c(L)`$ vs $`\mathrm{log}_{10}L`$ as shown in Fig. 1. Note that where error bars are not visible they are smaller than the symbol height. As is always the case when trying to ascertain asymptotic behavior from simulations at finite $`L`$, there can be some debate as to which, if any, data points should be omitted from the fit because of corrections to scaling. Here and elsewhere in our analysis, we proceed by dropping points one at a time in order of increasing L until either 1) a reasonably good fit is obtained, 2) the fit ceases to improve significantly with further cuts, or 3) we are left with only three data points. We employ standard, weighted $`\chi ^2`$ fitting, using the confidence level <sup>*</sup><sup>*</sup>*The confidence level is the probability that a $`\chi ^2`$ as poor as the measured value would occur, assuming that the underlying model is correct and that the measurement errors are normally distributed . (Note that confidence level is denoted by the symbol $`Q`$ in Ref. .) (CL) as our goodness-of-fit measure, and consider a fit to be “reasonably good” if $`\mathrm{CL}10\%`$.

For the Ising model, a fit to the last four data points ($`L128`$) yields $`w=1.020\pm .003`$ ($`\mathrm{CL}=16\%`$) in agreement with the observation $`w1`$ reported in Ref. . In the case of the 3-state Potts model, a fit to the last three points gives $`w=1.313\pm .008`$, but, because of the poor confidence level (2%) and the upward curvature visible in the data, this value should probably just be regarded as a lower bound on $`w`$ for $`q=3`$. We emphasize that the error bars on these and subsequent exponent estimates are purely statistical in nature and do not reflect the uncertainty of extrapolating to infinite system size.

Turning now to the dynamic variables (see Table II), we plot the logarithms of the autocorrelation times $`\tau _\epsilon `$, $`\tau _f`$, $`\tau _m`$, and $`\tau _\chi `$ vs. $`\mathrm{log}_{10}L`$ for $`q=2`$ in Fig. 2 and for $`q=3`$ in Fig. 3. For $`q=2`$ we find that $`\tau _\epsilon `$ and $`\tau _f`$ decrease with $`L`$ (perhaps in a rather complicated fashion) and $`\tau _m`$ and $`\tau _\chi `$ remain approximately constant for the range of $`L`$ values used in our simulations. These results are in agreement with initial observations  that led to speculation of no critical slowing. For $`q=3`$, however, we observe critical slowing in the behavior of $`\tau _m`$ and $`\tau _\chi `$. Fits to the data for $`L64`$ yield $`z_m=.191\pm .004`$ ($`\mathrm{CL}=25\%`$) and $`z_\chi =.206\pm .005`$ ($`\mathrm{CL}=20\%`$), but it seems likely that $`\tau _m\tau _\chi `$ in the asymptotic limit. We note that the value of $`z_m`$ is somewhat smaller than the previous estimate $`z_m.28`$ .

Next we consider $`\tau _\epsilon (l,L)`$, the integrated autocorrelation time for the energy per spin $`\epsilon (l,L)`$ measured on an $`l\times l`$ sub-lattice of the whole $`L\times L`$ lattice. As expected from the discussion in Sec. III, we see that as $`l`$ is increased for a given $`L`$, $`\tau _\epsilon (l,L)`$ increases, reaches a maximum, and then decreases, as shown in Fig. 4 for $`q=2`$ and in Fig. 5 for $`q=3`$.

To find $`\tau _{\epsilon ,\mathrm{max}}`$ and its location $`l_{\mathrm{max}}`$ for a given $`L`$ (see Table III), we fit a parabola to the region of the curve near the maximum. In order to do this objectively, we began, for each $`L`$, by omitting the data point with the smallest value of $`\tau _\epsilon `$ and performing the fit. We then dropped the point with the next smallest $`\tau _\epsilon `$, re-fit, and continued in this fashion until 1) a CL of 50% or greater was obtained and 2) the values of $`\tau _{\epsilon ,\mathrm{max}}`$ and $`l_{\mathrm{max}}`$ that were obtained by dropping an additional point remained within error bars of the current best-fit values.

We then attempted to determine the dynamic exponent $`z^{\mathrm{IC}}`$, defined in Eq. (3), by fitting a line to a plot of $`\mathrm{log}_{10}\tau _{\epsilon ,\mathrm{max}}(L)`$ vs. $`\mathrm{log}_{10}L`$ for $`q=2`$ and $`q=3`$ . The results are shown in Fig. 6 along with results for $`\tau _\epsilon ^{\mathrm{SW}}`$ taken from Baillie and Coddington  for $`q=2`$ and from Salas and Sokal  for $`q=3`$. We note that the increase of $`\tau _{\epsilon ,\mathrm{max}}(L)`$ with $`L`$ for $`q=2`$ is the first observation of critical slowing for the IC algorithm in the case of the Ising model.

For $`q=2`$ the autocorrelation times for the IC algorithm are nearly the same as for the Wolff algorithm and smaller than those for the SW algorithm by a factor of about 2 (see Table III), but the $`L`$-dependence is similarly obscure. A good fit to a power law could not be obtained for the IC data for $`q=2`$. The line shown in the figure, having slope $`.21`$ is the best fit for $`L64`$, but it clearly does not describe the data very well. The best power-law fit to the SW data for $`64L512`$ (also shown in the figure) yields $`z^{\mathrm{IC}}.25`$ as reported in Ref. . Although the fit is considerably better than that for the IC data, the CL is still poor ($`<0.1\%`$), and the better fit might be primarily due to the absence of data for $`L=1024`$ in the SW case.

Since it has been suggested  that $`\tau _\epsilon ^{\mathrm{SW}}`$ increases logarithmically with $`L`$ rather than as a power of $`L`$, Baillie and Coddington also fit their data to a logarithm with somewhat better results ($`\mathrm{CL}=13\%`$). For the IC algorithm, a logarithmic fit is still atrocious, albeit somewhat better than the power law. Later in this section we consider the possibility that we have underestimated the error bars on $`\tau _{\epsilon ,\mathrm{max}}`$, which, of course, could result in a poor fit even if the underlying model had been correctly identified. Still, even the general trend in the data is difficult to discern, indicating that corrections to scaling are probably significant for the system sizes studied here. Therefore, we conclude that high precision data for larger lattices are needed before a more definitive statement can be made concerning the asymptotic behavior of $`\tau _{\epsilon ,\mathrm{max}}`$.

For the 3-state model, however, the picture appears to be somewhat clearer. Although a slight downward curvature in the data is visible in Fig. 6, a good fit ($`\mathrm{CL}=69\%`$) to a power law is obtained for $`256L1024`$, yielding $`z^{\mathrm{IC}}=.346\pm .002`$. This is to be compared with the value of $`z_\epsilon ^{\mathrm{SW}}=.515\pm .006`$ obtained by Salas and Sokal for the SW algorithm with $`128L1024`$ ($`\mathrm{CL}=80\%`$). For the Wolff algorithm $`z_\epsilon ^{\mathrm{Wolff}}=.57\pm .01`$ was reported in Ref. . While there is no guarantee that $`\tau _{\epsilon ,\mathrm{max}}`$ is the system’s longest relaxation time, it is interesting that $`z^{\mathrm{IC}}`$ is significantly smaller than $`z_\epsilon ^{\mathrm{SW}}`$ and also less than the Li-Sokal bound ($`z^{\mathrm{SW}}\alpha /\nu =2/5`$).

Now we proceed to test the scaling hypothesis (Eq. (2)) presented in Sec. III. In order to do this, we first need to determine the exponent $`b`$ defined by $`\sigma (f)L^b`$. We plot $`\mathrm{log}_{10}\sigma (f)`$ vs. $`\mathrm{log}_{10}L`$ in Fig. 7 (the data are listed in Table I) for $`q=2`$ and $`q=3`$ along with the best-fit lines. A fit to all six $`q=2`$ points yields $`b=.4781\pm .0006`$ ($`\mathrm{CL}=54\%`$), while, for $`q=3`$ a good fit ($`\mathrm{CL}=68\%`$) is obtiained for the last four points ($`L128`$), resulting in $`b=.3252\pm .0009`$. These results are consistent with previous estimates .

We now test Eq. (2) for the variable $`c^{\mathrm{IC}}(l,L)l^d\mathrm{var}(\epsilon (l,L))`$. Since $`c_0^{\mathrm{SW}}(l)\mathrm{log}l`$ and we expect $`c_{\mathrm{}}^{\mathrm{SW}}(l)c_0^{\mathrm{SW}}(l)`$ (recall $`c_0^{\mathrm{SW}}(l)`$ is the specific heat for an $`l\times l`$ lattice with periodic boundary conditions, and $`c_{\mathrm{}}^{\mathrm{SW}}(l)`$ is the specific heat for an $`l\times l`$ sub-lattice immersed in an infinite system), we plot $`c^{\mathrm{IC}}(l,L)/\mathrm{log}_{10}l`$ vs. $`l/L^y`$ in Fig. 8, where $`y=b=.4781`$ for the Ising model ($`\nu =1`$). The observed data collapse provides strong support for Eq. (2). For $`q=3`$, however, it was found in Ref.  that the asymptotic form $`c_0^{\mathrm{SW}}(l)l^{\alpha /\nu }`$ does not describe the SW data very well for the range of lattice sizes considered here. Therefore, we cannot expect Eq. (2) to provide a good description of the IC data if the asymptotic form is used for $`c_{\mathrm{}}^{\mathrm{SW}}(l)`$. Nevertheless, if we plot $`c^{\mathrm{IC}}(l,L)/c_0^{\mathrm{SW}}(l)`$ vs. $`l/L^y`$, using the measured values of $`c_0^{\mathrm{SW}}(l)`$ from Ref.  and $`y=b\nu =.2711`$ for $`q=3`$, data collapse is apparent in Fig. 9, although perhaps a bit less convincing than for the Ising case.

We note that the curve in Fig. 9 extrapolates to about $`.85`$ on the vertical axis for $`l/L^y=0`$. As previously mentioned, we would expect this value to be 1 if SW data were collected for sub-systems immersed in larger systems so as to reproduce the boundary conditions applied to the IC sub-systems. Thus we conclude that our scaling hypothesis (Eq. (2)) does appear to be valid for the variable $`c^{\mathrm{IC}}(l,L)`$ for $`q=2`$ and $`q=3`$.

Although the static quantity $`c^{\mathrm{IC}}(l,L)`$ seems to be well described by Eq. (2) with $`y=b\nu `$, the same is not true for the dynamic quantity $`\tau _\epsilon (l,L)`$. This is easy to see, since Eq. (2) predicts that the location $`l_{\mathrm{max}}`$ of the maximum in Figs. 4 and 5 should scale as $`L^y`$; however, the plots of $`\mathrm{log}_{10}l_{\mathrm{max}}`$ vs. $`\mathrm{log}_{10}L`$ shown in Fig. 10 reveal that this is not the case—at least not if $`y=b\nu `$ is required. For both $`q=2`$ and $`q=3`$ the slope of the best-fit line for $`64L1024`$, shown in Fig. 10, is approximately $`.62`$, although the confidence levels are poor. Unlike the situation encountered earlier in this section, when fitting $`\tau _{\epsilon ,\mathrm{max}}`$ to a power law in $`L`$, the points seem to be scattered randomly about the best-fit line. Therefore, we suspect that $`l_{\mathrm{max}}`$ does scale as a power of $`L`$, but that our error bars on $`l_{\mathrm{max}}`$ are somewhat underestimated.

There are three aspects of our analysis that could lead to underestimates in the error bars on $`\tau _{\epsilon ,\mathrm{max}}`$ and $`l_{\mathrm{max}}`$. First, there always exists the possibility that the assumption of normally distributed measurement errors is not valid. Second, the blocking method, used to calculate error bars on values of $`\tau _\epsilon (l,L)`$, treats successive blocks as if they were independent runs, an approximation that may not be entirely justified even though the block length was greater than $`10^3\tau _{\epsilon ,\mathrm{max}}`$ in all cases. Finally, and probably most importantly, the error bars on $`\tau _{\epsilon ,\mathrm{max}}`$ and $`l_{\mathrm{max}}`$ resulting from the weighted $`\chi ^2`$ fit to a parabola are calculated by assuming that the measurements of $`\tau _\epsilon (l,L)`$ at different values of $`l`$ for a given $`L`$ are independent. This is clearly not a good approximation, since all of the sub-lattices extend outward from the same corner of the $`L\times L`$ lattice. Therefore, all the spins in a given $`l\times l`$ sub-system are also contained in every larger sub-system, and thus sub-systems for comparable values of $`l`$ are highly correlated.

In any case, it is still clear that Eq. (2) with $`y=b\nu `$ cannot explain the data for $`\tau _\epsilon `$. Since the logic leading to the scaling form seems sound, we hypothesize that Eq. (2) does hold for $`\tau _\epsilon `$ but that the crossover length for $`\tau _\epsilon `$ is different from that for $`c`$ so that $`yb\nu `$ in the case of $`\tau _\epsilon `$. This seems plausible, since there is no a priori reason why the thermodynamic argument by which we arrived at $`y=b\nu `$ must apply to the dynamic quantity $`\tau _\epsilon `$. Nevertheless, if Eq. (2) still holds for $`\tau _\epsilon `$, we can obtain the crossover exponent $`y`$ from Fig. 10 as described above.

To test our scaling hypothesis for $`\tau _\epsilon `$, we plot $`\tau _\epsilon (l,L)/\tau _{\epsilon ,0}^{\mathrm{SW}}(l)`$ vs. $`l/L^y`$ with $`y=.6176`$ ($`y=.626`$) for $`q=2`$ ($`q=3`$) in Fig. 11 (Fig. 12). The values of $`\tau _{\epsilon ,0}^{\mathrm{SW}}(l)`$ are taken from Ref.  for $`q=2`$ and from Ref.  for $`q=3`$. The data collapse is not terribly convincing in either case, but seems too good to completely rule out Eq. (2) as the correct asymptotic form. The fact that both curves extrapolate to about 1 for $`l/L^y=0`$ provides further support for the scaling hypothesis. Still, it appears that additional tests are needed to confirm or disprove Eq. (2) for $`\tau _\epsilon `$.

## VI Conclusion

Using the invaded cluster (IC) algorithm with the topological spanning rule, we simulated the critical Ising and 3-state Potts models in two dimensions for systems ranging in size from $`L=32`$ to $`L=1024`$. In accord with previous results , we find that the $`L`$-dependence of several static and dynamic quantities is very different from that observed for the Swendsen-Wang (SW) algorithm which samples from the canonical ensemble. In particular, the quantity $`c(L)L^d\mathrm{var}(\epsilon )`$ is not proportional to the specific heat and the integrated autocorrelation time $`\tau _\epsilon `$ for the energy per spin decreases with $`L`$. However, we find that the corresponding quantities $`\tau _\epsilon (l,L)`$ and $`c(l,L)`$, measured for a sub-system of size $`l`$, exhibit a crossover to SW behavior for $`l`$ sufficiently less than $`L`$.

To describe the crossover phenomenon, we propose the scaling form $`A^{\mathrm{IC}}(l,L)=A_{\mathrm{}}^{\mathrm{SW}}(l)F_A(l/L^y)`$, where $`A_{\mathrm{}}^{\mathrm{SW}}(l)`$ is any observable measured for the SW algorithm on an $`l\times l`$ sub-lattice immersed in an infinite system, $`A^{\mathrm{IC}}(l,L)`$ is the same observable measured for the IC algorithm on an $`l\times l`$ sub-lattice of an $`L\times L`$ lattice, and $`F_A`$ is a scaling function with the property that $`F_A(x)1`$ as $`x0`$. We have argued that the crossover exponent $`y`$ should equal $`b\nu `$, where $`\nu `$ is the usual correlation-length exponent and $`b`$ is defined by $`\sigma (f)L^b`$ with $`f`$ the ratio of occupied to satisfied bonds. We find that the proposed scaling form with $`y=b\nu `$ provides a good description of our data for the static variable $`c(l,L)`$, but is less successful for the dynamic variable $`\tau _\epsilon (l,L)`$, even if the possibility $`yb\nu `$ is admitted.

In addition, we define the dynamic exponent $`z^{\mathrm{IC}}`$ for the invaded cluster algorithm in terms of the maximum value $`\tau _{\epsilon ,\mathrm{max}}`$ attained for a given $`L`$ according to $`\tau _{\epsilon ,\mathrm{max}}L^{z^{\mathrm{IC}}}`$. For $`q=3`$ we find that $`z^{\mathrm{IC}}=.346\pm .002`$, which is smaller than recent numerical estimates of $`z^{\mathrm{SW}}`$ and $`z^{\mathrm{Wolff}}`$ and also less than the Li-Sokal bound on $`z^{\mathrm{SW}}`$. For $`q=2`$ we also observe critical slowing, but the $`L`$-dependence of $`\tau _{\epsilon ,\mathrm{max}}`$ is less clear. It appears from our simulations that $`z^{\mathrm{IC}}<.21`$ and possibly that $`\tau _{\epsilon ,\mathrm{max}}\mathrm{log}L`$ ($`z^{\mathrm{IC}}=0`$), but high precision data for larger lattices are needed before a more definitive statement can be made concerning the asymptotic behavior of $`\tau _{\epsilon ,\mathrm{max}}`$ for the Ising model.

## Acknowledgements

We are grateful to Yongsoo Choi for providing the original code for the invaded cluster algorithm and to Robert Guyer for useful discussions concerning the analysis. This work was supported in part by NSF Grant DMR-9632898 and NSA Grant MDA 904-98-1-0518.