main obstacle to sampling is that the appropriate normalisation factor, the sum of the weights of all colourings, is hard to compute. To be precise: for a graph $G=(V,E)$, a *(spin) configuration* $\sigma$ is a function which assigns to each vertex $i$ a colour $\sigma_i \in \{1, \ldots, q \}$ (also called states or spins). The probability of finding the system in a given configuration $\sigma$ is given by the Gibbs distribution: $$\pi(\sigma) = Z^{-1} e^{\beta \sum_{(i,j)\in E} 
   J \delta(\sigma_i, \sigma_j)},$$ where $\delta(\sigma_i, \sigma_j)$ is the Kronecker-$\delta$ (taking value 1 $\sigma_i = \sigma_j$, and taking value 0 otherwise); $\beta = (kT)^{-1}>0$ is the inverse temperature (here $k$ is Boltzman’s constant and $T$ is temperature); and $Z=Z(G,\beta,J,q)$, is the *partition function* i.e. the appropriate normalisation factor to make this a probability distribution. The strength of the interaction between neighbouring vertices is given by the coupling constant $J$. If $J>0$ then the bias is towards having many edges with like colours at the endpoints; this is the ferromagnetic region. If $J<0$ then the bias is towards few edges with like colours at the endpoints: this is the anti-ferromagnetic region.

Our results concern only the ferromagnetic region, where $J > 0$, although we 