Critical behavior of 3D Z(N) lattice gauge theories at zero temperature

Three-dimensional $Z(N)$ lattice gauge theories at zero temperature are studied for various values of $N$. Using a modified phenomenological renormalization group, we explore the critical behavior of the generalized $Z(N)$ model for $N=2,3,4,5,6,8$. Numerical computations are used to simulate vector models for $N=2,3,4,5,6,8,13,20$ for lattices with linear extension up to $L=96$. We locate the critical points of phase transitions and establish their scaling with $N$. The values of the critical indices indicate that the models with $N>4$ belong to the universality class of the three-dimensional $XY$ model. However, the exponent $\alpha$ derived from the heat capacity is consistent with the Ising universality class. We discuss a possible resolution of this puzzle. We also demonstrate the existence of a rotationally symmetric region within the ordered phase for all $N\geq 5$ at least in the finite volume.


Introduction
Models possessing global and/or local discrete Z(N) symmetry play an important role in many branches of physics ranging from the solid state physics to the description of the universality features of the deconfining transition in SU(N) gauge theories. In this paper we are interested in the Z(N) lattice gauge theory (LGT) at zero temperature. Assigning the gauge fields s n (x) = 0, 1, · · · , N − 1 to the links of a simple hypercubic lattice, the most general action of the isotropic Z (N) LGT can be written as where e n , n = 1, 2, 3, denotes the unit vector in the n-th direction. Similarly, the most general action of the Z(N) spin model is given by In both cases we used the convention The standard Potts model corresponds to the choice when all β k are equal. Then, the sum over k reduces to a delta-function on the Z(N) group. The conventional vector model corresponds to β k = 0 for all k = 1, N − 1. For N = 2, 3 the Potts and vector models are equivalent. Two-dimensional (2D) standard and vector Z(N) LGTs are exactly solvable both in the finite volume and in the thermodynamic limit. They exhibit no phase transition at any finite value of the coupling constant β. In particular, the rectangular R × T Wilson loop in the representation k obeys the area law W k (S) = exp (−σ k (N) R T ) , thus implying a permanent confinement of static charges. For example, the string tension of the vector model in the thermodynamic limit reads Here, I k (x) is the modified Bessel function.
No exact solution has been found for any Z(N) model in 3D, where the phase structure becomes highly non-trivial. While the phase structure of the general model defined by (1) remains unknown, it is well established that Potts models and vector models with only β 1 non-vanishing have one phase transition from a confining phase to a phase with vanishing string tension [1,2,3]. Via duality, Z(N) gauge models can be exactly related to 3D Z(N) spin models. In particular, a Potts gauge theory is mapped to a Potts spin model, and such a relation allows to establish the order of the phase transition. Hence, Potts LGTs with N = 2 have second order phase transition, while for N ≥ 3 one finds a first order phase transition. Since Z (2) LGT is equivalent to the Ising model, its critical behavior is well known (see Refs. [4] and references therein). Generally, the Z(N) global symmetry of the finite-temperature 4D SU(N) gauge theory motivated thorough investigations, both analytical and numerical, of the 3D spin models, especially for N = 2, 3 [5,6] (for more recent studies, see [7] and references therein). The 3D Potts models for N > 3 have been simulated in [8] and studied by means of the high-temperature expansion in [9].
Surprisingly, much less is known about the critical behavior of Z(N) vector LGTs when N > 4. They have been studied numerically in [10] up to N = 20 on symmetric lattices with size L ∈ [4 − 16]. It was confirmed that zero-temperature models possess a single phase transition which disappears in the limit N → ∞. A scaling formula proposed in [10] shows that the critical coupling diverges like N 2 for large N. Thus, the U (1) LGT has a single confined phase in agreement with theoretical results [11]. We are not aware, however, of any detailed study of the critical behavior of the vector models with N ≥ 4 in the vicinity of this single phase transition. Slightly more is known about the critical properties of Z(N) vector spin models. In particular, it has been suggested that all vector spin models exhibit a single second order phase transition [12]. An especially detailed study was performed on the Z(6) model, because the Z(6) global symmetry appears as an effective symmetry of the Z(3) antiferromagnetic Potts model [13,14]. The computed critical indices suggest that the Z(6) vector model belongs to the universality class of the 3D XY model. An interesting feature of the Z(6) model and, possibly, of all vector models with N > 4, is the appearance of an intermediate rotationally symmetric region below the critical temperature of the second order phase transition. The mass gap was however found to be rather small, but non-vanishing in this region [13]. Combined with a renormalization group (RG) study, the analysis concluded that this intermediate region presents a crossover to a low-temperature massive phase, where the discreteness of Z(6) plays an essential role [14].
The main goal of the present work is to fill the gap in our knowledge about the critical behavior of the 3D Z(N) LGTs. Another motivation comes from our recent studies of the deconfinement transition in the Z(N) vector LGT for N > 4 at finite temperatures [15,16,17]. The major findings of these papers was the demonstration of two phase transitions of the Berezinskii-Kosterlitz-Thouless type and the existence of an intermediate massless phase. The critical indices at these transitions have been found to coincide with the indices of the 2D vector spin models. An interesting question then arises regarding the construction of the continuum limit of the finite-temperature models in the vicinity of the critical points. For this to accomplish it might be useful, and even necessary, to know the scaling of quantities such as string tension, correlation length, etc. near the critical points of the corresponding zero temperature models.
In this work we are going to: • perform an analytical study of the general 3D Z(N) LGT for various N using a "phenomenological" RG; • locate critical points of the vector models via Monte Carlo simulations of the dual of the Z(N) LGT and determine their scaling with N; • compute some critical indices and establish the universality class of the models; • illustrate the existence of a rotationally symmetric region within the ordered phase for all N ≥ 5 at least in the finite volume.
This paper is organized as follows. In Section 2 we formulate our model and recall the exact duality relation with a generalized 3D Z(N) spin model. Section 3 is devoted to a RG study of the models. Within a modified version of the phenomenological RG, we explore the space of coupling constants, find fixed points and calculate the critical index ν. In Section 4 we present the setup of Monte Carlo simulations, define the observables used in this work and present the numerical results. In particular, we locate the position of critical points and compute various critical indices at these points. As a cross-check, we also simulated Z(N = 2, 3) models for which high-precision results. The same Section deals also with the computation of the average action and the heat capacity in the vicinity of critical points and the derivation of the index α from a finite size scaling (FSS) analysis of the heat capacity. In Section 5 we discuss some findings regarding the symmetric region below the critical point. All results are finally summarized in Section 6.
The 3D Z(N) LGT on an isotropic lattice can generally be defined as The most general Z(N)-invariant Boltzmann weight is The U(1) LGT is defined as the limit N → ∞ of the above expressions.
To study the phase structure of 3D Z(N) LGTs one can map the gauge model to a generalized 3D Z(N) spin model with the action given by Eq. (2). The relation between spin β s k and gauge β g k couplings can be computed exactly and reads where the dual Boltzmann weight is defined as In what follows we are going to simulate the Z(N) vector LGT. In this case we use β 1 = β N −1 = β and the dual weight becomes An example of the explicit relations between couplings in this case together with some further important comments on the relations can be found in [17].
With the goal of performing RG transformations, it is somewhat more convenient to use a different but equivalent representation for the Boltzmann weight, namely The set of coupling constants {t k } can be chosen to satisfy The coupling constants t k and β k can be connected to each other via the Fourier transform on the Z(N) group. The dual of the partition function (5) can then be presented as Here, the summations over j n , n = 1, 2, 3, enforce the global Bianchi constraints due to the periodic BC and η n = j n on a set of links dual to any fixed closed surface wrapping the original lattice in the directions perpendicular to n, otherwise η n = 0.

Results of the RG study
As a first step we consider the general 3D Z (N) LGT with the Boltzmann weight (10) and study its phase structure with the help of the phenomenological renormalization group (PH RG) [18]. For the 3D Ising model this RG was used in [19]. As is well known, the PH RG gives in many cases not only the qualitatively correct phase diagram of a model, but also a good quantitative approximation for the observables near the critical points, and this approximation can be systematically improved. We follow the general strategy described in the original papers [18] and only briefly outline our modifications (for a detailed account of these modifications, together with application to other models, see [20]). As is common for this type of RG, we preserve the mass gap during the RG steps. We proceed as follows.
1. As a starting point we always use the dual formulations, i.e. Eq. (11) in the present case.
2. As a simplest 3D strip, we take a lattice Λ s = (2 × 2 × L) with periodic BC in both transverse directions and periodic or free BC in the longitudinal direction. Correlation functions are computed on Λ s and on a one-dimensional lattice for all representations.
3. The requirement of mass gap preservation leads to equations for the fixed points where λ j ({t k }) is the eigenvalue of the transfer matrix corresponding to the correlation function in the representation j. This set of equations is treated as the recursion relations for the renormalized coupling constants. 4. The summation over the spins lying in one plane perpendicular to the longitudinal direction introduces all possible interactions between the remaining spins. The transfer matrix is constructed for the evolution of all independent couplings of this general interaction. This substantially reduces the size of the matrix.

5.
Combining the PH RG in the form described above with the cluster decimation approximation of Ref. [21], we can construct new partition and correlation functions on a decimated lattice with double lattice spacing. The next iteration is performed with newly computed constants t (1) j .
In the framework of this approach we have explored the phase structure of the generalized Z(N) model. Only in the case of the standard Potts model we can restrict ourselves to one iteration, since the RG steps do not generate new interactions. For the general case one should perform many iterations to locate precisely the critical points. The critical indices are calculated in the fixed point of the iterations where the critical points of the models of a given universality class are flowing to. Thorough discussion of the general case will be given in [20]. Here we report on the results for the vector LGTs relevant for this paper.
In Table 1 we present the estimates of the critical points β c , both in standard Potts models and in vector models, coming from the PH RG, and compare them with the results of Monte Carlo numerical simulations. In the case of standard Potts models, these simulations were performed in Ref. [8], while for the vector models, they were carried out in this paper. One can see that the PH RG based on the smallest possible strip and combined with the cluster decimation approximation gives quite accurate predictions both for the critical coupling and, as we will see later on, also for the index ν.

Setup of the Monte Carlo simulation
The model under exam in this work is described by the action given in Eqs. (5) and (6), with couplings β 1 = β N −1 ≡ β, β n = 0, n = 2, . . . , N − 2. To study the phase transitions it turns out to be more convenient to simulate the dual spin model, whose action is given in (2) with the coupling constants computed according to Eq. (7). Simulations were performed by means of a cluster algorithm on symmetric lattices L 3 with periodic BC and L in the range 8 -96. For each Monte Carlo run the typical number of generated configurations was 2.5 · 10 6 , the first 10 5 of them being discarded to ensure thermalization. Measurements were taken after every 10 updatings and error bars were estimated by the jackknife method combined with binning.
We considered the following observables: • population S L , where n i is number of s(x) equal to i; • real part of the rotated magnetization M R = |M L | cos(Nψ) and normalized rotated magnetization m ψ = cos(Nψ); • susceptibilities of M L , S L and M R : • Binder cumulants U We computed also the average action and the heat capacity in the vicinity of the critical points.

Critical couplings and their scaling with N
We obtained the critical couplings using the Binder cumulant crossing method described in [22]. In particular, we computed by Monte Carlo simulations the Binder cumulant U (M ) L and its first three derivatives with respect to β for the different lattice sizes, thus allowing to build the function U The values of β c are quoted in the column five of Table 1. The error bars take into account some of the systematic effects, the main one being the dependence of the estimated value of β c on the set of lattice sizes considered in the analysis. Another, less relevant, source of systematics is the uncertainty in the analytic dependence of U (M ) L on β near the critical point.
In [10] the dependence of the 3D Z(N) critical couplings on N was described as With our new data we checked if the value 1.5 is exact and provided the next-order correction to it. Using different forms for the next-order corrections, we found that our data exclude a correction of O(N 0 ), but reveal corrections of the order 1/N 2 , which we take in the form C 1 − cos 2π N . The critical coupling values for the 3D Z(N > 4) vector models were then fitted with the formula giving the following results: A = 1.50122 (7), C = 0.0096(5), χ 2 /d.o.f. = 13.1 (see Fig 1). Despite the large χ 2 , probably due to the underestimation of the error bars of critical couplings, the proposed function nicely interpolates data over a large interval of values of N.

Critical indices and hyperscaling relation
The procedure to determine the critical index ν is also inspired by Ref. [22]: for each lattice size L the known function U The best estimate of ν is found by minimizing the deviation of dU  Table 2, do not differ within error bars.
The critical indices β/ν and γ/ν can be extracted from the FSS analysis of the magnetization M L and its susceptibility χ (M ) L , according to the following fitting functions, which include also the first subleading corrections. The critical index η will then be given by 2 − γ/ν and the hyperscaling relation d = 2β/ν + γ/ν must be satisfied with d = 3.
In Tables 3 and 4 we summarize the results for the critical indices, when the FSS analysis is performed at β fixed at the central value of the determination of β c . Results are given for the cases when the subleading corrections, depending on the exponent δ, are included in neither fit, in both fits or only in one of them. When considered, the exponent δ has been fixed to the value 0.53/ν, as in the 3D XY model. Varying δ in a wide interval does not give any significant change in the results.
In Tables 5, 6, 7 we summarize the values of the critical indices obtained by two alternative methods: (i) performing the fit with the functions given in (20) not at β c , but at the pseudocritical β pc which maximizes the susceptibility χ

Heat capacity and the index α
The critical index α, determined from the ν values obtained in the previous subsection by means of the relation α = 2 − dν, gets negative values for all N ≥ 5, meaning that the transition is of order higher than two. In fact, these negative values are very close to that of the 3D XY model [22]. However, the plots of the heat capacity (see Fig. 2(left)) clearly show that it diverges in the vicinity of the critical point. Moreover, the maxima of the heat capacity and of the susceptibility χ This suggests that a different value for the index ν can be found if a FSS analysis is done on the peak values of the heat capacity, using as fitting function where the possible inclusion of subleading corrections has been taken into account. After α/ν is extracted, from the relation α = 2 − dν the value of ν can be obtained.
In Table 8 we summarize the results for ν determined in the described way, with and without the inclusion of the subleading correction term. When considered, the exponent δ has been fixed to the value 0.5/ν. The error estimates given in the table do not include systematic uncertainties brought by the localization of the maximum of the heat capacity by using analytic continuation, which can be unreliable in regions where the heat capacity changes fast. The constant B in front of the subleading correction appears to be of order unity, B ∼ 0.4 − 0.7, and does not seem to depend on N.
We see that while for N = 2, 4 the resulting ν agrees with the value of ν in the 3D Ising model, and the agreement improves if we include the subleading correction, for N > 4 this is not the case. For N > 4 the difference between the ν values obtained with and without inclusion of subleading corrections is much smaller than for N = 2, 4. The most important fact is, however, that in all cases the ν indices obtained in this way are close to ν ≈ 0.63the critical index for the Ising model. The difference between ν indices obtained from the U (M ) L cumulants and from the heat capacity leads us to conclude that we have two kinds of singularity depending on whether one approaches the critical coupling from above (3D XY model-like singularity) or from below (3D Ising universality class), for N > 4.

Symmetric phase
Another interesting phenomenon we have encountered during our study is the appearance of a symmetric phase just below β c for all N ≥ 5. That such phase exists in the vector Z(6) spin model has been known for a long time [12,13,14]. Here we confirm its existence for all vector Z(N ≥ 5) LGTs. The phase exhibits itself, e.g., in the behavior of magnetization. As an example, we give in Fig 3 the scatter plots of magnetization and rotated magnetization, together with the histogram of the magnetization angle for Z(6) on a 64 3 lattice. One sees that below the critical coupling at β = 2.97 the symmetry is not broken on this lattice. Only starting from approximately β = 2.7 one can observe the appearance of a symmetry-broken phase. In addition, we have studied the behavior of the population susceptibility below β c = 3.00683. Fig 4 shows that it has a second broad maximum, which slowly moves to β c = 3.00683 with increasing lattice size. A similar picture is observed for all N ≥ 5. While, however, for N = 5 the peak of the population susceptibility moves rather fast and practically collapses with the peak at the critical coupling on the largest available lattice L = 96, for larger N the peak stays rather far from the corresponding critical coupling, even for L = 96 (with our data we cannot even exclude a situation when the convergence of the second maximum is logarithmic). We can imagine two scenarios to explain such behavior: 1. This symmetric phase exists only in finite volume. When the lattice size increases, the second maximum approaches the critical coupling and, eventually, the symmetric region shrinks and disappears. The explanation proposed in [13,14] might work in this case, too. Namely, the symmetric phase on the finite lattice is a phase with a very small mass gap and describes a crossover region to the symmetry-broken phase.
2. For N > 5 the second maximum of the population susceptibility stays away from the critical couplings even in the infinite volume limit. In this case it might correspond to some higher order phase transition and the symmetric phase with tiny or even vanishing mass gap exists also in the thermodynamic limit.
In both cases it is tempting to speculate that this symmetric region is reminiscent of the massless phase which appears in these models at finite temperature [17]. Whichever scenario of the above two is realized, one needs to study the models on much larger lattices to uncover it.

Summary
In this paper we have studied the 3D Z (N) LGT at zero temperature aiming at shedding light on the nature of phase transitions in these models for N ≥ 4. This study was based on the exact duality transformations of the gauge models to generalized 3D Z(N) spin models. In Section 2 we presented an overview of the exact relation between couplings of these two models. In Section 3 we have studied the models analytically using a version of the phenomenological RG based on the preservation of the mass gap combined with a cluster decimation approximation. This study provided us with an approximate location of the critical couplings as well as with the value of the index ν. These calculations show that the index ν is approximately the same for N = 2, 4, ν ≈ 0.616, indicating that these two models might belong to the same universality class. We find ν > 2/3 and approximately equal for N > 4. Hence, Z(N > 4) vector LGTs belong to a different universality class. The numerical part of the work has been devoted to the localization of the critical couplings, computation of the various critical indices and check of the hyperscaling relation. The main results can be shortly summarized as follows: • We have determined numerically the position of the critical couplings for various Z(N) models. For N = 2, 3 we find a reasonable agreement with the values quoted in the literature. For larger N, we have significantly improved the values given in [10]. This allowed us to improve the scaling formula for the critical couplings with N.
• The critical indices ν and η derived here for N = 2, 4 suggest that these models are in the universality class of the 3D Ising model, while our results for all N > 4 hint all vector Z(N ≥ 5) LGTs belong to the universality class of the 3D XY model. This is especially evident from the value of the index ν, which stays very close to the XY value, ν ≈ 0.6716, given in [22]. The index α in this case takes a small negative value. It thus follows that a third order phase transition takes place for N ≥ 5.
• A careful investigation of the specific heat and of the index α extracted from it suggests however a more complicated picture of the critical behavior. In this case we find a value which roughly agrees with the value of 3D Ising model for all N studied. The fact that we observe two different values of the index α dependently on whether we approach the critical point from below or from above leads to the conclusion that the first derivative of the free energy could exhibit a cusp in the thermodynamic limit if N > 4.
• Our data also revealed the existence of a symmetric phase for all Z(N) vector LGTs if N > 4. However, substantially larger lattices are required to see if this phase survives the transition to the thermodynamic limit.

Acknowledgments
The work of Ukrainian     Table 3: Critical indices β/ν and γ/ν of 3D Z(N) vector models with N = 2, 4, 5, 6, determined by the fits in Eqs. (20), for different choices of the minimum lattice size L min . The χ 2 of the two fits, given in columns four and six, are the reduced one. Column seven gives the dimension d derived from the hyperscaling relation d = 2β/ν +γ/ν, while column eight contains the values of η = 2 − γ/ν. The three sets of parameters corresponding to the same L min refer to the cases of subleading corrections (term with the exponent δ in Eqs. (20)) included (i) in neither fit, (ii) only in the fit for γ/ν, (iii) in both fits. Table 4: Same as Table 3 for 3D Z(N) vector models with N = 8, 13,20. Table 5: Critical indices of 3D Z(N) vector models obtained from two alternative methods (see the text); the results of the first alternative are given in column two, those of the second alternative (three variants) in columns from three to five. The legenda for the entries in each large cell is the following: η = 2 − γ/ν (first line), β/ν with its χ 2 (second line), γ/ν with its χ 2 (third line) and the dimension d = 2γ/ν + β/ν (fourth line). Subleading corrections (term with the exponent δ) were not included in either fits.   Table 5, with subleading corrections (term with the exponent δ) included only in the fits for γ/ν.   Table 5, with subleading corrections (term with the exponent δ) included in both the fits for γ/ν and for β/ν.  (11) (10) 3.032 (19) 3.02 (2) 3.027 (13)