Dynamic Critical Behavior of the Swendsen-Wang Algorithm: The Two-Dimensional 3-State Potts Model Revisited

TL;DR

This study investigates the dynamic critical behavior of the Swendsen-Wang algorithm applied to the 2D 3-state Potts model, revealing that the Li-Sokal bound is nearly sharp with a divergence in the ratio of autocorrelation time to specific heat.

Contribution

The paper provides high-precision Monte Carlo results showing the near-sharpness of the Li-Sokal bound and the divergence behavior of the autocorrelation time ratio.

Findings

Li-Sokal bound is almost sharp but not quite

The ratio τ_{int,E} / C_H diverges as a small power or logarithm

Autocorrelation times grow faster than specific heat near criticality

Abstract

We have performed a high-precision Monte Carlo study of the dynamic critical behavior of the Swendsen-Wang algorithm for the two-dimensional 3-state Potts model. We find that the Li-Sokal bound ($\tau_{int,E} \geq const \times C_H$) is almost but not quite sharp. The ratio $\tau_{int,E} / C_H$ seems to diverge either as a small power ($\approx 0.08$) or as a logarithm.

Tables

Table 1. Table 1: Static data from the MC simulations at the critical point of the 3-state Potts model. For each lattice size ( L 𝐿 L ), we include the number of performed measurements ( M ​ C ​ S 𝑀 𝐶 𝑆 MCS ) in units of 10 6 superscript 10 6 10^{6} , the susceptibility ( χ 𝜒 \chi ), the specific heat ( C H subscript 𝐶 𝐻 C_{H} ), the second-moment correlation length ( ξ 𝜉 \xi ), and the energy ( E 𝐸 E ). The quoted errors correspond to one standard deviation (i.e. confidence level ≈ \approx 68%). The first row (“exact”) gives the exact results for the 4 × 4 4 4 4\times 4 lattice, and the last row (“ L = ∞ 𝐿 L=\infty ”) gives the exact energy in the limit L → ∞ → 𝐿 L\rightarrow\infty .

$L$	$MCS$	$\chi$		$C_H$		$\xi$		$E$
4	exact	111112.204711		1.496719		3.862380		1.710062
4	4.9	12.2093$\pm$	0.0054	1.4952$\pm$	0.0021	3.8635$\pm$	0.0043	1.71051$\pm$	0.00039
8	6.9	41.3527$\pm$	0.0204	2.4712$\pm$	0.0031	7.5419$\pm$	0.0072	1.65391$\pm$	0.00026
16	9.9	138.4491$\pm$	0.0723	3.7162$\pm$	0.0046	14.9267$\pm$	0.0134	1.62070$\pm$	0.00016
32	14.9	462.7014$\pm$	0.2408	5.2941$\pm$	0.0064	29.8509$\pm$	0.0250	1.60223$\pm$	0.00010
64	19.9	1542.6921$\pm$	0.8428	7.3839$\pm$	0.0094	59.6964$\pm$	0.0501	1.59163$\pm$	0.00006
128	29.9	5135.9512$\pm$	2.7481	10.1007$\pm$	0.0127	119.3401$\pm$	0.0952	1.58552$\pm$	0.00003
256	40.8	17082.8221$\pm$	9.2802	13.7116$\pm$	0.0176	238.4565$\pm$	0.1888	1.58202$\pm$	0.00002
512	12.9	56760.2838$\pm$	65.0770	18.4752$\pm$	0.0507	475.9494$\pm$	0.7779	1.58001$\pm$	0.00003
1024	5.5	189676.5530$\pm$	387.3500	24.5281$\pm$	0.1247	958.8929$\pm$	2.7875	1.57893$\pm$	0.00003
$\infty$	exact							1.577350

Table 2. Table 2: Autocorrelation times for the runs performed at the critical point of the 3-state Potts model. For each lattice size ( L 𝐿 L ), we include the integrated autocorrelation times for the squared magnetization ( τ int , ℳ 2 subscript 𝜏 int superscript ℳ 2 \tau_{{\rm int},{\cal M}^{2}} ), the bond occupation ( τ int , 𝒩 subscript 𝜏 int 𝒩 \tau_{{\rm int},{\cal N}} ), the energy ( τ int , ℰ subscript 𝜏 int ℰ \tau_{{\rm int},{\cal E}} ), and the nearest-neighbor connectivity ( τ int , ℰ ′ subscript 𝜏 int superscript ℰ ′ \tau_{{\rm int},{\cal E}^{\prime}} ).

$L$	${\tau }_{\mathrm{int},{ℳ}^2}$		${\tau }_{\mathrm{int},𝒩}$		${\tau }_{\mathrm{int},ℰ}$		${\tau }_{\mathrm{int},{ℰ}^{\prime }}$
4	4.013$\pm$	0.018	3.205$\pm$	0.013	4.023$\pm$	0.018	4.034$\pm$	0.018
8	5.983$\pm$	0.028	5.194$\pm$	0.023	6.033$\pm$	0.028	6.080$\pm$	0.028
16	8.816$\pm$	0.041	8.084$\pm$	0.036	9.025$\pm$	0.043	9.117$\pm$	0.043
32	12.648$\pm$	0.057	12.200$\pm$	0.055	13.280$\pm$	0.062	13.438$\pm$	0.063
64	18.101$\pm$	0.085	18.313$\pm$	0.086	19.549$\pm$	0.095	19.769$\pm$	0.097
128	25.665$\pm$	0.117	27.094$\pm$	0.127	28.525$\pm$	0.137	28.820$\pm$	0.139
256	35.691$\pm$	0.164	39.204$\pm$	0.189	40.824$\pm$	0.200	41.212$\pm$	0.203
512	49.775$\pm$	0.480	56.667$\pm$	0.583	58.511$\pm$	0.612	58.974$\pm$	0.619
1024	68.387$\pm$	1.183	80.489$\pm$	1.511	82.488$\pm$	1.567	83.030$\pm$	1.583

Table 3. Table 3: Ratios τ int , ℰ / C H subscript 𝜏 int ℰ subscript 𝐶 𝐻 \tau_{{\rm int},{\cal E}}/C_{H} for the runs performed at the critical point of the 3-state Potts model. In parentheses we give our error-bar estimates. The first number shows the error bar coming from performing the “bunch” method with 100 bunches. The second number is obtained by using the triangular inequality with the numerical results coming from the standard method (see Section  4 ).

$L$	${\tau }_{\mathrm{int},ℰ}/C_H$
4	2.6906	$\pm$ (0.0111,	$\le$ 0.0161)
8	2.4413	$\pm$ (0.0092,	$\le$ 0.0145)
16	2.4285	$\pm$ (0.0090,	$\le$ 0.0145)
32	2.5084	$\pm$ (0.0084,	$\le$ 0.0147)
64	2.6476	$\pm$ (0.0111,	$\le$ 0.0163)
128	2.8241	$\pm$ (0.0111,	$\le$ 0.0171)
256	2.9773	$\pm$ (0.0113,	$\le$ 0.0184)
512	3.1670	$\pm$ (0.0281,	$\le$ 0.0418)
1024	3.3630	$\pm$ (0.0442,	$\le$ 0.0810)

Table 4. Table 4: Numerical estimates for the static and dynamic critical exponents of the 3-state Potts model (second column). We also include the exact results for comparison (third column).

Exponent	numerical	exact
$\gamma /\nu$	$1.73444\pm 0.00043$	26/15
$\alpha /\nu$	$0.4240\pm 0.0030$	2/5
$z_{\mathrm{int},ℰ}$	$0.515\pm 0.006$	$\ge 2/5$
$z_{\mathrm{int},{ℰ}^{\prime }}$	$0.514\pm 0.006$	$\ge 2/5$
$z_{\mathrm{int},𝒩}$	$0.529\pm 0.006$	$\ge 2/5$
$z_{\mathrm{int},{ℳ}^2}$	$0.475\pm 0.006$

Equations

$${\tau }_{\mathrm{int},\mathcal{N}},{\tau }_{\mathrm{int},\mathcal{E}},{\tau }_{\exp }\ge \mathrm{const}\times C_H.$$

$$z_{\mathrm{int},\mathcal{N}},z_{\mathrm{int},\mathcal{E}},z_{\exp }\ge \frac{\alpha }{\nu },$$

$${\mathcal{H}}_{\mathrm{Potts}}=-\beta \sum _{⟨ij⟩}({\delta }_{{\sigma }_i,{\sigma }_j}-1),$$

$$Z=\sum _{\{\sigma \}}{e}^{-\mathcal{H}}=\sum _{\{\sigma \}}\exp \left[\beta \sum _b({\delta }_{{\sigma }_b}-1)\right].$$

$$W_{\mathrm{Potts}}(\{\sigma \})=\frac{1}{Z}\prod _b{e}^{\beta ({\delta }_{{\sigma }_b}-1)}=\frac{1}{Z}\prod _b(1-p+p{\delta }_{{\sigma }_b})$$

$$W_{\mathrm{joint}}(\{\sigma \},\{n\})=\frac{1}{Z}\prod _b\left[(1-p){\delta }_{n_b,0}+p{\delta }_{{\sigma }_b}{\delta }_{n_b,1}\right].$$

$$W_{\mathrm{RC}}(\{n\})=\frac{1}{Z}\left[\prod _{b:n_b=1}p\right]\left[\prod _{b:n_b=0}(1-p)\right]q^{\mathcal{C}(\{n\})}$$

$$P_{SW}=P_{\mathrm{bond}}{P}_{\mathrm{spin}},$$

$$C_{\mathcal{OO}}(t)\equiv ⟨\mathcal{O}(s)\mathcal{O}(s+t)⟩-⟨\mathcal{O}{⟩}^2$$

$${\rho }_{\mathcal{OO}}(t)\equiv \frac{C_{\mathcal{OO}}(t)}{C_{\mathcal{OO}}(0)}=\frac{C_{\mathcal{OO}}(t)}{var(\mathcal{O})}.$$

$$\mathcal{N}\equiv \sum _b{n}_b.$$

$$E(n_{b}|\{\sigma\})&=p\delta_{\sigma_{b}}\\ E(n_{b}n_{b^{\prime}}|\{\sigma\})=\left\{\begin{array}[]{ll}p^{2}\delta_{\sigma_{b}}\delta_{\sigma_{b^{\prime}}}&\quad\hbox{for $b\neq b^{\prime}$}\\ p\delta_{\sigma_{b}}&\quad\hbox{for $b=b^{\prime}$}\end{array}\right.$$

$$\langle{\cal N}\rangle&=p\langle{\cal E}\rangle{}\\ \langle{\cal N}^{2}\rangle=p^{2}\langle{\cal E}^{2}\rangle+p(1-p)\langle{\cal E}\rangle\\ C_{\cal NN}(0)\equiv\hbox{var}({\cal N})=p^{2}\hbox{var}({\cal E})+p(1-p)\langle{\cal E}\rangle$$

$$\mathcal{E}\equiv \sum _b{\delta }_{{\sigma }_b}.$$

$$C_{\mathcal{NN}}(1)\equiv ⟨\mathcal{N}(0)\mathcal{N}(1)⟩-⟨\mathcal{N}{⟩}^2=var(P_{\mathrm{bond}}\mathcal{N})=var(E(\mathcal{N}|\{\sigma \})).$$

$$P_{\mathrm{bond}}\mathcal{N}\equiv E(\mathcal{N}|\{\sigma \})=p\sum _b{\delta }_{{\sigma }_b}=p\mathcal{E}.$$

$$C_{\mathcal{NN}}(1)=p^{2}var(\mathcal{E})$$

$${\rho }_{\mathcal{NN}}(1)\equiv \frac{C_{\mathcal{NN}}(1)}{C_{\mathcal{NN}}(0)}=\frac{p{C}_H}{p{C}_H+(1-p)E}=1-\frac{(1-p)E}{p{C}_H+(1-p)E},$$

$$E$$

$$C_H$$

$${\rho }_{\mathcal{NN}}(t)={\int }_0^1{\lambda }^{|t|}d\nu (\lambda )$$

$${\rho }_{\mathcal{NN}}(t)\ge {\rho }_{\mathcal{NN}}(1{)}^{|t|}.$$

$${\tau }_{\mathrm{int},\mathcal{N}}$$

$${\tau }_{\exp ,\mathcal{N}}$$

$${\tau }_{\mathrm{int},\mathcal{N}}$$

$${\tau }_{\exp ,\mathcal{N}}$$

$$C_{\mathcal{EE}}(t)=p^{-2}{C}_{\mathcal{NN}}(t+1)$$

$${\rho }_{\mathcal{EE}}(t)=\frac{{\rho }_{\mathcal{NN}}(t+1)}{{\rho }_{\mathcal{NN}}(1)}\ge {\rho }_{\mathcal{NN}}(t).$$

$${\tau }_{\exp ,\mathcal{E}}={\tau }_{\exp ,\mathcal{N}},$$

$${\tau }_{\mathrm{int},\mathcal{N}}\le {\tau }_{\mathrm{int},\mathcal{E}}=\frac{{\tau }_{\mathrm{int},\mathcal{N}}-\frac{1}{2}}{{\rho }_{\mathcal{NN}}(1)}-\frac{1}{2}\le \frac{{\tau }_{\mathrm{int},\mathcal{N}}}{{\rho }_{\mathcal{NN}}(1)}.$$

$$z_{\mathrm{int},\mathcal{E}}=z_{\mathrm{int},\mathcal{N}}.$$

$${\mathcal{E}}^{\prime }\equiv \sum _b{\gamma }_b,$$

Full text

Key Words: Potts model; Swendsen–Wang algorithm; cluster algorithm; Monte Carlo; Li–Sokal bound; dynamic critical phenomena.

1 Introduction

Monte Carlo (MC) simulations [1, 2, 3, 4] have become a standard and powerful tool for gaining new insights into statistical-mechanical systems and lattice field theories. However, their practical success is severely limited by critical slowing-down: the autocorrelation time $\tau$ — that is, roughly speaking, the time needed to produce one “statistically independent” configuration — diverges near a critical point. More precisely, for a finite system of linear size $L$ at criticality, we expect a behavior $\tau\sim L^{z}$ for large $L$. The power $z$ is a dynamic critical exponent, and it depends on both the system and the algorithm.

Single-site MC algorithms (such as single-site Metropolis or heat bath) have a dynamic critical exponent z\mathrel{\hbox to0.0pt{\lower 3.0pt\hbox{\mathchar 536\relax}\hss}\raise 2.0pt\hbox{\mathchar 318\relax}}2. This makes it very hard to get high-precision data very close to the critical point on large lattices.

In some cases, a much better dynamical behavior is obtained by allowing non-local moves, such as cluster flips.111 See [4, 5, 6] for reviews of collective-mode Monte Carlo methods.

The Swendsen–Wang (SW) cluster algorithm [7] for the $q$-state ferromagnetic Potts model achieves a significant reduction in $z$ compared to the local algorithms: one has $z$ between 0 and $\approx\!1$, where the exact value depends on $q$ and on the dimensionality of the lattice [6]. The most favorable case is the two-dimensional (2D) Ising model ($q=2$), in which z\mathrel{\hbox to0.0pt{\lower 3.0pt\hbox{\mathchar 536\relax}\hss}\raise 2.0pt\hbox{\mathchar 316\relax}}0.3 [7, 8, 9, 10] (see below). In other cases, the performance of the SW algorithm is less impressive (though still quite good): e.g., $z=0.55\pm 0.03$ for the 2D 3-state Potts model [11], $z\approx 1$ for the 2D 4-state Potts model [11, 12, 13, 14], and $z\approx 1$ for the 4D Ising model [15, 16]. Clearly, we would like to understand why this algorithm works so well in some cases and not in others. We hope in this way to obtain new insights into the dynamics of non-local Monte Carlo algorithms, with the ultimate aim of devising new and more efficient algorithms.

There is at present no adequate theory for predicting the dynamic critical behavior of an SW-type algorithm. However, there is one rigorous lower bound on $z$ due to Li and Sokal [11]. The autocorrelation times of the standard (multi-cluster) SW algorithm for the ferromagnetic $q$-state Potts model are bounded below by a multiple of the specific heat:

[TABLE]

Here ${\cal N}$ is the bond density in the SW algorithm, ${\cal E}$ is the energy, and $C_{H}$ is the specific heat; $\tau_{\rm int}$ and $\tau_{\rm exp}$ denote the integrated and exponential autocorrelation times, respectively [4, 6]. As a result one has

[TABLE]

where $\alpha$ and $\nu$ are the standard static critical exponents. Thus, the SW algorithm for the $q$-state Potts model cannot completely eliminate the critical slowing-down in any model in which the specific heat is divergent at criticality. The bound (1.1)/(1.2) has also been proven to hold [13] for the direct SW-type algorithm [17] for the Ashkin-Teller (AT) model [18, 19].

The important question is the following: Is the Li–Sokal bound (1.1)/(1.2) sharp or not? An affirmative answer would imply that we could use the bound to predict the dynamic critical exponent(s) $z$ given only the static critical exponents of the system. There are three possibilities:

1. i)

The bound (1.1) is sharp (i.e., the ratio $\tau/C_{H}$ is bounded), so that (1.2) is an equality. 2. ii)

The bound is sharp modulo a logarithm (i.e., $\tau/C_{H}\sim\log^{p}L$ for $p>0$). 3. iii)

The bound is not sharp (i.e., $\tau/C_{H}\sim L^{p}$ for $p>0$), so that (1.2) is a strict inequality.

Unfortunately, the empirical situation, even for the simplest cases, is far from clear. Let us review the status of this problem for the 2D Potts models. For the Ising case, the numerical data [7, 8, 9, 10] are consistent with a power-law behavior with $z_{{\rm int},{\cal E}}$ ranging from $0.35\pm 0.01$ [7] to $0.25\pm 0.01$ [9, 10]. However, the data are also consistent with a logarithmic behavior $z_{{\rm int},{\cal E}}=0\times\log$ [8]. In [13] we reanalyzed the high-precision data of Baillie and Coddington [9, 10]222 The exact value of the specific heat at finite $L$ was taken from the paper by Ferdinand and Fisher [20]. , and we found that the ratio $\tau_{\rm int,{\cal E}}/C_{H}$ behaves most likely as a pure power law $\sim L^{p}$ with $p=0.060\pm 0.004$ (statistical error only) or as a logarithm $\sim A+B\log L$. This means that the bound (1.1) is either non-sharp by a small power, or else it is sharp modulo a logarithm (in the latter case, the leading term would be $\tau_{{\rm int},{\cal E}}\sim\log^{2}L$). It is extremely difficult to distinguish between these two scenarios with lattice sizes up to only $L=512$.

The 3-state Potts model was first considered under this perspective in [11]. The dynamic critical exponent was found to be $z_{{\rm int},{\cal E}}=0.55\pm 0.03$, which is significantly larger than the exact result $\alpha/\nu={2\over 5}=0.4$ [21]. So, it seemed that the bound (1.2) was not sharp at all.

The 4-state Potts model is rather peculiar: the naive fit to the data, $z_{{\rm int},{\cal E}}=0.89\pm 0.05$ [11], is smaller than the (exactly known) value of $\alpha/\nu=1$ [22]. The explanation of this paradox is that the true leading term in the specific heat has a multiplicative logarithmic correction, $C_{H}\sim L\log^{-3/2}L$ [14, 23, 24, 25]. Indeed, a naive power-law fit to the specific heat yielded $\alpha/\nu=0.75\pm 0.01$, consistent with the bound (1.2). A high-precision study of this model was carried out in [12, 13].333 In [12, 13] the “embedding” version of an SW-type algorithm for the AT model was used. This algorithm reduces to the standard SW algorithm at the Ising subspace, but not at the 4-state Potts subspace. However, it was shown numerically that both algorithms belong (as expected) to the same dynamic universality class at the 4-state Potts subspace.

Naive power-law fits to the data showed that the bound (1.2) was satisfied: $z_{{\rm int},{\cal E}}=0.876\pm 0.012$ and $\alpha/\nu=0.768\pm 0.009$. On the other hand, the behavior of the ratio $\tau_{{\rm int},{\cal E}}/C_{H}$ was consistent only with two scenarios: a power law $\sim L^{p}$ with $p=0.118\pm 0.012$ (statistical error only), or a logarithmic growth $\sim A+B\log L$. This means that (1.1) fails to be sharp either by a small power or by a logarithm. In conclusion, there are two likely behaviors for the autocorrelation time: either $\tau_{{\rm int},{\cal E}}\sim L^{q}\log^{-3/2}L$ with $q\approx 1.12$, or else $\tau_{{\rm int},{\cal E}}\sim L\log^{-1/2}L$. In both cases, we find multiplicative logarithmic corrections to the autocorrelation time, which make the numerical analysis extremely difficult. These two scenarios for the ratio $\tau_{{\rm int},{\cal E}}/C_{H}$ coincide with those obtained for the Ising case.

In summary, for the 2D Potts models with $q=2$ and $q=4$, the Li–Sokal bound might be either sharp modulo a logarithm or else not sharp with a very small power $p\equiv z-\alpha/\nu$ (0.04\mathrel{\hbox to0.0pt{\lower 3.0pt\hbox{\mathchar 536\relax}\hss}\raise 2.0pt\hbox{\mathchar 316\relax}}p\mathrel{\hbox to0.0pt{\lower 3.0pt\hbox{\mathchar 536\relax}\hss}\raise 2.0pt\hbox{\mathchar 316\relax}}0.12). Moreover, if we interpolate between these models by following the self-dual (critical) curve of the AT model [18, 19], we obtain [12, 13] the same two scenarios: the ratio $\tau_{{\rm int},{\cal E}}/C_{H}$ either grows like a power law with a very small power $p$ (which increases slightly as we move from the Ising model to the 4-state Potts model along the AT self-dual curve), or else grows like a logarithm. Thus, there is some kind of continuity along the AT self-dual curve for the dynamic critical behavior of the SW algorithm.444 A similar study was carried out in [17] for the single-cluster version of the algorithm, but the lattices were not very large ($L\leq 256$).

There results have led us to reappraise the status of the Li–Sokal bound for the 2D 3-state Potts model. In [11] the bound was declared not sharp on the basis of numerical evidence suggesting that $p\equiv z-\alpha/\nu=0.15\pm 0.03$. However, this value is not much larger than that obtained [12, 13] for the 4-state Potts model; this suggests that the data for the 3-state model might also be consistent with a logarithm. Indeed, a closer look at the results of [11] concerning the 3-state model reveals that the lattices studied were not very large ($L\leq 256$), and the statistics were not very high (the number of iterations for the $L=256$ lattice was at most $25\times 10^{3}\tau_{{\rm int},{\cal E}}$). This motivated us to reconsider this model and extend the results of [11] to larger lattices, with higher statistics.

This paper is organized as follows: Section 2 reviews the basics of the Swendsen–Wang algorithm for the Potts models, as well as the proof of the Li–Sokal bound for these models. In Section 3 we describe our Monte Carlo simulations. In Section 4 we discuss in detail our methods of statistical data analysis. Finally, in Section 5 we present the analysis of our numerical (static and dynamic) data, culminating in a discussion of the sharpness of the Li–Sokal bound.

2 Basic set-up and notation

2.1 Potts model and Swendsen–Wang algorithm

The $q$-state Potts model assigns to each lattice site $i$ a spin variable $\sigma_{i}$ taking values in the set $\{1,2,\ldots,q\}$; these spins interact through the reduced Hamiltonian

[TABLE]

where the sum runs over all the nearest-neighbor pairs $\langle ij\rangle$. To simplify the notation we shall henceforth write $\delta_{\sigma_{i},\sigma_{j}}\equiv\delta_{\sigma_{b}}$ for a bond $b=\langle ij\rangle$. The ferromagnetic case corresponds to $\beta>0$. The partition function is defined as

[TABLE]

Finally, the Boltzmann weight of a configuration $\{\sigma\}$ is given by

[TABLE]

where $p=1-e^{-\beta}$.

The idea behind the Swendsen–Wang algorithm [4, 6, 26] is to decompose the Boltzmann weight by introducing new dynamical variables $n_{b}=0,1$ (living on the bonds of the lattice), and to simulate the joint model of old and new variables by alternately updating one set of variables conditional on the other set. The Boltzmann weight of the joint model is

[TABLE]

The marginal distribution of (2.4) with respect to the spin variables reproduces the Potts-model Boltzmann weight (2.3). The marginal distribution of (2.4) with respect to the bond variables is the Fortuin–Kasteleyn [27, 28, 29] random-cluster model with parameter $q$:

[TABLE]

where ${\cal C}(\{n\})$ is the number of connected components (including one-site components) in the graph whose edges are the bonds with $n_{b}=1$.

We can also consider the conditional probabilities of the joint distribution (2.4). The conditional distribution of the $\{n\}$ given the $\{\sigma\}$ is as follows: Independently for each bond $b=\langle ij\rangle$, one sets $n_{b}=0$ when $\sigma_{i}\neq\sigma_{j}$, and sets $n_{b}=0$ and 1 with probabilities $1-p$ and $p$ when $\sigma_{i}=\sigma_{j}$. Finally, the conditional distribution of the $\{\sigma\}$ given the $\{n\}$ is as follows: Independently for each connected cluster, one sets all the spins $\sigma_{i}$ in that cluster equal to the same value, chosen with uniform probability from the set $\{1,2,\ldots,q\}$.

The Swendsen–Wang algorithm simulates the joint probability distribution (2.4) by alternately applying the two conditional distributions just described. That is, we first erase the current $\{n\}$ configuration, and generate a new $\{n\}$ configuration from the conditional distribution given $\{\sigma\}$; we then erase the current $\{\sigma\}$ configuration, and generate a new $\{\sigma\}$ configuration from the conditional distribution given $\{n\}$.

2.2 Li–Sokal bound

To prove the Li–Sokal bound we first notice that the transition matrix $P_{SW}$ of the Swendsen–Wang algorithm can be written as a product

[TABLE]

where $P_{\rm bond}$ (the update of the bond variables) and $P_{\rm spin}$ (the update of the spin variables) are given by the conditional expectation operators $E(\,\cdot\,|\{\sigma\})$ and $E(\,\cdot\,|\{n\})$, respectively.

The strategy behind the proof is to compute explicitly the autocorrelation function at time lags 0 and 1 for a suitable observable ${\cal O}$ (which we will choose to be a “slow mode” of the algorithm). The unnormalized autocorrelation function is defined as

[TABLE]

(where the expectations are taken in equilibrium), and the normalized autocorrelation function as

[TABLE]

Then, using some general properties of reversible Markov chains, we will deduce lower bounds for the autocorrelation times $\tau_{{\rm int},{\cal O}}$ and $\tau_{{\rm exp},{\cal O}}$. These will in turn imply lower bounds on the dynamic critical exponents $z_{{\rm int},{\cal O}}$ and $z_{{\rm exp},{\cal O}}$.

For the observable ${\cal O}$, we shall use the bond occupation

[TABLE]

From the joint Boltzmann weight (2.4) it is easy to compute the following bond expectation values conditional on the spin configuration $\{\sigma\}$:

[TABLE]

From these equations it is easy to compute the mean values $\langle{\cal N}\rangle$ and $\langle{\cal N}^{2}\rangle$, and thus $C_{\cal NN}(0)\equiv\hbox{var}({\cal N})\equiv\langle{\cal N}^{2}\rangle-\langle{\cal N}\rangle^{2}$:

[TABLE]

where the energy is defined as

[TABLE]

The unnormalized autocorrelation function at time lag 1 is given by

[TABLE]

Now, $P_{\rm bond}{\cal N}$ is equal to

[TABLE]

Therefore,

[TABLE]

and

[TABLE]

where the energy density $E$ and the specific heat $C_{H}$ have been defined as

[TABLE]

and $V$ is the number of lattice sites.

The correlation functions of ${\cal N}$ under $P_{SW}\equiv P_{\rm bond}P_{\rm spin}$ are the same as under the positive-semidefinite self-adjoint operator $P^{\prime}_{SW}\equiv P_{\rm spin}P_{\rm bond}P_{\rm spin}$. This implies (see e.g., [4]) that we have a spectral representation

[TABLE]

with a positive measure $d\nu$. From this spectral representation we conclude that

[TABLE]

If we now recall the definitions of the integrated and exponential autocorrelation times

[TABLE]

we conclude from (2.16)/(2.20) that

[TABLE]

These are precisely the bounds (1.1). If we take into account the expected behavior close to the critical point of the specific heat and the autocorrelation times, we deduce immediately the bounds (1.2).

Similar bounds hold for the autocorrelation times of the energy ${\cal E}$. This can be seen from the fact that $P_{\rm bond}{\cal N}=p{\cal E}$, which implies that

[TABLE]

and hence

[TABLE]

This equation proves that

[TABLE]

and furthermore allows us to bound $\tau_{{\rm int},{\cal E}}$ above and below in terms of $\tau_{{\rm int},{\cal N}}$:

[TABLE]

If the critical slowing-down is not completely eliminated, we expect the factor $\rho_{\cal NN}(1)$ to approach 1 (from below) as $L\rightarrow\infty$. Moreover, irrespective of the presence or absence of critical slowing-down, we expect $\rho_{\cal NN}(1)$ to be bounded away from zero as $L\rightarrow\infty$. Modulo this very weak hypothesis, (2.28) implies the equality of the dynamic critical exponents for the energy and bond-occupation observables:

[TABLE]

Unfortunately, we do not know how to rigorously rule out “exotic” behaviors, in which $\rho_{\cal NN}(1)$ tends to zero as $L\rightarrow\infty$ and yet $\tau_{{\rm int},{\cal N}}$ diverges because the autocorrelation function has an extremely long tail.

Finally, we can define a new (“energy-like”) bond observable: the nearest-neighbor connectivity

[TABLE]

where $\gamma_{b}$ equals 1 if both ends of the bond $b$ belong to the same cluster, and 0 otherwise. More generally, the connectivity $\gamma_{ij}$ can be defined for an arbitrary pair $i,j$ of sites:

[TABLE]

The interest in ${\cal E}^{\prime}$ stems from the fact that the conditional expectation of $\delta_{\sigma_{b}}$ given the bond configuration $\{n\}$ is essentially $\gamma_{b}$:

[TABLE]

This implies the following relation between the energy and the connectivity densities:

[TABLE]

where $k$ is the coordination number of the lattice (i.e., the number of nearest neighbors of any given site), and the connectivity density $E^{\prime}$ is defined as

[TABLE]

Furthermore, (2.32) tells us that

[TABLE]

This implies that

[TABLE]

and hence

[TABLE]

Thus, the bounds (1.1)/(1.2) hold also for the autocorrelation times of ${\cal E}^{\prime}$. Furthermore, we can obtain bounds analogous to (2.27)–(2.29) for the observable ${\cal E}^{\prime}$:

[TABLE]

and

[TABLE]

Using again (2.28) and (2.26) we arrive at

[TABLE]

If $\rho_{\cal NN}(2)$ is bounded away from zero as $L\rightarrow\infty$, we can conclude that the dynamic critical exponents of ${\cal N}$, ${\cal E}$ and ${\cal E}^{\prime}$ are all equal:

[TABLE]

3 Description of the simulations

3.1 Observables to be measured

We have made simulations of the two-dimensional 3-state Potts model at criticality,

[TABLE]

on a periodic square lattice of linear size $L$.

We have measured five basic observables. Three of them have been already defined: the energy ${\cal E}$ (2.12), the bond occupation ${\cal N}$ (2.9), and the nearest-neighbor connectivity ${\cal E}^{\prime}$ (2.30). The other two are

[TABLE]

and

[TABLE]

where \mbox{\boldmath\sigma}_{x}\in\hbox{\specialroman R}^{q-1} is the Potts spin in the hypertetrahedral representation [for $q=3$ this means \mbox{\boldmath\sigma}_{x}=(\cos(2\pi\sigma_{x}/3),\sin(2\pi\sigma_{x}/3))], $V=L^{2}$ is the number of lattice sites, and $(x_{1},x_{2})$ are the Cartesian coordinates of point $x$. The observable ${\cal F}$ can be regarded as the square of the Fourier transform of the spin variable at the smallest allowed non-zero momenta [i.e., $(\pm 2\pi/L,0)$ and $(0,\pm 2\pi/L)$ for the square lattice]. It is normalized to be comparable to its zero-momentum analogue ${\cal M}^{2}$.

From these observables we compute the following expectation values: the energy density $E$ (2.17), the specific heat $C_{H}$ (2.18), the connectivity density $E^{\prime}$ (2.34), the bond density

[TABLE]

the magnetic susceptibility

[TABLE]

the correlation function at momentum $(2\pi/L,0)$

[TABLE]

and the second-moment correlation length

[TABLE]

This definition of the correlation length is not equal to the exponential correlation length (=1/mass gap); but it is expected that both correlation lengths scale in the same way as we approach the critical point.

Remark: As a check we have also computed the mean-square size of the clusters,

[TABLE]

where the sum is over all the clusters ${\cal C}$ of activated bonds (i.e., with $n_{b}=1$) and $\#({\cal C})$ is the number of sites of the cluster ${\cal C}$. Using the Fortuin–Kasteleyn identities [4, 27, 28, 29], it is not difficult to show that

[TABLE]

For each observable ${\cal O}$ discussed above we have measured its autocorrelation function $\rho_{\cal OO}(t)$ (2.8); and from this function we have estimated the corresponding integrated autocorrelation time $\tau_{{\rm int},{\cal O}}$ (2.21). In Section 4 we explain in detail how we derived estimates of the mean values and the error bars for both static and dynamic quantities.

3.2 Summary of the simulations

We have run our Monte Carlo program on lattices with $L$ ranging from 4 to 1024. In all cases the initial configuration was random, and we have discarded the first $10^{5}$ iterations to allow the system to reach equilibrium; this discard interval is always greater than $10^{3}\,\tau_{{\rm int},{\cal E}}$.555 We expect that $\tau_{{\rm int},{\cal E}}/\tau_{{\rm exp},{\cal E}}\approx 0.96$ and $\tau_{{\rm exp},{\cal E}}=\tau_{{\rm exp}}$ for this algorithm [13] (see Section 5.5). So the discard interval is always greater than $10^{3}\tau_{{\rm exp}}$, which is more than sufficient.

For $4\leq L\leq 256$ the total run length was approximately $10^{6}\,\tau_{{\rm int},{\cal E}}$; for $L=512$, it was $2.2\times 10^{5}\,\tau_{{\rm int},{\cal E}}$; and for $L=1024$, it was $6.7\times 10^{4}\,\tau_{{\rm int},{\cal E}}$. In all cases, the statistics are high enough to permit a high accuracy in our estimates of the static (error $\sim$ 0.1–0.5%) and dynamic (error $\sim$ 0.5–2%) quantities.

For $4\leq L\leq 128$ our data were obtained from a single long run. For $L=256$ we made two independent runs (with different random-number-generator seeds); for $L=512$ we made three independent runs; and for $L=1024$ we made four independent runs. In each case, we discarded the first $10^{5}$ iterations of each run. The individual runs are all of length \mathrel{\hbox to0.0pt{\lower 3.0pt\hbox{\mathchar 536\relax}\hss}\raise 2.0pt\hbox{\mathchar 318\relax}}10^{4}\tau_{{\rm int},{\cal E}}, which is long enough to allow a good determination of the dynamic quantities.

To test the program, we compared the static results on the $4\times 4$ lattice to the exact solution (obtained by enumerating all the possible configurations on this lattice). The agreement was excellent (see Table 1). In addition, for all lattice sizes we checked the relations (LABEL:eq2.11a), (2.33) and (3.9) between the mean values of static observables; the relations (2.25)/(2.36) among the autocorrelation functions $C_{\cal NN}$, $C_{\cal EE}$ and $C_{\cal E^{\prime}E^{\prime}}$; and the relation (2.16) between the autocorrelation function $\rho_{\cal NN}(1)$ and the static observables $E$ and $C_{H}$. In all cases the agreement was also good.

The CPU time required by our program is approximately 7.2 $L^{2}$ $\mu$s/iteration on an IBM RS-6000/370, and 3.6 $L^{2}$ $\mu$s/iteration on a single processor of an IBM RS/6000 SP2. The total CPU time used in the project was approximately 1.5 years on an IBM RS/6000 SP2. The smallest lattices were run on a IBM RS-6000/370 at NYU, while the largest ones ($L\geq 128$) were run on the IBM SP2 cluster at the Cornell Theory Center.

4 Statistical analysis of the Monte Carlo data

In this paper we are aiming at extremely high precision for both static and dynamic quantities; and furthermore we need to disentangle the effects of statistical errors from the effects of systematic errors due to corrections to scaling. For this, it is essential to obtain accurate estimates not only of the static and dynamic quantities of interest, but also of their error bars: in this way we will be able (see Section 5) to perform $\chi^{2}$ tests which provide an objective measure of the goodness of fit in each scaling Ansatz.

In this section we discuss in some detail how we performed the statistical analysis of our raw Monte Carlo data. In particular, we describe how to compute the estimators for the mean value and the variance of both static and dynamic quantities. These methods are based on well-known results of time-series analysis [30, 31], which we review briefly in Section 4.1.666 A review of time-series-analysis methods as applied to MC simulations can be found in [32, Appendix C] and in [4].

Then (Section 4.2) we describe an alternative analysis method, based on independent “bunches”, and report the results of detailed cross-checks that confirm (with one slight exception) the validity and reliability of the standard time-series-analysis method.

4.1 “Standard” time-series-analysis method

Let us consider a generic observable ${\cal O}$, whose mean is equal to $\mu_{\cal O}$. Its corresponding unnormalized and normalized autocorrelation functions are denoted by $C_{\cal OO}(t)\equiv\langle{\cal O}(0){\cal O}(t)\rangle-\langle{\cal O}\rangle^{2}$ and $\rho_{\cal OO}(t)\equiv C_{\cal OO}(t)/C_{\cal OO}(0)$, respectively. We also define the integrated autocorrelation time

[TABLE]

Given a sequence of $n$ Monte Carlo measurements of the observable ${\cal O}$ — call them $\{{\cal O}_{1},\ldots,{\cal O}_{n}\}$ — the natural estimator of the mean $\mu_{\cal O}$ is the sample mean

[TABLE]

This estimator is unbiased and has a variance

[TABLE]

This means that the variance is a factor $2\tau_{{\rm int},{\cal O}}$ larger than it would be if the measurements were uncorrelated. It is, therefore, very important to estimate the autocorrelation time $\tau_{{\rm int},{\cal O}}$ in order to ensure a correct determination of the error bar on the (static) quantity $\mu_{\cal O}$.

The natural estimator for the unnormalized autocorrelation function $C_{\cal OO}(t)$ is

[TABLE]

if the mean $\mu_{\cal O}$ is known, and

[TABLE]

if the mean $\mu_{\cal O}$ is unknown. We emphasize that, for each $t$, the estimators $\widehat{C}_{\cal OO}(t)$ and $\widehat{\widehat{C}}_{\cal OO}(t)$ are random variables [in contrast to $C_{\cal OO}(t)$, which is a number]. The estimator $\widehat{C}_{\cal OO}(t)$ is unbiased, and $\widehat{\widehat{C}}_{\cal OO}(t)$ is biased by terms of order $1/n$. The covariance matrices of $\widehat{C}_{\cal OO}$ and $\widehat{\widehat{C}}_{\cal OO}$ are the same to leading order in the large-$n$ limit (i.e., $n\gg\tau_{{\rm int},{\cal O}}$), and we have [30, 31]

[TABLE]

where $t,u\geq 0$ and $\kappa$ is the connected 4-point autocorrelation function

[TABLE]

The natural estimator for the normalized autocorrelation function $\rho_{\cal OO}(t)$ is

[TABLE]

if the mean $\mu_{\cal O}$ is known, and

[TABLE]

if the mean $\mu_{\cal O}$ is unknown. The estimators $\widehat{\rho}_{\cal OO}(t)$ and $\widehat{\widehat{\rho}}_{\cal OO}(t)$ are biased by terms of order $1/n$, as a result of the ratios of random variables in (4.8)/(4.9). The covariance matrices of $\widehat{\rho}_{\cal OO}$ and $\widehat{\widehat{\rho}}_{\cal OO}$ are the same to leading order in $1/n$. If the process is Gaussian, this covariance matrix is given in the large-$n$ limit by [31]

[TABLE]

If the process is not Gaussian, then there are additional terms proportional to the fourth cumulant $\kappa(m,t,t-u)$; these terms are, like those in (4.10), of order $1/n$. The simplest assumption is to consider the stochastic process to be “not too far from Gaussian”, and drop all the terms involving $\kappa$. If this assumption is not justified, then we are introducing a bias in the estimate of this covariance.

Finally, we shall take the estimator for the integrated autocorrelation time to be [32]

[TABLE]

[or the same thing with $\widehat{\widehat{\rho}}_{\cal OO}(t)$] where $M$ is a suitably chosen number. The reason behind the cutoff $M$ is the following: if we were to make the “obvious” choice $M=n+1$, then the resulting estimator would have a variance of order 1 even in the limit $n\rightarrow\infty$; this is because the terms $\widehat{\rho}_{\cal OO}(t)$ with large $t$ have errors (of order $1/n$) that do not vanish as $t$ grows [cf. (4.10)], and their number is also large ($\sim n$). Taking $M\ll n$ restores the good behavior of the estimator as $n\rightarrow\infty$. The bias introduced by this rectangular cutoff777 We could use more general cutoff functions, but this rectangular cutoff is the most convenient for the present purposes.

is given by

[TABLE]

The variance of the estimator $\widehat{\tau}_{{\rm int},{\cal O}}$ can be computed from the covariance (4.10); the final result is [32]

[TABLE]

where the approximation $\tau_{{\rm int},{\cal O}}\ll M\ll n$ has been made. A good (self-consistent) choice of $M$ is the following [32]: let $M$ be the smallest integer such that $M\geq c\widehat{\tau}_{{\rm int},{\cal O}}(M)$, where $c$ is a suitable constant. If the normalized autocorrelation function is roughly a pure exponential888 This has been found empirically to be true in the SW algorithm for the 2D Ising [33] and 4-state Potts models [13], and is confirmed here for the 3-state Potts model (see Section 5.5). , then the choice $c\approx 6$ is reasonable. Indeed, if we take $\rho_{\cal OO}(t)=e^{-t/\tau}$ and minimize the mean-square error

[TABLE]

using (4.12)/(4.13), we find that the optimal window width is

[TABLE]

For $n/\tau\approx 10^{6}$ (resp. $10^{4}$), we have $M_{\rm opt}/\tau\approx 6.56$ (resp. 4.26).

As noted above, we expect the estimator $\widehat{\tau}_{{\rm int},{\cal O}}$ to have a bias of order $\tau_{{\rm int},{\cal O}}/n$, due to the nonlinearities in (4.8)/(4.9).999 The bias on the estimator $\widehat{\tau}_{{\rm int},{\cal O}}$ also induces a bias on the estimated variance (4.3) of the sample mean $\overline{\cal O}$. This bias is of order $1/n^{2}$, i.e. a factor $1/n$ down from the variance (4.3) itself.

To make this bias negligible we need long runs. It has been shown empirically that this procedure works fairly well when n\mathrel{\hbox to0.0pt{\lower 3.0pt\hbox{\mathchar 536\relax}\hss}\raise 2.0pt\hbox{\mathchar 318\relax}}10^{4}\widehat{\tau}_{{\rm int},{\cal O}} [4].

Remarks: 1. The estimation of the error bar for the specific heat is a little bit more complicated. One can obtain var($C_{H}$) by computing var(${\cal E}$), var(${\cal E}^{2}$), and cov(${\cal E},{\cal E}^{2}$). This procedure has a numerical drawback: sometimes the covariance matrix for the observables ${\cal E}$ and ${\cal E}^{2}$ is nearly singular; then a small statistical fluctuation can cause the estimator of this matrix to be non-positive-definite. We are not aware of any procedure that ensures that the estimator of a covariance matrix is also positive-definite. To overcome this difficulty, we considered the observable ${\cal O}=({\cal E}-\mu_{\cal E})^{2}$, which can be studied using the standard method. As we do not know exactly the value of $\mu_{\cal E}$, we can use instead the sample mean $\overline{\cal E}$ (which should be computed first). To leading order in $1/n$ this procedure gives the right error bar for the specific heat.

2. We have a similar problem when computing the error bar for the second-moment correlation length $\xi$, which is a function of the two quantities $\chi$ and $F$ [cf. (3.6)]. In this case we considered the random variable

[TABLE]

which has automatically a zero mean. The error bar for the second-moment correlation length can be written easily as a function of the susceptibility $\chi$, the quantity $F$, and the variance of the above-mentioned observable ${\cal O}^{\prime}$. With this trick we take into account the cross-correlation between ${\cal M}^{2}$ and ${\cal F}$. This method needs the mean values $\mu_{{\cal M}^{2}}$ and $\mu_{\cal F}$; in practical situations they are substituted by the corresponding sample means $\overline{{\cal M}^{2}}$ and $\overline{\cal F}$ (which must be computed first).

This is the standard procedure we have used to analyze each of our Monte Carlo runs. Estimates coming from multiple independent runs (for $L\geq 256$) were merged using the standard formulae for statistically independent data. The results are reported in Table 1 (static quantities) and Table 2 (dynamic quantities). The results on $\tau_{{\rm int},{\cal E}}$ for $16\leq L\leq 256$ can be compared with those reported in [11]; the agreement is excellent ($\chi^{2}=6.17$, 5 degrees of freedom, confidence level = 56%).

4.2 Method based on independent bunches

Unfortunately, this standard method does not provide an easy way to compute the error bar for complicated quantities such as the ratio $\tau_{{\rm int},{\cal E}}/C_{H}$, which will play a central role in our analysis of the sharpness of the Li–Sokal bound (see Section 5.3). We can, of course, give an upper bound on the actual error bar by using the triangle inequality; but this upper bound will be a significant overestimate of the true value. If we were then to use this overestimate in fits, we would find artificially small values of $\chi^{2}$; as a result, the confidence levels would be artificially high, and useless for distinguishing good from bad fits. (At best we could distinguish better versus worse fits, by looking at the relative values of $\chi^{2}$.)

This fact motivated us to look for an alternative method to compute the error bars. There is also another advantage in having an alternate method: we can independently check the assumptions and approximations made in the standard procedure.

The second method works as follows. First, we split the whole sample of $n$ MC measurements $\{{\cal O}_{1},{\cal O}_{2},\ldots,{\cal O}_{n}\}$ into $m$ bunches of equal (or almost equal) length $\ell\equiv n/m$. For each bunch $i$ ($1\leq i\leq m$) we can compute the sample mean of the observable ${\cal O}$ and an estimate of its integrated autocorrelation time ($\overline{\cal O}_{i}$ and $\widehat{\tau}_{{\rm int},{\cal O},i}$, respectively). Indeed, we can also estimate the corresponding variances, but these estimates do not play any role in what follows. When computing the autocorrelation functions within each bunch, we used the whole sample mean $\overline{\cal O}$ as our estimate of $\mu_{\cal O}$ (instead of the bunch sample mean $\overline{\cal O}_{i}$); this trick reduces the bias in the estimates of the autocorrelation functions.

In this way we obtain two sequences of single-bunch estimates $\{\overline{\cal O}_{i}\}$ and $\{\widehat{\tau}_{{\rm int},{\cal O},i}\}$. If the bunches are long enough (i.e. $\ell\gg\tau_{{\rm int},{\cal O}}$), then the estimates from distinct bunches are almost statistically independent. Thus, we can define our estimates as follows:

[TABLE]

The quality of the results depends on the total number of measurements $n$ and the number $m$ of bunches we use. The merging of data coming from different runs is trivial in this case.

For $4\leq L\leq 256$ we have extremely good statistics ($n\sim 10^{6}\,\tau_{{\rm int},{\cal E}}$). This allows us to vary $m$ over at least one order of magnitude, and thereby to provide a cross-check on the standard time-series analysis. In this discussion it is useful to divide the observables into three categories: linear static, non-linear static, and dynamic. The first category includes ${\cal E}$, ${\cal E}^{\prime}$, ${\cal N}$, ${\cal M}^{2}$ and ${\cal F}$, whose sample mean values are linear in the raw MC data. The second category includes the specific heat and the second-moment correlation length, whose mean values are non-linear functions of the raw MC data. Finally, the third category contains all the autocorrelation times, which are also non-linear functions of the raw data.

For $4\leq L\leq 256$ we first divided the whole sample into $m=100$ bunches, each of them with a length $\ell\sim 10^{4}\,\tau_{{\rm int},{\cal E}}$. For $4\leq L\leq 64$ we repeated the analysis using $m=1000$ bunches, each of them with a length $\ell\sim 10^{3}\,\tau_{{\rm int},{\cal E}}$.

For the linear static observables, the “standard” and “bunch” methods always give identical mean values; this is a trivial identity, provided that the bunch lengths are exactly equal. As for the error bars on these observables, we find unsystematic discrepancies between the estimates given by the two methods: for $m=100$ the discrepancies are of order 10%, and for $m=1000$ they are of order 2%. In other words, the size of these discrepancies is roughly of order $\sim 1/\sqrt{m}$, and the sign is random; this is exactly what one expects on theoretical grounds for the statistical fluctuations in the estimators (LABEL:O_error_bunch) and (LABEL:tau_int_error_bunch).

For the non-linear static observables, the agreement between the mean values coming from the standard and the bunch methods is exact only in the case of the specific heat; this is because we have used the same estimator $\overline{\cal E}$ in both methods. For the correlation length, the mean values show small systematic discrepancies between the two methods, of order 0.05–0.1 standard deviations when $m=100$; the bunch method always gives a larger estimate than the standard method. In absolute value, these discrepancies range from $6\times 10^{-4}$ (for $L=4$) to $2\times 10^{-2}$ (for $L=256$). The same qualitative behavior is found when we repeat the analysis with $m=1000$ bunches. Now the differences in the mean value of the correlation length are 0.6–1.6 standard deviations. In absolute value, they range from $8\times 10^{-3}$ (for $L=4$) to $4\times 10^{-2}$ (for $L=64$), i.e. they are roughly one order of magnitude larger than in the $m=100$ case. Thus, the discrepancies in the mean value of the correlation length are systematic with a size of order $\sim m$; this is exactly what one expects for the mean value of a non-linear observable, which is afflicted by a bias of order $1/n$ in the standard method versus $1/\ell=m/n$ in the bunch method. Regarding the error bars, we find unsystematic discrepancies of order 10% when $m=100$, and of order 2% when $m=1000$. They behave as the error bars of the static linear observables.

For the autocorrelation-time estimators, we find systematic differences between the two methods: the estimates coming from the standard method are consistently smaller than those coming from the bunch method. For $m=100$, these differences are rather small compared to the statistical error bars (\mathrel{\hbox to0.0pt{\lower 3.0pt\hbox{\mathchar 536\relax}\hss}\raise 2.0pt\hbox{\mathchar 316\relax}}0.5 standard deviations). When $m=1000$ these discrepancies are much more relevant: their size is roughly two standard deviations. In absolute value, the discrepancies when $m=1000$ are one order of magnitude larger than when $m=100$. Thus, we find a systematic bias of size $\sim m$, again as expected for a quantity which exhibits a bias of order $1/n$ versus $1/\ell$. The error bars for the autocorrelation times are consistently larger for the standard method than for the bunch method, except for the specific heat where the behavior is consistently the opposite one. For $m=100$ bunches, these discrepancies are of order 15%; and for $m=1000$ the size remains at the same level. Thus, the discrepancies for the error bars of dynamic observables do not depend much on $m$; rather they are of order $\sim 1$.

From the above discussion, we conclude that for the linear static observables the two methods show excellent agreement, both for the mean values (trivial equality) and for the error bars (random discrepancies of order $1/\sqrt{m}$). The same holds for the error bars of the non-linear static observables. This confirms that the standard method is giving accurate estimates of the error bars, at least when the run length $n$ is large enough to provide a good determination of the autocorrelation time (which largely determines the static-quantity error bar). This $1/\sqrt{m}$ dependence also confirms our theoretical prediction that the standard method gives a more accurate estimate of the error bars than the bunch method. Indeed, the bunch method can be considered roughly equivalent to employing the standard method with the window width $M$ taken of order the bunch length $\ell$; but if $\ell\gg\tau_{{\rm int},{\cal O}}$ (as it must be in the bunch method), this is an unnecessarily large window width, and thus leads to unnecessarily large statistical fluctuations.

For the mean values of the non-linear observables (both static and dynamic), we likewise confirm the validity of the standard method. Once again the standard method is more reliable than the bunch method, but for a different reason: the bias of order $1/n$ is much smaller than the bias of order $1/\ell=m/n$. The latter bias becomes particularly serious when the number $m$ of bunches is large (as it must be in order to get good estimates of the error bars!).

Finally, the least understood piece is the determination of the error bar of the autocorrelation time: our data from the bunch method suggest that the standard method may be making a systematic error of order $15\%$. Perhaps this systematic error (if indeed it is real) arises from our neglect of the contributions of the fourth-order cumulant $\kappa$ to the covariance (4.10). This point definitely merits further investigation.

In summary, we think that the bunch method provides a good confirmation of the estimates given by the standard method. We shall hereafter consider the values given by the standard method to be the definitive ones, except for the ratio $\tau_{{\rm int},{\cal E}}/C_{H}$ where the standard method does not yield any correct error bar. In this latter case we shall use the central value coming from the standard method (which in fact agrees with the bunch-method value to within 0.1–1%; the discrepancies are of order 0.2 standard deviations), and the error bars coming from the bunch method with 100 bunches (for $4\leq L\leq 256$), 26 bunches (for $L=512$) or 55 bunches (for $L=1024$). We shall also compute the upper bound on the error bar coming from the standard method combined with the triangle inequality. These results for $\tau_{{\rm int},{\cal E}}/C_{H}$ are shown in Table 3.

5 Data analysis

For each quantity ${\cal O}$, we carry out a fit to the power-law Ansatz ${\cal O}=AL^{p}$ using the standard weighted least-squares method. As a precaution against corrections to scaling, we impose a lower cutoff $L\geq L_{min}$ on the data points admitted in the fit, and we study systematically the effects of varying $L_{min}$ on the estimates $A$ and $p$ and on the $\chi^{2}$ value. In general, our preferred fit corresponds to the smallest $L_{min}$ for which the goodness of fit is reasonable (e.g., the confidence level101010 “Confidence level” is the probability that $\chi^{2}$ would exceed the observed value, assuming that the underlying statistical model is correct. An unusually low confidence level (e.g., less than 5%) thus suggests that the underlying statistical model is incorrect — the most likely cause of which would be corrections to scaling.

is \mathrel{\hbox to0.0pt{\lower 3.0pt\hbox{\mathchar 536\relax}\hss}\raise 2.0pt\hbox{\mathchar 318\relax}} 10–20%), and for which subsequent increases in $L_{min}$ do not cause the $\chi^{2}$ to drop vastly more than one unit per degree of freedom.

Our final estimates for static and dynamic critical exponents are collected in Table 4.

5.1 Static quantities

There are a few exactly known results concerning the 2D 3-state Potts model. We know all the critical exponents [19, 21] and, in particular, the ratios

[TABLE]

which are the relevant quantities we can directly estimate from our Monte Carlo data. The leading correction-to-scaling exponent is also known [34, 35]:

[TABLE]

using the notation in which the correction term is $|\beta-\beta_{c}|^{\theta}$ or $L^{-\theta/\nu}$.

We can write the singular part of the free energy as a function of the thermal field $t=\beta-\beta_{c}$, the ordering field $h$, the leading irrelevant field $u$, and the linear size of the system $L$ as [36]

[TABLE]

where $d$ is the dimensionality of the lattice. If we differentiate (5.4) twice with respect to the thermal field $t$ and then take the limit $t=h=0$, we get the specific heat at criticality on a finite lattice:

[TABLE]

where $\alpha/\nu=2/\nu-d$, $G$ is the second derivative of $F$ with respect to its first argument, and the dots indicate subdominant corrections. Thus, the corrections to the specific heat at criticality are given by $L^{-\Delta}$ with $\Delta=\theta/\nu=4/5=0.8$. A similar analysis can be carried out for the magnetization and the susceptibility, giving again corrections proportional to $L^{-\Delta}$.

The energy $E$ is obtained by differentiating of the full free energy $f=f_{s}+f_{ns}$ with respect to the thermal field $t$. The contribution of the non-singular piece is believed to be trivial: there is numerical evidence that $f_{ns}(t,L)=f_{ns}(t,\infty)$ [36]. In other words, this contribution has no $L$-dependence, and gives merely the infinite-volume value of the energy at the given temperature, $E(\beta,\infty)$. The contribution of the singular piece can be obtained by differentiating (5.4) once with respect to $t$; it goes to zero like $L^{-d+1/\nu}$, which thus gives the leading correction to scaling for the energy. This correction is of order $L^{-4/5}$, which (by pure coincidence as far as we can tell) is exactly the same order as the correction $L^{-\Delta}$ for the divergent static observables. Finally, we note that the energy of the 2D $q=3$ Potts model at criticality is also exactly known [19] to be $E(\beta_{c},\infty)=1+1/\sqrt{3}\approx 1.577350$.

We can check these predictions by performing the fit $E-E(\beta_{c},\infty)=AL^{-w}$. The quality of the fit is very good already for $L_{min}=16$:

[TABLE]

with $\chi^{2}=3.23$ (5 DF, level = 66%). The agreement with the predicted exponent $2-1/\nu=4/5=0.8$ is truly spectacular. (Indeed, it is probably a coincidence that the agreement is so good. In general, one can’t expect to obtain anywhere near this accuracy for correction-to-scaling exponents.)

The fits of the susceptibility to a power law $AL^{\gamma/\nu}$ are quite stable. Our preferred fit corresponds to $L_{min}=64$:

[TABLE]

with $\chi^{2}=4.09$ (3 DF, confidence level = 25%). This result is 2.6 standard deviations away from the exact value $\gamma/\nu=1.73333$. This discrepancy could be due to corrections to scaling. We can try to fit our data to $AL^{26/15}(1+BL^{-\Delta})$ with various choices for $\Delta=1.1,1.0,\ldots,0.1$ as well as $0\times\log$ (i.e., a correction $1/\log L$). We find that the best fits correspond to $\Delta\approx 0.8$, in agreement with the theoretical prediction; for $\Delta=0.8$ and $L_{min}=8$ we obtain $\chi^{2}=6.55$ (6 DF, level = 36%). Surprisingly, in all these fits we find that the $\chi^{2}$ remains almost constant when $L_{min}$ is increased beyond our “preferred” value (and the confidence level consequently deteriorates): for example, with $L_{min}=256$ we get $\chi^{2}=3.37$ (1 DF, level = 7%) for the pure power-law behavior and $\chi^{2}=3.91$ (1 DF, level = 5%) for the fit with $\gamma/\nu=26/15$ and $\Delta=0.8$. This suggests that the point with $L=1024$ is off by about two standard deviations (possibly because the error bar is underestimated111111 This is quite possible: though the total run length at $L=1024$ is $7\times 10^{4}\tau_{{\rm int},{\cal E}}$, the individual runs (on which the time-series analysis was performed) ranged in length from only $10^{4}\tau_{{\rm int},{\cal E}}$ to $2.5\times 10^{4}\tau_{{\rm int},{\cal E}}$; and with such short runs the time-series analysis may not be completely reliable. ). As a matter of fact, if we drop this point, we obtain a good power-law fit for $L_{min}=128$:

[TABLE]

with $\chi^{2}=0.39$ (1 DF, level = 53%). This result agrees excellently with the exact value. If we impose the right power $\gamma/\nu=26/15$ and try to fit the first correction-to-scaling exponent $\Delta$, we again find that the best fits correspond to $\Delta\approx 0.8$. For $L_{min}=8$ we get $\chi^{2}=3.84$ (5 DF, level = 57%).

For the specific heat we find that the fits to power law $AL^{\alpha/\nu}$ are not stable at all: the confidence levels are horrible, and there is a clear trend towards smaller values of $\alpha/\nu$ as $L_{min}$ is increased. The least bad fit is obtained for $L_{min}=256$:

[TABLE]

with $\chi^{2}=3.80$ (1 DF, level = 5%). This value is eight standard deviations away from the exact result $\alpha/\nu=2/5=0.4$. Unlike the 4-state Potts model [23, 24], we do not expect multiplicative corrections to the leading term of the specific heat. We do, however, expect additive corrections to scaling of the form $AL^{2/5}(1+BL^{-\Delta})$. If we try the same exponents $\Delta$ as in the susceptibility, we find a decent fit for $\Delta\approx 0.6$: for $L_{min}=128$ we get $\chi^{2}=1.21$ (2 DF, level = 55%). This value of the exponent $\Delta$ is not far from the expected value 0.8, but the cause of the discrepancy is unknown; perhaps there is a large next-to-leading correction to scaling.

Finally, the second-moment correlation length $\xi$ is expected to behave linearly in $L$ as $L\rightarrow\infty$. In particular, the ratio $x\equiv\xi/L$ should approach a constant $x^{\star}$. We have tested this behavior. Already for $L_{min}=16$ our data are consistent with a constant value

[TABLE]

with $\chi^{2}=8.08$ (6 DF, level = 23%). The fact that a good fit can be obtained with such a small $L_{min}$ implies that the corrections to scaling are very small for this observable. [If we fit to $\xi(L)/L=x^{\star}+BL^{-\Delta}$ with 0.1\mathrel{\hbox to0.0pt{\lower 3.0pt\hbox{\mathchar 536\relax}\hss}\raise 2.0pt\hbox{\mathchar 316\relax}}\Delta\mathrel{\hbox to0.0pt{\lower 3.0pt\hbox{\mathchar 536\relax}\hss}\raise 2.0pt\hbox{\mathchar 316\relax}}1.5, we find a very slight improvement in the goodness of fit, and the estimated $x^{\star}$ decreases somewhat (by less than 0.0008 if \Delta\mathrel{\hbox to0.0pt{\lower 3.0pt\hbox{\mathchar 536\relax}\hss}\raise 2.0pt\hbox{\mathchar 318\relax}}0.6). But the estimated coefficient $B$ is consistent with zero at the $1.5\sigma$ level.] As in the case of the susceptibility, we also find that the goodness of fit deteriorates as $L_{min}$ is increased (e.g., for $L_{min}=512$ we have $\chi^{2}=4.80$, 1 DF, level = 3%) This might be due to the fact that the value for $L=1024$ is a little bit off (or its error bar is underestimated). If we drop this point, we obtain a good fit again for the same $L_{min}=16$:

[TABLE]

with $\chi^{2}=5.82$ (5 DF, level = 32%). But now the fit with $L_{min}=256$ is reasonable ($\chi^{2}=1.24$, 1 DF, level = 27%). We remark that the value $x^{\star}$ should in principle be calculable by conformal field theory; we hope that someone will perform this calculation and test our results (5.10)/(5.11).

5.2 Dynamic quantities

In this section we are going to fit the autocorrelation times for the observables ${\cal O}={\cal E}$, ${\cal E}^{\prime}$, ${\cal N}$ and ${\cal M}^{2}$ to a simple power law $\tau_{{\rm int},{\cal O}}=AL^{z_{{\rm int},{\cal O}}}$.

Let us start with the energy ${\cal E}$. The fit $\tau_{{\rm int},{\cal E}}=AL^{z_{{\rm int},{\cal E}}}$ is not very stable: the estimate of the power decreases systematically as $L_{min}$ is increased, and the $\chi^{2}$ is poor until $L_{min}=128$ where the estimate stabilizes within errors (this behavior suggests that there are strong corrections to scaling). Our preferred fit corresponds to $L_{min}=128$:

[TABLE]

with $\chi^{2}=0.44$ (2 DF, level = 80%). This value is greater than our estimate for $\alpha/\nu$; thus, the Li–Sokal bound (1.2) holds, though apparently not as a strict equality. To compare our result with the one reported in [11], we redid the fit using only the data with $L\leq 256$. The fit again shows a systematic decrease of the exponent as $L_{min}$ is increased, as well as a poor $\chi^{2}$, until $L_{min}=64$ where we get $z_{{\rm int},{\cal E}}=0.531\pm 0.005$ with $\chi^{2}=2.67$ (1 DF, level = 10%). This is consistent with the result $z_{{\rm int},{\cal E}}=0.55\pm 0.02$ reported in [11], but our error bar is one-fourth of theirs. The slightly higher estimate of $z_{{\rm int},{\cal E}}$ in [11] seems to arise from corrections to scaling induced by their choice $L_{min}=16$. Actually, the best fit to the data in [11] corresponds to $L_{min}=32$ and gives $z_{{\rm int},{\cal E}}=0.52\pm 0.02$ with $\chi^{2}=3.93$ (2 DF, level = 14%).

The fit of the autocorrelation time for the bond occupation ${\cal N}$ follows the same pattern: the estimate of the power decreases strongly as $L_{min}$ increases, and the $\chi^{2}$ is initially horrible; eventually the power stabilizes within errors, and the $\chi^{2}$ becomes reasonable. Our preferred fit is $L_{min}=128$:

[TABLE]

with $\chi^{2}=0.71$ (2 DF, level = 70%). The difference with respect to $z_{{\rm int},{\cal E}}$ is only 2.3 standard deviations, consistent with the theoretical prediction (2.29) that $z_{{\rm int},{\cal N}}=z_{{\rm int},{\cal E}}$. To check this result, we studied the ratio $\tau_{{\rm int},{\cal N}}/\tau_{{\rm int},{\cal E}}$. Since the standard time-series-analysis method would give only a upper bound on the right error bar for this quantity, we used instead the error bar provided by the “bunch” method. Among all the Ansätze we tried, only one gave a good fit for $L_{min}=128$: asymptotically constant with corrections $\sim L^{-1/4}$ ($\chi^{2}=0.75$, 2 DF, level = 69%). Thus, we conclude that in the SW algorithm the dynamic critical behavior of the energy and the bond density are the same: $z_{{\rm int},{\cal N}}=z_{{\rm int},{\cal E}}$.

In the same way, we considered the nearest-neighbor connectivity ${\cal E}^{\prime}$. The pattern of the fit is the same as above, and our preferred fit again corresponds to $L_{min}=128$:

[TABLE]

with $\chi^{2}=0.47$ (2 DF, level = 79%). This time, the agreement with $z_{{\rm int},{\cal E}}$ is extremely good, confirming the theoretical prediction (2.41) that $z_{{\rm int},{\cal E^{\prime}}}=z_{{\rm int},{\cal E}}$.

Finally, a similar behavior is observed in the fit for the autocorrelation time of the squared magnetization ${\cal M}^{2}$: the estimate of the power shows a clear trend towards smaller values as $L_{min}$ is increased, and the $\chi^{2}$ is initially poor. Our preferred fit is $L_{min}=128$, for which we get

[TABLE]

with $\chi^{2}=0.36$ (2 DF, level = 84%). This power is slightly smaller than $z_{{\rm int},{\cal E}}$ (the difference is seven standard deviations). It would be interesting to know whether this difference is significant or not. Let us consider the ratio $\tau_{{\rm int},{\cal M}^{2}}/\tau_{{\rm int},{\cal E}}$, again using the error bar provided by the “bunch” method. We tried to fit this ratio to various Ansätze, but only two gave good results for $L_{min}=32$: a pure power-law behavior $AL^{p}$ with $p=-0.0409\pm 0.0007$ ($\chi^{2}=0.57$, 5 DF, level = 97%); and a constant plus corrections of the type $\sim L^{-1/16}$ ($\chi^{2}=0.41$, 5 DF, level = 98%). We are therefore unable to resolve whether $z_{{\rm int},{\cal M}^{2}}$ is exactly equal to $z_{{\rm int},{\cal E}}$ or not; but if it is not equal, then it is only very slightly smaller ($z_{{\rm int},{\cal M}^{2}}-z_{{\rm int},{\cal E}}\approx-0.04$).

In conclusion, we have shown numerically that the energy ${\cal E}$, the bond occupation ${\cal N}$, and the nearest-neighbor connectivity ${\cal E^{\prime}}$ all have the same dynamic critical exponent $z_{\rm int}$, as was “almost proved” in Section 2. On the other hand, the observable ${\cal M}^{2}$ has a similar but perhaps not identical dynamic critical behavior: $z_{{\rm int},{\cal M}^{2}}$ may coincide with $z_{{\rm int},{\cal E}}$, or it may be slightly smaller.

5.3 Analysis of the Li–Sokal bound

In the previous subsection we have seen that all the observables considered have the same (or, in the case of ${\cal M}^{2}$, almost the same) dynamic critical exponent, and the common $z_{\rm int}$ is strictly larger than $\alpha/\nu$. This implies that the Li–Sokal bound (1.2) is not sharp. However, the difference $z_{\rm int}-\alpha/\nu\approx 0.115$ is actually not very large. There are a few arguments in favor of a more detailed analysis:

1. i)

The power-law fit to the specific heat was not very good: the estimated value of $\alpha/\nu$ decreased as $L_{min}$ increased, and we were unable to reach (within statistical errors) the exact value $\alpha/\nu=2/5$. If we compare the observed values of $z_{\rm int}$ and $\alpha/\nu$, we find that $z_{\rm int}-\alpha/\nu$ is only $\approx 0.09$. 2. ii)

The power-law fits to the autocorrelation times also exhibited this trend to smaller values as $L_{min}$ increased. Even though the fits seemed stable for $L_{min}=128$, this might well be an artifact due to the large error bars associated with the largest lattices ($L\geq 512$) compared to the smallest ones ($L\leq 256$). 3. iii)

In [12, 13] it was shown that differences $z_{\rm int}-\alpha/\nu$ of order 0.1 could actually be due to multiplicative logarithmic corrections. Indeed, such a multiplicative logarithmic correction was found to be a likely scenario even for models not having such a logarithmic correction in the specific heat.

Here, we will follow the approach of Refs. [12, 13], and consider the ratio $\tau_{{\rm int},{\cal E}}/C_{H}$; we shall use the error bars obtained from the “bunch” method (see Table 3).121212 In this section we only consider the energy ${\cal E}$, as the other observables have the same (or very slightly smaller) dynamic critical exponent.

We have tried to fit this ratio to different Ansätze:

A pure power-law behavior $AL^{p}$.

• 2)

A logarithmic growth, either as $A+B\log L$ or as $A\log^{p}L$.

• 3)

Asymptotically constant with additive corrections to scaling $A+BL^{-\Delta}$. We have considered the cases $\Delta=2,1,{1\over 2},{1\over 4},{1\over 8}$, and $0\times\log$ (i.e., $\tau_{{\rm int},{\cal E}}/C_{H}=A+B/\log L$).

Among all these Ansätze, only two were reasonably good. The best one is the simple power-law behavior $AL^{p}$. For $L_{min}=32$ we obtain

[TABLE]

with $\chi^{2}=1.72$ (4 DF, level = 79%). When $L_{min}>32$, the value of the power $p$ stays stable, and the $\chi^{2}$ decreases slowly and consistently.

The second-best fit corresponds to the logarithmic growth $A+B\log L$. For $L_{min}=64$, we get $\chi^{2}=1.93$ (3 DF, level = 59%). Again, the estimates are stable for $L_{min}>64$ and the $\chi^{2}$ is reasonable. However, this Ansatz seems to be slightly inferior to the power-law fit: both $L_{min}$ and the $\chi^{2}$ are greater than in the power-law case. To test the logarithmic Ansatz, we can impose the known critical behavior of the specific heat and perform the fit $\tau_{{\rm int},{\cal E}}=L^{2/5}(A+B\log L)$. A reasonably good result is obtained for $L_{min}=32$, giving $\chi^{2}=1.54$ (4 DF, level = 82%).

The asymptotically constant fits are always horrible, unless one take $L_{min}=256$ or larger. The only semi-exception is $\Delta={1\over 8}$, which gives a tolerable fit ($\chi^{2}=5.54$, 3 DF, level = 14%) already for $L_{min}=64$.

Of course for $L_{min}=256$ we obtain reasonably good fits for all Ansätze. But this is because the error bars for $L>256$ are so large that we are unable to distinguish between the various scenarios.

In summary, we have shown that there are only two likely scenarios for the ratio $\tau_{{\rm int},{\cal E}}/C_{H}$:

[TABLE]

or equivalently

[TABLE]

The first scenario implies that the Li–Sokal bound (1.2) is a strict inequality, with $p\equiv z_{{\rm int},{\cal E}}-\alpha/\nu\approx 0.08$–0.12; while the second scenario means that this bound is violated only by a logarithm (i.e., the bound is sharp modulo logarithms). Thus, the same “dynamic continuity” we found interpolating between the 2-state and the 4-state Potts models along the AT self-dual line is also found along the $q$-state Potts model critical line. Moreover, the power $p$ found is midway between the power found for the Ising model ($p\approx 0.05$) and the one found for the 4-state Potts model ($p\approx 0.12$) [13].

5.4 Further discussion of the Li–Sokal bound

The proof of the Li–Sokal bound (Section 2.2) is based on computing the autocorrelation function for the observable ${\cal N}$ at time lags 0 and 1 in terms of static observables, and then exploiting the inequality $\rho_{\cal NN}(t)\geq\rho_{\cal NN}(1)^{|t|}$. The apparent non-sharpness of this bound indicates that the large-$t$ behavior of $\rho_{\cal NN}(t)$ is not fully predicted by its behavior at $t=1$. Otherwise put, if we define the initial autocorrelation time

[TABLE]

then the Li–Sokal bound

[TABLE]

is apparently not sharp: that is, $\tau_{{\rm int},{\cal N}}$ and $\tau_{{\rm init},{\cal N}}$ diverge with different critical exponents $z_{{\rm int},{\cal N}}$ and $z_{{\rm init},{\cal N}}=\alpha/\nu$.

One might ask whether the situation could be improved by computing $C_{\cal NN}(t)$ exactly at (for example) time lag $t=2$ or $t=3$. Because of the identities (2.25)/(2.36), this would be equivalent to carrying out the Li–Sokal proof using ${\cal E}$ or ${\cal E}^{\prime}$ as the test observable in place of ${\cal N}$. At present we do not know how carry out this computation — the trouble is that we do not know how to express $\langle\gamma_{ij}\gamma_{kl}\rangle$ in terms of spin observables — but we can nevertheless test numerically whether we would obtain a better bound on $z_{{\rm int},{\cal E}}$ and/or $z_{{\rm int},{\cal E}^{\prime}}$ if we could carry out this computation.

We thus computed $\tau_{{\rm init},{\cal O}}$ for ${\cal O}={\cal N}$, ${\cal E}$, ${\cal E}^{\prime}$, with error bars given by (4.10)131313 Since the autocorrelation function for these observables is almost exactly a pure exponential, it suffices to perform the sums in (4.10) analytically and then insert the appropriate value of $\tau$ (see [13]). ; we then tried various fits, among others the fit $\tau_{{\rm init},{\cal O}}/C_{H}=AL^{r}$.141414 For the error bars on $\tau_{{\rm init},{\cal O}}/C_{H}$, we used for simplicity the triangle inequality rather than the (more correct) bunch method.

For ${\cal O}={\cal N}$ we have the identity

[TABLE]

where $E(\beta,\infty)$ is the infinite-volume energy at inverse temperature $\beta$; and we have taken into account the leading terms and the corrections to scaling discussed in Section 5.1. At criticality $\beta=\beta_{c}=\log(1+\sqrt{3})$ [hence $p=p_{c}=\sqrt{3}/(1+\sqrt{3})$] we know [19] the exact constant term:

[TABLE]

Indeed, if we fit our data to the Ansatz $\tau_{{\rm init},{\cal N}}/C_{H}-3/(1+\sqrt{3})=bL^{-s}$ we obtain a good fit for $L_{min}=16$, and the estimate is

[TABLE]

($\chi^{2}=1.92$, 3 DF, level = 86%), in excellent agreement with the theoretical prediction. Of course, we could try also the more primitive fit $\tau_{{\rm init},{\cal N}}/C_{H}=AL^{r}$, and we get a good fit for $L_{min}=128$:

[TABLE]

($\chi^{2}=0.90$, 2 DF, level = 64%), which is very close to the exact result $r=0$. All this is trivially to be expected, since our raw data for $\rho_{\cal NN}(1)$, $C_{H}$ and $E$ are in good agreement with the theoretical identity (2.16).

The analysis becomes non-trivial when we look at ${\cal O}={\cal E}$ and ${\cal E^{\prime}}$. For ${\cal O}={\cal E}$ we fit $\tau_{{\rm init},{\cal E}}/C_{H}=AL^{r}$, and get a good fit for $L_{min}=64$:

[TABLE]

($\chi^{2}=0.91$, 3 DF, level = 82%). For ${\cal O}={\cal E}^{\prime}$ the analogous fit is good for $L_{min}=32$:

[TABLE]

($\chi^{2}=0.51$, 4 DF, level = 97%). These exponents $r$ are a factor of $\approx 5$ smaller than the exponent estimates for $p\equiv z-\alpha/\nu$ found in Section 5.3. This suggests that $\tau_{{\rm init},{\cal E}}/C_{H}$ and $\tau_{{\rm init},{\cal E}^{\prime}}/C_{H}$ in fact tend to finite constants as $L\rightarrow\infty$ (as we know rigorously to be the case for $\tau_{{\rm init},{\cal N}}/C_{H}$), and that the apparent non-sharpness of the Li–Sokal bound cannot be remedied by using ${\cal E}$ or ${\cal E}^{\prime}$ as a test observable in place of ${\cal N}$. Rather, the non-sharpness of the bound arises from the fact that the long-time behavior of $\rho_{\cal NN}(t)$ is not sufficiently well predicted by its behavior at any small time.

5.5 Exponential autocorrelation time

The exponential autocorrelation time for an observable ${\cal O}$ is defined as151515 For a general Markov chain, the “$\lim$” should strictly speaking be replaced by “$\limsup$”, and $\rho_{\cal OO}(t)$ should be replaced by its absolute value. But here it can be proven that the limit really exists, and that $\rho_{\cal OO}(t)\geq 0$ for all $t$; this follows from the spectral representation (2.19) [or rather its analogue for ${\cal O}$].

[TABLE]

This autocorrelation time measures the decay rate of the “slowest mode” of the system, provided that this mode is not orthogonal to $\cal O$.

The critical behavior of $\tau_{{\rm exp},{\cal O}}$ is, in general, different from the behavior of $\tau_{{\rm int},{\cal O}}$. This fact can be seen from the standard dynamic finite-size-scaling Ansatz for the autocorrelation function $\rho_{\cal OO}(t)$:

[TABLE]

(Here the dependence on the coupling constants has been suppressed for notational simplicity.) Summing (5.28) over $t$, it follows that

[TABLE]

or equivalently,

[TABLE]

Thus, only when $p_{\cal O}=0$ do we have $z_{{\rm int},{\cal O}}=z_{{\rm exp},{\cal O}}$ [4]. In this latter case the Ansatz (5.28) can be rewritten in the equivalent form

[TABLE]

To test this latter Ansatz, we have plotted $\log\rho_{\cal OO}(t)$ versus $t/\tau_{{\rm int},{\cal O}}$ for the observables ${\cal O}=$ ${\cal N}$ (Figure 1), ${\cal E}$ (Figure 2) and ${\cal E}^{\prime}$ (Figure 3). On each figure we have plotted the data coming from different lattice sizes ($4\leq L\leq 128$) with different symbols; the error bars are computed from (4.10), using for simplicity the approximation [13] that the decay is a pure exponential (which is here almost exact). On each graph, we have also depicted for reference a line corresponding to the pure exponential $\rho_{\cal OO}(t)=e^{-t/\tau_{{\rm int},{\cal O}}}$. In these plots we have omitted the data for $L\geq 256$ for the sake of visual clarity; these data agree well with the curve found for smaller $L$, but with huge statistical error bars.

For ${\cal O}={\cal N}$ we see that the data coming from $16\leq L\leq 128$ collapse well onto a single curve. The lattices $L=4,8$ show slight systematic deviations from this limiting curve: these deviations are negative for t/\tau_{{\rm int},{\cal N}}\mathrel{\hbox to0.0pt{\lower 3.0pt\hbox{\mathchar 536\relax}\hss}\raise 2.0pt\hbox{\mathchar 316\relax}}1.5 and positive for t/\tau_{{\rm int},{\cal N}}\mathrel{\hbox to0.0pt{\lower 3.0pt\hbox{\mathchar 536\relax}\hss}\raise 2.0pt\hbox{\mathchar 318\relax}}1.5. This trend continues for $16\leq L\leq 128$, but the deviations are in most cases smaller than the error bars (especially for the larger lattices).

A similar behavior is found for ${\cal O}={\cal E}$ (Figure 2) and ${\cal O}={\cal E}^{\prime}$ (Figure 3), but the corrections to scaling are much weaker, and their sign is opposite from those seen for ${\cal N}$.

Thus, within statistical errors, we have found that the Ansatz (5.31) is satisfied. This implies that the integrated and exponential autocorrelation times for these three observables have exactly the same dynamic critical exponent, i.e. $z_{{\rm int},{\cal O}}=z_{{\rm exp},{\cal O}}$ for ${\cal O}={\cal N},{\cal E},{\cal E^{\prime}}$. This equality does not hold as a general rule in the theory of dynamic critical phenomena [4, 6], but it does appear to hold for algorithms of Swendsen–Wang type.

Finally, we find that $\rho_{\cal NN}(t)$ differs slightly but noticeably from the pure exponential $e^{-t/\tau_{{\rm int},{\cal N}}}$. The discrepancy from a pure exponential becomes smaller when we consider the other two observables ${\cal E}$ and ${\cal E}^{\prime}$. This is to be expected theoretically from the relations ${\cal E}\sim P_{\rm bond}{\cal N}$ and ${\cal E^{\prime}}\sim P_{\rm spin}{\cal E}$ [cf. (2.14)/(2.35)]: each action of $P_{\rm bond}$ or $P_{\rm spin}$ helps to “purify” the slowest mode, so that the autocorrelation function becomes closer to a pure exponential. On the other hand, the identities (2.26)/(2.37) imply161616 Using the fact that $\rho_{\cal NN}(1)\to 1$ and $\tau_{{\rm int},{\cal N}}\to\infty$ as $L\to\infty$, as follows from the Li–Sokal identity (2.16) and the fact that the specific heat $C_{H}$ is divergent.

that the limiting (scaling) functions $\widehat{h}_{\cal N}$, $\widehat{h}_{\cal E}$ and $\widehat{h}_{\cal E^{\prime}}$ are identical (assuming they exist at all); this means that at least some of the curves in Figures 1–3 have not yet reached their scaling limit.

Finally, a crude fit suggests that $\tau_{{\rm int},{\cal E}}/\tau_{{\rm exp},{\cal E}}\approx 0.96$, in agreement with the idea that $\rho_{\cal EE}(t)$ is almost but not quite a pure exponential.

Acknowledgments

We wish to give our big thanks to Wolfhard Janke for pointing out an error in an earlier version of this paper. We also thank the Cornell Theory Center for the computer time needed to complete this project. The authors’ research was supported in part by U.S. National Science Foundation grants DMS-9200719 and PHY-9520978, and by Metacenter grant MCA94P032P.