hamiltonian can be written in the form $$H =
  -\sum_{( \vec r, \vec r')}
  \left( \sum_{p=1}^{q-1} M^{[\vec r]}_{x,p} M^{[\vec r']}_{x,q-p} \right)
    -\lambda \sum_{\vec r} M^{[\vec r]}_z,$$ where $M^{[\vec r]}_{x, p}$ and $M^{[\vec r]}_{z, p}$, with $p = 1,\ldots,q-1$, are spin matrices of size $q \times q$ acting at site $\vec r$. The parameter $\lambda$ is the analogue of the transverse magnetic field in the Ising case. In this formulation the spin matrices acting at each site are given by [@qPotts] $$M_{x,1} = \left(\begin{array}{cc} 0 & I_{q-1} \\ 1 & 0 \end{array} \right), \quad
  M_z = \left( \begin{array}{cc} q-1 & 0 \\ 0 & -I_{q-1} \end{array}
  \right),$$ where $I_{q-1}$ is the $(q-1)\times (q-1)$ identity matrix with $M_{x,p} = \left(M_{x,1}\right)^p$ and $\left(M_{x,1}\right)^q = I_q$.

![(color online) Numerical values for the local order parameter $\langle M_{x_1}\rangle$ and the GE per site ${\mathcal E}(\lambda)$ for the 2D quantum $q$-state Potts model for (a) $q=3$, (b) $q=4$ and (c) $q=5$ as a function of the transverse field strength $\lambda$ for the indicated values of the iPEPS truncation dimension $D$. For each case the curves for both the local order parameter and the GE show a discontinuous behavior, indicative of a discontinuous phase transition. []{data-label="fig2"}](q3.eps "fig:"){width="1.0\linewidth"} 