1. Introduction
Phase transitions are widely observed in nature and represent a fundamental topic in classical statistical physics and thermodynamics. They exemplify some of the most prominent collective behaviors in many-body systems [1] and may also be characterized by abrupt changes in the topology of the N-body configuration space [2, 3]. These transitions often occur suddenly and can lead to catastrophic consequences, such as irreversible climate change [4, 5] or ecological collapse [6, 7]. Consequently, accurately identifying critical points and establishing effective early-warning mechanisms has emerged as a central objective of contemporary research. Previous studies have shown that systems near critical transitions often display distinct spatiotemporal signatures [8-13]. Furthermore, recent efforts have investigated transition warning signal detection using thermodynamic indicators [14] and machine learning techniques [15]. Among these developments, the identification of critical transition signals in spontaneous symmetry-breaking (SSB) phase transitions at finite temperatures has drawn considerable attention. Researchers have systematically examined critical transition warning indicators in classical spin systems undergoing conventional phase transitions. Notably, through the application of microcanonical inflection-point analysis to introduce third-order independent and dependent transitions, in conjunction with geometric order parameters, robust and reliable transition warning signals have been established in the Ising model [16], Baxter-Wu model [17], Potts model [18, 19], and Blume-Capel model [20].
The microcanonical inflection-point analysis identifies phase transition signals by locating the least-sensitive inflection points of derivatives of various orders of the microcanonical entropy [21]. This least-sensitive point exhibits remarkable consistency in both finite-size systems and the thermodynamic limit. In [16], the microcanonical analysis results for the Ising model in the thermodynamic limit demonstrate the same behavior as in the finite-size case. This consistency forms the theoretical basis for applying microcanonical inflection-point analysis to investigate critical transition warning signals in finite-size systems. A third-order dependent transition occurs only in the presence of a lower-order phase transition. Specifically, for a third-order dependent transition to occur, the system must undergo either a first-order or a continuous transition; however, experiencing a first-order or continuous transition does not necessarily entail the presence of a third-order dependent transition. In other words, when a system exhibits a third-order dependent transition, further cooling will ultimately result in either a first-order or a continuous transition. By contrast, a third-order independent transition is independent of any lower-order phase transition, may occur on its own, and generally arises before the onset of global phase transition behavior. Cluster analysis results indicate that both third-order independent and dependent transitions are associated with local changes within the system, whereas first-order and continuous transitions are associated with global changes [16, 19, 20]. Third-order transitions have been shown to lack finite-size effects and do not exhibit divergent behavior in the thermodynamic limit [16-20]. Therefore, they are referred to as a 'pseudo transition'. Our objective is to predict global system behavior by tracking changes in such local behaviors, which has been thoroughly validated in previous studies concerning SSB phase transitions. However, another important class of phase transitions also warrants further investigation——those induced by topological defects, namely the Berezinskii-Kosterlitz-Thouless (BKT) transitions. This study aims to explore transition warning signals associated with such topological phase transitions.
The BKT transition, first proposed in 1972, is a topological phase transition characterized by the unbinding of vortex-antivortex pairs [22, 23]. The two-dimensional XY model on a square lattice undergoes such a transition, which has consequently attracted significant research attention. In recent years, numerous studies have further investigated this phenomenon [24-28]. The two-dimensional q-state clock model is often regarded as a discrete analog of the XY model [29]. Due to its rich and intricate phase behavior, this model has been extensively explored through various analytical and numerical techniques [29-33]. It is well established that for q ≤ 4, the q-state clock model exhibits only a single second-order transition, resembling that of the Ising model, representing a direct transition from an ordered ferromagnetic phase to a disordered paramagnetic phase. In contrast, for q ≥ 5, the system exhibits two distinct BKT-type transitions: one from a phase with true long-range order (TLRO) to a quasi-long-range ordered (QLRO) phase at the critical temperature ${T}_{{\rm{C}}}^{2}$, and another from the QLRO phase to a disordered paramagnetic phase at ${T}_{{\rm{C}}}^{1}$ [30, 34]. Hereafter, G2 and G1 are used to denote these two phase transitions, respectively. As q asymptotically approaches infinity, the q-state clock model converges to the continuous two-dimensional XY model [35-37]. Recent findings have revealed that the G2 in the six-state clock model corresponds to a higher-order SSB transition, whereas the G1 belongs to the BKT transition [33]. Simultaneously, a single dominant clock state emerges beyond the phase transition point in both the Q-state Potts model [38] and the clock model [39], revealing an intriguing phenomenon. Inspired by these insights, we selected the six-state clock model as the basis for our investigation of higher-order phase transitions. This choice enables the exploration of both spontaneous symmetry breaking and BKT transitions, thereby extending the theoretical framework to encompass topological phase transitions that involve the BKT mechanism.
We employ both the microcanonical inflection-point analysis and the canonical ensemble approach to examine the phase-transition behavior of the six-state clock model under periodic boundary conditions on a two-dimensional square lattice and to evaluate the universality of previously proposed critical-transition warning indicators. This approach broadens the applicability of our findings and paves the way for future applications in more realistic systems. The density of states (DOS) for the six-state clock model is obtained via the efficient Wang-Landau (WL) algorithm, while canonical order parameters are calculated using Metropolis Monte Carlo simulations. Based on these data, we examine the phase transition behavior within both microcanonical and canonical frameworks and locate higher-order transitions using canonical order parameters. Section 2 introduces the model and analytical methods; section 3 presents and discusses the main results; and the final section concludes the study and proposes avenues for future investigation.
2. The model and the method
2.1. Six-state clock model
The Hamiltonian of the six-state clock model, defined on a two-dimensional L × L lattice, is expressed as:
$\begin{eqnarray}H=-\displaystyle \sum _{\langle ij\rangle }J{\,{\boldsymbol{S}}}_{i}\cdot {{\boldsymbol{S}}}_{j}=-\displaystyle \sum _{\langle ij\rangle }J\cos ({\theta }_{i}-{\theta }_{j}).\end{eqnarray}$
In this expression, Si and Sj represent the spins at lattice sites i and j, respectively, while J denotes the interaction strength between neighboring spins. In this study, we restrict our analysis to the case of J = 1. The spin orientations are discrete, defined as θi = 2πq/6, where q = 0, 1, …, 5.2.2. Wang-Landau sampling
The WL algorithm [40] is both efficient and broadly applicable for directly estimating the DOS, g(E), of a system. Utilizing g(E) obtained through this method enables the direct investigation of higher-order phase transitions while eliminating the need for post-processing, thereby minimizing errors and improving the accuracy of the results. The DOS is estimated by performing a random walk in energy space based on the criterion of histogram flatness, using the WL sampling method. As the exact form of g(E) is initially unknown, it is initialized as g(E) = 1 for all energy levels and is subsequently refined through iterative updates guided by a transition probability:
$\begin{eqnarray}p({E}_{i}\to {E}_{j})=\min \left(\frac{g({E}_{i})}{g({E}_{j})},1\right),\end{eqnarray}$
where Ei and Ej are the system energies before and after a spin flip, respectively. Upon acceptance of a new energy state, the DOS is modified as: $\begin{eqnarray}g(E)\to g(E)f,\end{eqnarray}$
where f is the modification factor, initially set to f0 = e ≈ 2.71 828, and the corresponding energy histogram H(E) is updated by one. The random walk proceeds by repeatedly flipping spins in the lattice to explore the energy landscape until the histogram meets the flatness condition. A histogram is considered flat when the minimum count across all visited energy levels exceeds 80% of the average count. Once flatness is achieved, the modification factor is reduced according to ${f}_{i+1}=\sqrt{{f}_{i}}$, and the histogram is reset for the next iteration. The simulation terminates when the modification factor is reduced below the threshold ffinal = 1 + 10-8.2.3. Microcanonical infection-point analysis method
Qi and Bachmann integrated microcanonical analysis with the principle of minimal sensitivity to successfully identify and classify first-order and higher-order phase transitions in complex systems [21]. This approach has since been widely applied to various spin models, capturing a broad spectrum of complex phase transition behaviors, including first-order, continuous (second-order), and higher-order transitions.
In physical systems, macroscopic properties are governed collectively by entropy and energy, with the microcanonical entropy serving as a comprehensive descriptor of the system's phase transition behavior. It is defined as
$\begin{eqnarray}S(E)={k}_{{\rm{B}}}{\mathrm{ln}}\,g(E),\end{eqnarray}$
where kB denotes the Boltzmann constant and S(E) is the entropy. Within a single-phase region, the entropy function remains smoothly concave. However, this concavity is characteristically disrupted during a phase transition, resulting in a distinct inflection point on the entropy curve. According to the principle of minimal sensitivity, only the minimally sensitive inflection points are considered physically meaningful.Microcanonical inflection-point analysis identifies significant inflection points of physical relevance by computing higher-order derivatives of the entropy function with respect to energy E, thereby enabling the detection of phase transitions. When a phase transition occurs, the first derivative of entropy, β(E), typically exhibits a positive local minimum, which indicates the transition point in energy space. At this point, the second derivative, $\gamma$(E), generally displays a pronounced peak: a negative peak suggests a continuous (second-order) phase transition, while a positive peak suggests a first-order transition. Furthermore, the third derivative, δ(E), may exhibit a positive local minimum or a negative local maximum, which corresponds to an independent or dependent third-order phase transition, respectively. To detect fourth-order transitions, we analyze the fourth derivative of entropy, denoted by ε(E), which can be used to identify both independent and dependent fourth-order transitions. A summary of this classification framework is provided in table 1.
Table 1. Signal of the order of the transitions. |
| Categories | Even order transitions | Odd order transitions |
|---|---|---|
| Independent | $\frac{{{\rm{d}}}^{2k}S(E)}{{\rm{d}}{E}^{2k}}\lt 0$ | $\tfrac{{{\rm{d}}}^{2k-1}S(E)}{{\rm{d}}{E}^{2k-1}}\gt 0$ |
| Negative maximum | Positive minimum | |
| | ||
| Dependent | $\tfrac{{{\rm{d}}}^{2k}S(E)}{{\rm{d}}{E}^{2k}}\gt 0$ | $\tfrac{{{\rm{d}}}^{2k-1}S(E)}{{\rm{d}}{E}^{2k-1}}\lt 0$ |
| Positive minimum | Negative maximum | |
When computing derivatives directly from discrete microcanonical entropy data, numerical noise can substantially compromise the accuracy of the results. To mitigate this issue and improve data reliability, we employed a multi-step data smoothing strategy. Specifically, we conducted five independent estimations of the DOS and averaged them to obtain a statistically smoothed profile. Based on this averaged DOS, we further applied the Bézier algorithm to construct a continuous and differentiable entropy curve [41], following the procedure described in our previous work. Finally, we repeated the entire process four times to obtain independent datasets and performed error analysis by evaluating the standard deviations across these datasets.
2.4. Canonical analysis method
To gain a more comprehensive understanding of the model, we performed a canonical ensemble analysis. On the one hand, we investigated the system's phase behavior using geometric order parameters obtained from spin configurations generated using Metropolis Monte Carlo sampling [42]. This geometric indicator, the average cluster perimeter, is defined consistently with our previous studies [19, 20]. On the other hand, we additionally employed an order parameter associated with the ${{\mathbb{Z}}}_{6}$ symmetry-breaking transition to characterize the nature of the phase transitions and detect potential higher-order behavior [43].
The average cluster perimeter, P, is defined for a given spin configuration X as the average perimeter of all clusters containing more than one spin:
$\begin{eqnarray}P=\frac{1}{n}\displaystyle \sum _{l}{C}_{l},\end{eqnarray}$
where n is the number of such clusters in configuration X, l indexes each cluster, and Cl is the perimeter of cluster l. The thermal average of this quantity is given by: $\begin{eqnarray}\langle P\rangle =\displaystyle \frac{1}{Z}\displaystyle \sum _{X}P(X){{\rm{e}}}^{-E(X)/{k}_{{\rm{B}}}T},\end{eqnarray}$
where T is the temperature, and $Z={\sum }_{X}\exp \left[-E(X)/{k}_{{\rm{B}}}T\right]$ is the canonical partition function. Our primary interest lies in the rate of change of this quantity with temperature, which we define as $\begin{eqnarray}{D}_{\langle P\rangle }=\displaystyle \frac{{\rm{d}}\langle P\rangle }{{\rm{d}}T}.\end{eqnarray}$
A peak in D‹P› indicates the location of the critical transition, whereas a local minimum beyond this point corresponds to the higher-order dependent transition. This criterion is consistent with the one proposed by Sitarachu and Bachmann in their analysis of average cluster size [16].The phase transition associated with the ${{\mathbb{Z}}}_{6}$ symmetry breaking is characterized by the following order parameter [43]:
$\begin{eqnarray}{M}_{{{\mathbb{Z}}}_{6}}=\frac{1}{N}{\left\{{\left[\displaystyle \sum _{i}\cos ({\varphi }_{i})\right]}^{2}+{\left[\displaystyle \sum _{i}\sin ({\varphi }_{i})\right]}^{2}\right\}}^{1/2}.\end{eqnarray}$
Our analysis focuses on the temperature dependence of the ${{\mathbb{Z}}}_{6}$ symmetry. As the temperature increases, the degree of ${{\mathbb{Z}}}_{6}$ symmetry gradually diminishes. The low-temperature phase transition is characterized by spontaneous ${{\mathbb{Z}}}_{6}$ symmetry breaking, which is signaled by a pronounced peak in the derivative of the corresponding order parameter. The high-temperature phase transition, associated with the breaking of continuous U(1) symmetry, is primarily driven by the continuous symmetry but also gives rise to a secondary, more pronounced peak in the rate of ${{\mathbb{Z}}}_{6}$ symmetry breaking. Between these two peaks, there exists an intermediate regime in which the change in ${{\mathbb{Z}}}_{6}$ symmetry is relatively smooth. We hypothesize that this intermediate regime is indicative of a higher-order dependent transition. We use $\tfrac{{\rm{d}}\langle {M}_{{{\rm{{\mathbb{Z}}}}}_{6}}\rangle }{{\rm{d}}T}$ to represent the rate of change of ${M}_{{{\mathbb{Z}}}_{6}}$ (the first derivative of ${M}_{{{\mathbb{Z}}}_{6}}$) as a diagnostic tool for identifying and characterizing phase transitions. It exhibits a negative local minimum at the phase transition point and a negative local maximum at the higher-order transition point.3. Results
3.1. Traditional transitions
In this study, the DOS of the system was estimated using the WL sampling method. To balance computational efficiency and accuracy, we restricted our simulations to system sizes of L = 20, 30, 40, 50, and 60. The phase transition points were determined via two independent and complementary approaches: microcanonical inflection-point analysis and canonical ensemble analysis by analyzing the temperature dependence of the specific heat. The specific heat is calculated according to the expression [43]:
$\begin{eqnarray}{C}_{V}=\displaystyle \frac{\langle {E}^{2}\rangle -{\langle E\rangle }^{2}}{N{k}_{{\rm{B}}}{T}^{2}},\end{eqnarray}$
where CV is normalized by the total number of spins N, which allows for consistent comparison across system sizes in a single plot.In figure 1(a), we present the specific heat curves obtained in the canonical ensemble for identifying the phase transition points. Clearly, two distinct peaks emerge in the specific heat, indicating the presence of two separate phase transitions in the system. Vertical blue and green lines are used to indicate the locations of the G2 and G1, respectively. As shown in figure 1(a), the specific heat peak near the G2 increases with system size and eventually approaches a system-size-independent value. In contrast, the specific heat in the G1 regime exhibits a non-monotonic behavior: it first increases and then decreases as the system size increases. This trend is evident in the inset (c), where the peak value at L = 30 exceeds those at L = 40, L = 50 and L = 60. Such behavior has also been reported in previous studies [27, 43], suggesting that it represents a distinct signature of a BKT-type transition. In [44], it was reported that in the XY model, the specific-heat peak temperature is approximately 17% higher than the BKT transition point. Furthermore, we estimated the 8-state clock model using data from [31] and found that the specific-heat peak temperature is approximately 12% higher than the BKT transition point. In our subsequent analysis, the phase transition temperature in the limit of an infinite system size was determined using ${M}_{{{\mathbb{Z}}}_{6}}$, and a red arrow was added in figure 1(a) to indicate this value. We found that the specific-heat peak temperature at the BKT transition is about 6% higher than the actual BKT transition point. This behavior is consistent with the characteristic features of the BKT transition, thereby demonstrating the reliability of our results.
Figure 1. (Specific heat and finite-size scaling analysis). (a) shows the specific heat curves computed for various system sizes L, where the presence of two pronounced peaks reveals evidence for two distinct phase transitions in the system. The temperature indicated by the red arrow is the phase transition temperature determined by ${M}_{{{\mathbb{Z}}}_{6}}$ in the limit of an infinite system size. Inset (c) provides a close-up view of the specific heat in the vicinity of the G1. (b) presents the Tc(1/L) plot of the six-state clock model for L = 20-60. The Tc(1/L) values are derived from finite-size scaling analysis of the peak positions of the specific heat. In the limit L → ∞, the results yield ${T}_{{\rm{C}}}^{1}=1.0070(20)$ and ${T}_{{\rm{C}}}^{2}=0.6065(67)$. |
From the positions of the maxima of the specific heat, we obtain the critical temperatures ${T}_{{\rm{C}}}^{1}=1.0070$ and ${T}_{{\rm{C}}}^{2}=0.6065$ in the limit L → ∞, as shown in figure 1(b). Although it is known that phase transitions do not occur precisely at the peaks of Cv, these positions serve as a reference. As shown in figure 3(b) of [32], where the phase transition points were determined from the specific heat peaks as ${T}_{{\rm{C}}}^{1}=1.017$ and ${T}_{{\rm{C}}}^{2}=0.6068$, our results are in close agreement. This suggests that the DOS obtained from our WL simulations is accurate, further confirming that the results from the microcanonical analysis using this DOS are reliable. In a subsequent table, we summarize the phase transition points reported in [32], which differ slightly from the transition temperatures presented here, as the transition points in the table were determined using the Fisher zero method.
In addition, we conducted a microcanonical inflection-point analysis of the data. figure 2(a) shows the S(E) for the six-state clock model. To facilitate comparison across different system sizes, the S(E) has been normalized by the system size N, enabling direct comparison across different system sizes in a single plot. Figure 2(b) shows the first derivative of the entropy, β(E), which corresponds to the inverse temperature. Figures 2(c) and (d) present the second and third derivatives of the entropy, $\gamma$(E) and δ(E), respectively.
Figure 2. (Microcanonical inflection point analysis of the six-State clock model). The figure presents the results of microcanonical analysis for the six-state clock model. (a) depicts the entropy profile as a function of energy E, while (b) presents the profile of its first derivative, β(E). (c) shows the evolution of the second derivative, $\gamma$(E), where two negative maxima are identified as indicators of two separate phase transitions. (d) presents the evolution of the third derivative, δ(E). Inset (e) shows a zoomed-in view of the fourth derivative of the density of states, ε(E), revealing a signature consistent with a fourth-order dependent transition. |
In figure 2(c), two phase transitions are identified, which correspond to two distinct negative peaks in the curve. By locating the associated energy levels at which the transitions occur and referencing these energies in panel (b), the corresponding inverse temperatures are obtained, from which the transition temperatures are subsequently determined. In figure 2(b), these transition points are indicated with vertical blue and green lines, which mark the G2 and G1, respectively. The corresponding transition temperatures are annotated using the same color scheme. According to [33], the G2 is identified as a higher-order SSB phase transition, whereas the G1 belongs to the BKT class. Based on the microcanonical inflection-point analysis, both the G2 and G1 are classified as continuous transitions. This interpretation is further supported by all the order parameters employed in our subsequent analysis, each of which exhibits smooth and continuous behavior in the vicinity of the transition points.
In table 2, we summarize the phase transition points determined from both the microcanonical inflection-point analysis and the peak positions of the specific heat. $1/{\beta }_{{\rm{C}}}^{1}$ and $1/{\beta }_{{\rm{C}}}^{2}$ represent the transition temperatures of G1 and G2 extracted from specific heat calculations within the canonical ensemble. The numbers in parentheses indicate the associated uncertainty in the numerical estimation. Similarly, $1/{\beta }_{{\rm{m}}}^{1}$ and $1/{\beta }_{{\rm{m}}}^{2}$ denote the transition temperatures of G1 and G2 identified via the microcanonical method. E/N denotes the energy per site corresponding to the phase transition point in the microcanonical ensemble. All data have been compiled and presented in table 2.
Table 2. Phase transition point information determined by microcanonical inflection-point analysis and canonical specific heat analysis. |
| L = 20 | L = 30 | L = 40 | L = 50 | L = 60 | |
|---|---|---|---|---|---|
| $1/{\beta }_{{\rm{C}}}^{1}$ | 1.0914(3) | 1.0652(10) | 1.0514(15) | 1.0411(4) | 1.0336(8) |
| $1/{\beta }_{{\rm{C}}}^{2}$ | 0.5927(14) | 0.5969(13) | 0.6002(17) | 0.6013(12) | 0.6013(15) |
| $1/{\beta }_{{\rm{m}}}^{1}$ | 1.1226 | 1.0910 | 1.0767 | 1.0848 | 1.1306 |
| E1/N | -1.17125 | -1.20222 | -1.2175 | -1.2026 | -1.1374 |
| $1/{\beta }_{{\rm{m}}}^{2}$ | 0.6195 | 0.6231 | 0.6245 | 0.6266 | 0.6277 |
| E2/N | -1.75125 | -1.73833 | -1.73375 | -1.7296 | -1.7275 |
To enable a more thorough investigation of the phase transition behavior and facilitate comparison with the microcanonical transition energies, we show in figure 3 the canonical probability distribution $P(E,T)\sim g(E){{\rm{e}}}^{-E/{k}_{{\rm{B}}}T}$ near the transition temperatures. Figure 3(a) illustrates the distribution P around the G2 point, while figure 3(b) corresponds to the G1 region. In figure 3, the red data points indicate the probability distribution evaluated precisely at the phase transition temperature, while the black and blue data points correspond to temperatures just below and above the transition point, respectively. In each distribution, the peak position signifies the energy level associated with the phase transition. The transition energy levels are marked with vertical blue and green lines, with color-coded labels denoting the corresponding energy values.
Figure 3. (Plot of the canonical distribution, $P(E,T)\sim g(E){{\rm{e}}}^{-E/{k}_{{\rm{B}}}T}$). For the six-state clock model with system size L = 60, the canonical distribution $P(E,T)\sim g(E){{\rm{e}}}^{-E/{k}_{{\rm{B}}}T}$ is examined at the phase transition temperature T. (a) illustrates the distribution of the order parameter near the G2 at ${T}_{{\rm{C}}}^{2}$, while (b) depicts the distribution near the G1 at ${T}_{{\rm{C}}}^{1}$. Red points indicate the distribution corresponding to the transition temperature, black points represent data slightly below it, and blue points represent data slightly above it. The energy levels at which the phase transitions occur are indicated by blue and green markers. |
In figure 3, the transition energy levels of the two phase transition points in the canonical ensemble for system size L = 60 are presented. In table 2, the corresponding transition energy levels of the two phase transition points in the microcanonical ensemble for L = 60 are also listed. A straightforward calculation shows that, in the G2 regime, the transition energy levels obtained from the canonical and microcanonical methods are in excellent agreement, with an energy difference of only -0.0356, which is nearly negligible. In contrast, under the same conditions, a larger discrepancy of -0.1405 is observed in the G1 regime. This deviation can be attributed primarily to two factors: first, the inherent theoretical and methodological discrepancies between the canonical and microcanonical analyzes; and second, the fact that the G1 is driven by the unbinding of vortex-antivortex pairs, which leads to a more intricate energy landscape characterized by extensive level merging. Since the microcanonical inflection-point method is highly sensitive to the detailed structure of the energy spectrum, such complexity can result in larger deviations and potential information loss in the G1 region. The results obtained here are thus deemed accurate and reliable. Furthermore, a comparison with previously reported transition points in the literature, as summarized in table 4, reveals that, although there is some numerical discrepancy, it can be attributed to the smaller system sizes used in our simulations. When compared with transition points obtained from WL simulations, the difference is relatively small.
Table 3. Transition points obtained from different canonical order parameters. |
| ACV | ${M}_{{{\mathbb{Z}}}_{6}}$ | ||||
|---|---|---|---|---|---|
| L | Tde | ${T}_{{\rm{C}}}^{2}$ | Tde | ${T}_{{\rm{C}}}^{2}$ | ${T}_{{\rm{C}}}^{1}$ |
| 20 | 0.901 | 0.603 | 0.789 | 0.615 | 1.098 |
| 30 | 0.871 | 0.611 | 0.786 | 0.623 | 1.065 |
| 40 | 0.859 | 0.614 | 0.788 | 0.627 | 1.046 |
| 50 | 0.851 | 0.616 | 0.786 | 0.630 | 1.033 |
| 60 | 0.845 | 0.617 | 0.787 | 0.634 | 1.025 |
| 70 | 0.843 | 0.618 | 0.789 | 0.636 | 1.016 |
| 80 | 0.840 | 0.619 | 0.791 | 0.638 | 1.010 |
| 90 | 0.838 | 0.620 | 0.783 | 0.640 | 1.005 |
| 100 | 0.837 | 0.619 | 0.784 | 0.640 | 1.002 |
| 110 | 0.836 | 0.620 | 0.793 | 0.642 | 0.998 |
| 120 | 0.836 | 0.620 | 0.789 | 0.643 | 0.996 |
Table 4. Comparison of estimated critical temperatures ${T}_{{\rm{C}}}^{2}$ and ${T}_{{\rm{C}}}^{1}$ by various methods. |
| References | Method | L or m | ${T}_{{\rm{C}}}^{2}$ | ${T}_{{\rm{C}}}^{1}$ |
|---|---|---|---|---|
| Challa and Landau [29] (1986) | MC | L = 72 | 0.68(2) | 0.92(1) |
| Hwang [32] (2009) | Wang-Landau MC | L = 28 | 0.632(2) | 0.997(2) |
| Baek et al [45, 46] (2010) | Wolff MC | L = 512 | —— | 0.9020(5) |
| Kumano et al [47] (2013) | Boundary-flip MC | L = 256 | 0.700(4) | —— |
| Chen et al [48] (2017) | HOTRG | m = 15 | 0.6658(5) | 0.8804(2) |
| Suruengan et al [49] (2019) | Swendsen-Wang MC | L = 512 | 0.705(8) | 0.898(5) |
| Hong and Kim [50] (2020) | HOTRG | L = 128 | 0.693 | 0.904 |
| Li et al [37] (2020) | VUMPS | m = 250 | 0.694(1) | 0.9127(5) |
| Ueda and Okunishi et al [51] (2020) | CTMRG (Correlation length, etc) | m = 768 | 0.694(3) | 0.908(3) |
| CTMRG (Entanglement spectrum) | m = 768 | 0.693 | 0.900 | |
| Shiina and Mori [52] (2020) | Machine-Learning | L = 64 | 0.66(1) | 0.93(1) |
| Tseng and Jiang [53] (2023) | Supervised neural network | L = 256 | 0.66(2) | 0.893(15) |
| This work | MIPA | L = 60 | 0.6277 | 1.1306 |
| ACV | L = 120 | 0.620 | —— | |
| ${M}_{{{\mathbb{Z}}}_{6}}$ | L = 120 | 0.643 | 0.996 |
3.2. Higher order transitions
The δ(E) curve in figure 2(d) is employed to identify third-order phase transitions. Identifying such higher-order transitions facilitates a more detailed analysis and efficient utilization of the data, thereby enhancing our overall understanding of the system's behavior.
In figure 2(d), we investigated potential indicators of third-order phase transitions in the system; however, no conclusive evidence for such a transition was observed. In inset (e), we analyze the behavior of the local fourth derivative of the DOS. A positive local minimum is detected in the deep-blue curve near E/N = -1.56 (T ≈ 0.8). However, the error in this region is considerable, suggesting the possibility of overfitting. A closer examination demonstrates that this positive minimum also appears in smaller systems, where the associated error is significantly reduced. This discrepancy can be tentatively attributed to the amplified numerical uncertainties inherent in microcanonical inflection-point analysis for larger systems. To further investigate the system's phase behavior, we employ two representative order parameters——the average cluster perimeter and the MZ6 symmetry. Our subsequent analysis provides evidence for a higher-order transition at this location, which, when combined with the microcanonical inflection-point analysis, can be identified as a fourth-order transition.
In figures 4(a) and (b), we show the variation of the average cluster perimeter and the ${M}_{{{\mathbb{Z}}}_{6}}$ symmetry with temperature, respectively. Figures 4(c) and (d) display the corresponding first derivatives with respect to temperature. Our primary interest lies in the rate of change of these order parameters, and thus we focus primarily on their derivative analysis. In the figure, vertical blue and green lines mark the locations of the G2 and G1, respectively, while the vertical orange line indicates the location of the higher-order dependent transition.
Figure 4. (Analysis of the Average Cluster Perimeter and ${M}_{{{\mathbb{Z}}}_{6}}$ Symmetry ). (a) and (b) show the temperature dependence of the average cluster perimeter and the ${M}_{{{\mathbb{Z}}}_{6}}$ symmetry, respectively, while (c) and (d) display the corresponding first derivatives of these order parameters. In the figures, the low-temperature phase transition point, high-temperature phase transition point, and the dependent transition point are indicated by blue, green, and orange markers, respectively, and the corresponding temperatures are labeled using the same colors. |
Figures 4(a) and (c) display the temperature dependence of the average cluster perimeter and its first-order derivative, respectively. The first derivative exhibits a local maximum followed by a negative local minimum. The local maximum corresponds to the G2, while the negative minimum marks the location of a higher-order dependent transition. In conjunction with the microcanonical inflection-point analysis, this dependent transition can be classified as a fourth-order dependent transition. At the G2, the average cluster perimeter increases most rapidly, indicating that the system's order is being rapidly disrupted. This behavior reflects the fragmentation of the dominant cluster, signaling the onset of a phase transition from an ordered to a quasi-ordered phase. As the largest cluster disintegrates, single-spin clusters and small clusters begin to emerge. Driven by spin-spin coupling interactions, these 'outlier' tend to either recombine or attract neighboring spins, forming new cluster structures. As a result, just beyond the phase transition point, the average cluster perimeter initially increases and then decreases with increasing temperature.With a further increase in temperature, the number of small clusters and single-spin domains continues to rise, which slows the rate of decrease in the average perimeter. Consequently, the first derivative of the average cluster perimeter exhibits a negative local minimum, indicating the location of the fourth-order dependent transition.
However, no corresponding signature is observed near the G1. This is because the G1 belongs to the BKT universality class, characterized by the unbinding of vortex-antivortex pairs. These vortex structures typically differ from their nearest neighbors and, when included in perimeter calculations, their unbinding leads to only subtle variations in the perimeter. As temperature increases further, the number of clusters steadily grows, leading to a gradual decrease in the average perimeter. However, the rate of this decrease becomes less pronounced, and the first derivative does not display a distinct peak. Therefore, no clear signal of a phase transition or higher-order transition is detected in the G1 regime.
To gain further insight into the temperature-dependent behavior of the average cluster perimeter, we analyzed the spin configurations at different temperatures, as shown in figure 5. Figures 5(a), (b), and (c) depict the system in the ordered, quasi-ordered, and high-temperature disordered phases, respectively. In figure 5(b), the observed vortices consist of several small clusters with gradually varying spin orientations, forming distinct vortex structures. This observation is consistent with our earlier description of 'outlier' clusters, which tend to recombine or attract neighboring spins to form new cluster configurations. These small clusters emerge, single-spin clusters appear first. These single-spin clusters subsequently attract neighboring spins, aligning them with their own orientation or reducing angular disparities, thereby generating small clusters with gradually varying spin directions, which ultimately culminate in the formation of vortices.
Figure 5. (Configuration distribution map of the six-state clock model). (a), (b), and (c) present representative snapshots of the six-state clock model of size L = 60 under varying temperatures. (a) shows the system in a low-temperature ordered phase, while (b) represents a quasi-ordered phase, with red boxes indicating the locations of vortices. (c) illustrates the system in a high-temperature disordered phase. |
Meanwhile, the location of the fourth-order dependent transition is further supported by the evolution of the ${{\mathbb{Z}}}_{6}$ symmetry-breaking order parameter ${M}_{{{\mathbb{Z}}}_{6}}$. As shown in figures 4(b) and (d), the ${M}_{{{\mathbb{Z}}}_{6}}$ symmetry and its first-order derivative are plotted as functions of temperature. In figure 4(d), we observe that the first derivative of ${M}_{{{\mathbb{Z}}}_{6}}$ displays two negative local minima and one negative local maximum. These extrema correspond, respectively, to the two phase transition points and the location of the fourth-order dependent transition associated with the G2. According to [33], the G2 is characterized by the spontaneous breaking of the ${M}_{{{\mathbb{Z}}}_{6}}$ symmetry, which manifests as a sharp decline in the ${M}_{{{\mathbb{Z}}}_{6}}$ order parameter. In contrast, the G1 results from the breaking of the U(1) symmetry, which likewise induces a marked suppression of the ${M}_{{{\mathbb{Z}}}_{6}}$ order parameter. This reduction is even more pronounced than the suppression observed during the spontaneous ${M}_{{{\mathbb{Z}}}_{6}}$ symmetry breaking. As a result, the negative minimum of the first derivative of ${M}_{{{\mathbb{Z}}}_{6}}$ at the G1 is more pronounced than the one observed at the G2 as shown in figure 4(d). To observe two negative local minima, a local maximum must necessarily occur between them. In figure 4(d), we clearly observe this negative local maximum, and its mechanism is similar to the explanation based on cluster perimeter: as the temperature increases, the system's symmetry gradually weakens. However, when spins reaggregate due to coupling interactions, the rate of symmetry breaking slows down, leading to an increase in $\tfrac{{\rm{d}}\langle {M}_{{{\rm{{\mathbb{Z}}}}}_{6}}\rangle }{{\rm{d}}T}$. Subsequently, as the temperature continues to rise to the second phase transition point, symmetry rapidly decreases again, resulting in the appearance of the negative local maximum. This feature corresponds to the position of the fourth-order dependent transition, and this position is close to the location determined by the average cluster perimeter analysis of the fourth-order transition.
Figure 6 presents the results of our analysis for the system with L = 40 clusters. We conducted 200 000 Monte Carlo simulations, followed by statistical analysis of the final 100 000 MC steps. Figure 6(a) illustrates the cluster distribution near ${T}_{{\rm{C}}}^{2}$. It is evident that at this temperature, the giant cluster dominates the system. However, as the temperature surpasses ${T}_{{\rm{C}}}^{2}$, the number of giant clusters significantly decreases, and the size of the largest cluster is reduced to approximately 1200. Simultaneously, clusters with sizes ranging from 400 to 800 (larger clusters) begin to emerge, as shown in figure 6(b). This indicates a transition from an ordered to a quasi-ordered phase at ${T}_{{\rm{C}}}^{2}$. Figures 6(c) and (d) illustrate the cluster distributions near the fourth-order dependent transition. It is evident that at this stage, the number of larger clusters remains relatively high. These larger clusters result from the disintegration of giant clusters and the merging of 'outlier' clusters due to spin-spin coupling interactions. This is consistent with the previous explanation of the average cluster perimeter and ${M}_{{{\mathbb{Z}}}_{6}}$ symmetry, reflecting a localized phenomenon. As the temperature increases, the number of larger clusters gradually decreases and eventually disappears, as shown in figures 6(e) and (f). Figure 6(f), it is evident that at ${T}_{{\rm{C}}}^{1}$, the larger clusters almost completely vanish, indicating a transition from a quasi-ordered to a disordered phase.
Figure 6. (Non-normalized occurrence frequency of clusters with different sizes). (a), (b), (c), (d), (e), and (f) represent the cluster distributions of a system with size L = 40 at temperatures T = 0.6270, T = 0.6500, T = 0.7880, T = 0.8000, T = 0.9000, and T = 1.0460, respectively. A total of 200 000 Monte Carlo steps were performed, with data from only the last 100 000 steps recorded for the cluster distribution. |
Why can the fourth-order dependent transition serve as an early-warning indicator for critical transitions? As shown in figure 4(d), a negative local maximum appears at the fourth-order dependent transition in $\tfrac{{\rm{d}}\langle {M}_{{{\rm{{\mathbb{Z}}}}}_{6}}\rangle }{{\rm{d}}T}$. With further temperature reduction, a subsequent negative local minimum is expected to emerge at G2. This occurs due to the system's rapid recovery of ${M}_{{{\mathbb{Z}}}_{6}}$ symmetry as the temperature continues to decrease. As the temperature increases, the ${M}_{{{\mathbb{Z}}}_{6}}$ symmetry inevitably reaches its minimum value. Near G1, the rate of decrease of ${M}_{{{\mathbb{Z}}}_{6}}$ is fastest, resulting in a pronounced negative local minimum in its first derivative. As shown in figure 6, the number of larger clusters remains relatively high near the fourth-order dependent transition. As the temperature decreases, the number of larger clusters declines, and giant clusters gradually dominate, occurring at the G2 phase transition. As the temperature increases, the number of larger clusters also decreases, and smaller clusters gradually dominate, occurring at the G1 phase transition. This suggests that the fourth-order dependent transition is a key point in the quasi-ordered region, and detecting it can help predict the impending transitions. Therefore, we use the fourth-order dependent transition as an early-warning indicator for the critical transitions in the six-state clock model.
In summary, the location of the fourth-order dependent transition serves as a reliable indicator for the critical transition warning of the two primary phase transition points. The transition temperatures obtained from the canonical order parameters——including the two primary phase transitions and the fourth-order dependent transition——are compiled in table 3. Here, ${T}_{{\rm{C}}}^{1}$ denotes the G1 phase temperature, ${T}_{{\rm{C}}}^{2}$ denotes the G2 phase temperature, and Tde indicates the position of the fourth-order dependent transition. As shown in the results in table 3, the locations of the fourth-order dependent transition determined by the two methods are not exactly the same. Upon observing the difference between them, we note that as the system size increases, the difference decreases from approximately 0.11 to about 0.04. This discrepancy primarily arises from the simulation noise in the two types of order parameters, which we consider to correspond to the same location. Additionally, we conducted a comparison between the transition points obtained from the canonical order parameters and those reported in the literature. This comparative analysis is summarized in table 4, where MC refers to Monte Carlo methods, HOTRG indicates the tensor renormalization group method with higher-order singular value decomposition, and VUMPS refers to the variational uniform matrix product state method. The parameter m represents the bond (truncation) dimension.
In table 4, it can be seen that the transition points obtained using three different methods are consistent. However, a comparison with the relevant literature reveals that the reported transition points are approximately ${T}_{{\rm{C}}}^{1}=0.94$ and ${T}_{{\rm{C}}}^{2}=0.67$, which show some discrepancy with our results. This difference is primarily due to the fact that our simulations were performed on finite-size systems. In table 4, the largest system size used to determine the transition points via microcanonical inflection-point analysis and specific heat analysis is L = 60, whereas the analyzes based on the average cluster perimeter and the ${M}_{{{\mathbb{Z}}}_{6}}$ symmetry employ system sizes up to L = 120. We selected representative order parameters, such as ${M}_{{{\mathbb{Z}}}_{6}}$, and performed extrapolations to estimate the transition temperatures in the limit L → ∞, as shown in figure 7. In this limit, the G1 and G2 transition temperature are ${T}_{{\rm{C}}}^{1}=0.9690$ and ${T}_{{\rm{C}}}^{2}=0.6516$, which corroborate the results reported in the literature in table 4, further confirming the reliability of our findings.
Figure 7. (Finite-size scaling analysis using ${M}_{{{\mathbb{Z}}}_{6}}$ symmetry). The figure presents the Tc versus (1/L) plot of the six-state clock model for L = 50-120. The Tc values are derived from finite-size scaling analysis using the ${M}_{{{\mathbb{Z}}}_{6}}$ symmetry, yielding ${T}_{{\rm{C}}}^{1}=0.9690(11)$ and ${T}_{{\rm{C}}}^{2}=0.6516(56)$ in the limit L → ∞. |
4. Summary
In this paper, we employed the WL sampling technique to estimate the DOS for the six-state clock model and applied the Metropolis algorithm to generate spin configurations, from which key observables——such as isolated spins, the average cluster perimeter, and the ${M}_{{{\mathbb{Z}}}_{6}}$ symmetry——were extracted. Based on the data obtained by both methods, we performed analyzes within the frameworks of both microcanonical and canonical ensembles, facilitating the accurate identification of the phase transition points in the G2 and G1. The close agreement between the results from the two sampling approaches further substantiates the validity and robustness of our findings. Moreover, using the microcanonical inflection-point analysis, we classified the G2 and G1 as continuous in nature. This conclusion is further supported by the behavior of the canonical order parameters.
A fourth-order dependent transition was also identified in the G2 based on two canonical order parameters: the average cluster perimeter and the ${M}_{{{\mathbb{Z}}}_{6}}$ symmetry. The transition points obtained from these two measures are in close agreement, supporting both the reliability of the data and the robustness of the analytical framework. More importantly, the analysis of the ${M}_{{{\mathbb{Z}}}_{6}}$ symmetry suggests that this fourth-order dependent transition can potentially serve as an early indicator of both the G2 and G1.
In general, the occurrence of a phase transition necessitates that the system approach the thermodynamic limit, in which physical quantities display non-analytic behavior. In our study, however, simulations were conducted on finite-size systems, which inherently do not exhibit true singularities. Nevertheless, we observed abrupt or sharp variations resembling genuine phase transitions; these phenomena are commonly referred to as pseudo phase transitions. Furthermore, the higher-order (third- and fourth-order) transitions, identified via microcanonical inflection-point analysis, are largely unaffected by finite-size effects; rather, they manifest primarily as localized phenomena. Therefore, strictly speaking, these transitions cannot be classified as true phase transitions. Nevertheless, cluster-based analyzes reveal that such higher-order transitions can indeed partition the system into two distinct states; therefore, they are referred to as pseudo transitions.
These findings contribute to a deeper understanding of phase transition behavior in the clock model and offer valuable insights for future investigations of phase transitions. Moreover, they provide a solid theoretical foundation for the study of critical transition warning signals. We plan to further develop and apply this approach in future investigations, with continued efforts to identify order parameters that may serve as reliable early-warning indicators of critical transitions.