Clock model interpolation and symmetry breaking in O(2) models

Motivated by recent attempts to quantum simulate lattice models with continuous Abelian symmetries using discrete approximations, we define an extended-O(2) model by adding a $\gamma \cos(q\varphi)$ term to the ordinary O(2) model with angular values restricted to a $2\pi$ interval. In the $\gamma \rightarrow \infty$ limit, the model becomes an extended $q$-state clock model that reduces to the ordinary $q$-state clock model when $q$ is an integer and otherwise is a continuation of the clock model for noninteger $q$. By shifting the $2\pi$ integration interval, the number of angles selected can change discontinuously and two cases need to be considered. What we call case $1$ has one more angle than what we call case $2$. We investigate this class of clock models in two space-time dimensions using Monte Carlo and tensor renormalization group methods. Both the specific heat and the magnetic susceptibility show a double-peak structure for fractional $q$. In case $1$, the small-$\beta$ peak is associated with a crossover, and the large-$\beta$ peak is associated with an Ising critical point, while both peaks are crossovers in case $2$. When $q$ is close to an integer by an amount $\Delta q$ and the system is close to the small-$\beta$ Berezinskii-Kosterlitz-Thouless transition, the system has a magnetic susceptibility that scales as $\sim 1 / (\Delta q)^{1 - 1/\delta'}$ with $\delta'$ estimates consistent with the magnetic critical exponent $\delta = 15$. The crossover peak and the Ising critical point move to Berezinskii-Kosterlitz-Thouless transition points with the same power-law scaling. A phase diagram for this model in the $(\beta, q)$ plane is sketched. These results are possibly relevant for configurable Rydberg-atom arrays where the interpolations among phases with discrete symmetries can be achieved by varying continuously the distances among atoms and the detuning frequency.


I. INTRODUCTION
In recent years, the idea of using quantum computers or quantum simulation experiments to approach the real-time evolution or the finite-density behavior of lattice models of interest for high-energy physics has gained considerable interest [1][2][3][4][5][6][7][8][9][10]. As the current noisy intermediate scale quantum (NISQ) devices that are available to implement this research program have a very limited number of quantum computing units, such as qubits, trapped ions or Rydberg atoms, it is essential to optimize the discretization procedure. Starting from the standard Lagrangian formulation of lattice field theory models with continuous field variables, one can either discretize the field variables [11][12][13] used in the path integral, expand the Boltzmann weights using character expansions [3,14,15], or use the quantum link method [16,17].
Models with continuous Abelian symmetries are of great physical interest. Besides the electromagnetic interactions of charged particles in 3+1 dimensions, this also includes models where a mass gap is dynamically generated [18,19] or a Berezinskii-Kosterlitz-Thouless (BKT) transition [20][21][22] occurs. For models with a U (1) symmetry, the character expansion mentioned above is simply the Fourier series. It has been shown [23,24] * jin-zhang@uiowa.edu that the truncation of these series preserves the original symmetry. On the other hand, the Z q clock approximation of the integration over the circle only preserves the Z q discrete subgroup. A recent proposal applies the Z q clock approximation to the simulation of the Abelian gauge theory in 2 + 1 dimensions, where transformations between the electric representation and the magnetic representation can significantly reduce the required computational resources [25]. In order to decide how good the Z q approximation is in a variety of situations, it is useful to build a continuous family of models interpolating among the various possibilities.
In this article, we focus on the case of the O(2) nonlinear sigma model in 1+1 dimensions. This model was key to understanding the BKT transition [20][21][22]26] and the corresponding Z q clock model has been studied extensively . We propose to interpolate among these models by starting with the standard O(2) action and introducing a symmetry-breaking term, ∆S(γ, q) = −γ x cos(qϕ x ). (1) When q is an integer, if we take the limit γ → ∞, we recover the Z q clock model. For the rest of the discussion, it is important to realize that the O(2)-symmetric action is 2π-periodic for all the ϕ x variables. In contrast, ∆S has a 2π/q periodicity. When q is an integer, if we apply the shift q times we obtain the periodicity of the O(2) action. In order to interpolate among the clock models, we will consider noninteger values of q while keeping a fixed ϕ interval of length 2π. The model and the effect of the symmetry breaking are discussed in Sec. II both in the standard Lagrangian and tensor formulations. The idea of having a doubly continuous set of models is interesting from a theoretical point of view but also from a quantum simulation point of view. If we attempt to quantum simulate these models using Rydberg atoms as in Refs. [48][49][50], it is possible to tune the ratio R b /a of the radius for the Rydberg blockade and the lattice spacing, as well as local chemical potentials continuously. This allowed interpolations among Z q phases for different integer values of q [50]. Sequences of clock models also appear in models for nuclear matter when the number of colors is varied [51].
It is often a difficult task to detect BKT transitions in the quantum Hamiltonian approach, as it is hard to find a good indicator of BKT transitions that has a clear discontinuity, peak or dip. The equivalence of the path integral formulation and statistical mechanics can be used to access universal features and detect phase transitions in statistical mechanics based on Monte Carlo (MC) simulations and tensor renormalization group (TRG) calculations. The Markov chain MC (MCMC) simulations efficiently explore the typical set of the physical configurations. MC calculations use the universal jump in the helicity modulus [52] as an indicator for BKT transitions. But ambiguities in the definition of the helicity modulus in the Z 5 clock model can result in controversial conclusions [39,53]. The TRG [54,55] calculations provide a coarse-grained theory where the size of the lattice spacing doubles at each step. If the truncations performed are under control, one can go to the thermodynamic limit quickly. Calculation of the magnetic susceptibility in the presence of a weak external field is a universal method to detect critical points that can be easily implemented in the TRG. However, it does not show a peak to indicate the large-β BKT transition in the five-state clock model [44,56]. The study of the γ → ∞ limit with fractional q not only provides us a clear picture of what phases the symmetry-breaking term will drive the XY model to, paving the way to discussions for the full phase diagram at finite γ, but also brings us a new tool to detect BKT transitions in Z n models. In contrast, the calculation of the specific heat at increasing volume allows us to discriminate between a second-order phase transitionwhere it diverges logarithmically with the volume in the Ising case-and a BKT transition or a crossover. This paper is organized as follows. Section II introduces the definition of the extended-O(2) model, the extended q-state clock model and thermodynamic quantities. The MC and TRG methods are introduced in Sec. III. The MC method is used to validate the TRG at small volume. The symmetry breaking is discussed in tensorial language. We discuss the behaviors of thermodynamic quantities and point out the crossover peak and the Ising peak in both the specific heat and the magnetic susceptibility in Sec. IV A. We analyze the crossover peak and the Ising critical point in Secs. IV B and IV C respec-tively. We change the integration interval and discuss a new case in Sec. IV D. The phase diagram in the (β, q) plane is sketched in Sec. IV E. We summarize our results and give outlooks in Sec. V.

II. THE MODEL
To define the models that we consider in the following we start with the two-dimensional classical O(2) nonlinear sigma model, or XY model, where the spin degrees of freedom σ are unit vectors whose possible directions are confined to a plane. They reside on the sites of a two-dimensional lattice of volume V = N x × N t (we prefer to label the second dimension as t to maintain the connection between two-dimensional classical models and 1+1-dimensional quantum field theories). The action is where the sum on x is over the sites of the twodimensional lattice, and on µ = 1, 2 over the directions. The field, h, is a uniform constant external magnetic field. It is convenient to parametrize the spins σ with a single angle ϕ ∈ [ϕ 0 , ϕ 0 + 2π). The action then takes the form where h = | h| and ϕ h is the direction of the external field that in the absence of other symmetry-breaking terms can be set to zero for convenience. Next, let us extend the model by introducing a term that can favor certain values of the angle: We call the model with the action (4) the "extended-O(2)" model. For integer q the limit γ → ∞ forces the spin angles to take the values ϕ (k) x = 2πk/q with k ∈ Z. Thus, while γ = 0 corresponds to the O(2) model, γ → ∞ corresponds to the q-state clock model. The action defined in Eq. (4) is also valid for noninteger q and we therefore consider Eq. (4) as our definition of the extension of the q-state clock model to noninteger q in the γ → ∞ limit. In this case the angle ϕ with k ∈ Z and some choice of domain [ϕ 0 , ϕ 0 + 2π). By varying ϕ 0 , we can obtain different sets of angles 2π/q +φ ϕ0 = −π (case 2) Figure 1. Arrows indicate the allowed spin orientations for the extended q-state clock model with the choice ϕ0 = 0 (left) and ϕ0 = −π (right). In this example, q = 5.5.
that are equivalent to either k = 0, 1, . . . , q (case 1) or k = 0, 1, . . . , q − 1 (case 2), since in the γ → ∞ limit-in the absence of an external field-the action only depends on the relative angle between nearest-neighbor sites ∆ϕ x (see Appendix A). Case 2 just has one fewer angle than case 1. As shown in Fig. 1, the angular distance between two adjacent values of ϕ x on a circle takes two values: 2π(q − q )/q < 2π/q for case 1, and 2π/q < 2π(1 + q − q )/q for case 2. These values including 0 have the largest Boltzmann weights in the partition function. The small angular distance is in case 1 and 2π/q in case 2. With the choice ϕ 0 = 0, we have case 1, while choosing ϕ 0 = −π is equivalent to case 2 (1) for odd (even) q . At noninteger q the Z q symmetry is explicitly broken since the action is not invariant under the operation k → mod (k + 1, q ). But there is still a Z 2 symmetry because the action is invariant under the operation k → q − k.
We also consider the limit γ → ∞ directly by simply restricting the values of the originally continuous angle ϕ to the values given in Eq. (5): We call the model (7) the "extended q-state" clock model for all values of q and the "fractional-q-state" clock model for fractional values of q. For integer q the extended qstate clock model reduces to the ordinary q-state clock model. Numerical results presented in later sections are from the extended q-state clock model.
The partition function is where the corresponds to ϕ0+2π ϕ0 for the continuous angles in the O(2) and extended-O(2) models and to ϕ (k) for the discrete angles in the extended q-state clock model.
With the models defined we turn to observables. The main observables that we compute to study the critical behavior are the internal energy, magnetization, and their corresponding susceptibilities. These quantities are defined in the same way for both the continuous and discrete angle cases. The internal energy is defined as where · · · denotes the ensemble average. The specific heat is In addition we consider the magnetization. The magnetization of a given spin configuration is and the ensemble average is then The magnetic susceptibility defined in a manifestly O(2)invariant way is We note that in Monte Carlo simulations at zero external field in a finite system in the absence of explicit symmetry breaking terms the definition of the spontaneous magnetization (12) gives M = 0. In such situations one often resorts to using a proxy observable [29,31,32,35,44,45] in place of M . The corresponding susceptibility is While one expects that | M | indicates the same critical behavior, in general, | M | is numerically different from M except deep in the ordered phase. Nevertheless, we expect both definitions of the magnetic susceptibility-Eqs. (13) and (15)-possess the same critical behavior, and can be relied upon to extract universal features. In the next section we detail the methods used to study the observables defined above.

III. METHODS
The allowed spin orientations in the extended q-state clock model, given by Eq. (5), are discrete, and the model can be studied using a heatbath algorithm. The heatbath algorithm is a MCMC algorithm that drives the lattice toward equilibrium configurations by choosing the new spin at each update according to the probability distribution defined by its neighboring spins. We adapted Fortran code developed by Bernd Berg for the standard Potts model [57].
Initial exploration of the extended q-state clock model was performed via MC on a 4 × 4 lattice with zero external magnetic field. For this model, the heatbath approach suffers from a slowdown that makes it difficult to study the large-β regime already on very small lattices. An alternative approach, the TRG, which does not suffer from this slowdown, was used to study the model on much larger lattices and in the thermodynamic limit. This allows us to perform finite-size scaling and characterize the phase transitions in the extended q-state clock models. The TRG results are validated by comparison with exact and Monte Carlo results on small lattices (see Appendix C). The TRG methods are not exact. Because the bond dimension, D bond , of the coarse-grained tensor increases exponentially with the renormalization-group (RG) steps, a truncation for D bond must be applied to avoid uncontrolled growth of memory needs on classical computers. For noncritical phases, the fixed-point tensor of the RG flow has small D bond , but a larger bond dimension is needed to have the correct RG flow near the critical point. It has been reported that TRG methods using D bond = 40 can locate the phase transition point with an error of order 10 −4 for the Ising model [56], and of order 10 −3 for the O(2) model [56,58] and the clock models [44,56].
To perform TRG calculations, we need to express the partition function as a contraction of a tensor network. We rewrite the weight of each link by a singular value decomposition (SVD): Then we sum over the original k indices and the partition function can be expressed in the dual space from the expansion in terms of n indices where l = n x−ŝ,s , r = n x,s , d = n x−τ ,τ , u = n x,τ for each site x, and the local rank-four tensor is defined as Notice that for integral q, the matrix U = V −1 can be chosen as U kn = exp(i2πkn/q), then if h = 0, the tensor C lrdu becomes a δ-function that gives a Z q selection rule for values of n: mod (n x−ŝ,s + n x−τ ,τ − n x,s − n x,τ , q) = 0. (19) The tensor reformulation of the expectation value of a local observable can be obtained in the same way. For example, the first component of the magnetization is equal to the expectation value of cos (ϕ) at an arbitrary site x 0 , m 1 = cos(ϕ (k) x0 ) , which can be expressed as where ,lrdu is an impure tensor residing at site x 0 , and To compute the internal energy, we need to calculate the expectation values of link interactions µ = cos(ϕ ) . Taking µ = s as an example, we perform another SVD for the target link and introduce two impure tensorsT i x0+ŝ,lrdu ,T i x0,lrdu residing at nearest neighbor sites x 0 +ŝ, x 0 , by replacing U k x 0 +ŝ ,l , G l with U i k x 0 +ŝ ,l , G i l and replacing V kx 0 ,r , G r with V i kx 0 ,r , G i r in Eq. (18), respectively. Thus the tensor reformulation of the expectation value of the link interaction is written as In the following, we use TRG and higher-order TRG (HOTRG) to contract tenser networks with impure tensors [59,60] up to a volume V = L 2 = 2 24 ×2 24 , calculate the first component of the magnetization m = (m 1 , m 2 ) and internal energy E = s + τ , and take derivatives of m and E with respect to h and β, respectively, to find the magnetic susceptibility and specific heat. 1,2 The locations and heights of the peaks of χ M and C v are obtained via a spline interpolation on datasets with ∆β = 10 −3 . The tensor contraction in HOTRG is performed with ITensors Julia Library [61].

A. Thermodynamics
In the extended q-state clock model, there is a Z q symmetry when q ∈ Z. When q / ∈ Z, this symmetry is explicitly broken. We choose ϕ 0 = 0, so the allowed spin orientations divide the unit circle into q arcs of which q −1 have measure 2π/q. The remainder has measureφ given in Eq. (6) and illustrated in Fig. 1. There remains a Z 2 symmetry and an approximate Z q symmetry.
Monte Carlo results obtained with a heatbath algorithm on a 4 × 4 lattice with zero external field are shown for 4.1 ≤ q ≤ 5.0 in Fig. 2. The four panels show the energy density and the specific heat defined in Eqs. (9) and (10) as well as the proxy magnetization and susceptibility defined in Eqs. (14) and (15). For q = 5, the energy density is zero at β = 0 because there is no linear term in the series expansion of the partition function due to the Z 5 symmetry. The nonzero energy density at β = 0 for q < 5 is consistent with the explicit Z 5 symmetry breaking. As q → 4 + , the stronger symmetry breaking results in a more negative energy density. There is a double-peak structure in the specific heat, where the large-β peak moves toward β = ∞ as q is decreased. The proxy magnetization also increases for smaller values of q with stronger symmetry breaking, and the peak of the magnetic susceptibility moves toward smaller β values. Note that the double-peak structure of this proxy magnetic susceptibility will appear at larger system sizes L ≥ 12 [31]. The true magnetic susceptibility with an external field at large volumes will show the double-peak structure (see below). More MC results and additional details are given in Appendix B.
In Fig. 3, we show the logarithm of the energy density of this model from TRG for q = 4.1, 4.3, 4.5, 4.7, and 4.9 with 0 ≤ β ≤ 16 and V = 1024 × 1024. For large enough β, we have so that the energy density converges exponentially with β. The results in Fig. 3 confirm this behavior. We also notice that for smaller q, there is a larger range of β where the energy density does not change much. In this range, 1 − cosφ is close to zero for q close to 4 from above, and the inverse temperature β is large enough that terms containing larger angular distances are negligible, but β is still not large enough to change the values of exp[−β(1 − cosφ)] significantly from 1 and ignore higher orders. The result for q = 4.1 shown in Fig. 3 indicates that the specific heat is almost zero for β > 2.5, which is confirmed in Fig. 4.
In Fig. 4, we show the specific heat for q = 4.1, 4.5, 4.9, and 5.0 at volumes ranging from 4 × 4 to 128 × 128. For generic q, there are two peaks in the specific heat. 3 In Fig. 4 we see only a single peak for q = 4.1 since the second peak is at much larger β. For q / ∈ Z and not too close to 5, the first peak shows little or no dependence on volume. The second peak grows logarithmically with volume, as shown in the insets for q = 4.5, 4.9. This is in contrast to the integer case q = 5 where there are two BKT transitions and both peaks show little dependence on volume for lattice sizes larger than 32 × 32. Because the specific heat is the second-order derivative of free energy, the results in Fig. 4 indicate that the first peak is associated with either a crossover or a phase transition with an order larger than 2, and the second peak is associated with a second-order phase transition. To conclusively characterize the phase transitions, if any, associated with these two peaks in the fractional-q-state clock model, we study the magnetic susceptibility in the next two subsections.
We find that the thermodynamic curves vary smoothly for n < q ≤ n + 1 where n is an integer. When q is taken slightly larger than n from below, these curves change abruptly since an additional degree of freedom is introduced. The specific heat exhibits a double-peak structure with the second peak at very large β. As q is increased further, this second peak moves toward small β, until at q = n + 1, the thermodynamic curves of the integer-(n + 1)-state clock model are recovered.
In the small-β (high temperature) regime, all allowed angles are nearly equally accessible, and the model behaves approximately like a q -state clock model. The model is dominated by the approximate Z q symmetry, and there is a peak in the specific heat. In Fig. 5, we show that at intermediate β, an explosion of the integrated autocorrelation time of the energy is observed in the MC simulation as the model quickly reduces to a rescaled Ising model. At large beta, the configuration space separates into thermodynamically distinct sectors, and the Markov chain has trouble adequately sampling both sectors. This is discussed further in Appendix B. At large-β, spin flips across the small angular distancẽ φ are strongly favored relative to spin flips across the other distances. Thus, in the large-β regime, the model behaves as a rescaled Ising model. The existence of an Ising critical point is conclusively established via TRG in Sec. IV C.
We next present our TRG results and discuss the phase transitions in the fractional-q-state clock model in the rest of this section. We first present the results for the magnetic susceptibility without an external field at small volumes in Fig. 6. For q < 5, there is a small-β peak converging quickly with volume, which means the peak is associated with a crossover. As q is increased, there is  a high plateau moving toward small β. The height of the plateau increases with volume as a power law (notice the logarithmic scale in the y-axis). The divergent plateau signals a phase transition. As there is no spontaneous symmetry breaking in small volumes, the response of the system to an external field remains high at low temperatures, so we cannot see the transition from a critical phase to the symmetry-breaking phase where there is a small magnetic susceptibility based on these results for h = 0 and small volumes. But notice that for fractional q like q = 4.9, there is a higher plateau after the first one for each volume. This is due to an "approximate Z q symmetry breaking" after a Z 2 symmetry breaking. This "approximate Z q symmetry breaking" is not a true phase transition so it moves to β = ∞ in the thermodynamic limit as shown in the results for q = 4.9. There is only one Z 5 symmetry breaking for q = 5 so there is a single plateau for each volume. We will use the magnetic susceptibility with a weak external field in the thermodynamic limit to detect the phase transitions.
In Fig. 7 we present the magnetic susceptibility in the thermodynamic limit as a function of β for q = 4.995, 4.999, 4.9999 with a small external field h = 4 × 10 −5 , 2 × 10 −5 . The height of the large-β peak is large for all three values of q, indicating a phase transition near this peak. The small-β peak for q = 4.995, h = 4 × 10 −5 is invisible because it is very small. For q = 4.999, h = 4 × 10 −5 , closer to 5, the small-β peak is higher, and the large-β peak moves toward the large-β BKT transition point β BKT q=5,c2 for q = 5. When the external field is decreased to h = 2×10 −5 , the large-β peak height is almost doubled, while the small-β peak height does not change, which means there is no phase transition near the smallβ peak. For q = 4.9999 (closer to 5) with h = 4 × 10 −5 , one can see that the small-β peak becomes higher than the large-β one, and the large-β peak is fading away, which is consistent with the results in Refs. [44,56] for the five-state clock model. When the external field is de-  creased to h = 2 × 10 −5 , the large-β peak grows a much larger amount of height than the small-β one does and becomes higher than the small-β one. We have confirmed that the small-β peak will eventually converge for small enough h. All these behaviors are evidence that for all fractional q > 4, the small-β peak in χ M does not diverge, and the large-β peak diverges at h = 0 indicating a phase transition.
In the following, we discuss the behavior of the two peaks in the thermodynamic limit in detail. The main observations are that for ϕ 0 = 0 and q > 4, there are two peaks in the specific heat and the magnetic susceptibility, the small-β one is finite and is associated with a crossover, and the large-β one diverges which is characteristic of an Ising critical point. When q is approaching an integer from below, the height of the crossover peak of the magnetic susceptibility diverges as q − q goes to zero with a power law: We can formulate the scaling hypothesis with ∆q = q − q and h, where λ parametrizes a scale transformation, p and r are the scaling dimensions, and d = 2 in two-dimensional space. Notice that the reduced temperature should not enter the homogeneous function independently because there is an essential singularity in the correlation length as a function of temperature for BKT transitions. We assume Eq. (26) holds for any critical temperature, in particular, the BKT crossover peak position β BKT (∆q, h), which is a power-law function of ∆q and h, considered in the following calculations. Then the magnetization and  the magnetic susceptibility satisfy the following relations: from which one can obtain where 1/δ = (d − r)/r, 1/δ = (d − r)/p. We then have δ = (δ − 1)/y = 14/y, where we have used the fact that δ = 15 for BKT transitions. Note that the expansion of the action to the first order in ∆q is where S Z q is the action for the integer-q clock model. The perturbation term that breaks the Z q symmetry has a very different form from the one for the external field h cos(ϕ x ). We numerically show in the following that δ is equal to the magnetic critical exponent δ = 15 for the BKT transitions and Ising critical points in two dimensions. Thus if we define a new susceptibility as ∂M/∂∆q, which should scale as (∆q) −1+1/δ , the exponent y = 1 − 1/δ = 14/15, still the same as y. However, calculating δ from y requires higher accuracy since δ = 1/(1−y ) where one significant digit is subtracted in the denominator. Both the peak position of the crossover and the Ising critical point go to the two BKT transition points at integer q with the same power law, which provides us a new way to locate the phase transitions in models with Z n symmetry. However, for ϕ 0 = −π and odd q , both peaks of χ M are finite for fractional q so there are no critical points. For q → 5 + , only the small-β peak can be used to extract the BKT transition of q = 5 since the large-β peak fades away.

B. Small-β peak: Crossover
We have shown in Fig. 7 that for q = 4.999, the height of the small-β peak converges to about 400 for small enough external field. The dependence of the peak height   The log-log plot of the maximal value of the magnetic susceptibility χM for the small-β peak as a function of 5 − q. The peak height diverges with a power law when q → 5 − . on the external field is larger for values of q closer to an integer from below. In Fig. 8, we show another example for q = 4.3. One can see that the h dependence of the peak height is much smaller. The peak height at h = 10 −2 differs from the value in the h → 0 limit by only about 0.04. As the external field is decreased, the peak height converges to a constant χ M ≈ 1.2, implying that there is no phase transition around this peak. This is true for all fractional q. Thus, for fractional q, the first peak in the specific heat is associated with a crossover rather than a true phase transition. As q approaches 5 from below, we expect the small-β peak height to diverge because there is a BKT transition for q = 5, and we expect the location of the peak to go to the small-β BKT transition point.
To check this, we calculate the converged peak height  Figure 10. Same as Fig. 9, but for q → 6 − . dh = 10 −10 in the numerical differentiation. We plot χ * M versus 5 − q in Fig. 9, where we see that the dependence of the peak height on D bond is very small for D bond ≥ 40. Applying a linear fit to ln(χ * M ) versus ln(5−q) shows that the peak height diverges as q → 5 − with a power law from which we obtain δ = 14.77(9) for D bond = 40 and δ = 14.911(10) for D bond = 50. The value of δ is close to the magnetic critical exponent δ = 15 for BKT transitions and Ising critical points in two dimensions, but now 5 − q (rather than the external field h) is playing the role of the symmetry-breaking parameter. Since both h and 5 − q break the Z 5 symmetry to a Z 2 symmetry, the agreement on the critical exponents is reasonable. The value of δ is also checked for q → 6 − in Fig. 10, where a larger D bond is applied in HOTRG. In this case, the linear fit gives The location of the converged peak β p as a function of 5 − q is depicted in Fig. 11. The discrepancy of the peak positions between D bond = 40 and D bond = 50 is invisible at 5 − q = 10 −4 and increases as 5 − q is decreased, which is reasonable because a larger bond dimension is needed for systems closer to a critical point. The overall discrepancy is small. We extrapolate the peak position to q = 5 with a power law and obtain the small-β BKT transition point β BKT q=5,c1 = 1.0494(6) with D bond = 40 and β BKT q=5,c1 = 1.0506(4) with D bond = 50 for the five-state clock model. The results are consistent with 1.0503(2) in Ref. [47]. The same procedure is performed for q → 6 −  in Fig. 12. The peak positions β p for D bond = 50 and those for D bond = 60 have little discrepancy for 6 − q ≥ 7 × 10 −7 . For 6 − q < 7 × 10 −7 , β p for D bond = 60 becomes larger than that for D bond = 50. The powerlaw fit gives us β BKT q=6,c1 = 1.0983(4) with D bond = 50 and β BKT q=6,c1 = 1.1019(5) for D bond = 60. The results are consistent with 1.101(4) in Ref. [62].
We have shown that q − q plays the same role as an external magnetic field for the small-β BKT transition in integer-q-state clock models. The magnetic susceptibility is always finite for fractional q < q , and diverges as q → q − with a critical exponent y = 14/15. This provides us an alternative way to extract the locations of BKT transitions in clock models. However, the situation is very different for the large-β peak of the magnetic susceptibility.

C. Large-β peak: Ising criticality
To understand the large-β peak in the specific heat, we again study the magnetic susceptibility χ M . For a fixed q, the critical point, if any, is given by the location of the peak of χ M in the limit h → 0 where the height χ * M of the peak is infinite. A power-law extrapolation to h = 0 is performed on peak positions of χ M for small values of h. In Fig. 13, we present the peak height χ * M of the susceptibility as a function of the external field h. The linear fit of ln(χ * M ) versus ln(h) gives χ * M ∼ h −0.93318(16) , from which we obtain the magnetic critical exponent δ = 14.97(4). This value is consistent with the value δ = 15 of the BKT transitions and Ising critical points. We have shown before that there is no phase transition around the small-β peak of χ M . A BKT transition should be accompanied with a continuous critical region, so the divergent large-β peak of χ M must be an Ising critical point. In the Ising universality class, where ν e = 1, β e = 1/8, γ e = 7/4 are the universal critical exponents. There is a universal function relating χ M /L 1.75 and L(β − β c ) with fixed hL 15/8 . In Fig. 15, we plot χ M /L 1.75 versus L(β − β c ) for various lattice sizes around the large-β peak of χ M , where the value of β c is obtained from the Ising approximation in Eq. (37) described below. One can see that all the data collapse onto a single curve, which gives strong evidence that this is a critical point of the Ising universality class. In Fig. 14, we show the logarithm of χ M as a function of log 2 (L) for β = 1.1. For q = 5, β = 1.1 is between two BKT points so is a critical point, and χ M keeps increasing with a positive power of L as expected. When q < 5, even  for a very small 5 − q = 10 −6 , χ M always saturates at a large enough volume. These results again prove that there are no BKT transitions in fractional q, no matter how close to q the q is. Because the maximal value of χ M should diverge as a negative power of 5 − q when q → 5 − , the increment of the height of the plateaus between 5 − q = 10 −n and 5 − q = 10 −n−1 should be a constant for any integer n, which is also confirmed in Fig. 14.
We next obtain the Ising critical point by extrapolating the peak position of χ M to h = 0. An example for q = 4.9 is shown in Fig. 16. The power-law extrapolation gives β c = 1.44614 (2). We can repeat the same procedure for other values of q. But notice that we need a larger bond dimension in TRG when q is very close to an integer from below, because more degrees of freedom become important, and the critical point is close to a BKT transition point. The phase transition in the Ising universality class is a transition from a disordered phase to a symmetrybreaking phase. The structure of the fixed-point tensor in TRG can easily characterize this phase transition. As proposed in Ref. [63], the symmetry-breaking indicator should be 1 in the disordered phase and 2 in the Z 2 symmetry breaking phase. Thus the discontinuity in X for the fixed point tensor as a function of β can be used to locate the phase transition. An example for q = 4.9 is shown in Fig. 17, where the value of X changes from 1 to 2 at 1.4461 < β < 1.4462, consistent with the result from the extrapolation of the peak position of χ M in Fig. 16. The advantage of this method is that we only need to contract a single tensor network for each value of β and scan a β range once to locate the phase transition point. This saves us a lot of computational effort when we extrapolate the Ising critical point to q = q from below. Now we use the value of X to locate the Ising critical point and find the large-β BKT transitions for q = 5, 6. Notice that for q = 5, although both a small external field and a small deviation of q from an integer break the Z q symmetry to a Z 2 symmetry, the magnetic susceptibility with a weak external field does not have a peak around the large-β BKT transition, so it fails to predict the location of the phase transition [44,56], but here we always have an Ising critical point for fractional q. In Fig. 18(a), we calculate the Ising critical point for q → 5 − with D bond = 40 and extrapolate the result to q = 5 with a power law. The value of 5 − q is between 6 × 10 −4 and 10 −3 , where the dependence of β c on D bond is small. The extrapolated BKT transition point for q = 5 is β BKT q=5,c2 = 1.1027 (14), consistent with the result 1.1039(2) obtained in Ref. [47]. As a comparison, we also present the extrapolation of the crossover peak with the same D bond in the same figure. The exponents of the two power-law scalings are the same within uncer-  Figure 17. The β dependence of X from the fixed point tensor for q = 4.9. The discontinuity is located between 1.4461 and 1.4462, consistent with the result from the extrapolation of the peak position of χM to zero external field in Fig. 16. The tensorial bond dimension is 40. tainties, and the values of the exponents are consistent with 0.2677(84) obtained in Ref. [44] for magnetic susceptibility with an external field h ≤ 10 −3 .
The results for q → 6 − are shown in Fig. 18(b), where we use a larger bond dimension D bond = 60 and data with 6 − q ≤ 10 −4 . The extrapolated BKT transition point is β BKT q=6,c2 = 1.435 (3), consistent with 1.441(6) in Ref. [62]. Comparing the Ising critical points and the crossover peak positions, we find that the power-law scaling parts are exactly the same except for a minus sign within uncertainties, which means they approach the two BKT transition points in the same manner. We believe this behavior can be seen for all q → n − for integer n ≥ 5.
To determine the exponent accurately, we need to use a larger D bond and data closer to an integer, which is beyond the scope of this work. However, this exponent should be the same as the power-law scaling of β p of χ M in a weak external field for all clock models with integer q ≥ 5, where there is always an emergent O(2) symmetry [42,64]. In the limit of O(2) model, this exponent is found to be around 0.162 [58,65,66].
At large β, the fractional-q-state clock model is a rescaled Ising (q = 2) model because the link interactions for the two angles 0 andφ dominate the weights in the partition function. There are two peaks in the specific heat, and if these peaks are sufficiently separated, then the second peak is that of an Ising model where β is rescaled as β → αβ, with a rescaling factor where the small angular distanceφ in the model depends on q and is defined in Eq. rescaled Ising model, In Table I we list some of these critical points for different values of q. These points give approximately the location of the large-β peak in the specific heat for the extended q clock model. The critical point β rIsing of the rescaled Ising model (see Table I) is a good approximation for the true critical point β c for values of q that are not too close to an integer from below. For example, in Fig. 19 we compare the specific heat from TRG with the specific heat of the rescaled Ising model for q = 4.5 and several different lattice sizes. For lattices 16 × 16 and larger, the specific heat for the rescaled Ising model accurately captures the large-β peak of the fractional-q-state clock model. As q approaches an integer from below, the approximation begins to fail. In Fig. 20, we compare the rescaled Ising critical points with the true critical points as q → 5 from below. Most of the true critical points are obtained from X defined in Eq. (35), and five points are from χ M and agree perfectly with those from X. For 4.0 < q 4.7, the difference between the true critical point and the Ising approximation is less than 0.01. As q → 5 − , the difference becomes larger and is around 0.175 for q = 5.

D. Integration interval
For the extended q-state clock model, the allowed angles are restricted to the integration interval ϕ ∈ [ϕ 0 , ϕ 0 + 2π). All results presented so far used ϕ 0 = 0, which is in the so-called case 1. We also considered the possibility ϕ 0 = −π. As shown in Fig. 21, for q = 4.5, the number of allowed spin orientations and their relative sizes are the same. In fact, the models with ϕ 0 = 0 and ϕ 0 = −π are equivalent and in case 1 for all even q . For odd q , the choice of ϕ 0 = 0 or ϕ 0 = −π results in a different number  of allowed spin orientations and different thermodynamic behaviors. For example, when 5 < q < 6, there are six allowed orientations when ϕ 0 = 0, but only five allowed orientations when ϕ 0 = −π (see Fig. 21 for q = 5.5). One can show that the model with ϕ 0 = −π is in case 2 for all odd q . For integer values of q, the extended q-state clock model reduces to the ordinary q-state clock model for both ϕ 0 = 0 and ϕ 0 = −π. We consider 5 < q < 6 and ϕ 0 = −π. In Fig. 22, we show the magnetic susceptibility at h = 0 as a function of β. First of all, one can see that the magnetic susceptibility is finite for all values of q presented here, which means there are no phase transitions for any 5 < q < 6 and ϕ 0 = −π. The maximal values of χ M are larger for q closer to 5. For q = 5.1, we see a double-peak structure q = 4.5  The susceptibility does not diverge with the volume, which implies there is no phase transition here. This is different from the model with ϕ0 = 0, which has a divergent peak in the susceptibility corresponding to a second-order phase transition. in χ M , with the large-β peak higher than the small-β peak. As q increases toward 6, the first crossover peak fades away and disappears around q = 5.3. The magnetic susceptibility eventually converges to a single-peak structure with the peak position around β = 1.5 and the peak height around 4. In order to see the behavior of χ M for q → 5 + , we plot χ M versus β for q = 5.001 and q = 5.00001 in Fig. 23. We see that the maximal value of χ M is much larger than those in Fig. 22, because we are approaching a Z 5 clock model from above. Unlike q = 5.1, the small-β peak becomes higher than the largeβ peak for q = 5.001, and the large-β peak fades away as we move closer to q = 5. The small-β peak moves towards the small-β BKT transition for q = 5 so it can be used to extrapolate the value of β BKT q=5,c1 , while the large-β peak moves across the large-β BKT transition for q = 5 from right to left so it fails to predict the value of β BKT q=5,c2 . This means that the magnetic susceptibility cannot capture the crossover going to the large-β BKT transition point for q → 5 + . But the crossover should go to the BKT transition point as the Z q symmetry is recovered. The cross derivative of the free energy should be able to capture this [56], but we will just focus on the magnetic susceptibility.
We then check the divergent behavior of the small-β peak of χ M as q → 5 + using D bond = 50. Figure 24 shows the linear fit for ln(χ * M ) versus ln(q − 5). We see that the peak height of χ M diverges as a power law χ * M ∼ 1/(q−5) 0.9275 (15) , which again gives a value of the critical exponent δ = 15.09(2), close to the expected value 15. We extrapolate the peak position to q = 5 from above in Fig. 25. One can see that the power-law scaling is the same as case 1 where q → 5 − , ϕ 0 = 0, and gives β BKT q=5,c1 = 1.0514(5), consistent with the result in Fig. 11.
phase and a Z q symmetry-breaking phase separated by a second-order phase transition. For q ≥ 5, there is a disordered phase at small-β and a Z q symmetry-breaking phase at large-β with a critical phase for intermediate β between them. The boundaries of the critical phase are two BKT transition points of infinite order [42]. In our extended q-state clock model, one must make a choice of the integration interval ϕ ∈ [ϕ 0 , ϕ 0 + 2π). For the choice ϕ 0 = 0 and fractional q, both the specific heat and the magnetic susceptibility have a double-peak structure. We have shown that the small-β peak is associated with a crossover, and the large-β peak is associated with a phase transition of the Ising universality class. For the choice ϕ 0 = −π and fractional q, the phase structure is a little more complicated. For even q , we get the same behavior as with ϕ 0 = 0, but for odd q , we get a trivial case with no critical point.
The phase diagrams for both ϕ 0 = 0 and ϕ 0 = −π are shown in Fig. 26. For ϕ 0 = 0, as q → q − , the Z q symmetry is recovered. For q > 4 and ϕ 0 = 0, both the small-β crossover line and the large-β Ising critical point are smoothly connected to the small-β BKT transition point and the large-β BKT transition point for integer q from below, respectively. Notice that for q < 4 and ϕ 0 = 0, only the large-β Ising critical point is smoothly connected to the second-order phase transition for integer q from below, while the crossover peak fades away for q close enough (around 3.9 for 3 < q < 4) to an integer from below. When q → q + and ϕ 0 = 0, the Ising critical point goes to infinity, while the crossover line goes to a smaller value than the phase transition point for the q -state clock model because there is 1 more degree of freedom than the Z q clock model. For even q and ϕ 0 = −π, the phase diagram is the same as that for even q and ϕ 0 = 0. For odd q and ϕ 0 = −π, the Z q symmetry is recovered when q → q + . For 5 < q < 6, both crossover lines are smoothly connected to the two BKT transition points for integer q from above. When q is increased toward an even q and ϕ 0 = −π, the smallβ crossover line fades away and the large-β crossover line goes to a larger value than the large-β BKT transition for integer q because there is 1 fewer degree of freedom than the q -state clock model. For 3 < q < 4, there is only one crossover line that is smoothly connected to the second-order phase transition point for q = 3 and goes to around 0.77 when q → 4 − .

V. SUMMARY AND OUTLOOK
Interpolations among Z n clock models have been realized experimentally using a simple Rydberg simulator, where Z n (n ≥ 2) symmetries emerge by tuning continuous parameters, the detuning and Rabi frequency of the laser coupling, and the interaction strength between Rydberg atoms [50]. This paves the way to quantum simulation of lattice field theory with discretized field variables. We are interested in a theory that can interpolate among the O(2) model and Z n clock models. We define an extended-O(2) model by adding a symmetry breaking term γ cos(qϕ x ) to the action of the two-dimensional O(2) model. For integer q, Z q clock models emerge for large enough γ. For fractional q, we believe there exists a much more interesting phase structure. The first step to graph the full phase diagram in the (γ, q, β) cube is to consider the limit γ → ∞. In this work, we studied the fractional-q-state clock model as the γ → ∞ limit of the extended-O(2) model with angular variables in the domain [ϕ 0 , ϕ 0 + 2π). In this limit, the angular variable takes discrete values ϕ x,k = 2πk/q with k integral. By varying ϕ 0 , the set of values of integer k take either case 1 (0, 1, . . . , q ) or case 2 (0, 1, . . . , q − 1). In case 1, Z q symmetry is recovered as q → q − , while Z q symmetry is recovered as q → q + in case 2.
For the integer-q-state clock model, there is a single second-order phase transition when q = 2, 3, 4. When q ≥ 5, there are two BKT transitions with a critical phase between them. We studied the fractional-q-state clock model using Monte Carlo (MC) and tensor renormaliza-  [65], and a critical phase at large β. In the γ = ∞ plane, it is the extended q-state clock model, which has the phase diagram shown in Fig. 26. In this example we have ϕ0 = 0. Establishing the phase diagram at finite-γ will be addressed in future work.
tion group (TRG) methods. We establish the phase diagram of the model for both ϕ 0 = 0 and ϕ 0 = −π. When ϕ 0 = 0, we are in case 1, and analysis of the finite-size scaling shows a crossover and a phase transition of the Ising universality class. When ϕ 0 = −π, we are in case 1 for even q and in case 2 for odd q . There are no critical points for case 2.
In case 1, we found that there are two peaks in both the specific heat and the magnetic susceptibility. The height of the small-β peak is always finite for fractional q. The large-β peak diverges and characterizes an Ising critical point. When q → q − and q < 4, the large-β Ising critical point is smoothly connected to the second-order phase transition point for Z q clock models, while the small-β peak fades away. When q → q − and q > 4, the large-β Ising critical point and the position of the small-β peak are smoothly connected, with the same power-law scaling ∼ ( q − q) b , to the large and small BKT points respectively for Z q clock models. We also found that the height of the small-β peak of the magnetic susceptibility diverges as a power law 1/( q − q) 14/15 , from which we obtain a critical exponent δ = 15 in the ansatz of the scaling of the magnetization M ∼ ( q − q) 1/δ . This critical exponent is equal to δ associated with the magnetization with an external field M ∼ h 1/δ . In case 2, there are no critical points. When q → q + and q > 5, the small-β peak also goes to the small BKT point with the same power-law scaling and the same δ exponent as case 1, while the large-β peak fades away and cannot be used to extrapolate the large BKT point of Z q clock models.
To use the magnetic susceptibility to locate a critical point, a weak external field must be applied for the magnetic susceptibility to be finite, and extrapolate the peak position to h = 0. This method works in most cases, but the peak fades away near the large-β BKT point of integer-q-state clock models. Our procedure provides an alternative approach to locate the BKT transitions of clock models, by breaking the Z q symmetry to a Z 2 symmetry in the q direction instead of h direction. This procedure creates an Ising critical point that can be used to extrapolate the large-β BKT point for clock models.
Our results clarify what phases the symmetry-breaking term γ cos(qϕ x ) will drive the system to. These phases should have boundaries in the finite-γ direction. For small enough γ, the extended-O(2) model should go back to the same universality class as the ordinary XY model, which has been studied extensively [58,65,[67][68][69][70][71][72]. For the ordinary XY model, there is a single BKT transition from a disordered phase to a quasi-long-range-ordered critical phase at β c = 1.11995(6) [65]. Figure 27 shows the work that remains to be done to figure out the phase diagram in the (β, γ, q) space interpolating between the known phase diagram at γ = 0 and the phase diagram at γ = ∞ discussed here. There should be a rich phase diagram in the finite-γ region, which is beyond the scope of this work and will be discussed in future work.
It is interesting to note that the BKT critical point found in the O(2) model-and here at the limit of the extended phase diagram-can also be reached through a completely different interpolation. By considering the O(3) nonlinear sigma model with an additional symmetry breaking term which breaks the O(3) symmetry down to an O(2) symmetry, one can interpolate between Z 2 to O(3), and from O(3) to O(2) by tuning the sign, and magnitude, of the additional symmetry-breaking term [73]. Further additional symmetry breaking terms could be interesting. Positive-definite worm algorithms have been constructed for the O(3) nonlinear sigma model, and could be used to simulate the model efficiently [74,75].
Topics currently under study include the autocorrelations at different volumes, dynamical critical exponents, spatial correlations, vortices, density of states, and zeros of the partition function. interpolation in n < q ≤ n + 1 for integer n. For any irrational q, this model without cutoff becomes equivalent to the ∞-state clock model, that is, the XY model.
In Fig. 28 we explore the effect of increasing the cutoff for the model with q = 3.141592654 ≈ π. For example, for the case ϕ ∈ [0, 2π), the allowed angles are 2πk/π = 2k with k = 0, 1, 2, 3, whereas for the case ϕ ∈ [0, 4π), we have k = 0, 1, . . . , 6. In all of these cases there remains a Z 2 symmetry, and we expect the model to have an Ising transition at very large β for certain cutoffs. As the upper limit of the domain is moved to infinity (i.e. as the cutoff is removed), this Ising transition moves to infinity and the model becomes the XY model. A minor detail is that 3.141592654 = 3141592654/1000000000 is in fact a rational number, so one would actually get the 3141592654-state clock model if the cutoff were removed completely. However, it would be indistinguishable from the XY model for practical purposes.

Appendix B: Monte Carlo Results
We used a heatbath algorithm to study the extended qstate clock model on a 4×4 lattice. The general structure of the Monte Carlo algorithm is as follows: Parameters used in the primary data production for the extended q-state clock model are given in Table II. The total number of heatbath sweeps performed (not including equilibration) is nrpt × nmeas × ndisc. The number of measurements taken is nrpt × nmeas. Ensemble averages and error bars for the energy and magnetization were calculated after binning with nrpt bins each of size nmeas. For specific heat and magnetic susceptibility, the measurements were binned and jackknifed. We studied the extended q-state clock model on a 4×4 lattice with zero external magnetic field. For each β, we initialized to a random lattice (hot start) then we performed 2 15 equilibrating sweeps followed by 2 22 measurement sweeps. Each measurement sweep was followed by 2 8 discarded sweeps. We calculated the energy density and specific heat as defined in Eqs. (9) and (10). We calculated the proxy magnetization and susceptibility as defined in Eqs. (14) and (15)  In the large-β (low temperature) regime, the heatbath algorithm has difficulty appropriately sampling the configuration space. At large-β, the lattice freezes along a particular magnetization direction. When q ∈ Z all magnetization directions are equivalent, however, when q / ∈ Z, the discrete rotational symmetry is broken and the direction of magnetization matters. The configuration space is split into two sectors with different thermodynamic properties. In one sector, the magnetization is in the direction 0 or −φ, defined in Eq. (6), where relatively large fluctuations may still be possible. In the other sector, the magnetization is in one of the other directions where fluctuations are less likely. To appropriately sample the configuration space one has to use very large statistics or run multiple heatbath streams at the same parameters but with different seeds for the random number generator. In Fig. 29, we show an example of this phenomenon for q = 4.5 and β = 2.5 on a 4 × 4 lattice. In this example, we need 2 32 heatbath sweeps to adequately sample both sectors. This Monte Carlo slowdown makes it difficult to study larger lattices, and is a strong motivation for using TRG.
The Monte Carlo slowdown is illustrated by an explosion of the integrated autocorrelation time in the intermediate-β regime. For an observable O, an estimator of the integrated autocorrelation time is given bỹ where τ O,int is nearly independent of t. The integrated autocorrelation time for 4 < q ≤ 5 is shown in the main document in Fig. 5. The values of T needed to extract these points are given in Fig. 30.
To mitigate the effect of autocorrelation in our results, we discarded 2 8 heatbath sweeps between each saved measurement. The saved measurements were then binned (i.e. preaveraged) with bin size 2 16 before calculating the means and variances.
Appendix C: Validating TRG with MC Whereas Monte Carlo methods are well understood in the context of classical spin models, TRG is a relatively new approach. We validate TRG results at small β and small volume using exact and Monte Carlo results. Exact results can be computed for q = 2 (Ising model) and q = 4 (two coupled Ising models). We use Monte Carlo to validate the TRG results for other (including fractional) values of q.
In Fig. 35, we show that the energy density from TRG agrees very well with the exact calculation, and that the tiny discrepancy appears only around the critical point. In TRG, the specific heat is calculated by taking a finite difference derivative of the energy. In Fig. 36, we compare the specific heat from TRG with the exact values for q = 4.0. TRG deviates from the exact results near the peak of the specific heat, but this deviation is mostly due to the discretization error from the derivative.
Exact solutions are not known for fractional-q, so we       Figure 36. Same as Fig. 35, but for the specific heat. TRG shows deviations from exact results near the peak, however, this is mostly due to discretization error from the derivative.
validate TRG by comparing with results from Monte Carlo at small β and small volume. For example, Fig. 37 shows that the discrepancy between the energy density from TRG and that from Monte Carlo is only of order 10 −4 . In Fig. 38, we compare the specific heat from TRG with that of Monte Carlo for q = 4.9. TRG deviates from  Figure 38. Same as Fig. 37, but for the specific heat. TRG shows deviations from Monte Carlo results near the peak, however, this is mostly due to discretization error from the derivative.
the Monte Carlo results near the peak of the specific heat. However, this deviation is again almost entirely due to discretization error from the derivative.