Monte Carlo Renormalization Flows in the Space of Relevant and Irrelevant Operators: Application to Three-Dimensional Clock Models

TL;DR

This paper introduces a Monte Carlo-based method to visualize and analyze renormalization group flows in three-dimensional clock models, effectively identifying the relevance of operators and the crossover between different fixed points.

Contribution

The paper presents a novel Monte Carlo approach for visualizing and quantifying RG flows, especially useful for complex flows involving dangerously irrelevant operators.

Findings

Method outperforms standard techniques in determining scaling dimensions.

Successfully analyzes the crossover from U(1) to Z_q fixed points.

Applicable to quantum-critical points with irrelevant operators.

Abstract

We present a way to visualize and quantify renormalization group flows in a space of observables computed using Monte Carlo simulations. We apply the method to classical three-dimensional clock models, i.e., the planar (XY) spin model perturbed by a $Z_q$ symmetric anisotropy field. The method performs significantly better than standard techniques for determining the scaling dimension $y_q$ of the $Z_q$ field at the critical point if it is irrelevant ($q\ge4$). Furthermore, we analyze all stages of the complex renormalization flow, including the cross-over from the U(1) Nambu-Goldstone fixed point to the ultimate $Z_q$ symmetry-breaking fixed point due to the relevance of the $Z_q$ field inside the ordered phase. We expect our method to be particularly useful in the context of quantum-critical points with inherent dangerously irrelevant operators that cannot be tuned away microscopically but whose renormalization flows can be analyzed exactly as we do here for the clock models.

Figures

Tables

Table 1. Table 1: Scaling dimensions y q subscript 𝑦 𝑞 y_{q} of the Z q subscript 𝑍 𝑞 Z_{q} field for q = 4 , 5 , 6 𝑞 4 5 6 q=4,5,6 . The numbers within parenthesis indicate the statistical errors (one standard deviation) of the preceding digit.

	4	5	6
Ref. Oshikawa00	0.2	1.5	3.0
Ref. Leonard15	0.114	1.16	2.29
Refs. Hasenbusch11 ; Okubo15	0.108(6)	1.25	2.5
Ref. Banerjee18	0.128(6)	1.265(6)	2.509(7)
This work	0.114(2)	1.27(1)	2.55(6)

Table 2. Table SI: Critical temperatures for various q 𝑞 q and h ℎ h values. The underlying analysis is presented in Fig.  S1 .

	4	5	6
1.0	2.20465(1)
2.0		2.20239(1)
5.0		2.20357(1)
$\infty$		2.20502(1)	2.20201(1)

Table 3. Table SII: Scaling dimension of the Z q subscript 𝑍 𝑞 Z_{q} field based on the fits in Fig.  S4 and compared with previous numerical results.

$q$	1	2	3
$y_q$	2.481(1)	1.7677(4)	0.876(13)
	2.4810(3) Campostrini06	1.7639(11) Hasenbusch11	0.8915(20) Hasenbusch11

Equations

$$f_s(t,h,L)=L^{-d}{F}_s(t{L}^{1/\nu },t{L}^{1/{\nu }^{\prime }},h{L}^y,\lambda L^{-\omega }),$$

$$H=-\sum _{⟨i,j⟩}\cos ({\theta }_i-{\theta }_j)-h\sum _i\cos (q{\theta }_i),$$

$${\varphi }_q=⟨\cos (q\Theta )⟩,$$

$${\varphi }_q=L^{y_q}\Phi (t{L}^{1/\nu },t{L}^{1/{\nu }_q^{\prime }}),$$

$${\varphi }_q\propto L^{y_q}(1+t{L}^{1/\nu }),$$

$$U=U(t{L}^{1/\nu })=U_{XY}+t{L}^{1/\nu }+L^{-\omega },$$

$$d_1\propto \sqrt{(t{L}^{1/\nu }+L^{-\omega }{)}^2+L^{2y_q}(1+t{L}^{1/\nu }{)}^2}.$$

$$D_1\propto t^{\frac{\omega }{1/\nu +\omega }}=t^{0.345(6)},L_1\propto t^{-\frac{1}{1/\nu +\omega }}=t^{-0.440(4)},$$

$$D_2\propto t^{-y_6\nu }=t^{1.71(4)},L_2\propto t^{-\nu }=t^{-0.6717(1)}.$$

$${\varphi }_q=L^{y_q}(t{L}^{1/\nu }{)}^{a}g(t{L}^{1/{\nu }_q^{\prime }}),$$

$${\varphi }_q=L^p{t}^{\nu (p-y_q)}g(t{L}^{1/{\nu }_q^{\prime }}).$$

$${\nu }_q^{\prime }=\nu (1-y_q/p)=\nu (1+|y_q|/p),$$

$${\varphi }_q\to 1-k(t{L}^{1/{\nu }_q^{\prime }}),$$

$$1-U\propto (t{L}^{1/\nu }{)}^{-r},$$

$$d_3=\sqrt{L^{-2r/\nu }{t}^{-2r}+L^4{t}^{4{\nu }_q^{\prime }}},$$

$$D_3\propto \sqrt{t^{2r(R-1)}+t^{4({\nu }_q^{\prime }-R\nu )}},L_3\propto t^{-\nu R},$$

$$U=2-\frac{⟨M^4⟩}{⟨M^2{⟩}^2},$$

$$M_x=\sum _{i=1}^N\cos ({\theta }_i),M_y=\sum _{i=1}^N\sin ({\theta }_i).$$

$$L_3\propto t^{-\nu R}=t^{-\nu (r+2{\nu }_q^{\prime })/(r+2\nu )},$$

$$t\propto L_3^{-(r/\nu +2)/(r+2{\nu }_q^{\prime })}.$$

$$t{L}_3^{1/{\nu }_q^{\prime }}$$

$$=$$

$$g(t{L}_3^{1/{\nu }_q^{\prime }})\to g(0)=1,$$

$$m(q,r_i)=\cos (q{\theta }_i),$$

$$C(q,r)=⟨m(q,r_i)m(q,r_j)⟩=⟨\cos (q{\theta }_i-q{\theta }_j)⟩,$$

$$C(q,r_{\text{m}})\sim a{L}^{-2{\Delta }_q}(1+b{L}^{-\omega }),$$

$$A(t,L)=L^{\sigma }g(t{L}^{1/\nu }),$$

Full text

The renormalization group (RG) is a powerful framework both for conceptual understanding of phase transitions and for calculations Wilson71a ; Wilson71b ; Fisher72 . A key concept is that a universal critical point can be stable or unstable in the presence of perturbations, depending on their scaling dimensions. Similarly, an ordered state can also be stable or unstable under the influence of perturbations. Under an RG process, a system flows in a space of couplings which change as the length scale is increased under coarse graining of the microscopic interactions, until finally reaching a fixed point corresponding to a phase or phase transition. At this point, all the initially present irrelevant couplings have decayed to zero.

RG flows can also be defined of physical observables obtained by Monte Carlo (MC) simulations, allowing controlled finite-size scaling analysis—some times referred to as phenomenological renormalization Fisher72 ; Binder81 ; Luck85 ; Wolff09 . Here we extend the standard finite-size scaling of a single observable to an entire flow in a space of two observables associated with relevant or irrelevant couplings. The method is particularly useful for quantifying dangerously irrelevant perturbations (DIPs)—those that are irrelevant at a critical point but become relevant upon coarse graining inside an adjacent ordered phase Amit82 .

Scaling and RG flows.—Consider a $d$-dimensional lattice model of length $L$ which can be tuned to a critical point by a relevant field $t$, e.g., the temperature ($t=T_{c}-T$). With a local operator $m_{i}$ and its conjugate field $h$, we add $h\sum_{i}m_{i}\equiv hM\equiv hL^{d}m$ to the Hamiltonian $H$. In a conventional RG calculation, a flowing field $h^{\prime}$ is computed under a scale transformation. Here we will instead vary the system size, which effectively lowers the energy scale, and calculate the response $\langle m\rangle$ using MC simulations. Together with some quantity $Q$ characterizing the critical point and phases of the system, we can trace out curves (MC RG flows) $(Q,\langle m\rangle)_{L}$ as $L$ increases for fixed values of $h$ and $T$. These flows are very similar to conventional RG flows in the space $(t,h^{\prime})$.

The singular part of the free-energy density takes the form $f_{s}(t,h,L)=L^{-d}F_{s}(tL^{1/\nu},hL^{y})$. At $t=0$, the leading $h$ dependent part is $f_{s}\propto hL^{y-d}$, while the statistical mechanics of $H$ gives a contribution $h\langle m\rangle\propto hL^{-\Delta}$ from the internal energy. Thus, we obtain the well known relation $y=d-\Delta$. The perturbation is irrelevant at the critical point if $y<0$, but, in the case of a DIP, it eventually becomes relevant as $L$ increases in the ordered phase. It has been known for some time that this cross-over is associated with a length scale $\xi^{\prime}\propto t^{-\nu^{\prime}}$ which may diverge faster than the correlation length $\xi\propto|t|^{-\nu}$ Oshikawa00 .

To take both divergent length scales properly into account, i.e., to reach the regime where $tL^{1/\nu^{\prime}}$ is large, we adopt the two-length scaling hypothesis Shao16 and write

[TABLE]

where we have also included a generic scaling correction with exponent $\omega>0$. The exponents $\nu^{\prime}$ and $y$ arise from the same DIP and there is a relationship between them that has been the subject of controversy Oshikawa00 ; Lou07 ; Okubo15 ; Leonard15 . Here we will derive the relationship from Eq. (1) and show how the entire RG flow of two observables can be explained.

Models and observables.—We study three-dimensional (3D) classical clock models on the simple cubic lattice,

[TABLE]

with $\theta\in[0,2\pi)$. Based on previous studies Oshikawa00 ; Lou07 ; Okubo15 ; Leonard15 ; Hove2003 ; Hasenbusch11 ; Pujari15 ; Ding16 , for $q\geq 4$ the phase transition for fixed $h$ at $T=T_{c}$ belongs to the 3D U(1) universality class, i.e., the clock field $h$ is irrelevant. However, for $T<T_{c}$ it is relevant, reducing the order parameter symmetry from U(1) to $Z_{q}$ when observed above the DIP length scale $\xi^{\prime}_{q}$.

In our MC simulations Wolff89 , for a given spin configuration we compute $M_{x}=\sum_{i}\cos(\theta_{i})$ and $M_{y}=\sum_{i}\sin(\theta_{i})$. With $M=(M_{x}^{2}+M_{y}^{2})^{1/2}$ and $\Theta=\arccos(M_{x}/M)$, an angular order parameter can be defined as

[TABLE]

which becomes non-zero in response to the $Z_{q}$ field. This quantity was used to study the length scale $\xi^{\prime}_{q}$ Lou07 ; Pujari15 ; Okubo15 (with a slightly different definition in Refs. Lou07 ; Pujari15 ), but here we will use it in a different way. For $T\geq T_{c}$, $\phi_{q}\to 0$ when $L\to\infty$, while $\phi_{q}\to 1$ for $T<T_{c}$. We will use $\phi_{q}$ in combination with the Binder cumulant $U=2-\langle M^{4}\rangle/\langle M^{2}\rangle^{2}$, which takes the limiting forms $U\to 0$ ($T>T_{c}$), $U\to 1$ ($T<T_{c}$) and $U\to U_{\rm XY}\approx 0.757$ (at $T=T_{c}$ with 3D XY universality Campostrini06 ).

MC RG Flows.—Fig. 1 shows flows of $(U,\phi_{q})_{L}$ for the $q=6$ ”hard” model, i.e., $h\to\infty$ in Eq. (2). Results for $q=4,5$ are discussed in Supplemental Material (SM) sm , where we also determine $T_{c}(h)$ for $q=4,5,6$. The RG process is manifested in the flows with increasing $L$ of the two observables at fixed $T$. The high-$T$ Gaussian fixed point (G) is at $(U,\phi_{q})=(0,0)$; the XY critical point at $(U_{\rm XY},0)$, the U(1) symmetry-breaking Nambu-Goldstone (NG) point at $(1,0)$, and the $Z_{q}$ symmetry-breaking point at $(1,1)$. For $T\geq T_{c}$, we observe simple flows to the fixed points, while for $T<T_{c}$ there are two stages in the flow away from the XY point; first toward the NG point and then an NG to $Z_{q}$ crossover. While this multi-stage flow is expected based on previous RG results Oshikawa00 ; Okubo15 ; Leonard15 , our description with a phenomenological scaling function for accessible observables provides a more practical and intuitive framework for numerical simulations.

Scaling dimensions.—We first study the scaling dimension $y_{q}$ of the $Z_{q}$ field, following the red curve that tends to the XY fixed point in Fig. 1. Previous MC estimates used $Z_{q}$ anisotropy correlators in the pure XY model for $q=4$ Hasenbusch11 . Since the $Z_{q}$ field is irrelevant for $q\geq 4$, the decay power $2\Delta_{q}$ of the correlation function is larger than $6$, which makes it difficult to determine $\Delta_{q}$ accurately (see SM sm for some results). The decay of the induced $\phi_{q}$ is analyzed in Fig. 2 for $q=4,5,6$ at selected $h$ values. The results listed in Table 1 demonstrate that $\phi_{q}$ scales as $M=L^{d}m$ in the general discussion above, i.e., $\phi_{q}\propto L^{-\Delta_{q}+d}=L^{-|y_{q}|}$.

For $q=4$ the $Z_{q}$ field may only be irrelevant for small $h$; the hard model ($h=\infty$) is equivalent to two decoupled Ising models, and for $h=2$ the transition already seems to not be in the XY universality class Pujari15 . Here we use $h=1$. Our simulations extend up to $L=120$ for $q=4$ but smaller for larger $q$ because of the long runs needed to obtain sufficiently small error bars on $\phi_{q}$. To reduce effects of scaling corrections we have excluded small systems until a good fit obtains. Our result $y_{4}=-0.114(2)$ agrees well with the best previous numerical result Hasenbusch11 , but the error bar is smaller. It also matches a high-order nonperturbative expansion Leonard15 . For $q=5$, we have used joint fit to data for several $h$ values, with a common exponent but different prefactors. Our result $y_{5}=-1.27(1)$ is close to an extrapolated value from simulations for smaller $q$ Okubo15 but differs significantly from the field-theory expansions Oshikawa00 ; Leonard15 . For $q=6$ we obtain $y_{6}=-2.55(6)$, which again agrees well with the extrapolated value Okubo15 but differs from those in Refs. Oshikawa00 ; Leonard15 . For all the $q$ values studied, our results show that the first-order $\epsilon$-expansion Oshikawa00 overestimates $y_{6}$, while the nonperturbative expansion Leonard15 underestimates it for $q>4$. All results agree well with a very recent MC calculation of an optimized correlation function Banerjee18 .

Having determined the scaling dimensions, the $Z_{q}$ order parameter in the ordered phase takes the form

[TABLE]

where we neglect the irrelevant arguments in Eq. (1) as they merely produce corrections here. We apply this form to curves such as those shown in Fig. 1, primarily by defining distances to the various fixed points. We study $q=6$ specifically but keep the general-$q$ notation.

Scaling near the XY point.—Though the critical point is well known, it is still useful to study the flows in the two-dimensional space in Fig. 1. We analyze the minimum distances of the $T<T_{c}$ curves to $(U_{\rm XY},0)$. Here $tL^{1/\nu^{\prime}_{q}}\ll tL^{1/\nu}\ll 1$ in Eq. (4), and to leading order

[TABLE]

where we do not include unimportant factors for simplicity. The Binder cumulant scales as

[TABLE]

where $\omega$ is the smallest correction exponent affecting $U$. The scaling form (i.e., without unimportant factors) of the distance $d_{1}$ to the XY fixed point is

[TABLE]

Since $\omega\ll|y_{6}|$, the first term in the square-root dominates; $d_{1}\propto tL^{1/\nu}+L^{-\omega}$, i.e., $d_{1}\to U-U_{XY}$ here (but not necessarily in general). Minimizing for fixed $t$ gives the distance $D_{1}$ and the corresponding system size $L_{1}$

[TABLE]

where we have used $\nu=0.6717(1)$ and $\omega=0.785(20)$ Campostrini06 . Fig. 3(a) shows $d_{1}$ versus $L$ and Fig. 3(b) shows power-law fits to $D_{1}(t)$ and $L_{1}(t)$, where the exponents are $0.372(1)$ and $-0.404(4)$, respectively. These values are in reasonable agreement with Eq. (8) considering scaling corrections for the rather small sizes sm and the neglected subleading $\phi_{6}$ contribution in Eq. (7). The error bars reflect only statistical fluctuations.

Another characteristic of the $T<T_{c}$ curves in Fig. 1 is the minimum distance to the horizontal axis. This RG stage between the XY and NG fixed points is still governed by the XY criticality because $tL^{1/\nu}$ and $tL^{1/\nu^{\prime}_{q}}$ are both small. Since $tL^{1/\nu^{\prime}_{q}}\ll tL^{1/\nu}$, $\phi_{q}$ is given by Eq. (5) and the minimum value $D_{2}$ and corresponding system size therefore scale with $t$ as (for $q=6$)

[TABLE]

The expected exponents indicated above agree reasonably well with our fits in Fig. 4, where the exponents are $1.88(2)$ and $-0.60(3)$, respectively. The deviations are again likely due to scaling corrections.

*Cross-over exponent $\nu^{\prime}_{q}$.—*When $tL^{1/\nu}\gg 1$ but $tL^{1/\nu^{\prime}_{q}}$ is arbitrary, Eq. (4) must reduce to

[TABLE]

where the exponent $a$ follows from the physics of the clock model. Specifically, we can ask how $\phi_{q}$ depends on $L$ at fixed $t$ when the $U(1)$ symmetry is barely broken down to $Z_{q}$, i.e., when $\phi_{q}\ll 1$. This is a subtle issue at the heart of the long-standing controversy regarding the symmetry cross-over Ueno91 ; Oshikawa00 ; Lou07 ; Okubo15 ; Leonard15 . Instead of invoking physical arguments, we will here simply posit that $\phi_{q}\propto L^{p}$ in the regime where $tL^{1/\nu}$ is large but $tL^{1/\nu^{\prime}_{q}}$ remains small [hence $g\approx 1$ in Eq. (10)], and later show how $p$ can be consistently determined from the MC RG flows. Thus, we have $a=\nu(p-y_{q})$ in Eq. (10);

[TABLE]

This form should apply also when $\phi_{q}\to 1$, demanding $g\to(tL^{1/\nu^{\prime}_{q}})^{b}$ with $b=-\nu(p-y_{q})$ and $\nu^{\prime}_{q}=-b/p$. Then

[TABLE]

which for $p=3$ agrees with Ref. Lou07 , while for $p=2$ it agrees with Refs. Okubo15 ; Leonard15 . When $\phi_{q}$ deviates from $1$, $g\to(tL^{1/\nu^{\prime}_{q}})^{b}[1-k(tL^{1/\nu^{\prime}_{q}})]$, so that for large $tL^{1/\nu^{\prime}_{q}}$

[TABLE]

where the function $k$ must be dimensionless.

The exponent $\nu^{\prime}_{q}$ in Eq. (13) can be determined by a standard data-collapse procedure Okubo15 ; Lou07 . Here we proceed in a different way: The function $k(x)$ can be Taylor expanded around some arbitrary point $x_{0}$ where $\phi_{q}=y_{0}$; $\phi_{q}=y_{0}+a(x-x_{0})$, or $\phi_{q}=ax+b$ for some $b$. For fixed $t$, we consider $L=L_{c}$ for which $\phi_{q}(L_{c})=e$ for some $e$, which gives $L_{c}\propto t^{-\nu^{\prime}_{q}}$. In Fig. 5(a) we extract $L_{c}$ for $e=0.5,0.55$, and $0.6$. Analyzing the scaling behavior with $t$ in Fig. 5(b), we find $\nu^{\prime}_{6}=1.52(4)$. Thus, Eq. (12) with $|y_{6}|=2.55(6)$ is satisfied if $p=2$, in agreement with Refs. Okubo15 ; Leonard15 . From Eq. (11), the initial growth of $\phi_{q}$ with $L$ is then $\phi_{q}\propto L^{2}$; not $\propto L^{3}$ Lou07 .

Near the NG fixed point.—Finally we consider the distance to the NG fixed point $(1,0)$, where Eq. (11) applies with $g\approx 1$ ($L\ll\xi^{\prime}_{q}$ can be tested self-consistently sm ). $U$ is close to $1$, but should remain of the form $U(tL^{1/\nu})$ because, as we will see, $L$ and $t$ for a given curve in the region of interest are related such that $t\to 0$ when $L\to\infty$. We need $1-U$, which has a non-trivial scaling form

[TABLE]

where it has been argued that, in some cases, $r=d\nu=3\nu$ Privman84 . However, this result is based on subtle assumptions and may not be generic Privman90 . As shown in SM sm , $r=1.52(2)\not=3\nu$ for the XY model.

The distance to the NG fixed point is, from Eq. (14) and Eq. (11) with $\nu(2-y_{q})=2\nu^{\prime}_{q}$ and $g\approx 1$;

[TABLE]

and minimizing with respect to $L$ leads to

[TABLE]

where $R=(r+2\nu^{\prime}_{q})/(r+2\nu)$. For the $q=6$ case we then have $D_{3}\propto t^{0.9(1)}$ and $L_{3}\propto t^{-1.07(3)}$. From the analysis in Fig. 6 the exponents are $1.19(3)$ and $-1.14(2)$, respectively, in reasonable agreement with the prediction, again considering that we have not included any scaling corrections. The cross-over behavior around the NG point is also the most intricate of all the regions in the way the two length scales intermingle.

Discussion.—The standard finite-size scaling hypothesis in the presence of a DIP (see, e.g., Ref. Kenna13 ) includes only $tL^{1/\nu}$ and the irrelevant field $hL^{y}$ in Eq. (1), which is sufficient for extracting the critical exponents close to $T_{c}$; up to $|T-T_{c}|\propto L^{-1/\nu}$. As we have shown here with the clock model, the other relevant variable $tL^{1/\nu^{\prime}_{q}}$ is necessary for describing the symmetry cross-over from U(1) to $Z_{q}$. By considering different necessary (for scaling) limiting forms when the arguments are small or large, we have quantitatively explained the entire MC RG flows.

The controversial relationship between $\nu^{\prime}_{q}$ and the scaling dimension $y_{q}$ Ueno91 ; Oshikawa00 ; Lou07 ; Okubo15 ; Leonard15 involves an exponent $p$ associated with the initial formation of an effective $Z_{q}$ symmetric potential for the order parameter. Analytical RG methods for related problems, e.g., the Sine-Gordon model with a weak potential are indeed highly non-trivial and sensitive to the type of approximation used Nandori04 . In our approach, $p$ for a given system is obtained from numerical data and can then be used to further understanding of the subtle physics of the DIP. We have here confirmed numerically that $p=2$ in the clock model Okubo15 ; Leonard15 , but this exponent is not necessarily universal—it may depend on a combination of the finite-size properties of the fixed point with the higher symmetry (here the well-understood NG point Hasenfratz90 ; Dimitrovic91 ) and the mechanisms of the DIP causing the lowering of the symmetry.

Our method should be useful in the context of deconfined quantum criticality Senthil04 ; Sandvik07 ; Lou09 , where a scaling ansatz with two relevant arguments was introduced to account for anomalous scaling in 2D quantum magnets Shao16 . There the DIP cannot be tuned away (unlike some fermionic models Liu18 ), because it is connected to the lattice itself. Thus, the method introduced here of studying scaling and RG flows in the presence of a finite DIP is ideal.

I Supplementary Information

Monte Carlo Renormalization Flows in the Space of Relevant and Irrelevant Operators: Application to Three-Dimensional Clock Models

Hui Shao,1,2,∗ Wenan Guo,3,2,† Anders W. Sandvik,4,5,3,‡

1*Center for Advanced Quantum Studies, Department of Physics,

Beijing Normal University, Beijing 100875, China

2Beijing Computational Science Research Center, Beijing 100193, China

3Department of Physics, Beijing Normal University, Beijing 100875, China

4Department of Physics, Boston University, 590 Commonwealth Avenue, Boston, Massachusetts 02215, USA* 5Beijing National Laboratory for Condensed Matter Physics and Institute of Physics,

Chinese Academy of Sciences, Beijing 100190, China

e-mail: ∗huishao@bnu.edu.cn, †waguo@bnu.edu.cn, ‡sandvik@bu.edu

We discuss further results that were used in the main text. In Sec. 1 we determine $T_{c}$ for the $q=4,5$ and $6$ clock models. In Sec. 2 we show MC RG flow diagrams for $q=4$ and $5$, complementing the $q=6$ results in Fig. 1 in the main paper. In Sec. 3 we determine the scaling dimensions of the $Z_{q}$ perturbations using the conventional correlation-function method for $q=1$-$4$. In Sec. 4 we determine the exponent $r$ governing the asymptotic form of the Binder cumulant $U(x)$ in Eq. (14), by MC calculations for large values of $x=tL^{1/\nu}$ in the ordered phase.

I.1 1. Determination of critical temperatures

To extract the critical temperatures for the clock models with different $q$ and $h$, we calculate the Binder cumulant of the two-component vector order parameter,

[TABLE]

with $M=\sqrt{M_{x}^{2}+M_{y}^{2}}$, where

[TABLE]

In a standard $(L,2L)$ crossing-point analysis Luck85 (described in detail and tested, e.g., in the Supplemental Information of Ref. Shao16 ), we have computed the cumulant for a series of system sizes around the critical point in each model and used cubic polynomials to interpolate and extract the crossing points defining the flowing critical temperature $T^{*}(L,2L)$ and the associated cumulant value $U^{*}(L,2L)$. In Fig. S1(a,b) we analyze the size dependence of these quantities for all the models studied in the main text. The infinite-size extrapolated $T_{c}$ values are summarized in Table. SI. We have also tested the consistency of the critical exponent $\nu$ of the correlation length (obtained from the derivatives $dU/dT$ at the crossing points) and the universal value of the Binder cumulant $U_{c}$ with the 3D O(2) universality class Campostrini06 ; the results with increasing $L$ tend to values fully consistent with the known numbers, as shown in Fig. S1(b) and (c), though the error bars of the $1/\nu$ estimates are large.

For the $q=6$ case, we also present results for the exponents of the scaling corrections in Fig. S1(a) and (b). In (a), the fit to a power-law correction gives $1/\nu+\omega=2.45(3)$ and in (b) we similarily find $\omega=0.94(3)$ from the correction to $U_{c}$. These results do not agree fully with the known values $1/\nu=1.4890(6)$ and $\omega=0.78(2)$ Campostrini06 , but if we fix $1/\nu$ to its known value in the estimate for $1/\nu+\omega$, then the value of $\omega$ is statistically consistent with the value from the $U^{*}$ fit. Since we have only included one correction here, and influence from the higher-order corrections may be significant still at these system sizes, the exponent $\omega$ should be considered as an “effective exponent”, whose value should approach the true value for system sizes larger than those used here.

I.2 2. MC RG Flows for the $q=4,5$ clock models.

In addition to the $q=6$ MC RG flows discussed in the main paper, we have performed more limited simulations of the cases $q=4$ and $5$. Results for the $q=5$ hard-constrained model is shown in Fig. S2, with data distributed mostly near the XY and NG fixed points. For the system sizes available, there is no $T$ for which we can observe both the flow toward the NG fixed point and the cross-over away from this point toward the $Z_{5}$ fixed point. However, we can see these parts of the flows separately for suitably chosen temperatures (the two groups of curves in Fig. S2). On a qualitative level the flows are very similar to the $q=6$ case.

Figure S3 shows results for the case $q=4$, $h=1$. At first sight, the flows here appear to be very different from the $q=5,6$ cases. However, this should just be due to the small scaling dimension of the $Z_{4}$ field, $|y_{4}|\approx 0.11$ (Table I in the main paper). This means that $\phi_{4}$ decays very slowly with increasing system size, as is clear both from Fig. 2 in the main paper and the red set of $T_{c}$ data in Fig. S3. For the system sizes available, there are not yet any sign of flows toward the NG point before the ultimate flow toward the $Z_{4}$ fixed point. For very large system sizes we expect that such cross-over behavior should be manifested also in this case, but to observe it requires a clear separation of the length scales $\xi\propto t^{-\nu}$ and $\xi^{\prime}_{q}\propto t^{-\nu^{\prime}_{q}}$. Since in this case the difference between the exponents is very small, $\nu^{\prime}_{4}-\nu=\nu|y_{4}|/2\approx 0.04$, if we would like to have, say, $\xi^{\prime}/\xi=10$, we need $t\approx 10^{-25}$ (assuming all proportionality factors are of order one). From our analysis of the flow away from the NG fixed point, summarized as Eq. (16), we then have roughly $L_{c}\propto t^{-\nu R}\approx 10^{18}$ for the system size where the cross-over will occur. This length scale is clearly beyond any current or future MC calculations.

It is also important to check whether it is always true that $tL^{1/\nu^{\prime}}$ asymptotically vanishes when the cross-over in the neighborhood of the NG point takes place. This was the assumption under which we derived the cross-over point with minimum distance $D_{3}$ to the NG point and the associated length $L_{3}$, because we set $g=1$ in the scaling form Eq. (10) of $\phi_{q}$. Rewriting the scaling form of $L_{3}$ in Eq. (16) as

[TABLE]

we have that $t$ scales with $L_{3}$ as

[TABLE]

Thus, the relevant scaling argument corresponding to the second length scale depends on $L_{3}$ as

[TABLE]

where that the exponent on $L_{3}$ is always negative because $\nu^{\prime}_{q}>\nu$ for a DIP. Therefore,

[TABLE]

and the self-consistency of the assumption is confirmed for any $q\geq 4$ in the neighborhood of the NG cross-over.

I.3 3. Scaling dimensions $y_{q}$ from correlation functions in the XY model

The standard way to obtain the scaling dimension of an irrelevant or relevant operator is to compute the related correlation function at the critical point in the model without the perturbation. In the case of the 3D XY model, the best MC calculation of the scaling dimension of the $Z_{4}$ clock perturbation is in Ref. Hasenbusch11 . Because of the rapid decay of the correlation functions for larger $q$, no MC results based on the conventional method are available for $q>4$, as far as we are aware. Our method presented in the main paper can reach larger $q$ because of the slower decay of the induced operator expectation value in the presence of the perturbation.

Here we contrast the conventional and new method by considering the $q=4$ case, computing the $Z_{q}$ correlator with MC simulations at the 3D XY critical point, using $T_{c}=2.20184$ Campostrini06 . The local operator corresponding to the $Z_{q}$ field can be taken as:

[TABLE]

and we study the corresponding correlation function

[TABLE]

where $r=r_{i}-r_{j}$ and the global rotational symmetry has been taken into consideration.

In Fig. S4 we analyze the long-distance correlation function $C(q,r_{\text{m}})$ in the three different lattice directions (i.e., $r_{\text{m}}$ is half the system length in the respective directions), as indicated in the inset of the figure. The asymptotic form should be

[TABLE]

where $\Delta_{q}=3-y_{q}$, with $y_{q}$ being the scaling dimension of the $Z_{q}$ field, and we have also included a scaling correction with exponent $\omega$. We perform joint fit to Eq. (S10) with the MC data along all three directions, where same exponents but different prefactors $a$ are used.

In Fig. S4 we present results for $q=1,2,3$, i.e., the cases in which the $Z_{q}$ fields are relevant. The results for the scaling dimensions are summarized in Table SII and compared with previous MC studies Campostrini06 ; Hasenbusch11 . The agreement is good, and in the case of $q=2$ we improve on the statistical error. We should note here that the previous study used a system-volume integrated correlator, for which the statistical errors of the correlations are smaller but the corrections may be larger.

When the $Z_{q}$ field becomes irrelevant, the decay exponent of the correlation function grows larger than $6$, and it becomes extremely hard to extract the scaling dimension in this way. We show our $q=4$ data in Fig. S5. Here we do not report any results of fitting, but only indicate the expected decay power $2\Delta=(3+|y_{4}|)\approx 6.23$ based on the scaling dimension $y_{q}\approx-0.114$ extracted with the alternative method in the main paper.

Here we should again note that the previous MC study Hasenbusch11 used a system-integrated correlator, for which the decay exponent is $2(3-\Delta_{q})=2y_{q}$. With the larger exponent due to summation over the system volume, the error bars are significantly reduced and the results were therefore considerably less noisy than in the data presented here. The long-distance correlator is possibly less affected by scaling corrections, though we have not tested this. Our approach of explicitly including the field still appears to work better, having a decay exponent of just $y_{q}$. Our main purpose of studying the $Z_{q}$ correlation functions here was mainly to establish the consistency between the two approaches.

I.4 4. Asymptotic form of the Binder cumulant

Recall that, in the critical finite-size scaling form of some singular quantity $A$,

[TABLE]

the exponent $\sigma$ must be compatible with the asymptotic form of the scaling function $g(x)$, $x=tL^{1/\nu}$. This behavior is connected to the size-independent scaling form in the thermodynamic limit, $A\propto t^{\kappa}$ (where $\kappa$ is a generic notation for the critical exponent for the quantity in question), which is obtained if $g\to x^{\kappa}$ when $x\to\infty$ (i.e., $L\to\infty$ for fixed small $t$). Then, to eliminate the $L$ dependence we must have $\sigma=-\kappa/\nu$.

In the case of the dimensionless Binder cumulant $U$, $\sigma=0$ and, accordingly, the corresponding scaling function $g(x)$ in Eq. (S11) must take the form $g\to c$, where $c$ is a constant which we know takes the value $c=1$ in the ordered phase (while $c=0$ in the disordered phase). The scaling form does not immediately tell us how $g$ approaches $1$, however, which is what we need in the analysis of the flow close to the NG fixed point in the main paper. It should be noted that the scaling regime of interest here does not yet correspond to Gaussian fluctuations in the ordered phase, because $t$ approaches zero with increasing length-scale, as shown in the main paper. A natural assumption is that $1-g(x)$ takes a power-law form, $1-g(x)\propto x^{-r}$, corresponding to the form of $1-U$ in Eq. (14). The exponent $r$ should presumably also be related to the critical exponents of the universality class in question.

Surprisingly, while the Binder cumulant is one of the most important quantities used to characterize critical points in numerical studies Binder81 ; Luck85 , the asymptotic form of $1-U$ has not been extensively studied—the focus has naturally been on the behavior for small arguments; $x=0$ and $x\approx 0$. We are only aware of Privman’s work on the asymptotic $x\to\infty$ behavior Privman84 ; Privman90 . He argued that $r=d\nu$ but also pointed out that the assumptions underlying this conclusion are somewhat speculative and untested.

To investigate the scaling behavior, we have carried out systematic MC calculations of the 3D XY model and the $q=6$ clock model inside their ordered phase in order to extract the exponent $r$ independently. Our results for the XY model are shown in Fig. S6(a). We performed dedicated simulations targeting $1-U$ for system sizes up to $L=256$ for a wide range of the scaling variable $tL^{1/\nu}$, sufficient to reliably observe data collapse and an asymptotic power-law form. A fit to the $L=256$ data gives the exponent $r=1.52(2)$, which is clearly different from Privman’s prediction $r=3\nu\approx 2.02$ Privman84 ; Privman90 . As Privman pointed out, there are subtle assumptions made in the derivation of his result, and the behavior may not be generic. In the case here, the exponent is consistent within statistical errors with the exponent $1/\nu$, but we see no obvious reason for this value.

In the case of the clock model, results for which are shown in Fig. S6(b), we have just plotted the same data that we used in the main paper, going up only to $L=64$. The data forming a group in the range $tL^{1/\nu}\approx 50\sim 100$ are fully consistent with the same exponent as in the XY model, and for lower values of the scaling variable the behaviors are also very similar. For small and moderate values of $tL^{1/\nu}$ it is clear that the clock and XY models should behave very similarly in this regard, since the clock field close to $T_{c}$ is irrelevant. However, when $tL^{1/\nu}$ is larger, e.g., when $tL^{1/\nu}\approx 100$ in in Fig. S6(b), there could in principle be a cross-over behavior also in $U$, where $tL^{1/\nu^{\prime}_{q}}$ may impact the scaling behavior (perhaps as a correction) when it also reaches large values. We do not see any evidence of a break-down of the $tL^{1/\nu}$ scaling, however.

It would be interesting to study $1-U$ also for other models, to test the generality of the results found here.