Deconfined Quantum Criticality, Scaling Violations, and Classical Loop Models

Numerical studies of the N\'eel to valence-bond solid phase transition in 2D quantum antiferromagnets give strong evidence for the remarkable scenario of deconfined criticality, but display strong violations of finite-size scaling that are not yet understood. We show how to realise the universal physics of the Neel-VBS transition in a 3D classical loop model (this includes the interference effect that suppresses N\'eel hedgehogs). We use this model to simulate unprecedentedly large systems (of size $L\leq 512$). Our results are compatible with a direct continuous transition at which both order parameters are critical, and we do not see conventional signs of first-order behaviour. However, we find that the scaling violations are stronger than previously realised and are incompatible with conventional finite-size scaling over the size range studied, even if allowance is made for a weakly/marginally irrelevant scaling variable. In particular, different determinations of the anomalous dimensions $\eta_\text{VBS}$ and $\eta_\text{N\'eel}$ yield very different results. The assumption of conventional finite-size scaling gives estimates which drift to negative values at large $L$, in violation of unitarity bounds. In contrast, the behaviour of correlators on scales much smaller than $L$ is consistent with large positive anomalous dimensions. Barring an unexpected reversal in behaviour at still larger sizes, this implies that the transition, if continuous, must show unconventional finite-size scaling, e.g. from a dangerously irrelevant scaling variable. Another possibility is an anomalously weak first-order transition. By analysing the renormalisation group flows for the non-compact $CP^{n-1}$ model ($n$-component Abelian Higgs model) between two and four dimensions, we give the simplest scenario by which an anomalously weak first-order transition can arise without fine-tuning of the Hamiltonian.


I. INTRODUCTION
The paradigmatic "deconfined" quantum phase transition is that separating the Néel antiferromagnet from the columnar valence bond solid (VBS) for a square lattice of spin-1=2s. The theoretical arguments of Refs. [1][2][3] indicate that the Néel-VBS phase transition is described by the noncompact CP 1 (NCCP 1 ) model [4], a field theory with bosonic spinons z ¼ ðz 1 ; z 2 Þ coupled to a noncompact U(1) gauge field: This theory is defined in three-dimensional Euclidean spacetime; the Néel order parameter is proportional to z †σ z, whereσ are the Pauli matrices.
Numerical results for the J-Q model (the Heisenberg antiferromagnet supplemented with a four-spin interaction [5]) support the validity of this continuum description [6][7][8][9][10][11], as does work on the SUðnÞ generalization of the problem at large n [12][13][14][15]. Unfortunately, though, the existence of a continuous phase transition in both the NCCP 1 model and the SU(2) lattice magnets remains a vexed question. While simulations of the J-Q model are compatible with a direct continuous transition, they show strong violations of finite-size scaling [9,10,16,17]. These persist up to the largest system sizes studied so far and hamper the extraction of meaningful critical exponents [10]. Additionally, direct numerical studies of the lattice NCCP 1 field theory have disagreed as to whether the transition is continuous [4,18] or whether scaling violations similar to those seen in the lattice magnets should be interpreted as the initial stages of runaway flow to a firstorder transition [19].
Are the scaling violations seen at the Néel-VBS transition indeed signs of a first-order transition, with an anomalously large correlation length [16,17,19], are they due to the critical theory possessing a weakly irrelevant scaling variable [9,20,21], or do they indicate something more exotic? This issue remains controversial. Its relevance extends beyond quantum magnets since the critical behavior of the NCCP 1 model is important for various other fundamental problems in statistical mechanics. For example, this field theory is believed to describe the threedimensional classical O(3) model when hedgehog defects are disallowed [4], as well as the columnar ordering transition in the classical dimer model on the cubic lattice [22][23][24][25][26][27]. [In the latter example, SU(2) symmetry is absent microscopically but argued to emerge at the critical point.] There is also numerical evidence that similar scaling violations afflict the SU(3) and SU(4) generalizations of the deconfined transition [10,28].
In this paper, we introduce a new model which is ideally suited for studying the universal features of the Néel-VBS transition and perform simulations on very large systems (of linear size up to 512 lattice spacings, and 640 for a few selected observables). We verify that the model shows the basic features expected from the NCCP 1 field theory [Eq. (1)]: an apparently continuous direct transition, with emergent U(1) symmetry for rotations of the VBS order parameter at the critical point. However, we show that scaling violations are even stronger than previously appreciated. Conventional finitesize scaling assumptions are not obeyed: The data cannot be made to show scaling collapse, and quantities that would normally be expected to be universal instead drift with system size. The larger sizes considered here show that these drifts are stronger than the logarithmic form conjectured previously [8,9].
In common with Ref. [10], we see a drift in finite-size estimates of critical exponents. We show that this is more drastic than previously apparent. Estimates of the anomalous dimensions of both the Néel and VBS order parameters, as extracted from the correlation functions GðrÞ at distances r comparable with the system size (e.g., r ¼ L=2), yield negative values at large sizes. Negative anomalous dimensions are ruled out for a conformally invariant critical point by the unitarity bounds [29,30]. This, together with the form of the drifts mentioned above, rules out a conventional continuous transition with conventional finite-size scaling.
On the other hand, the decay of GðrÞ with r for r ≪ L appears consistent with the large positive anomalous dimensions suggested for a deconfined critical point. It is conceivable that the transition could be continuous but that conventional finite-size scaling could fail as a result of a dangerously irrelevant variable [28]. As we discuss, in this scenario correlators GðrÞ with 1 ≪ r ≪ L would be expected to show the true positive anomalous dimensions, while correlators with r of order L would behave anomalously (as in, e.g., ϕ 4 theory above 4D [31]). The hypothetical dangerously irrelevant variable discussed here should not be confused with the much discussed Z 4 anisotropy for the VBS order parameter (Sec. IV B), which is dangerously irrelevant in a different sense.
Therefore, unless there is a reversal of the drift in exponents at still larger sizes, which seems unlikely, there are two possibilities: Either the transition is continuous with unconventional finite-size scaling behavior (for example, as a result of a dangerously irrelevant scaling variable), or it is first order. We will discuss both possibilities but cannot rule out either. We do not see the conventional signs of a first-order transition, such as double-peaked probability distributions for the energy and other quantities.
On the other hand, an alternative hypothesis put forward previously-that the scaling violations are due simply to a weakly or marginally irrelevant scaling variable [8,9,20]is not supported by our data. We also rule out any explanation in terms of unconventional dynamic scaling, i.e., deviations from dynamical exponent z ¼ 1: Our model has z ¼ 1 by construction since it is isotropic in three dimensions. This isotropy is also a convenient feature from the point of view of simulations.
Turning to theory, we analyze the topology of the renormalization group (RG) flows in the NCCP n−1 model between two and four dimensions, in order to assess the possibility of an anomalously weak first-order transition. This analysis unifies what is known about this field theory in 4 − ϵ dimensions, in 2 þ ϵ dimensions, and at large n, and extends previous partial results [32]. Treating n as continuously varying, we argue that in 3D, there is a universal value n Ã below which the deconfined critical point disappears by merging with a tricritical point. The argument does not fix the value of n Ã , which could be greater or smaller than 2, so it does not tell us whether the NCCP 1 model has a continuous transition. However, it does have the following consequences. If n Ã happens to be greater than 2, there is a possible mechanism by which a very weak first-order transition can appear for a range of n values (i.e., a large correlation length can be obtained without the need for finetuning of the Hamiltonian). The argument also shows that n Ã is greater than 1. This means that the inverted XY transition in the model with n ¼ 1 is not analytically connected to the critical point in the large-n regime of the NCCP n−1 model, contrary to assumptions made in previous work.
This picture for the RG flows also clarifies that the usual 2 þ ϵ expansion of the O(3) sigma model does not describe the conventional O(3) transition in 3D but rather the deconfined critical point (if it exists at n ¼ 2). This is natural: Hedgehogs are crucial in determining the critical behavior of the O(3) model in 3D [4], and the 2 þ ϵ expansion presumably fails to account for them. The conclusion is also in line with the RG result that the 2 þ ϵ approach to the 3D OðMÞ model should fail when M is less than a universal value M c , conjectured to be above 3, as a result of neglecting the topology of the sphere [33]. Interestingly, our best estimates for the correlation-length exponent at the deconfined transition (Sec. IV E) are close to ν ¼ 1=2, smaller than most previous estimates but in good agreement with the 2 þ ϵ predictions (Sec. VI D).
Returning to lattice models, our numerical strategy is, instead of focusing on a simple two-dimensional quantum Hamiltonian, to construct a simple 3D classical model that is well adapted to large-scale Monte Carlo simulations. While the correspondence with classical lattice models in one dimension higher is a standard tool for studying quantum phase transitions, one might, at first glance, think that this tool is not available for deconfined criticality. This is because deconfinement relies crucially on the fact that the Euclidean action for the spins in 2 þ 1 dimensionsunlike the energy functional for a classical spin model in three dimensions-contains imaginary terms (Berry phases). The effect of these terms is to endow hedgehogs in the Néel order parameter with position-dependent complex fugacities [34][35][36]. After coarse-graining, this leads to a phase cancellation effect that suppresses hedgehogs [1][2][3].
Contrary to the naive expectation above, we show that the remarkable physics of deconfinement, including the suppression of hedgehogs by phase cancellation, is present in our 3D classical model. This model is formulated in terms of configurations of loops on a lattice and is a variant of the models of Refs. [37,38]. The loop configurations have positive Boltzmann weights, so they define a conventional classical statistical mechanics problem. However, the partition function can also be mapped to a lattice field theory for CP n−1 spins, and in this representation, the Boltzmann weights are not necessarily real. We show by a direct calculation that they include the complex hedgehog fugacities necessary for deconfined criticality.
The loop model introduced here has qualitative features in common with loop ensembles arising in worldline quantum Monte Carlo techniques for sign-problem free Hamiltonians such as the J-Q model [39] (see also Ref. [40]). However, direct simulation of a quantum Hamiltonian leads to an ensemble of worldlines in continuous imaginary time, whereas the loop model is an isotropic three-dimensional lattice model. This is a desirable feature for numerical simulations as it fixes an otherwise unknown velocity and eliminates a potentially significant source of corrections to scaling [41]. The geometric form of the model also motivates new observables-for example, we find it useful to consider some percolationlike observables such as the number of system-spanning strands and the fractal dimension of the loops.

II. LOOP MODEL
An astonishing variety of critical phenomena can be studied using classical loop gases. The present lattice model involves two species (colors) of loops, or n colors in the SUðnÞ generalization. It has a phase in which infinite loops proliferate and one in which all loops are short. The short-loop phase spontaneously breaks lattice symmetry because the system must choose between four symmetryrelated ways to pack the short loops.
The transfer matrix for loop models of this kind gives a correspondence with a 2D quantum magnet on the square lattice [38]. The color of a strand is related to the state of the spin (at a given point in Euclidean space-time) in the quantum problem. The infinite-loop phase corresponds to the Néel phase: The presence of infinite loops is equivalent to the presence of long-range spin correlations. The four degenerate short-loop phases map to the four equivalent columnar VBS patterns on the square lattice. The schematic correspondence between the loop model and the continuum field theory [Eq. (1)] is that the two species of loops are worldlines of the two species of bosonic spinons (z 1 , z 2 ). (See Sec. V for more details on the continuum limit.) The loop model is a modification of those studied in Refs. [37,38], with an additional interaction chosen to drive the model through a transition without explicitly breaking the symmetry between the four short-loop (VBS) states. The model lives on a four-coordinated lattice with cubic symmetry proposed by Cardy [42]. This "3D L lattice" is shown in Fig. 1 (left). Formally, it can be defined by starting with two interpenetrating cubic lattices, C 1 and C 2 , with lattice spacing 2, displaced from each other by (1,1,1): The faces of C 1 intersect the faces of C 2 along lines: These define the links of the L lattice. The L lattice is bipartite; its two sublattices are marked in yellow and black in Fig. 1.
We orient the links of the L lattice such that each node has two incoming and two outgoing links, with the two incoming links parallel and the two outgoing links parallel, as in Figs. 1 and 2. This assignment is unique up to a reversal of all orientations.
Breaking up each L lattice node by pairing the links in one of the two ways shown in Fig. 2 gives a completely packed loop configuration. In the simplest case, the partition function is just the equal-weight sum over all such configurations, with one of n colors assigned to each loop (the case of main interest here is n ¼ 2): We can, of course, perform the sum over colors explicitly to give Z ¼ P loop configs n no: loops . Note that the loops are automatically consistently oriented since the pairing of links at a node is always between an incoming and an outgoing link.
With the above Boltzmann weight, the system is (for n ≤ 4) in a phase where infinite loops proliferate [37,38]. We wish to add an interaction that drives the system into a phase with only short loops. First, consider the extreme limit of such a phase, which is a configuration in which every loop has the minimal possible length (which is six links). There are only four such configurations, and they are related by lattice symmetry. One is shown in Fig. 1, right.
A general loop configuration is determined by the binary choice of link pairing at each of the nodes. We denote this binary degree of freedom at node r by an Ising-like variable σ r ¼ AE1. To fix the sign convention for σ, let us pick one of the four minimal-length configurations as a reference and declare that all σ r are equal to þ1 in it. The four minimallength configurations are then those in which all the σs on the same sublattice (A or B) have the same sign. We can define an order parameter of the schematic formφ ¼ ðσ A ; σ B Þ which distinguishes the four short-loop phases. This is the analogue of the VBS order parameter in the quantum problem.
We introduce an interaction between nearest-neighbor σs on the same sublattice (i.e., between nodes of like color in Fig. 1): (The sum over uncolored loop configurations is equivalent to a sum over the σs.) With this choice, there is a direct transition at J c between a phase that has extended loops and hφi ¼ 0, and one that has only short loops and hφi ≠ 0.
As we would expect from the quantum correspondence [38], the continuum description of the above model is the NCCP n−1 model. In Sec. V, we show this directly by mapping the loop model to a lattice CP n−1 model and coarse graining, paying special attention to the fate of hedgehogs.

III. OVERVIEW OF RESULTS
We first summarize the salient results of our simulations. At the most basic level, they confirm that the loop model shows the central features of the deconfined Néel-VBS transition and that it probes the same universal physics as the J-Q [5] and related quantum models.
We find a direct and apparently continuous transition between Néel and VBS phases. Figure 3 shows the order parameters for these phases, for various system sizes L, very close to the critical point (see details in Sec. IVA). The data suggest a single transition. This is confirmed by examining finite-size pseudocritical couplings J c ðLÞ determined from various observables (inset to Fig. 3); all extrapolate to the same value as L → ∞ within error bars, so we are confident that there is a single transition at At small sizes, the estimates of critical exponents are compatible with those found in the J-Q model at similar sizes and in direct simulations of the NCCP 1 model [6,7,9,10,18]. As expected [1,5,7], we see an emergent U(1) symmetry for rotations of the VBS order parameterφ close to this critical point. field in the continuum action Eq. (1) [1,2].] Within the VBS phase, the emergent U(1) symmetry survives up to a length scale ξ VBS that is parametrically larger than the correlation length ξ. As for the J-Q model [7], the U(1) symmetry is apparent in the probability distribution forφ; see Fig. 4. Despite the above features, finite-size scaling properties at the transition are anomalous in various ways. For example, an appropriately defined stiffness for the Néel vector-which would be a universal constant at a conventional critical point-increases slowly with system size, and the critical exponent estimates also drift as the size is increased. Similar features were seen in previous numerical work on the J-Q model [5,9,10,16,17], but the larger sizes considered here show that the scaling violations are stronger than previously apparent. For a detailed picture of the transition, we analyze a variety of observables. Violations of finite-size scaling are visible in almost all quantities and do not decrease as L is increased. For this reason, we are unable to fit the size dependence of the data near the critical point assuming either scaling corrections coming from an irrelevant scaling variable (even if it is very weakly irrelevant) or logarithmic corrections similar to those considered in Refs. [8,9]. (See Secs. 3, IV C, IV E.) Figure 5 shows the "spanning number" N s versus the coupling J for various system sizes. N s is the average number of strands which span the system in a given direction, and it is a measure of the stiffness of the Néel order parameter. Instead of tending to a universal value as dictated by standard finite-size scaling, the crossing points N Ã s drift upwards as a power of L: We calculate the correlation-length exponent ν and the anomalous dimensions η Néel and η VBS using several observables. In the text, results will be presented with statistical errors in the last significant digit shown in brackets in the usual way; for reasons that will be apparent, we are not generally able to estimate systematic errors.
Estimates of ν obtained from finite-size scaling analyses of different quantities are in reasonable agreement but drift significantly, from ν ≳ 0.6 at small sizes to values around ν ∼ 0.46 for the largest sizes. In contrast, values of ν obtained from the variation of the correlation length with distance from the critical point lie in the range 0.45-0.5 with less dependence on size (Sec. IV E).
Strikingly, the behavior of the correlation functions at the critical point suggests different values of the anomalous dimensions η Néel and η VBS depending on the range of r used to extract them. Values obtained from correlation functions at separation L=2 both drift from values above 0.2 to values below zero at large sizes. Negative ηs violate the unitarity bound η ≥ 0 [29,30]. In contrast, there is evidence that behavior for r ≪ L is consistent with positive values for the anomalous dimensions. We note that the use of correlators at separation L=2 to determine η assumes finite-size scaling, which is a stronger assumption than that of the continuity of the transition, as we discuss in Sec. IVA 1.
In view of the above, the transition can only be continuous if some subtlety invalidates the usual finitesize-scaling expectations. (Of course, in principle, there could be a drastic change in behavior at still larger sizes L ≫ 512, but the data give no reason to expect this.) Therefore, it is natural to ask whether the transition is first order, with an anomalously large correlation length. But while the probability distributions of various quantities show violations of finite-size scaling, we do not see the standard signs of an incipient first-order transitiondouble-peaked probability distributions, etc. (Secs. IVA 4, IV C 1, IV D). Figure 6 shows the heat capacity C, which quantifies the fluctuations of the energy. This has a diverging peak at large sizes, as expected for a critical point with a positive heat capacity exponent (α ¼ 2 − 3ν > 0). The peak only emerges from the background at relatively large L. Surprisingly though, the peak fits well to a power law after subtracting a constant to account for the background, C max − C 0 ∼ AL α=ν (Fig. 6, upper inset). This gives ν ≃ 0.44, corresponding to a divergence ∼L 1.52 . This divergence is much slower than the L 3 expected asymptotically at a first-order transition. For a more intuitive picture, the lower panel of Fig. 6 shows the standard deviation of the energy, divided by the volume, at J ¼ 0.08850. For a first-order transition, this should saturate to a constant (proportional to the square of the difference in energy density between the two phases), while here there is no sign of saturation. (More details can be found in Sec. IV D.) To shed light on these perplexing observations, we analyze the topology of the RG flows in the NCCP n−1 theory in Sec. VI. The topology we find allows for a scenario with an anomalously weak first-order transition for a range of n, as a result of a coupling which "walks" (runs slowly) in the proximity of a fixed point located at a spatial dimension slightly below three. This is one possible reconciliation of the above numerical observations.
A more radical possibility is that the transition is continuous but disobeys finite-size scaling because of a dangerously irrelevant variable. This was hypothesized for the SU(3) and SU(4) cases in Ref. [28]. As pointed out below, in this scenario we expect scaling violations in correlation functions when the separation r of the points is comparable with L, but not when r is fixed and L → ∞. In Secs. IVA 1 and VII we consider this possibility in the light of the data.
At present, we cannot rule out either scenario; we sum up the situation in Sec. VII.

A. Néel and VBS order parameters and correlators
The deconfined transition separates phases that break different symmetries. In the VBS (short-loop) phase, lattice symmetry is broken: This is quantified by the order parameterφ introduced in Sec. II, whose spatially uniform part isφ This is normalized so jφj 2 ¼ 1 for perfect VBS order (there are N sites =2 sites on each sublattice). In the Néel (infiniteloop) phase, SU(2) spin-rotation symmetry is broken. In the loop representation, the magnitude of the Néel order parameter is the probability N that a given link lies on an infinite loop [38]. For a finite system, one may define N to be the probability that a link lies on a strand that spans the system in the z direction. If the transition is second order, we expect finite-size scaling forms [44] forφ and N, where δJ ¼ J − J c . However, attempting a scaling collapse using these forms gives negative ηs and very poor collapse. Raw data for the order parameters were shown above in Fig. 3.

Correlation functions
Next, we examine the critical two-point correlation functions forφ and the Néel vector. In the loop representation, the Néel correlator is simply the probability that two links lie on the same loop [38]. We denote these correlators G VBS ðr; LÞ and G Néel ðr; LÞ, where r is the separation of the points (taken parallel to a coordinate axis) and L is the system size. Raw data are shown in Fig. 7.
Conventionally, at J c one would expect with different ηs and different scaling functions c for each of the two observables. This would imply a collapse when plotting L 1þη Gðr; LÞ against r=L. Here, this collapse fails because the effective values of η at small and at large distances differ, as we now quantify. The full correlation function is relatively complicated because it depends on two length scales, r and L. Therefore, a standard approach is to examine G Néel ðL=2; LÞ and The effective values η Néel ðLÞ and η VBS ðLÞ determined from the slope are shown in Fig. 8 (lower inset). Note that for large L, the estimates for both exponents reach negative values. As mentioned above, negative values of η VBS or η Néel are ruled out for a continuous phase transition governed by a conformally invariant fixed point (though see below).
Another way to quantify the violation of finite-size scaling is via the ratios which should be universal according to Eq. (10) but instead drift significantly with L; see Fig. 8, upper inset. At certain critical points-for example, in ϕ 4 theory above 4D-a dangerously irrelevant variable invalidates standard finite-size scaling for the correlators. In this scenario, it may happen that the correlator is conventional in the limit L → ∞ [i.e., Gðr; ∞Þ ∼ r −1−η with η ≥ 0] but anomalous when r is comparable with L, or even with a smaller power of L (see below). Although a priori there is no theoretical reason to expect this phenomenon here, it suggests examining correlators in the regime r ≪ L. From   Fig. 7, it is conceivable that a well-defined power law will emerge in the limit L → ∞, although if so, the convergence in L is rather slow.
For further insight, we examine the derivatives of the correlators, G 0 Néel ðrÞ ¼ dG Néel ðrÞ=dr and G 0 VBS ðrÞ ¼ dG VBS ðrÞ=dr, in Fig. 9. Remarkably, these quantities show quite clean power-law behavior up to at least r ∼ 100, with and quite good scaling collapse (insets to Fig. 9). Reasonable scaling collapse is also obtained for "subtracted" correlators, defined as Gðr; LÞ − GðL=2; LÞ.
Straight lines corresponding to the above exponent values are shown in Fig. 7 for comparison with the raw data.
These results indicate that the strongest effect of the scaling violations is on the zero modes of the fields, which have anomalously large fluctuations. This is suggestive because in ϕ 4 theory above 4D, where the quartic term is dangerously irrelevant, the violation of finite-size scaling is caused by anomalously large fluctuations of the field's zero mode [31]. (Unlike the modes at nonzero momentum, which appear in the kinetic term, the critical fluctuations of the zero mode are controlled only by the irrelevant quartic term.) In that simple theory, we can easily separate out the contribution of the zero mode to the critical correlation function of the scalar field. With periodic boundary conditions, this gives (d is the dimension, and u is proportional to the bare coefficient of the quartic term) The second term, which violates conventional finite-size scaling, dominates as soon as r ≳ L d=ð2d−4Þ . Since the contribution of this mode to the two-point function depends on L but not on r, scaling can be repaired in this case by differentiation or subtraction.
The fact that this works perfectly in ϕ 4 theory is expected to be a special feature of the fixed point being free (and of the choice of correlator). Nevertheless, the good scaling of G 0 Néel and G 0 VBS at the deconfined transition is striking given the strong violation of scaling for the correlators themselves (Fig. 8), and may possibly indicate that a dangerously irrelevant variable is playing a role in the scaling violations. The example of vector ðφ 2 Þ 2 theory also makes it clear that a dangerously irrelevant variable of this typeassociated with controlling anomalously large fluctuations of zero modes-would lead to a spin stiffness that diverges at the critical point instead of taking a universal value [45].
In an appropriate loop gas picture for ðφ 2 Þ 2 theory, this diverging spin stiffness is associated with the appearance of "anomalously long" loops at the critical point.

Fractal structure of loops
The geometrical interpretation of the anomalous dimension η Néel is in terms of the fractal dimension of the loops, which according to conventional scaling relations, is given by d f ¼ ð5 − η Néel Þ=2, and determines the power-law relation between the root-mean-square end-to-end distance R of a strand and its length (see Ref. [38] for details). This again gives a large positive η Néel , in contrast to the drift towards negative values seen in the estimate from G Néel ðL=2Þ. The simplest fit, taking strands with R ≲ 100 to minimize effects of finite R=L, gives η Néel ¼ 0.42ð6Þ (data not shown). We note that this is considerably larger than Eq. (12). However, attempting to include finite R corrections in the fit gives smaller values in the range 0.25 ≲ η Néel ≲ 0.42 [47].

Susceptibilities
To compare with the estimates above, we calculate η Néel and η VBS from the Néel [38] and VBS [48] susceptibilities. These are shown in Fig. 10. According to finite-size scaling, the peaks should diverge as L 2−η Néel;VBS . The insets show log-log plots of the peak heights against L. The slopes indicate a downwards drift from η VBS ¼ 0.164ð13Þ to η VBS ¼ −0.35ð10Þ and from η Néel ¼ 0.355ð9Þ to η Néel ¼ 0.126ð3Þ. Figure 11 shows the Binder cumulant for the VBS order parameter, defined as

Binder cumulant
At a first-order transition, there should be a dip in U VBS which diverges with the system size [49]. In our case, there is no sign of this dip. In the inset to Fig. 11, we plot the maximum value of the slope dU VBS =dJ for each L. For a second-order transition, this value diverges as L 1=ν at the critical point. From the inset, we see that there is different behavior for small and large system sizes, giving ν ¼ 0.62ð1Þ for sizes L ≤ 64 and ν ¼ 0.476ð18Þ for L ≥ 256. (See Sec. IV E for other estimates of ν.)

B. Emergent symmetries
The deconfined criticality scenario assumes that the leading operator that breaks the symmetry for rotations ofφ from U(1) down to Z 4 is dangerously irrelevant: irrelevant at the critical point but relevant within the VBS phase [2]. This leads to the prediction of a crossover between U(1) and Z 4 symmetry within the VBS phase, on a length scale ξ VBS which is parametrically larger than the correlation length: ξ VBS ∼ ξ 1þjy 4 j=3 , where y 4 < 0 is the RG eigenvalue of the fourfold anisotropy [50]. This has been confirmed in the J-Q model [7,51]. Figure 4 gave visual evidence for the emergent U(1) symmetry in the loop model. A quantitative measure of Z 4 anisotropy is hcos 4θi, whereφ ¼ jφjðcos θ; sin θÞ. Figure 12 shows data for sizes up to L ¼ 200. Ignoring scaling violations, the anisotropy should behave as [50] hcos 4θi ¼ f(L 1=ν 4 δJF 2 ðδJÞ); ð15Þ where F 2 ðxÞ ¼ 1 þ ax þ bx 2 takes into account nonlinear dependence of the scaling variable on J (which is needed here because of the larger range of δJ studied for this observable) and ν 4 ¼ νð1 þ jy 4 j=3Þ. The inset to Fig. 12 shows the attempted scaling collapse using Eq. (15). The exponent ν 4 ¼ 1.09ð6Þ is obtained from the fit. This confirms the irrelevance of fourfold anisotropy to the behavior at the transition and to the explanation of the scaling violations. The corresponding value of jy 4 j is dependent on the assumed value of ν but is considerably larger than the estimate in Ref. [7] for a variant of the J-Q model. The closeness of the finite-size effective values of η VBS and η Néel in Fig. 8 and Eq. (12) makes it tempting to speculate about a much larger emergent symmetry-an SO(5) symmetry relating the Néel and VBS vectors. This can be incorporated into an alternative field theory for the deconfined critical point [52,53], which was argued to be equivalent to Eq. (1) [54]. This symmetry enhancement would be analogous to the emergent SO(4) symmetry of the 1D spin-1=2 chain, which relates the spin-Peierls order parameter and the Néel vector [55]. In the future, it would be interesting to check explicitly for SO(5) symmetry. [This has now been addressed in [57].]

C. Néel stiffness and spanning strands
A useful observable is the spanning number N s , defined as the number of strands that span the system in (say) the z direction. Its mean value hN s i may be taken as a definition of the stiffness of the Néel order parameter [58]. At a conventional critical point, hN s i has scaling dimension zero and the scaling form Therefore, hN s i should be a universal constant at a critical point, modulo corrections due to irrelevant scaling variables: Plots of hN s i versus J for different L should cross at J c . In the VBS phase, hN s i tends to zero exponentially in L, and in the Néel ordered phase, it grows as L.  The mean spanning number was shown in Fig. 5. Contrary to the above expectation, hN s i appears to diverge slowly with system size at the critical point. This is manifested in the upwards drift of the crossing points in the main panel. In the inset, we show pseudocritical values N Ã s ðLÞ defined in two different ways. The data cannot be fitted with conventional scaling corrections from an irrelevant variable, i.e., N Ã s ¼ N crit s − AL y with negative y: Attempting such a fit leads to a positive (relevant) y.
Previous work on the SU(2) J-Q model found a drift in a closely related winding number and proposed that this indicated logarithmic corrections to scaling [9]. Similar drifts were found for the SU(3) and SU(4) J-Q models [28], fitting slightly better to a power law than a logarithm. The larger sizes considered here for the SU(2) case show that the divergence is certainly stronger than logarithmic. On attempting to represent it by a pure power law N Ã s ∼ AL a , we find that the exponent drifts upwards, but a power law plus constant fits the results for all L. This divergence is of course still slower than the linear behavior expected asymptotically at a first-order transition.

Drift of critical probability distribution
In addition to the mean hN s i, we examine the full probability distribution of N s . Let P k be the probability that N s is equal to 2k, meaning that k oriented strands span the system in a specified direction and k in the reverse direction. This again has scaling dimension zero, so conventionally we would expect the scaling form By contrast, at a first-order transition, where the short-loop and infinite-loop phases coexist, P k would have a peak at k ¼ 0 from the short-loop phase and a peak at k ∝ L from the infinite-loop phase. The distribution P k obtained numerically is shown in Fig. 13, for various L. To compare different sizes, we tune J for each L so that P 0 ¼ 0.3 (using the Ferrenberg method [43]). For comparison, the inset shows the distribution at fixed J ¼ 0.0885, very close to the critical point. Contrary to Eq. (18), the data show no sign of tending to a universal distribution. On the other hand, we do not see a double-peaked structure developing either.
The scaling form Eq. (18) would imply that P k is a universal function of hN s i for each k. [Explicitly, P k ¼ ðg k ∘h −1 ÞðhN s iÞ.] Therefore, a plot of, e.g., P 1 against hN s i would show scaling collapse without the need to adjust any parameters. See Ref. [38] for successful examples of such scaling collapse for the compact CP 1 [i.e., Oð3Þ] and compact CP 2 models. It is clear from Fig. 13 that such a collapse will not work here. Figure 14 shows this for P 1 (interpolating curves are obtained with the Ferrenberg method [43]). The dramatic failure to collapse is quantified in the inset, which shows the maximum value of P 1 as a function of L. At a conventional critical point, this would reach a finite constant as L → ∞, while at a first-order transition, it should tend to zero. There is no evidence for saturation over this size range, though it cannot be ruled out.
We note that the maximum of P 1 decreases below the universal value for the O(3) universality class [38], which is close to 0.4 [a stiffness in the J-Q model also drifts beyond the O(3) value [17]]. This is further confirmation that we are dealing with a direct transition rather than two separate transitions that are too close to be resolved.

D. Energy distribution
The behavior of the heat capacity, proportional to the variance in the energy E, has been discussed in connection with Fig. 6. Further information is contained in the full probability distribution for E [which is defined in Eq. (4)]. The top panels of Fig. 15 show how this distribution evolves as J is varied. The bottom panel shows the distribution at J ≃ J c for various L. We do not see a double-peaked distribution. The width of the critical distribution also decreases with increasing system size, contrary to the expectation for a first-order transition (Fig. 6, lower inset).
The Binder parameter plotted in the inset to Fig. 16, is an alternative quantity for analysis. The data show a peak near the transition, with a height V max that decreases with L (it is necessary to use the Ferrenberg interpolation method [43] for accurate estimates of V max ). At a first-order transition, V max should saturate to a constant as a result of the double-peaked energy distribution, while at a continuous transition with α > 0, the peak height should tend to zero as V max ∼ L 2=ν−6 . Here, a direct estimate of ν using the slope gives a value that drifts from ν ∼ 0.621ð5Þ for system sizes 32 ≤ L ≤ 64 to ν ∼ 0.481ð12Þ for L ≥ 256. However, in addition to the peak, V has a large background contribution scaling as L −3 (see inset). It is natural to subtract such a correction. This allows V max to be fitted to the power-law form corresponding to ν ¼ 0.468ð6Þ-see Fig. 16, main panel.

E. Correlation-length exponent
In this section, we discuss two different approaches to determining the correlation-length exponent ν, one relying on finite-size scaling and one not (see Sec. IVA for anomalous dimensions η Néel , η VBS ).
First, we obtain estimates using the standard finite-size scaling forms for various observables. For example, the Binder cumulant for the VBS order parameterφ would naively scale as U VBS ¼ f U ðL 1=ν δJÞ, so the maximum value of dU VBS =dJ should grow as L 1=ν . Therefore, we can define an effective exponent via ν eff ðLÞ −1 ¼ d logð½U 0 VBS max Þ=d log L. We calculate such numerical derivatives using four consecutive system sizes. Figure 17 shows the resulting estimates ν eff ðLÞ from U VBS , from the probability P 0 of having no spanning strands, and the order parameters N and jφj. ν eff ðLÞ drifts from large values, around 0.62 (in accordance with previous studies [5][6][7][8]10]) to values around 0.46. The latter is in agreement with the estimate from the heat capacity with  the background subtracted (Sec. III). Despite the drift, ν eff ðLÞ remains above the unitarity bound 2=5. A drift in ν was also identified in Ref. [10].
The above estimates all rely on finite-size scaling forms for observables on the scale of the system size. To avoid the assumption of conventional finite-size scaling, we also estimate ν directly from the correlation length ξ in the regime where ξ ≪ L. For values of J in the VBS (short-loop) phase, we determine ξ by fitting the L dependence of the spanning number to the expected form N s ∝ L 2 e −L=ξðJÞ . In the Néel phase, the spanning number is expected to grow as N s ∼ AL=ξðJÞ. This allows us to determine ξðJÞ up to the overall constant A.
The results are shown in Fig. 18. The power-law fits shown give ν ¼ 0.477ð4Þ for the data in the Néel phase and ν ¼ 0.503ð9Þ for the data in the VBS phase. These values are close to the estimates in Fig. 17 at the largest sizes. But it is remarkable that here this behavior sets in at much smaller length scales. For the VBS phase [where we can determine ξðJÞ without the complication of an overall constant], the above exponent fits the data well starting from scales as small as ξ ∼ 15.
We believe that if the transition is indeed continuous, the correlation-length exponent is close to ν ¼ 0.5, thus considerably smaller than most earlier estimates from J-Q and NCCP 1 models. This value is close to 2 þ ϵ expansion results for the CP 1 nonlinear sigma model, which should apply to the deconfined critical point, assuming the transition is continuous (see Sec. VI D).

V. FIELD THEORY FOR LOOP MODEL; HEDGEHOG FUGACITIES
Models for completely packed, oriented loops can be mapped to lattice CP n−1 models with an unconventional but simple form [37,38]. At first sight, the continuum limit of these models is simply the (compact) CP n−1 sigma model. Here, we discuss this continuum limit in more detail and show that hedgehog defects in the CP n−1 spin configuration contribute imaginary terms to the action that are analogous to the Berry phases in the Euclidean action for the 2D quantum Heisenberg model [34,35]. These terms are crucial for the present model. By the reasoning of Refs. [1][2][3][4], they change the effective continuum description from the compact CP n−1 model to NCCP n−1 . (By contrast, the imaginary terms were unimportant for the transitions discussed in Refs. [37,38] as a result of the lower lattice symmetry there.) The quantities of interest are determined by symmetry, so it is enough to consider the case J ¼ 0, where the lattice field theory for the loop model is simplest. The CP n−1 spins are placed on the links l of the lattice. They are complex vectors z l ¼ ðz 1 l ; …; z n l Þ, with fixed length jzj 2 ¼ 1 and the gauge redundancy z l ∼ e iφ l z l . In a loose notation where the incoming links at a given node are denoted i and i 0 and the outgoing links o and o 0 , the partition function is Here, "Tr" is the integral over the zs. Under a gauge transformation of z l , the terms for the two nodes adjacent to l pick up opposite phases, so the Boltzmann weight is invariant. The mapping between Eq. (20)  model follows from a straightforward graphical expansion which is described in Ref. [38]. Let us consider the continuum description of Eq. (20). To begin with, take a configuration in which z is slowly varying. Each term in the product over nodes is then close to 1, and we may obtain a continuum sigma-model Lagrangian by a derivative expansion. In 3D, the only term with two derivatives allowed by global, gauge, and lattice symmetries is the standard sigma-model kinetic term. Let us focus on the n ¼ 2 case and parameterize CP 1 (which is simply the sphere) using the Néel vector instead of the gauge-redundant field z. Then, A crude way to estimate a bare value of K is to calculate the Boltzmann weight in Eq. (20) for a spin configuration with a uniform twist, giving K ¼ 1=16. (For general n, L σ may be written in terms of the matrix Q ¼ zz † − 1=n.) The Lagrangian obtained by the derivative expansion can fail to capture the true scaling behavior in two ways. It fails in a trivial way whenÑ varies strongly at a node. This will of course be the case in the lattice model, and it leads to an order-one renormalization of the stiffness.
More importantly, the phase of z can vary rapidly even if N is slowly varying. Nodes whereÑ is approximately constant but where this phase varies abruptly contribute imaginary terms to the action. For smooth configurations of trivial topology, these phases cancel. However, in the presence of hedgehog defects, there remain nontrivial phases that are missed by the derivative expansion. This is because in a configuration with a hedgehog, it is impossible to find a gauge in which z is everywhere slowly varying, even far from the hedgehog core. This follows from the fact that topological flux density B μ , which when integrated over a closed surface gives the signed number of hedgehogs inside, is a total derivative when written in terms of z: B μ ¼ ð1=iÞϵ μνλ ∇ μ ðz † ∇ ν zÞ. If z were continuous, integrating the topological density over a large sphere would give zero. Therefore, z must be discontinuous somewhere on the sphere if the sphere encloses a hedgehog.
A simple calculation is required to determine what effect the imaginary terms have on the weight for a configuration with a hedgehog. We do this calculation in Appendix B. For specificity, we take the hedgehog to be centered on a site of C 1 or a site of C 2 -for example, at the center of the cube in Fig. 1. These locations form a body-centred cubic (bcc) lattice, with four sublattices. We find that the weight of a configuration with a hedgehog acquires a fugacity proportional to depending on which sublattice it sits on. (More precisely, only the relative phase between different locations is meaningful [60].) This result also generalizes immediately to larger n. The result matches nicely what is found for the 2D quantum Heisenberg model [34][35][36], where the fugacity for instantons-hedgehogs in spacetime-takes the same set of values as above depending on which of the four sublattices of the square lattice the instanton occurs on.
By symmetry, we infer that the coarse-grained hedgehog fugacity vanishes. In other words, it vanishes as a result of phase cancellation between configurations in which the hedgehog is centered on nearby sites on different sublattices. Thus, the arguments of Refs. [1][2][3][4] apply, giving the NCCP 1 model [Eq. (1)] as the continuum description.
An alternative argument for the NCCP 1 description of square lattice spin-1=2 antiferromagnets was given in Ref. [61], focusing on the VBS order parameterφ and its vortex defects rather than the Néel order parameterÑ and its hedgehog defects. The key point of this alternative argument is that a vortex in the VBS order parameter carries a single unpaired spin at its center. In spacetime, this corresponds to an extended spinon worldline running along the vortex core. This has a direct interpretation in the loop model: A vortex line in the node order parameterφ has a single extended loop running along its core. We have confirmed this explicitly by constructing such configurations.
In previous work, we have considered transitions in a different version of the loop model which does not preserve the full lattice symmetry [37,38]. The transitions in that less-symmetric model are described by the compact CP n−1 model, unlike the present loop model whose transition is described by NCCP n−1 . This is because breaking lattice symmetry spoils the cancellation between the values in Eq. (23), leaving a nonzero hedgehog fugacity. This result allows the standard critical behavior of the compact CP n−1 model [i.e., of the usual O(3) model when n ¼ 2]. We note that in the case n ¼ 3, the compact CP 2 model appears to show an interesting continuous transition that is naively forbidden by Landau theory [37,38,62]. A RG explanation for why the expectation from Landau theory breaks down in this case was given in Ref. [38].

VI. RG FLOWS IN THE NCCP n−1 MODEL (n-COMPONENT ABELIAN HIGGS MODEL)
In this section, we make a conjecture for the topology of the RG flows in the NCCP n−1 model, the n-component generalization of Eq. (1) with z ¼ ðz 1 ; …; z n Þ:

041048-13
This generalization is also known as the n-component Abelian Higgs model. We treat both n and the spatial dimension d (between 2 and 4) as continuously varying. Scaling violations are seen in a wide variety of different lattice models that are related to this field theory (at n ¼ 2) and persist to very large length scales [63,64], so we believe a plausible explanation for them should appeal to universal physics of the NCCP n−1 model and not to accidental features of specific Hamiltonians. Results for SU(3)-and SU(4)-symmetric models (n ¼ 3, 4 [10,28]) suggest that a satisfactory explanation should also account for scaling violations across a range of n. Figure 19 shows the basic topology of the RG flows we find. This is a sheet of RG fixed points projected on the space of n, d, and a scaling variable λ. [Close to 4D, λ is the quartic coupling in Eq. (1), but in lower dimensions, one cannot make this identification.] The RG flow is parallel or antiparallel to the λ axis since n and d do not flow. λ is irrelevant on the critical sheet and relevant on the tricritical sheet. The strongly relevant coupling that drives the transition (i.e., the mass) is not shown since we consider the theory at the critical value. For any fixed value of n between zero and n Ã4 ≃ 183, the critical point exists so long as we are sufficiently close to two dimensions. When d is increased, the critical point disappears at a universal value d Ã ðnÞ, by merging with the tricritical point.
To begin with, consider three limits in which the NCCP n−1 model is solvable. First, it is tractable by saddle-point at large n, where it yields a nontrivial critical point for 2 < d < 4. This critical point describes a direct transition between a Higgs phase, where z is condensed and SUðnÞ symmetry is broken, and a Coulomb phase where SUðnÞ symmetry is unbroken and the gauge field A is massless.
The field theory is also tractable in a 4 − ϵ expansion [65]. For infinitesimal ϵ, a weak-coupling critical point exists only if n is greater than or equal to a value that we denote n Ã4 . This value is quite large, n Ã4 ≃ 183. In fact, in the regime n > n Ã4 where the critical point exists, the 4 − ϵ expansion also yields a tricritical point at a smaller value of the quartic coupling λ. As n approaches n Ã4 from above, these two fixed points approach each other, and they annihilate when n reaches n Ã4 . For n < n Ã4 , there is no nontrivial fixed point: The theory is expected to flow to a discontinuity fixed point at large negative λ representing a first-order transition.
Finally, the NCCP n−1 model can be studied in 2 þ ϵ dimensions by switching from a soft-spin formulation to a nonlinear sigma model [66][67][68]. In this regime, a continuous phase transition is found for all values of n greater than zero. (The "replicalike" regime n ≤ 1 is meaningful and describes certain classical loop models [37,59,69].) In the 2 þ ϵ approach, it does not matter whether the nonlinear sigma model is formulated with a dynamical U(1) gauge field or as a pure nonlinear sigma model with target space CP n−1 : The two formulations give identical results [66]. (When the dynamical gauge field is included, its coupling flows to infinity, so it can be integrated out, leaving the usual CP n−1 nonlinear sigma model.) A crucial point, made in Ref. [32], is that the fixed point found in the large n approach is the same as that found in both the 2 þ ϵ and 4 − ϵ expansions. This can be seen by comparing the results for the critical exponents in the regions of overlap of the expansions. Viewing n and d as continuous variables, there is therefore a continuous family of fixed points in a region of the (n, d) plane for 2 < d < 4 and sufficiently large n [32]. This region is defined by n > n Ã ðdÞ, where n Ã ðdÞ is the d-dependent value of n at which the fixed point disappears. From the 4 − ϵ approach, we know that the limiting behavior as d → 4 is n Ã ðdÞ → n Ã4 ≃ 183 [70], and from the 2 þ ϵ expansion, we know that n Ã ðdÞ → 0 as d → 2. Using the sigma model, we may also argue that the slope of n Ã ðdÞ is finite as d → 2 (see the endnote [72]). The value of n Ã ð3Þ is not known, but it is possible to show that n Ã ð3Þ > 1 (Sec. VI C). We assume that n Ã ðdÞ increases monotonically with d.
We can go further by noting that in the 4 − ϵ expansion, the way in which the nontrivial fixed point disappears at n Ã ðdÞ is by annihilation with a tricritical point. On the basis of continuity, we expect that the mechanism for the disappearance of the fixed point at n Ã ðdÞ is the same for all d. This leads directly to the topology in Fig. 19. It will be convenient to denote the value of d where the merging of the critical and tricritical points happens, for a given n, by d Ã ðnÞ.
A. RG flows close to the merging line is a broad range of n values where the line lies close to three dimensions, i.e., where jd Ã ðnÞ − 3j is small. So it is worth studying the RG flows in this regime.
On the line, λ is marginal. After rescaling and shifting λ by a constant, its RG equation is Moving slightly away from the line, the RG equation becomes, to lowest order in d Ã ðnÞ − d, , which is zero at n ¼ n Ã ð3Þ (where the critical point disappears in 3D) and small over the range of n where the merging line lies close to d ¼ 3: When n < n Ã ð3Þ, ΔðnÞ is negative and there is no fixed point in 3D. Instead, the RG flows go to large negative λ, suggesting a first-order transition. However, if we are close to the line, so jΔðnÞj is small, the RG flow becomes very slow. Integrating Eq. (27) shows that this implies an exponentially large correlation length at this transition:

B. Interpretation
At first sight, the topology we have found for the RG flows suggests two possible explanations for the scaling violations.
First, if we (speculatively) assume that n Ã ð3Þ > 4, but that the merging line lies close to d ¼ 3 over the range 2 ≤ n ≤ 4, then we obtain an anomalously weak first-order transition for n ¼ 2, 3, 4. In this scenario, there is pseudocritical behavior (with drifting exponents [75]) up to the exponentially large length scale of Eq. (28), thanks to the "nearby" fixed point at slightly smaller spatial dimension. The virtues of this scenario are that it appeals to universal features of the RG flow and so may explain why numerous different lattice models see very similar scaling violations, and that it can produce scaling violations over a range of n.
(Note that if jΔj is very small, there exists a range of sizes where λ appears to be a conventional marginally irrelevant variable. However, this range ends at a size L Ã which is parametrically smaller than ξ [76] when ξ is large. The basic point is that the stages of the RG flow with λ > 0 and with λ < 0 take roughly equal amounts of RG time, but the latter corresponds to a vastly larger length scale because of the logarithmic relationship between length and RG time.) Second, we might try to explain the scaling violations differently by postulating that in 3D the values n ¼ 2, 3, 4 lie just below the merging line [n Ã ð3Þ < 2] so that ΔðnÞ is small and positive for n ¼ 2, 3, 4. This would give a true critical point with large (but conventional) scaling corrections due to a small irrelevant exponent y irr . However, our numerical results strongly indicate that a small y irr is not sufficient, on its own, to explain what we see.
We emphasize that since our argument fixes the topology of the RG flows but not the numerical value of n Ã ðdÞ, the scenario above for a weak first-order transition is speculative. In Sec. VII, we discuss another possible conjecture.

C. Bound on n Ã ð3Þ
The value of n Ã ð3Þ is not known, but we can argue that This may look surprising at first glance since the singlecomponent Abelian Higgs model, n ¼ 1, certainly has a continuous phase transition in 3D. This transition is related by duality to that of the XY model [77]. Equation (30) means that this "inverted XY" phase transition does not lie on the sheet we are considering: It is not analytically connected to the deconfined critical point at large n [78]. Formally, one can see this as follows. If the inverted XY transition did lie on the critical sheet of Fig. 19, we could describe it by setting n ¼ 1 in the 2 þ ϵ expansion of the CP n−1 nonlinear sigma model. But this is evidently not the case. The inverted XY transition is a conventional thermodynamic phase transition with nontrivial signatures in the free energy. In contrast, the sigma model at n ¼ 1 is a replicalike theory, in which the number of degrees of freedom becomes zero and the free energy vanishes identically. The same reasoning implies that n Ã ð3Þ is strictly greater than 1. Otherwise, the n ¼ 1 model would have a Higgs transition with an unphysical replicalike continuum description.
Instead of being connected to the critical points in Fig. 19, the inverted XY transition is the n ¼ 1 limit of a much simpler transition-namely, that which (for n > 1) separates the Higgs phase, where z is condensed, from a pair condensed phase where the bilinear zz † is condensed but z itself is not. (This phase appears for appropriate couplings [19,79].) SUðnÞ symmetry is broken on both sides of this transition, so it is not like the critical points in Fig. 19. One can check that the transition is in the inverted XY universality class for all n [as the interactions between the critical sector and the Goldstone modes of the broken SUðnÞ symmetry are irrelevant], so the n dependence of the critical behavior is trivial [80].

D. "Failure" of the 2þϵ expansion of the O(3) model
The CP 1 nonlinear sigma model is the O(3) nonlinear sigma model by another name, as the target space is the sphere (CP 1 ¼ S 2 ). Therefore, the topology of the flows in Fig. 19  . The applicability of the 2 þ ϵ expansion to the Abelian Higgs model was also argued in Ref. [32].
Although this is contrary to what is often assumed, it should not be surprising in the light of knowledge about hedgehogs in the 3D O(3) model [4,81]. Suppressing these topological defects has been convincingly argued to change the critical behavior [4]. It is natural that the 2 þ ϵ expansion, which considers only spin waves, describes the behavior in the absence of hedgehogs [82].
The conclusion is also supported by a RG approach to the OðMÞ model that employs a double expansion in (d − 2) and (M − 2) [33]. This shows that the standard 2 þ ϵ expansion fails to capture the critical behavior of the OðMÞ model when M is smaller than a d-dependent critical value M c . To first order in d − 2, the relationship is ðM c − 2Þ ≃ ðπ 2 =4Þðd − 2Þ, which gives M c ∼ 4.5 in 3D. While higher-order corrections in (d − 2) could be significant, this result suggests M c > 3 [33].
Reference [33] also notes the poor agreement between the 2 þ ϵ exponents in 3D and the exponents of the O(3) model in 3D. For example, ν is equal to ν ¼ 1=2 at order ϵ 2 , or to 2=5 at order ϵ 3 [67,84]. In the light of the preceding, we should instead apply these exponents to the SU(2) deconfined transition. Indeed, the ϵ-expansion values ν ¼ 1=2 and ν ¼ 2=5 are remarkably close to our best estimates of ν (see Sec. IV E).

VII. INTERPRETATION AND CONCLUSIONS
We have argued that the scaling violations at the deconfined critical point are too severe to be explained as corrections to scaling from a weakly or marginally irrelevant scaling variable. This clearly excludes the most conventional scenario for the first time. We have also sharpened the possible alternatives as follows.
Some of our numerical results suggest that estimates of the exponents ν, η Néel , and η VBS may be better defined if we avoid using finite-size scaling to obtain them (see, in particular, Secs. IVA 1 and IV E)-although, of course, abandoning finite-size scaling restricts us to length scales much smaller than the system size. The most radical conjecture would therefore be to attribute the scaling violations to a dangerously irrelevant variable (DIV) which leaves critical behavior intact but modifies finite-size scaling (see also Ref. [28]). In the simplest picture, the role of this DIV would be to cut off fluctuations of some zero mode(s) of the fields that are unbounded in the pure fixed-point theory.
The main examples we know of this phenomenon are in free theories (such as ϕ 4 theory above 4D [31] or the quantum Lifshitz theory [85]); another type of example is in theories that are dual to free theories (see endnote [86]). Further work is required to determine whether it is a plausible possibility in an interacting theory like Eq. (1).
The alternative possibility of an anomalously weak firstorder transition has been discussed in detail in Sec. VI. We have shown that, in principle, there is a mechanism by which a very large correlation length can appear without the need for fine-tuning of the Hamiltonian and that in this scenario there would be pseudocritical behavior over a large range of scales, with (for example) drifting critical exponents. Such a possibility is hard to exclude numerically. We note, however, that we do not see the usual signs of an incipient first-order transition, despite studying much larger scales than the early simulations used to argue for first-order behavior in the J-Q model [16]. The first-order scenario would also leave the good scaling of, for example, the derivatives of the correlators (Sec. IVA 1) a mystery.
Intriguing questions therefore remain for the future. The loop model is an ideal platform for further work on the deconfined transition, since it provides an intuitive geometrical picture and since isotropy in three dimensions is a convenient feature. It would also be interesting to perform simulations at other values of n [which in the formulation just after Eq. (3) need not be an integer] in order to probe the scenario of Sec. VI.
Various modifications to the loop model are possible. For example, one may allow a third node configuration (see Fig. 2) in which the two incoming links are joined and the two outgoing links are joined, so that the loops are no longer consistently oriented. The symmetry is then broken from SUðnÞ to SOðnÞ, and the relevant field theories are RP n−1 models [37,38,46,59] which can show a Z 2 spin liquid phase [87,88]. One may also study the effects of various anisotropies and symmetry-breaking perturbations. Finally, the loop model generalizes to four dimensions, where it may be a useful tool to search for new types of critical behavior arising in 3 þ 1D quantum magnets. from the Gordon and Betty Moore Foundation under the EPiQS initiative (Grant No. GBMF4303). M. O. thanks KITP for hospitality. We thank I. Klebanov for bringing Ref. [71] to our attention.

APPENDIX A: METHODS
The numerical procedure is as follows. An initial state is constructed by choosing at random one of the two possible configurations of a node with equal probabilities. Loops are then formed by following turning instructions at each node. We generate the fugacity n from a sum over loop colors, assuming n is an integer. So a color is associated with each loop chosen with equal probability from n alternatives.
Subsequent states are generated using three kinds of parallelized Monte Carlo moves to ensure that at equilibrium configurations are distributed according to the partition function, Eq. (4). The first one updates the state of the nodes using a checkerboardlike algorithm. Each node, if its four links have the same color, changes its state with probability min fexp ½−2JσðiÞ P j nn σðjÞ; 1g. The updates of all the nodes in the lattice are done in three stages. In the L lattice, there are two sublattices A and B. Each sublattice is tripartite, and the three sublattices are simple cubic types. At each time, one sublattice of A and one sublattice of B are updated. The second type of Monte Carlo move chooses a link at random and changes the color of all the links of the associated loop to a different color, chosen with uniform probability from the n − 1 possibilities. This move is also parallelized by letting each thread choose a link and change the color of the loop if it is not already visited. The third type of move is to recolor all loops in the system, with the new colors selected independently and at random for each loop. It is designed to ensure that the colors of short loops equilibrate efficiently.
A number of these moves are combined to form a composite update which we term a Monte Carlo sweep. The ingredients in a single sweep are as follows. First, we iterate 20 times a sequence in which moves of the first two types are intercalated and repeated for each of the three sublattices. Then, we apply the third type of move. Measurements are performed once every Monte Carlo sweep. The autocorrelation function of the energy is used to estimate a correlation time. The blocking and bootstrap methods [89] are used to estimate errors.
As an interpolation scheme, we have used Ferrenberg's multiple histogram method [43] to obtain a continuous set of values as a function of J of the different quantities. This method is also employed whenever a derivative has to be calculated and a maximum or minimum has to be obtained. Errors related to this technique are calculated using the bootstrap method.
We consider system sizes of up to 3.9 × 10 8 links or lateral size L ¼ 640, with extensive results for L ≤ 512. The minimum number of Monte Carlo sweeps used is 10 5 for any J and L, and it increases with decreasing L.

APPENDIX B: CALCULATION OF HEDGEHOG FUGACITY
Regarded as a lattice magnet for classical CP n−1 spins, the partition function (20) has the peculiar feature that although it is both local and gauge invariant, it is not simply expressed in terms of the gauge-invariant quantityÑ (for n ¼ 2) or Q (for larger n). The consequence of this, as noted in Sec. V, is that the sigma-model action arising from a derivative expansion may need to be supplemented by purely imaginary terms from nodes at which the phase of z changes abruptly. In the presence of hedgehogs, such nodes are inevitable, even far from the hedgehog core. As a result, the hedgehog fugacity acquires a spatially varying phase.
For simplicity, we take the hedgehogs to sit at the center of a cube of C 1 or C 2 [Eq. (2)], for example, at the center of the cube in Fig. 1 (left). These locations form a bcc lattice. We take the origin at one such bcc site and the coordinate axes parallel to the links.
Focusing on the case n ¼ 2 (the generalization to larger n is immediate), let us first consider the representative configuration in which the hedgehog is centered at the origin and the Néel vectorÑ (defined in Sec. V) points in the radial direction. In polar coordinates, 0 ≤ θ ≤ π, 0 ≤ φ < 2π, this is N ¼ ðsin θ cos φ; sin θ sin φ; cos θÞ; ðB1Þ where the coordinates of a link are those of its midpoint. To write this in terms ofz, we must pick a gauge. It is convenient to choose one in which the Boltzmann weight is approximately equal to 1 for as many nodes as possible.
For the links with positive and negative z coordinates (θ < π=2 and θ > π=2), we take, respectively, z ¼ðcosθ=2;e iφ sinθ=2Þ; z ¼ðe −iφ cosθ=2;sinθ=2Þ: ðB3Þ We see from Fig. 1 that there are also links in the equatorial plane. For these links, we take Now consider e −S node . In fact, it will suffice to consider only nodes far from the core and to treat z as constant except for the discontinuities in our gauge choice: Other contributions to the action are either included in the spatially independent amplitude of the hedgehog fugacity or are already captured by the naive derivative expansion. With the above gauge choice, the nodes at which z varies abruptly all lie in the equatorial plane and have two of their links within this plane, one above it and one below. Figure 1 shows four such nodes (all black). From Eqs. (B2)-(B4), we find that most of these nodes have e −S node ≃ 1, despite the variation in the phase of z. However, there is a string of nodes along the positive x axis (φ ¼ 0, θ ¼ π=2), each of which contributes a minus sign to the Boltzmann weight (i.e., e −S node ≃ −1 for these nodes). If we translate the core of the hedgehog by the vector (2, 0, 0), while keeping the configuration far from the core fixed, we change the number of nodes on this string by 1. Therefore this translation changes the sign of the Boltzmann weight.
Let the phase term in the hedgehog fugacity be denoted expðiηðrÞÞ; where the spatial vector r lies on a bcc site. [An "antihedgehog" of negative topological charge has phase factor e −iηðrÞ , as we see from Eq. (B2) and the fact that complex conjugating z exchanges hedgehogs and antihedgehogs.] The phase ηðrÞ is defined only up to a constant: For example, in a closed system with periodic boundary conditions, there are equal numbers of hedgehogs and antihedgehogs, so the constant part of η drops out. It may be seen from Fig. 1 that the translational symmetry between bcc sites is not spoiled by the link orientations. Using this, we may argue that ηðrÞ is of the form ηðrÞ ¼ k:r for some momentum k. By the above calculation, e iηðrÞ ¼ −e iηðrþ2xÞ , wherex ¼ ð1; 0; 0Þ. By symmetry, we have similar results for translations in the y and z directions. This is enough to fix k up to a sign: One of these signs applies to the hedgehog and one to the antihedgehog. We have not fixed which is which, but it does not matter. With this k, the hedgehog fugacity takes four distinct values on the four sublattices of the bcc lattice, proportional to AE1 and AEi. This is the result quoted in Sec. V.