Dualities and the phase diagram of the $p$-clock model

A new"bond-algebraic"approach to duality transformations provides a very powerful technique to analyze elementary excitations in the classical two-dimensional XY and $p$-clock models. By combining duality and Peierls arguments, we establish the existence of non-Abelian symmetries, the phase structure, and transitions of these models, unveil the nature of their topological excitations, and explicitly show that a continuous U(1) symmetry emerges when $p \geq 5$. This latter symmetry is associated with the appearance of discrete vortices and Berezinskii-Kosterlitz-Thouless-type transitions. We derive a correlation inequality to prove that the intermediate phase, appearing for $p\geq 5$, is critical (massless) with decaying power-law correlations.


Introduction
In this article we investigate, via the use of dualities, two-dimensional (D = 2) classical systems, such as the XY and clock models [1,2], that display Berezinskii-Kosterlitz-Thouless (BKT)-type transitions [3,4,5,6]. BKT transitions, notably characterized by essential singularities in the free energy, emerge in many physical situations including screening in Coulomb gases [7], surface roughening [8], melting in D = 2 solids [4,9], and many other classical and quantum problems, such as deconfinement in D = 3 + 1 lattice gauge theories [10,11,12,13]. We study these models by invoking a bond-algebraic approach that we have recently developed [12,13,14,15]. Within our duality-based method, one relates singularities in the free energy at one temperature (or coupling constants) to those at a dual temperature (or dual coupling constants).
Specifically, we investigate exact dualities of the D = 2 XY and p-clock models, and exploit those dualities to unravel their phase structures. Those transformations are exact even for finite systems after appropriate boundary terms are included. It is noteworthy that unlike nearly all of the analytical work done to date, our dualities do not rely on the approximation scheme of Villain [16,10,17,18], yet they can be related to the exact dualities of the Villain model in appropriate limits. Furthermore, our analysis leads to exact dualities for general p-clock models and yields a better understanding of the appearance of two transitions in systems with p ≥ 5 states (the XY model with only one transition is recovered in the p → ∞ limit). By fusing our duality results with the Peierls argument, we will be able to (1) prove that p ≥ 5 clock systems can be made to be self-dual; (2) prove by a Peierls argument that there exists a lower ordering temperature T (1) ∼ 1/p 2 , associated to domain-wall excitations, below which the global Z p symmetry is broken; (3) demonstrate that a second transition occurs at a temperature T (2) ∼ O(1) when p ≥ 5 (in the self-dual case, it follows that if T (1) is not the self-dual temperature T * , then there must necessarily be a second phase transition at T (2) of an identical character); and (4) characterize the nature of the topological excitations, and further explain that at p = 5 a new type of topological excitation, with an associated discrete winding number, appears. Our considerations suggest that these discrete vortices may proliferate above the temperature T (2) . We will also (5) determine an analytic expression for the self-dual temperature T * , which relates and clarifies temperature scales discussed in [19], and most importantly, (6) establish the non-Abelian polyhedral symmetry group P (2, 2, p) of the p-clock and related models, and explicitly unveil the U(1) continuous symmetry that emerges when p ≥ 5. Indeed, the latter is intimately tied to the existence of the BKT transition. Finally, (7) we derive a correlation inequality to prove that the intermediate phase, appearing for p ≥ 5, is critical (massless) with decaying power-law correlations.
Despite several analytic [10,20,21] and numerical calculations [22,19,23], the precise nature of the two phase transitions (p ≥ 5) is not completely understood. It was proven [21] that for large enough p, clock models exhibit a BKT-type transition (actually, it has only been proved that there exists an intermediate critical phase with power-law correlations). The question whether that still holds when p = 5 remains open [22]. By relying on exact results, we shed light on the character of the BKT transitions in these systems. We will relate the BKT transition to a continuous U(1) emergent symmetry of an usual type. Although BKT transitions are often discussed in terms of specific anomalous exponents and jumps in the helicity modulus, we will not address such non-universal issues.
Our treatment of classical dualities is based on a new approach developed in Refs. [14,13] that relies on the transfer matrix or operator formalism [11]. In statistical mechanics, two models a and b are dual if their partition functions Z a = tr [T N a ] and Z b = tr [T N b ] are related as (N is the linear size of the system in D = 2) with A some analytic function of the set of couplings K of model a, and dual couplings K * . In principle, Eq. (1) establishes an extremely broad relationship that could be achieved through many transformation schemes, including the standard one based on taking the Fourier transform of individual Boltzmann weights [17,24]. However, it was discovered in Refs. [14,13] that low-temperature(strong coupling)/high-temperature(weak coupling) dualities correspond to a unitary equivalence of transfer matrices or operators, T a and T b , with U d a unitary operator. This observation is extremely insightful because there is a simple and systematic way to look for unitary equivalences between physical operators, based on the notion of bond algebra, or algebra of interactions [14,13]. The outline of this article is as follows. In Section 2 we define the classical XY model, and then in 2.1 establish its transfer operator. In Section 2.2 we discuss the form of the exact one-dimensional quantum analogue whose partition function is that of the D = 2 classical XY model with coupling constants K 1 and K 2 . In the limit of large coupling K 2 along columns, the quantum model is the O(2) quantum rotor model. In Section 2.3, we establish the duality of the D = 2 XY model to a solid-on-solid-like and also to a lattice Coulomb gas-like models and, moreover, determine the disorder variables. These dualities do not rely on the Villain approximation scheme but are exact dualities obtained by our bond-algebraic method [13].
We next proceed to analyze in Section 3 the p-clock model [1]. This model provides a particular controlled limit to the XY model (the p → ∞ limit). We replicate the same steps undertaken in the analysis of the XY model of Section 2, but now the Weyl algebra [25] and the theory of circulant matrices [26] play a key role. We construct in Section 3.1 its transfer matrix, and proceed in 3.2 to establish the corresponding one-dimensional quantum Hamiltonian, that is not self-dual for p ≥ 5. We study the dualities of these systems in Section 3.3. The system is exactly self-dual for p = 2, 3, 4, and becomes approximately self-dual for K 2 ≫ K 1 when p ≥ 5. In Section 3.4 we introduce a variant of the classical p-clock model that is exactly self-dual for all p. We examine, in Section 3.5, the exact and emergent symmetries of these systems and, notably, unveil the U(1) symmetry that emerges when p ≥ 5. This continuous emergent symmetry is responsible for the existence of the intermediate critical (massless) phase.
Finally, in Section 4 we utilize our previous findings to better understand the phase diagram of the p-clock model. Here we present an analytic expression for the self-dual temperature T * , an important scale in the problem, and advance a Peierls argument. We also introduce a topological invariant, that we call the discrete winding number k, to unravel the nature of the topological excitations. Starting at p ≥ 5 a new type of topological excitation appears with a non-zero value of k that one may call discrete vortex, and which is responsible for the phase transition to a disordered state. Domain-wall topological excitations are key at low-temperatures and their energy cost depends on p, and on the relative spin configurations (except for 2 ≤ p ≤ 4). By using both duality and energy-versus-entropy balance considerations we show that the transition from the critical to the disordered phase scales as T (2) ∼ O(1), while the one from the broken Z p symmetry to the critical phase goes as T (1) ∼ 1/p 2 for large p. We derive a correlation inequality, i.e., show that the two-point correlation function G is a monotonically decreasing function of temperature, allowing us to prove that, for p ≥ 5, the intermediate phase is critical (massless) with decaying power-law correlations. The appendices provide technical developments including a duality of the XY model to qdeformed bosons, illustrating the key physical difference between compact and non-compact degrees of freedom.

The XY model: A paradigm of BKT phenomenology
The D = 2 classical XY model is the paradigmatic example of a system displaying a BKT transition at a finite temperature T (2) BKT > 0. This model, also known as planar rotator or planar O(2), consists of an N × N array of classical two-component spins S r located at the vertices r = i e 1 + j e 2 (i, j being integers) of a square lattice with unit vectors e µ , µ = 1, 2, as indicated in Fig. 1. Its partition function is where the spin S r = S x r e 1 + S y r e 2 , the coupling K µ = βJ µ is the product of the inverse temperature β = 1/k B T and the exchange coupling J µ , and h is a temperature-rescaled external magnetic field. Fixing the magnitude of the spin variable S 2 r to one allows us to re-write the partition function as where the continuous angle variables θ r take values in the interval θ r = θ i,j ∈ [0, 2π), i.e., it is a compact variable. We assume, without loss of generality, that h = h e 1 . The sum over configurations represents an integral In the remainder of this article, we will concentrate on the case with zero external magnetic field, i.e., h = 0. The XY model displays a (global) continuous U(1) symmetry, which amounts to the invariance of the model under a simultaneous rotation of every spin in the lattice by the same angle. In low dimensional systems with continuous symmetries, such as the XY model above, long-range order is more fragile. Thermal fluctuations may induce instabilities with the end result that long-range order is actually non-existent in D ≤ 2 dimensions. Spin-wave excitations are responsible for destroying such an order. The Mermin-Wagner theorem formalizes this qualitative picture. In the context of the D = 2 XY model, the theorem states that this system does not display spontaneous magnetization at finite temperatures. Figure 1: The two-dimensional classical XY model. On each vertex r = i e 1 + j e 2 of the square lattice there is a classical spin S r = S x r e 1 + S y r e 2 of magnitude 1, i.e., S r · S r = 1, and S x r = cos θ r , S y r = sin θ r with θ r ∈ [0, 2π). Nearest neighbor spins interact with an exchange constant value J 1 or J 2 depending on the spatial direction.
A common physical mechanism behind the formal proofs for both classical and quantum versions of the XY model, can be found in Ref. [11].
A phase transition is said to occur whenever a thermodynamic function of the system under study displays a non-analyticity. The latter may occur even when the ground state is unique, so that there is no spontaneous symmetry breakdown. This is the case for the D = 2 XY model, that is known to have a special phase transition at a finite, non-zero temperature T (2) BKT . This BKT transition is characterized by an essential singularity in the free energy and correlation length at T (2) BKT . If T > T (2) BKT , the correlators of the XY model decay exponentially with distance, as is typical of a disordered, paramagnetic phase. In the low-temperature phase, T < T (2) BKT , the correlators decay algebraically with distance, just as if every temperature below T (2) BKT represented an ordinary critical point. The fact that this power-law behaviour extends over the finite temperature range 0 < T < T (2) BKT , without long-range order, is known as quasi-long-range order.

A transfer operator for the XY model
In this section we set up a transfer operator for the XY model, in preparation for the detailed study of its duality properties and symmetries. We assume open boundary conditions in the e 1 -direction and periodic ones in the e 2 -direction.
Since we are considering the XY model on a square lattice of size N × N we will need the operators satisfying the following commutation relations The eigenstates of the unitary operators e ±iθ i , The plane wave eigenstates |n i of L z,i form an orthonormal basis of the Hilbert space of square integrable functions on the cirle L 2 (U(1)) = H i , and are related to the states |θ i via θ i |n i = e iθ i n i / √ 2π. On the one hand, e ±iθ i represent position operators for the spin at site i, since their simultaneous eigenstates |θ i specify one, and only one point on the unit circle. On the other hand, L z,i represents their canonically conjugate momentum, the infinitesimal generator of translations This last equation follows from Eq. (7), since The product states, that are simultaneous eigenstates of all the position operators e ±iθ i , form an orthonormal basis of the total Hilbert space We have now all the ingredients needed to write down a transfer operator for the XY model. Consider for concreteness the following row-to-row (j to j + 1) matrix elements of the transfer operator and the diagonal operator Both matrices are defined in the basis of states introduced in Eq. (12). It is straightforward to check that if we set T XY ≡ T 2 T 1 , then recovers the partition function for the XY model of Eq. (3), provided that we set the external magnetic field h to zero. Next, we rewrite the operators T 1 , T 2 in terms of the operators introduced in Eq. (7). The result reads as can be checked by taking matrix elements of T 2 T 1 in the basis of Eq. (12). Notice that T 1 factors into a product of two-body operators that involves only nearest neighbors, while T 2 factors into a product of one-body operators. This important simplification is a direct reflection of locality. The relevant symmetries of the classical XY model translate into unitary transformations that commute with T XY ≡ T 2 T 1 . Besides the obvious geometrical symmetries of the lattice, T XY commute with two operators that represent internal, global symmetries. The continuous global U(1) symmetry under global rotations of the classical spin direction θ r → θ r + α, ∀r, guarantees that [L z , T XY ] = 0, where L z = N i=1 L z,i is the total angular momentum. There is also a discrete symmetry C 0 = N i=1 C 0i that is alluded to, somewhat inaccurately, as "charge conjugation" [27]. The operator C 0i acts on position eigenstates as C 0i |θ i = |2π −θ i , and on angular momentum eigenstates as C 0i |n i = | − n i . Thus, C 2 0i = 1 and C † 0i = C 0i . Since C 0 does not commute but rather anticommutes with L z , we notice that the full group of internal symmetries of the XY model is non-Abelian.

Hamiltonian form of the XY model
Models that can be written in terms of an Hermitian transfer matrix or operator, such as the XY model, can be translated into quantum-mechanical problems [11] by simply defining a quantum Hamiltonian according to While this is a powerful tool, often used in the literature, its actual value is diminished by the technical problem of computing ln(T XY ) and, perhaps more importantly, because H XY turns out to be a highly non-local operator. The standard way out of this difficulty is to make approximations that solve both of these problems. The qualitative picture that emerges is intuitively appealing but problematic if the approximations are not reasonably controlled. We will not determine H XY in closed form, but rather we will compute in closed form, and then exploit the Baker-Campbell-Haussdorf (BCH) formula to obtain the expansion Since T XY = T 2 T 1 , as defined in the previous section, is not Hermitian, we have to set H XY = (H + H † )/2. This only affects terms quadratic and higher order in the commutators. We will also study the conditions for the nondiagonal (kinetic) part of the Hamiltonian H 2 to reduce to the intuitively appealing form H 2 ∝ i L 2 z,i /2. Referring back to the operator form of T 1 , T 2 , Eq. (16), we see that it is straightforward to compute H 1 . Since T 1 is already in diagonal form, we obtain On the other hand, T 2 is not diagonal. Computing H 2 is further simplified by the fact that T 2 factors into the product of N commuting one-body operators. It follows that Next we notice that, since e −iθ 1 L z,i e −iθ 2 L z,i = e −i(θ 1 +θ 2 )L z,i , H 2,i should be of the form (see Appendix A) By combining this expression with Eq. (21), we get an equality between functions of L z,i that can be evaluated in that operator's diagonal basis. The result is an infinite set of equations that one can use to determine a K 2 (θ) by taking the Fourier transform of the left-hand side. The modified Bessel functions of the first kind and of integer order n, I n (K 2 ), satisfy [28] e K 2 cos θ = n∈Z I n (K 2 )e iθn = I 0 (K 2 ) + 2 ∞ n=1 I n (K 2 ) cos(θn), after using the relation I −n (K) = I n (K) [28]. It follows from Eq. (23) that a K 2 (2π − θ) = a K 2 (θ), so that we can make the substitution e −iθL z,i → cos(θL z,i ) in Eq. (22). Then, we can Taylorexpand the cosine function to get with a m (K 2 ) = (−1) m 2π 0 dθ θ 2m a K 2 (θ)/(2m)!. This equation provides a very convenient representation of H 2,i , especially if we are allowed to discard terms beyond m = 1. To investigate this possibility, consider Eqs. (21) and (25), and evaluate those expressions in the L z,i 's diagonal basis to get ∞ m=0 a m (K 2 ) n 2m = ln(2πI n (K 2 )).
For large K 2 , the functions I n (K 2 ) have the following asymptotic expansion [28] with (28) Notice that this can be rearranged into an expansion in n 2 that can be compared to the left-hand side of Eq. (26). Keeping, for each m, only the leading order in 1/K 2 , and expanding the logarithm accordingly (ln(1 + x) ∼ x), we obtain ln(2πI n (K 2 )) ∼ ln 2π Comparing with Eq. (26), so that a m+1 /a m ∼ −1/(2(m + 1)K 2 ). In summary, in the large K 2 limit, which is the one-dimensional O(2) quantum rotor model.

Duality of the XY model without the Villain approximation
We now exploit the detailed understanding we have gained on the exact operator structure of the XY model to look for its dual representations. The standard approach to the dualities of the XY model starts by replacing it with the Villain model (see Appendix C), then mapping the latter to the solid-on-solid (SoS) model, and finally mapping the SoS model to a lattice Coulomb gas [11]. In contrast, in this section we establish directly exact dual representations of the XY model. In Appendix B, we establish a duality to a q-deformed boson Hamiltonian which illustrates the fact that non-canonical bosons need to emerge because of the compact nature of the degrees of freedom of the XY model.
Our methodology starts with the transfer operators T 1 , T 2 introduced in Eq. (16). The algebra of interactions, or bond algebra in the language of Refs. [12,13,14,15], underlies their basic structure. In this case, this is the von Neumann algebra A XY generated by the bonds Notice that T 1 , T 2 ∈ A XY , since these operators are expressible as sums of products of the bonds listed in Eq. (32). A XY reflects the interactions present in the XY model and is at the same time easy to characterize in terms of relations. Then we can look for other dual realizations A D XY that are isomorphic images of A XY . By the general properties of von Neumann algebras, it must be that provided both algebras act on state Hilbert spaces of the same dimensionality. This is all we need to establish a duality for the XY model. The dual partition function is determined from the dual transfer operator The goal of this section is to look for a dual representation of the XY model in terms of integer-valued degrees of freedom, so we can expect the dual bond algebra A D XY to act on the state space N i=1 L 2 (Z). Let us introduce the states |n = N i=1 |n i , and the operators X i , R i that satisfy the algebra Then, the operators generate an isomorphic dual representation A D XY of the bond algebra of the XY model. The isomorphism Φ d connecting the two bond algebras is illustrated in Fig. 2, for N = 3 sites. The resulting dual transfer operators are The next and last step is to compute the dual partition function in the product basis states of Eq. (32). T D 2 is already diagonal in that basis, leading simply to where The evaluation of n ′ |T D 1 |n factors into the computation of the one-body matrix elements, This is simplified by noticing that the Fourier transform operator of Eq.
thus putting the exponential in diagonal form. It follows that If we now put all the pieces together, we arrive at the conclusion that we have obtained the exact duality (43) This illustrates a typical characteristic of dualities: the coupling K 1 (K 2 ) in the e 1 (e 2 )-direction of the XY model regulates the interaction in the orthogonal e 2 (e 1 )-direction of the dual model.
It is standard to argue that the Villain model is an excellent approximation to the XY model, specially at low temperatures. As shown in Appendix C, the Villain model is dual to the SoS model. Thus, it must be that the SoS model is approximatly related to Z D XY defined by the right-hand side of Eq. (43), at least for low temperatures. Consider then the limit of large K 1 , K 2 (i.e., low temperatures). We can then use the asymptotic expansion of Eq. (27) to obtain an asymptotic form of the dual potential of Eq. (39), where c(K) is independent of n and can be computed from Eqs. (39) and (27). It follows that (45) to the same level of approximation. Thus we have recovered the well known result that the XY model at very low temperatures (strong coupling) is well represented by the (approximately dual) SoS model at very high temperatures (weak coupling).
The action of the duality of Eq. (35) can be extended to act on the operator e iθ N as Φ d (e −iθ N ) = R N . It follows that Φ d generates the following set of N dual variables (below we distinguish a dual variable by an overtilde), The dual variables satisfy the algebra of Eq. (7), confirming that Φ d defines an algebra isomorphism. Since this is also the algebra of the variables that should be compared to n m+r,n+s n m,n = tr (T The classical dual variables e −iθm,n are difficult to define directly, but they are well defined in the sense that any correlator in the ensemble Z D XY can be computed by a straightforward generalization of (47).
The duality to a lattice Coulomb gas is of a very special nature (Poisson duality) and quite different from every other duality discussed in this article (or in literature on dualities in general). Its general features are discussed in Ref. [13]. Here, we briefly summarize the exact Coulomb gas-like dual model for the exact dual partition function Z D XY computed above. According to Ref. [13], the lattice gas representation of Z D XY is defined by with interaction energy E D defined through the global Fourier transform The interaction energy E D can be determined in closed-form in the limit in which Z D XY reduces to the SoS model, and corresponds to the standard lattice Coulomb gas result [11]. We see from Eq. (51) that Poisson dualities are in general only of practical use for models with Gaussian energy functionals.

The p-clock model: A close relative of XY
The p-clock model, also known as the vector Potts or Z p model, represents a hierarchy of approximations to the XY model as a function of the positive integer p, and is a test ground for rich critical behavior. In D = 2 dimensions, the configurations of the classical p-clock model are described by a set of p discretized angles θ r The partition function is given by This is, in appearance, identical to the classical XY model (h = 0), except for the essential fact that the degrees of freedom {θ r } are now discrete and countable. In the large p limit (p → ∞), however, one supposedly recovers the XY model. Since the p points e iθr close a Z p subgroup of U(1) (the group of p'th roots of unity), the p-clock model manages to provide an approximation to the XY model that features a finite number of states per site, without sacrificing the XY's natural group structure. Also like the XY model, the p-clock has interesting but hard to uncover duality properties [20]. We will address this problem by the same methods applied to the XY model. In fact, we will follow the discussion of the XY model as closely as possible, to highlight the connections between the two models.

A transfer matrix for the p-clock model
To introduce a transfer matrix for the p-clock model, we need to define a suitable Hilbert space and a set of basic kinematical operators. Let Z p be defined on an N × N square lattice with open boundary conditions on the e 1 -direction and periodic ones on the e 2 -direction. The states on each site i = 1, · · · , N of a row can be described by orthonormal vectors such that s i represents the discrete angle 2πs i /p. They span the state space |s i for elements of the product basis of H p , then we can define a matrix related to any pair of adjacent rows j, j + 1, and the diagonal matrix These definitions guarantee that The degrees of freedom of the p-clock model (at any one site of the lattice) can take any value out of a discrete, equidistant subset of points of the unit circle. To proceed in analogy to Section 2.1, we need to introduce position operators and their conjugate momenta in this discrete setting. The formalism that emerges was used extensively by Schwinger in his work on the foundations of quantum mechanics [25]. In what follows, we consider only one site (one degree of freedom), for the sake of clarity. We will consider all N sites again near the end of the section.
It is easy to restrict the position operators e ±iθ used for the XY model to the subset of configurations available to a clock handle in the p-clock model. The result is the operator U satisfying with ω ≡ e i2π/p representing a pth root of unity. The position operator U and its Hermitian-conjugate U † satisfy UU † = 1 = U p . The momentum operator V conjugate to U rotates any state counter-clockwise to its nearest-neighbor Momentum and position operators are represented, in (p × p) matrix form, as It follows that V † implements a clockwise rotation, and that The fundamental algebraic relation follows directly from the definitions of U and V .
As is well known from quantum mechanics, the ordinary position operator x and its conjugate momentum operatorp are related by a Fourier transform F , a unitary transformation in the space of wave functions. Essentially the same holds for the operators U, U † and V, V † . The appropriate unitary transformation in this context is the discrete Fourier transform F , that in matrix form reads This is also known as Schur matrix [29]. It follows that [25] and so the eigenvectors of V ,s = 0, 1, · · · , p − 1, are easily determined via a Fourier transform of the eigenvectors of U.
In the mathematical literature V is known as the fundamental circulant matrix. This is so as it generates the algebra of circulant matrices [26] (meaning that any circulant matrix C is of the form C = p−1 m=0 a m V m , a m ∈ C). Together, U and V generate the full algebra of (p × p) complex matrices [25], that we continue to call the Weyl group algebra, to emphasize that we are working with a distinguished set of generators. This shows that they constitute a convenient basis set of kinematic operators, because we can write any other operator in terms of them.
We need to re-introduce the row spatial index i to apply the technology just developed and rewrite the transfer matrices of Eqs. (56) and (55) , N, will be our basic set of operators. They commute at different sites, satisfy the relation (60) at any one site i, and act on the state space H p = N i=1 H p,i . One then obtains This last expression for T 2 follows from the fact that ). It should be compared to the analogous expression for the continuum circle, Eq. (16).

Hamiltonian form of the p-clock model
In this section we compute the Hamiltonian form of the p-clock model following the strategy of Section 2.2. We start by computing H µ = − ln T µ , µ = 1, 2, with T 1 , T 2 as defined in Eq. (64).
Since T 1 is diagonal, we can write As explained in Appendix A, we can solve this equation to obtain with Then, the Hamiltonian H p for the p-clock model follows provided we truncate the BCH expansion of ln T p to linear order (see the discussion in Section 2.2). We notice for future reference that the discrete As discussed in Appendix A, the coefficients a m (K 2 ) have simple asymptotic forms in the limit K 2 → ∞. The corresponding approximation to H p reads with and (see Eq. (A.12)) Equation (71) shows a boundary term (−K 1 (U N + U † N )/2) not present in Eq. (69), and that we include to make this approximation to H p exactly self-dual [13].
The approximation made in going from Eq. (69) to Eq. (71), that keeps only V † i and V †(p−1) i = V i , is reminiscent of the one introduced in Section 2.2 based on Eq. (26), whereby we replaced the operator cos(θL z,i ) for the simpler L 2 z,i in Eq. (31). Indeed, the two aproximations coincide in the p → ∞ limit. A simple way to see this is to notice that we can realize the operator V † i directly in the Hilbert space of the XY model as V † i → e −i2πL z,i /p . Then, in the limit p → ∞, Figure 3: The duality isomorphic mapping Φ d of Eq. (76), for N = 3 sites.

Dualities of the p-clock model
The dualities of the p-clock model appear as isomorphic representations of the bond algebras associated to the transfer matrices defined in Eq. (64). The bond algebra A p generated by is simple to work with and adequate to our purposes. It has a dual (isomorphic) representation A D p generated by the same bonds listed in Eq. (75), except for V N , V † N that have to be removed from the set of generators, and replaced by U 1 , U † 1 . The duality isomorphism Φ d : together with the corresponding Hermitian-conjugate entries (since it must Fig. 3, and should be compared with the duality of Section 2.3 for the XY model, illustrated in Fig. 2. of the p-clock model that follows from Eq. (76) is defined by the dual transfer matrices Clearly, the p-clock model is not self-dual for arbitrary p and arbitrary couplings. However, the model is approximately self-dual in the extreme anisotropic limit with K 2 ≫ K 1 , and it is exactly self-dual for p = 2, 3, 4, and any coupling. We study these aspects of the p-clock model in the next sections and in Appendix D.
We mention, in closing, that there is a p-state model that approximates the p-clock in the same sense in which the Villain model approximates the XY model. This Z p Villain model [10] is exactly self-dual for any p, but is otherwise quite different from the self-dual p-clock model to be introduced next.

Self-dual classical p-clock model
In this section we introduce a classical model Z sdp that we call the self-dual p-clock. It is closely related to the p-clock model (identical for 2 ≤ p ≤ 4), yet it is exactly self-dual for any value of p, and has the distinct advantage over the Z p Villain model [10] that its transfer matrix is remarkably simple. To define the self-dual p-clock, we introduce the transfer matrices (78) Then, Z sdp [a 0 , a 1 , a 1 ], and where a 0 and a 1 are free parameters of the model to be determined, for instance, by the requirement that the approximation be as good as possible for arbitrary K 2 . Z sdp is self-dual due to the existence of a unitary transformation that maps This fact results from a bond-algebraic analysis, but we omit the details which can be found in Ref. [13]. It follows from Eq. (80) that the self-dual line is specified by K 1 = 2a 1 . The next issue then is to understand the structure of Z sdp in terms of classical variables. On one hand, it is clear that the interaction energy in the e 1 -direction is still of the form K 1 cos(2π(s i+1,j − s i,j )/p). On the other hand, the interaction energy in the e 2 -direction, u(s i,j − s i,j+1 ), is determined by the relation is exactly self-dual under the exchange 2a 1 ↔ K 1 .
One can see directly from Eq. (82) that u(p − m) = u(m). It follows that  A.13). The point to notice is that to make the p-clock self-dual, we need to add higherorder harmonics (terms cos(2πm(s i,j − s i,j+1 )/p), with m = 2, · · · , [p/2]) to the basic cosine interaction. We show next that Z sdp becomes a very good approximation to the standard p-clock model in a suitable limit.
A comparison of the self-dual p-clock to the standard p-clock model shows that the latter has an approximate self-duality for p ≥ 5. As explained in Appendix A, in the limit in which K 2 is large, Eq. (79) becomes almost exact, with a 0 ≈ K 2 , a 1 ≈ e K 2 (cos 2π (see Eq. (A.12)), so that (see Eq. (72)), with dual couplings We emphasize that this approximate self-duality, in the extreme anisotropic limit, is valid for any value of p. We consider exact self-dualities for the particular cases p = 2, 3, 4 in Appendix D.

Exact and emergent symmetries of the p-clock model
Non-Abelian, discrete symmetries. The representation of the transfer matrix T p = T 2 T 1 , Eq. (64), is very convenient for understanding the internal, global symmetries of the p-clock model. It is apparent that the model has an Abelian Z p symmetry, but, as it turns out, its full group of symmetries is considerably larger and non-Abelian, provided p ≥ 3. To prove this we will show that there are two Hermitian operators that commute with T p and satisfy These relations show that, if p ≥ 3, C 0 and C 1 generate a unitary representation of the so called polyhedral group P (2, 2, p) [30] of order 2p, and so the group of internal symmetries of the p-clock model is at least as big as this non-Abelian group. Notice that C 0 C 1 ≡Q, known as Z p charge, generates a Z p subgroup of P (2, 2, p). This is the standard Abelian symmetry of the p-clock model that gets broken in the low-temperature ordered phase (see Section 4). It becomes a U(1) symmetry in the limit p → ∞, and corresponds to the usual continuous symmetry of the classical XY model θ r → θ r + α, with α an arbitrary real number. As C 0 and C 1 are products of one-site operators, let us focus on a single site i for now. We define the operators C 0i and C 1i by specifying their action on the basis of Eq. (54), The arithmetic in these definitions is modular, mod(p). For example, if p = 5, then C 0i |0 = | − 0 = |0 , C 0i |1 = | − 1 = |4 , and so on. Keeping this in mind, one can check that C 0 and C 1 are Hermitian, C † 0i = C 0i , C † 1i = C 1i , and satisfy the relations listed in Eq. (90) for C 0 and C 1 . In particular, as C 0i C 1i = V i , (C 0i C 1i ) p = 1. If p = 2, then C 0 = 1, and C 1 generates the Abelian Z 2 symmetry of the Ising model which is different from the non-Abelian P (2, 2, 2).
A routine calculation shows the action of C 0i , C 1i on the discrete position U i and momentum V i operators, It is easy to check that C 0 and C 1 commute with T p . The operator C 0 is known in the literature as the "charge-conjugation" operator [27]. However, as we alluded to earlier, this is something of a misnomer. Geometrically speaking, C 0 is the exact analogue of the parity operator P|x = | − x on the real line.
In fact, C 0 is related to the discrete Fourier transform asF 2 = (F † ) 2 = C 0 , just as its counterpart on the real line F is connected to the parity operator as F 2 = P. We wish to emphasize that these non-Abelian symmetries are shared by a large number of classical and quantum p-state models besides the p-clock, including the self-dual p-clock introduced in Section 3.4 and the Z p Villain models.
Emergent U(1) symmetry. For p ≥ 5 the discrete charge symmetryQ gets enhanced into a continuous U(1) symmetry. In reality, this is not an exact symmetry it is an emergent one [31], but it is essential to establish the intermediate BKT critical (massless) phase (see Section 4). Let us derive this emergent symmetry.
Given the generators of the SU(2) algebra in the spin S = (p − 1)/2 , and S − = (S + ) † , one may study the transformation properties of the Weyl's group generators U and V under the U(1) mapping Since U = ω p−1 2 U 2π/p , it commutes with U φ . The transformation of V requires some thinking: Let us rewrite V , Eq. (59), as the sum of two operators i.e., the matrix that has only a 1 in the lower-left corner. Then, We are interested in analyzing how the transfer matrices T 1 and T 2 of Eq.
It is indeed easier to analyze the Fourier transform transfer matrices,T µ =F † T µF , Clearly,Û φ commutes withT 2 but does not withT 1 unless φ = 2π/p, which not surprisingly corresponds to the (Fourier transform) discreteQ symmetry. However,Û φ is an exact continuous symmetry of the modified transfer matrix and becomes the usual U(1) symmetry of the XY model when p → ∞. This emergent symmetry may allow for the construction of spin-wave excitations in the critical region. Note that in the original transfer matrix representation T 1 , T 2 , the continuous emergent symmetry is represented byFÛ φF † . Moreover, it is an emergent symmetry of both p-clock and self-dual p-clock models.

Phase diagram: From the p-clock to the XY model
We are thus left with the task of establishing the phase diagram of the p-clock model, the nature of its phase transitions and excitations, and its behavior as p → ∞. One may argue that the phase structure of the model is well understood [10] (see Fig. 4). At very low temperatures, there is a ferromagnetic phase characterized by long-range order of the two-point, spinspin, correlation function G(|r − r ′ |) = cos(θ r − θ r ′ ) , and the breakdown of the Z p symmetryQ. At very high temperatures, the system is in a disordered phase with G(|r − r ′ |) decaying as an exponential function of the distance. For 2 ≤ p ≤ 4, these two phases are separated by a continuous second-order phase transition of the Ising (p = 2, 4) or Potts (p = 3) type. (It is very easy to prove that the p = 4 case is identical to two uncoupled p = 2 Ising models [32], see Appendix D.) For p 5 there is an additional intermediate critical phase separating the ferromagnetic from the disordered phase. It is characterized by a power-law behavior of G(|r − r ′ |) ∼ |r − r ′ | −η with a nonuniversal exponent η, and by the absence of symmetry breakdown and quasilong-range order. In the p → ∞ limit the broken-symmetry phase disappears as one recovers the D = 2 XY model with a continous U(1) symmetry. This qualitative picture leaves several issues unresolved that numerical simulations have not been able to resolve either: • What is the nature of the two phase transitions for p 5?
• What is the nature of the relevant topological excitations in each phase?
• What is the physical origin of the critical (massless) phase?
• What is the minimum p after which the transitions are of the BKT type? Figure 4: Phase diagram of the p-clock model. For p ≥ 5 there are three phases, the broken Z p (low-temperature) phase disappearing in the limit p → ∞ (XY limit). A transition is of BKT-type whenever it is associated to an essential singularity of the free energy. The critical phase is characterized by power-law correlations, i.e., quasi-long-range order, with non-universal exponents.
To understand qualitatively the nature of the phases of the model, consider the ground state of the self-dual quantum Hamiltonian H p defined in Eq. (71), in the large (low-temperature) and small (high-temperature) K 1 limits (a 0 (K 2 ) = 0). Let us start with the broken Z p symmetry, low temperature sector that corresponds to the line (K 1 , a 1 (K 2 ) = 0). Then, the p-fold degenerate subspace of ground states is trivial to describe in terms of the simultaneous eigenvectors of the U i of Eq. (57), The ground state energy is E 0 = −K 1 N for periodic boundary conditions (E 0 = −K 1 (N − 1) for open boundary conditions), and Ψ r 0 |Ψ s 0 = δ rs [33].
The fully disordered, high-temperature phase is defined by the sector (K 1 = 0, a 1 (K 2 )). The ground state (E 0 = −2a 1 (K 2 )N) is unique and given by (in terms of the eigenstates of V i , Eq. (63)), and satisfies C 0 |Φ 0 = +|Φ 0 . It is difficult to obtain exact results for arbitrary couplings. There is, however, an interesting exact relation that holds at the self-dual line K 1 = 2a 1 (K 2 ) ≡ K * and follows from the fact that the self-duality unitary U d becomes a new symmetry of the problem on that line. It is clear from Eq. (71) that [13]. Since [H p [K * ], U d ] = 0, we can choose the energy eigenstates |Ψ n , n = 0, 1, · · ·, to be also eigenstates of U d , U d |Ψ n = e iφn |Ψ n . Then For 2 ≤ p ≤ 4, the p-clock model is exactly self-dual and the transition from the ferromagnetic to the disordered phase happens at the self-dual point K 1 = 2a 1 (K 2 ). For p ≥ 5, Eq. (69) shows that the p-clock model is no longer exactly self-dual, but the self-dual approximation of Eq. (71) or (78) allow us to establish the following self-dual equation for arbitrary p From the self-dual condition K 1 = 2a 1 (K 2 ) one can determine the self-dual temperature T * . The self-dual point is a point of non-analyticity of the free energy for 2 ≤ p ≤ 4, but for p ≥ 5 it is analytic. Some values are indicated in Table 1, assuming isotropic couplings K 1 = K 2 = J/(k B T ). It follows from very general considerations (see Section 8 of Ref. [13]) that the two critical points c 1 and c 2 bounding the self-dual point when p ≥ 5 are exactly related by It is interesting to analyze the large-p limit of the self-dual Eq. (102). In that limit BKT , goes to zero when p → ∞, as T BKT ∼ 1/p 2 , and the highest critical temperature T · · · 1/(2 ln(2 cos( π 9 ))) large p · · · 2π/p and from the asymptotic expansion of the modified Bessel functions I 0,1 , Eq. (27), one gets the relation between the transfer matrix and direct couplings 1 One can then use Eq. (103) to obtain a relation between the two critical temperatures T (1) BKT and T (2) The Peierls argument developed in Appendix E, on the other hand, provides a rigorous scaling for the lowest transition temperature, T , as p → ∞. Thus this lowest critical temperature vanishes for the classical XY model (where a broken-symmetry phase is not allowed) and, according to Eq. (106), T BKT ∼ O(1) for large p. The self-dual temperature T * has its own "intermediate" scaling with p, k B T * = (2π/p)J, so that it also vanishes in the p → ∞ limit. This is to be expected since the XY model is not self-dual, nor has a natural self-dual approximation.
To understand what makes p ≥ 5 different from p < 5 and explain the ap- two is that domain walls exist for any p, while vortex-like excitations do not exist for 2 ≤ p ≤ 4, becoming manifest only for p ≥ 5. Also, if 2 ≤ p ≤ 4 the energy cost to create a domain wall is independent of |∆s| = |s in − s out |, with s in(out) indicating the angular configurations at the two sides of the wall. This changes for p ≥ 5, allowing for twists of the spin of size |∆s| = 2, · · · , (p − 2).
In the Peierls argument provided in Appendix E, the upper bound for the probability of having a domain wall corresponds to a change of orientation between two domains of (±2π/p) (or equivalently |∆s| ∼ O(1)). Such a change always appears in all domain walls in systems with p = 2, 3, 4. This, however, is not the case for p ≥ 5 where there exist general domain wall topologies that do not allow for a uniform twist of the angle between neighboring domains. In such instances, |∆s| can be O(p) and, as shown in Fig. 5, two types of excitations are generally possible. More precisely, in p ≥ 5 systems, a vorticity arises. The topological invariant characterizing these configurations, that we call discrete winding number k, is given by the circulation sum taken around an oriented loop Γ, with the argument ∆θ rr ′ ∈ [−π, π) given by In contrast, the sum in Eq. (107) should not be taken mod(2π) but just as an ordinary sum of real numbers (otherwise it would vanish identically). To make this definition lucid, consider a (p = 5) configuration such as the one shown in Fig. 5, with θ = 2πs/p, and a loop Γ. Herein, the circulation sum explicitly reads (∆θ 01 +∆θ 12 +∆θ 23 +∆θ 30 ) = 2π/5+2π/5+2π/5+4π/5 = 2π. Thus, the configuration in Fig. 5 has a vortex of strength k = 1 at its origin. As θ 0 − θ 3 = −6π/5, we set ∆θ 30 = 4π/5. This shift in the value of ∆θ rr ′ (that must lie in the interval [−π, π)) leads to the non-zero value of k in this case.
We may now use energy-versus-entropy balance considerations to argue for the relevance of these topological excitations in establishing the two phase transitions. The Peierls argument presented in Appendix E rigorously establishes that domain walls oriented relative to one another by the minimal energy cost (i.e., twists of (±2π/p)) are responsible for the existence of a low-temperature ferromagnetic broken Z p symmetry phase. Both the energy penalties and entropic costs associated with such minimal cost domain walls scale with ℓ (the domain wall length), and the analysis leads to a transition temperature that behaves as k B T  Note that this energy-versus-entropy balance argument does not rely on the existence or non-existence of the self-dual property of the model.
It is important to mention that while the physics of the low-temperature phase is associated to the exact discrete Z p symmetry, the existence of vortexlike excitations is directly related to the emergence of the continuous sym-metryÛ φ unveiled in Section 3.5. This U(1) symmetry becomes more exact at high temperatures (T T * or 2a 1 K 1 ), for a fixed p ≥ 5, or it is exact at any temperature when p → ∞. Thus, the physical origin of the critical phase and the extended universality concept introduced in Ref. [19] is simply our emergentÛ φ continuous symmetry.
What is the nature of the phase transitions when p ≥ 5? Given the current debates [22], it is important to say what we mean by a "BKT-type phase transition". We simply mean a transition characterized by an essential singularity in the free energy (or the ground state energy in a quantum model). This includes those cases where there is an essential singularity but, for instance, the correlation function exponent η = 1/4 (1/4 is the exponent for the XY model [11]). It is very difficult to prove analytically the existence of an essential singularity but numerical simulations seem to indicate that for p ≥ 5 the two phase transitions are continuous with continuous derivatives [22], supporting the BKT scenario. Moreover, our new self-duality argument proves that the two transitions must be of the same nature [13]. In other words, if there is an essential singularity (the function and all its derivatives remain continuous) in the free energy at T (1) BKT , then, there should be the same type of singularity at T (2) BKT [13]. This, of course, does not mean that the self-duality fixes the value of, for instance, η, to be the same at the two transition points.
In what follows, we will show that when p ≥ 5, the temperature region T (1) BKT must be critical (massless) with algebraic correlations. The assumptions are: (1) The phase transitions at T (1) BKT and T (2) BKT are continuous, and (2) for large separations |r−r ′ |, the two-point correlation function G(|r− r ′ |, T ) = cos(θ r − θ r ′ ) is of the canonical Ornstein-Zernike-like form with M representing the order parameter, i.e., magnetization, ξ the correlation length, η an anomalous exponent, and A an amplitude. In principle, all of these quantities (M, ξ, η, and A) are functions of temperature T . We start by demonstrating that G(|r − r ′ |, T ) is a monotonically decreasing function of T for any p. To this end, we derive in Appendix F a Griffiths'-type inequality as in general ferromagnetic systems [34]. Assume, for simplicity, the uniform case K µ = K > 0. Then, it is straightforward to show (see Appendix F) that or, equivalently, ∂ T G(|r −r ′ |, T ) ≤ 0, which proves that G is a monotonically decreasing function of temperature (0 ≤ G(|r − r ′ |, T ) ≤ 1). Now, if G(|r − r ′ |, T ) is monotonically decreasing with temperature T for any fixed separation |r − r ′ | then, in the absence of magnetization (i.e., when M = 0) the correlation length ξ must be a monotonically decreasing function of temperature. The proof of this assertion is trivial. Consider two rather general temperatures in this region T a , T b , such that T BKT , then the ratio of the corresponding asymptotic correlation functions is because of the monotonicity property of G. This is only possible if ξ b < ξ a for otherwise for large |r − r ′ |, the ratio of the two correlation functions would diverge exponentially in |r − r ′ |. As T a and T b were rather general temperatures, it follows that ξ(T ) is a monotonically decreasing function of the temperature T so long as the magnetization M = 0 (as it is for T > T BKT ). Therefore, a divergence of the correlation length at T (1) BKT and T (2) BKT implies that within the entire interval T (1) BKT , the correlator G is algebraic in |r − r ′ |. We thus proved that there must exist a power law, critical phase, between T (1) BKT and T (2) BKT .

Acknowledgements
This work is partially supported by the National Science Foundation under Grant No. 1066293 and the hospitality of the Aspen Center for Physics. Fruitful discussions with C. D. Batista are acknowledged.

Appendix A. Exponential of shift operators
In this appendix we collect some useful formulas to compute exponentials and logarithms of shift operators. Let us start with the operator L z , the infinitesimal generator of translations on the circle. We have the general relation This, together with the orthogonality relation 2πδ(θ ′ − θ) = n∈Z e i(θ ′ −θ)n , leads to In Section 3.1 we introduced a diagonal matrix U, and a shift operator V † that describes translations in a p-points discretization of the circle. This shift operator plays a role similar to that of L z (actually, e −iθLz ). Now we have the general relation These are the expressions that are most useful in physical applications. Suppose next that b m = e K u(m) , where K is a positive constant, and u(m) is a real function of m = 0, · · · , p − 1 (for example, u(m) = cos( 2πm p ) for the classical p-clock model). We would like to study the behaviour of the a m to next-to-leading order in K, in the limit that K grows very large (this could happen at low temperature). Notice that in this limit to next-to-leading order, assuming that the inequalities In this section we study a duality that illustrates the essential differences between compact,θ, and non-compact,x, degrees of freedom. The algebraic tool of choice is the q-oscillator algebra [35], specified by a positive real number q, a creation operator a † , its Hermitian conjugate a, and a Hermitian operatorn, satisfying To compute Z D XY , one should take the trace in the eigenbasis ofn i . This basis is described in Ref. [35].
This description of the XY model in terms of q-deformed bosons with q = e −2 suggests that the algebra of Eq. (B.1) affords a continuous interpolation between the XY model and ordinary phonons (characterized by q = 1), but this is not the case: The XY model belongs to a representation of the algebra of Eq. (B.1) that is inequivalent to that describing phonons (i.e., canonical bosons). The reason is that Eq. (7) is not enough to specify the algebra of translations in the circle. We must also have that e ±iθ e ∓iθ = 1. (B.4) The mapping of Eq. (B.2) will respect this constraint only if a, a † satisfy at least for q = e −2 , including the relations listed in Eq. (B.1). But the resulting set of four relations becomes inconsistent at q = 1. This shows that the q-oscillator algebra cannot interpolate continuously between canonical bosons and compact excitations.
was introduced in Ref. [16] to provide a Gaussian approximation to the XY model that preserves the essential property of compacticity, and is a good approximation at sufficiently low temperatures. We now show, by using our bond-algebraic approach, that it is dual to the solid-on-solid (SoS) model of the roughening transition, characterized by integer-valued degrees of freedom m r ∈ Z [11]. We work directly in the thermodynamic limit, N → ∞, to avoid dealing with boundary terms.
The transfer operator T SoS = T 2 T 1 for the SoS model can we written as in terms of the operators X i , R i , R † i defined in Eq. (32). Now, however, i ∈ Z labels the sites of an infinite straight line. The duality of bond algebras affords a dual representation of T SoS , in terms of compact degrees of freedom. The next step is to compute Z D SoS = tr [(T D 2 T D 1 ) N ] in the basis introduced in Eq. (9). T D 2 is already diagonal in that basis At this point we could proceed by analogy to previous sections and rewrite this expression in terms of an interaction potential V K (θ) ≡ − ln m e K 2 2 m 2 e −iθm , but this will not turn out be the most convenient approach. Instead, let us proceed to compute the matrix elements of T D 1 . This task reduces to computing the matrix elements of a one-body operator, which results from recalling that the orthonormal states of L z,i are the plane waves θ i |n i = e iθ i n i / √ 2π. Notice that the function e K 2 x 2 is the Fourier transform of e x 2 2K / √ K. It then follows that we can use Poisson's summation formula to write Putting all the pieces together, we obtain The last expression is exactly (2π/K 2 ) N Z V [K µ ], and thus the Villain model is dual to the SoS model. Notice the reciprocal relation between the couplings: The Villain model is strongly coupled only if its dual SoS representation is weakly coupled.
Appendix D. The p-clock model for p = 2, 3, and 4 Let us start with the simplest p = 2 case. Then, U i = U † i = σ z i and V i = V † i = σ x i , and the transfer matrix T p = T 2 T 1 of Eq. (64) reduces to This finite Ising model is self-dual up to boundary corrections. The substitution T 1 → e K 1 σ z N T 1 renders the model exactly self-dual for any N [13]. If p = 3, then V †2 i = V i , and T 2 becomes Moreover, it is easy to check that the operator 2C i = (1l + σ z (known as a controlled-NOT gate in quantum computation) maps and thus it follows that C = N i=1 C i maps that clearly defines two decoupled Ising models, with couplings that are half of those of the p = 2 model. In particular, the p = 4 clock model is exactly self-dual provided T 1 → e K 1 2 (σ z 1,N +σ z 2,N ) T 1 .

Appendix E. Peierls argument for the p-plock model
We now use the Peierls argument to prove that there should be a broken symmetry phase (low-temperature ordered phase) in the p-clock model on the square lattice. The proof establishes the existence of a phase transition at a temperature T (1) below which global Z p symmetry is broken. For large p ≫ 1, Specifically, our objective is to show that if uniform boundary conditions pertaining to one of the clock states θ = 2πs/p, with 0 ≤ s ≤ p − 1 a fixed integer, are applied on the boundary of the square lattice, then there provably exists a temperature T Peierls > 0 such that for temperatures T < T Peierls , spontaneous symmetry breaking (SSB) of the global Z p symmetry arises. (In the context of our discussions thus far, T Peierls < T (1) ; asymptotically, for large p, both temperatures scale as 1/p 2 .) By SSB in this context, we refer to the lifting of the symmetry triggered by applying the uniform boundary conditions at spatial infinity. That is, when the aforementioned boundary conditions are introduced then, for T < T Peierls , the probability distribution P(θ 0 ) for the angular orientation of the spin at the origin S 0 is not symmetric between the p possible values of θ 0 . In particular, we will demonstrate that P(θ 0 ) is maximal when θ 0 has an orientation that matches that on the boundary, θ ∞ . In other words, for temperatures T < T Peierls , To prove this inequality, we note that where Prob(X) denotes the probability of the set of events X. The domain wall is defined as the boundary between differently oriented spins. The logic underlying Eq. (E.2) is clear: if θ 0 = θ ∞ then, by its very definition, at least one domain wall must separate the spin at the origin from the spins on the boundaries of the lattice. We now will bound the probability of having a particular domain wall. Specifically, let us denote by {C α } the set of configurations that have Γ as the outer-most domain wall surrounding the origin. That is, Γ separates spins with an orientation θ = θ ∞ from those having another (uniform) orientation θ in . [Note that, generally, more than one domain wall may be present and thus θ in need not be the same as θ 0 .] The upper bound on the probabilities in Eq. (E.2) is a sum over the probabilities of having such different outer-most domain walls Γ. We furthermore define the partition function where E α ≡ E(C α ) is the energy of the spin configuration C α . Thus, Z Γ is smaller than the total partition function Z p of the system. This is so as Z Γ contains only a subset of the Boltzmann weights appearing in Z p . That is, in Eq. (E.3) we sum only over spin configurations with at least one domain wall surrounding the origin. We next define a new configurationC α formed by rotating all of the spins inside Γ by a uniform angle ∆θ such that the outermost domain wall that surrounds the origin is removed. That is, for all spins S r that (i) lie inside the region bounded by the domain wall Γ, we perform the transformation θ r → (θ r + ∆θ) with an angle of rotation where ∆s is an integer. (ii) All spins lying outside the domain wall Γ have an orientation θ = θ ∞ ; these spins are not rotated. When present, any other more internal domain walls will remain unchanged by this uniform rotation of all the spins inside Γ. In order to bound, from above, the probability of having an outermost domain wall Γ, we now consider (Eᾱ = E(C α )) The probability of having the domain wall Γ is fixed by the ratio of the sum of Boltzmann weights associated with having the domain wall Γ divided by the sum of Boltzmann weights associated with all spin configurations (i.e., the partition function Z p ). As ZΓ contains a sum only over a subset of all Boltzmann weights that appear in Z p , we have The smallest energy difference between a configuration C α andC α is bounded by where ℓ is the length of the domain wall Γ (exchange constants are set to unity, J µ = 1, µ = 1, 2  with x = e β(1−cos 2π p ) /3. In performing the summation in Eq. (E.12), we assumed a sufficiently low temperature so that x > 1, andw is a monotonically decreasing function of β. Notably,w can be made arbitrarily close to zero for large enough β. Let us denote by β Peierls the solution to the equationw(β Peierls , p) = p−1 p . Then, for β > β Peierls , the probability of the spin at the origin being the same as that on the boundary is P(θ 0 = θ ∞ ) > 1/p. In other words, for T < T Peierls , we clearly have SSB. It is important to emphasize that T Peierls is only a lower bound to the transition temperature, and the actual SSB occurs for T (1) > T Peierls . An estimate for T Peierls resulting from this analysis is T Peierls ≈ (1 − cos 2π p )/ ln 6 (∼ O(1/p 2 ) for large p). As in the lower bound derived herein, the energy cost for a domain wall is anticipated to determine the actual ordering temperature. This bound is rigorous. Physically, in Eq. (E.10), the logarithms of the two terms N ℓ and D ℓ capture, respectively, bounds on the entropy and energy costs associated with domain walls of length ℓ.
Discrete vortices such as the one shown in Fig. 5 with a typical change of angle ∆θ = O(1) (or |∆s| = O(p)) across the intersecting domain walls that extend over a linear distance ℓ may entail, for all p ≥ 5, an energy cost that scales as ℓ. This is to be contrasted with the minimal energy penalty associated with a difference in angle of |∆θ| = 2π/p for which the corresponding energy penalty as ℓ/p 2 (and that physically sets the bounds that we derived in the Peierls argument above). Thus, from energy-versus-entropy balance considerations, the temperature below which it is unfavorable to have vortices is T Vortex ∼ O(1) (or of order J): the energy for such domain walls scales as ℓ as does the entropy associated with a network of possible intersecting domain walls that have a total length ℓ.

Appendix F. Proof of monotonicity of the correlation function G
We now prove Eq. (110) for large (yet finite) lattices. In finite size systems (no matter how large), there are no thermodynamic phase transitions. Thus, the free energy and all its derivatives (including, in particular, the two-point correlation function G(|r − r ′ |, T )) are analytic for all values of β. We will prove monotonicity by expanding (the analytic) G as a power series in β, which is everywhere convergent, and illustrate that the coefficients multiplying each power of β are non-negative. The uniformity of the sign of all contributions to the series coefficients follows from repeated applications of the identity with the Kronecker delta above defined mod(p). That is, δ 0,0 = δ ±p,0 = δ ±2p,0 = · · · = 1, otherwise it vanishes.
Longhand, the correlator G(|r − r ′ |, T ) is given by (K = βJ ≥ 0) with θ r = 2πs r /p. We Taylor expand the argument of the exponential, i.e, exp[βA] = 1 + βA + (βA) 2 /2! + · · · , in both the numerator and the denominator of Eq. (F.2), and may represent pictorially the expansion terms by Feynman-type diagrams. In this scheme, each appearance of cos(θ x ′ − θ x ) relates to an internal propagator linking sites x and x ′ . Similar to perturbative schemes in the continuum, where the number of Feynman-type "bubble diagrams" in Z p [K] with nearest-neighbor links that share common vertices with any given "connected" diagram which contains both the sites r and r ′ is negligible, there is a near cancellation, in Eq. (F.2), of all "bubble diagrams". What remains from the ratio of Eq. (F.2) is the sum of all diagrams in which sites r and r ′ are linked to each other by the line stemming from cos(θ r − θ r ′ ) as well as internal lines which pass through the points x = 1, 2, · · · , m. Up to (inherently positive) symmetry factors, the numerical value of any such resulting diagram is given by a sum of the form p−1 sr,s 1 ,s 2 ,··· ,sm,s r ′ =0 β t r1 +t r2 +···+t 12 +···+t mr ′ cos(θ r − θ r ′ )(cos(θ r − θ 1 )) t r1 ×(cos(θ r − θ 2 )) t r2 · · · (cos(θ 1 − θ 2 )) t 12 · · · (cos(θ m − θ r ′ )) t mr ′ , (F. 3) where the integers t ab > 0 represent the number of lines linking sites a and b. It is simple to show that sums of the form of Eq. (F.3) are manifestly non-negative. Replacing cos(θ a − θ b ) by (exp(i(θ a − θ b )) + exp(i(θ b − θ a )))/2, a sum of exponentials with positive weights, Eq. (F.3) reduces to a sum of individual products with positive weights, each being of the form x=r,1,2,··· ,m,r ′ S p (n x ) ≥ 0, (F. 4) with n a set by sums of the powers t ab in Eq. (F.3). Then, it follows that when the ratio of Eq. (F.2) is expanded in β, i.e., G(|r − r ′ |, T ) = t a t β t , the prefactor a t multiplying each individual power β t is non-negative. As such, G(|r − r ′ |, T ) is manifestly monotonic in β, i.e., ∂ β G(|r − r ′ |, T ) ≥ 0.
which shows that the p = 4 clock model is equivalent to two decoupled Ising systems.