(The simulation data in Fig. 3 are obtained by extrapolating moment ratio crossings [@moments] for system sizes $L=200$, 400, 800 and 1600). The extrapolated value of $m$ at the crossings is not particularly good (for $\mu=1$ we find $m_c = 1.110$, compared with the best estimate of 1.1736(1) [@moments]). Although we expect that this would improve using larger systems, our objective here is to find a relatively fast and simple method to predict the phase boundary. (The cpu time required to converge to the QS distribution is comparable to that required to integrate the equations numerically in the $n=6$ cluster approximation.)

To apply the QS method to the two-dimensional 2SCP, we devised an algorithm that enumerates configurations and transitions for a general graph of $N$ vertices; the graph structure is specified by the set of bonds ${\cal B} = \{(i_1,j_1), (i_2,j_2), ..., (i_m,j_m)\}$ linking pairs of vertices $i_k$ and $j_k$. To represent a portion of the square lattice, with periodic boundaries, each vertex must be linked to four others. This can be achieved rather naturally for a square ($m \times m$) or rectangle ($m \times (m+1)$); for other values of $N$ we use a cluster close to a square, and 