Potts model: step in computing the Gibbs free energy using the variational approach To be clear: This is a pure computational question, it is not about the physics of the Potts model.
In the Potts model the Hamiltonian is given by:
$\mathcal{H} = \sum_{i,j}(1-\delta_{\sigma_i\sigma_j})$,
where $\delta_{ij}$ is the Kronecker delta and $\sigma_{i,j}, \tau \in \{1,\dots,q\}$ and $i \in \{1,\dots,N\}$ where the sum $\sum_{i,j}$ is the sum over nearest neighbour pairs. I want to compute the Gibbs free energy by using the variational approach. Eventually one has to solve some averages, and this is where my question comes in. Given a probability distribution $\mathcal{Q} = \Pi_{i=1}^{N} \left( \sum_{\tau}\rho_{\tau}\delta_{\sigma_i \tau}\right)$ where $\rho_{\tau}$ is the fraction of spins in state $\tau$, how do you compute:
\begin{equation}
\langle \sum_{i,j}\delta_{\sigma_i\sigma_j}\rangle_{\mathcal{Q}}
\end{equation}
?
This comes from the need of calculating $-\beta\langle \mathcal{H}\rangle_{\mathcal{Q}}$.
 A: First, we can simply rewrite the sum more explicitly as pairs of spins interacting :
\begin{equation}
   \left\langle \sum_{\langle i, j \rangle} \delta_{\sigma_i, \sigma_j}
   \right\rangle = \sum_{\langle i, j \rangle} \langle \delta_{\sigma_i,
   \sigma_j} \rangle = \sum_{\langle i, j \rangle} \sum_{\tau} \langle
   \delta_{\sigma_i, \tau} \delta_{\sigma_j, \tau} \rangle
\end{equation}
Then, you have to realize that the distribution $\mathcal{Q}$ is a product of independent spin distributions, and thus actually describes decoupled spins, so that
\begin{equation}
   \langle \delta_{\sigma_i, \tau} \cdot \delta_{\sigma_j, \tau}
   \rangle_{\mathcal{Q}} = \langle \delta_{\sigma_i, \tau}
   \rangle_{\mathcal{Q}} \cdot \langle \delta_{\sigma_j, \tau}
   \rangle_{\mathcal{Q}} = \rho_{\tau}^2
\end{equation}
Thus,
\begin{equation}
   \left\langle \sum_{\langle i, j \rangle} \delta_{\sigma_i, \sigma_j}
   \right\rangle_{\mathcal{Q}} = \sum_{\langle i, j \rangle} \sum_{\tau}
   \rho_{\tau}^2 = \frac{Nz}{2} \sum_{\tau} \rho_{\tau}^2
\end{equation}
where $z$ is the coordinance of the lattice.
This is actually a mean-field description of the Potts model, because spins are assumed independent. The free energy you will obtain in the Bragg-Wiliams free energy. Even if this is very approximate, especially at low dimensions and low number of spin directions, it stills yields important insights on the behavior of the Potts model and the phase transition which occurs in this model.
