Sweeny dynamics for the random-cluster model with small $Q$

The Hamiltonian of the $Q$-state Potts model [1] reads $\mathcal{H}=-J\sum_{\langle ij\rangle}\delta_{\sigma_{i},\sigma_{j}}$, where $J>0$ is the ferromagnetic coupling constant, $\langle ij\rangle$ stands for a pair of nearest-neighboring sites, and each lattice site $i$ has a spin $\sigma_{i}$, taking one of the $Q$ states $\sigma_{i}=1,2,...,Q$. For $Q=2$, the Potts model corresponds to the famous Ising model. Given a spin configuration, under the Fortuin-Kasteleyn transformation [2, 3], bonds between nearest-neighboring sites with the same spin value are independently occupied with the probability $p=1-e^{-J/kT}$ ($k$ is the Boltzmann constant, $T$ is the temperature, and, for convenience, we set $kT=1$), and the resulting graph is a configuration of the random-cluster model [4]. For bond configurations $A\subseteq{E}$ of a given graph $G=(V,E)$, the random-cluster model is defined by the partition function

Here $A$ is the set of “occupied bonds”, $|A|$ is the number of occupied bonds, $k(A)$ is the number of clusters (connected components), real positive numbers $Q$ and $v=p/(1-p)$ are the cluster and bond weights, respectively. For integer $Q\geq 2$, the random-cluster model serves as a graphical representation of the $Q$-state Potts model. As $Q\to 1$, the random-cluster model reduces to independent bond percolation [5].

In equilibrium, the weight ratio between two configurations are $W(\beta)/W(\alpha)=Q^{\delta k(A)}v^{\delta|A|}$, where $\delta k(A)$ ($\delta|A|$) represents the changes of $k(A)$ ($|A|$) from configuration $\alpha$ to $\beta$. As $Q\rightarrow 0$, when $v/Q$ is finite, a configuration consists of a spanning forest [6, 4]. This is because in this case $W(\beta)/W(\alpha)\sim Q^{\delta[k(A)+|A|]}$, and spanning forests have the minimum $k(A)+|A|$. From a spanning forest $\alpha$, occupying any internal-perimeter bond (defined later in Fig. 2) leads to a nonforest $\beta$ with $\delta[k(A)+|A|]=1$, such that $W(\beta)/W(\alpha)\sim Q\rightarrow 0$. In the case $Q\rightarrow 0$ with $v/Q^{\sigma}(0<\sigma<1)$ being finite, similarly the model consists of uniform spanning trees which have the smallest $k(A)+\sigma|A|$ [6, 4].

Markov Chain Monte Carlo (MCMC) algorithms [7, 8] are effective in simulating the random-cluster model, such as the Swendsen-Wang-Chayes-Machta algorithm [9, 10] for $Q>1$ with non-local cluster updates, and the Sweeny algorithm [11] for $Q>0$ with local bond updates. To assess the efficiency of the MCMC algorithms, autocorrelation times are commonly used as a metric [12]. For any observable ${\cal O}$, one first defines the autocorrelation function

Here expectations are taken after equilibrium is attained, and the time unit is one “hit”, in which only one bond is visited (a Monte Carlo sweep consists of $|E|$ hits, where $|E|$ is the number of all bonds). Then one defines the exponential autocorrelation time for an observable $\cal O$

and the overall exponential autocorrelation time

The integrated autocorrelation time for an observable $\cal O$ is defined as

Observables typically have $\tau_{{\rm exp},{\cal O}}=\tau_{\rm exp}$ and $\tau_{{\rm int},\cal O}\leq\tau_{\rm exp}$. The autocorrelation times $\tau_{\rm exp}$ and $\tau_{{\rm int},\cal O}$ play different roles in MCMC simulations [12]: $\tau_{{\rm exp}}$ is the relaxation time of the slowest mode in the system, and $\tau_{{\rm int},\cal O}$ controls the statistical error in the measurement of $\langle\cal O\rangle$. Thus $\tau_{{\rm int},\cal O}$ mainly determines the efficiency of MCMC sampling. For a time series of total length $M$, the sample mean

has variance

The above Eq. (7) can be used to measure $\tau_{{\rm int},{\cal O}}$.

In MCMC simulations, the dynamic processes usually undergo critical slowing-down [13]: the autocorrelation time diverges as $\tau/{|E|}\sim\xi^{z}$ when the critical point is being approached, where $\xi$ is the spatial correlation length and $z$ is called the dynamic exponent. For two-dimensional Ising model, Metropolis-Hastings [14, 15] and Glauber [16, 8] algorithms have $z\approx 2$.  To overcome the $critical$ $slowing$-$down$, the Swendsen-Wang-Chayes-Machta algorithm [9, 10] was proposed for the random-cluster model with $Q>1$. The latter algorithm has $z\approx 0.2$ for the two-dimensional Ising model, which represents a huge improvement. The gained efficiency is associated with nonlocal cluster updates.

However, some local algorithms may still be very efficient, such as the worm algorithm [17] and the Sweeny algorithm [18]. The worm algorithm outperforms the cluster algorithms for the three-dimensional Ising model in measuring quantities such as the susceptibility and correlation length. With local single-bond updates, the Sweeny algorithm can surprisingly exhibit critical speeding-up for a suitable range of $Q$ [18]: the relaxation time of some observables decreases significantly at time scales much shorter than one sweep as the lattice size $L$ increases, and this phenomenon becomes even more pronounced as $Q$ becomes smaller. We illustrate the latter property of the Sweeny algorithm by observing the dynamical behavior of two quantities, i.e., the number of occupied bonds $\mathcal{N}=|A|$ and the second moment of cluster sizes $\mathcal{S}_{2}=\sum{|\mathcal{C}|^{2}}$. Hereafter, if not specified, all simulations were performed on the square lattice at the critical point $v_{c}(Q)=\sqrt{Q}$. The average $\langle\mathcal{S}_{2}\rangle/|V|$ is a susceptibility-like quantity of the random-cluster model [2, 3], with $|V|$ being the volume. For $Q\geq 2$, it is found [18] that the dynamic exponent $z_{exp}\approx z_{{\rm int},\mathcal{N}}$, and that both $z_{{\rm int},\mathcal{N}}$ and $z_{{\rm int},\mathcal{S}_{2}}$ are only slightly larger than the lower bound $\alpha/\nu$ [19], and could possibly be equal to it (here $\alpha$ and $\nu$ are critical exponents for the specific heat and correlation length, respectively). For $Q<2$, the autocorrelation function of $\mathcal{N}$ follows an almost purely exponential function, as demonstrated for $Q=0.005$ and $Q=0.5$ in Fig. 1(a) and  1(c), respectively. However, the dynamic behavior of $\mathcal{S}_{2}$ is richer [18]: (1) as shown in Fig. 1(b) and 1(d), the autocorrelation function $\rho_{\mathcal{S}_{2}}(t)$ decays very rapidly in much less than one sweep, as a prelude to the final exponential decay $\sim e^{-t/\tau_{exp}}$; (2) the exponent $z_{{\rm int},\mathcal{S}_{2}}$ is a negative number, and it decreases with the increasing of lattice size $L$; (3) from Fig. 1(d) to 1(b), it is seen that the critical speeding-up is more pronounced for smaller $Q$.

In this article, we are interested in the dynamic behavior of the Sweeny algorithm for very small $Q$ values. A $Q$-dependent critical slowing-down is found for the random-cluster model in two dimensions. The phenomenon can be observed in the dynamics of a specific local pattern on the square lattice, i.e., a unit square with configuration or , where solids lines represent occupied bonds. We call the local pattern as a half face. We find that, when $Q\rightarrow 0$, the integrated autocorrelation time of the number of half faces in terms of sweeps (i.e., $\tau_{{\rm int},\mathcal{N}_{\rm HF}}/2L^{2}$, and hereafter $\tau_{\mathcal{N}_{\rm HF}}/L^{2}$ is used for brevity) diverges as $Q^{-1/2}$. This indicates that the Sweeny algorithm becomes not ergodic as $Q\to 0$.

The above $Q$-dependent critical slowing-down can be explained by the acceptance probability of updates in the Sweeny algorithm. As illustrated in Fig. 2, considering different occupation and connectivity, the bonds can be classified into four types — the $backbone$, $bridge$, $internal$-$perimeter$ and $external$-$perimeter$ bonds. At $v_{c}=\sqrt{Q}$, when transforming an internal-perimeter bond into a backbone bond, one has $\delta k(A)=0$ and $\delta|A|=1$, which yields the weight ratio $W(\beta)/W(\alpha)=\sqrt{Q}$; and transforming a bridge bond into an external-perimeter bond leads to $\delta k(A)=1$ and $\delta|A|=-1$, which also gives $W(\beta)/W(\alpha)=\sqrt{Q}$. It is shown in Sec. II that these weight ratios are associated with vanishing acceptance probabilities $~{}\sqrt{Q}$ as $Q\rightarrow 0$. To overcome the $Q$-dependent slowing-down, we first use the knowledge of four kinds of bonds in the detail balance condition and propose an improved Sweeny algorithm. The improved Sweeny algorithm largely alleviates the slowing-down but does not fully solve the problem. For example, at $Q=0$ where the graphs become spanning trees, the improved Sweeny algorithm is still invalid due to zero probability of deleting or adding a bond. We then employ the Kawasaki algorithm [20], which proposes the random exchange of an empty and an occupied bond. At $Q=0$ on spanning trees, the Kawasaki algorithm still has an $L$-dependent divergence $\tau_{\mathcal{N}_{\rm HF}}/L^{2}\sim L^{3/2}$ which can also be related to the acceptance probability of updates. The latter divergence can be eliminated through an improved Kawasaki algorithm which uses the categories of bonds. Finally, we combine the use of the improved Sweeny and improved Kawasaki algorithms, and find no critical slowing-down for small $Q$, i.e., the autocorrelation time $\tau/L^{2}$ is of order ${\cal O}(1)$, being independent of $Q$.

The remainder of this paper is organized as follows. Section II describes the Sweeny algorithm and the $Q$-dependent slowing-down. Section III presents an improved Sweeny algorithm which can greatly improve simulation efficiency. Section IV introduces an improvement of the Kawasaki algorithm and its combination with the improved Sweeny algorithm, which can completely eliminate the $Q$-dependent slowing-down. Section V contains a discussion and conclusion.

As mentioned in the Introduction, the Sweeny algorithm is a local single-bond update algorithm for the random-cluster model. Actually, it is the only known MCMC algorithm for the model with $0<Q<1$. The original Sweeny algorithm [11] is presented in Alg. 1 . A basic step of the Sweeny algorithm is proposing to change the status of a randomly selected bond, and accepting the proposal according to the Metropolis criterion.

For an efficient implementation of the Sweeny algorithm, the main challenge is in checking the connectivity between the two endpoints of a chosen bond. There are three typical methods with different asymptotic run-time scaling: breadth-first search, union-and-find, and dynamic connectivity algorithms [21, 22]. We use the method of simultaneous breadth-first searches starting at both endpoints of a bond $xy$, since it has very simple code and is effective enough for small and medium size $L$ [23]. We note that a simple modification of Alg. 1 can save some time by not requiring connectivity checking for a fraction of updates. The modification for $Q\leq 1$ at $v_{c}=\sqrt{Q}$ is shown in Alg. 2. The basic idea is to consider together several probabilities in Alg. 1: if the generated random number is less than the minimum of these probabilities, the update must be accepted, and if it is greater than the maximum of these probabilities, the update must be rejected. For $Q\to 1$, this small change in the algorithm can lead to a significant improvement of efficiency.

As described in the Introduction, the Sweeny algorithm exhibits rich dynamics, including both critical slowing-down and critical speeding-up. However, its dynamics at small $Q$ values are not uncovered fully. For Alg. 1 and Alg. 2, one can observe that they encounter a problem when $Q$ is small: at $v_{c}=\sqrt{Q}$, probabilities of some updates (occupying an internal-perimeter bond and deleting a bridge bond) are proportional to $\sqrt{Q}$, which causes the algorithm being not ergodic at $Q\to 0$. The problem was not severe in previous studies, since the studied $Q$ values are not extremely small, and conventional observables such as $\mathcal{N}$ and $\mathcal{S}_{2}$ could not detect any slowing-down effect due to the $\sqrt{Q}$-form updating probabilities. Actually on the square lattice in the limit $Q\to 0$, a critical equilibrium configuration consists of a spanning tree, thus the above observables remain fixed as $\mathcal{N}=L^{2}-1$ and $\mathcal{S}_{2}=L^{4}$ for a given $L$.

We find that the above problem can become dominant in sampling other quantities, such as the number of half faces $\mathcal{N}_{\rm HF}$. Using the original Sweeny algorithm, we calculated $\tau_{\mathcal{N}_{\rm HF}}$ as plotted in Fig. 3. It can be seen that $\tau_{\mathcal{N}_{\rm HF}}/L^{2}$ diverges as $Q^{-1/2}$, independent of the system size $L$. We refer to this as $Q$-dependent critical slowing-down. This newly discovered phenomenon at small $Q$ enriches the dynamics of the Sweeny algorithm. It is important to note that, away from $Q=0$, the original Sweeny algorithm is not inherently flawed for small $Q$. For example, when $Q=10^{-6}$, the autocorrelation time $\tau_{\mathcal{N}_{\rm HF}}/L^{2}$ is of order $O(10^{2})$, which means that $\mathcal{N}_{\rm HF}$ can still be measured rather efficiently.

To understand the $Q$-dependent critical slowing-down, we associate the above factor $Q^{-1/2}$ to the updating probabilities in the Sweeny algorithm. In the following we measure the numbers of the aforementioned four kinds of bonds, and then give these updating probabilities.

On the square lattice in the limit $Q\to 0$, the number of bridge bonds ($\mathcal{N}_{\rm BR}$) and of internal-perimeter bonds ($\mathcal{N}_{\rm IP}$) are $\mathcal{N}_{\rm BR}\to L^{2}-1$ and $\mathcal{N}_{\rm IP}\to|E|-\mathcal{N}_{\rm BR}$, respectively. These numbers can be obtained from the characteristics of a spanning tree: since every node (i.e., site) corresponds to an occupied bond except the root, there are $L^{2}-1$ occupied bonds for a tree with $L^{2}$ nodes; the $L^{2}-1$ occupied bonds are all bridge bonds, and other (empty) bonds are internal-perimeter bonds. For finite small $Q$, changes of the numbers of four kinds of bonds with the decrease of $Q$ are shown in Fig. 4, with fitting results for two of them in Table 1. For $L=16$, we find that for small $Q$, the average number of backbone bonds ($\mathcal{N}_{\rm BB}$) and external-perimeter bonds ($\mathcal{N}_{\rm EP}$) are respectively

both being proportional to $\sqrt{Q}$, as plotted in Fig. 4(a). The average number of occupied bonds is found to be $\langle\mathcal{N}_{B}\rangle\simeq L^{2}-1$ for small $Q$, which is consistent with $\mathcal{N}_{\rm BR}\to L^{2}-1$ in the limit $Q\to 0$. In the simulations, for convenience of a loop construction which will be described in Section III, we use periodic boundary conditions in one direction of the square lattice and open boundary conditions in the other direction. The total number of bonds on the lattice is $|E|=2L^{2}-L$. Thus, considering $\mathcal{N}_{B}=\mathcal{N}_{\rm BB}+\mathcal{N}_{\rm BR}$ and $|E|-\mathcal{N}_{B}=\mathcal{N}_{\rm EP}+\mathcal{N}_{\rm IP}$, we derive

The above scaling behaviors are demonstrated by a plot of $1-\langle\mathcal{N}_{\rm BR(IP)}\rangle/L^{2}$ as shown in Fig. 4(b). In the thermodynamic limit $L\rightarrow\infty$, from the self-duality of the square lattice, the critical numbers of four kinds of bonds can be obtained [24] (an alternative derivation is given in the Appendix of the present paper) as

In the limit $Q\rightarrow 0$, the factor $0.94$ in Eqs. (8) and (9) deviates from $1$ in Eq. (10) due to finite-size effects.

With the above bond numbers for small $Q$, we can get the updating probabilities in the Sweeny algorithm as following: the probability of selecting a bridge bond or internal-perimeter bond is $\sim 1$, and the acceptance probability ($Q/v$ or $v$) is $\sqrt{Q}$, then the associated updating probability is $\sim\sqrt{Q}$; the probability of selecting a backbone bond or external-perimeter bond is $\sim\sqrt{Q}$, and the acceptance probability is 1, then the associated updating probability is also $\sim\sqrt{Q}$. Thus the $Q^{-1/2}$-form critical slowing-down is associated with the $\sqrt{Q}$-form updating probabilities.

To overcome the $Q$-dependent critical slowing-down, we introduce an improved Sweeny algorithm as shown in Alg. 3. Instead of selecting a bond at random, the improved algorithm randomly proposes one of the following four operations: deleting a backbone bond (delete-backbone), deleting a bridge bond (delete-bridge), adding a backbone bond (add-backbone), and adding a bridge bond (add-bridge). Acceptance probabilities of the operations shall be given later. The algorithm consists of recording and dynamically updating the information of four kinds of bonds, which are directly used in proposing MCMC updates. To the end of this section, we shall see that the improved Sweeny algorithm can largely solve the $Q$-dependent critical slowing-down.

For Alg. 3, a list is created to dynamically keep the information for each kind of bonds. The lists record the indices of corresponding bonds, and the numbers of four types of bonds, i.e., $\mathcal{N}_{\rm BB},\mathcal{N}_{\rm BR},\mathcal{N}_{\rm EP},\mathcal{N}_{\rm IP}$. The bonds are classified by using the Baxter-Kelland-Wu loops [25] on the medial lattice, as exemplified in Fig. 5. The loops are constructed by assigning each occupied bond two loop arcs to its two sides, and assigning each unoccupied bond two loop arcs which cross the bond. The types of bonds are determined as following: an occupied bond is a backbone (bridge) bond iff the loop arcs on its two sides do not (do) belong to the same loop, and an unoccupied bond is an external (internal) -perimeter bond if the two arcs crossing the bond do not (do) belong to the same loop. To avoid the case that an occupied bond is a backbone bond but has its two loop arcs in the same loop, the boundary conditions are chosen as periodic in one direction and open in the other direction such that the finite-$L$ square lattice is planar. During the update of a bond state, the connectivity checking and labelling of loop arcs in the loop configuration are performed by simultaneous breadth-first searches starting from two loop arcs of the bond [23, 21]. When a loop is split into two or two loops are merged into one during the update, the simultaneous searching time is only proportional to the size of the smaller loop.

The acceptance probabilities in Alg. 3 are derived as follows. Using $\alpha$ and $\beta$ to represent two configurations that can transform into each other in an update, the detailed balance condition has the form

where $\pi$ is the configuration weight, $P_{pro}$ and $P_{acc}$ are the proposal and acceptance probabilities, respectively. The detailed balance condition associated with the deletion or addition of a backbone bond can be written as

and that associated with the deletion or addition of a bridge bond as

Following the Metropolis-Hastings criterion, the acceptance probabilities are calculated as

Our results of $\tau_{\mathcal{N}_{\rm HF}}$ at $v_{c}=\sqrt{Q}$ for the improved Sweeny algorithm is shown in Fig. 6. When fixing the size $L$, it can be seen that, as $Q$ decreases, $\tau_{\mathcal{N}_{\rm HF}}/L^{2}$ is almost independent of $Q$ before a turning point $Q_{c}(L)$, after which $\tau_{\mathcal{N}_{\rm HF}}/L^{2}$ tends to increase as $\sim Q^{-1/2}$. The value of $Q_{c}(L)$ decreases with increasing system size $L$, and it is expected that $Q_{c}(L)\rightarrow 0$ in the limit $L\rightarrow\infty$. Thus, comparing with results of the original Sweeny algorithm (i.e., Fig. 3), the improved Sweeny algorithm greatly improves simulation efficiency. For $Q=10^{-6}$ and $L=32$, the autocorrelation time for the improved algorithm is of order $O(10)$, and it decreases with $L$, making the algorithm sufficient for most studies of the random-cluster model with small $Q$.

The performance of the improved Sweeny algorithm can be understood from the acceptance probabilities: $P^{(+)}_{BB}$ and $P^{(-)}_{BR}$ are of order $O(1)$ for $Q\gg Q_{c}(L)$; $P^{(-)}_{BB}$ and $P^{(+)}_{BR}$ are of order $O(1)$ for all $Q$. The residual $Q^{-1/2}$ behavior is attributed to that $P^{(+)}_{BB}$ and $P^{(-)}_{BR}$ are of order $O(\sqrt{Q})$ for $Q\ll Q_{c}(L)$. Detailed analyses of these acceptance probabilities will be presented in Section IV. In the limit $Q\rightarrow 0$, the improved Sweeny algorithm fails since $P^{(-)}_{BR}=P^{(+)}_{BB}=0$, i.e., no bridge bond can be deleted from a uniform spanning tree and no backbone bond can be added to the tree.

In the limit $Q\rightarrow 0$, since both the original and improved Sweeny algorithms are invalid, one has to employ other methods. When $Q\rightarrow 0$, though no bond can be erased or occupied, bonds may exchange their states. Thus, to simulate models with $Q\rightarrow 0$, we can use the Kawasaki algorithm [20], in which a randomly selected empty bond is proposed to switch with a randomly selected occupied bond.

We applied the Kawasaki algorithm to the uniform spanning trees, which appear in the limit of $Q\rightarrow 0$ with $v_{c}/Q^{1/2}$ being finite, for the critical square-lattice model. The results are shown in Fig. 7(a), where it can be seen that $\tau_{\mathcal{N}_{\rm HF}}/L^{2}$ diverges as $\sim L^{3/2}$. We understand this $L$-dependent slowing-down as follows.

A Kawasaki update can be performed in two steps: first deleting a randomly chosen occupied bond so that the spanning tree is split into two trees, then occupying a randomly chosen empty bond. We call the configuration between the two steps an intermediate state. For a uniform spanning tree, the overall update would be accepted if and only if the chosen empty bond is an external-perimeter bond in the intermediate state.

The number of external-perimeter bonds in the intermediate states $\mathcal{N}^{*}_{\rm EP}$ can be counted by observing changes of the loops: on a spanning tree, the deletion of an occupied bond decompose a large loop into two loops, and $\mathcal{N}^{*}_{\rm EP}$ is proportional to the size of the external perimeter associated with the smaller loop. The loop surrounding a cluster defines its hull, and the fraction of loop arcs in the hull corresponding to external-perimeter bonds is called the external perimeter of the cluster. The size of the smaller loop is determined as $\sim L^{d_{\rm H}-x_{2}}$ in Ref. [23], where $d_{\rm H}$ is the fractal dimension of the hull and $x_{2}$ is the two-arm exponent [26]. Analogously, we conjecture that the size of the external perimeter associated with the smaller loop is $\sim L^{d_{\rm EP}-x_{2}}$, where $d_{\rm EP}$ is the fractal dimension of the external perimeter. According to the literature [26, 27, 28, 23], the fractal dimensions take values $d_{\rm H}=1+2/g$ and $d_{\rm EP}=1+g/8$; the parameter $g$ and $k$-arm exponents $x_{k}$ are given by

For $Q\rightarrow 0$, one gets $g=2$, $x_{2}=3/4$, and $d_{\rm EP}=5/4$. Then we expect $\langle\mathcal{N}^{*}_{\rm EP}\rangle\sim L^{5/4-3/4}\sim L^{1/2}$, which is numerically confirmed in Fig. 7(b). Thus, the probability that a randomly selected empty bond is an external-perimeter bond is $\langle\mathcal{N}^{*}_{\rm EP}\rangle/2L^{2}\sim L^{-3/2}$, which leads to $\tau_{\mathcal{N}_{\rm HF}}/L^{2}\sim L^{3/2}$ on the uniform spanning trees.

To eliminate the above $L$-dependent slowing-down, similar to the idea in improving the Sweeny algorithm, we make use of the real-time information of four kinds of bonds and propose an improved Kawasaki algorithm, as shown in Alg. 4. This improved algorithm either exchange a bridge bond with an external-perimeter bond, or exchange a backbone bond with an internal-perimeter bond. Since the proposed update always exist and the acceptance probability is $1$, the improved Kawasaki algorithm does not have the above $L$-dependent slowing-down. This is confirmed by our simulation results which have $\tau_{\mathcal{N}_{\rm HF}}/L^{2}\simeq 0.55$ for $L=4$ to $128$.

For nonzero small $Q$, since updates in the (improved) Kawasaki algorithm cannot change the total number of occupied bonds, one has to employ the (improved) Sweeny algorithm. In Section III, for the improved Sweeny algorithm at fixed $L$, we have shown the appearance of a $Q$-dependent slowing-down after an turning point $Q_{c}(L)$. This can be understood by observing the acceptance probabilities as follows. First, for updates of deleting a backbone bond or adding a bridge bond, substituting Eqs. (8) and (9) into Eqs. (14) and (17), one has

where $\mathcal{N}_{\rm IP}(\alpha)+1=\mathcal{N}_{\rm IP}(\beta)$ and $\mathcal{N}_{\rm BR}(\beta)+1=\mathcal{N}_{\rm BR}(\alpha)$ are used. Then, for updates of deleting a bridge bond or adding a backbone bond, more considerations are needed. For example, after deleting a bridge from configuration $\alpha$, there is no simple relation between $\mathcal{N}_{\rm EP}(\beta)$ and $\mathcal{N}_{\rm EP}(\alpha)$. Instead, one has $\mathcal{N}_{\rm EP}(\beta)=\mathcal{N}^{*}_{\rm EP}$, which scales as $\langle\mathcal{N}^{*}_{\rm EP}\rangle\sim L^{d_{\rm EP}-x_{2}}$ in the limit $Q\rightarrow 0$ as previously discussed. For nonzero small $Q$ and large $L$, we expect that Eq. (8) still holds, thus suppose that $\langle\mathcal{N}^{*}_{\rm EP}\rangle\sim L^{2}\sqrt{Q}+aL^{d_{\rm EP}-x_{2}}$. Substituting this result and Eq. (9) into Eq. (16), we get the acceptance probability

Thus, we expect that the autocorrelation time scales as the inverse of the above acceptance probability

From the above equation, for a fixed $Q$, as $L$ increases, since $d_{\rm EP}-x_{2}-2<0$, one has $\tau_{\mathcal{N}_{\rm HF}}/L^{2}\to{\cal O}(1)$; and for a fixed $L$, as $Q$ decreases, $\tau_{\mathcal{N}_{\rm HF}}/L^{2}\to Q^{-1/2}$. These imply the turning near $Q_{c}(L)$ and explain the results of the improved Sweeny algorithm, i.e., Fig. 6.

Finally, since the acceptance probabilities for operations in the improved Kawasaki algorithm are also $1$ for nonzero $Q$, we can combine the improved Sweeny algorithm with the improved Kawasaki algorithm. Our simulation results for this combined method are shown in Fig. 8, where $\tau_{\mathcal{N}_{\rm HF}}/L^{2}$ versus $L$ is almost a straight line. Thus the improved Kawasaki method eliminates the $\sim Q^{-1/2}$ form slowing-down of the improved Sweeny algorithm at fixed $L$, and the combined improved Sweeny-Kawasaki method serves as a complete solution for simulating the random-cluster model with small $Q$.

To summarize, by introducing the number of half faces, we observe a $Q$-dependent slowing-down in two dimensions, which further enriches the dynamics of the Sweeny algorithm for the random-cluster model. The slowing-down is understood by observing the updating probabilities, which contains the critical bond weight $v_{c}=\sqrt{Q}$ on the square lattice. By classifying the bonds into four types and dynamically updating their type information, we propose an improved Sweeny algorithm, which (at fixed $L$) pushes the slowing-down to much smaller $Q$-region. Since the $Q$-dependent slowing-down renders the (improved) Sweeny algorithm being non-ergodic for $Q\rightarrow 0$, we resort to the Kawasaki algorithm. The latter algorithm is found to exhibit an $L$-dependent slowing-down, and this slowing-down can be fully suppressed in an improved version of the algorithm by making use of the bond-type classification. Finally, we formulate an improved Sweeny-Kawasaki method which combines advantages of the two improved algorithms and eliminates the above critical slowing-down.

We note that updates in the improved algorithms proceed with the aid of Baxter-Kelland-Wu loops, as it is convenient to classify the bonds using the loops, and it is also faster to perform connectivity checking in the loops than in the clusters. The computational cost of the algorithms primarily stems from the connectivity checking: simultaneous breadth-first search [23, 21] in the clusters has a runtime scaling of $\sim L^{d_{\rm F}-x_{2}}$, while that in the loops has a runtime scaling of $\sim L^{d_{\rm H}-x_{2}}$ with $d_{\rm H}\leq d_{\rm F}$. For not too large $L$, the improved Sweeny-Kawasaki method in this work should be adequate for simulating the random-cluster model. To speedup the simulations for large $L$, it should be worth to study how to incorporate dynamic connectivity algorithms which have runtime scaling at most logarithmic in the system size [21, 22].

It is also noted that, solely for the uniform spanning trees, there exist algorithms which have excellent performance. These include the celebrated Wilson’s algorithm [29], which uniformly constructs a spanning tree from scratch by a loop-erased random walk, and a smart Kawasaki algorithm which uniformly generates spanning trees by using a fast link-cut tree data structure [30]. The latter is different from our improved Kawasaki algorithm, which constructs loops to classify the bonds and applies also to models with $Q\neq 0$.

Furthermore, though our simulations were performed at $v_{c}=\sqrt{Q}$ on the square lattice, for $v=Q^{\sigma}$ with general values of $0<\sigma<1$, as $Q\rightarrow 0$ there is also $Q$-dependent slowing-down for the Sweeny algorithm. This is because in the latter case the updating probabilities tend to zero as $Q^{\sigma}$ or $Q^{1-\sigma}$. The improved Sweeny and Kawasaki algorithms can also solve this $Q$-dependent slowing-down, due to reasons similar to those for $v=\sqrt{Q}$. Also, the $Q^{-1/2}$-form critical slowing-down for the random-cluster model is expected to exist on other 2D lattices. For example, on the triangular ($\mathbb{T}$) and honeycomb ($\mathbb{H}$) lattices, when $Q\leq 4$, it is known [31, 1, 4] that

These critical bond weights $v^{\mathbb{T}}_{c}$ and $v^{\mathbb{H}}_{c}$ also vanish as $\sim\sqrt{Q}$ when $Q\rightarrow 0$. However, when the spatial dimension $d\geq 3$, for small $Q$ the critical point behaves as $v_{c}\sim Q$ (hence the limit $Q\to 0$ corresponds to the spanning-forest model) [32]. Thus, the $Q$-dependent critical slowing-down does not occur for $d\geq 3$. Nevertheless, the improved Sweeny and Kawashaki algorithms may still help to improve the simulation efficiency.

Finally, the number of half faces proves to be very useful in demonstrating the slowing-down of the unimproved algorithms. Further work may explore more quantities which have similar behavior. Despite the mechanism discussed in this work, we cannot rule out any other property which may cause more severe slowing-down.