Damage spreading in the random cluster model

We investigate the damage spreading effect in the Fortuin-Kasteleyn random cluster model for 2- and 3-dimensional grids with periodic boundary. For 2D the damage function has a global maximum at $p=\sqrt{q}/(1+\sqrt{q})$ for all $q>0$ and also local maxima at $p=1/2$ and $p=q/(1+q)$ for $q\lesssim 0.75$. For 3D we observe a local maximum at $p=q/(1+q)$ for $q\lesssim 0.46$ and a global maximum at $p=1/2$ for $q\lesssim 4.5$. The chaotic phase of the model's $(p,q)$-parameter space is where the coupling time is of exponential order and we locate points on its boundary. For 3-dimensional grids the lower bound of this phase may be equal to the corresponding critical point of the $q$-state Potts model for $q\ge 3$.


I. INTRODUCTION
The study of damage spreading has its early origin in statistical systems where binary vectors are iteratively mapped by randomly chosen Boolean functions. If two independently chosen vectors then are exposed to the same sequence of functions, how does their Manhattan distance behave with time and when do the vectors become equal? (see for example Ref. [1] and references therein). In the 1980s similar questions were asked for Ising spin systems [2]. Here a replica of an equilibrated 2dimensional spin system is made with one spin changed, giving it an initial damage. Updating the two systems with the Glauber rule while using the same sequence of random numbers to guide the updating, the difference between the two systems, the damage D t , then evolves with time t. For low temperatures T < T c the two systems quickly become identical (D t = 0 at t = τ ) but for high temperatures T > T c the damage spreads through the replica, reaching a positive fraction of the spins (D t > 0).
The damage spreading phenomenon has, since then, mainly been studied for Ising spin glass systems but also in many other contexts [3][4][5][6][7][8][9][10]. Let us describe the phenomenon in the spin-glass case in slightly more detail. Here two independent and randomly chosen ±1-spin vectors of length N , having the same interaction strengths J ij , are exposed to the same sequence of random numbers for updating the spins, this time with heat-bath dynamics. Now the situation becomes the opposite of that just described for the Ising case [2,3]: above a certain temperature T D , in the ordered ("frozen" or "stable") phase the two spin vectors become identical after just a few time steps (say, τ = O(N )) and then stays identical, i.e., the two systems have coupled. Below this temperature, in the chaotic (or disordered) phase, the coupling time τ instead diverges exponentially with N and the two systems never seem to meet. In this phase, the distance D t fluctuates around some temperature-dependent mean value * per.hakan.lundow@math.umu.se D = D t . In short, at T D the coupling time τ changes (with decreasing T ) from staying bounded, through logarithmic and polynomial growth, to exponential growth, see Fig. 3 of Ref. [8].
The opposite behaviour, compared to the Ising case, is due to the different dynamics: heat-bath versus Glauber updating [2]. Recall that with heat-bath dynamics a spin is set to its new value with some probability, whereas for Glauber dynamics a spin is instead flipped with some probability. Note also that for spin glasses, T D is not the same as the critical temperature T c , as in the 2D Ising case. For example, for 2D systems where J ij = ±1 with equal probability, we have T D ≈ 1.69 [7] even though T c = 0.
Here we instead investigate damage spreading in the setting of the random-cluster model where two parameters p and q control the distribution of states [11,12]. With this model it is the state of edges (on/off) rather than the state of spins (±1) that is updated. The underlying graphs will be 2-and 3-dimensional grids of linear order L. Recall that the random cluster model contains the percolation model (q = 1) and the q-state Potts model (integer q ≥ 2) which in turn generalises the Ising model without an external field.
For smaller L we obtain an support set S L of the damage function D(p, q) by simply updating a system pair for a long time, depending on L, to see if D → 0 in that time. The chaotic phase S is obtained from estimating the growth rate of the coupling time for many choices of parameters p, q. Effectively the chaotic phase is a limit of the support sets, S L → S. However, determining S L for larger L can be computationally demanding near its boundary, but the small L tells us where to look for the boundary of S. Thus it can actually be easier to estimate the boundary of S by simply measuring growth rates of the coupling time.
We have computed the damage function and estimated the support set for 2D and 3D systems by scanning a grid of p, q-parameters and for a range of L. Based on this we have approximated the boundary of the chaotic phase. When q < 1 for 2D grids, this boundary has a smooth, almost triangular shape, but is considerably more complex for 3D grids. We also guess simple formulas describing the boundary. For 3D grids with q ≥ 3 the lower bound of the chaotic phase seems to be closely related to the critical point of the q-state Potts model.

II. DEFINITIONS AND DETAILS
The underlying graph G = (V, E) is the d-dimensional grid graph of order L with periodic boundary conditions thus having N = L d vertices and m = dL d edges. The partition function of the Fortuin-Kasteleyn (FK) model [11], for parameters 0 < p < 1 and q > 0 is where = (A) = |A| is the number of edges in A and κ = κ(A) is the number of (connected) components, or clusters, in the spanning subgraph G(V, A). For convenience, we will overload the notation and write (p, q) = /m and κ(p, q) = κ /N for the normalised mean values under the distribution of states. The parameter p controls the edge-density, in fact (p, q) is nondecreasing in p for fixed q > 0 [12]. Parameter q controls the bias towards clusters; q > 1 means we prefer many clusters, q < 1 means we prefer less of them, and q = 1 gives a random subset of the edges E without any such bias.
We will describe the process of updating A in some detail. The state of an edge e ∈ E is either e ∈ A or e / ∈ A. We wish to update this state by choosing at random whether it should be reversed (flipped). Reversing the state of e = {u, v} changes by ∆ = ±1. However, ∆κ depends both on the state of e and whether there is a uv-path in A not using e, here denoted u v. The four possible cases are then An edge e ∈ E selected at random will then belong to one of these cases and we denote the probability for these cases P 1 , P 2 , P 3 , P 4 . Alternatively, the underlying graph being edge transitive, as the present grid graphs, we might also say that for a fixed edge, the probability that this edge belongs to case i is P i . Note, by the way, that P 1 +P 2 is the probability that e ∈ A so that (p, q) = P 1 + P 2 .
As usual, briefly put, our update rule should maintain so-called detailed balance between the current set A and the next set A , differing only at the edge e, The probability ratio between the states is With F (x) = Pr(A → A ) we have F (1/x) = Pr(A → A) and thus we want a function satisfying 0 ≤ F (x) ≤ 1 and F (x) = xF (1/x). Choosing a random number U uniformly distributed over the interval [0, 1), the update rule is now simply: reverse the state of edge e if U ≤ F (x). We have mainly used the traditional F (x) = min(1, x) of Metropolis fame, but we will also briefly use F (x) = x 1+x (heat-bath) for comparison. Functions like F (x) = x+x 2 +...+x k 1+x+x 2 +...+x k work equally well but do not seem to be used in the literature. See Ref. [12] for a more complete description of Metropolis/heat-bath in the random cluster context and Ref. [13] for more on update rules.
The process starts at time step t = 0 with a random subset A ⊆ E where Pr(e ∈ A) = Pr(e / ∈ A) = 1/2. A second set B ⊆ E is initially chosen as the complement of A, i.e., B = E \ A. The damage can now be defined where I(s) is an indicator function returning 1 if s is true and 0 otherwise. At t = 0 the damage is thus always D = 1. Both sets are now exposed to exactly the same sequence of random numbers and the edges in both sets are updated according to the protocol described above. A time step consists of performing a sweep where we update m randomly selected edges from E. The same sequence of edges are then candidates for addition/deletion at both sets and they are tried with the same random numbers. If we at some time step t receive D t = 0, the coupling time τ has been reached and the damage will stay zero at all future time steps. The damage will eventually go to zero but the coupling time divides the p, q-parameter space into a chaotic phase and an ordered phase. In the ordered phase, the coupling time is of at most polynomial order, τ = O(L c ) for some finite c, thus potentially D t → 0 comparatively fast. In the chaotic phase the coupling time is instead super-polynomial, typically exponential. Except for small L we will in practice never obtain D t = 0 in this phase. Instead, the damage stabilizes after some equilibration, around a value D(p, q) = D t .
When moving along a line in the p, q-parameter space, from the chaotic to the ordered phase, we typically find τ = O(L c ) at a single point. We will take a simple view on the transition between the two phases, going from exponential, via polynomial at a point, to subpolynomial, say logarithmic or bounded τ . However, we do not wish to rule out that this transition may be more complex than this.
For a given L and after some practical number of equilibration steps, we can treat D(p, q) as, effectively, a function with support set S L = {(p, q) : D(p, q) > 0}. The chaotic phase will then be the set of (p, q)-points where τ is of superpolynomial order, i.e., τ = ω(L c ) for all constants c, using little-omega notation. At least for small L, where τ is still not too large, it is practical to view S as a limit set of S L . We will denote by S(q) the set for fixed q. Denote by S(q) andS(q) the lower and upper bounds for fixed q.
For d = 2 we define which for q ≥ 1 is known to be the critical point of the random cluster model [12,14], though we will use this notation also for q < 1. For d = 3 no general formula for p c (q) is known but estimates exists of course for specific integer q-values in the context of Ising-or Potts-model estimates. Recall that the translation between the temperaure parameter β of the q-state Potts model and p of the random cluster model is through the simple relation p = 1 − e −β , see for example Refs. [12,15]. For d = 3 this gives p c (2) = 0.358091335(6) (via Ising model β c -estimate [16]), p c (3) = 0.42330(5) (Potts model estimates, see for example Refs. [17,18]), p c (4) = 0.4666(3) (see for example Refs. [19,20]). For q = 5 the phase transition becomes first order and β c is harder to estimate. We are unfortunately not aware of any published estimates of β c for q ≥ 5.
In the present work we have collected random cluster data at a large number of p, q-points for L = 6, 7, . . . , 32 when d = 2. and for L = 3, 4, . . . , 20 when d = 3. For both dimensions we also have data for some larger L but at much fewer p, q-points. We also have gathered coupling data at many points around the border of S, sometimes for larger L. All data were collected on a MacBook Pro laptop computer and a PowerMac desktop computer, both with an 8-thread processor, running programs written in Fortran. All code can of course be obtained from the author.

III. TRAJECTORIES AND COUPLING TIME
In Fig. 1 we show the median damage trajectories for a 32 × 32-grid for three different p-values with q = 1/2. With Metropolis updating for p = 0.27 the damage reaches 0 at τ = 94(2), for p = 0.28 at τ = 415(10) while for p = 0.29 the damage does not reach zero within this time frame. Using instead heat-bath updating the coupling times change to τ = 150(2) and τ = 620 (20) respectively. Thus we expect larger coupling times with heat-bath updating. As the inset of Fig. 1 shows, for p = 0.29 using the Metropolis rule the damage stabilises very quickly to D ≈ 0.0924 after roughly 200 sweeps while for the heat-bath rule the equilibration takes considerably longer, ca 600 sweeps and then ends up at a lower value, D ≈ 0.0868. These examples suggest that the Metropolis rule leads to a faster convergence both in the ordered and the chaotic phase. It should here be mentioned that at p c (1/2) ≈ 0.4142 we see no difference between the respective D, when t 50, though the heatbath curve is still slower to equilibrate. Let us now instead look at the coupling time τ plotted versus L. In Fig. 2 we show the median τ taken over many trajectories (error bars from bootstrap method) for a range of L and for the same p, q-values as in Fig. 1, for both heat-bath and Metropolis updating. Clearly the coupling time differ roughly by a factor for these p-values and we will henceforth only focus on the Metropolis updating.
In principle one should be able to locate a p-value on the boundary of S(q) as a straight line in a log-log plot of τ versus L. This is often true, but in this case we see rather strong finite-size effects of τ , which complicates matters. At this point we can only say for certain that p = 0.29 is in the chaotic phase since the coupling time clearly grows exponentially and is well-fitted by the estimate τ = 0.59(3) exp(1.06(1)L) for the Metropolis rule.
At p = 0.27 the coupling time grows logarithmically and is well-fitted over L ≥ 24 by τ = −68(2)+46(1) log L. For p = 0.28 the small-L effects are too strong to obtain a useful fit. Thus, in the presence of such strong finite-size effects we receive a rough estimate, S(1/2) = 0.28(1).
Let us instead look at the damage function D(p, 1/2) near p = 0.28 as shown in Fig. 3. Damage curves for a range of L agree on a common limit curve and the curve stay positive closer to the boundary for larger L. From the data points in the plot we estimate the limit curve D = 3.80(p − 0.2806) 0.80 , giving an excellent fit. This approach gives the estimate S(1/2) = 0.2806 (2). One can thus resort to this approach when τ has strong finite-size effects if high accuracy is important. Usually, however, we will settle for the accuracy in the first approach.
We will also show another example, see  In the right panel of Fig. 4, the damage D(p, 1/2) is shown. Again, from standard log-log fits, we find the simple formula D = 3.78(0.5625 − p) 0.845 which fits these curves very well. Unfortunately, due to slow equilibration, it is much more difficult to obtain good D-data at this end of S(1/2) and it would require larger systems to improve significantly on the accuracy in our estimate ofS(1/2). Hence, locating the boundary of S by searching for polynomial order τ is actually practical in this case. This will in fact be the main approach used to locate the boundary of S. The accuracy of the estimates are indicated by error bars but often the error is smaller than the points in the figures. Having established our approach we can now show the chaotic phase for d = 2, 3 in the coming sections.
IV. THE CHAOTIC PHASE, d = 2, 0 < q < 1 A rough estimate of S can be obtained by equilibrating D for as many sweeps as possible, over a fine grid of p, q-points, for some L. We show these for L = 8 and 32 in Fig. 5. The points indicate estimated boundary points of S, including also the points (p, q) = (0, 0) and (1, 1/2). The upper bound is closely approximated by the simple formulaS(q) = 0.50 + 0.086(2) √ 1 − q, but the lower bound S(q) has a more complicated behaviour for which we have no suggested formula.
The dotted lines in Fig. 5 indicate the location of three local maxima of D(p, q) for fixed q, namely: p = 1/2, p = p c (q) and p = p c (q 2 ) = q/(1 + q), see Eq. (6). This is demonstrated in Fig. 6 for q = 1/2 where the peaks are located exactly, as far as we can tell, at p = 1/2, p = p c (1/2) = √ 2 − 1 and p = p c (1/4) = 1/3. This picture is quite representative for q 0.75, while for q 0.75 we see instead a triangle-shaped function with only one cusp-like maximum at p = p c (q).
The global maximum of D at p = p c (q) supports that Eq. (6) indeed is a critical p-value also for q < 1 as well, seemingly without finite-size corrections. It is far from clear why p = 1/2 and p = q/(1 + q) also would give local maxima. Note that both p = 1/2 and p = q/(1 + q) gives cusp-like singularities in the edge-flip probability P f . This is not surprising since P f is a linear combination of P i and F (x) with x depending on the cases of Eq. (2). With F (x) = min(1, x) (Metropolis) both of these p-values give this effect on P f , see Fig. 7. We find this effect in the damage function remarkable but can only note the apparent relationship with the edge-flip probability.
The damage function differs very little when using heat-bath updating but the edge-flip probability is instead quite smooth without any cusps as in Fig. 6.
For q > 1 the chaotic phase surrounds p c (q) in a narrow band as shown in Fig. 8 for 1 < q ≤ 2. Clearly the finitesize S L is growing outwards, approaching S. Estimating the boundary of S proceeds as for q < 1. We note that the lower bound of S is sometimes difficult to locate due to substantial finite-size effects in τ . However, almost no finite-size effects are seen for the upper bound of S. An example is shown in the left panel of Fig. 9 where we show the behaviour of τ nearS (2).
At the boundary, estimated toS(2) = 0.6325(25), we  The right panel of Fig. 9 shows the damage function D at q = 2 for a range of L. There are stronger finitesize effects away from the maximum. Just as for q < 1, the damage maximum seems to be located at exactly p c (q) = √ 2/(1 + √ 2), for all L and for all 1 < q ≤ 6 that we have checked. For 3-dimensional systems the chaotic phase has a quite different shape from that of d = 2. In Fig. 10 we show the support sets (gray) for L = 6, 10, 12 after several thousand sweeps. For L ≤ 9 these sets are disconnected, splitting into two parts along p ≈ 0.42. The points in the figure (with error bars) show our attempt at pinpointing the boundary of the chaotic phase. It turns out that a good approximation for the lower bound is the simple S(q) = q/(2.805(1) + q), for q 0.75, while the upper bound is excellently fitted byS(q) = 0.50 + 0.0300(7)(1 − q) 0.75 for q < 1. Both these curves are shown in Fig. 10.
A complicating feature of this region is the wedgeshaped lower right boundary which stretches up to roughly (p, q) ≈ (0.21, 0.75). Then it climbs upwards to the left in the figure to about (p, q) ≈ (0.42, 0.40) before finally taking a more spike-shaped form. The error bars are often rather large due to quite significant finitesize effects in τ . Determining the correct shape here will require a more ambitious computational investment. In most other cases around the boundary, the finite-size effects are negligible, thus giving very small error bars.
The damage function D(p, q) has a considerably less regular behaviour for d = 3 than for d = 2. We show an example for q = 0.1 and 0.3 in Fig. 11. For fixed q there is a cusp-like maximum at p = 1/2 for all 0 < q < 1. There is also a cusp-like maximum at p = q/(1 + q), though it is no longer a maximum for q 0.46. More peculiar is that there is a cusp also at p = √ q/(1 + √ q) of Eq. (6), but again this vanishes for q 0.46. We have no explanation for why a critical point for 2D would occur in the 3D-case. For 0.46 q 0.66 a new local maximum appears to the left of p = q/(1 + q), but we do not have a conjecture for this point. For 3-dimensional systems the chaotic phase for q > 1 is interesting as well. In Fig. 12 we show the support sets for L = 6, 8, 12 together with some estimated boundary points of the chaotic phase. Again we find a very good approximation for the upper boundary,S(q) = 0.50 + 0.0333(1) √ q − 1, for q 20 (!). For q 20 the formula does not give a good prediction of the upper bound and we do not know the shape of S beyond this point.
The lower boundary points S(q) agree very well with the critical inverse temperature of the q-state Potts model, after using the relation p c (q) = 1 − e −βc(q) . Let us begin with q = 2, where p c (2) ≈ 0.358091. In Fig. 13 we show τ versus L around this value and the damage function along q = 2. In this case it appears that τ grows exponentially at p = 0.358091, approximately τ  S(2) < p c (2) though the difference is only about 0.5%.
In somewhat less detail, let us also consider q = 3, recall that p c (3) = 0.42330 (5). In Fig. 14  For q = 4, as shown in Fig. 15, the finite-size effects in τ complicate the picture and it is difficult to estimate S(4) with good accuracy. At p = 0.467, just above the estimate p c (4) = 0.4666(3) we see a distinctly exponential growth rate in τ . At p = 0.465, just below the estimate, τ actually appears to be bounded. We show also τ for p = 0.466 but it would take larger L to decide if it has polynomial growth. This gives us the estimate S(4) = 0.466(1). This estimate of S(4) then overlaps the estimate of p c (4).
We continue with q = 5 and give an estimate of S(5). In the left panel of Fig. 16 we see that p = 0.49625 gives distinct exponential growth in τ and at p = 0.49375 the growth is logarithmic. We also show τ for p = 0.495 but  the small-L effects make it difficult to say which category it belongs to. We obtain the estimate S(5) = 0.4950 (13) which translates to β c (5) = 0.6832(26), if we assume that p c and the lower bound are equal. Analogously, we find for q = 6 that the lower bound of the chaotic phase is S(6) = 0.518(1) suggesting β c (6) = 0.730(1) for the 6state Potts model.
In the right panel of Fig. 16, showing D(p, 5), we note that the maximum is no longer located at p = 1/2, but rather at p ≈ 0.51. In fact, our data suggest that the maximum is located at p = 1/2 for q ≤ 4.5 but at some p larger than 1/2 for q > 4.5.
The phenomenon of an emerging discontinuity in D(p, q) is also reflected in other fundamental quantities of the random-cluster model. In Fig. 17 we show (p, q) and  κ(p, q) versus p for q = 2, 3, 4, 5, 6 and L = 16. A discontinuity is clearly emerging with increasing q, becoming more distinct at q = 4, and it is located at exactly the same point where the damage drops to zero. Thus it is perhaps not surprising that the critical temperature of the q-Potts model is connected to this point. On the other hand, the upper bound of S does not seem to have a corresponding effect on and κ. For rigorous results on the discontinuity of , see Section 7.5 of Ref. [12])

VIII. CONCLUSIONS
We have investigated the phenomenon of damage spreading in the random-cluster model for grids of dimen-sions two and three with periodic boundary conditions. The chaotic phase, where the coupling time τ grows exponentially with L, appears as an island in the p, q-plane and here the damage function takes positive values and its average remains stable for many time steps.
For 2-dimensional grids the damage function has a global cusp-shaped maximum at exactly p = p c (q) = √ q/(1 + √ q). When q 0.75 there are also local maxima, again cusp-shaped, at p = q/(1 + q) and p = 1/2. For comparison, we observe that the edge-flip probability, under Metropolis updating, also has local cusp-like maxima at p = q/(1 + q) and p = 1/2. We have estimated lower and upper bounds for the chaotic phase and suggested an approximation formula forS(q) when q < 1.
For 3-dimensional grids the damage function has a richer behaviour, depending on q. Notably, for q 0.46 the damage function has a cusp-shaped local maximum at p = q/(1 + q) and we observe a cusp also at p = √ q/(1 + √ q). The chaotic phase has a more complex shape for q < 1 but we give good approximation formulas for the upper bound when q < 20 and for the lower bound when q < 0.75. When q > 2 is integer, the lower bound appears remarkably close to the critical point p c (q) of the q-state Potts model. They differ slightly for q = 2 but estimates of p c (q) and S(q) overlap for q = 3, 4. For q > 2 the damage function (with fixed q) develops a discontinuity at S(q) and the same effect is observed for both and κ at this point.