Consistent Scaling Exponents at the Deconfined Quantum-Critical Point^*

Among the many proposed exotic quantum states and quantum phase transitions beyond the Landau–Ginzburg–Wilson (LGW) paradigm in two dimensions,^[1–3] the deconfined quantum critical point (DQCP)^[4] is special because it has concrete lattice realizations in sign-free "designer models" accessible to quantum Monte Carlo (QMC) simulations.^[5] Indeed, the first hints of an LGW-forbidden continuous transition between antiferromagnetic (AFM) and spontaneously dimerized valence-bond solid (VBS) ground states came from QMC simulations,^[6] and following the DQCP concept (which builds on previous works on VBS phases and topological defects^[7–13]) numerous additional studies have been reported. The most compelling results for DQCP physics have been obtained with J–Q models,^[14–27] which are Heisenberg antiferromagnets in which the exchange of strength J is supplemented by multi-spin couplings of strength Q that induce singlet correlations and eventually cause spontaneous dimerization. The space-time loop structure employed in QMC simulations^[28–30] can also be used to formulate analogous classical three-dimensional loop models, which exhibit behaviors very similar to the J–Q models.^[31,32]

Even though very large lattices have been studied, with linear size L up to 256 for the J–Q model^[18,26] and twice as large for the loop model,^[32] it has not yet been possible to draw definite conclusions on the nature of the AFM–VBS transition. While no explicit signs of a first-order transition have been detected in the best DQCP candidate models (in contrast to intriguing discontinuous transitions with emergent symmetry in related models^[33–35]), some observables exhibit scaling behaviors incompatible with conventional quantum criticality. Such behaviors have been interpreted as runaway flows toward what would eventually become a first-order transition on lattices even larger than those studied so far.^[16,22,36] Another proposal is that the DQCP is even more exotic than initially anticipated, with novel relationships between critical exponents originating from the presence of two divergent length scales^[26]—in addition to the standard correlation length ξ, there is a larger scale ξ' associated with a "dangerously irrelevant" perturbation and emergent U(1) symmetry of the near-critical VBS fluctuations in the DQCP scenario.^[4,37] The weak first-order scenario has attracted attention in the context of non-unitary conformal field theories (CFTs), which have critical points in the complex plane.^[38–41] In this scenario, the AFM–VBS transition is a "walking" first-order transition,^[42,43] where the renormalization-group flow (which is also manifested in finite-size scaling) is affected by the inaccessible nearby critical point and only slowly "walks" to a first-order instability.

In support of the weak first-order scenario, the J–Q and loop models are often invoked as supporting evidence, though there are no concrete predictions that have been compared with the numerical results. In the absence of any quantitative tests or clear signs of discontinuities or coexistence state in the lattice models, the walking scenario should not be accepted as the final word on the fate of the DQCP. Here we will show that the cited^[39–41] large scaling corrections affecting estimates of the correlation length exponent ν^[21,26,32] are not precursors to a first-order transition. We reach this conclusion by extracting ν from the scaling dimension of the relevant field of the model. The corresponding correlation function exhibits only small scaling corrections and delivers an exponent compatible with results based on Binder cumulants; ν = 0.455(2). Given the well behaved estimators of ν, a continuous transition is the most likely scenario.

To set the stage, we briefly summarize some standard facts on critical scaling. Consider a Hamiltonian H_c tuned to a quantum critical point to which a perturbation is added that maintains all the symmetries of H_c;

$$\begin{eqnarray}&&H={H}_{{\rm{c}}}+\delta \displaystyle \sum _{r}D(r),\end{eqnarray} \tag{ 1 }$$

where r denotes the lattice coordinates and D(r) are local operators. Normally H is written in a form with some tunable parameter g such that, for some critical value g = g_c, H(g_c) = H_c and δ = g – g_c. We assume that the system develops long-range order when δ > 0, with an order parameter m(r) such that 〈m〉 = 〈m(r)〉 ∝ δ^β for small δ > 0 and m = 0 for δ < 0. The critical exponent β depends on the universality class of H_c in the thermodynamic limit. On either side of the phase transition, the exponential decay of the correlation function C_m(r) = 〈m(0)m(r)〉 – 〈m〉² defines the divergent correlation length, ξ ∝ |δ|^−ν. At δ = 0, the correlation function takes the critical form C_m(r) ∝ r^−2Δ_m, where Δ_m = β/ν is the scaling dimension of the operator m.

In QMC calculations ν is typically extracted using finite-size scaling of some dimensionless quantity, such as the Binder ratio R = 〈M⁴〉/〈M²〉², where M = ∑_r m(r). Neglecting scaling corrections, in the neighborhood of the critical point we have R(δ,L) = R(δL^1/ν), by which ν (and the critical point g_c if it is not known) can be obtained from data for different values of δ and L. A less common method is to use the relation 1/ν = d – Δ_D, where d is the space-time dimensionality (here d = 3) and Δ_D is the scaling dimension of the perturbing operator D in Eq. (1). The scaling dimension can be obtained from the power-law decay C_D(r) ∝ r^−2Δ_D of the correlation function C_D(r) = 〈D(0)D(r)〉 – 〈D〉² at g_c.

It is not clear to us why ν is not commonly extracted from C_D(r), but there are two potential drawbacks: (i) Often Δ_D is rather large, e.g., in the case of the O(3) universality class (of which we will show an example below) Δ_D ≈ 1.6, so that the correlation function decays rapidly and is difficult to compute precisely (with small relative statistical errors) at large r. (ii) The operator D is often off-diagonal and may appear to be technically difficult to compute. However, although the latter issue is absent in simulations of classical systems, the scaling dimension Δ_D is still normally not computed.

Here we will take advantage of the fact that existing estimates of ν at the DQCP (ν ≈ 0.45 in both the J–Q^[26] and loop^[32] models) correspond to a rather small value of the scaling dimension, Δ_D ≈ 0.8, and therefore it may be possible to compute it reliably in this case (as was done recently for the transverse-field Ising chain, where, in the notation used here, Δ_D = 1^[44]). Furthermore, we point out that off-diagonal correlation functions of operators that are terms of the Hamiltonian have very simple estimators within the Stochastic Series Expansion (SSE) QMC method.^{[29,30,45,46]} The quantum fluctuations are here represented by a string of length n of terms H_i of H, with mean length 〈n〉 = |〈H〉|/T, where T is the temperature. A connected correlation function of any two terms is given by^[46]

$$\begin{eqnarray}\begin{array}{ll}{C}_{ab} & \equiv \langle {H}_{a}{H}_{b}\rangle -\langle {H}_{a}\rangle \langle {H}_{b}\rangle \\ & ={T}^{2}(\langle {n}_{ab}(n-1)\rangle -\langle {n}_{a}\rangle \langle {n}_{b}\rangle ),\end{array}\end{eqnarray} \tag{ 2 }$$

where n_a is the number of operators H_a in the string and n_ab is the number of times that H_a and H_b appear adjacent to each other. This expression can be easily applied to all location pairs (a,b) in a single scan of the operator string, and translational invariance can be exploited at no additional cost to improve the statistics.

As a demonstration of the method, we first consider the S = 1/2 bilayer Heisenberg Hamiltonian

$$\begin{eqnarray}&&H={J}_{1}\displaystyle \sum _{a=1,2}\displaystyle \sum _{\langle ij\rangle }{S}_{a,i}\cdot {S}_{a,j}+{J}_{2}\displaystyle \sum _{i=1}^{N}{S}_{1,i}\cdot {S}_{2,i},\end{eqnarray} \tag{ 3 }$$

where 〈ij〉 denotes nearest neighbors on a square periodic lattice with N = L² sites and a is the layer index. This system has an AFM ground state for g ≡ J₂/J₁ < g_c and is a quantum paramagnet dominated by inter-layer singlet formation for g > g_c. The O(3) quantum phase transition has been investigated in many previous works. Here we take g_c = 2.52205 for the critical point^[47,48] and study a correlation function corresponding to the perturbation D in Eq. (1). Since both the J₁ and J₂ interactions drive the system away from the critical point, we can study correlations between either type of terms (i.e., they have the same scaling dimension). We use the J₂ terms, which form a simple square lattice, and define

$$\begin{eqnarray}&&{C}_{2}({r}_{ij})\equiv \langle ({S}_{1,i}\cdot {S}_{2,i})({S}_{1,j}\cdot {S}_{2,j})\rangle -{\langle {S}_{1,i}\cdot {S}_{2,i}\rangle }^{2},\end{eqnarray} \tag{ 4 }$$

where r_ij denotes the separation of the sites i and j.

Investigating the decay of the correlations, we can either study large lattices and focus on r ≪ L to eliminate finite-size effects or take r of order L and study the size dependence. Here we opt for the latter method with r = (L/2 – 1,0), for which there are more equivalent points for averaging than for the high-symmetry points (L/2,0) and (L/2,L/2). For the expected O(3) universality class in 2+1 dimensions ν ≈ 0.71,^[49] corresponding to a scaling dimension Δ₂ ≈ 1.594 of the J₂ interaction. As shown in because of the rapid decay we can access only rather modest distances, but the results still show a remarkably good agreement with the expected form C₂(r) ∝ r^−2Δ₂ starting from r = 4 (L = 10). In the SSE simulations we have used T = c/L (in units with J₁ = 1), reflecting the emergent Lorentz invariance of the system (i.e., the dynamic exponent z = 1), with two different proportionality factors; c = 2 and 1/8. Apart from the different amplitudes of the correlations, both data sets exhibit the same decay.

Turning now to the J–Q model, we express the AFM Heisenberg interaction as a singlet projector, –P_ij, on S = 1/2 spins; P_ij = 1/4 – S_i ⋅ S_j. To simplify the notation, we use a bond index b to implicitly refer to two nearest-neighbor spins 〈i,j〉_b; P_b ≡ P_ij. We also use an index p to refer to a 2 × 2 plaquette with sites in the arrangement ${\left(\begin{array}{cc}i & j\\ k & l\end{array}\right)}_{p}$ and define Q_p ≡ P_ij P_kl + P_ik P_jl. With these definitions the J–Q Hamiltonian is^[14]

$$\begin{eqnarray}&&H=-J\displaystyle \sum _{b}{P}_{b}-Q\displaystyle \sum _{p}{Q}_{p}.\end{eqnarray} \tag{ 5 }$$

We define the coupling ratio g ≡ J/Q and use the SSE method to compute the z component of the staggered magnetization (the AFM order parameter)

$$\begin{eqnarray}&&{m}_{z}=\displaystyle \frac{1}{N}\displaystyle \sum _{r}{S}_{r}^{z}{(-1)}^{{r}_{x}+{r}_{y}},\end{eqnarray} \tag{ 6 }$$

and the two-component dimer (VBS) order parameter, also defined with the z spin components,

$$\begin{eqnarray}&&{d}_{\alpha }=\displaystyle \frac{1}{N}\displaystyle \sum _{r}{S}_{r}^{z}{S}_{r+\hat{\alpha }}^{z}{(-1)}^{{r}_{\alpha }},\end{eqnarray} \tag{ 7 }$$

where α stands for the x or y lattice direction. We scale the temperature in units of Q as T = c/L, with c = 2.38 being the estimated critical velocity of excitations^[25] (i.e., the system is in the "cubic" scaling regime,^[48,50] as in the case 1/T = L/2 for the bilayer model in Fig. 1).

Zoom In Zoom Out Reset image size

Early QMC studies placed the VBS–AFM transition at g_c ≈ 0.040,^[14–16] while more recent works show a somewhat larger value, g_c ≈ 0.045,^{[18,25,26,30]} as a consequence of significant finite-size corrections. We now have data for system sizes up to L = 512 and present the Binder cumulants U_z and U_d defined in the standard way such that U_x → 1 with increasing system sizes if there is order of type x and U_x → 0 otherwise;

$$\begin{eqnarray}&&{U}_{{\rm{z}}}=\displaystyle \frac{5}{2}-\displaystyle \frac{5}{6}\displaystyle \frac{\langle {m}_{z}^{4}\rangle }{{\langle {m}_{z}^{2}\rangle }^{2}},\,\,\,{U}_{{\rm{d}}}=2-\displaystyle \frac{\langle {({d}_{x}^{2}+{d}_{y}^{2})}^{2}\rangle }{{\langle {d}_{x}^{2}+{d}_{y}^{2}\rangle }^{2}}.\end{eqnarray} \tag{ 8 }$$

Results for several system sizes are shown in Fig. 2(a).

Zoom In Zoom Out Reset image size

To improve the g_c estimate, we analyze crossing points g = g*, where U_z(g*,L)=U_d(g*,L) and also where (for different g*) U_x(g*,L/2) = U_x(g*,L) with x = z or x = d. As shown in Fig. 2(b), these crossing points flow to g_c = 0.04510(2) as L → ∞. The extrapolation is based on a fit to two power laws for each data set, with a common g_c. Unconstrained fits also result in consistent g_c values. We have excluded small systems until a statistically sound fit is obtained, with L ≥ 64 included in the final analysis. From now on we fix the coupling ratio to g = 0.0451 ≈ g_c.

We here examine the correlation function of the Q-terms in the Hamiltonian (5),

$$\begin{eqnarray}&&{C}_{Q}({r}_{ij})=\langle {Q}_{i}{Q}_{j}\rangle -{\langle {Q}_{i}\rangle }^{2},\end{eqnarray} \tag{ 9 }$$

which is less noisy than the J-energy correlator. As shown in Fig. 3(a), the correlations exhibit strong even-odd oscillations, with amplitude decaying with the distance. The reason for the oscillating behavior is that the columnar VBS correlations are also detected by the plaquette correlation function C_Q(r) (for a detailed general discussion of this, see Ref. [35]). In a columnar state with x-oriented dimers, C_Q(0,y) will be small while C_Q(x,0) will have signs (–1)^x due to the dimerization. In an ergodic QMC simulation, C_Q(x,y) will reflect averaging over states with x- and y-oriented dimers. The contributions from the VBS order parameter then cancel in C_Q(x,0) for odd x, while C_Q(x,x) retains the VBS contributions with (–1)^x signs. These behaviors are seen in Fig. 3(a), where the amplitude decay is due to the system being a critical VBS. Since the system has emergent U(1) symmetry of the order parameter,^[14,16] we should consider C_Q(r) as averaged over an angle ϕ ∈ [0,2π) corresponding to a circular-symmetric distribution P(d_x,d_y). The above-mentioned behaviors of C_Q(r) along the lines r = (x,0) and r = (x,x) will also hold in this case.

Zoom In Zoom Out Reset image size

In addition to the large contributions to C_Q(r) from the VBS order parameter, there should be a uniform component reflecting the scaling dimension of the full Q operator. Since the VBS contributions are absent at (x,0) with odd x, examining the correlations at these distances is a good way to access the uniform component. In Fig. 3(a), small rapidly decaying values are indeed seen, and in Fig. 3(b) the functional form is analyzed on a log–log plot. We use a large system, L = 256, with x ≪ L, as well as x = L/2 – 1 for smaller sizes. In both cases we observe the same algebraic asymptotic decay, and a power-law fit to the x = L/2 – 1 data for x > 12 gives Δ_Q = 0.800(4). This scaling dimension corresponds to 1/ν = 2.200(4), in good agreement with the previous (less precise) results for the J–Q^[26] and loop^[32] models.

Next we consider the cumulant slopes S_x ≡ dU_x/dQ, x = d,z, computed with direct SSE estimators as previously carried out for S_z with L ≤ 160 in Ref.[26]. Here we present the results for L up to L = 448 (our L = 512 results are too noisy). The slopes should scale asymptotically as L^1/ν. In order to account for the leading correction we also include a second power-law term with smaller exponent, and exclude small systems until good fits are obtained. The results are shown in Fig. 4. The inset shows the same data sets and fits converted into 1/ν* ≡ ln[S(L)/S(L/2)]ln⁻¹(2), which flows to 1/ν as L → ∞. We note that: (i) 1/ν = 2.23(2) is fully consistent with the previous result from smaller systems,^[26] and (ii) the value also agrees with the above result from the scaling dimension of the Q terms (with a difference less than 1.5 standard deviations).

Zoom In Zoom Out Reset image size

While the finite-size corrections in 1/ν obtained from the cumulant slopes in Fig. 4 are substantial, the corrections to the r^−2Δ_Q form of the correlation function in Fig. 3 are very small. The good agreement of the extracted exponents with the relationship 1/ν = 3 – Δ_Q should alleviate any concerns of 1/ν eventually flowing to the value 3 ( = d) expected at a conventional first-order transition (or to d + 1, as found at an unconventional transition in Ref.[33]). Weak first-order transitions are often most clearly manifested in 1/ν,^[51] and the results presented here simply do not indicate anything unusual.

Previously, anomalous scaling was found of the order parameters and the spin stiffness,^{[16,18,22,32]} and it was argued that the standard finite-size scaling hypothesis must be replaced by a form taking into account two divergent length scales in a new way.^[26] Though this interpretation has not been independently confirmed, the results presented here reinforce the notion that anomalies are not present in the magnitude L^−1/ν of the critical window. The value ν=0.455(2) is still puzzling in the sense that it violates a CFT bound from the bootstrap method.^[52] This disagreement suggests that the transition is either not described by a CFT or that the arguments underlying the bound have some loophole. It would be interesting to compute ν from the relevant critical correlator for fermionic DQCP candidate models,^[53] where the standard finite-size analysis is difficult because of the rather small accessible lattices.

We would like to thank Jun Takahashi for useful discussions. The numerical calculations were carried out on the Shared Computing Cluster managed by Boston University's Research Computing Services.