{'fraction_non_alphanumeric': 0.035531972058033316, 'fraction_numerical': 0.012667920472864051, 'mean_word_length': 4.395448285859245, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Correlated valence-bond states

We study generalizations of the singlet-sector amplitude-product (AP) states in the valence-bond basis of S=1/2 quantum spin systems. In the standard AP states, the weight of a tiling of the system into valence bonds (singlets of two spins) is a product of amplitudes depending on the length of the bonds. We here introduce correlated AP (CAP) states, in which the amplitude product is further multiplied by factors depending on two bonds connected to a pair of sites (here nearest neighbors). While the standard AP states can describe a phase transition between an antiferromagnetic (Neel) state and a valence-bond solid (VBS) in one dimension (which we also study here), in two dimensions it cannot describe VBS order. With the CAP states, Neel-VBS transitions are realized as a function of some parameter describing the bond correlations. We here study such phase transitions of CAP wave-functions on the square lattice. We find examples of direct first-order Neel-VBS transitions, as well as cases where there is an extended U(1) spin liquid phase intervening between the Neel and VBS states. In the latter case the transitions are continuous and we extract critical exponents and address the issue of a possible emergent U(1) symmetry in the near-critical VBS. We also consider variationally optimized CAP states for the standard Heisenberg model in one and two dimensions and the J-Q model in two dimensions, with the latter including four-spin interactions (Q) in addition to the Heisenberg exchange (J) and harboring VBS order for large Q/J. The optimized CAP states lead to significantly lower variational energies than the simple AP states for these models.


I. INTRODUCTION
The valence-bond (VB) basis 1-6 is ideally suited for describing many different types of ground states and low-energy excitations of quantum spin models. [7][8][9][10][11][12][13][14][15][16][17] In the case of S = 1/2 spins in the singlet sector, a basis state corresponds to a tiling of the lattice into bonds connecting pairs of sites forming singlets, such that each spin belongs to one bond. This basis is overcomplete if bonds of all lengths are included. To describe the ground state of a Hamiltonian with bipartite interactions, only bonds connecting sites on different sublattices have to be included-this restricted VB basis exactly reproduces Marshall's sign rule 18 for the ground state of such a system. Thus, in this basis the wave function is positive definite and can be sampled using Monte Carlo (MC) techniques. We here investigate a class of bipartite correlated VB wave functions which can exhibit valence-bondsolid (VBS) order and related interesting quantum phase transitions in one and two dimensions.
In this introductory section we provide some further background and motivation for studying VB states. We review the definition and properties of the well-studied Liang-Doucot-Anderson amplitude-product states 4 and introduce their more versatile generalizations-the correlated AP (CAP) states that we focus on in this paper. We discuss reasons to study such states in the context of quantum phase transitions from the antiferromagnetically ordered Néel state into non-magnetic VBS and spin liquid states.

A. Valence-bond states and Marshall's sign rule
While some analytical work has been carried out in the VB basis, 8,12 in most quantitative calculations MC sampling of the bonds must normally be used to reliably evaluate expectation values. Since the basis is overcomplete, the non-negative definiteness of the wave function is a requirement to avoid problems due to negative sampling weights (the sign problem). Thus, in most cases VB MC calculations are restricted to bipartite (non-frustrated) systems. A two-spin singlet (VB) connecting sites a and b on sublattices A and B is then defined according to the following phase convention: Marshall's sign rule is then incorporated for any tiling of an even number N of spins into N/2 singlets, i.e., when expressed in the standard basis of z spin components, |Z = |S z 1 , . . . , S z N , the sign of a non-zero coefficient ψ V (Z), i.e., for states with antiparallel spins on each bond, is given by where n A↓ is the number of ↓ spins on sublattice A. The wave function ψ 0 (V ) of the ground state of such a system expressed in the VB basis, is therefore non-negative. When using MC simulations, e.g., with a variational wave function |Ψ approximating |Ψ 0 , this is essential, because the overcompleteness implies that the sampling weigh is not |ψ(V )| 2 , as it would be in a standard orthogonal basis, but ψ(V )ψ(V ′ ) V ′ |V for simultaneously sampled non-orthogonal bra and ket configurations |V and V ′ |. The overlap V ′ |V and matrix elements of operators of interest in characterizing the states can be easily calculated, 3-5 as we will discuss later below.

B. Amplitude-product states
The most commonly used variational states in this context are the AP states introduced by Liang, Doucot, and Anderson. 4 Here one associates a bond connecting two sites (a, b) with an amplitude h(a, b), which in the case of a translationally-invariant system is only a function of the lattice vector r ab separating the two sites; h(a, b) = h(r ab ). The wave function coefficient for a VB configuration V is then where n r is the number of bonds of "shape" r in the configuration.
The amplitudes h(r) can be used as variational parameters. In the original work with the AP states to describe the ground state of the two-dimensional (2D) Heisenberg model, 4 only the amplitudes for a small number of short bonds were optimized, and different functional forms (exponentially or power-law decaying with the distance r) were tested. In later work all the amplitudes (on finite lattices) were optimized, leading to relative energy errors (deviations from results from unbiased quantum Monte Carlo, QMC, calculations) of less than 0.1%. 19,20 In the optimal state the amplitudes decay asymptotically as 1/r 3 , which is also the result of a mean-field VB approach. 12 In some cases, if one is just interested in the properties of some class of states without reference to a specific Hamiltonian, the optimization step is not needed. This approach has been taken in recent studies of the prototypical resonating VB (RVB) spin-liquid state consisting of the superposition of all configurations of the shortest (nearest-neighbor) bonds on the square lattice, [21][22][23] and also in the presence of some fraction of the second bipartite bond (fourth-neighbor). 21 These wave functions, for which the parent hamiltonian was recently identified (in the case of nearest-neighbor bonds only), 24 has exponentially decaying spin correlations but power-law decaying c 0 c 1 c 2 c 3 c 4 FIG. 1: (Color online) Configurations of short valence bonds (shown in red) connected to a lattice link b (indicated by thick black bars). Their associated CAP weights, Eq. (7), are C b (r1, r2) with r1 = r2 = 1 and b here being a horizontal link. We will later use a notation with weights ci for CAP states including only short-bond correlations, with i corresponding to (r1, r2) according to the labeling above.
VBS correlations. A phase transition from the Néel state into this kind of spin liquid can be achieved by using amplitudes of the form h(r) ∝ 1/r κ and tuning the exponent κ to a critical value. 12,25 C. AP states with bond correlations (CAPs) One of the motivations of the work reported in the present paper is to obtain a variational description of the 2D Néel-VBS transition. For this, we need a class of wave functions beyond the AP states, since they do not exhibit VBS order (while the 1D variants do, as we will discuss in Sec. III). The 2D non-magnetic AP states are beleived to always be spin liquids, with exponentially decaying spin correlations and power-law VBS correlations, similar the prototypical short-bond RVB states. 12, 21 We study a class of generalized AP states defined with bond-correlation factors multiplying the AP wave function (6). We take these factors to be of the form C b [r 1 (b), r 2 (b)], where b denotes a nearest-neighbor link on a 1D chain or 2D square lattice (or, more generally, any lattice with some imposed bipartition), and r 1 (b), r 2 (b) are the shapes of the two VBs connected to this bond (with the case of there being just a single bond connecting the two sites being a special case). Thus, the wave-function coefficient is For a translationally invariant system C b (r 1 , r 2 ) for given (r 1 , r 2 ) depends only on the orientation (horizontal or vertical in two dimensions) of the bond b, and these weights also should obey applicable lattice symmetries. The number of different correlation factors is then ∝ N 2 for a system of N spins. For simplicity of the notation we hereafter suppress the subscript b.
In principle, in variational CAP calculations all correlation factors can be optimized, along with the amplitudes h(r), but one can also opt to consider only those factors C(r 1 , r 2 ) for which r 1 , r 2 ≤ r max , with some maximum bond length r max , and set the remaining weights to unity. The possible two-bond configurations with r max = 1 are illustrated in Fig. 1. In variational cal-culations one would expect the energy to decrease monotonically with r max , which we will demonstrate explicitly in Sec. V.
Beyond improving the energy in variational calculations, the correlation factors also play an important qualitative role in 2D systems-without bond correlations, the standard AP states are either long-range Néel ordered (although they do not, by construction, break the spin-rotation symmetry, they can still develop the magnitude of the sublattice magnetization) or are RVB spin liquids with critical VBS correlations (as discussed above in the context of the short-bond RVBs). They cannot form VBS order. In contrast, the trivial 1D AP state with only short bonds is an extreme case of a two-fold degenerate VBS state with alternating links with or without a VB. This kind of long-range order remains stable also in the presence of some fraction of longer bonds, as we will discuss below in Sec. III. The generalized CAP states (7) can exhibit 2D VBS order if the correlation factors favor such correlations strongly enough. This is true even with correlations only involving only the shortest bonds on the square lattice, illustrated in Fig. 1. The CAP states open the possibility to study the Néel-VBS transition in classes of wave functions with the bond correlations parametrized in some way, and to carry out improved variational calculations for systems with VBS order or Néel states with significant VBS fluctuations.

D. Purpose and outline of the paper
One of our main reasons to study the CAP states here is to investigate their abilities to describe the 2D Néel-VBS transition. This transition has received considerable interest recently in the context of "deconfined" quantum-criticality (DQC). 26,27 Following earlier work on VBS states and quantum-critical phenomena in antiferromagnets 10,28-31 and topological aspects of phase transitions in 3D classical spin systems, 32 Senthil et al. proposed 26,27 that the 2D Néel-VBS transition is of an unusual kind where the standard Landau rule stipulating a generically first-order transition between the two ordered states is violated. The two order parameters in the DQC scenario are both consequences of the same underlying more fundamental objects-spinons interacting with an emergent U(1) symmetric gauge field. The spinons condense in the Néel state and confine as pairs in the VBS state. Unbiased QMC studies of J-Q models, 33-38 which include certain multi-spin interactions (Q) in addition to the standard Heisenberg exchange (J), are in general in good agreement with the theoretical predictions. Among the most interesting features observed is an emergent U(1) symmetry of the VBS order parameter (presumably reflecting the emergent gauge field-the "photon") as the critical point is approached from the VBS side. Moreover, studies of SU(N ) generalizations of the J-Q model 36 and other spin models 39 have allowed direct connections with analytical large-N calculations for the non-compact CP N −1 field theory argued 26,27 to describe the transition.

Scope of the paper
We here investigate whether the 2D Néel-VBS transition can be correctly captured with a simple ansatz wave-function of the form (7) with fixed single-bond amplitudes (of a power-law form) and continuously varying short-bond correlation weights of the form in Fig. 1. The result so far is negative, in the sense that we do not observe the same kind of continuous VBS transition as in the J-Q model. Instead, with parameters chosen such that there is direct Néel-VBS transition, we find strong discontinuities. In other cases we find an RVB spin liquid intervening between the ordered phases. Thus, it still remains an interesting challenge to find a simple VB description of the DQC point.
Looking at the VBS order-parameter distribution, we do not observe any emergent U(1) symmetry of the VBS at the continuous VBS to RVB transition, impliying that this is not a "deconfined" transition in this case. Nevertheless, we find interesting scaling properties of the angular VBS fluctuations, although the length-scale associated with them are not divergent.
In addition to the 2D studies, we also closely examine the Néel-VBS transition within the standard AP states in one dimension, with amplitudes h(r) of the form 1/r α . This transition, which occurs at a critical value of α (which is not universal but depends on the detailed form of the amplitudes for small r) was previously studied by Beach, 12 but only the spin correlations were computed. Here we extract also the VBS correlations and confirm that there is a single critical point versus α. The exponents are continuously varying, depending on the shortbond amplitudes.
We also report variational calculations including optimization with the CAP states, minimizing the energy for 1D and 2D Heisenberg and J-Q models. Naturally, bond correlations have a significant improving effect in VBS phases, but they help also to improve the 2D Néel state and the critical ground state of the 1D Heisenberg chain.

Outline of the paper
In Sec. II we briefly describe the technical aspects of MC calculations with AP and CAP states. In Sec. III we study the Néel-VBS transition in 1D AP states, and in Sec. IV we study the more rich set of states and quantum phase transitions in 2D CAP states. In Sec. V we present some 1D and 2D AP and CAP variational calculations (minimizing the energy as a function of the amplitudes and correlation factors) for prototypical model Hamiltonians with Néel and VBS order (the Heisenberg and J-Q models), showing how bond correlations improve the states. We conclude in Sec. VI with a summary and discussion of future prospects.

II. MC SAMPLING OF CAP STATES AND CALCULATION OF OBSERVABLES
In an AP or CAP state with the wave function of the form (7) the expectation value of an observableÔ can be written for the purpose of importance sampling as where the configuration weight is and the estimator corresponding toÔ given by Here the overlap V β |V α is evaluated by counting the number l αβ of loops in the transition graph of V α and V β ; 3,4 An example with two loops is shown in Fig. 2. Below we briefly discuss MC sampling of the VB configurations and estimators for some important observables.

A. Bond and spin sampling schemes
For a given functional form for h(r) and bond correlation factors C(r 1 , r 2 ), an expectation value Ô can be computed stochastically by importance sampling according to the weight W αβ . Using some random reconfiguration of bonds in either the state V α or V β or both of them, the standard Metropolis acceptance probability for a modified configuration (V α ′ , V β ′ ) is where the weight ratio is Even for the simplest update involving changes in just two bonds, 4 the calculation of the change in the new number of loops (l α ′ β ′ −l αβ ) can require a computational time up to ∝ N for each new update proposal, since a Néel state has extensive loops (while magnetically disordered states have only short loops). Since the number of such updates in each MC sweeps should also be proportional to N , this type of update leads to a total computational time O(N 2 ) for a full sweep in a Néel state, while in a non-magnetic state the scaling is O(N ). The unfavorable scaling in the Néel state can be avoided by working in a combined space of both spins and bonds, 19 where the VBs are also sampled, by randomly selecting either ↑ a ↓ b or ↓ a ↑ b for each singlet (a, b). Since the spin basis is orthogonal, all spins in the bra and the ket have to be the same, and a consistent assignment for both V α and V β , thus, implies that the spins on each loop in the transition graph follow a staggered, ↑↓↑↓ . . ., pattern. The overlap (11) in the pure VB basis then follows, since there are two possible staggered configurations on each loop. The spins are periodically updated by flipping all spins in randomly selected loops. By such spin sampling, the weighting by the number of loops is accounted for automatically, due to the entropic effect of favoring configurations with large numbers of loops, without the need for actually counting the loops. For a detailed description of the combined spin-bond basis and simulations in it, we refer to Refs. 19 and 21. Here we just note two different ways of updating the bond configurations: (i) In the two-bond reconfiguration scheme, an elementary MC move consists of choosing two sites on the same sublattice at random (typically the two bonds on a randomly chosen pair of next-nearest-neighbor spins), and exchanging the bonds connected to these two sites for the other possible bipartite configuration. Such a reconfiguration is only possible if the spin states on the two selected sites are the same. If that is the case, the acceptance probability (12) is applied, where in the combined spin-bond basis the weight is instead of Eq. (9), and with W α ′ β ′ /W αβ evaluated using only the bond amplitudes and correlation factors in (7) affected by the change.
(ii) In a loop update, we start by removing a dimer randomly from two connected sites, creating two defects ("holes"). We keep one defect stationary and move the second one by connecting one end (the one on the same sublattice) of a chosen bond to it (hence moving the hole to the previous location of that end of the bond). The bond to move should be chosen probabilistically in such a way as to satisfy detailed balance, which is relatively straight-forward in the case of AP states 19,40 but more complicated when bond correlations are included. In the present work we have used loop updates only for pure AP states, while we use two-bond updates for CAP states. The latter are also efficient enough to study relatively large lattices (with thousands of spins).
It should be noted here that VB configurations can be classified according to topological "winding numbers". 43 In AP or CAP states defined with only short bonds, the two-bond update conserves the winding number, but with no restriction on the bond-length such updates can change the winding number. In practice, if the bond probability (which depends on the single-bond amplitudes as well as the correlation factors in CAP states) decays very rapidly with the length, a simulation for a large system may still be confined to the sector of zero winding number.

B. Spin and dimer correlations
In order to characterize the different phases realized by the CAP states, we evaluate order parameters for detecting antiferromagnetic (Néel) order and VBS order. Néel order can be characterized using the standard two-spin correlation function, where we use r ij to denote the vector separating the lattice sites i and j and the phase factor cancels the signs of the staggered spin correlations obtaining in the systems we study. Alternatively, one can study the full sublattice magnetization averaged over the whole system; where φ i = +1 on sublattice A and φ i = −1 on sublattice B. Since the singlet AP and CAP states manifestly cannot break the spin-rotation symmetry, order must be detected in the squared order parameter, m 2 s , which in the limit of large system size will be identical to the long-distance spin correlation (15).
To accurately locate an antiferromagnetic phase transition, the Binder cumulant is very useful. It is defined according to 42 where the factors are chosen for the 3-component Néel order parameter such that U (N → ∞) = 0 in the disordered phase (where the order-parameter distribution is a Gaussian with zero average) and U (N → ∞) = 1 in the ordered phase (where the radial distribution is peaked at non-zero m s ). Typically crossing points of U graphed versus a control parameter for different system sizes approach the critical point vary rapidly as a function of increasing system size.
To characterize VBS order we use the dimer correlation function, defined as in terms of the bond operators directed along the unit lattice vectorsx andŷ. We will not need the mixed x-y correlations here. In some cases we will characterize VBS order by the long-distance behavior of (18). The states we will be studying have a 2-site VBS unit cell, forming a staggered weak-strongweak-strong pattern in one dimension and an analogous columnar pattern in two dimensions. In both cases we can extract the dominant component of the correlations, corresponding to the squared order parameter, by taking the appropriate difference of (18) evaluated at nearby distances. We here use a symmetric version of this difference; and a function D * yy (r) for y-oriented dimers defined analogously. We will also study the full order parameter, which in two dimensions can be defined using the q = (π, 0) and q = (0, π) Fourier transforms of the nearest-neighbor bond correlations (19); The magnitude D of the order parameter can be computed as the square-root of the average squared operator, D 2 = D 2 x + D 2 y . In addition to the expectation values, we will also investigate the probability distribution P (d x , d y ), in which emergent U(1) symmetry can be detected. Here d x and d y are the expectation values of the corresponding operators (21) evaluated in a given sampled configuration based on the transition graph. We refer to Ref. 52 for further details on this quantity, which is not a conventional quantum mechanical expectation value but still very useful for characterising VBS states.
All of the above two-and four-spin correlations are related to the transition-graph loops generated in the VB MC sampling process. For instance, the estimator for the two-spin correlation is given by 3,4 where L ℓ is the size (the number of sites) of loop ℓ.
Both the dimer correlation function and the fourth power of the staggered magnetization involve four-spin correlations. Detailed descriptions on how to calculate these based on the transition graph of two VB configurations can be found in Refs. 5 and 21. Here we only write down the expression for the fourth power of the staggered magnetization, needed for the Binder cumulant (17), which is also solely determined by the sizes of all loops formed in the transition graph. We note that the Binder cumulant of the VBS order parameter is much more difficult to evaluate, since its definition in analogy with (17) requires eight-spin correlations. While these also in principle can be evaluated in terms of the transition-graph loops, 5 the expressions are quite complicated to implement in practice and we have not done so.

III. NÉEL TO VBS TRANSITION IN ONE DIMENSION
In one dimension, the standard AP states given in Eq. (6) are able to reproduce a Néel-VBS transition without correlation factors. We will study this 1D transition carefully in this section, using the very efficient loop update of the VB configurations.
It is natural to study the evolution of the state as a function of some parameter governing the long-distance behavior of the amplitudes, e.g., using the power law h(r) = 1/r κ with tunable κ or an exponential form. Here we will use the power-law. However, it is already known that the nature of the state is not just determined by the asymptotic behavior of h(r), but also depends on details of the short-bond weights. 12 In addition to the exponent κ we here tune the shortest-bond amplitude h(r = 1) = λ. The wave function is, thus, explicitly given by where n r (V ) again refers to the number of bonds of length r in the bond configuration V . It is clear that for λ > 0 and large κ this AP state is a VBS, since in the limit κ → ∞ only two configurations contribute; those with r = 1 bonds on alternating links. For small κ there is instead Néel order but no VBS order. 12 Note that long-range order corresponding to broken SU(2) symmetry is possible in this kind of 1D system, since viewed as a classical statistical-mechanics problem there are long-range interactions (since the bonds have unbounded length), and the Mermin-Wagner theorem 41 prohibiting 1D Néel order does not apply. Note also again that the AP wave function is a singlet and, thus, the SU(2) symmetry is not actually broken (as in any calculation targeting the singlet ground state). The magnitude of the Néel order measured by m 2 s , Eq. (16), or the longdistance correlation function (15) can still evolve toward a non-zero value as the system size grows, tending to the square of the symmetry-broken value of m s in the corresponding thermodynamic-limit state with no constraint on the total spin. Beach has previously studied Néel ordering in this class of wave functions (with a somewhat different parametrization of the short-bond amplitudes). 12 He found a continuous transition between the Néel state and the non-magnetic state. Here we also investigate the VBS correlations and find a single transition point where both the spin and dimer correlations are critical. We study the evolution of the transition in the plane (κ, λ).
For fixed λ, in order to find the critical value of κ c of the AP state we study the Néel Binder cumulant (17). The behavior of curves for different system sizes L crossing each other as a function of κ is illustrated in the upper panel of Fig. 3. The crossing points do not fall exactly on a single point due to subleading size corrections. We observe a systematic smooth drift of the crossing points as the system size is increased. In order to eliminate this size effect and determine the critical point from data such as those in Fig. 3, we extract κ-values corresponding to crossing points of (L, 2L) size pairs, and plot these points against 1/L, as shown in the lower panel of Fig. 3. We then extrapolate these values to L → ∞ and obtain κ c . The fitting function we use here for extrapolation is the standard power-law; 42 The extrapolated κ c values versus λ are plotted in Fig. 4; the phase diagram of 1D AP states with the two tuning parameters λ and κ. The inset of this figure demonstrates the qualitatively different behaviors of the spin correlation functions in the two phases and at the critical point, using λ = 1 results as an example. At κ = 1.6 the correlations decay faster than power-law, as is expected for a non-magnetic VBS ordered state. In contrast, at κ = 1.4, the correlations for small L first decay somewhat but then converge to a non-zero value for larger L (even increasing somewhat for large systems), demonstrating the presence of long-range Néel order. At the critical value κ c extracted using the Binder crossings as explained above, the decay of the correlations are consistent with a critical, power-law form.
To determine whether the VBS correlations are also critical at the κ c points extracted from the Néel Binder cumulant, we further study both the spin and dimer correlations at these points. The results confirm the expectation of a common critical Néel and VBS point. By studying chains as large as L = 4096, we can extract the exponents governing the critical correlation functions with relatively small error bars (thanks to the powerful VB MC loop update discussed in Sec. II). The analysis of the power laws is presented in Fig. 5. Note that in order to avoid boundary modifications of the power-law correlation functions as a function of the distance r in systems of fixed L, we study the long-distance correlations versus the system size, with r = L/2 for the spin correlations and the staggered component of the dimer correlations extracted based on r = L/2 and L/2 − 1 As λ increases, larger system sizes are needed to observe the asymptotic critical forms. Especially for the largest λ studied, λ = 8, one can observe in Fig. 5 (upper panel) a clear cross-over of the spin correlation function from a rapidly decaying short-distance form to the asymptotic power-law form. The straight lines in Fig. 5 are fits to the simple asymptotic form aL −b . We have also tried to include shorter chains in an analysis including corrections, by fitting to the form aL b + cL d . This form is, however, not capable of describing the small size effect in this model (in contrast to 2D critical spin liquid RVB states, where this form works very well 21 ). In any case, the large-L behaviors appear to be reasonably well converged to the simple power law and the exponents extracted should be reliable. An exception is λ = 0, for which the dimer correlations decay very rapidly and are too noisy to allow the exponent β to be reliably determined (and we have therefore not graphed these correlations in Fig. 5). It is even possible that the VBS state for λ = 0 is of a different kind than for λ > 0. Further studies will be needed to settle this issue.
We plot the extracted critical exponents as a function of λ in Fig. 6. The exponents vary continuously with λ, with the dimer exponent decreasing monotonically and the spin exponent increasing. An interesting conclusion that can be drawn from these results is that the critical state becomes increasingly "quasi VBS ordered" with increasing λ, with the decay exponent of the dimer correlations perhaps vanishing as λ → ∞, although this is difficult to confirm definitely (because the simulations become increasingly difficult for large λ). The behavior is in line with the expectation that a large lambda favors VBS ordering because of the predominance of the very shortest bonds, i.e., when moving on the critical line toward higher λ the density of short bonds increases, and this leads to a strengthening of the VBS quasi-order. At the same time, the exponent of the spin correlations appear to approach 1. However, Néel order still exists for large λ when reducing κ from the critical value. In terms of the transition graph estimators of the correlation functions, VBS correlations correspond to certain loop correlations, 5 while Néel order is related to the presence of long (∝ L) loops. While long-range Néel and VBS orders are mutually exclusive in these states, the Néel state in the neighborhood of the critical curve for large λ approaches a coexistence situation. Here the magnitude of the Néel order parameter also becomes very small, however, and the coexistence is therefore not robust.
The VB formulation of the ground state can be viewed as a 1D classical statistical mechanics problem, but at the same time it should also correspond to a path-integral formulation in 1 + 1 dimensions (with some underlying parent Hamiltonian). One may then expect the system to be classifiable according to the standard 2D conformal field theories by a central charge c. Varying critical exponents, as we have found here, normally imply c ≥ 1, but the fact that the system includes long-range interactions may invalidate this requirement, although it is not clear how the power-law bond length translates into effective interactions in an underlying parent Hamiltonian (and the interactions in it may well be short-ranged). One notable aspect of the AP states is that they are not able to reproduce the ground state of the critical Heisenberg chain, where α = β = 1. 44 We will address this issue further in Sec. V with variational AP calculations for the Heisenberg chain.
It would be interesting in the future to compute the bipartite entanglement entropy of the 1D AP states, to test its system size scaling and consistency between c extracted from it 45,46 and from the correlation functions. Such calculations can also be carried out using the VB MC sampling scheme used here. 23,47

IV. NÉEL TO VBS TRANSITION IN TWO DIMENSIONS
The Néel state is known to be the ground state of the square lattice Heisenberg antiferromagnet with homogeneous nearest-neighbor couplings. There is convincing numerical evidence 20 as well as mean-field arguments 12 showing that the standard AP states with power-law decaying amplitudes of the asymptotic form h(r) = 1/r 3 is an optimal variational wave function for the 2D Néel state (for any finite size, where the ground state is a singlet with no explicitly broke symmetry).
The AP states have no Néel order for rapidly decaying (exponentially or according to a power law 1/r κ with large κ) bonds. 4 An extreme case is the state that contains only nearest neighbor bonds (dimers). Such a shortbond VB state on the square lattice normally corresponds to an U(1) spin liquid with critical VBS correlations, 21 in contrast to the 1D AP VBS state discussed in the previous section. One should expect the 2D AP state to turn into a long-range ordered VBS when appropriate bond correlations are included. On the square lattice, which we will consider here, the simplest kinds of VBS states (with 2-site unit cell) form columnar, staggered, or plaquette ordering patterns. We here study the Néel to VBS transition on a square lattice within the CAP states by imposing bond correlations that favor or suppress only certain types of shortdimer alignments. All possible configurations of short dimers connected to a pair of nearest-neighbor sites on a square lattice are shown in Fig. 1. We assign a weight to each of those two-dimer configurations according to Eq. (7), with all weights with r 1 , r 2 = 1 set to 1. To simplify the notation we here use c i for the short-bond correlation factors, instead of C(r 1 , r 2 ), with the correspondence between the two shown in Fig. 1. two reference sites we set c 0 = 1 as a normalization factor for the correlations.
To reduce the number of control parameters in our simulations we introduce a single parameter p such that c i = p > 1 for favored two-dimer configurations i, while c i = 1/p < 1 for unfavored configurations and p = 1 for cases that are considered "neutral". If all dimers are uncorrelated, i.e., p = 1, the state reduces to the standard AP state, for which we choose the single-bond amplitudes to be h(r) = 1/r 3 . This choice, which we keep also for the CAP states, is motivated by the fact that this gives the correct description of the Néel state. The generalized CAP states we use in our simulations are therefore characterized by the single parameter p controlling the bond correlations.
Below we investigate two different parametrizations of the bond correlations. In both cases we use c 1 = p > 1 to locally favor the columnar or plaquette VBS pattern (and whichever of these two VBS patterns that actually will be realized is not clear from the outset). Other types of dimer correlations are suppressed, by setting c 2 = c 3 = c 4 = 1/p in the first case-case (I)-while they are set to neutral, c 2 = c 3 = c 4 = 1, in case (II).
We destabilize Néel order by increasing the control parameter p and study the phase transition into a VBS. For case (I), we have found a first-order Néel to columnar VBS transition, while for case (II) we have found a continuous transition into a critical U(1) spin liquid, followed by a second continuous transition to the columnar VBS. We discuss the two cases in order.

A. A first-order Néel to VBS transition
Case (I) again corresponds to favoring parallel VB bond configuration by setting c 1 = p > 1 in Fig. 1 and suppressing fluctuations by setting c i>1 = 1/p. As in the 1D case discussed in Sec. III, we here first use the Binder cumulant of the staggered magnetization U to detect the transition of the Néel order. This quantity is also useful for distinguishing between a first-order and continuous phase transitions. As shown in Fig. 7(a), the Binder cumulant as a function of the control parameter p exhibits a minimum separating the Néel phase, where U → 1, and a non-magnetic phase, where U → 0. The minimum value of U is negative for all system sizes we studied, and the negative peak becomes narrower and deeper as the system size increases. In fact, the negative peak diverges as −L 2 when L → ∞, as plotted in Fig. 8, which provides strong evidence 48 for a first-order phase transition.
Note that the divergence of U min of the form L d expected for a classical d-dimensional system could in principle change to L d+z for a quantum system, where z is a dynamical exponent. 49 However, our definition of the Binder cumulant is purely a real-space definition and does not include integration over the imaginary time dimension (which we do not even have access to because it relies on a path-integral formulated using the underlying, unknown parent Hamiltonian). The form L 2 therefore is expected.
To check whether the non-magnetic phase exhibits VBS order, we next compute the columnar component of the dimer correlation defined according to (20). Fig. 7(b) shows D * (r) for the largest distance, r = (L/2, L/2). The correlation indeed converges to a non-zero value in the non-magnetic phase, tending to a step function at p c as L increases. The location of the discontinuity coincides with the point where the Binder cumulant reaches its minimum in Fig. 7(a). Note also that the curves for different L cross each other. This size-independence of the order parameter (as supposed to size-independence after multiplying with some power of L corresponding to an exponent of critical correlations) at the transition also supports a first-order scenario.
The location of the transition point p c in the thermodynamic limit can be determined, e.g., by extrapolating the U minimum location p min (Fig. 7) to the infinite-L limit. The finite-size scaling plot in Fig. 9(a) shows that the finite-size shift of the transition point defined in this way is consistent with ∝ L −2 for large L, where the exponent 2 again is the one expected based on scaling at a first-order transition, as discussed above. We estimate p c ≈ 1.500 from an extrapolation to L → ∞.
In the regime for p < p c the cumulant for different system sizes exhibits crossing points versus p. We expect that the crossing points should coincide with the minimum location when L = ∞. By finite-size extrapolation of crossing points p * for pairs of two system sizes L/2 and L, shown in Fig. 9(b), we estimate p * c ≈ 1.500 in the thermodynamic limit, in perfect agreement with the result obtained from the cumulant minimum.
In addition to the Binder cumulant signaling the transition of the Néel order, we also estimate the transition point of the VBS order from the scaling of the crossing points p ′ of the long-distance dimer correlation function. The inset of Fig. 7(b) shows a zoom-in of the region where the crossings occur. As shown in Fig. 9(c), these crossing points also appears to shift as p ′ ∼ L −2 for large L, and the extrapolation to L → ∞ yields an estimated location of the transition point p ′ c ≈ 1.503. This is marginally above the two other estimates discussed above, but given the very small range of data points for which the L −2 fits work well (we have used the points for the three largest systems in all cases, but the data still show some nonasymptotic curvature here), this result is still consistent with a single Néel-VBS transition point. Larger system sizes would be required to extract the location of this point more precisely.
Finally for case (I) we examine the histogram P (d x , d d ) of the order parameters d x and d y , corresponding to VBS order with x-and y-oriented dimers, Eq.  a circular shape with a central peak. At p = 1.48 in the transition region, the distribution shows coexistence of the Néel order (characterized by the central circular region) and the columnar VBS order [characterized by the four narrow peaks at angles φ = arctan(d y /d x ) = 0, π/2, π, 3π/2]; this again provide clear evidence for a first-order Néel-VBS phase transition. At p = 1.50, only barely inside the columnar VBS phase of this finite system, the distribution exhibits only the four VBS peaks.
From the many consistent results discussed in this subsection, we can conclude that the CAP states with favored parallel dimer pairs and suppressed flips of such pairs can characterize a first-order phase transition between the Néel and the columnar VBS phases. In such a CAP state, the Néel order is destroyed by formation of parallel dimers as the weight p increases. A first-order transition is also of course what would normally be expected for an order-order transition involving two unrelated order parameters. Note also that the system sizes studied in this section were rather modest; up to L = 36 (while much larger systems, up to L = 128, will be considered in the next section). For larger systems it becomes very difficult to obtain good statistics and smooth curves versus the control parameter, because of hysteresis effects related to the first-order nature of the transition. Still, as we have shown, the system sizes studied are sufficient to study the asymptotic finite-size scalimng forms.
In the context of VBS ground states of Hamiltonians, a first-order transition was previously observed with a J-Q model with the multi-spin interaction Q arranged to favor a staggered state. 50 In that case, the first-order transition was expected, because local dimer fluctuations are strongly suppressed with this kind of bond order. Alternatively, one can make an argument based on the nature of vortex-like defects in the VBS. 51 In contrast, in a columnar state parallel short-dimer pairs can fluctuate by 90 • rotation, unless such fluctuations are energetically expensive. In the DQC theory, these fluctuations are essential and correspond to an emergent U(1) symmetry of the VBS order parameter, which has been confirmed in J-Q models with plaquette VBS ground states. 33,35,36,52 In the CAP state considered here, we suppressed the fluctuations out of the perfect columnar state by the use of correlation weights and the observed first-order transition then is in line with expectations based on the DQC theory and the earlier studies of the Néel-VBS transitions in various models.

B. A critical U(1) spin liquid and a second-order transition to the VBS
Given the findings and discussion in the previous section, we now investigate whether the removal of the correlation factors suppressing dimer fluctuations can change the nature of the Néel-VBS transition of the CAP states. In this case (II) we only set c 1 = p > 1 and keep the other correlation factors in Fig. 1 as neutral; c 2 = c 3 = c 4 = 1. Fig. 11 shows the Néel Binder cumulant and the staggered dimer correlation function against the control parameter p. Like in case (I), the Néel order, characterized by the Binder cumulant tending to 1 as L grows, survives in the small p region up to a certain value of p, and a substantial columnar dimer correlation D * sets in when the Néel order is destroyed.  Fig. 11) increases only logarithmically with the system size L, indicating a continuous transition.
We notice that the negative peak of the Binder cumulant, which occurs only for large systems, is less pronounced than the cumulant peak in case (I). A negative Binder cumulant is often taken as evidence for a first-order transition, 53 but there are now known examples of rigorously understood continuous classical phase transitions associated with this behavior, e.g., the 2D 4state Potts model and the related Ashkin-Teller model in the neighborhood of this special point. 54 In such cases of "pseudo-first-order" scaling the minimum U min of the Binder cumulant diverges much slower than the expected L d form at a first-order transition (or, in some cases, possibly even converges to a finite value) .
The finite-size scaling plot Fig. 12 shows that −U min for CAP in case (II) grows only logarithmically with the system size. Thus, we conclude that the Néel order here vanishes in a continuous transition associated with pseudo-first-order behavior (related to anomalies in the critical order-parameter distribution).
Next we determine the phase boundaries more precisely. For a continuous phase transition, a frequently used method to determine the critical point is to find the unique asymptotic crossing point of the Binder cumulant for different system sizes. Due to finite size effects, the cumulant curves will intersect at a single point only when the system sizes are sufficiently large. For our case, we find the intersection point of the cumulant for L ≥ 80, and it is located at p ≈ 1.276 [ Fig. 13(a)], which can be identified as the Néel-spin liquid transition point p 1c . As the cumulant for a large system size exhibits a negative peak before it vanishes in the spin liquid phase, we also extrapolate the location of the peak to L → ∞ to determine the critical point from another route. By doing so, we estimate p 1c ≈ 1.276, in perfect agreement with the cumulant intersection point.
Beyond the order of the transition, another major dif- Finite size scaling of the location pmin at which the cumulant reaches its minimum. The critical point of the Néel-spin liquid transition in the thermodynamics limit is also here estimated as p1c ≈ 1.276, from an extrapolation to L = ∞ ference from case (I) is the behavior of the dimer correlations. As seen in Fig. 11(b), there is a wide intermediate region between the Néel phase and the VBS phase for larger p (where the correlation clearly converges to a non-zero constant for large systems); in this intermediate region, the dimer order parameter still decays versus the system size and, for the largest systems, the curve versus p develops a behavior suggestive of a second phase transition between a disordered state and the VBS above p = 1.6. The dimer correlations decay in the intermediate region as a power law, D * (r) ∼ r −β , with the exponent β depending on p, as shown in Fig. 14. Algebraically decaying dimer correlations were previously found in shortbond resonating valence bond (RVB) spin liquids investigated in Refs. 21,22. We thus tentatively identify this intermediate p-region with critical dimer correlations as a spin liquid in the same class of RVB states, for which it is known that the exponent of the dimer correlations depends on details of the bond fugacities. 21 We next study the distribution of the VBS order parameter P (d x , d y ). The examples of distributions shown in Fig. 15 for L = 64 are ring shaped in the intermediate region (at p = 1.4 and p = 1.5) before evolving into the expected Z 4 -symmetric shape in the larger-p regime where columnar VBS order is formed (as seen in the location of the peaks; a 45 • rotated distribution would correspond to a plaquette VBS). The same kind of ring-shaped P (d x , d y ) distribution was also found in the prototypical short-bond RVB spin liquid (although only when the bond configurations are restricted to the dominant topological sector of zero winding number). 21,22 An important issue here is whether the VBS hosts emergent U(1) symmetry when the critical point is approached. To investigate this, we also need to determine the location of the point where the VBS becomes longrange ordered. This is not easy to do just based on the scaling behavior versus the system size in Fig. 14, because there is a whole critical VBS phase and the change from the power-law decay to convergence to a samall non-zero constant is subtle. We can instead characterize the VBS state by the quantity cos(4φ), where φ = arctan(d y /d x ) is the angle in the VBS order parameter (d x , d y ) computed for the individual VB configurations, i.e., based on the histogram, P (d x , d y ). This expectation value measures the degree of the developed Z 4 symmetry of the VBS order parameter. In the spin liquid (as well as in the Néel state), where the U(1) symmetry is preserved as L → ∞, we have cos(4φ) = 0. In the VBS states in which the distribution for large systems develops Z 4 symmetry, we have cos(4φ) → 1 as L → ∞ for a columnar VBS (while it approaches −1 for a plaquette VBS). As cos(4φ) is dimensionless, one expects cos(4φ) -curves for different system sizes to asymptotically become size-independent at the transition point into the VBS state.
As seen in Fig. 16(a), for the case-(II) CAP, a crossing point develops at a non-zero value of cos(4φ) , at p ≈ 1.65, which we thus identify as the liquid-VBS transition point. The fact that cos(4φ) > 0 at this point shows that there is no emergent U(1) symmetry at the VBS-liquid transition, since the order parameter remains Z 4 symmetric exactly at the transition point. For ref- , respectively], the distributions are ring-shaped with weight at all angles (bright regions). As p increases to p = 1.6 (c) and p = 1.7 (d), the U(1) symmetric distributions evolve into Z4-symmetric ones, with higher densities at the angles 0, π/2, π, 3π/2 corresponding to columnar VBS order. erence, in Fig. 17 we show QMC results for the J-Q model with 6-spin columnar Q-interactions, for which the order-parameter symmetry was previously analyzed in a slightly different way. 36 The results for cos(4φ) here show crossing points decaying toward 0 in the vertical direction. The system sizes are not yet sufficiently large to see that the crossings tend toward the critical point, (J/Q) c ≈ 0.66. The contrast with the CAP states in Fig. 16(a) is stark, however, with the absence of four-fold symmetry-presence of emergent U(1) symmetry-at the transition being very plausible.
Since the results shown in Fig. 16(a) appear to give a rather precise estimate for the transition point, p 2c = 1.650(5), without any scaling needed of the vertical axis, we now use this result to investigate scaling of the VBS order parameter. In Fig. 18, in order to achieve data collapse, we have rescaled the long-distance dimer correlation functions in Fig. 13(b) both vertically, multiplying by L β , with β = 0.22, and horizontally, multiplying (p − p 2c ) by L 1/ν , with 1/ν = 0.44. The correlationlength exponent is, thus, ν ≈ 2.3, which is anomalously large and definitely rules out a first-order transition (in which case ν = 1/d = 1/2 would be expected).
The cos(4φ) curves in Fig. 16(a) can also be scaled in the horizontal direction to collapse all the data onto a single curve, as shown in panel (b). In cases where there is emergent U(1) symmetry, this procedure, using the control parameter scaled as L 1/aν (p − p c ), gives the correlation-length exponent ν multiplied by a number a > 1, 36 reflecting the faster divergence of the lengthscale Λ controlling the emergent symmetry; Λ ∼ ξ a . In the case at hand, we have already concluded that there is no emergent U(1) symmetry, since cos(4φ) remains non-zero as p → p 2c , and the exponent aν, with a < 1 should instead reflect a shorter length-scale (irrelevant operator), governing the reduction of the angular VBS fluctuations ( cos(4φ) → 1) in the VBS and the growth of these fluctuations ( cos(4φ) → 0) in the liquid. We find very good scaling, and indeed the factor a ≈ 0.8 is clearly less than one. We finally investigate the nature of the spin correlations in the spin liquid phase. Results for the squared sublattice magnetization for several points representing the three different phases are shown in Fig. 19. Here we plot the results on a log-log scale, in order to study power-law correlations. In the Néel state the sublattice magnetization approaches a non-zero constant, while in the liquid and VBS states we observe a clear 1/L 2 decay. This is the form expected with exponentially decaying spin-spin correlation functions. We see a power-law behavior with a non-trivial exponent, ∼ L −α , with α = 1. only at the Néel-liquid critical point. These results again confirm that the liquid is of the same type as the prototypical AP RVB states (i.e., with no correlation factors), where the varying power-law for the critical VBS correlations corresponds to a varying stiffness constant in a mapping to a height model. 21,55 V. VARIATIONAL CALCULATIONS In this section we explore variational optimization of CAP states, carrying out energy minimization based on derivatives along the same lines as in Refs. 19,20. We consider two models. First, the standard Heisenberg model, 20 50 100

defined by the Hamiltonian
where ij denotes a pair of nearest-neighbor sites. We will consider both the 1D chain and the 2D square lattice (in both cases adopting periodic boundary conditions for systems an even number of spins). We also consider the J-Q model, 33 which includes four-spin interactions in addition to the exchange J: Here ij and kl form opposite edges on an elementary 2 × 2 plaquette on the square lattice and the summation includes both horizontal and vertical orientations of these edges on all plaquettes (i.e., the Hamiltonian obeys all the symmetries of the square lattice). With the negative prefactor of the Q term, this interaction clearly is related to enhancement of the parallel-dimer weight c 1 , in the notation of Fig. 1, in the ground state wave function (although the state is still of course not expected to be exactly reproduced by the CAP ansatz).

A. Optimization method
We start a variational calculation from some initial value of the parameters in the CAP state (7), typically a power-law form of the amplitudes h and all the correlation constants C = 1. When optimizing states for different values of some parameter (e.g., J/Q), we also normally start the calculation for a new parameter value from a previous calculation for some nearby value. One can also use this approach for different system sizes, although when increasing the system size initial values for the parameters corresponding to the longest bonds are of course not available and have to be set to some suitable values based on the longest previous bonds. In general, we have found that the starting point does not play an important role in optimization of AP and CAP states, indicating that the energy landscape is relatively smooth.
To minimize the energy as a function of all parameters, we compute the energy and its derivatives. We then apply either (a) the steeped-decent method or (b) a stochastic variant of it where only the signs of the derivatives are used, as discussed in Ref. 20. A generic parameter p is in these two cases updated according to Here, in (a) R ∈ [0, 1) is a random number and in (b) max(|dE|) denotes the derivative that is the largest in magnitude among all the derivatives considered. The maximum shift δ is gradually reduced, so that the variational parameters eventually will converge. If δ is reduced sufficiently slowly, then one will reach a minimum. This minimum is not necessarily the global one, however. The stochastic scheme (a) should be better in avoiding local minimums, although one can of course never be completely guaranteed to find the global minimum. For the case at hand, the energy landscape appears to be relatively smooth, with no serious problems in consistently reaching the same minimum energy (regardless of the starting point, as mentioned above). Occasionally the method does fail, with independent runs of the same system leading to different final results. Typically, in a set of several runs, most of them will be consistent with each other, with only a small fraction of them deviating significantly from the majority value. As expected, the failure rate decreases when increasing the number of MC sweeps for sampling the VB configurations (leading to smaller error bars on the derivatives) and when reducing δ at a slower rate. We typically use a protocol based on an iteration number k. For each k, the parameters are adjusted some number M of times (e.g., M = 100) based on derivatives obtained in MC simulations with ∝ k 2 steps. This way, the derivatives become more precisely determined as the solution approaches the minimum (where the derivatives decrease in magnitude, thus necessitating a larger number of MC sampling steps to obtain statistically useful information). The maximum parameter shift δ in Eqs. (29) is of the form δ 0 /k α , with α = 3/4 a suitable exponent in practice, based on experience. In practice, one can quite easily converge calculations for small systems essentially completely, while large systems require very long runs. There is still some room for improvement of the energy minimization protocol, as what we have explored so far are essentially schemes based on trial-and-error approaches. Nevertheless, the results to be presented next can be considered as almost optimized and we do not anticipate that our conclusions would change based on more complete optimizations.

B. Heisenberg chain
An interesting question in one dimension is whether AP or CAP states can describe the critical ground state of the Heisenberg chain. As we saw in Sec. III, with the simple parametrized 1D AP wave function (25) the correct critical decay exponents (α = β = 1) 44 corresponding to this system cannot be achieved. In a variationally optimized state, we are not tied to any particular form of the amplitudes, and a sufficiently flexible variational wave function should then be able to capture the cor- rect criticality, including the logarithmic corrections that arise in the field-theory language due to a marginally irrelevant operator. The question then is whether the AP or CAP states have this kind of flexibility, to possibly capture even such a subtle effect as the logarithmic corrections to the correlation functions. 57 To answer the above question, we have carried out energy minimizations with the simple AP state (with all amplitudes as adjustable parameters) as well as with two types of CAP states. To include only the minimum amount of bond correlations beyond the AP state we include in (7) only the two 1D bond-pair configurations with length-1 bonds; C(r 1 , r 2 ) with r 1 = ±1 and r 2 = ±1. The case r 1 = 1, r 2 = −1 corresponds to a single bond on a nearest-neighbor link and we can regard this as a normalization for the correlation factors, C(1, −1) = 1, in the same way as we also set h(r = 1) = 1. There is then only one other correlation factor, C(−1, 1) to optimize at this level. We also consider the extreme case of optimizing all C(r 1 , r 2 ), r max = L/2 − 1, again with C(1, −1) = 1.
Let us first discuss the energy. As an example, for L = 256 the exact energy per site is E/L = −0.44316 while for the AP state we obtained −0.44184. For the CAP the best energy when r max=1 is E/L = −0.44272, and with r max = L/2 − 1 it decreases to −0.44306. As discussed in the previous section, it is difficult to completely optimize long chains, so the optimal variational energies may still be somewhat lower. Following the trends as a function of system size, the relative energy error with the CAP state seems to remain at about 0.05%.
Turning now to the spin correlations, in Fig. 21 the staggered spin structure factor, defined according to is graphed versus the system size for all the cases discussed above. Since the exact C(r) decays as 1/r with a multiplicative logarithmic correction, the exact S(π) grows slightly faster with L than ln(L), as demonstrated with unbiased QMC results in Fig. 21. All the AP and CAP results exhibit a faster growth with L. When graphed on a log-log scale (instead of the lin-log scale used in Fig. 21), the AP behavior is consistent with a power law, S(π) ∼ L 1−α with α ≈ 0.70. With the CAP states, the data move closer to the exact points, but even with the maximally correlated CAP state the divergence is still somewhat too fast. It is also interesting to examine the optimized amplitudes of the AP state. Fig. 22 shows results for L = 256. Interestingly, a power law applies here for short and moderate bond lengths, with the deviations (enhancements) at large lengths likely related to the periodic boundary conditions (and some jaggedness of the large-r data due to imperfect optimization, reflecting the total energy not being very sensitive to these "noise" features). Even the r = 1 amplitude falls on the common power-law line in Fig. 22, i.e., in the notation of Sec. III the optimized state has λ = 1. Looking at Fig. 6, when λ = 1 the exponent α ≈ 0.75, quite close to α ≈ 0.70 obtained above with the optimized amplitudes. Thus, the boundary effects on h(r) seen in Fig. 22 appear to have only minor effects on the critical behavior. The conclusion for the optimized AP state is, thus, that a critical behavior is reproduced, but with the wrong exponents for the correlation functions. Note, however, that α ≈ β for the applicable power-law obtained here, which is also the case for the true Heisenberg correlations (but with larger values, α = β = 1).

C. Two dimensions
We next systematically investigate the improvement of the energy with the inclusion of bond correlations in two  (7) is optimized (i.e., can be different from 1) are those marked by n.
dimensions, using several choices for the maximum bondlength r max in the correlation factors C(r 1 , r 2 ). Fig. 23 illustrates all the bond shapes (r 1 , r 2 ) at three correlation levels, with r max = 1, √ 5, and 3 for correlation levels 1,2, and 3, respectively.

Heisenberg model
For the 2D Heisenberg model, previous variational AP calculations have shown that the energy error within this class of state is < 0.1% for large systems, and the spin correlations are reproduces to within 1% or better. 19,20 Although the system is strongly Néel-ordered and only has rapidly decaying short-range VBS correlations, including bond correlations with CAP states can still significantly improve the energy further. Fig. 24 shows results for L × L systems with L = 16, 32, and 64 at different correlation levels. The deviation from unbiased QMC calculations decreases with increasing correlation level. For L = 16 with r max = 3 the relative error is as small as ≈ 4 × 10 −5 , while for the larger systems it is somewhat larger, about 10 −4 .
Going further and optimizing all correlation weights C(r 1 , r 2 ) with r ≤ L/2 − 1, one should in principle be able to further improve the energy and obtain the best possible CAP state (with the kind of correlations included here) when L → ∞. The energy only improves marginally on the r max = 3 results, however. Fig. 25 shows results versus the system size for the energy as well as the sublattice magnetization. On the scale of the graphs, one can barely see any differences between the CAP and unbiased QMC results for L ≤ 20, while for the larger systems there are some visible deviations.
Here it should again be noted that the results for large systems are likely not completely optimized. As discussed above in Sec. V A, the energy depends only very weakly on the long-bond statistics, which implies that MC evaluations of the corresponding derivatives are affected by relatively large fluctuations, leading to slow convergence. The sublattice magnetization is more sensitive to the long bonds, however, and this makes it very difficult to obtain completely unbiased results for large systems. For example, five independent optimizations for L = 32 with r max = 3 all gave the same energy within statistical errors, but the sublattice magnetization showed significant fluctuations, with the following results: m 2 = 0.1131, 0.1113, 0.1132, 0.1129, 0.1094, with the error bar approximately equal to 2 in the last digit.
Here it can be noted that three of the results agree well, while two of them are clearly off. One may then conclude that the best optimized results should be around 0.1131, although to confirm this one should carry out a much larger number of independent runs. The correct results based on unbiased QMC calculations 19 is m 2 = 0.1128, less that 0.3% below the average of the above 3 consistent points.

J-Q model
As discussed in Sec. I, the J-Q model (28) exhibits a Néel-VBS transition at a critical value of the coupling ratio J/Q, with most precise estimate so far being (J/Q) c = 0.0447 (2). 38 An interesting question is whether this transition can be described by the CAP states. Here we consider the case J = 0, where the ground state is a columnar VBS. This VBS is very complex, however, since the order parameter is only about 20% of the maximum possible value (i.e., for states with length-1 singlets forming columns and no fluctuations around this configuration). 52 The fluctuations are significant and have U(1) character up to a very large length scale (larger than what can currently be studied). As it turns out, the fully optimized CAP state (using r max = L/2−1) for this J = 0 system does not reproduce the VBS order. Instead, as we will see below, the system is still on the Néel side of the quantum phase transition. First, let us again investigate the impact of including longer bonds in the correlation factor in Eq. (7). The energy is improved very dramatically with increasing r max .
For example, for L = 32 the best optimized AP state has an energy E 0 /N = −0.8013. The CAP states have E/N = −0.8215, −0.8229, and −0.8232, at correlation levels 1, 2, 3 in the scheme of Fig. 23. Going to r max = 5 there is only a marginal energy improvement to −0.8233, which differs by about 0.1% from an unbiased QMC result; E/N = −0.8240.
Although an energy deviation of 0.1% would normally be considered excellent in variational calculations, the order parameters are still not well reproduced. Fig. 26 shows the dependence of both the Néel and VBS order parameters on r max . With increasing r max , the sublattice magnetization is reduced and the VBS order parameter increases, as would be expected with CAP states in a VBS state. However, the VBS order parameter is still more than 30% too small at r max = 5, and it appears to be essentially converged at that point. Accordingly, the Néel order parameter is instead too large. Fig. 27 shows both order parameters calculated with r max = 3 as a function of the inverse system size, along with unbiased QMC results. Here one can observe that the agreement between the two calculations is very good for small systems, but the agreement becomes worse for increasing L. Asymptotically, the variational calculations tend toward a weakly Néel ordered state, with no VBS long-range order, as opposed to the actual VBS ground state. This calculation serves to illustrate the insensitivity of the energy to long-distance correlations and the related difficulty in using the energy as a reliable measure of the quality of a state obtained by variational means (see Ref. 58 for a different example of this issue).
While this result could be seen as a failure of the variational CAP states, it should be noted again that the actual Néel-VBS transition takes place at a very small coupling ratio, (J/Q) c ≈ 0.045, and one cannot expect a variational calculation to reproduce the critical point exactly. For the CAPs considered here we have confirmed that the transition is pushed to a small negative J/Q, but we leave more detailed studies of the transition (which requires very well optimized states) for a future study.

VI. SUMMARY AND DISCUSSION
In conclusion, we have discussed VBS states and associated quantum phase transitions in 1D and 2D wave functions in the VB basis. VBS states appear naturally within the standard 1D AP states, and we have here characterized the continuous Néel-VBS transition in such a class of states with amplitudes decaying as a power-law of the bond length. To stabilize a 2D VBS requires explicit bond correlations beyond the AP states. We have introduced the CAP states, where correlations are enforced through factors corresponding to bond pairs, as in the wave function (7). We have shown how tuning of parameters in 2D CAP states can lead to transitions from the standard Néel antiferromagnet to a VBS, in some cases with an intervening spin liquid phase. With the parametrization considered here, the direct Néel-VBS transition is first-order, while the Néel-liquid and liquid-VBS transitions are continuous. The 2D Néel-liquid transition is of the same kind as in the pure AP states, although the short-distance VBS fluctuations are enhanced. Interestingly, the liquid-VBS transition is not associated with an emergent U(1) symmetry, although the columnar VBS (which is the VBS variant stabilized with the CAPs studied here) in principle supports this phenomenon 26 and has been observed in QMC studies of J-Q models at the Néel-VBS transition. 33,35,52 Thus, in the states we have studied here the VBS should be induced by a relevant operator, instead of the dangerously irrelevant operator associated with the emergent U(1) symmetry in the "deconfined" criticality scenario.
It remains an interesting challenge to find a parametrization of the CAP states such that a DQC point is obtained. In the study with short-bond correlations in Sec. IV we kept the single-bond amplitudes fixed with the form h(r) = 1/r 3 and varied only a correlation parameter p. In principle we could also use h(r) = 1/r κ and also vary κ. It would be interesting to study the full phase diagram in the plane (p, κ) for the two parametrizations of the bond correlations, and also with different choises of the correlation factors.
We note a recent study of AP states, with a certain parametrization of the amplitudes, to describe the Néel to quantum-paramagnetic transition in the bilayer Heisenberg model. 59 The correct 3D classical Heisenberg exponents were obtained for this transition, i.e., the AP states contain effectively long-range interactions that allow (2+1)-dimensional criticality to be reproduced within a 2D configuration space. In principle it seems that this should be possible to achieve also for the DQC transition, and if such a program to construct a simple CAP wave function is successful, it would likely lead to further useful insights into the mechanism of spinon deconfinement and emergence of the effective U(1) gauge field.
It is possible that CAP states with all parameters adjusted to minimize the energy of a model such as the J-Q model could lead to a DQC point. Here we carried out such variational QMC calculations for the standard J-Q model with four-spin interactions. This model has a critical point very close to J/Q = 0, however, and in the variational calculation the J = 0 system is still inside the Néel phase. In principle one can still study the phase transition by going to negative J, but in that case Marshall's sign rule cannot be proven rigorously (although most likely it should still hold when J/Q is small and negative). We will investigate this case in future studies, and also consider the model with sixspin interactions, 36 where the transition point is at much larger J/Q and should remain well within the range of positive J/Q values within optimized CAP states. If indeed the correct type of criticality can be achieved, then by examining the bond amplitudes and correlation factors in the optimized state it may also be possible to construct a class of CAP states with a single tunable parameter (instead of all the parameters changing as a function of J/Q in the variationally optimized states) to drive this type of criticality-which the parametrization used in the present paper was not capable of.
Beyond the case (I) and (II) parametrizations of the CAP states that we considered here, we have also explored other cases. In particular, with c 1 = c 3 = p and c 2 = c 4 = 1/p, in the notation of Fig. 1, we have found a plaquette VBS, in contrast to the columnar VBS states obtaining in cases (I) and (II). In future studies it will also be interesting to study the nature of the transition from a Néel antiferromagnet to a VBS in this case.
It is in principle possible to further improve on the CAP states, by including more correlations. We here considered pairs of bonds connected to a nearest-neighbor link only. We have verified that this gives a better variational energy (for the Heisenberg and J-Q models) than next-nearest-neighbor links, and therefore, most likely, the nearest-neighbor links are optimal for introducing correlations in this kind of CAP states. One could also combine several types of correlation factors, and in-clude also factors for correlations between more than two bonds. In practice this may not be worth the effort, however, as the main utility of CAP states should be (i) to have a simple class of states to capture the Néel-VBS transition and (ii) to use them as "trial states" for projector QMC calculations in the VB basis. 19,[60][61][62] While the variational states can be improved, in practice it is better to project out the ground state exactly using QMC if completely unbiased results are needed, and too many parameters in a CAP state defeats the purpose of (i).
A very interesting question is whether the phase transitions we have discussed here can be realized in ground states of reasonable Hamiltonians, with only local interactions. In particular, the Néel-liquid-VBS series of phases would be of interest in this regard. We already know from the worh of Cano and Fendley on the shortbond RVB that there is a local parent Hamiltonian for that state. It is then also appears very plausible that some local pertutbations of this Hamiltonian will effect the stiffness constant characterizing the RVB 21,55 and governing its critical VBS correlations. Thus, a class of local Hamiltonians should be able to capture the whole spin liquid phase in our case (II). Then, it also seems plausible that other local perturbations can drive the system into a Néel or VBS states, e.g., the Q term of the J-Q model should do this. Since the Cano-Fendley Hamiltonian has a sign problem in QMC calculations, some other methods would be needed to study phase transitions in perturbations of it.