<figure><img src="image_1.png" /></figure>

HAL
open science

# Phase transitions in q-state clock model

Arpita Goswami, Ravi Kumar, Monikana Gope, Shaon Sahoo

► To cite this version:

Arpita Goswami, Ravi Kumar, Monikana Gope, Shaon Sahoo. Phase transitions in q-state clock model. 2024. hal-04653617v2

HAL Id: hal-04653617

https://hal.science/hal-04653617v2

Preprint submitted on 2 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Phase transitions in $q$-state clock model

Arpita Goswami, Ravi Kumar, Monikana Gope, and Shaon Sahoo*

Indian Institute of Technology Tirupati, India, 517619

The $q$-state clock model, sometimes called the discrete XY model, is known to show a second-order (symmetry breaking) phase transition in two-dimension (2D) for $q \le 4$ ($q = 2$ corresponds to the Ising model). On the other hand, the $q \to \infty$ limit of the model corresponds to the XY model, which shows the infinite order (non-symmetry breaking) Berezinskii-Kosterlitz-Thouless (BKT) phase transition in 2D. Interestingly, the 2D clock model with $q \ge 5$ is predicted to show three different phases and two associated phase transitions. There are varying opinions about the actual characters of phases and the associated transitions. In this work, we develop the basic and higher-order mean-field (MF) theories to study the $q$-state clock model systematically. Our MF calculations reaffirm that, for large $q$, there are three phases: (broken) $\mathbb{Z}_q$ symmetric ferromagnetic phase at the low temperature, emergent $U(1)$ symmetric BKT phase at the intermediate temperature, and paramagnetic (disordered) phase at the high temperature. The phase transition at the higher temperature is found to be of the BKT type, and the other transition at the lower temperature is argued to be a large-order spontaneous symmetry-breaking (SSB) type (the largeness of transition order yields the possibility of having some of the numerical characteristics of a BKT transition). The higher-order MF theory developed here better characterizes phases by estimating the spin-spin correlation between two neighbors.

# I. INTRODUCTION

The two-dimensional (2D) Ising model shows a spontaneous symmetry-breaking (SSB) phase transition at a finite temperature [1, 2]. The 2D XY model does not show this SSB transition due to the restriction imposed by the Mermin-Wagner theorem [3]. This model, however, shows a different type of phase transition, known as Berezinskii-Kosterlitz-Thouless or BKT transition [4–8], due to the unbinding of bound vortexes-antivortexes pairs (topological defects). The $q$-state clock model, a generalization of the Potts model [9, 10] and one of the different possible discretizations of the XY model [11, 12], represents the Ising and the XY models in two different limits. The clock model becomes the $\mathbb{Z}_2$ symmetric Ising model when $q = 2$ and it represents the $U(1)$ symmetric XY model in the $q \to \infty$ limit [13–15]. The clock model exhibits rich physics in terms of the nature of its phases and the associated critical phenomena.

Many questions arise while studying the 2D clock model. For example, since the Ising and the XY model show two different types of phase transitions, it is of interest to know what happens to the 2D clock model when $q$ is finite and $q > 2$. It is generally agreed that there is some critical value $q_c$; the 2D clock model shows only one phase transition (of the SSB type) when $q < q_c$. On the other hand, when $q \ge q_c$, the model is predicted to show two phase transitions [13–16]. It is now generally accepted that $q_c = 5$ [17, 18].

The next issue is about the nature of phases and the associated phase transitions. For a finite value of $q$, the Mermin-Wagner theorem [3] does not apply to the clock model.

In addition, the non-symmetry breaking BKT phase transition is originally shown to exist for the XY model where spins have continuous rotational symmetry $U(1)$. This raises the question of whether the phase transitions observed in the 2D clock model are of SSB type, especially for small values of $q$ (say, 5 or 6). It is then somewhat counterintuitive that a large number of the recent studies point out that the natures of the phase transitions in the $\mathbb{Z}_q$ symmetric clock model are of the BKT type [15, 19–22]. The counterclaims about the nature of the transitions do also exist, especially for the small values of $q$ [23–25].

Although there are varying opinions on the nature of the phase transitions, it is generally accepted that there are three phases for $q \ge 5$ [15, 20, 21, 26], which include an ordered ferromagnetic phase at low temperature and a disordered paramagnetic phase at high temperature. There is an intermediate phase which is predicted to be a quasiliquid (sometimes called the BKT phase). It is important to investigate how the $U(1)$ symmetry arises from a $\mathbb{Z}_q$ symmetric spin model and how the BKT phase forms (in terms of the vortexes-anti vortexes pairs) with this emergent symmetry. Since many of the issues in this field are not yet completely resolved and some results require more clarity, it is worthwhile to re-investigate and gain more insight into the nature of the phases and associated transitions. To this effort, we develop a mean-field (MF) theory (both basic and higher order) for the $q$-state clock model. To the best of our knowledge, ours is the first systematic mean-field study of the $q$-state clock model.

The MF theory is an important approach to analyze properties of model systems [2]. Its success in analyzing classical or quantum models is well known. The MF theories work better in the higher dimensions and/or when the coordination number is large. Although, in this

<footnote>* Corresponding author: shaon@iittp.ac.in</footnote>
work, we are mainly interested in the two-dimensional systems, we expect to get reasonable results in the large $q$ limit. This is because the change in the mean field will be much less due to elementary excitation or low-temperature fluctuation in the neighbors. Additionally, since the MF theory works better in higher dimensions, our MF theory can help us analyze properties of higher dimensional (say, 3-dimensional) $q$-state clock model.

Our main contributions in this paper include the following. We develop the basic and higher order MF theories for the $q$-state clock model and our analysis reaffirms the existence of three separate phases in 2D clock model (for large $q$). We demonstrate how the $U(1)$ symmetry arises in the intermediate temperature regime for a Hamiltonian with $\mathbb{Z}_q$ symmetry. We find that the transition at the higher temperature is of BKT type, and the transition temperature is almost $q$-independent, i.e., $T_{BKT} = O(1)$. We argue that the transition at the lower temperature is a large-order SSB transition with possibility of having some numerical characteristics of a BKT transition. This transition temperature decreases as $1/q^2$, i.e., $T_{SSB} = O(q^{-2})$.

This article is arranged in the following way. In Sec. II, we develop the basic (zeroth order) MF theory and analyze the 2D clock model using the theory. Next, in Sec. III, we construct the first-order MF theory which gives us better insight into the phases of the model. Using this higher-order theory, we estimate the nearest-neighbor spin-spin correlation and find a slightly better estimation of the BKT transition temperature. We conclude our work in Sec. IV.

## II. MEAN-FIELD THEORY

The $q$-state clock model represents a spin system where the spins lie in the $XY$ plane and can orient only in $q$ possible directions. This model, without a magnetic field, is represented by the following Hamiltonian:

$$H = -J \sum_{\langle ij \rangle} \cos(\theta_i - \theta_j), \qquad (1)$$

where $J$ is the nearest-neighbor exchange constant (taken to be positive) and $\langle ij \rangle$ represents summation over the nearest neighbors. Here $\theta_i$ is the angle between the $i$-th spin and the reference, say $X$, axis. This angle takes $q$ possible values $\theta_i = 2\pi k_i/q$ where $k_i$ can take any integer value from 1 to $q$. As the exact results for this model are scarce, it is important to develop an appropriate mean-field approach to understand the behavior of this model.

As a first step in developing the (basic or zeroth order) mean-field theory for this model, we define a complex spin variable $S_k$ corresponding to the $k$-th spin: $S_k = e^{i\theta_k}$, where $\theta_k = 2\pi j_k/q$. Here $j_k$ takes any integer value from 1 to $q$. With this complex spin variable, we can rewrite the Hamiltonian in Eq. 1 as,

$$H = -JRe \sum S_i S_j^*. \qquad (2)$$

We now perform the basic (zero-th order) MF theory over the complex spin variables. In the MF approach, each variable ($S_i$) is replaced by its mean value plus a small fluctuation term ($\delta S_i$),

$$S_i = \langle S_i \rangle + \delta S_i.$$

If we define the complex mean-field parameter as $l = \sum_{i=1}^N \frac{\langle S_i \rangle}{N}$, for translationally invariant system, we can write

$$S_i S_j^* = (l + \delta S_i)(l + \delta S_j)^* \approx |l|^2 + l \delta S_j^* + l^* \delta S_i,$$

after neglecting the small (second order) term $\delta S_i \delta S_j^*$. After replacing $\delta S_i$ by $S_i - l$, we get the following mean-field Hamiltonian from Eq. 2,

$$H_{MF} = -JRe \sum_{\langle ij \rangle} (|l|^2 + l^*(S_i - l) + l(S_j - l)^*).$$

We here note that $\sum_{\langle ij \rangle} S_i = \sum_{\langle ij \rangle} S_j$ and $\sum_{\langle ij \rangle} = \frac{1}{2} \sum_{i=1}^N \sum_{j \in (nn)}$, where $\sum_{j \in (nn)}$ represents the sum over the nearest neighbors, $N$ is the total number of lattice sites. With these, we get

$$H_{MF} = -J \frac{1}{2} Re(-N r |l|^2 + r \sum_{i=1}^{N} (l^* S_i + l S_i^*)), \quad (3)$$

where $r$ is the number of nearest neighbors. As $l^* S_i + l S_i^* = 2Re(l^* S_i)$, we can write

$$H_{MF} = \frac{1}{2} N J r |l|^2 - J r \sum_{i=1}^{N} Re(l^* S_i). \qquad (4)$$

The mean-field parameter $l$ is a complex quantity, which we now write in terms of two real parameters: $l = l_1 + il_2$. Physically, $l_1$ and $l_2$ represent the average values of the $X$ and the $Y$ components of a spin. With $S_j = \cos \frac{2\pi j_k}{q} + i \sin \frac{2\pi j_k}{q}$, we can recast Eq. 4 as

$$H_{MF} = \frac{1}{2} N J r (l_1^2 + l_2^2) - J r \sum_{i=1}^{N} \left( l_1 \cos \frac{2\pi i_k}{q} + l_2 \sin \frac{2\pi i_k}{q} \right). \qquad (5)$$

It may be mentioned again here that $i_k$ can take any of the $q$ values from 1 to $q$. If we write $H_{MF} = \sum_{i}^N h_i$, we get the following effective mean-field Hamiltonian corresponding to the $i$-th site,

$$h_i = \frac{1}{2} J r (l_1^2 + l_2^2) - J r \left( l_1 \cos \frac{2\pi i_k}{q} + l_2 \sin \frac{2\pi i_k}{q} \right). \qquad (6)$$

The $q$ possible values of $h_i$ are given by

$$\epsilon(k) = \frac{1}{2} J r (l_1^2 + l_2^2) - r (l_1 \cos \frac{2\pi k}{q} + l_2 \sin \frac{2\pi k}{q}), \quad (7)$$

where $k$ takes any integer value between 1 and $q$. The energy $\epsilon(k)$ does not depend on the site index $i$ for a transitionally invariant system.
In the following, we will calculate the free energy to determine $l_1$ and $l_2$ and to analyze the equilibrium properties from the free energy minima at a given temperature.

The partition function of each spin is given by,

$$z = \sum_{k=1}^{q} e^{-\beta\epsilon(k)} \qquad (8)$$

Since the spins in the mean-field approach are non-interacting, the total partition function of the system is given by $Z = z^N$. The free energy corresponding to a single spin is

$$F = -\frac{\ln z}{\beta}, \qquad (9)$$

where $\beta = \frac{1}{k_B T}$ with $k_B$ being the Boltzmann constant and $T$ being the absolute temperature. Written explicitly, we get

$$F = \frac{Jr(l_1^2 + l_2^2)}{2} - \frac{1}{\beta} \ln \sum_{k=1}^{q} e^{\beta J r(l_1 \cos \frac{2\pi k}{q} + l_2 \sin \frac{2\pi k}{q})} \qquad (10)$$

### A. Free energy plots at different temperatures

The free energy $F$ in Eq. 10 is a function of $l_1$, $l_2$ and $\beta$. The system equilibrium for a given temperature corresponds to the free energy minimum (minima) with respect to $l_1$ and $l_2$. Thus, to first gain insight into the problem, we plot $F$ as a function of $l_1$ and $l_2$ for different fixed temperatures. The plots can be seen in Fig. 1.

#### 1. Three different phases

We first note that, at high temperatures, we have only one minimum at $l_1 = l_2 = 0$. On the other hand, at low temperatures, we have a local maximum at the origin ($l_1 = l_2 = 0$) but $q$ number of minima around the origin. Interestingly, there is an intermediate range of temperature where we get a ring of minima around the origin (in this case the free energy looks like a $U(1)$ symmetric “Mexican hat”).

As mentioned earlier, the parameters $l_1$ and $l_2$ respectively represent the average values of the X and Y components of a spin. Keeping this in mind, we can have the following interpretations of the above observations. At high temperatures, we have one free energy minimum at the origin ($|l|^2 = l_1^2 + l_2^2 = 0$). In this phase, an individual spin can orient along any possible direction. This corresponds to the disordered (paramagnetic) phase. As we lower the temperature, we see the minimum turns into a local maximum, and a ring of minima forms around the maxima. This is the emergence of the $U(1)$ symmetry from a $Z_q$ symmetric model. This is further discussed in Sec. II B. Unlike in the paramagnetic phase, in this

new phase, an individual spin can orient along a fixed direction ($|l|^2 = l_1^2 + l_2^2 \neq 0$), but different spins can orient in different directions. In fact, we later show using the higher order MF theory (Sec. IIIC) that the nearest-neighbor correlation takes an intermediate value between 0 and 1 (these two values, respectively, correspond to the disordered phase and ferromagnetic phase). We identify this phase with the bound vortex-antivortex or Berezinskii-Kosterlitz-Thouless (BKT) phase. It may be mentioned here that the individual spins will break the emergent $U(1)$ symmetry in this phase, as was conjectured earlier for higher dimension in Ref. [11].

If we lower the temperature further, we get $q$ number of minima around the origin ($l$ takes $q$ possible non-zero (complex) values corresponding to the minima). Although all spins can orient along any of the $q$ possible directions at the low temperature, the spins will spontaneously choose one of them. This phase is identified as the (broken) $Z_q$ symmetric ordered (ferromagnetic) phase.

#### 2. Nature of two phase transitions

The transition between the paramagnetic phase and the emergent $U(1)$ symmetric BKT phase is called the BKT phase transition, and the corresponding temperature would be here denoted by $T_{BKT}$. The other transition between the ferromagnetic and the BKT phases needs more careful analysis; we deliberate on this below.

Starting from the BKT phase, if we lower the temperature further, the ring of minima will disappear, and $q$ distinct localized minima will form (Fig. 1). Although all $q$ minima are equally probable, the system will spontaneously break the $Z_q$ symmetry, and the spins will pick up a particular direction corresponding to one of the minima. Since the symmetry of the model ($Z_q$) is spontaneously broken as we go across the transition point, we call this transition a spontaneously symmetry breaking (SSB) transition and the corresponding temperature is denoted by $T_{SSB}$.

It is worth mentioning here that many recent studies found this transition between the ferromagnetic and BKT phases to be of BKT type [15, 20, 21]. There are also some studies where the BKT nature of this transition is ruled out [23, 25]. As explained, at this transition, the symmetry is certainly spontaneously broken, hence it is a SSB transition but it may not be a second order or any small-order transition. We therefore conjecture that the transition at the lower temperature is a large-order SSB transition; because of the large transition order, it would exhibit some of the numerical characteristics of a BKT transition. Since no finite order derivative of the mean-field free energy (Eq. 10) show any singularity, it is not possible to directly determine the actual order of this transition using our MF approach. But since there are $q$ free energy minima below $T_{SSB}$, the Landau-style free energy functional will be a large-order polynomial
<figure><img src="image_3.png" /><figcaption>FIG. 1: Basic mean-field results: free energy plots at low, intermediate and high temperatures are shown in (a)-(c), (d)-(f) and (g)-(i), respectively. Results in the left, middle and right columns are for $q = 5, 6$ and $7$, respectively.</figcaption></figure>

in the MF parameters (the functional may contain some other non-algebraic terms). The order of such polynomial can be at-least $q + 2$. We know from the theory of higher-order (greater than second order) phase transition that one requires a large-degree polynomial in the energy functional to describe the transition [27, 28]. This leads to our aforementioned conjecture on the nature of the phase transition at the lower temperature in the 2D clock model for large $q$.

## B. Understanding emergence of $U(1)$ symmetry

As seen in Fig. 1, also as discussed in some earlier works (see, e.g., Refs. [11, 20, 29]), there is an emergence of continuous $U(1)$ symmetry in the $\mathbb{Z}_q$ symmetric 2D clock model in the intermediate temperature regime. This emergence can be understood in a straightforward way using mean-field calculations.

To appreciate how a continuous $U(1)$ symmetry emerges in a $\mathbb{Z}_q$ symmetric model, we first express the free energy, found in Eq. 10, as a function of parameters $|l|$ and $\theta$, where $l = l_1 + il_2 = |l|e^{i\theta}$. We find that

$$F(|l|, \theta) = \frac{Jr|l|^2}{2} - \frac{1}{\beta} \ln \sum_{k=1}^{q} e^{\beta J r |l| \cos\left(\frac{2\pi k}{q} - \theta\right)}. \quad (11)$$

In the high or in the intermediate temperature limit (where $\beta$ is small and $|l| < 1$), we can effectively expand $e^{\beta J r |l| \cos\left(\frac{2\pi k}{q} - \theta\right)}$ in a series. Noting that $\sum_{k=1}^q \cos^n\left(\frac{2\pi k}{q} - \theta\right) = 0$ or $\frac{q}{2^n}C(n, \frac{n}{2})$, respectively, for odd or even $n$ (here $C(n, m)$ denotes a binomial coefficient), we conclude that the free energy $F(|l|, \theta)$ does not depend on $\theta$ in the high or intermediate temperature limits. This result shows that there is an emergence of the continuous $U(1)$ symmetry in the aforementioned temperature limits. This symmetry disappears in the low-temperature limit, and the free energy $F(|l|, \theta)$ starts depending on the parameter $\theta$. As our numerical results show (in Fig. 1), the free energy will have $q$ number of minima, and the system will recover its $\mathbb{Z}_q$ symmetry (albeit the fact that there will be spontaneous breaking of this discrete symmetry in the low-temperature limit).

To see how the $U(1)$ symmetry disappears and the free energy starts depending on $\theta$, we note that in this limit, $\beta$ is very large (and $|l| \approx 1$). In this limit, it is not advisable to do calculations using a series expansion of an exponential function. To understand what happens in this limit, we first replace the sum in Eq. 11 by a more symmetric one in the following way:
$$F(|l|, \delta\theta) = \frac{Jr|l|^2}{2} - \frac{1}{\beta} \ln \sum_{k=-\frac{q}{2}+1}^{\frac{q}{2}} e^{\beta J r |l| \cos(\frac{2\pi k}{q} - \delta\theta)}. \quad (12)$$

Although we assume here $q$ to be even, the results followed are valid for the odd $q$ as well. Since we want to see how the free energy starts depending on the parameter $\theta$, we here assume a small value ($\delta\theta$) for the parameter for mathematical convenience. In the following we would assume $\delta\theta < \frac{2\pi}{q}$. Noting that, in the low temperature limit, $e^{\beta J r |l| \cos(\frac{2\pi k}{q} - \delta\theta)}$ is significant only for those values of $k$ for which $|\cos(\frac{2\pi k}{q} - \delta\theta)| \approx 1$, we approximate $F(|l|, \delta\theta)$ in the following way:

$$F(|l|, \delta\theta) \approx \frac{Jr|l|^2}{2} - \frac{1}{\beta} \ln Q(t, \delta\theta), \quad (13)$$

where $Q(t, \delta\theta) = e^{t \cos(\frac{2\pi}{q} + \delta\theta)} + e^{t \cos(\frac{2\pi}{q} - \delta\theta)} + e^{t \cos(\delta\theta)}$ with $t = \beta J r |l|$. Here we only consider the three most significant values for $k$, i.e., $k = 0, \pm 1$ in the sum in Eq. 12. This approximation is valid as long as we are in a low-temperature regime and $\delta\theta < \frac{2\pi}{q}$. Now if $q$ is large enough so that $\cos(\frac{2\pi}{q}) \approx 1 - \frac{1}{2}(\frac{2\pi}{q})^2$ and $\sin(\frac{2\pi}{q}) \approx \frac{2\pi}{q}$, then,

$$Q(t, \delta\theta) \approx e^{t \cos \delta\theta} \left[ 1 + 2e^{-\frac{1}{2}\left(\frac{2\pi}{q}\right)^2 t \cos \delta\theta} \cosh\left(\frac{2\pi}{q} t \sin \delta\theta\right) \right]. \quad (14)$$

Now to proceed further, we work in the low-temperature regime where $t(\frac{2\pi}{q})^2 < 1$. This assumption is clearly not valid when $\beta \to \infty$ (or $t \to \infty$), but can be appropriate near the phase transition point at the lower temperature

(i.e. near $T_{SSB}$). Noting that $\delta\theta < \frac{2\pi}{q}$, we get the following value of $Q(t, \delta\theta)$ retaining the terms up to $(\delta\theta)^2$:

$$Q(t, \delta\theta) \approx 3e^t \left[ 1 - \frac{1}{3}t\left(\frac{2\pi}{q}\right)^2 - \frac{1}{2}t(\delta\theta)^2 \right]. \quad (15)$$

The above expression for $Q(t, \delta\theta)$ tells us how the free energy in Eq. 13 starts depending on the parameter $\theta$ as we lower the temperature, starting from the phase with the emergent $U(1)$ symmetry. The temperature at which this transition happens can be estimated from Eq. 15. At the transition temperature (denoted here by $T_{SSB}$), the $\theta$-dependent term becomes significant in comparison to the $\theta$-independent terms in Eq. 15, i.e., $1 - \frac{2}{3}t(\frac{2\pi}{q})^2 = O(t(\delta\theta)^2)$ at the transition temperature. We also note that a physically significant change in $\theta$ is of the order of $\frac{2\pi}{q}$, i.e., we have $\delta\theta = O(\frac{1}{q})$ in Eq. 15. This gives us: $1 - \frac{2}{3}t(\frac{2\pi}{q})^2 = O(\frac{t}{q^2})$ at the transition temperature. Solving this we get that $t = O(q^2)$ or $T_{SSB} = O(\frac{1}{q^2})$.

### C. Calculation of $T_{BKT}$

In this part we estimate $T_{BKT}$, the temperature at which transition happens between the BKT and the paramagnetic phases. This transition corresponds to the free energy minimum at the origin ($l_1 = l_2 = 0$) turning into a local maximum. Thus, we will get an inflection of the free energy $F(l_1, l_2)$ at the temperature $T_{BKT}$.

In practice, we get the transition temperature from either of the equations $\frac{\partial^2 F}{\partial l_1^2}|_{(0,0)} = 0$ or $\frac{\partial^2 F}{\partial l_2^2}|_{(0,0)} = 0$. From Eq. 10, we get

$$\left(\frac{\partial^2 F}{\partial l_1^2}\right)_{(0,0)} = Jr - \frac{(Jr)^2 \beta \sum_{j=1}^q 1 \sum_{j=1}^q \cos^2 \frac{2\pi j}{q}}{(\sum_{j=1}^q 1)^2} + \frac{\beta(Jr)^2 (\sum_{j=1}^q \cos \frac{2\pi j}{q})^2}{(\sum_{j=1}^q 1)^2} = 0, \text{ i.e.,} \\ Jr - \frac{(Jr)^2 \beta \sum_{j=1}^q \cos^2 \frac{2\pi j}{q}}{q} + \frac{\beta(Jr)^2 (\sum_{j=1}^q \cos \frac{2\pi j}{q})^2}{q^2} = 0 \quad (16)$$

We now note the following trigonometric identities,

$$\sum_{j=1}^{q} \cos \frac{2\pi j}{q} = 0, \text{ and}$$

$$\sum_{j=1}^{q} \cos^2 \frac{2\pi j}{q} = \frac{\sum_{j=1}^{q} (1 + \cos \frac{4\pi j}{q})}{2} = \frac{q}{2}.$$

We get the second identity by using $\sum_{j=1}^{q} \cos \frac{4\pi j}{q} = 0$.

Thus, we get from Eq. 16

$$Jr - \frac{(Jr)^2 \beta}{2} = 0.$$

Noting that, $\frac{1}{\beta} = k_B T_{BKT}$, we finally get

$$T_{BKT} = \frac{Jr}{2k_B} \quad (17)$$

For a two-dimensional square lattice $r = 4$, hence, $T_{BKT} = \frac{2J}{k_B}$. We note here that the BKT transition
temperature is independent of $q$. We again estimate this temperature using the first-order MF theory in Sec. III.

### D. Calculation of $T_{SSB}$

As we already have seen in Sec. II B, there is a transition temperature at which the emergent $U(1)$ symmetry breaks. This temperature, denoted here by $T_{SSB}$, decreases quadratically over the variable $q$, i.e., $T_{SSB} = O(\frac{1}{q^2})$. This can also be understood in the following way. At very low temperatures, all the spins align in a particular direction. This ferromagnetic order is broken due to the low-lying excitations as we increase the temperature. The temperature of transition ($T_{SSB}$) corresponds to the minimum energy required to break the ferromagnetic order. We note that the larger the value of $q$, the less energy we need to break the ferromagnetic order. This shows that $T_{SSB}$ depends on $q$, decreasing as $q$ increases. To know how $T_{SSB}$ functionally depends on $q$, we must first calculate the energy required for a spin to make the smallest possible angle with its original direction in the ferromagnetic phase. Eq. 1 shows that this energy is proportional to $1 - \cos(2\pi/q)$, i.e., proportional to $1/q^2$ for large $q$. Thus, we infer that $T_{SSB} \propto 1/q^2$.

In the following we go on to estimate $T_{SSB}$ using yet another approach. We recall that, see Fig. 1, at low temperature (when the system is in ferromagnetic phase), all the spins spontaneously choose one of the $q$ free energy minima points. As we increase temperature, this ferromagnetic order will be broken when spins move from one minimum to the next. We note that while moving from one minimum to the next, the system needs to cross a barrier (local maximum). So we set the following criterion to estimate $T_{SSB}$: this is the temperature at which thermal energy at any of the minima equals the barrier height between two consecutive minima.

The order of the thermal energy is calculated from $TS = U - F$; for this purpose, we first separately calculate internal energy ($U$) and the free energy $F$) at one of the $q$ global minima. While calculating the energy barrier, we calculate the difference between one of the global minima and the local maximum between two consecutive global minima.

Calculating the energy barrier requires different strategies for even and odd values of $q$. We describe them separately below.

#### 1. Calculation of $T_{SSB}$ when $q$ is odd

For odd $q$, in the low-temperature regime (Fig. 1), we have one global minimum on the positive $l_1$-axis and one local maximum (between two consecutive global minima) on the negative $l_1$-axis. For illustration, we plot the low-temperature free energy as a function of $l_1$ keeping $l_2 = 0$ for the $q = 3$ model, as shown Fig. 2(a)). In this plot, we see a global minimum and a local minimum, which is

actually the local maximum while going from one to the next global minimum. We can calculate their positions along the $l_1$-axis by taking $\frac{\partial F(l_1, 0)}{\partial l_1} = 0$, where $F(l_1, 0)$ represents the free energy as function of $l_1$ keeping $l_2 = 0$.

For $q = 3$, according to Eq. 10, the free energy function $F(l_1, 0)$ is

$$F(l_1, 0) = \frac{Jl_1^2 r}{2} - \frac{\ln \sum_{i=1}^{3} e^{\beta J r(l_1 \cos \frac{2\pi i}{q})}}{\beta}. \quad (18)$$

Taking the partial derivative of $F(l_1, 0)$ w.r.t $l_1$, and equating it to zero, we get

$$l_1 = \frac{\sum_{i=1}^{3} \cos \frac{2\pi i}{3} e^{\beta l_1 J r \cos \frac{2\pi i}{3}}}{\sum_{i=1}^{3} e^{\beta l_1 J r \cos \frac{2\pi i}{3}}}. \quad (19)$$

Considering $r = 4$ and $J = 1$, we obtain

$$l_1 = \frac{e^{4\beta l_1} - e^{-2\beta l_1}}{e^{4\beta l_1} + 2e^{-2\beta l_1}}. \quad (20)$$

This is a transcendental equation that can not be solved directly. We need to do some numerical calculations to solve this equation. Let

$$g(l_1) = \frac{e^{4\beta l_1} - e^{-2\beta l_1}}{e^{4\beta l_1} + 2e^{-2\beta l_1}} - l_1. \quad (21)$$

For a given $\beta$, the values of $l_1$ corresponding to the zeros of the function $g(l_1)$ are the solutions of Eq. 20. These solutions can easily be found from the plot of $g(l_1)$, as can be seen in Fig. 2(b). We clearly see that there are three solutions: one corresponds to the global minimum ($l_1 > 0$), the next one corresponds to the local maximum ($l_1 = 0$), and the last one corresponds to the local maximum between two consecutive minima ($l_1 < 0$). Now, the energy barrier can be found by taking the difference between the global minimum and local maximum calculated from $F(l_1, 0)$ using appropriate $l_1$ values (solutions of Eq. 20). The energy barrier as a function of $\beta$ can be found in Fig. 2(c).

Now, for the calculation of thermal energy, we first calculate the internal energy of the spin by using the relation

$$U = \frac{\sum_{k=1}^{3} \epsilon(k) e^{-\beta \epsilon(k)}}{\sum_{k=1}^{3} e^{-\beta \epsilon(k)}}$$

where the expression of energy $\epsilon(k)$ can be found in Eq. 7 with $q = 3$. Written explicitly, we have

$$U = \frac{2(2l_1^2 + 2l_1)e^{-2l_1\beta} + (2l_1^2 - 4l_1)e^{4\beta l_1}}{2e^{-2\beta l_1} + e^{4\beta l_1}}. \quad (22)$$

The order of thermal energy ($U-F$) at a global minimum is calculated using Eqs. 18 and 22 with the appropriate $l_1$ value coming from solving Eq. 20. The thermal energy as a function of $\beta$ can be found in Fig. 2(c). From this plot, we can easily check the temperature at which the thermal energy and the energy barrier match.
<figure><img src="image_4.png" /><figcaption>FIG. 2: Plots of (a) $F(l_1, 0)$ vs. $l_1$ and (b) $g(l_1)$ vs. $l_1$ for $q = 3$ clock model with $\beta = 0.82$. (c) The thermal energy and the energy barrier are plotted as function of $\beta$ for $q = 3$. $T_{SSB}$ is estimated from where the two lines cross.</figcaption></figure>

We find that around $\beta = 0.65$, they are equal; thus, for $q = 3$, the symmetry-breaking phase transition temperature is $T_{SSB} = \frac{1}{0.65k_B} = \frac{1.538}{k_B}$. In the same way, we can calculate symmetry-breaking phase transition temperature for other odd $q$ values. For $q = 5$ and 7, calculated $T_{SSB}$ values can be found in Table I.

## 2. Calculation of $T_{SSB}$ when $q$ is even

While calculating $T_{SSB}$ for the $q$-state clock model with even $q$, the trick used for odd $q$ will not work. If we plot the function $F(l_1, 0)$ with $l_1$ or plot $F(0, l_2)$ with $l_2$, we will get two global minima, as can be seen in Fig. 3(b). However, we need the global minimum and the local maximum (between two successive global minima) to calculate the energy barrier between two successive global minima. To get the local maximum, we need to analyze the free energy along the line joining the origin and the midpoint between two consecutive global minima. For convenience, we can take the line joining the origin and the midpoint between the first two global minima (the first one being on the positive $l_1$-axis). The equation for this line is $l_2 = l_1 \tan \frac{\pi}{q}$. Along this line, the $F(l_1, l_2)$ plot can be seen in Fig. 3(a). Along this line, we can get the local maximum by solving the equation $\frac{\partial F}{\partial l_1} = 0$. Replacing $l_2$ by $l_1 \tan \frac{\pi}{q}$ in the expression of $F(l_1, l_2)$, then following the procedure described for odd $q$, we get the following transcendental equation for $l_1$ (taking $q = 4$ and $J = 1$):

$$l_1 = \frac{e^{4\beta l_1} - e^{-4\beta l_1}}{2(e^{4\beta l_1} + e^{-4\beta l_1})} \quad (23)$$

To solve this equation for $l_1$, we define the following function,

$$g_1(l_1) = \frac{e^{4\beta l_1} - e^{-4\beta l_1}}{2(e^{4\beta l_1} + e^{-4\beta l_1})} - l_1. \quad (24)$$

The roots of the equation $g_1(l_1) = 0$ are the solutions of Eq. 23. As can be seen from Fig. 3(c), there are three solutions. Two solutions for which $l_1 \neq 0$ correspond to the local maxima (between global minima).

There are two global minima on the $l_1$-axis, as can be seen in Fig. 3(b); we can get them by solving $\frac{\partial F(l_1, 0)}{\partial l_1} = 0$. This leads to the following transcendental equation for $l_1$ (again taking $q = 4$ and $J = 1$):

$$l_1 = \frac{e^{4\beta l_1} - e^{-4\beta l_1}}{2 + (e^{4\beta l_1} + e^{-4\beta l_1})} \quad (25)$$

We now define the following function $g_2(l_1)$ to find the solutions of the above equation,

$$g_2(l_1) = l_1 - \frac{e^{4\beta l_1} - e^{-4\beta l_1}}{2 + (e^{4\beta l_1} + e^{-4\beta l_1})} \quad (26)$$

We see from Fig. 3(d) that there are three solutions, two of them for which $l_1 \neq 0$ corresponds to the global minima. The energy barrier can be found by taking the difference between the global minimum and local maximum of the free energy as calculated above. The energy barrier as a function of $\beta$ can be seen in Fig. 4.

Next, to estimate the thermal energy ($U-F$) at the global minimum, we need to calculate the internal energy $U$. The expression of $U$, in this case, is given by

$$U = \frac{4l_1^2 + (2l_1^2 + 4l_1)e^{-4\beta l_1} + (2l_1^2 - 4l_1)e^{4\beta l_1}}{2 + e^{-4\beta l_1} + e^{4\beta l_1}} \quad (27)$$

where $l_1$ corresponds to a global minimum. The details of how the internal energy is calculated are already discussed for the case of odd $q$. The thermal energy at a global minimum as a function of $\beta$ can be found in Fig. 4. From this figure, we can easily check at which temperature the thermal energy matches with the energy barrier; we see that, for $q = 4$, two quantities match near $\beta = 0.93$. We, therefore, estimate that $T_{SSB} = 1.075/k_B$. We also calculated $T_{SSB}$ for $q = 6$, which can be found in Table I.

From this Table I, we can see that the spontaneous symmetry-breaking phase transition temperature ($T_{SSB}$) decreases with an increase in $q$; this is also evident in Fig. 5. Although these results help us understand how $T_{SSB}$ changes with $q$, they do not reveal the exact functional dependence of the transition temperature on $q$. For a better understanding of the functional dependence, in the following, we use our mean-field theory to analytically estimate the value of $T_{SSB}$ in the large-$q$ limit.
<figure><img src="image_5.png" /><figcaption>FIG. 3: Plots of (a) $F(l_1, l_2 = l_1)$ vs. $l_1$, (b) $F(l_1, 0)$ vs. $l_1$, (c) $g_1(l_1)$ vs. $l_1$ and (d) $g_2(l_1)$ vs. $l_1$ for $q = 4$ clock model with $\beta = 0.82$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 4: The thermal energy and the energy barrier are plotted as function of $\beta$ for $q = 4$. $T_{SSB}$ is estimated from where the two lines cross.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIG. 5: Plot of $T_{SSB}$ vs. $q$</figcaption></figure>

<table><caption>TABLE I: Numerically calculated $T_{SSB}$ values for different $q$ values</caption><thead><tr><td>$q$-value</td><td>$T_{SSB}$</td></tr></thead><tbody><tr><td>3</td><td>1.538/$k_B$</td></tr><tr><td>4</td><td>1.075/$k_B$</td></tr><tr><td>5</td><td>0.7692/$k_B$</td></tr><tr><td>6</td><td>0.556/$k_B$</td></tr><tr><td>7</td><td>0.4225/$k_B$</td></tr></tbody></table>

### E. Analytical calculation of $T_{SSB}$ in the large-$q$ limit

In this part, we first calculate the energy barrier, and then we calculate the thermal energy to estimate $T_{SSB}$ in the large-$q$ limit.

#### 1. Calculation of energy barrier

For the large-$q$ limit, we can convert the summation appearing in the free energy expression, as seen in Eq.10, into an integration. Consequently, we have the following expression of the free energy,

$$F = \frac{J(l_1^2 + l_2^2)r}{2} - \frac{1}{\beta} \ln \int_1^q e^{\beta J r (l_1 \cos \frac{2\pi j}{q} + l_2 \sin \frac{2\pi j}{q})} dj. \quad (28)$$

However, performing this integration is not straightforward. To go further, we replace the integrand function with a judiciously chosen simpler function. We first convert the integration limits to a more symmetric ones; ac-
cordingly,

$$F = \frac{J(l_1^2 + l_2^2)r}{2} - \frac{1}{\beta} \ln \int_{-\frac{q}{2}+1}^{\frac{q}{2}} e^{\beta J r(l_1 \cos \frac{2\pi j}{q} + l_2 \sin \frac{2\pi j}{q})} dj. \quad (29)$$

We only consider the even $q$ for the current calculations. If we are interested in calculating the global minima on the $l_1$-axis, we have to deal with the integrand function $e^{\beta J r l_1 \cos \frac{2\pi j}{q}}$. While replacing this function with a simpler function, we choose a function so that the extrema points of the original and the new function match. We note that the original function has a maximum at $j = 0$ and the function effectively vanishes when $j > \frac{q}{4}$ or $j < -\frac{q}{4}$, since $\cos \frac{2\pi j}{q}$ picks up a negative sign in the ranges. With these properties in mind, we now replace this integrand function with the following approximate function:

$$e^{\beta J r l_1 \cos \frac{2\pi j}{q}} \rightarrow (1 - c j^2) e^{\beta J r l_1}.$$

The value of $c$ is determined by the fact that the original function effectively vanishes at $j = \pm \frac{q}{4}$. This gives us $c = (\frac{4}{q})^2$. It may be mentioned here that the original and the approximate functions have the same maximum, and they effectively vanish at the same points. It is also worth commenting that the integration over the new approximate function should be performed only from $-\frac{q}{4}$ to $\frac{q}{4}$. With this judiciously chosen approximate function, it is now easy to estimate the value of the integral:

$$\int_{-\frac{q}{2}+1}^{\frac{q}{2}} e^{\beta J r (l_1 \cos \frac{2\pi j}{q})} dj \rightarrow \int_{-\frac{q}{4}}^{\frac{q}{4}} \left(1 - \left(\frac{4}{q}\right)^2 j^2\right) e^{\beta J r l_1} dj \\ = \frac{q}{3} e^{J r \beta l_1}.$$

Using this result, the expression for the free energy $F$ along the $l_1$-axis can be written as,

$$F = \frac{J l_1^2 r}{2} - J r l_1 - \frac{\ln\left(\frac{q}{3}\right)}{\beta}. \quad (30)$$

$$\int_{-\frac{q}{2}+1}^{\frac{q}{2}} e^{\beta J r l \cos \frac{\pi}{q} \cos \frac{(2\pi j - \pi)}{q}} dj \rightarrow \int_{-\frac{q}{4}+\frac{1}{2}}^{0} (1 - c j^2) e^{(J r \beta l \cos^2 \frac{\pi}{q})} dj + \int_{1}^{\frac{q}{4}+\frac{1}{2}} (1 - b(j-1)^2) e^{(J r \beta l \cos^2 \frac{\pi}{q})} dj. \quad (33)$$

After performing the integration, we get the following:

$$\int_{-\frac{q}{2}+1}^{\frac{q}{2}} e^{\beta J r l \cos \frac{\pi}{q} \cos \frac{(2\pi j - \pi)}{q}} dj \rightarrow \frac{q-2}{3} e^{(J r \beta l \cos^2 \frac{\pi}{q})}. \quad (34)$$

Thus, from Eq. 32, we get the following expression for

Now the global minimum can be found from $\frac{\partial F}{\partial l_1} = 0$. This gives us $l_1 = 1$ at the global minimum. So, the free energy at the global minimum is,

$$F_G = -\frac{J r}{2} - \frac{1}{\beta} \ln\left(\frac{q}{3}\right). \quad (31)$$

For the calculation of the energy barrier, we next calculate the local maximum (between consecutive global minima). We note that if we connect by straight lines all the global minima to the origin ($l_1=0, l_2=0$), the angle between consecutive lines would be $\frac{2\pi}{q}$. Let us consider that the coordinate of the local maximum, between the first two global minima, is ($l_1, l_2$) and the coordinate of the global minimum, which is on the $l_1$-axis, is given by ($l, 0$). Some trigonometric considerations show that $l_1 = l \cos^2 \frac{\pi}{q}$ and $l_2 = l \cos \frac{\pi}{q} \sin \frac{\pi}{q}$. Then, the expression for free energy, as seen in Eq. 29, at the local maximum will take the following form:

$$F = \frac{J r l^2 \cos^2 \frac{\pi}{q}}{2} - \frac{1}{\beta} \ln \int_{-\frac{q}{2}+1}^{\frac{q}{2}} e^{\beta J r l \cos \frac{\pi}{q} \cos \frac{(2\pi j - \pi)}{q}} dj. \quad (32)$$

Once again, we will replace the integrand function with a simpler function with the same extreme points. We note that the function $e^{\beta J r l \cos \frac{\pi}{q} \cos \frac{(2\pi j - \pi)}{q}}$ effectively vanishes when $j < -\frac{q}{4} + \frac{1}{2}$ or $j > \frac{q}{4} + \frac{1}{2}$, and has a maximum at $j = \frac{1}{2}$. But $j$ is always an integer, so we have to consider the nearest integer of $\frac{1}{2}$, which are 0 and 1, while approximating the function; accordingly, the new function will be defined in two regions: one is $(-\frac{q}{4} + \frac{1}{2})$ to 0, and the other one is from 1 to $(\frac{q}{4} + \frac{1}{2})$.

In the first interval, i.e. from $(-\frac{q}{4} + \frac{1}{2})$ to 0, we can approximate the function as $(1 - c j^2) e^{(J r \beta l \cos^2 \frac{\pi}{q})}$. Both the original function and this approximate function have the same maximum at $j = 0$. The constant $c$ can be determined from the fact that the approximate function should vanish at $j = -\frac{q}{4} + \frac{1}{2}$. From this condition we get $c = (\frac{4}{q-2})^2$.

Similarly, for the second interval from $j=1$ to $j=(\frac{q}{4} + \frac{1}{2})$, we approximate the function as $(1 - b(j-1)^2) e^{(J r \beta l \cos^2 \frac{\pi}{q})}$ with $b = (\frac{4}{q-2})^2$. Thus, we get the following:

the free energy at the local maximum,

$$F_L = \frac{J r l^2 \cos^2 \frac{\pi}{q}}{2} - \frac{\ln\left(\frac{q-2}{3} e^{(J r \beta l \cos^2 \frac{\pi}{q})}\right)}{\beta}, \quad (35)$$

with $l = 1$, as obtained earlier. Taking $\cos^2 \frac{\pi}{q} \approx 1 - \frac{\pi^2}{q^2}$
for large $q$, we get

$$F_L = -\frac{1}{2}Jr\left(1 - \frac{\pi^2}{q^2}\right) - \frac{1}{\beta}\ln\left(\frac{q-2}{3}\right). \quad (36)$$

From Eqs. 31 and 36, we now get an estimate for the energy barrier:

$$F_L - F_G = \frac{Jr\pi^2}{2q^2} - \frac{\ln \frac{q-2}{3}}{\beta} + \frac{\ln \frac{q}{3}}{\beta}. \quad (37)$$

## 2. Calculation of thermal energy

To calculate the thermal energy ($U-F$) at a global minimum, we use the following thermodynamic relation to first calculate the internal energy at the global minimum ($U_G$):

$$U_G = -T^2 \frac{\partial}{\partial T} \left( \frac{F_G}{T} \right), \quad (38)$$

where $F_G$ is given by Eq. 31. A simple calculation shows that

$$U_G = -\frac{Jr}{2}. \quad (39)$$

Hence, the thermal energy at the global minimum is

$$U_G - F_G = \frac{1}{\beta} \ln\left(\frac{q}{3}\right). \quad (40)$$

## 3. Estimation of $T_{SSB}$

As discussed earlier, an estimation of the spontaneous symmetry-breaking transition temperature ($T_{SSB}$) can be found by equating the energy barrier between two consecutive global minima and the thermal energy at a global minimum. Accordingly, from Eqs. 37 and 40, we get

$$\frac{Jr\pi^2}{2q^2} - \frac{\ln \frac{q-2}{3}}{\beta} + \frac{\ln \frac{q}{3}}{\beta} = \frac{1}{\beta} \ln\left(\frac{q}{3}\right).$$

Solution of this equation gives (after replacing $\beta$ by $\frac{1}{k_B T_{SSB}}$),

$$T_{SSB} = \frac{Jr\pi^2}{2k_B q^2 \ln \frac{q-2}{3}}. \quad (41)$$

For large-$q$, $q^2$ dominates over $\ln \frac{q-2}{3}$, so we see that $T_{SSB}$ decreases as $q^{-2}$, as expected earlier for large-$q$.

The zeroth order mean-field theory that we discussed provides many insights into the $q$-state clock model. However, this theory does not allow us to estimate the important quantities like spin-spin correlation. The spin-spin correlation helps us gain more insight into the phases

<figure><img src="image_8.png" /><figcaption>FIG. 6: First order mean-field configuration</figcaption></figure>

of the model. To calculate this quantity, we now develop a higher-order (first order) mean-field theory where two neighboring spins with exact interaction are targeted. This higher-order theory also helps us improve some of the results we already obtained using the zeroth order theory.

## III. HIGHER-ORDER MEAN-FIELD THEORY

Let us consider a two-dimensional square lattice and focus on the two spins ($S_0$ and $S_1$) that interact exactly. These two spins are in the mean field created by their neighboring spins ($S_2, S_3, S_4, S_5, S_6, S_7$), as shown in Fig. 6. The Hamiltonian for the system of two spins can be written as:

$$h_{01} = -J \mathrm{Re}(S_0 S_1^* + \frac{1}{2}(S_5 S_0^* + S_1 S_2^* + S_6 S_0^* + S_0 S_4^* + S_7 S_1^* + S_1 S_3^*)). \quad (42)$$

The factor 1/2 is to avoid double counting the interactions between the targeted spins and their neighbors. If there are $N$ number of spins, then there are $N/2$ such spin pairs. We recall here that the complex spin variable ($S_k$) associated with the $k$-th spin is by definition: $S_k = e^{i\theta_k}$, where $\theta_k = 2\pi j_k/q$. The variable $j_k$ takes any integer value from 1 to $q$. To develop the MF theory, we now replace a spin variable with a complex mean-field parameter ($l$) and a small fluctuation or deviation term: $S_k = l + \delta S_k$. We take a translationally invariant system for which the MF parameter $l$ is the same for all the spins. Using the parameter $l$, any interaction term can be approximated as: $S_i S_j^* \approx -|l|^2 + l^* S_i + l S_j^*$. With the help of this expression, we get the following mean-field version of Hamiltonian in Eq. 42,

$$\begin{aligned} h_{MF} = & -J \mathrm{Re}(S_0 S_1^* \\ & -\frac{3}{2}|l|^2 + \frac{3}{2}l^* S_0 + \frac{1}{2} \sum_{j \in N_0} l S_j^* \\ & -\frac{3}{2}|l|^2 + \frac{3}{2}l^* S_1 + \frac{1}{2} \sum_{k \in N_1} l S_k^*). \end{aligned} \quad (43)$$
Here, $N_0$ denotes the set of nearest neighbors of the 0-th spin except the 1st spin. Similarly, $N_1$ is the set of nearest neighbors of the 1st spin except the 0-th spin. It may be mentioned here that since we are interested in the real part of an interaction term, terms like $l^*S_0$ and $lS_0^*$ can be taken to be equivalent and can be written interchangeably. From Eq. 43, if we now write the mean-field Hamiltonian for the whole two-dimensional system, we will have exact interaction terms like $S_0S_1^*$ corresponding to each target spin pair and the terms like $l^*S_0$ appearing 3 times in the Hamiltonian. Therefore, the total mean-field Hamiltonian of the system can be written as the sum of $N/2$ mean-field pair Hamiltonians; this mean-field pair Hamiltonian corresponding to $i-j$ pair (nearest neighbors) is given by,

$$h_{ij} = -J \mathrm{Re}(S_i S_j^* - 3|l|^2 + 3l^* S_i + 3l S_j^*). \quad (44)$$

Now we recall that $S_k = \cos \theta_k + i \sin \theta_k$, where $\theta_k = 2\pi j_k/q$ with $j_k$ taking values from 1 to $q$. If we write $l = l_1 + il_2$, we get from Eq. 44:

$$h_{ij} = -J \cos(\theta_i - \theta_j) + 3J(l_1^2 + l_2^2) - 3Jl_1(\cos \theta_i + \cos \theta_j) - 3Jl_2(\sin \theta_i + \sin \theta_j). \quad (45)$$

$$F = 3J(l_1^2 + l_2^2) - \frac{1}{\beta} \ln \sum_{i=1}^{q} \sum_{j=1}^{q} e^{J\beta(\cos \frac{2\pi(i-j)}{q} + 3l_1(\cos \frac{2\pi i}{q} + \cos \frac{2\pi j}{q}) + 3l_2(\sin \frac{2\pi i}{q} + \sin \frac{2\pi j}{q})).} \quad (48)$$

### A. Free energy plots at different temperatures

Using this new approximation, the behavior of free energy can be analyzed at a high temperature, intermediate temperature, and low temperature in a similar way as we have done previously (within the zeroth order mean-field approach). The free energy plots in three different temperature ranges for $q=5, 6$, and $7$ can be found in Fig. 7. These new results are qualitatively the same as obtained using the zeroth order theory. There are some quantitative differences, which can be appreciated from the updated value of $T_{BKT}$, which we again calculate in the following using the new mean-field approach.

### B. Calculation of $T_{BKT}$ using new mean-field approximation

We can use the same scheme to calculate the BKT transition temperature as in the first part (Sec. II). We will find at which temperature the free energy function has an inflection point at $(0,0)$. From Eq. 48, we have

The possible energy values of the spin pair can be obtained from Eq. 45 by taking different possible values of $\theta_i$ and $\theta_j$. We can now calculate the partition function of the system as $Z = Z_{ij}^{N/2}$, where $Z_{ij}$ is the partition function corresponding to the pair $i-j$ and is given by,

$$Z_{ij} = \sum_{\theta_i, \theta_j} e^{-\beta h_{ij}}. \quad (46)$$

Here, the sum is over possible values of $\theta_i$ and $\theta_j$. The free energy corresponding to a spin pair is given by,

$$F = -\frac{1}{\beta} \ln Z_{ij}. \quad (47)$$

Written explicitly, we get the following expression for the free energy,

the free energy as function of $l_1$ ($l_2=0$) as follows:

$$F(l_1, 0) = 3J|l_1|^2 - \frac{1}{\beta} \ln \sum_{i=1}^{q} \sum_{j=1}^{q} e^{J\beta(\cos \frac{2\pi(i-j)}{q} + 3l_1(\cos \frac{2\pi i}{q} + \cos \frac{2\pi j}{q}))}. \quad (49)$$

Now taking the second order partial derivative of $F(l_1, 0)$ at $l_1 = 0$ and equating it to 0, we obtain:

$$\begin{aligned} \left(\frac{\partial^2 F}{\partial l_1^2}\right)_{l_1=0} &= 6J \\ &+ 9\beta J^2 \frac{\left(\sum_{i=1}^q \sum_{j=1}^q (\cos \frac{2\pi i}{q} + \cos \frac{2\pi j}{q}) e^{\beta J \cos \frac{2\pi(i-j)}{q}}\right)^2}{\left(\sum_{i=1}^q \sum_{j=1}^q e^{\beta J \cos \frac{2\pi(i-j)}{q}}\right)^2} \\ &- 9\beta J^2 \frac{\sum_{i=1}^q \sum_{j=1}^q (\cos \frac{2\pi i}{q} + \cos \frac{2\pi j}{q})^2 e^{\beta J \cos \frac{2\pi(i-j)}{q}}}{\sum_{i=1}^q \sum_{j=1}^q e^{\beta J \cos \frac{2\pi(i-j)}{q}}} = 0. \end{aligned} \quad (50)$$

It is a nontrivial task to find $\beta$ for which the above equation (Eq. 50) is satisfied. To go further, we perform the following approximation. As seen in the first part, in Sec. II, $\beta J \sim \frac{1}{2}$ near the BKT transition temperature. Therefore, near the transition temperature, $|\beta J \cos \frac{2\pi(i-j)}{q}| \lesssim$
<figure><img src="image_9.png" /><figcaption>FIG. 7: First order mean-field results: free energy plots at low, intermediate and high temperatures are shown in (a)-(c), (d)-(f) and (g)-(i), respectively. Results in the left, middle and right columns are for $q = 5$, 6 and 7, respectively.</figcaption></figure>

$\frac{1}{2}$. Accordingly, we take the following first order approximation: $e^{\beta J \cos \frac{2\pi(i-j)}{q}} \approx 1 + \beta J \cos \frac{2\pi(i-j)}{q}$. Noting that $\sum_{i=1}^{q} \cos \frac{2\pi i}{q} = 0$ and $\sum_{i=1}^{q} \cos^2 \frac{2\pi i}{q} = \frac{q}{2}$, we get the following:

$\beta$ by $\frac{1}{k_B T_{BKT}}$),

$$T_{BKT} = \frac{1.895J}{k_B}. \qquad (53)$$

$$ \begin{aligned} & \frac{\sum_{i=1}^{q} \sum_{j=1}^{q} (\cos \frac{2\pi i}{q} + \cos \frac{2\pi j}{q})^2 e^{\beta J \cos \frac{2\pi(i-j)}{q}}}{\sum_{i=1}^{q} \sum_{j=1}^{q} e^{\beta J \cos \frac{2\pi(i-j)}{q}}} \\ &= \frac{\sum_{i=1}^{q} \sum_{j=1}^{q} (\cos \frac{2\pi i}{q} + \cos \frac{2\pi j}{q})^2 (1 + \beta J \cos \frac{2\pi(i-j)}{q})}{\sum_{i=1}^{q} \sum_{j=1}^{q} (1 + \beta J \cos \frac{2\pi(i-j)}{q})} \\ &= \frac{q \cdot \frac{q}{2} + q \cdot \frac{q}{2} + 2\beta J \cdot \frac{q}{2} \cdot \frac{q}{2}}{q \cdot q} \\ &= 1 + \frac{\beta J}{2}. \end{aligned} \qquad (51) $$

Using this result, we get from Eq. 50,

$$ \begin{aligned} & 6J + 9\beta J^2 \times 0 - 9\beta J^2(1 + \frac{\beta J}{2}) = 0, \text{ i.e.,} \\ & 4 - 6\beta J - 3(\beta J)^2 = 0. \end{aligned} \qquad (52) $$

Taking only the positive root of Eq. 52, we get the following value for the transition temperature (after replacing

Instead of solving Eq. 50 analytically by taking some approximation, it is also possible to directly solve the equation numerically. For given $q$, we can numerically check at which $\beta$ the equation is satisfied. For $q = 6$ and 7, plots of $(\frac{\partial^2 F}{\partial l_1^2})_{l_1=0}$ vs. $\beta$ can be seen Fig. 8. Numerical value of $T_{BKT}$ calculated from such plots show that, for $q = 7$ or more, the value of the transition temperature does not change much with $q$ and this constant value is estimated to be $T_{BKT} = \frac{1.895J}{k_B}$. We note that, compared to the zeroth order result, the current value of $T_{BKT}$ is slightly closer to the reported value of $0.893J/k_B$ as found using Monte Carlo calculations [30]. Thus, we see that the consideration of a single exact interaction in the higher-order mean-field theory improves the results quantitatively.
<figure><img src="image_10.png" /><figcaption>FIG. 8: Plots of $(\partial^2 F / \partial l_1^2)_{l_1=0}$ vs. $\beta$ for (a) $q=6$ and (b) $q=7$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIG. 9: NN spin-spin correlation for (a) $q=6$, and (b) $q=7$.</figcaption></figure>

### C. Estimation of the spin-spin correlation

The spin-spin correlation is an important quantity in studying different phases of a system. The correlation shows different behaviors in different phases, and a change in its behavior normally indicates a phase transition. Mean-field theories, in general, do not predict exact correlation behavior. Our first-order mean-field theory considers one exact interaction between two neighboring spin pairs. As a result, it is possible to estimate the nearest-neighbor (NN) spin-spin correlation using this approach. The NN correlation can tell us whether the system is in an ordered or disordered phase at a given temperature. Using our first-order MF theory, we estimate the NN spin-spin correlation for the q-state clock model. We note that the NN spin-spin correlation is given by,

$$ \langle \mathrm{Re}(S_0 S_1^*) \rangle = \frac{\sum_{\theta_0} \sum_{\theta_1} \mathrm{Re}(S_0 S_1^*) e^{-\beta h_{01}}}{\sum_{\theta_0} \sum_{\theta_1} e^{-\beta h_{01}}}. \quad (54) $$

Here $h_{01}$ is given by Eq. 45 with indices i and j taking values 0 and 1 respectively. In terms of variables $\theta_0$ and $\theta_1$, we also have $\mathrm{Re}(S_0 S_1^*) = \cos(\theta_0 - \theta_1)$. Evaluating this correlation analytically is, in general, difficult. It is, though, possible to numerically evaluate the value of NN correlation at different values of $\beta$. In the following we first discuss how one can calculate NN correlation numerically, we then discuss how one can get analytical results in two different limits (high and low temperature limits).

### 1. Numerical calculation of NN correlation

To calculate NN correlation as function of $\beta$, we first need to find the appropriate numerical values of the parameters $l_1$ and $l_2$ at different $\beta$ values, since these parameters appear in Eq. 54 via $h_{01}$. We note from the free energy plots (Fig. 7) that there is always a minimum on the positive $l_1$-axis (including origin). In the high-temperature limit, the minimum is at the origin ($l_1 = l_2 = 0$). In the low-temperature limit, there are $q$ minima around the origin; one of these minima is on the positive $l_1$-axis. In the intermediate temperature, there is a ring of minima around the origin. This observation allows us to take $l_2 = 0$ and an appropriate value of $l_1$ (corresponding to the minimum on the positive $l_1$-axis) in $h_{01}$. To find this appropriate value of $l_1$ at a given $\beta$, we numerically find the free energy minimum on the positive $l_1$-axis from the equation $\partial F(l_1, 0) / \partial l_1 = 0$, where $F(l_1, 0)$ is given by Eq. 49.

For $q=6$ and $7$, we calculate $l_1$ corresponding to the minimum on the $l_1$-axis for different $\beta$ values and then use that to calculate NN correlation at those $\beta$ values. Results can be seen in Fig. 9. We observe from the plots of NN correlation that there is a kink near $T_{BKT} = 1.895J / k_B$. This result correctly predicts a change in phase at $T_{BKT}$. We, though, do not observe any noticeable change in correlation at $T_{SSB}$.
## 2. Analytical expression of NN correlation in high-temperature limit

We know from the free energy plots (Fig. 7) that, in the high-temperature limit, there is only one minimum corresponding to $l_1 = l_2 = 0$. Accordingly, we can evaluate Eq. 54 by taking $l_1 = l_2 = 0$ in the Hamiltonian $h_{01}$. Thus, in this temperature limit, $h_{01} = -J \cos(\theta_0 - \theta_1)$ and $e^{-\beta h_{01}} \approx 1 + \beta J \cos(\theta_0 - \theta_1)$, since $\beta$ is small. Written explicitly, we get from Eq. 54,

$$ \langle \mathrm{Re}(S_0 S_1^*) \rangle = \frac{\sum_{\theta_0} \sum_{\theta_1} \cos(\theta_0 - \theta_1)(1 + \beta J \cos(\theta_0 - \theta_1))}{\sum_{\theta_0} \sum_{\theta_1} (1 + \beta J \cos(\theta_0 - \theta_1))} \\ = \frac{\sum_{i=1}^q \sum_{j=1}^q \cos \frac{2\pi(i-j)}{q} (1 + \beta J \cos \frac{2\pi(i-j)}{q})}{\sum_{i=1}^q \sum_{j=1}^q (1 + \beta J \cos \frac{2\pi(i-j)}{q})}. \quad (55) $$

Noting the identities $\sum_{i=1}^q \cos \frac{2\pi i}{q} = 0$ and $\sum_{i=1}^q \cos^2 \frac{2\pi i}{q} = \frac{q}{2}$, we finally get

$$ \langle \mathrm{Re}(S_0 S_1^*) \rangle = \frac{\beta J}{2}. \quad (56) $$

We note here that, in the high-temperature limit, the NN spin-spin correlation decreases with increasing temperature. In accordance with our analytical result, we also note from Fig. 9 that the numerically calculated NN correlation is linear with $\beta$ near the origin.

## 3. Analytical expression of NN correlation in low-temperature limit

In the low-temperature limit ($\beta$ being large), $l_1 \approx 1$ corresponds to the free energy minimum on the positive $l_1$ axis. In this limit, it is not possible to approximate $e^{-\beta h_{01}}$ by retaining a couple of terms from its (Taylor) series. To go further, for large $q$, we transform the summation to integration. With $l_1 = 1$ and $l_2 = 0$ in Eq. 45, we have

$$ h_{01} = -J \cos(\theta_0 - \theta_1) + 3J - 3J(\cos \theta_0 + \cos \theta_1), \quad (57) $$

with the indices $i$ and $j$ taking their values 0 and 1, respectively. In the large $q$ limit and taking symmetric ranges for the integration variables $\theta_0$ and $\theta_1$, we get from Eq. 54,

$$ \langle \mathrm{Re}(S_0 S_1^*) \rangle = \frac{\int_{\theta_0=-\pi}^{\pi} \int_{\theta_1=-\pi}^{\pi} \cos(\theta_0 - \theta_1) e^{-\beta h_{01}} d\theta_0 d\theta_1}{\int_{\theta_0=-\pi}^{\pi} \int_{\theta_1=-\pi}^{\pi} e^{-\beta h_{01}} d\theta_0 d\theta_1}, \quad (58) $$

where $h_{01}$ is given by Eq. 57. This can be further written as:

$$ \langle \mathrm{Re}(S_0 S_1^*) \rangle = \lim_{x \to 1} \frac{1}{\beta J} \frac{\partial}{\partial x} \ln \left( \int_{-\pi}^{\pi} \int_{-\pi}^{\pi} E(\theta_0, \theta_1; x) d\theta_0 d\theta_1 \right), \quad (59) $$

where $E(\theta_0, \theta_1; x) = e^{J\beta(x \cos(\theta_0 - \theta_1) + 3(\cos\theta_0 + \cos\theta_1))}$. We now approximate the function $E(\theta_0, \theta_1; x)$ in a manner similar to how in the first part, we approximated an exponential function by a second-order polynomial. We note that function like $e^{t \cos \theta}$, with $t > 0$, has a maximum $e^t$ at $\theta = 0$. This exponential function effectively vanishes when $\theta \le -\frac{\pi}{2}$ or $\theta \ge \frac{\pi}{2}$. For large $t$, this allows us to approximate the function as $e^{t \cos \theta} \approx e^{t(1-\theta^2/2)} \to e^t(1 - tA\theta^2)$. If $t$ is a given constant, $A$ is determined by the fact that at $\theta = \pm \frac{\pi}{2}$, the approximate function should vanish. If $t$ is a variable, then $A$ is determined by replacing $t$ by its limiting value and making sure that the approximate function vanishes at $\theta = \pm \frac{\pi}{2}$. It may be noted here that the original exponential function and the approximate polynomial function both have the same maximum at $\theta = 0$.

Above discussion helps us to approximate $E(\theta_0, \theta_1; x)$ in the following way:

$$ E(\theta_0, \theta_1; x) \to e^{\beta J(x+6)} (1 - xA(\theta_0 - \theta_1)^2) (1 - B\theta_0^2) (1 - C\theta_1^2). \quad (60) $$

With this approximation, the integration ranges in Eq. 59 should be taken appropriately. First, we now determine the values of the constants $A$, $B$, and $C$. The constants $B$ and $C$ are determined by making sure that the approximate function should vanish at $\theta_0 = \pm \frac{\pi}{2}$ and $\theta_1 = \pm \frac{\pi}{2}$. This gives us $B = C = (\frac{2}{\pi})^2$. While determining the constant $A$, we note that while $\theta_0 = \theta_1$ corresponds to the maximum for both exponential and the approximate polynomial function, both functions should vanish when $|\theta_0 - \theta_1| \ge \frac{\pi}{2}$. Since finally we are going to take the limit $x \to 1$, we determine that $A = (\frac{2}{\pi})^2$.

To find the appropriate ranges of integrations in Eq. 59 after replacing $E(\theta_0, \theta_1; x)$ by its approximation, we first note that the function $e^{3\beta J \cos \theta_0}$ or $e^{3\beta J \cos \theta_1}$ effectively vanishes when $|\theta_0| \ge \frac{\pi}{2}$ or $|\theta_1| \ge \frac{\pi}{2}$. On the other hand, the function $e^{x\beta J \cos(\theta_0 - \theta_1)}$ effectively vanishes when $|\theta_0 - \theta_1| \ge \frac{\pi}{2}$. This tells us that, when $-\frac{\pi}{2} \le \theta_0 \le 0$, the range of the other variable should be $-\frac{\pi}{2} \le \theta_1 \le \frac{\pi}{2} + \theta_0$. Similarly, when $0 \le \theta_0 \le \frac{\pi}{2}$, the range of the other variable should be $-\frac{\pi}{2} + \theta_0 \le \theta_1 \le \frac{\pi}{2}$.

Let us denote the approximate function by $\tilde{E}(\theta_0, \theta_1; x)$, i.e.,

$$ \tilde{E}(\theta_0, \theta_1; x) = e^{\beta J(x+6)} (1 - \frac{4x}{\pi^2}(\theta_0 - \theta_1)^2) (1 - \frac{4\theta_0^2}{\pi^2}) (1 - \frac{4\theta_1^2}{\pi^2}). \quad (61) $$

Now to find the NN correlation using Eq. 59, we need to perform the following integration,

$$ I(x) = \int_{\theta_0=-\frac{\pi}{2}}^{0} \int_{\theta_1=-\frac{\pi}{2}}^{\frac{\pi}{2}+\theta_0} \tilde{E}(\theta_0, \theta_1; x) d\theta_0 d\theta_1 \\ + \int_{\theta_0=0}^{\frac{\pi}{2}} \int_{\theta_1=-\frac{\pi}{2}+\theta_0}^{\frac{\pi}{2}} \tilde{E}(\theta_0, \theta_1; x) d\theta_0 d\theta_1. \quad (62) $$

We find that $I(x) = (3.863 - 0.9678x)e^{\beta J(x+6)}$. We now
have,

$$ \begin{aligned} \langle \mathrm{Re}(S_0^* S_1) \rangle &= \lim_{x \to 1} \frac{1}{\beta J} \frac{\partial \ln I(x)}{\partial x} \\ &= 1 - \frac{1}{3\beta J}. \end{aligned} \qquad (63) $$

We see that the NN spin-spin correlation saturates to the value 1 as the temperature reduces (i.e., $\beta \to \infty$). This analytical result matches well with our numerical results shown in Fig. 9.

### 4. Characterizing phases using NN correlation

We noticed from our numerical as well as analytical results that $\lim_{\beta \to 0} \langle \mathrm{Re}(S_0^* S_1) \rangle = 0$. This result signifies that, at high-temperature limit, the neighboring spins are not correlated or they move randomly. This high-temperature phase is, then, a disordered phase, as expected.

Similarly, in the other limit, we find that $\lim_{\beta \to \infty} \langle \mathrm{Re}(S_0^* S_1) \rangle = 1$. This result says that, in low-temperature limit, the neighboring spins are maximally correlated. For a ferromagnetic model, this implies a low-temperature ordered phase (all spins orient along the same direction).

In general, in the intermediate temperature regime, $0 < \langle \mathrm{Re}(S_0^* S_1) \rangle < 1$. This result says that the neighboring spins are not parallel and they, on average, make some angle with each other. This finding, along with the free energy results as shown in Figs. 7 (d)-(f), implies that, in the intermediate temperature limit, an individual spin orient along a particular direction but different neighboring spins orient in such a way that we get the vortex-antivortex topological defects (as seen in XY model). Here we also refer to our discussion in Sec. II A.

we develop a mean-field theory (both basic and higher order) to study the q-state clock model. To the best of our knowledge, this is the first systematic mean-field study of the model with finite q.

Our work provides a comprehensive analysis of the q-state clock model in two-dimension. Our findings reaffirm the existence of three different phases: (broken) $Z_q$ symmetric ferromagnetic (ordered) phase at the low-temperature, emergent $U(1)$ symmetric BKT phase at the intermediate temperature and paramagnetic (disordered) phase at the high-temperature. We find that the transition at the higher temperature is of BKT type. We also argue here that the transition at the lower temperature is of the large-order SSB type. Our MF theory helps us confirm that the transition temperature $T_{BKT}$ does not depend on q value (for large q) and the transition temperature $T_{SSB}$ decreases quadratically with q.

The higher-order mean-field approach developed here refines our understanding by incorporating exact treatments of nearest-neighbor interactions. This approach enhances our ability to describe the model's behavior with better accuracy (resulting in, for example, slightly better estimation for $T_{BKT}$). Using this approach, we could estimate NN spin-spin correlation, which then helped us better characterize the phases.

We know that there exists a critical value of q (denoted as $q_c$) for the 2D clock model, below which the model only shows a SSB transition and above which the model shows both SSB and BKT transition. Although our mean-field theory did not determine this critical value precisely, it serves as a robust framework for exploring the complex phase structure of the q-state clock model. Future research could build upon these findings to investigate the model further: as example, our mean-field approach can be used in advancing our understanding of the clock model in higher dimensions and in studying the model under applied field.

## IV. CONCLUSION

The q-state clock model garnered a lot of interest due to the fact that it shows a double phase transitions for finite $q \ge 5$ in two dimension. There are varying opinions about the nature of phases and associated phase transitions, especially concerning the nature of the transition at lower temperature. Except for a few special values of q, the 2D model is not exactly solvable. Different numerical techniques are developed to gain insights into the properties of this model. To investigate some of the unanswered questions and to address some of the debates in this field,

## ACKNOWLEDGMENTS

We acknowledge Smitarani Mishra for her participation in discussions in the initial phase of this project.

**Author contributions:** RK and MG calculated the BKT transition temperature using the zeroth order MF theory, AG performed most of the remaining calculations using the zeroth and higher order MF theories, SS proposed and closely supervised the project besides completing some calculations, SS wrote the manuscript with the help from AG.

[1] L. Onsager, Crystal statistics. I. a two-dimensional model with an order-disorder transition, Phys. Rev. 65, 117 (1944).

[2] R. K. Pathria and P. D. Beale, Statistical Mechanics, 3rd ed. (Elsevier, 2011).

[3] N. D. Mermin and H. Wagner, Absence of ferromag-
netism or antiferromagnetism in one- or two-dimensional isotropic heisenberg models, Phys. Rev. Lett. 17, 1133 (1966).

[4] Berezinskii, Destruction of long-range order in one-dimensional and two-dimensional systems having a continuous symmetry group I. classical systems, Sov. Phys. JETP 32, 493 (1971).

[5] J. M. Kosterlitz, Kosterlitz-Thouless physics: a review of key issues, Reports on Progress in Physics 79, 026001 (2016).

[6] J. M. Kosterlitz, The critical properties of the two-dimensional XY model, J. Phys. C 7, 1046 (1974).

[7] J. M. Kosterlitz and D. J. Thouless, Ordering, metastability and phase transitions in two-dimensional systems, Journal of Physics C: Solid State Physics 6, 1181 (1973).

[8] R. Gupta, J. DeLapp, G. G. Batrouni, G. C. Fox, C. F. Baillie, and J. Apostolakis, Phase transition in the 2D XY model, Physical Review letters 61, 1996 (1988).

[9] R. J. Baxter, Potts model at the critical temperature, Journal of Physics C: Solid State Physics 6, L445 (1973).

[10] F. Y. Wu, The potts model, Rev. Mod. Phys. 54, 235 (1982).

[11] A. C. D. van Enter, C. Külske, and A. A. Opoku, Discrete approximations to vector spin models, J. Phys. A: Math. Theor. 44, 475002 (2011).

[12] N. Butt, X.-Y. Jin, J. C. Osborn, and Z. H. Saleem, Moving from continuous to discrete symmetry in the 2D XY model, Phys. Rev. D 108, 074511 (2023).

[13] J. Chen, H.-J. Liao, H.-D. Xie, X.-J. Han, R.-Z. Huang, S. Cheng, Z.-C. Wei, Z.-Y. Xie, and T. Xiang, Phase transition of the q-state clock model: Duality and tensor renormalization, Chinese Physics Letters 34, 050503 (2017).

[14] O. A. Negrete, P. Vargas, P. F. J. G. Saravia, and E. E. Vogel, Short-range Berezinskii-Kosterlitz-Thouless phase characterization for the q-state clock model, Entropy 23, 1019 (2021).

[15] Z.-Q. Li, L.-P. Yang, Z. Y. Xie, H.-H. Tu, H.-J. Liao, and T. Xiang, Critical properties of the two-dimensional q-state clock model, Phys. Rev. E 101, 060105 (2020).

[16] H. Chen, P. Hou, S. Fang, and Y. Deng, Monte Carlo study of duality and the Berezinskii-Kosterlitz-Thouless phase transitions of the two-dimensional q-state clock model in flow representations, Phys. Rev. E 106, 024106 (2022).

[17] S. Elitzur, R. B. Pearson, and J. Shigemitsu, Phase structure of discrete abelian spin and gauge systems, Phys. Rev. D 19, 3698 (1979).

[18] J. Tobochnik, Properties of the q-state clock model for q = 4, 5, and 6, Phys. Rev. B 26, 6201 (1982).

[19] S. Hong and D.-H. Kim, Logarithmic finite-size scaling correction to the leading fisher zeros in the p-state clock model: A higher-order tensor renormalization group study, Phys. Rev. E 101, 012124 (2020).

[20] G. Ortiz, E. Cobanera, and Z. Z. Nussinov, Dualities and the phase diagram of the p-clock model, Nuclear Physics B 854, 780 (2012).

[21] T. Surungan, S. Masuda, Y. Komura, and Y. Okabe, Berezinskii-kosterlitz-thouless transition on regular and villain types of q-state clock models, J. Phys. A: Math. Theor. 52, 275002 (2019).

[22] Y. Kumano, K. Hukushima, Y. Tomita, and M. Os-hikawa, Response to a twist in systems with $z_p$ symmetry: The two-dimensional p-state clock model, Phys. Rev. B 88, 104427 (2013).

[23] C. M. Lapilli, P. Pfeifer, and C. Wexler, Universality away from critical points in two-dimensional phase transitions, Phys. Rev. Lett. 96, 140603 (2006).

[24] S. K. Baek, H. Mäkelä, P. Minnhagen, and B. J. Kim, Residual discrete symmetry of the five-state clock model, Phys. Rev. E 88, 012125 (2013).

[25] C.-O. Hwang, Six-state clock model on the square lattice: Fisher zero approach with wang-landau sampling, Phys. Rev. E 80, 042103 (2009).

[26] F. C. Alcaraz and R. Koberle, Duality and the phases of z(n) spin systems, J. Phys. A: Math. Gen. 13, L153 (1980).

[27] P. Kumar, Theory of a higher-order phase transition: The superconducting transition in $Ba_{0.6}K_{0.4}BiO_3$, Phys. Rev. B 68, 064505 (2003).

[28] C. E. Ekuma, G. C. Asomba, and C. M. I. Okoye, Ginzburg-Landau theory for higher order phase transition, Physica C 472, 1 (2012).

[29] J. Fröhlich and T. Spencer, The Kosterlitz-Thouless transition in two-dimensional abelian spin systems and the coulomb gas, Commun. Math. Phys. 81, 527 (1981).

[30] Y.-D. Hsieh, Y.-J. Kao, and A. W. Sandvik, Finite-size scaling method for the Berezinskii-Kosterlitz-Thouless transition, J. Stat. Mech.: Theory Exp., P09001.