# Quantum criticality with two length scales

\clearpage\pagenumberingarabic\clearpage\pagenumberingarabic

The theory of âdeconfinedâ quantum critical points describes phase transitions at temperature outside the standard paradigm, predicting continuous transformations between certain ordered states where conventional theory requires discontinuities. Numerous computer simulations have offered no proof of such transitions, however, instead finding deviations from expected scaling relations that were neither predicted by the DQC theory nor conform to standard scenarios. Here we show that this enigma can be resolved by introducing a critical scaling form with two divergent length scales. Simulations of a quantum magnet with antiferromagnetic and dimerized ground states confirm the form, proving a continuous transition with deconfined excitations and also explaining anomalous scaling at . Our findings revise prevailing paradigms for quantum criticality, with potentially far-reaching implications for many strongly-correlated materials.

#### Introduction

In analogy with classical phase transitions driven by thermal fluctuations, condensed matter systems can undergo drastic changes as parameters regulating quantum fluctuations are tuned at low temperatures. Some of these quantum phase transitions can be theoretically understood as rather straight-forward generalizations of thermal phase transitions [1, 2], where, in the conventional Landau-Ginzburg-Wilson (LGW) paradigm, states of matter are characterized by order parameters. Many strongly-correlated quantum materials seem to defy such a description, however, and call for new ideas.

A promising proposal is the theory of deconfined quantum critical (DQC) points in certain two-dimensional (2D) quantum magnets [3, 4], where the order parameters of the antiferromagnetic (Néel) state and the competing dimerized state (the valence-bond-solid, VBS) are not fundamental variables but composites of fractional degrees of freedom carrying spin . These spinons are condensed and confined, respectively, in the Néel and VBS state, and become deconfined at the DQC point separating the two states. Establishing the applicability of the still controversial DQC scenario would be of great interest in condensed matter physics, where it may play an important role in strongly-correlated systems such as the cuprate superconductors [5]. There are also intriguing DQC analogues to quark confinement and other aspects of high-energy physics, e.g., an emergent gauge field and the Higgs mechanism and boson [6].

The DQC theory represents the culmination of a large body of field-theoretic works on VBS states and quantum phase transitions out of the Néel state [7, 8, 9, 2, 10]. The postulated SU() symmetric non-compact (NC) CP action can be solved when [11, 5, 12] but non-perturbative numerical simulations are required to study small . The most natural physical realizations of the Néel–VBS transition for electronic SU(2) spins are frustrated quantum magnets [9], which, however, are notoriously difficult to study numerically [13, 14]. Other models were therefore pursued. In the - model [15], the Heisenberg exchange between spins is supplemented by a VBS-inducing four-spin term which is amenable to efficient quantum Monte Carlo (QMC) simulations [15, 16, 17, 18, 19, 20, 21, 22, 23]. Although many results for the - model support the DQC scenario, it has not been possible to draw definite conclusions because of violations of expected scaling relations that affect many properties. Similar anomalies were later observed in three-dimensional loop [24] and dimer [25] models, which are also potential realizations of the DQC point. Simulations of the NCCP action as well have been hard to reconcile with the theory [26, 27, 21].

One interpretation of the unusual scaling behaviors is that the transitions are first-order, as generally required within the LGW framework for order–order transitions where unrelated symmetries are broken. The DQC theory would then not apply to any of the systems studied so far, thus casting doubts on the entire concept [17, 26, 21]. In other interpretations the transition is continuous but unknown mechanisms cause strong corrections to scaling [18, 27, 28] or modify the scaling more fundamentally in some yet unexplained way [19, 24]. The enigmatic current state of affairs is well summed up in the recent Ref. [24].

Here we show that the DQC puzzle can be resolved based on a finite-size scaling Ansatz including the two divergent length scales of the theory—the standard correlation length , which captures the growth of both order parameters (), and a faster-diverging length associated with the thickness of VBS domain walls and spinon confinement (the size of a spinon bound state). We show that, contrary to past assumptions, can govern the finite-size scaling even of magnetic properties that are sensitive only to in the thermodynamic limit. Our simulations of the - model at low temperatures and in the lowest (two-spinon) excited state demonstrate complete agreement with the two-length scaling hypothesis, with no other anomalous scaling corrections remaining.

#### Finite-size scaling forms

Consider first a system with a single divergent correlation length , where is the distance to a phase transition driven by quantum fluctuations arising from non-commuting interactions controlled by at and is the critical value of . In finite-size scaling theory [29], for a system of linear size (volume in dimensions), close to a singular quantity takes the form

(1) |

where the exponents are tied to the universality class, also depends on , and approaches a constant when (up to scaling corrections ). We assume (or, alternatively, ) so that scaling arguments depending on have been eliminated.

The form (1) fails for some properties of the - model [18, 19, 22] and other DQC candidate systems [26, 24, 25]. A prominent example is the spin stiffness, which for an infinite 2D system in the Néel phase should scale as with [3, 4, 1]. To eliminate the size dependence when and in Eq. (1), we must have and for large . Thus, and should be constant when if . However, at criticality instead appears to diverge slowly [17, 18, 21]. At first sight this might suggest , but other quantities, e.g., the magnetic susceptibility, instead behave as if [30]. Strong scaling corrections have been suggested as a way out of this paradox [18, 19, 28]. Claims of a weak first-order transition have also persisted [26, 27, 21], though the continuous DQC scenario is supported by the absence of any of the usual first-order signals, e.g., the Binder cumulant does not exhibit any negative peak [18, 24].

To explain the scaling anomalies phenomenologically, in the presence of a second length in the VBS, we propose that Eq (1) should be replaced by the form

(2) |

where, unlike what was assumed in the past, is not necessarily the same as the exponent which governs the behavior of most observables in the thermodynamic limit. Instead, we show that the criticality in the - model generically demands .

First assume . The correct thermodynamic limit with for can then be obtained from Eq. (2) if for large , and, as before, . This behavior can also be expressed using a scaling function where the second argument is the ratio of the two lengths; . If is constant when , then acts like just another irrelevant field, as in the standard scenario for dangerously irrelevant perturbations in classical clock models [31]. Our proposal is a different large- limit of Eq. (2) controlled by , which leads to concrete predictions of scaling anomalies. In the case of the stiffness, the correct thermodynamic limit is obtained with and if for large . Then , which we can also obtain with and for . A function behaving as a power of for was implicitly suggested in Ref. [19], though with no specific form.

This alternative scaling behavior corresponds to saturating at when upon approaching the critical point, in contrast to the standard scenario where grows until it also reaches [32]. The criticality at distances is conventional, whereas is governed by the unconventional power laws. Different behaviors for and were actually observed in the recent loop-model study [24] and a dangerously irrelevant field was proposed as a possible explanation, but with no quantitative predictions of the kind offered by our approach. The anomalous scaling law controlled by , which we confirm numerically below, is an unexpected feature of DQC physics and may also apply to other systems with two divergent lengths.

#### Quantum Monte Carlo Results

The - model [15] for spins is defined using singlet projectors as

(3) |

where denotes nearest-neighbor sites on a periodic square lattice with sites and the pairs and in form horizontal and vertical edges of plaquettes. This Hamiltonian has all symmetries of the square lattice and the VBS ground state existing for (with ) is columnar, breaking the translational and rotational symmetries spontaneously. The Néel state for breaks the spin-rotation symmetry.

We will study several quantities in the neighborhood of . Although we have argued that the asymptotic behavior when in Eq. (2) is controlled by the second argument of , the critical finite-size scaling close to (when is of order ) can still be governed by the first argument [32]. We will demonstrate that, depending on the quantity, either or is the relevant argument, and, therefore, and can be extracted using finite-size scaling with effectively single-parameter forms. We will do this for manifestly dimensionless quantities, in Eq. (2), before testing the anomalous powers of in other quantities.

If the effective one-parameter scaling holds close to , then Eq. (2) implies that at some and a crossing-point analysis (Fisher’s phenomenological renormalization) can be performed [29]. For a quantity, if and with constant, a Taylor expansion of shows that the crossing points approach as if is the relevant exponent (which we assume here for definiteness). approaches its limit as , and one can also show that

(4) |

converges to at the rate . In practice, simulation data can be generated on a grid of points close to the crossing values, with polynomials used for interpolation and derivatives. We present details and tests of such a scheme for the Ising model in Supplementary Material.

In the sector, spinon physics can be studied with projector QMC simulations in a basis of valence bonds (singlet pairs) and two unpaired spins [33, 34]. Previously the size of the spinon bound state in the - model was extrapolated to the thermodynamic limit [35], but the results were inconclusive as to the rate of divergence upon approaching the critical point. Here we study the critical finite-size behavior. We define the size of the spinon pair using the strings connecting the unpaired spins in valence-bond QMC simulations [33, 34], as illustrated in Fig. 1 and further discussed in Supplementary Material.

If when , then follows from our proposed limit of Eq. (2). If manifestly probes only the longer length scale also in a finite system, which we will test, then is the exponent controlling the crossing points of . Data and fits are presented in Fig. 2 (left). Unlike other quantities used previously to extract the critical point [18], the drift of with is monotonic in this case and the convergence is rapid. All points are consistent with the expected power-law correction, with and , where the number in parenthesis indicates the statistical uncertainty (one standard deviation) in the preceding digit. The critical point agrees well with earlier estimates [18]. The crossing value also clearly converges and a slope analysis according to Eq. (4) gives .

Next, in Fig. 2 (right) we analyze a Binder ratio, defined with the -component of the sublattice magnetization as and computed at as in Ref. [18]. Here the non-monotonic behavior of the crossing points necessitates several scaling corrections, unless only the largest sizes are used. In either case, the behavior of is fully consistent with obtained from . has an uncertainty of over because of the small value of the correction exponent; . The slope estimator (4) of the exponent is monotonic and requires only a single correction, also with a small exponent . The extrapolated exponent is close to the value obtained recently for the loop model [24].

The above results support a non-trivial deconfinement process where the size of the bound state diverges faster than the conventional correlation length. However, in the DQC theory the fundamental longer length scale is the thickness of a VBS domain wall. It can be extracted from the domain wall energy per unit length, , which in the thermodynamic limit should scale as [4]. In Supplementary Material we re-derive this form using a two-length scaling Ansatz and discuss simulations of domain walls in a 3D clock model and the - model. At criticality, in the conventional scenario (exemplified by the clock model) both and saturate at and . For the - model we instead find with for large , as illustrated in Fig. 3(a). Our interpretation of this unconventional scaling is that, when saturates at , also stops growing and remains at . Thus with , which agrees reasonably well with obtained from the quantities in Fig. 2. The large error bar on the latter ratio leaves open the possibility that the spinon confinement exponent is between and the domain-wall exponent [4].

We also calculated the critical spin stiffness and susceptibility for the smallest wave-number , using . Conventional quantum-critical scaling dictates that these quantities should decay as when [2]. Instead, Figs. 3(b) and (c) demonstrate clearly slower decays, and being slowly divergent, as had been found in earlier works as well[17, 18, 30, 19, 21]. The unconventional limit of the scaling function (2) requires and to diverge as . The behaviors are indeed consistent with this form and extracted from , with a correction with a small (close to the correction for in Fig. 2). The mutually consistent scaling behavior of the three different quantities lends strong support to the new type of criticality where the magnetic properties are not decoupled from the longer VBS length scale for finite . The results are incompatible with a first-order transition, where, , , .

#### Discussion

We have shown that the effects of the larger divergent length scale at the Néel–VBS transition are more dramatic than those caused by standard [31] dangerously-irrelevant perturbations, and we therefore propose the term super dangerous for this case. The universality class, in the sense of the normal critical exponents in the thermodynamic limit at , are not affected by such perturbations, but anomalous power laws of the system size appear generically in finite-size scaling. We have determined the value for the exponent ratio governing the anomalous scaling in the - spin model.

Loop and dimer models exhibit similar scaling anomalies [24, 25] and it would be interesting to test the consistency between different quantities as we have done here. In simulations of the NCCP action [21, 26, 27] one would at first sight not expect any effects related to the longer DQC length scale, because the monopoles responsible for the VBS condensation are not present in the continuum theory [3]. There could still be some other super dangerous operator present [see also Ref. [24]], perhaps related to lattice regularization.

The consequences of our findings extend also to quantum criticality in the thermodynamic limit, because is the thickness of an equivalent system in the path integral formulation [1, 2]. Anomalous finite- behaviors of the - model have already been observed [18, 30]. For instance, the spin correlation length at , which should be affected by deconfined spinons, grows faster than the normally expected form and the susceptibility vanishes slower than . Remarkably, the asymptotic forms and can account for the respective behaviors (Fig. S10 [32]). Thus, we find a strong rationale to revise the experimentally most important tenet of quantum criticality—the way scaling at the critical point is related to power laws in at .

We conclude that quantum criticality with two divergent length scales is much richer than previously anticipated. Our findings may apply to a wide range of strongly-correlated quantum systems with more than one length scale and may help to resolve the mysteries still surrounding scaling behaviors in materials such as the high- cuprate superconductors.

#### Acknowledgments

We would like to thank Kedar Damle for suggesting a definition of spinons using valence-bond strings [33, 34]. The research was supported by NSFC Grant No. 11175018 (WG), the Fundamental Research Funds for the Central Universities (WG), NSF Grant No. DMR-1410126 (AWS), and by the Simons Foundation (AWS). HS and WG gratefully acknowledge support from the Condensed Matter Theory Visitors Program at Boston University and AWS is grateful to the Institute of Physics, Chinese Academy of Sciences, and Beijing Normal University for their support. Some of the computations were carried out using Boston University’s Shared Computing Cluster.

## References and Notes

- [1] M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. S. Fisher, Phys. Rev. B 40, 546 (1989).
- [2] A. V. Chubukov, S. Sachdev, and J. Ye, Phys. Rev. B 49, 11919 (1994).
- [3] T. Senthil, A. Vishwanath, L. Balents, S. Sachdev, and M. P. A. Fisher, Science 303, 1490 (2004).
- [4] T. Senthil, L. Balents, S. Sachdev, A. Vishwanath, and M. P. A. Fisher, Phys. Rev. B 70, 144407 (2004).
- [5] R. K. Kaul and S. Sachdev, Phys. Rev. B 77, 155105 (2008).
- [6] Y. Huh, P. Strack, and S. Sachdev, Phys. Rev. Lett. 111, 166401 (2013).
- [7] N. Read and S. Sachdev, Phys. Rev. B 42, 4568 (1990).
- [8] G. Murthy and S. Sachdev, Nucl. Phys. B 344, 557 (1990).
- [9] N. Read and S. Sachdev, Phys. Rev. Lett. 66, 1773 (1991).
- [10] O. I. Motrunich and A. Vishwanath, Phys. Rev. B 70, 075104 (2004).
- [11] M. A. Metlitski, M. Hermele, T. Senthil and M. P. A. Fisher, Phys. Rev. B 78, 214418 (2008).
- [12] E. Dyer, M. Mezei, S. S. Pufu, and S. Sachdev, J. High Energy Phys. 2015, 37 (2015).
- [13] T. Li, F. Becca, W. Hu, and S. Sorella, Phys. Rev. B 86, 075111 (2012).
- [14] S.-S. Gong, W. Zhu, D. N. Sheng, O. I. Motrunich, and M. P. A. Fisher, Phys. Rev. Lett. 113. 027201 (2014).
- [15] A. W. Sandvik, Phys. Rev. Lett. 98, 227202 (2007).
- [16] R. G. Melko and R. K. Kaul, Phys. Rev. Lett. 100, 017203 (2008).
- [17] F. J. Jiang, M. Nyfeler, S. Chandrasekharan, and U. J. Wiese, J. Stat. Mech. (2008), P02009.
- [18] A. W. Sandvik, Phys. Rev. Lett. 104, 177201 (2010).
- [19] R. K. Kaul, Phys. Rev. B 84, 054407 (2011).
- [20] K. Harada, T. Suzuki, T. Okubo, H. Matsuo, J. Lou, H. Watanabe, S. Todo, and N. Kawashima, Phys. Rev. B. 88, 220408 (2013).
- [21] K. Chen. Y. Huang, Y. Deng, A. B. Kuklov, N. V. Prokof’ev, and B. V. Svistunov, Phys. Rev. Lett. 110, 185701 (2013).
- [22] M. S. Block, R. G. Melko, and R. K. Kaul, Phys. Rev. Lett. 111, 137202 (2013).
- [23] S. Pujari, F. Alet, and K. Damle, Phys. Rev. B 91, 104411 (2015).
- [24] A. Nahum, J. T. Chalker, P. Serna, M. Ortunõ, and A. M. Somoza, Phys. Rev. X 5, 041048 (2015).
- [25] G. J. Sreejith and S. Powell, Phys. Rev. B 89, 014404 (2014).
- [26] A. B. Kuklov, M. Matsumoto, N. V. Prokof’ev, B. V. Svistunov, and M. Troyer, Phys. Rev. Lett. 101, 050405 (2008).
- [27] O. I. Motrunich and A. Vishwanath, arXiv:0805.1494 (2008).
- [28] L. Bartosch, Phys. Rev. B 88, 195140 (2013).
- [29] J. M. Luck, Phys. Rev. B 31, 3069 (1985).
- [30] A. W. Sandvik, V. N. Kotov, and O. P. Sushkov, Phys. Rev. Let. 106, 207203 (2011).
- [31] F. Léonard and B. Delamotte, Phys. Rev. Lett 115, 200601 (2015).
- [32] See Supplementary Material.
- [33] Y. Tang and A. W. Sandvik, Phys. Rev. Lett. 107, 157201 (2011).
- [34] A. Banerjee and K. Damle, J. Stat. Mech. (2010) P08017.
- [35] Y. Tang and A. W. Sandvik, Phys. Rev. Lett. 110, 217213 (2013).
- [36] Contributions originating from the regular part of the free energy are also of the same form and in some cases, where the exponents of the irrelevant fields are large, they can cause the leading scaling corrections [37].
- [37] G. Kamieniarzt and H. W. J. Blöte, J. Phys. A: Math. Gen. 26, 201 (1993).
- [38] E. Marinari, Optimized Monte Carlo Methods, Advances in Computer Simulation, edited by J. Kertesz and I. Kondor (Springer, Berlin, 1997); arXiv:cond-mat/9612010.
- [39] A. W. Sandvik and H. G. Evertz, Phys. Rev. B. 82, 024407 (2010).
- [40] H. Shao, W. Guo, and A. W. Sandvik, Phys. Rev. B 91, 094426 (2015).
- [41] We thank Kedar Damle for suggesting a definition based on the entire spinon string.

Supplementary Material for

Quantum Criticality with Two Length Scales

Hui Shao, Wenan Guo, and Anders W. Sandvik

In Sec. A we discuss the crossing-point analysis (phenomenological renormalization) underlying the finite-size scaling studies summarized in Fig. 2. We derive the scaling properties of the crossing points as a function of the system size and use the 2D Ising model as a bench-mark case to demonstrate the unbiased nature of the method when all sources of statistical errors and scaling corrections are considered.

In Sec. B we present the scaling arguments underlying the analysis of the domain-wall energy of the critical - model in Fig. 3(a). We begin with the simpler case of a generic system with discrete symmetry-breaking at a critical point with a single divergent length scale, deriving the scaling form of in the thermodynamic limit and for finite size. We then generalize to the case when the thickness of the domain wall diverges faster than the correlation length and discuss the different scenarios for finite-size scaling at criticality. We use the 2D Ising model as an example to illustrate Monte Carlo (MC) procedures we have developed for computing the free-energy differences needed for at thermal phase transitions. We also demonstrate conventional finite-size scaling in a classical system with a dangerously irrelevant perturbation; the 3D six-state clock model. For the - model at , we discuss calculations of the ground state energy with and without domain walls and supplement the results shown in Fig. 3(a) with results for a different domain-wall boundary condition, demonstrating consistent anomalous scaling with the same exponent ratio in both cases.

In Sec. C, to further motivate the two-length scaling form (2) and its unconventional limiting behavior [demonstrated numerically in Figs. 3(b,c)], we present derivations of the quantum-critical scaling forms of the spin stiffness and magnetic susceptibility by generalizing the approach of Fisher et al. [1] to the case of two divergent length scales.

In Sec. D we re-analyze previously published [30] results for the temperature dependence of the spin correlation length and the magnetic susceptibility of the critical - model. We demonstrate that their scaling anomalies can also be explained by exponents modified by the exponent ratio .

In Sec. E we provide some more details of the and (ground-state projector) QMC methods used in the studies of the - model.

## Appendix A Crossing-point analysis

The crossing-point analysis employed in Fig. 2 is an extension of Fisher’s “phenomenological renormalization”. We follow essentially the formalism developed and tested with numerically exact transfer-matrix results for the Ising model in Ref. [29], but apply it to QMC data. In Sec. A.1 we discuss formalities and derivations of the exponents governing the drifts of crossing points in the standard case, when there is a single divergent length scale. In Sec. A.2 we discuss why the single-length scaling form (1) can still be used to analyze crossing points and extract the exponent controlling the shorter length scale, even in the case when the criticality is described by the two-length ansatz (2) with the anomalous limit controlled by the longer length scale. In Sec. A.3 we discuss several practical issues and potential error sources (statistical as well as systematical) that should be properly taken into account when analyzing crossing points. We illustrate the procedure with data for the 2D Ising model, demonstrating the unbiased nature of the approach by reproducing the exactly known critical temperature and critical exponents to within small statistical errors.

### a.1 Scaling corrections and crossing points

Consider first the standard case of a single divergent length scale (correlation length) as a function of the distance to a critical point (a classical transition driven by thermal fluctuations at or a quantum phase transition at ). For some other singular quantity with the behavior in the thermodynamic limit (valid for , , or both, depending on the quantity) the finite size scaling is governed by the form

(S1) |

where and the variables are irrelevant fields which in principle can be tuned by introducing some other interactions in the Hamiltonian [36]. Keeping only the most important irrelevant field, using the notation for convenience, and suppressing the dependence on the unknown value of , we have Eq. (1) in the main text. The scaling function is non-singular and we can Taylor expand it in the neighborhood of the critical point;

(S2) |

For two system sizes and (), the two curves and take the same value (cross each other) at the point

(S3) |

Thus, in general the finite-size value of the critical point defined using such curve-crossing points shifts with the system size as . However, if the quantity is asymptotically size-independent at the critical point, , the first term in Eq. (S3) vanishes and the shift is faster;

(S4) |

where the constant of proportionality depends on the chosen aspect ratio and the generally unknown coefficients of the Taylor expansion (S2). The value of the quantity at the crossing point is obtained by inserting into Eq. (S2), which for both the general case and the special case can be written as

(S5) |

with some constants and . Thus, in principle a crossing point analysis can be used to obtain the leading critical exponents and as well as the subleading exponent . However, it should be noted that the higher-order terms in Eq. (S2) can play a significant role for system sizes attainable in practice, and often and extracted from fitting to power laws according to Eqs. (S4) and (S5) should be considered only as “effective” exponents which change with the range of system sizes considered (with the correct exponents obtained only for very large system sizes where the subleading corrections become negligible). To extract the critical point, a dimensionless quantity () should be chosen as the convergence then is the the most rapid, given by Eq. (S4). The value of the critical point obtained from fitting to this functional form is normally not very sensitive to the imperfection of the power-law correction with the effective value of the exponent, as long as the fit is statistically sound.

There are many other ways of analyzing crossing points. For instance, the exponent can be obtained more directly than the difficult extraction based on the correction terms in the shift analysis above. Consider a dimensionless quantity , such as the Binder ratio (or the corresponding cumulant). We then have, including also some terms of higher order in Eq. (S2),

(S6) |

and from the derivative with respect to or we have

(S7) |

We will now assume that is positive in the region of interest, and if not we redefine it with a minus sign. At we then have

(S8) |

with some constants and . Thus, for large , at the critical point depends linearly on and the slope is the exponent . A drawback of this method for extracting is that the critical point has to be determined first, and a careful analysis should also take into account the uncertainties in the estimated value of .

To circumvent the requirement of having to determine first, we observe that, instead of evaluating the derivative (S7) exactly at the critical point, we can use the crossing point of the quantity for two system sizes [or, as in Ref. [29], one can use with a constant , which only modifies some unimportant prefactors of the results derived below]. Inserting the crossing value (S4) of into (S7) we obtain

(S9) | |||||

Having access to two different slopes at the crossing point, we can take the difference of the logarithms of these and obtain

(S10) |

with some constant . We can therefore define an exponent estimate corresponding to the crossing point,

(S11) |

and this estimate approaches the correct exponent at the rate for large ;

(S12) |

with some constant and various higher-order terms again left out.

With all the crossing-point quantities discussed above, the infinite-size values , , and can be obtained by fitting data for several system-size pairs , using Eqs. (S4), (S5), and (S12). One can either use the leading form as written with only the asymptotically dominant correction (in the case of ) or (for and ) if the system sizes are large enough for the higher-order terms to be safely neglected, or one can include higher-order terms explicitly and fit to a larger range of system sizes. The former method has the advantage of the optimum fit being easier to find, while fits with multiple power laws are some times challenging or affected by large fluctuations in the parameters unless the statistical errors are very small.

### a.2 The case of two length scales

We now turn to systems with two divergent lengths, where the critical scaling is governed by Eq. (2). When the thermodynamic limit corresponds to the scaling function being a power of the first argument for large and , the effect of the second argument is the same as in the standard case of a dangerously irrelevant field scaling as . The crossing-point analysis then remains the same as in the previous section. In the anomalous case, which we have termed the super dangerous perturbation, the second scaling argument (the longer length scale) generically controls the behavior and demands the modified powers of in front of the scaling function. This case requires some additional discussion.

In general, the scaling in this case is much more complex. In the main paper we have discussed how the correct thermodynamic limit is obtained when the scaling function is controlled by . This limit corresponds directly to the intuitive physical picture of the shorter length saturating at when the longer length reaches , and, therefore, should not be replaced by at criticality but instead by . This change imposes an anomalous power law at criticality for any observable which can be written as some nonzero power of the correlation length close to the critical point. It should be noted that, there are special non-generic observables, such as the Binder ratio, which by construction neither have any -dependent prefactors of the finite-size scaling function nor any dependence on in the thermodynamic limit (e.g., the Binder ratio takes constant values in the phases and a different value at the critical point). In such non-generic cases there are also no modified power laws, since there are no powers to be modified by the ratio in the first place. All other generic observables are expected to develop anomalous power laws.

We next note that, in the above large limit of , both the arguments and become large. When we are interested in crossing points close to , we are far from this limit, however. We can anticipate crossing points as in the single-length case when the first argument is of order one (i.e., is of order ), whence the second argument is very small, . There is no a priori reason to expect that this limit is controlled by . The most natural assumption, which can be tested, is that is irrelevant in this regime. Then we are back at a situation where the standard crossing-point analysis can be performed and the exponent delivered by such an analysis should generically be , not . An exception is an observable which is manifestly dependent only on the longer length scale, in which case the shorter length scale will play the role of an irrelevant correction. The simplest quantity of this kind is a length scale which is proportional to the longer length itself. In the main text we have analyzed the size of the spinon bound state and found its crossing points to be controlled by an exponent exponent which is indeed significantly larger than , and also holds in the neighborhood of the critical point, as expected from the scaling function controlled by when in the thermodynamic limit.

We now have concluded that the limits of when and are controlled by different exponents in the generic case; by both and in the former case and only by in the latter case. This implies an interesting cross-over behavior between these limits. In principle, such a cross-over can be tested explicitly by numerical data, by graphing results for a wide range of system sizes and couplings (in the case of the - model, that should be done inside the VBS phase) against both and . One should observe data collapse onto common scaling functions in both cases, but only in the relevant regimes controlled by the different scaling arguments; small or large . It would clearly be desirable to carry out such an analysis for the - model, which we have not yet done due to the large computational resources required to do this properly for sufficiently large system sizes. We anticipate the analysis of the cross-over to be complicated also by the small exponent of the leading scaling corrections, as demonstrated in Fig. 2 in the main paper.

Even if no tests of the cross-overs are available currently, the two limits and have already been confirmed in this work; the former by the scaling of the Binder cumulant with the exponent (the shorter length scale) and the latter more indirectly by the presence of anomalous powers of . An anomalous exponent which is very well converged as a function of the system size and completely inconsistent with any other previous scenario (neither large scaling corrections nor a first-order transition) is best provided by the domain-wall energy , which is analyzed in Fig. 3(a) of the main paper and also further below in Sec. B.4.

### a.3 Tests on the 2D Ising model

In order to demonstrate the reliability of the method of obtaining the critical point and exponents from crossing points, and to discuss practical issues in implementing it, we here present results based on the Binder cumulant of the standard 2D Ising model;

(S13) |

where is the magnetization

(S14) |

MC simulations were carried out on lattices of size with periodic boundary conditions, using a mix of Wulff and Swendsen-Wang (SW) cluster updates, with each sweep of Wulff updates (where on average spins are flipped) followed by an SW update where the system is decomposed into clusters, each of which is flipped with probability . The SW clusters are also used to measure and with improved estimators (after each SW update). We carried out simulations of sizes , at temperatures in the neighborhood of the relevant crossing points of the Binder cumulant for system-size pairs , i.e., using aspect ratio in the expressions of Sec. A.1. Up to measurements were collected for the smaller sizes and for the largest sizes.

Figure S1 shows examples of data for three different system sizes, where cubic polynomials have been fitted to the data. The crossing points can be extracted numerically using bisection. In order to analyze and in the thermodynamic limit, it suffices to consider a small number of points very close to each crossing point to be analyzed. To obtain from the slopes according to Eq. (S11), where the derivative in Eq. (S7) is taken of the fitted polynomials, it is better to have a more extended range of points. However, for a very large range a high order of the polynomial has to be used in order to obtain a good fit, and it is then better in practice to adapt the window size so that a relatively low order polynomial can be used. In the tests reported here, cubic polynomials were used and all fits were statistically sound.

In order to compute the statistical errors (error bars) a bootstrap method can be used, i.e., by generating a large number of random samples of the binned MC data. Each bootstrap sample is computed using randomly chosen bins for each system size and temperature, where is also the total number of data bins available from simulations at . The standard deviations of the values (the horizontal and vertical crossing points and the slope estimator for ) computed for these bootstrap samples correspond to the error bars, which later will be used in the fits to extrapolate to infinite size. In evaluating the cumulant (S13), for the full data set or a bootstrap sample, the individual expectation values and should be computed first based on all the bins included in the sample, after which the ratio is evaluated. If one instead uses ratios computed for each bin separately, a statistically significant systematical error can be introduced in the ratio, due to nonlinear contributions to the statistical error which do not vanish as the number of bins is increased (for fixed bin size) but do decrease properly in the bootstrap method when the sample size is increased.

We next fit crossing points for a series of system pairs to the expected forms, Eqs. (S4), (S5) with , and (S12), and compare with exact and previous numerical results for the 2D Ising model. Onsager’s rigorous analytical solution gives and . The value of at is not known exactly, but Blöte obtained by extrapolating exact numerical finite-size transfer-matrix data to infinite size [37]. For the Binder cumulant the dominant subleading correction has the exponent [37]. These results should all be obtained within statistical errors from the crossing point analysis of the MC data if sufficiently large systems are used and the data are analyzed using appropriate statistical methods. For small sizes the expected higher-order corrections will cause deviations beyond the statistical errors from the leading-order forms, which can be detected in the goodness of the fits to the leading forms (S4), (S5), and (S12). Our strategy is to remove small system sizes until a statistically sound fit is obtained for a given quantity.

The crossing points for the different size pairs , , are not all statistically independent, because the same system size can appear in two different pairs. One should therefore define the goodness of the fit, per degree of freedom (the number of data points minus the number of parameters of the fit), with the full covariance matrix instead of just its diagonal elements (which are the conventional variances). Using to denote some quantity defined based on the crossing point (the crossing temperature , the value of of at the crossing point, or obtained from the slopes evaluated using the fitted polynomial), we thus use

(S15) |

where is either the mean value obtained from all available bins or an average obtained from the bootstrap procedure (the two estimates should differ only by an amount much smaller than the standard deviation based on the bootstrap analysis), is the value of the quantity evaluated using the fitted function (here a power-law correction to the infinite-size value), and is the total number of system-size pairs used. The covariance matrix is defined as

(S16) |

where the expectation value for each pair for which is again evaluated using bootstrap sampling (as explained above for the error bars, which correspond to the square-roots of the diagonal elements ). We use of the order bins and generate several thousand bootstrap samples to obtain accurate estimates of the covariance matrix.

To compute error bars on the extracted quantities, we repeat the fits to Eqs. (S4), (S5), and (S12) several hundred times using the bootstrap method and define the final means and statistical errors (one standard deviation) using these bootstrap samples. When defining as in Eq. (S15) for data fits based on bootstrap samples, the covariance matrix (S16) should be multiplied by a factor , due to the two statistically equal sources of fluctuations; the original MC fluctuations and those in the bootstrap samples. Then, for a statistically sound fit, is expected for the bootstrap-averaged goodness of the fit.

To quantitatively define a criterion for an acceptable fit, we consider the standard deviation of the distribution. For degrees of freedom, the standard deviation of is . We systematically eliminate the small sizes until falls roughly within two standard deviations of its expected mean;

(S17) |

Clearly this criterion is sensitive to the quality of the data—if the elements of the covariance matrix are very small, even fits including only relatively large system sizes can detect the presence of higher-order corrections and not pass our test, while with noisy data also small system sizes can be included (but the error bar on the final extrapolated value will be larger).

If a fit satisfies the goodness-of-fit criterion (S17) it can still not be completely guaranteed that no effects of the higher-order corrections are present in the final result, but in general one would expect any remaining systematical errors to be small relative to the statistical error. In principle one can estimate the magnitude of the systematical error using the parameters obtained from the fit and some knowledge or estimate of the nature of the higher-order corrections. We will not attempt to do that here, because in general such knowledge will be very limited. To minimize possibly remaining systematical errors one can continue to exclude more system sizes even after the soundness criterion (S17) is satisfied, at the price of increasing the statistical errors of the parameters extracted from the fits.

The above method implies a ’curse of good data’, as less data points are actually included in the final fit when longer simulations are carried out for a fixed set of system sizes. However, the discarded data still contain valuable information on the convergence properties and can in principle be used to analyze higher-order scaling corrections (which we do not pursue here).

Results for the horizontal (temperature) and vertical (cumulant) crossing values of the 2D Ising model are shown in Fig. S2. For the horizontal points in (a), our fits start to satisfy the criterion (S17) when including sizes (the average goodness of the fit is then with ) and we show that case in the figure. The fit gives and the exponent combination . Thus, the critical temperature comes out correct within the remarkably small error bar, while is about twenty of its error bars outside the true (asymptotic) value . As discussed above, it is typical in finite-size scaling that corrections-to-scaling exponents do not come out reliably until very large systems are used, and we therefore do not consider the mismatch as a failure here, rather as a confirmation of the known fact that the exponent should be considered as an “effective exponent” which slowly changes as larger system sizes are included.

For the crossing value of the cumulant we find a similar trend. In this case a good fit requires that only the points are used, giving and , again with (). The value deviates by about an error bar from Blöte’s result quoted above, while the correction exponent again is relatively far (considering the size of the error bar) from its asymptotic value . Interestingly, extracted as the difference of the two exponents comes out close to the correct value , within the statistical error.

The insets of Fig. S2 show the differences between the data points and the fitted curves. Here it can be seen that the points are not quite randomly distributed around , as they should be if the fitted functions are of the correct form. The overall shape with noisy but discernible minimums and maximums suggests the presence of a correction which is barely detectable for the range of system sizes at this level of statistics. One can then conclude that the deviations of by two standard deviations from in these fits are not purely statistical fluctuations (which is not clear from the values alone), but due to the neglected higher-order corrections. Nevertheless, the most important extrapolated values and were not adversely affected statistically, thus demonstrating the ability of the effective exponent and the prefactor of the correction term in Eqs. (S4) and (S5) to reproduce the overall trend of the data sufficiently well for extrapolating to infinite size.

To illustrate the effect of excluding even more system sizes, with the minimum size we obtain , two error bars away from the correct value (still a statistically acceptable match), and , also about two error bars from the previous (Blöte’s) value. From the fit we obtain in this case and from the fit . These exponents are now correct to within statistical errors, but the error bars are about 10 times larger than before, while the error bars on and only doubled. The average value of is very close to for both these fits and the deviations from the fitted function look completely random. Upon excluding even more points, the error bars increase rapidly but the extracted parameters remain statistically in good agreement with their correct values.

Next, we extract the exponent using the log-slope formula (S11). Fig. S3 shows the results along with a fit including all the system sizes (). Remarkably, the fit is statistically perfect, with , already at this small minimum size and the inverse exponent extrapolates to , in excellent agreement with the exact result . The slope data are much more noisy than the underlying values and the error bars grow very rapidly with for the largest sizes. The fit is therefore dominated by the smaller sizes. Naturally, the large error bars mask the effects of higher-order corrections, as discussed above. It is nevertheless remarkable that the extracted exponent does not show any effects of the neglected corrections at all, even though, again, the leading correction exponent, which comes out to , is not very close to the correct value and its error bar is large. Again, the flexibility of the leading finite-size term allows it to mimic the effects of the correction terms without significant effects in the extrapolation of the fit.

It should be noted that the 2D Ising model has logarithmic corrections in addition to the higher-order scaling corrections that we have neglected here [37], which is not a generic feature of critical points (except for systems at their upper critical dimension). The logarithms of multiply powers of higher than those of the leading corrections and we therefore do not expect them to affect the procedures used above.

These results demonstrate the unbiased nature of the crossing-point analysis when it is carried out properly. We have used the same scheme to analyze the results for the - model in Fig. 2 of the main text. In the left column, the behavior of is similar to that of of the Ising model in Fig. S2, with a relatively large correction exponent which makes the fits and extrapolations to stable and and visually convincing. In the right column, it is clear that the leading correction exponent for is small, , and that there are other significant corrections present in the top two panels. The fact that the critical point nevertheless agrees perfectly to within small error bars with that extracted from the spinon bound state is very reassuring. As in the Ising model, the fit to only requires a single scaling correction, though it can not be excluded that this correction is an effective one, mimicking the collective effects of several corrections with the same sign. In any case, the extrapolations are stable, e.g., excluding some of the small- points does not dramatically change the extrapolation, though of course the error bar grows.

We advocate the systematic curve-crossing method as outlined above to determine the critical temperature (or critical coupling of a quantum phase transition) and the critical exponents, instead of often used [also in DQC studies [15, 20, 22]] data-collapse techniques where many choices have to be made, concerning the range of data included, use of corrections, etc. Although trends when increasing the system size can also be studied with data collapse [as done in Ref. [20])], the solid grounding of the present scheme directly to the finite-size scaling form (S1) makes it the preferred method.

## Appendix B Domain-wall energy

As we discussed in the main text, the fundamental longer length scale in the DQC theory is the thickness of a domain wall in the VBS. In Fig. S4 we illustrate a generic domain wall in a 2D system in which a discrete symmetry is broken. In the case of a broken continuous symmetry, e.g., the magnetization vector in the XY spin model, there is no domain wall but the order parameter (its direction) gradually twists uniformly over the entire width of the system. This case will be discussed in Sec. C in the context of a twist of the Néel order parameter of the - model. For a discrete broken symmetry it is energetically favorable for the system to instead restrict the size of the region (the domain wall) over which the order parameter deviates significantly from the values imposed at the boundaries. Note, however, that the domain wall is not strictly fixed at some location, and, e.g., in an MC simulation the local order parameter will not detect the intrinsic width of a domain wall, because averaging is performed over all locations of the wall. Therefore, other means have to be employed to detect the intrinsic domain-wall thickness, e.g., using suitably defined correlation functions.

As we showed in the main text, the length scale is conveniently present in the - model in the finite-size scaling of the energy density of a VBS domain wall. Here, in Sec. B.1 we derive the scaling form of , in the thermodynamic limit and for finite system size, using a simple Ansatz generalizing the treatment by Fisher et al. [1] in a different context (considered further in Sec. C) to the case of discrete symmetry breaking with two divergent length scales. The formalism applies both to classical and quantum systems. We present our MC procedures to compute at classical (thermal) phase transitions, using the 2D Ising model as a concrete example in Sec. B.2. We also present results for the 3D classical six-state clock model at its critical temperature in Sec. B.3, before describing the details of the QMC calculations of for the - model at in Sec. B.4.

### b.1 Scaling forms

Let us first consider the case of a -dimensional system with single divergent length scale . Following Fisher et al. [1], we consider the singular part of the free-energy density, which we can write for a classical system at finite temperature or a quantum system at (in which case the free energy is just the ground state energy) as

(S18) |

where formally the dynamic exponent for a classical system. Introducing a domain wall, the free-energy difference with respect to the system without domain wall should scale in a similar way but with a different size-dependent function [1];

(S19) |

This density should be understood as a quantity averaged over the inhomogeneous system (or, equivalently, in a finite system the domain wall location is not fixed and all properties are averages over all locations of the domain wall), and the total free-energy difference is

(S20) |

where is the volume of the system.

We can also write down a different expression for the free-energy difference, by explicitly considering the cost of twisting the order parameter. If the domain wall has width and the total twist of the order parameter across the wall is , then the cost per lattice link inside the wall is , which also defines the stiffness constant . Outside the wall region the local energy cost vanishes, and, since the total volume occupied by the domain wall is we have

(S21) |

Consistency in the dependence between this expression and Eq. (S20) requires that the scaling function has the form , and therefore

(S22) |

The domain wall energy per generalized cross-section area of the wall (its length for , area for , etc.) is then

(S23) |

which no longer has any dependence and, thus, represents the behavior in the thermodynamic limit. We can also read off the scaling of the stiffness constant,

(S24) |

Since we have written all expressions in terms of the correlation length, we can now switch to finite-size scaling at a critical point by simply making the substitution . For the domain wall energy (S23) of interest here we obtain

(S25) |

Now consider a system with two length scales, with a conventional correlation length and a domain wall thickness , with . A simple generalization of Eq. (S20) suggests that

(S26) |

Note that only the shorter length scale should appear in front of the size-dependent scaling function because the free energy in the thermodynamic limit should only depend on the two lengths in an additive way, , in order for the specific-heat exponent () relation to hold, i.e., for hyper-scaling to apply (which we thus assume). Since diverges slower than , is asymptotically dominated by the term, and (S26) should then describe the leading singular behavior.

We can also easily generalize Eq. (S21) to a domain wall of thickness ;

(S27) |

Now consistency between Eqs. (S26) and (S27) for both the dependence and the dependence requires that , and we arrive at

(S28) |

for the scaling of in the thermodynamic limit. Note the consistency of this form and the single-length form (S23) when . In the particular case of a DQC point (, ), Eq. (S28) reduces to , which was derived in a different way by Senthil et al. [4].

To convert Eq. (S28) to finite-size scaling, in the standard treatment of two length scales arising from a dangerously irrelevant perturbation [31], the longer scale is not present in the leading finite-size scaling behavior. This can be understood physically as follows: Upon approaching the critical point from the ordered phase, when reaches we simply replace by . However, continues to grow and controls the scaling behavior until it reaches . At the critical point also is replaced by , and the critical finite-size scaling of obtained from (S28) is, thus, identical to the single-length form (S25). Since neither nor appear here, there is no information on these exponents in the finite-size scaling of in the standard scenario.

As we argued in the main text, there is also another possibility, namely, the growth of in Eq. (S28) is halted when reaches . Then , leading to the finite-size scaling

(S29) |

In the case of DQC, this reduces to . It is very interesting that the ratio appears here in a simple way and can be extracted using critical finite-size scaling. The result in Fig. 3(a) leaves little doubt that , which represents unambiguous evidence for anomalous scaling in the - model. Below, in Sec. B.4, we will present details of these calculations, along with additional results showing that for the - model.

### b.2 2D Ising model

It is instructive to first test the domain-wall scaling using a simple system such as the 2D Ising model. A domain wall in the ferromagnet can be enforced in different ways using suitable boundary conditions. Here we use systems with periodic boundaries in the -direction and compare two different boundaries, as illustrated in Fig. S5. The boundaries are open, with the edge columns coupled with the same strength as the bulk coupling to fixed spins and , equivalent to boundary fields of strength . Here the domain-wall imposing column of spins to the right extends only partially through the system, to illustrate the mechanism we use for computing the required free-energy difference.

It is not easy to compute the free energy in MC simulations, but it is relatively easy to compute a free-energy difference, if the two systems of interest, let us call them “1” and “2”, can be simulated collectively as a partition function . If there are updates switching the simulation between system states 1 and 2 with detailed balance satisfied, then the free-energy difference , where are the probabilities of the simulation “visiting” the respective states. Such multi-canonical simulations [38] can be extended to an arbitrary number of systems , and any can then be accessed, provided that the simulation can easily transition between the different states .

In the studies of domain walls considered here, the different systems correspond to boundary conditions fluctuating between the normal periodic boundaries and the domain-wall boundaries. To enhance the ability of the system to fluctuate between these boundary conditions of interest, the whole boundary is not changed at once, but in small steps where the right boundary has a change from to at some vertical location , as illustrated in Fig. S5. Thus, corresponds to the normal periodic boundaries (no domain wall) and corresponds to the boundary enforcing a full vertical domain wall. For the domain wall does not extend vertically through the whole system and instead has a horizontal part connecting to the location where the boundary changes. MC updates are used to move this location, , using heat-bath acceptance probabilities.

We find that the probability of the boundary conditions generated is the highest, as expected, for . There is also a local maximum at , and a minimum around . To further increase the efficiency of the boundary moves, a weight factor is multiplied with the Boltzmann probability for the spins and gradually adjusted until the histogram of the relative number of times the boundary is at becomes almost flat. Then, the actual probability without the re-weighting factor is , and the free-energy difference between the systems with and without domain wall is (leaving out the unimportant temperature factor),

(S30) |

MC results for are shown in Fig. S6. The inset shows the running exponent extracted on the basis of size pairs by postulating and , whence . The results are fully compatible with when , as predicted by Eq. (S25) when , with a correction . We have also carried out simulations of the 3D Ising model and confirmed that .

### b.3 3D clock model

The existence of two length scales in the DQC theory relies heavily [3, 4] on an analogy with the classical 3D clock model, where the standard XY model is deformed by an external potential for all the angles . This term is known to act as a dangerously-irrelevant perturbation, leading to a domain-wall thickness . It is therefore natural to also test the scaling of the domain-wall energy in this case. Here we use the standard XY interaction between nearest neighbors on the 3D simple cubic lattice

(S31) |

where the angles are constrained to the clock angles, , . The hard constraint is equivalent to the limit with the cosine perturbation.

The exponent should be independent of (including the fully-constrained limit considered here) but depends on , diverging as . There has been some controversy regarding methods to compute the exponent in MC simulations, as summarized in the recent Ref. [31], but for small several calculations are nevertheless in good agreement with each other and we can use them as reference points.

In order for the exponent ratio to be significantly different from one we here use , in which case and, since the 3D XY exponent , the ratio . Results for the domain-wall energy scaling at the critical point are shown in Fig. S7. The results are completely consistent with the form (S25) with , corresponding to the expected standard scenario where finite-size scaling is obtained from the thermodynamic-limit form by replacing both divergent length scales by . The results are completely inconsistent with the alternative scenario (S29), where the decay exponent should approach . This result reinforces the unusual form of the scaling of in the - model, Fig. 3(a) of the main text, which we will discuss in more detail in the next section.

We also comment on the applicability of the generic two-length scaling form (2) in the main paper to in the clock model. Using the finite-size scaling we found above, we should have

(S32) |

To obtain the correct thermodynamic limit, when , we must have , which is also natural because, given the form in the thermodynamic limit, should be separable, , where the two factors just correspond to the expected scaling forms for the length scales and themselves. In contrast, in the - model we have argued for an anomalous form which corresponds to a generally non-separable scaling function with the thermodynamic limit controlled only by the second argument.

### b.4 J-Q model

In the - model we are interested in ground state energies of systems with and without domain walls and these can be computed in standard QMC simulations. The multi-canonical approach employed in the previous section, developed to circumvent the difficulties of MC calculations of individual free energies at , are therefore neither useful nor needed. We use the projector QMC approach with applied to a valence-bond trial state of the amplitude-product type [39, 40], choosing the “projection time” sufficiently large, up to , to converge the ground-state energy. Domain walls are introduced by boundary conditions in two different ways, schematically illustrated in Fig. S8.

The VBS order parameter is a vector , where the operators corresponding to the two components can be defined as