{'fraction_non_alphanumeric': 0.05105356281164549, 'fraction_numerical': 0.020620878237091845, 'mean_word_length': 4.046428571428572, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 16, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dynamic critical behavior of cluster algorithms for 2D Ashkin-Teller and Potts models

We study the dynamic critical behavior of two algorithms: the Swendsen-Wang algorithm for the two-dimensional Potts model with q=2,3,4 and a Swendsen-Wang-type algorithm for the two-dimensional symmetric Ashkin-Teller model on the self-dual curve. We find that the Li--Sokal bound on the autocorrelation time \tau_{{\rm int},{\cal E}} \geq const \times C_H is almost, but not quite sharp. The ratio \tau_{{\rm int},{\cal E}}/C_H appears to tend to infinity either as a logarithm or as a small power (0.05 \ltapprox p \ltapprox 0.12). We also show that the exponential autocorrelation time \tau_{{\rm exp},{\cal E}} is proportional to the integrated autocorrelation time \tau_{{\rm int},{\cal E}}.


Introduction
Monte Carlo (MC) simulations [4,5,25] 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 limitated by critical slowing-down: the autocorrelation time τ , which roughly measures the MC time needed to produce a "statistically independent" configuration, diverges near a critical point. More precisely, for a finite system of linear size L at criticality, we expect a behavior τ ∼ L z for large L. The power z is a dynamic critical exponent and it depends on both the system and the algorithm.
In single-site MC algorithms (such as single-site Metropolis or heat bath), critical slowing-down arises fundamentally from the fact that the updates are local; and have a dynamic critical exponent z ∼ > 2. This is a severe critical slowing-down: the total amount of computer work needed to study a d-dimensional lattice of linear size L grows L 2 times faster than the naive geometrical factor L d . To study static and dynamic critical behavior we need high-precision data (run lengths ∼ < 10 4 τ ), so in practice, it is very difficult to obtain high-precision data for large lattices with local algorithms. The geometrical factor L d is unavoidable for the usual MC simulations, so the elimination (or the reduction) of the critical slowing down is the only way to make MC simulations feasible close to a critical point.
In some cases, a much better dynamic behavior is obtained by allowing nonlocal moves, such as cluster flips. In particular, the Swendsen-Wang (SW) algorithm for the q-state ferromagnetic Potts model [26] achieves a significant reduction in z compared to the local algorithms: numerical experiments suggest that 0.22 ∼ < z ∼ < 1, where the exact value depends on the dimensionality of the lattice d and on q (See Table 1). The most favorable case is the two-dimensional (2D) Ising model (q = 2), for which the numerical data suggest z = 0.2186 ± 0.0068 [24] but may also be compatible with a logarithmic growth [12,1,2]. In other cases, the performance of the SW algorithm is less impressive (though still quite good): e.g., z = 0.515 ± 0.006 for the 2D 3-state Potts model [22], z = 0.876 ± 0.011 for the 2D 4-state Potts model [21,23], and z ≈ 1 for the 4D Ising model [16,20]. 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 MC algorithms, with the ultimate aim of devising new and more efficient algorithms.
where α and ν 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 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: i) The bound (1.1) is sharp (i.e., the ratio τ /C H is bounded), so that (1.2) is an equality.
ii) The bound is sharp modulo a logarithm (i.e., τ /C H ∼ log p L for p > 0).
iii) The bound is not sharp (i.e., τ /C H ∼ 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. From Table 1 we see that the Li-Sokal bound is apparently not sharp for the Ising model in d > 2. However, the differences z int,E −α/ν are much smaller in d = 2. A first look at Table 1 reveals that the Li-Sokal bound for the Ising model could only be sharp if the autocorrelation time grows like a logarithm [12,1,2]. Otherwise, the bound (1.2) would be non-sharp. The Li-Sokal bound for the 3-state Potts model seems to be non-sharp: the difference z int,E − α/ν = 0.115 ± 0.006 is clearly not consistent with zero. Finally, the 4-state Potts model is rather peculiar: the naive fit to the data, z int,E = 0.876 ± 0.011 [23], is smaller than the (exactly known) value of α/ν = 1. The explanation of this paradox is that the true leading term in the specific heat has a multiplicative logarithmic correction, C H ∼ L log −3/2 L [19,7,6,23]. Indeed, a naive power-law fit to the specific heat yields α/ν = 0.770 ± 0.008 [23], consistent with the bound (1.2). The Li-Sokal bound would be sharp modulo a logarithm if τ ∼ L log p L with p ≥ −3/2. Thus, in d = 2, q = 2 and q = 4 are candidates for a sharp (perhaps modulo a logarithm) Li-Sokal bound; but for q = 3 this bound is apparently nonsharp, at least if the numerical estimates are to be trusted. It is important, however, to proceed cautiously in interpreting these estimates, which are after all only fits to data from a limited range of L values (typically L ∼ < 1024) and which may be significantly biased by corrections to scaling. In particular, it is extremely difficult to distinguish numerically between a logarithm and a small power. 1 The above empirical situation motivated in part the study of the 2D Ashkin-Teller (AT) model using a SW-type algorithm [21]. This model interpolates between the Ising and the 4-state Potts models. The study of the AT model proved to be useful to obtain new insights on the sharpness of the Li-Sokal bound.
In Section 2 we will define more rigorously the concept of autocorrelation time. In Section 3 we will review the SW algorithm for the q-state ferromagnetic Potts model and the Li-Sokal bound for such algorithm. In Section 4 we will introduce two different SW-type algorithms for the AT model. In Section 5 we will analyze the sharpness of the Li-Sokal bound for both cases. Finally, in Section 6 we will consider the dynamic critical behavior of the exponential autocorrelation time.

Autocorrelation times
The goal of a MC simulation is to generate random samples of the spins σ = {σ i } distributed according to a probability distribution Π. This measure is usually Gibbsian Π = e −H(σ) /Z, where H is the reduced Hamiltonian of the system and Z is the partition function (which is generally unknown).
We then simulate the Markov process defined by the transition matrix P , starting at some arbitrary initial configuration σ (0) . In this way we generate a random sequence of configurations σ (1) , σ (2) , . . . σ (t) , . . .. Given the two properties above (e.g., ergodicity and stationarity of Π), it is guaranteed that this Markov chain converges to the equilibrium distribution Π, irrespective of the initial configuration.
In every MC simulation two different autocorrelation times can be defined. Given an arbitrary observable of the spins A(σ), we can define its equilibrium normalized autocorrelation function: Typically, one expects that this function decays exponentially (e.g., ρ AA (t) ∼ e −|t|/τ ). This motivates the following definition of the exponential autocorrelation time for the observable A τ exp,A = lim sup The autocorrelation time τ exp,A measures the relaxation time of the slowest mode coupled to the observable A. The exponential autocorrelation time of the system is the relaxation time of the slowest mode of the system This number can be infinite for some perfectly legal algorithms! [25]. The integrated autocorrelation time of the observable A is also defined in terms of ρ AA : There are two fundamental (and distinct) issues in dynamical MC simulations: 1) Initialization bias. If the initial configuration σ (0) is not "characteristic" of the equilibrium probability distribution Π, the first MC configurations are not distributed according to Π; they will introduce a systematic error in our MC estimates. To eliminate this bias, we have to discard the first N 0 iterations such that when t > N 0 the system has reached the equilibrium distribution Π. But how large should N 0 be? One answer is given by the exponential autocorrelation time τ exp (2.3). This autocorrelation time is a worst-case bound on relaxation to equilibrium. For most systems, it suffices to discard the first N 0 ∼ > 20τ exp iterations. In this case, the deviation from equilibrium (in the l 2 sense) will be at most e −20 (≈ 2 × 10 −9 ) times the initial deviation from equilibrium.
2) Autocorrelation in equilibrium. Even if the system is in equilibrium, the MC configurations are not statistically independent, but highly correlated. The MC estimate of the mean value A is given byĀ where N is the total number of measurements. The variance of such estimator can be written as when N ≫ τ int,A . Thus, var(Ā) is 2τ int,A times larger than in independent sampling. This means we need a good determination of the integrated autocorrelation time (2.4) in order to obtain a reliable estimate of var(Ā). In practice, the condition N ≫ τ int,A is replaced by N ∼ > 10 4 τ int,A [18].
Close to a second-order phase transition both autocorrelation times (2.2)/(2.4) diverge as a power law 2 The dynamic critical exponents z exp,A and z int,A are in general different [25] (See Section 6).

Swendsen-Wang algorithm for the Potts model
The SW algorithm [26] is defined for the ferromagnetic q-state Potts model. This model assigns to each lattice site i a spin variable σ i taking values in the set {1, 2, . . . , q}. These spins interact through the reduced Hamiltonian where the sum runs over all the nearest-neighbor pairs ij , and the all the coupling constants are positive J ij ≥ 0 ∀ ij . The Boltzmann weight of a configuration {σ} is given by where p ij = 1 − e −J ij , and the partition function is given by Z = {σ} e −H Potts . The idea behind the SW algorithm [26,9,25] is to decompose the Boltzmann weight by introducing new dynamical variables n ij = 0, 1 (living on the bonds of the lattice), and to simulate the joint model of spin and bond variables by alternately updating one set of variables conditional on the other set. The Boltzmann weight of the joint model is The marginal distribution of (3.3) with respect to the spin variables reproduces the Potts-model Boltzmann weight (3.2). The marginal distribution of (3.3) with respect to the bond variables is the Fortuin-Kasteleyn [15,11] random-cluster model with parameter q: where C({n}) is the number of connected components (including one-site components) in the graph whose edges are the bonds with n ij = 1. We can also consider the conditional probabilities of the joint distribution (3.3). The conditional distribution P bond = E(·|{σ}) of the {n} given the {σ} is as follows: independently for each bond ij , one sets n ij = 0 when σ i = σ j , and sets n ij = 0 and 1 with probabilities 1 − p ij and p ij when σ i = σ j . The conditional distribution P spin = E(·|{n}) of the {σ} given the {n} is as follows: independently for each connected cluster, one sets all the spins σ i in that cluster equal to the same value, chosen with uniform probability from the set {1, 2, . . . , q}.
The SW algorithm simulates the joint probability distribution (3.3) by alternately applying the two conditional distributions just described (i.e., P SW = P bond P spin ).
Step 1. Generate a new bond configuration {n} from P bond .
Step 2. Generate a new spin configuration {σ} from P spin .
The performance of this algorithm has been discussed in Section 1 (See also Table 1). Let us now briefly review the proof of the Li-Sokal bound (1.1)/(1.2) for the homogeneous Potts model p ij = p [17]. The strategy of the proof is simple: first, we choose two "slow" observables: the bond and energy densities and compute their autocorrelation functions at time lags 0 and 1. Then, using some general properties of reversible Markov chains, we will deduce lower bounds for the autocorrelation times τ int,A (for A = N , E) and τ exp . These will in turn imply lower bounds on the dynamic critical exponents z int,A and z exp . From (3.3) we can compute the expectation values of bond variables n ij conditional to the spin variables. In particular, we can compute the autocorrelation function of the bond occupation N at time lag 1 [17,22] where the energy density E and the specific heat C H are defined by: At criticality, p → p c > 0 and E → E c > 0, so the constant A does not vanish. The correlation functions of N under P SW = P bond P spin are the same as under the positive-semidefinite self-adjoint operator P ′ SW = P spin P bond P spin . This implies [17] that we have a spectral representation with a positive measure dν. It follows that If we introduce the inequalities (3.6)/ (3.9) in the definition of the autocorrelation times (2.2)/(2.4) for the bond occupation we obtain the corresponding Li-Sokal bounds (1.1) [17]. Similar bounds for the energy E are obtained by using the identity E(N |{σ}) = pE [17,22]. This equation implies the following relation between the autocorrelation functions of the bond occupation and energy: This implies in particular, the following set of equalities and inequalities [22]: From (3.11) the Li-Sokal bounds (1.1) for the energy follow (and also for the exponential autocorrelation time τ exp ). We can also deduce from Eq. (3.11a) that if ρ N N (1) = 0 as L → ∞, then [22].
4 Swendsen-Wang-like algorithm for the Ashkin-Teller model The AT model is a generalization of the Ising model to a four-state model. To each site i of the lattice we assign two Ising spins σ i = ±1 and τ i = ±1; they interact through the Hamiltonian It can be regarded as two Ising models (with nearest-neighbor couplings J and J ′ ) interacting via a four-spin coupling K. This model enjoys several symmetries: 1) It is symmetric under permutations of (σ, τ, στ ). This implies that (4.1) is symmetric under permutations of the couplings (J, J ′ , K). 2) On any bipartite lattice, (4.1) is symmetric under sublattice flip of two of the spins (σ, τ, στ ). This means that (4.1) is symmetric under sign flip of two of the couplings (J, J ′ , K).
In this paper we are mainly interested in one particular case of the Hamiltonian (4.1): the symmetric Ashkin-Teller (SAT) model characterized by J = J ′ . Although we do not know how to solve this model analyticly, we have a fairly good understanding of its phase diagram on a square lattice [3] (See Figure 1). The line K = 0 (dotted line in Figure 1) corresponds to two decoupled Ising models; its ferromagnetic critical point is denoted by DIs. The line J = K corresponds to the ferromagnetic 4-state Potts model, and its critical point is denoted by P (J = K = 1 4 log 3). The line J = −K corresponds to the 4-state Potts antiferromagnet, which in known to be non-critical at all temperatures T ≥ 0 [10]. The 4-state Potts subspace is depicted with dash-dotted lines in Figure 1. The line J = 0 corresponds to an Ising model in the variable στ ; its ferromagnetic critical point is denoted by Is, and its antiferromagnetic counterpart by AFIs. Finally, when K = ∞ we have another Ising model with σ = τ , and coupling 2J. Its ferromagnetic critical point is denoted by Is'. There are several critical curves in the SAT model. The most important one is the self-dual curve e −2K = sinh 2J. This curve is critical only for K ≤ 1 4 log 3 (solid curve in Figure 1), and is noncritical for K > 1 4 log 3 (dashed curve in Figure 1). The critical part belongs to the universality classes of the conformal field theories with conformal charge c = 1. Along this line the critical exponents vary continuously and they are known exactly [21, and references therein]. At the point P two more critical curves emerge: one goes to the critical point Is and the other one goes to the critical point Is' at K = ∞. Finally, there is another critical curve emerging from the critical point AFIs and pointing towards K = −∞. The exact location of these three curves is unknown, although there are good numerical approximations [14]. These lines are believed to belong to the Ising universality class. We have considered two different cluster algorithms for the AT model. In all cases we can assume that J, J ′ ≥ |K| (4.2) This condition is always achievable via the symmetries of the AT Hamiltonian (4.1) for a bipartite lattice.
Direct algorithm: The idea behind this algorithm is the same as in the standard SW algorithm for the Potts model: we decompose the Boltzmann weight by introducing new dynamical variables, and we then simulate the joint model by alternately updating one set of variables conditional on the other set. As we have two distinct sets of Ising spins, we expect to introduce two distinct sets of auxiliary variables: m ij , n ij = 0, 1 (associated to the spins σ, τ respectively). The Boltzmann weight of the joint model is Indeed, the Boltzmann weight of the original AT model is the marginal distribution of (4.3) with respect to the spin variables. From (4.3) we can read off the conditional probabilities needed in the SW procedure (See Ref. [21] for more details): Step 1. Update the bond variables {m, n} using the conditional distribution P bond = E(·|{σ, τ }).
This algorithm reduces on the ferromagnetic 4-state Potts subspace J = K > 0 to the standard SW algorithm. Wiseman and Domany [29] have introduced essentially this same decomposition of the Boltzmann weight using a different derivation as ours. They then studied numerically the single-cluster version of this algorithm. Here we study the many-cluster ("SW") version.
Embedding algorithm: The direct algorithm is perfectly legal, but it is somewhat complicated to write the computer code for its Step 1 in an efficient way. This motivated the introduction of a variant algorithm in which we deal with only one kind of spin (σ or τ ) at a time. Consider the Boltzmann weight of a given bond ij , conditional on the {τ } configuration (i.e., the τ spins are kept fixed): it is We can simulate this system of σ spins using a standard SW algorithm. The effective nearest-neighbor coupling is no longer translation-invariant, but this does not matter. The key point is that the effective coupling is always ferromagnetic, due to the condition (4.2). An exactly analogous argument applies to the {τ } spins when the {σ} spins are held fixed. The embedding algorithm for the AT model has therefore two parts: Step 1: Given the {τ } configuration (which we hold fixed), we perform a standard SW iteration on the σ spins. The probability p ij arising in the SW algorithm takes the value p ij = 1 − exp[−2(J + Kτ i τ j )].
Step 2: Given the {σ} configuration (which we hold fixed), we perform a standard SW iteration on the τ spins. In this case, Wiseman and Domany [29] also constructed an embedding version of their singlecluster algorithm. Furthermore, they showed that, in the single-cluster context, the direct and embedding algorithms define the same dynamics 3 ; only the computer implementation is different. However, this equivalence does not hold for our manycluster algorithm. In the direct algorithm we have independent clusters of σ spins and τ spins that could be flipped simultaneously. In the embedding algorithm we have at each step only one of the two types of clusters.
The embedding algorithm, due to its simplicity, is the one used in our MC study of the square-lattice SAT model. We expect that it lies in the same dynamic universality class as the direct algorithm, on the grounds that one SW hit of {σ} followed by one SW hit of {τ } should be roughly equivalent to one joint hit of {σ, τ }. Of course, we do not expect the autocorrelation times for the two algorithms to be equal, but we do expect them to be asymptotically proportional as the critical point is approached.
The embedding algorithm does not reduce to the standard SW algorithm on the 4-state Potts subspace. However, we expect that they do belong to the same dynamic universality class. This fact has been numerically verified by comparing the dynamic data for the 4-state Potts model with the embedding algorithm to those quoted in Ref. [17] (which correspond to the standard SW algorithm). We see that the ratio of the two autocorrelation times is more or less constant within statistical errors, and that there is no systematic trend as L grows. Thus, we conclude that they are proportional in the limit L → ∞ [21]: From Table 2 we see that the performance of the embedding algorithm is reasonably good on the self-dual curve of the SAT model compared to the local algorithms. The direct algorithm for the AT model satisfies the Li-Sokal bound (1.1). The proof is a straightforward generalization of the one given in Section 3 for the Potts case [21, Appendix A]. As we expect that both the direct and the embedding algorithms belong to the same universality class, we expect that there is a Li-Sokal bound also for the embedding algorithm. In Table 2 we show the dynamic critical exponents and the values of the ratio α/ν for several points on the self-dual curve of the SAT model. We have studied two points on the SAT self-dual curve between the critical Ising and 4-state Potts models; they are denoted by X2 and ZF (See Figure 1). X2 = (J ≈ 0.344132, K ≈ 0.147920) is one of the points considered in [29], and ZF = Model z int,E α/ν X2 0.477 ± 0.028 [21] ≈ 0.4183 ZF 0.733 ± 0.014 [21] 2/3 ≈ 0.666667 Table 2: Numerical estimates of the dynamical critical exponent z int,E for two points on the self-dual curve of the square-lattice symmetric Ashkin-Teller model. These points interpolate between the Ising and 4-state Potts critical models (See Figure 1). We also show the values of the ratio α/ν.
(J ≈ 0.302923, K ≈ 0.220343) corresponds to a model studied by Zamolodchikov and Fateev [30]. In both cases the Li-Sokal bound is satisfied and the difference z int,E − α/ν is not very large.

Testing the sharpness of the Li-Sokal bound
The most straightforward (and naive) method to study the sharpness of the Li-Sokal bound (1.1)/(1.2) is to fit the autocorrelation time to a power-law τ int = AL z int and compare the dynamic critical exponent z int to the ratio α/ν. This procedure has a clear disadvantage: in some models the leading term of the specific heat is not a pure power law (e.g., in the 2D q = 2, 4 Potts models). And there is no reason why this should not happen for the autocorrelation times as well. This motivates the search for a different criterium: The Li-Sokal bound (1.1) tells us that the ratio τ int,E /C H ≥ const., so we may try to fit it to different Ansätze. 4 In this section we will consider only the energy autocorrelation times, as this is one of the slowest modes of the system. We have considered three scenarios: 1. The Li-Sokal bound is sharp: z int,E = α/ν. In this case we expect that the ratio τ int,E /C H will converge to a constant as L → ∞. In particular, we have tried three different Ansätze compatible with this scenario The second (resp. third) Ansatz corresponds to power-law (resp. additive logarithmic) corrections to scaling.
2. The Li-Sokal bound is sharp modulo a logarithm: z int,E = α/ν ×log p with p > 0. In this case, the ratio τ int,E /C H will diverge like log p L. We have considered two Ansätze compatible with this scenario 3. The Li-Sokal bound is not sharp: z int,E = α/ν + p with p > 0. In this case the ratio τ int,E /C H will diverge like a power law. This scenario is represented by the Ansatz: Our numerical results [21,22,23,24] show that there are only two Ansätze that are likely to the six models we have considered here (See Tables 3-4). These are Thus, the Li-Sokal bound seems to be either sharp modulo a logarithm or non-sharp by a small power 0.05 ∼ < p ∼ < 0.12.  Table 3: We show the fits to the ratio τ int,E /C H for the different models we consider in this paper. For each Ansatz, we give the value of L min , the value of the χ 2 , the number of degrees of freedom (DF ), and the confidence level. For the power-law Ansatz AL p we also give the estimate of the power p. Note that for the 4-state Potts, X2 and ZF models, the error bar of the ratio τ int,E /C H was computed using the triangle inequality. Thus, the corresponding χ 2 values are somewhat overestimated.

Model
In the Ising model we have performed a high-precision MC simulation on lattices ranging from L = 4 to L = 512. In all cases we measured 10 6 τ int,E iterations [24]. In  Table 4: We show the fits to the integrated autocorrelation times for the different models we consider in this paper. For each Ansatz, we give the value of L min , the value of the χ 2 , the number of degrees of freedom (DF ), and the confidence level.
For the power-law Ansätze we also give the estimates of the powers (p or z int,E ).
order to increase the statistics we merged our data to that of Baillie and Coddington [1,2]. From Table 3 we see that the sharp Ansatz (5.1) is clearly disfavored: with a larger L min it achieves a poorer confidence level. On the other hand, there is no significant difference between the other two Ansätze (5.2)/(5.3). This can be seen clearly in Figure 2(a). If we look at Table 4, we see that the pure-power-law Ansatz is disfavored, as it is the logarithmic Ansatz proposed by several authors [12,1,2]. Both of them need a larger value of L min than in the other two Ansätze. Finally, the non-sharp-by-a-power Ansatz (5.3) is now slightly favored over the non-sharpby-a-logarithm Ansatz (5.2), but the difference is again probably not significant. In Figure 2(b) we see clearly that both fits are clearly favored over the simple power-law and the logarithmic growth.
In the 3-state Potts model we have performed a MC simulation on lattices ranging from L = 4 to L = 1024. In all cases we measured 10 6 τ int,E iterations [22]. From Table 3 we see that the power-law Ansatz (5.3) to the ratio τ int,E /C H is slightly favored over the logarithmic Ansatz (5.2). However, if we look at Table 4 we observe that the Ansatz L 2/5 (A log L + B) is favored over the non-sharp-by-a-power Ansatz: the former achieves a similar value of the confidence level with a smaller value of L min . In this model the finite-size corrections in the specific heat are very important [22] and this is possibly the cause that one scenario is favored when we study τ int,E , and a different one is favored when we consider the ratio τ int,E /C H .
In the 4-state Potts model we performed a MC simulation using the embedding algorithm for the SAT model on lattices ranging from L = 16 to L = 1024. In all cases (but L = 1024), we measured at least 10 4 τ int,E iterations; for L = 1024 we could only measure 1500τ int,E iterations [23]. If we look at Table 3 we see that the power-law fit is favored over the logarithmic fit: for the same value of L min the latter gets twice as much χ 2 than the former. However, if we look at the fits of τ int,E (Table 4), we conclude that the non-sharp-by-a-logarithm Ansatz is favored over the other two (i.e., the non-sharp-by-a-power and the pure power-law scenarios). With a smaller value of L min it achieves the best confidence level of all three.
Finally, the models X2 and ZF were studied using the embedding algorithm for SAT model. The simulations were performed on lattices ranging from L = 16 to L = 512. We performed at least 10 4 τ int,E measurements for the ZF model, and 3 × 10 4 τ int,E for the X2 model. In Tables 3-4 we see almost no difference between the fits corresponding to the non-sharp-by-a-power 5.3 and the non-sharp-by-a-logarithm 5.2 Ansätze.
In all cases we conclude that the Li-Sokal bound (1.1)/(1.2) is not sharp by a small quantity: either a logarithm or a small power p. Furthermore, there is some kind of continuity among the values of the power p: it grows from q ≈ 0.05 in the Ising model to q ≈ 0.12 in the 4-state model, irrespective if we go through the 3-state Potts model (p ≈ 0.08) or through the SAT self-dual curve (p(X2) ≈ 0.05 and p(ZF) ≈ 0.08). It is extremely difficult to distinguish between these two scenarios with lattices up to L = 1024. We need larger lattices L ≫ 1024 in order to solve this question.
It is intriguing why the Li-Sokal bound is so close to be sharp in two dimensions, and clearly non-sharp for d > 2. The explanation is unknown; but a reasonable guess is that in d > 2 the energy and the bond occupation are not among the slowest modes of the system (as in d = 2). It would be very interesting to characterize such slow modes for d > 2.

Exponential autocorrelation times
So far we have analyzed the behavior of the integrated autocorrelation time. However, it is not guaranteed that this behavior coincides with that of the exponential autocorrelation time. In general, we expect that z int,A = z int,A for most observables A [25]. This issue can be explained using the scaling relation for the autocorrelation function (2.1) where p A is an exponent, h A is a scaling function, and ξ is the correlation length of the system. Summing over the time t, we obtain where h A is another scaling function. This Ansatz can be numerically verified by plotting ρ EE as a function of t/τ int,E and testing if the data coming from different values of L collapse well onto a single curve. If so, then p E = 0.
We have plotted ρ EE versus t/τ int,E for the Ising model in Figure 3. Similar plots for the 3-and 4-state Potts models can be found in Refs. [22,21]. In all cases we find that ρ EE is close to a pure exponential and that the Ansatz (6.3) is satisfied within statistical errors (usually we find small corrections to scaling for the smaller values of L; e.g., in the Ising case the Ansatz (6.3) is satisfied for L ≥ 64). Thus, the conclusion is that the exponential and integrated autocorrelation times for the energy are likely to have the same critical exponent z int,E = z exp,E (6.4) for the SW algorithm (Potts model) and for the embedding algorithm (SAT model).
The value of τ exp,E can in principle be obtained by fitting ρ EE to an exponential. However, this is not an easy task as the MC estimates of ρ EE for different t are highly correlated [21]. One should compute the full covariance matrix for these random variables. In Ref. [21] a method was proposed to perform this computation. However, high-precision data is needed in order to obtain meaningful results: we need at least 6 × 10 4 τ int,E measurements. For all the models considered here we can nevertheless obtain crude estimates of the ratio τ int,E /τ exp,E , which are slightly smaller than 1.