# References

The discrete $`_6`$ spin models describe the behaviour of a collection of spins $`S_i=\mathrm{exp}(i\theta _i)`$, where the $`\theta _i`$ are integer multiples of $`\pi /3`$. These spins live on a square two-dimensional lattice and interact according to a reduced Hamiltonian of the form

$$H\frac{1}{k_BT}=\underset{<ij>}{}V(\theta _i\theta _j),$$
(1)   
the sum running over nearest-neighbour pairs of sites $`<ij>`$. Imposing $`V(\theta )=V(\theta )`$, a particular system is characterised by the three numbers $`V_r=V(\pi r/3)V(0)`$, $`r=1,2,3`$. A duality transformation maps the space of such models to itself: with the couplings parametrised by variables $`x_r=\mathrm{exp}(V_r)`$, a sum over the $`\theta _i`$ is equivalent to one over dual variables $`\stackrel{~}{\theta }_{\stackrel{~}{ı}}`$, with the triplet of couplings $`\{x_r\}`$ replaced by the dual set $`\{\stackrel{~}{x}_r\}`$:

$`\stackrel{~}{x}_1`$ $`=`$ $`(1+x_1x_2x_3)/\mathrm{\Delta }`$   
$`\stackrel{~}{x}_2`$ $`=`$ $`(1x_1x_2+x_3)/\mathrm{\Delta }`$   
$`\stackrel{~}{x}_3`$ $`=`$ $`(12x_1+2x_2x_3)/\mathrm{\Delta }`$ (2)   
where $`\mathrm{\Delta }=1+2x_1+2x_2+x_3`$. The transformation leaves invariant points on the line

$$(x_1,x_2,x_3)=(t,\beta \alpha t,\alpha 2\beta t)$$
(3)   
where $`\alpha =3\sqrt{6}`$ and $`\beta =\sqrt{6}2`$. This self-dual line intersects the line of 6-state Potts models at the point $`x_1=x_2=x_3=(\sqrt{6}1)/5`$, denoted $`P`$ below. There is also a distinguished surface on which the models decouple into independent three-state Potts and Ising models. This is revealed on writing each spin $`S_i`$ as a product of a $`_2`$ and a $`_3`$ valued variable:

$$S_i=e^{i\theta _i}=\mathrm{\Sigma }_i\sigma _i;\mathrm{\Sigma }_i=1,e^{\pm 2\pi i/3},\sigma =\pm 1$$
(4)   
The pairwise interaction energy can then be written as

$`V(\theta )`$ $`=`$ $`J_1[1\mathrm{cos}(\theta )]+J_2[1\mathrm{cos}(2\theta )]+J_3[1\mathrm{cos}(3\theta )]`$ (5)   
$`=`$ $`J_1[1\frac{1}{2}(S+S^1)]+J_2[1\frac{1}{2}(\mathrm{\Sigma }+\mathrm{\Sigma }^1)]+J_3[1\frac{1}{2}(\sigma +\sigma ^1)]`$   
$`=`$ $`3J_1[1\delta _{\sigma _i\sigma _j}\delta _{\mathrm{\Sigma }_i\mathrm{\Sigma }_j}]+\frac{3}{2}(J_2J_1)[1\delta _{\mathrm{\Sigma }_i\mathrm{\Sigma }_j}]+(2J_3J_1)[1\delta _{\sigma _i\sigma _j}]`$   
where $`\theta =\theta _i\theta _j`$, $`S=S_i/S_j`$, $`\mathrm{\Sigma }=\mathrm{\Sigma }_i/\mathrm{\Sigma }_j`$, $`\sigma =\sigma _i/\sigma _j`$, and

$$J_1=\frac{1}{3}\mathrm{ln}\frac{x_1}{x_2x_3}J_2=\frac{1}{3}\mathrm{ln}\frac{x_3}{x_1x_2}J_3=\frac{1}{6}\mathrm{ln}\frac{x_2^2}{x_1^2x_3}.$$
(6)   
The spins $`\mathrm{\Sigma }_i`$ and $`\sigma _i`$ decouple on the surface $`J_1=0`$, $`x_1=x_2x_3`$, on which $`\frac{3}{2}J_2`$ is the 3-state Potts coupling, and $`2J_3`$ the Ising coupling.

The more general phase structure of these models has been studied by various authors over the years (see for example ). The question is an interesting exercise in its own right, and is also of wider relevance – to, for example, the effect of hexagonal symmetry-breaking on the isotropic planar model , the behaviour of the cubic model , and the spectra of Heisenberg antiferromagnetic spin chains . An initial phase diagram was suggested by Domany and Riedel in , with the three-dimensional thermodynamic phase space partitioned into four domains, one disordered and the others exhibiting $`_2`$, $`_3`$ and $`_6`$ ordering. The existence of an additional massless phase was then demonstrated, first along the Villain line, and then throughout a whole three-dimensional region. This phase was incorporated into the Domany-Riedel diagram by Alcaraz and Koberle in, its end point on the self-dual line later being identified with a particular integrable point $`C`$ on the phase diagram (further support for such an identification can be found in).

However there are reasons to believe that the story is not yet complete. We give two examples.

(i) In the diagrams of and, the point $`P`$ touches surfaces of transition from the disordered phase into regions of $`_2`$ and $`_3`$ order. Such transitions are expected to be of second order, and so the correlation length on these surfaces should be infinite. This contradicts the known behaviour of the 6-state Potts model at $`P`$, where the correlation length remains finite, albeit large.

(ii) On the decoupling surface $`J_1=0`$, there is a line of Ising transitions, $`J_3=\frac{1}{2}\mathrm{ln}[1+\sqrt{2}]`$, and a line of 3-state Potts transitions, $`J_2=\frac{2}{3}\mathrm{ln}[1+\sqrt{3}]`$. These lines cross at a renormalisation group fixed point $`D`$, the product of a critical Ising model and a critical three-state Potts model. Three $`_6`$-invariant operators at this point are the Ising and three-state Potts energy densities $`ϵ`$ and $``$, and their product $`ϵ`$. Since their scaling dimensions are $`1`$, $`4/5`$ and $`9/5`$ respectively, all are relevant and the fixed point is triply unstable. This corrects the approximate (Migdal) renormalisation group result used in, which gave $`D`$ as being doubly-unstable.

We can take this second point a little further, using continuum field theory arguments which are valid in the scaling region around the point $`D`$. The fixed point itself is described by a product of $`c=1/2`$ and $`c=4/5`$ conformal field theories, and nearby points by perturbations of this product $`c=13/10`$ conformal field theory by combinations of the continuum operators $`ϵ`$, $``$ and $`ϵ`$. The first two are anti-self-dual under (2), while the third is self-dual and moves the model away from $`D`$ along the line (3). Minimal models coupled by local operators have received a fair amount of attention (see, for example, ) but this particular instance does not seem to have been studied in any detail. However, Zamolodchikov’s counting argument  can be used to show that the $`ϵ`$ perturbation preserves at the very least conserved charges of spins $`\pm 3`$ and $`\pm 5`$. The perturbed (continuum) theory should therefore be integrable  and we can hope to obtain information about the scaling region of the self-dual line near to $`D`$ via the thermodynamic Bethe ansatz (TBA) technique. This method expresses the finite-volume ‘effective central charge’ $`c(r)`$, $`r=mR`$, of a model with bulk length scale $`1/m`$ confined to a circle of circumference $`R`$ in terms of the solutions $`\epsilon _a(\theta )`$ to a set of coupled integral equations. A candidate system for this case has already been identified: in it was observed that the following TBA system for the functions $`\epsilon _1\mathrm{}\epsilon _7`$

$`\epsilon _a(\theta )`$ $`=`$ $`\delta _{a1}r\mathrm{cosh}\theta {\displaystyle \frac{1}{2\pi }}{\displaystyle \underset{b=1}{\overset{7}{}}}l_{ab}^{[E_7]}\varphi \mathrm{ln}(1+e^{\epsilon _b})(\theta ),`$   
$`c(r)`$ $`=`$ $`{\displaystyle \frac{3}{\pi ^2}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\theta r\mathrm{cosh}\theta \mathrm{ln}(1+e^{\epsilon _1(\theta )})`$ (7)   
predicts an ultraviolet central charge of $`13/10`$, and a small-$`r`$ behaviour of $`c(r)`$ compatible with an expansion in even powers of a coupling $`\lambda `$ to an operator of dimension $`9/5`$. (Here $``$ denotes the convolution, $`fg(\theta )=_{\mathrm{}}^{\mathrm{}}𝑑\theta ^{}f(\theta \theta ^{})g(\theta ^{})`$, $`\varphi (\theta )=1/\mathrm{cosh}\theta `$, and $`l_{ab}^{[E_7]}`$ is the incidence matrix of the $`E_7`$ Dynkin diagram, with $`1`$ labeling the node at the end of the middle-length arm.) This is the behaviour expected of the perturbation of the product Ising and three-state Potts conformal field theory by its $`ϵ`$ operator. (On dimensional grounds, $`m`$ will be related to $`\lambda `$ as $`m\lambda ^5`$, but we do not attempt to find the constant of proportionality here.) Assuming that the TBA is correct, it can be used to extract a non-trivial prediction about the vacuum structure of the perturbed model, using arguments described in section 4.3 of. The system (7) implies the following asymptotic for the ground-state energy $`E(m,R)`$:

$$E(m,R)=\frac{\pi }{6R}c(mR)\frac{\sqrt{6}m}{\pi }K_1(mR)$$
(8)   
where $`K_1`$ is the modified Bessel function of order one and the prefactor $`\sqrt{6}`$ can be interpreted as the Perron-Frobenius eigenvalue of the incidence matrix for single kinks interpolating degenerate vacua of the perturbed theory. These vacua must support a representation of the global $`_6`$ symmetry, and we find that the simplest incidence matrices compatible with both this fact and the eigenvalue implied by (8) are the following:

$$I^a=\left(\begin{array}{ccccccc}\text{0}& \text{1}& \text{1}& \text{1}& \text{1}& \text{1}& \text{1}\\ \text{1}& \text{0}& \text{0}& \text{0}& \text{0}& \text{0}& \text{0}\\ \text{1}& \text{0}& \text{0}& \text{0}& \text{0}& \text{0}& \text{0}\\ \text{1}& \text{0}& \text{0}& \text{0}& \text{0}& \text{0}& \text{0}\\ \text{1}& \text{0}& \text{0}& \text{0}& \text{0}& \text{0}& \text{0}\\ \text{1}& \text{0}& \text{0}& \text{0}& \text{0}& \text{0}& \text{0}\\ \text{1}& \text{0}& \text{0}& \text{0}& \text{0}& \text{0}& \text{0}\end{array}\right)I^b=\left(\begin{array}{ccccc}\text{0}& \text{0}& \text{1}& \text{1}& \text{1}\\ \text{0}& \text{0}& \text{1}& \text{1}& \text{1}\\ \text{1}& \text{1}& \text{0}& \text{0}& \text{0}\\ \text{1}& \text{1}& \text{0}& \text{0}& \text{0}\\ \text{1}& \text{1}& \text{0}& \text{0}& \text{0}\end{array}\right)$$
(9)   
The first is consistent with the coexistence of disordered and $`_6`$ ordered phases, the second with coexistence of $`_2`$ and $`_3`$ ordered phases. This leads us to the first part of our proposal for the resolution of points (i) and (ii) above: we suggest that the entire segment of the self-dual line from $`D`$ to $`C`$, which includes the 6-state Potts transition point $`P`$, lies on a surface of first-order transitions separating disordered and fully-ordered regions, and in the scaling region near to $`D`$ is described by an integrable scattering theory with kink incidence matrix $`I^a`$. (In the massive scaling region near to $`C`$, the vacuum structure comes as the $`N=6`$ case of the results of, and is also given by $`I^a`$.) It is natural to suppose that the matrix $`I^b`$ describes the vacuum structure on the opposite side of the point $`D`$, and support for this idea comes from the following argument.

Consider the behaviour of the model near to the decoupling surface $`J_1=0`$, taking either $`J_21`$ or $`J_31`$. In the first case, the effect will be to freeze out the three-state Potts spins $`\mathrm{\Sigma }_i`$, and all terms $`\delta _{\mathrm{\Sigma }_i\mathrm{\Sigma }_j}`$ can be approximated by $`1`$ in the final line of (5). The result is an effective interaction for the remaining unfrozen spins, equal to $`2(J_3+J_1)[1\delta _{\sigma _i\sigma _j}]`$. Thus the only effect of a non-zero $`J_1`$ in this ‘$`\mathrm{\Sigma }`$-frozen’ region is to replace $`J_3`$ with $`J_3+J_1`$. The line of Ising transitions on the decoupling surface was at $`J_3=J_3^c=\frac{1}{2}\mathrm{ln}[1+\sqrt{2}]`$; we now conclude that a small non-zero $`J_1`$ in the region $`J_21`$ simply shifts this line to $`J_3=J_3^cJ_1`$. A similar argument shows that the line of three-state Potts transitions, situated at $`J_2=J_2^c=\frac{2}{3}\mathrm{ln}[1+\sqrt{3}]`$ on the decoupling surface, is shifted to $`J_2=J_2^cJ_1`$ in the region $`J_31`$ where the $`\sigma `$ spins are frozen. Finally, duality can be used to see that in the opposite regimes $`J_21`$ and $`J_31`$, the critical lines move in the opposite senses when $`J_1`$ shifts away from zero. This line of argument has nothing to say directly about the central region where the Ising and three-state Potts spins are both near their critical points, but the simplest hypothesis, also consistent with the continuum results of the last paragraph, is to continue the lines in from the asymptotic regions as in Fig. 1.

Much of our numerical work has been devoted to the confirmation of this picture. Before we describe this, we pause to introduce a third set of coordinates on the thermodynamic phase space. First, notice that the parameter $`K`$ defined by

$$K(x_1,x_2,x_3)=(x_1x_2x_3)/\mathrm{\Delta }$$
(10)   
is mapped into itself by duality: $`\stackrel{~}{K}=K`$. Therefore surfaces of constant $`K`$ in phase space are mapped into themselves by duality; $`K=0`$ is the decoupling surface. Points with $`K>0`$ lie ‘below’ this surface, on the same side as the 6-state Potts point $`P`$. It is easy to check that the two parameters

$$y_2=\mathrm{ln}\frac{1+2x_2}{\sqrt{36K}}y_3=\mathrm{ln}\frac{1+x_3}{\sqrt{24K}}$$
(11)   
have the simple duality transformation rule

$$\stackrel{~}{y_2}=y_2\stackrel{~}{y_3}=y_3$$
(12)   
It is convenient to use $`(K,y_2,y_3)`$ as coordinates in phase space: on any fixed-$`K`$ surface the point $`(K,0,0)`$ is the self-dual point and duality is just the reflection in this point. The origin $`(0,0,0)`$ is the fixed point $`D`$ on the decoupling surface, and the 6-state Potts transition point $`P`$ is at $`(K_P,0,0)`$ with $`K_P=(\sqrt{6}1)^2/250.084`$. On the decoupling surface, $`y_2`$ and $`y_3`$ are functions of the Ising and Potts couplings respectively.

A cluster algorithm appears to be the ideal choice to simulate the relevant regions of the phase diagram, since it is reasonable to expect large correlation lengths even when the transitions are only first-order. Since cluster algorithms are especially simple to implement for $`_2`$ and $`_3`$ models, we exploited the possibility of writing the Hamiltonian as in (5): the algorithm we used performs alternate cluster updates of the $`\sigma _i`$ and $`\mathrm{\Sigma }_i`$ variables, where those which are not being updated provide effective, site-dependent couplings for the others. For example, suppose we are updating the $`_2`$ variables: the effective coupling to be used on the link $`<ij>`$ is then $`J_3+J_1\left(\mathrm{\Sigma }+\mathrm{\Sigma }^1\right)/2`$ with $`\mathrm{\Sigma }=\mathrm{\Sigma }_i/\mathrm{\Sigma }_j`$, while when updating the $`_3`$ variables the effective coupling is $`J_2+J_1\sigma `$ with $`\sigma =\sigma _i/\sigma _j`$. Similar algorithms were introduced in  for the Ashkin Teller model, defined as two coupled Ising models.

The fact that the effective couplings are site-dependent does not pose a problem for the cluster updates. However the effective $`_2`$ and $`_3`$ couplings can become negative, i.e. antiferromagnetic. This will happen for the effective $`_2`$ coupling whenever $`J_3+J_1<0`$ or $`J_3J_1/2<0`$, and for the $`_3`$ coupling whenever $`J_2+J_1<0`$ or $`J_2J_1<0`$. Negative couplings on some links can in turn lead to frustrations, which in principle can make the cluster algorithm highly ineffective. However it is easy to convince oneself that the $`_3`$ model is never actually frustrated: for every configuration of the $`_2`$ spins there exists a configuration of the $`_3`$ ones such that all the links are satisfied. Therefore the possibility of frustrations exists only for $`J_1<J_3`$ or $`J_1>2J_3`$. These relations are never satisfied in the regions we considered.

We explored several fixed-$`K`$ surfaces and mapped the various phases and their boundaries using the Binder cumulants method . In our case, we define two cumulants, one for each order parameter:

$$Q_\sigma =\frac{3}{2}\frac{1}{2}\frac{m_\sigma ^4}{m_\sigma ^2^2}Q_\mathrm{\Sigma }=2\frac{m_\mathrm{\Sigma }^4}{m_\mathrm{\Sigma }^2^2}$$
(13)   
where $`m_\sigma `$ and $`m_\mathrm{\Sigma }`$ are the $`_2`$ and $`_3`$ magnetizations per site. For each surface we measured the $`_2`$ and $`_3`$ Binder cumulants on a grid of points around the self-dual point $`(K,0,0)`$ for two different lattice sizes $`L_1\times L_1`$ and $`L_2\times L_2`$. For every fixed value of $`y_2`$ we estimated the value of $`y_3`$ where each Binder cumulant remains constant when the lattice size is increased, i.e. the transition point. This allows us to determine a set of points belonging to the transition lines for the $`_2`$ and $`_3`$ variables. These are the points plotted in Figs. 2,3,4 for the $`K=0.06,0.04,0.08`$ surfaces respectively. The consistency of the procedure was checked by repeating it with the roles of $`y_2`$ and $`y_3`$ inverted. Finite-size effects are signaled by violations of the duality symmetry, that is, in our coordinates, symmetry under reflection in the origin. Small violations are visible in our figures; we checked that they become smaller when the lattice sizes are increased. These effects turned out to be much larger in the $`K<0`$ region. In fact we observed that to control them satisfactorily in the $`K>0`$ region lattice sizes as small as $`L_1=15`$, $`L_2=30`$ sufficed, while for $`K<0`$ we had to use $`L_1=45`$, $`L_2=60`$.

The numerical results confirm the qualitative picture shown in Fig. 1, and our claim that the decoupling point $`D`$ on the self-dual line marks a change from the coexistence of disorder and $`_6`$ order to the coexistence of $`_2`$ and $`_3`$ order. The reduced length of the first-order segment at $`K=0.04`$ is consistent with scaling predictions, though we are not close enough to the point $`D`$ to see a complete collapse of data. In conclusion, we have proposed a significant modification to the previously-accepted phase diagram of the $`_6`$ spin models, and this has been supported by a detailed numerical study.

Acknowledgements –  We would like to thank John Cardy, Michele Caselle, Eytan Domany, Ferdinando Gliozzi, Martin Hasenbusch, Klaus Pinn and Jean-Bernard Zuber for useful conversations. The work was supported in part by a TMR grant of the European Commission, contract reference ERBFMRXCT960012. PED thanks the UK EPSRC for an Advanced Fellowship, and RT thanks SPhT Saclay for hospitality.