Hidden orders and phase transitions for the fully packed quantum loop model on the triangular lattice

Quantum loop and dimer models are prototypical correlated systems with local constraints, which are not only intimately connected to lattice gauge theories and topological orders but are also widely applicable to the broad research areas of quantum materials and quantum simulation. Employing our sweeping cluster quantum Monte Carlo algorithm, we reveal the complete phase diagram of the triangular-lattice fully packed quantum loop model. Apart from the known lattice nematic (LN) solid and the even $\mathbb{Z}_2$ quantum spin liquid (QSL) phases, we discover a hidden vison plaquette (VP) phase, which had been overlooked and misinterpreted as a QSL for more than a decade. Moreover, the VP-to-QSL continuous transition belongs to the $(2+1)$D cubic* universality class, which offers a lattice realization of the (fractionalized) cubic fixed point that had long been considered as irrelevant towards the O($3$) symmetry until corrected recently by conformal bootstrap calculations. Our results are therefore of relevance to recent developments in both experiments and theory, and facilitate further investigations of hidden phases and transitions.

A QDM where two dimers touch every lattice site is also dubbed the fully packed quantum loop model (QLM) [33,34,[48][49][50][51]].In the strongly interacting limit, the QLM shares the Hilbert space and low-energy properties of a hard-core boson model (or spin-1/2 XXZ model) on a kagome lattice at 1/3 filling and many other similar frustrated models with local constraints [33,34,52,53].The QLM is also related to the resonant-valence-loop phase in frustrated spin-1 models [54,55], in analogy to the relation between the one-dimer-per-site QDM and spin- 1  2 models.The triangular-lattice QLM has been investigated by field-theoretical analyses and various numerical approaches [33,34], but, so far, no consistent phase diagram has been obtained.At the RK point, the model realizes a Z 2 quantum spin liquid (QSL) phase, which has even Z 2 topological order and hosts fractionalized quasiparticles called visons.Tuning away from the RK point, the model also The triangles v1 and v2 represent two sublattices for the visons, and the red and grey colors in the triangles denote vison numbers (±1) with opposite sign.The first-order phase transition between the LN and VP states occurs at V = 0.05 (5).The first row of the middle subfigure is the kinetic energy correlation pattern and the second row is the vison plaquette (VP) pattern based on QMC simulation results, respectively (see Fig. S3(a,b) in Supplementary Note 3 for values in each triangle).VP is a hidden dimer solid phase without dimer order.The right subfigures illustrate the representative dimer coverings in the Z2 QSL phase.The continuous phase transition between the VP phase and the QSL occurs at Vc = 0.59 (2).
realizes various topologically trivial symmetry-breaking phases.The transitions from the QSL to the topologically trivial phases can be understood as condensation of the visons.Previous numerical studies [33,34] show that decreasing the repulsion between parallel dimers leads to a lattice nematic (LN) solid phase.Furthermore, the phase transition between the QSL and the LN phase, described by vison condensation, is suggested to have an emergent O(3) symmetry and therefore, belongs to the 3D O(3) * universality class; the * conveys that the transition is between the O(3) symmetry-breaking phase and a Z 2 topologically ordered phase (rather than a trivial symmetric phase).
These two works [Refs.[33,34]] employ different methods, namely, exact diagonalization (ED) of small systems and density-matrix renormalization group (DMRG) calculations [33], or variational quantum Monte Carlo (QMC) simulations with intermediate system sizes of 12 × 12 [34].The numerical evidence presented in these two works, however, disagree on the boundaries of the phase diagram.In particular, the former indicates that the LN-QSL transition occurs at V /t ∼ −0.3 whereas the latter estimated this point to be at V /t ∼ 0.7 and found that the LN order is vanishingly small between V /t ∈ (0.1, 0.7).As we will explain below, this is because there is a hidden ordered phase between the LN solid and the QSL.
Theoretically, the speculation of the LN-QSL transition being in the O(3) * universality class is based on the assumption that the cubic anisotropy is irrelevant at the 3D O(3) Wilson-Fisher fixed point [56][57][58][59][60].However, this assumption has been disproven by recent conformal bootstrap analyses [61,62] demonstrating that the 3D O(3) fixed point is unstable and will flow towards the nearby corner cubic fixed point, whose symmetry-breaking phase is the corner cubic phase.We will show later that there indeed emerge cubic anisotropies in this model which change the O(3) criticality suggested in previous works to a cubic one.The cubic phase transition occurs between the hidden order and the QSL.
Moreover, a previous study [52] on the honeycomblattice transverse-field Ising model (which is related to the current QLM by representing the vison excitations with Ising spins) shows a first-order transition instead of a continuous one, with clear face-centered cubic anisotropy at the transition point.Therefore, resolving the controversy in the phase diagram of the triangularlattice QLM requires an additional and consistent understanding that unifies the aforementioned recent developments [33,34,52,61], and provides the correct mechanism of the vison-condensation transition from the Z 2 QSL to the solid phases.The results in Ref. [34] and, in particular, the vanishing of the LN order between V /t ∈ (−0.3, 0.7), already hints at such a reconciliation.This first-order phase transition can also be clearly understood in our work, as we will show below, in that it occurs between the LN and the hidden-order phases.
In summary, our work addresses these longstanding questions and presents a comprehensive picture of the quantum phases and phase transitions of the QLM.Using both QMC and ED analyses, we find that between the proposed LN phase [63] and the Z 2 QSL phase, there exists a hidden vison plaquette (VP) solid state without apparent dimer order.The phase transition between the LN and VP phases is first-order, and that from the VP to the Z 2 QSL phase is continuous and belongs to the (2+1)D cubic * universality class.Moreover, our numerical discoveries of the LN-VP first-order transition and the VP-QSL continuous transition provide a unified understanding that is fully consistent with the recent conformal bootstrap analysis of the relation between the O(3) and cubic fixed points [56][57][58][59][60] in (2 + 1)D [61].Our results for the QLM on the triangular lattice provide a nontrivial realization of this scenario in a setting that is experimentally accessible in programmable quantum simulators based on highly tunable Rydberg atom arrays [20,21,28,[64][65][66][67][68][69].

The constrained model
The Hamiltonian of the QLM on a triangular lattice is defined as where α represents the three possible orientations of all plaquettes of the triangular lattice, as shown in Fig. 1(a) for an example of the LN state.The local constraint of our model requires two dimers to touch every triangularlattice site in any configuration, thereby forming the fully packed quantum loops [34].The kinetic term is controlled by t, which changes the dimer covering of every flippable plaquette while respecting the local constraint, and V is the repulsion (V > 0) or attraction (V < 0) between dimers facing each other.The special RK point is located at t = V and has an exact Z 2 QSL solution [39].We set t = 1 as the unit of energy in our simulations.
We employ the recently developed powerful sweeping cluster quantum Monte Carlo method [28,36,47,70,71] to solve this model.Our simulations are performed on the triangular lattice with periodic boundary conditions and system sizes N = 3L 2 for linear dimensions L = 8, 10, 12, 16, 20, while setting the inverse temperature β = L.Further description of the QMC implementation and its advantages can be found in Supplementary Note 4.

Physical observables
To explore the phase diagram of the Hamiltonian in (1), we first measure the order parameter describing the soft vison modes [33,[72][73][74], which is given by where r runs over all the unit cells of the triangular lattice, as shown in Fig. 1.The vector ϕ = (ϕ 1 , ϕ 2 , ϕ 3 ) encapsulates the 3D order parameters of the visons, which are obtained from the Fourier transforms of vison configurations v i v j = (−1) N P ij , with N Pij being the number of dimers cut along the path P between triangular plaquettes i and j [36].The label (v 1,r , v 2,r ) denotes visons living at the centers of the two sublattices of the triangular plaquette at r [see the left subfigure in Fig. 1(b)].In our simulations, we choose a reference vison on a certain triangular plaquette, e.g., by setting v 1,r=(0,0) = ±1 on the upper triangle of the first unit cell.The vison configurations are obtained with respect to such a reference, from which we can assign values of (v 1,r , v 2,r ) for all r.The three momenta ), and M 3 = (0, 2π/3) correspond to the low-energy modes of the vison dispersion and are invariant under the projective symmetry transformations [33,74] as we show in Supplementary Note 1.The associated u j are (1, 1) T for M 1 and M 2 , and (1, −1) T for M 3 .The derivation of (2), the symmetry analysis of the 3D order parameter of the visons, and representative vison configurations in the VP phase are presented in the Supplementary Note 1.
According to the analysis of the effective critical theory [33,74], the low-energy effective Hamiltonian of the quantum loop model is dual to the ferromagnetic transverse field Ising model on the honeycomb lattice.In the large-external-field limit, the spin in the z direction fluctuates, representing the fluctuation of the boson number on each site, which corresponds to the Z 2 spin liquid.As the Ising coupling increases, the energy dispersion acquires its minima at three M points, which correspond to three patterns of the LN phase.Therefore, the transition between the Z 2 QSL and the nontopological phase can be regarded as the process of vison condensation, which can be described by the vison order parameter (ϕ 1 , ϕ 2 , ϕ 3 ) [33,52].We take the magnitude |ϕ| = ϕ 2 1 + ϕ 2 2 + ϕ 2 3 , plotted in Fig. 2, as the order parameter to roughly detect the solid phases and their transitions to the Z 2 QSL phase.As illustrated in Fig. 2, the order parameter |ϕ| clearly vanishes via a twostep process that sharpens with increasing system sizes.When V < 0, the order parameter is finite, denoting long-range order in the vison pattern corresponding to the LN phase, as shown in the left subfigure in Fig. 1(b); when V > 0, ϕ drops to a finite plateau, signifying another symmetry-breaking phase.It is worth noting that the order vanishes around V = 0 if we measure it through dimer correlations (e.g., see Fig. S5 of SI).Ww thus ask: why does the nature of the order seem so different when viewed in terms of dimers and visons?
To further reveal the nature of the symmetry-breaking 1.0 0.8 0.6 0.4 0.2 0 phases, we plot the histogram of the order parameter ϕ in Fig. 3.The distribution of ϕ is indeed different in the two phases: it is peaked at the six face centers of the cube in the phase for V < 0, and at the eight corners of the cube in the phase for 0 < V < 0.6.According to the symmetry transformation of ϕ in Eqs.(S8) and (S9) of the SI, a peak at the face center corresponds to breaking threefold rotational symmetry while preserving translational and twofold rotational symmetries.From the broken symmetries, we recognise that this is the LN phase.On the other hand, a peak at the corner is associated with breaking the translational symmetry in a manner that doubles the unit cell in both directions while preserving the rotational symmetries.Therefore, the ground state actually belongs to a plaquette-ordered phase with a 2 × 2 unit cell, which we refer to as the VP phase.As we will discuss later, the ground state possesses emergent hidden order in that it does not exhibit any symmetry breaking in the dimer-dimer correlations but does so in other correlation functions.The LN and VP ordered phases break incompatible symmetries and thus cannot be connected via spontaneous symmetry breaking, which points to a first-order phase transition between them.The histogram in a region of phase coexistence, plotted in Fig. 3 (c), also indicates the first-order nature of the transition.More detailed results on the evolution from the LN to the VP solid phase can be seen in Supplementary Note 2.

Hidden order: the vison plaquette phase
It is only natural to ask why the vison plaquette phase had not been identified for more than a decade [33,34].This is likely because, except for the vison order parameter histogram in Fig. 3, one cannot find a corresponding dimer order in the VP phase, i.e., the order is hidden in the dimer basis.In order to uncover this hidden order in the VP phase, we perform the following analysis.
Motivated by the QMC data, we use exact diagonalization (ED) to further demonstrate that the hidden order cannot be measured in the dimer basis.As there is no clear peak in the structure factor of dimer correlations in the VP phase for all the system sizes simulated, we seek assistance from the vison histograms to amplify any signals of possible symmetry-breaking phases in the dimer patterns.In the VP phase, the eight cubic fixed points stay in eight octants in the histograms [such as in Figs.3(d), (e), and (f)], so every vison configuration can be classified into a certain octant except those on the boundaries; since two vison classes in opposite octants correspond to the same dimer configuration, there are only four classes of dimer configurations.We classify QMC dimer configurations into these four classes and average the ones in the same class (which would amplify the symmetry breaking in the dimers, if there is any) to obtain the dimer density on the strongest vison-ordered configurations.We collected about 500, 000 (300, 000) such configurations for an L = 4 (L = 8) system at V = 0.4 but found that the real-space dimer density is homogeneous, ∼ 1/3, on all the bonds of the lattice.Furthermore, as shown in Supplementary Note 6 (and Fig. S4), a similar analysis using the ground-state wavefunction from ED of a 4 × 4 lattice at V = 0.3 (which eliminates the statistical error in QMC simulations) also confirms that there is no translational symmetry breaking in the dimer density distribution.Therefore, despite evidence of translational symmetry breaking in Fig. 3, both QMC and ED results strongly suggest that the VP order is a hidden-order phase in the dimer basis, with a homogeneous dimer occupation.We note that such a unique hidden symmetry-breaking phase has not been reported earlier in the QDM literature.Whether there are intricate (emergent) symmetries that actually protects the homogeneity is an interesting question for future investigation.
This unsuccessful detection in the dimer basis can also be understood easily from the vison configurations.In fact, from the real-space vison correlation function in Fig. S3 (b) of Supplementary Note 3, the difference in correlation functions between two closest triangles is seen to be a constant (about 0.14 at V = 0.3).It is well known that this difference can be translated to the dimer occupation on the bond separating two neighboring triangles [75].Thus, the constant difference in vison correlations also points to a uniform dimer occupation.It is worth emphasizing that a homogeneous dimer density only requires a constant difference between closest visons The structure factor at the M point, CTq=M , as a function of V ; in the VP phase, CTq=M acquires long-range order.The inset shows that the VP phase (for V = 0.1, 0.15, and 0.3, the extrapolated value at L = ∞ is 0.006(2), 0.008(2), and 0.0097(3), respectively) persists in the thermodynamic limit.The error bars here represent the standard error of the mean, which is calculated as σ/ √ N , where σ is one standard deviation and N is the total number of independent samples.instead of a constant vison density.This gives rise to the possibility that this VP order exists separately from the QSL.
Considering that the vison is a fractional quasipaticle which cannot be measured in experiments while the QDM and QLM are widely investigated in frustrated quantum materials and blockaded ultracold atomic arrays, we now propose a measurable observable to identify this hidden order.We notice that the translational symmetry breaking indicated by the histograms in Fig. 3(d), (e) and (f) is indeed reflected in other correlation functions.Firstly, this can be seen from the vison correlations ⟨v i v j ⟩, as displayed in the phase diagram in Fig. 1 and discussed in detail in Supplementary Note 3. We find that although the vison correlation function appears to change sign under mirror reflections and sixfold rotations, the physical state actually preserves these symmetries, because of the two-to-one correspondence between vison and dimer configurations.Furthermore, the translational symmetry breaking can also be detected from correlation functions of local operators.
Therefore, we construct the following composite order parameter which can be measured in experiments: where the sum runs over the kinetic t-terms of the Hamiltonian on three rhombi with labels 1, 2, 3 containing the triangle i [Fig.4(a)].This is a natural choice of a composite order parameter because it preserves the threefold rotational symmetry.The correlation function of CT ≡ ⟨T i T j ⟩ shows a structure-factor peak at the M point of momentum space in the VP phase, as shown in Fig. 4(b).The corresponding real-space CT correlation function-plotted in Fig. S3(a) and sketched in the middle subfigure of Fig. 1(b)-clearly reveals the pattern of translational symmetry breaking and the associated 2×2 unit cell.We can also try to identify CT q=M with operators of the low-energy effective action (4).The operator should break translational symmetry while preserving the sixfold rotational symmetry.Since CT q=M is constructed from dimers in (3), it should also be gauge invariant.The natural candidate is therefore which is clearly invariant under the sixfold rotation operation defined in Eq. (S9) of Supplementary Note 1.After this identification, we can infer the critical behavior of CT q=M near the second-order phase transition at V ∼ 0.6 from the knowledge of the scaling dimension of X of the cubic conformal field theory (CFT).As mentioned previously, the cubic CFT and the O(3) CFT have similar critical exponents; since X belongs to the j = 2 representation of O(3) (the T representation in the notation of Ref. [61]), we obtain the critical exponent η * ≈ 1.42 [57,61,62], which is close to the critical exponent for the (2+1)D O(2) * transition observed before in a frustrated kagome magnet [14,32].
The cubic * phase transition From the histograms in Fig. 3 (e), (f) and (g), we note that there is a sizeable cubic anisotropy order parameter even close to the phase transition point between the VP and QSL phases.Theoretically, the cubic anisotropy was thought to be irrelevant at the 3D O(3) Wilson-Fisher fixed point [56][57][58][59][60].However, recent conformal bootstrap analyses [61,62] demonstrate that the 3D O(3) fixed point is unstable and will flow towards the nearby cubic fixed point, whose symmetry-breaking phase is the corner cubic phase.Thus, the phase transition point here is indeed of the cubic * criticality, with the * representing that the order parameter here is constructed by a fractional quasiparticle.
We now turn to the critical behavior of the VP-QSL transition shown in Fig. 5. From the plot of the vison order parameter |ϕ| as a function of V in Fig. 2 for different system sizes, the transition around V = 0.6 is clearly continuous, and-per the discussion abovebelongs to the cubic * CFT.However, the critical exponents of the cubic fixed point and those of the O(3) fixed point are not practically distinguishable with the resolution in our study (the scaling dimension of the cubic anisotropy differs by only ∼ 0.01 [61]).For example, the most recent classical Monte Carlo simulations give −0.00001 < ν O(3) − ν cubic < 0.00007 and    2), β = 0.33 (5), and ν = 0.75 (8).The values of ν and β are consistent with the standard O(3) value, ν = 0.7112 (5) and β = 0.3689 (3).(e) and (f) show the data collapse of the vison order parameter and its Binder ratio B2.Here, the values of Vc, β, and ν are determined independently from our fitting procedures.The error bars here represent the standard error of the mean, which is calculated as σ/ √ N , where σ is one standard deviation and N is the total number of independent samples.
η O(3) − η cubic = −0.00061(10)[76], confirming an earlier prediction of a perturbative calculation [77].The conformal bootstrap result in Ref. [61] also implies that the difference should be ∼ 10 −4 .Thus, the logic is that if we find a cubic anisotropy of the order parameter near the phase transition and the critical exponents are the same as the O(3) criticality up to QMC precision, then this amounts to the cubic phase transition being demonstrated numerically.Therefore, in Fig. 5, we perform stochastic data collapse with both the vison order parameter ϕ and its Binder cumulant [78,79] as shown in Figs.5(a) and  (b); the details can be found in Supplementary Note 7. By including more data points and a larger parameter space for the optimization process, our scheme, as shown in Figs.5(b) and (e), gives independent estimates of V c = 0.59(2), β = 0.33 (5), and ν = 0.75 (8), which are consistent with the O(3) critical exponents of ν = 0.7112 (5) and β = 0.3689 (3).Furthermore, in Figs.5(e) and (f), ϕ and its Binder ratio are illustrated with a rescaled x-axis to show the high-quality data collapse.Our numerical analyses firmly establish that the QLM in (1) realizes different limits of the effective action ((4) in the next section) and are consistent with the latest understanding of the fixed points of cubic symmetry.

Effective action and renormalization-group analysis
In order to fully understand the LN-VP first-order transition and the VP-QSL continuous transition, as well as their fundamental relation to the O(3) and cubic fixed points in (2 + 1)D, we will need to first explain the effective action and the renormalization-group (RG) analysis of the problem in a more general and up-to-date setting.The low-energy description of these transitions [33], with the O(3) order parameter ϕ in (2), can be cast into the action S = dtd 2 x L with the Lagrangian coupled to a Z 2 gauge theory.The . . .denote higherorder terms, whose explicit forms can be found in the SI.The Z 2 gauge theory does not change the dynamics of the system; its effect is to gauge out operators that are odd under the Z 2 symmetry which flips the sign of all three ϕ i fields.The first three terms in (4) preserve the O(3) symmetry.The ν 4 piece (ignoring other higher-order terms), on the other hand, breaks the O(3) symmetry to a threedimensional cubic symmetry.This cubic model has a long history since its first appearance in describing the structural phase transition of perovskites [56,59] (for a review, see Ref. [60]).More recently, the experimentally observed structural phase transitions of single-layer transition metal dichalcogenides [80][81][82], such as MoS 2 , have also been described by the effective action above [83].
Mean-field theory suggests that there are two phases: when ν 4 > 0, in the symmetry-broken phase (r < 0), the minima of the effective potential are located at Such a phase is commonly called the corner-cubic phase (i.e., the VP phase in our case), yielding the bright spots in the histograms in Figs.3(d)-(f); here, the eight possible states are in one-to-one correspondence with the eight vertices of a three-dimensional cube.When ν 4 < 0, the symmetry-broken phase is described by together with the other states where either ⟨ϕ 2 ⟩ or ⟨ϕ 3 ⟩ acquires an expectation value.The six symmetry-broken states in this case correspond to the face centers of a three-dimensional cube, thus defining a face-cubic phase (i.e., the LN phase in our case which corresponds to the bright spots in the histograms in Figs.3(a) and (b)], which has also been seen in the honeycomb-lattice transverse-field Ising model [52] that bears a low-energy action similar to (4).
To further demonstrate the first-order phase transition here, with the above background, we now move back to the QMC data in Figs. 3 and 6.It is important to notice that at the cubic fixed point, the coupling constant ν 4 is a small positive number.This leaves room for another weakly first-order phase transition when ν 4 changes sign.This is precisely the transition between the corner-and face-cubic phases that we observed at V /t ≈ 0.05 in the vison order-parameter histograms in Fig. 3.
To quantitatively describe the phase transition between the LN and VP phases, we consider the anisotropy parameter ν 4 in Eq.( 4).In the Monte Carlo simulations, it can be expressed as [52] (see derivations in Supplementary Note 5) where Y 0 4 and Y 0 6 are two spherical harmonics given by Y is defined in (2) for every component.The volume factor is vol = L 2 β.Notice the leading term in the above equation essentially measures the angular dependence of the expectation of ⟨ϕ(θ, ψ)⟩, projecting onto Y 0 4 , with the signs of this term different in the corner-and face-cubic phases.This term gives the leading contribution to ν 4 .As shown in Figs. 6, the sign of the anisotropy parameter ν 4 changes from negative to positive near the LN-VP transition.A similar approximate emergent continuous symmetry at a first-order transition has been reported in the checkerboard J-Q model in both 2D and 3D lattices [84,85]; the latter is in fact related to ongoing experimental efforts in understanding the thermodynamic data of the Shastry-Sutherland material SrCu 2 (BO 3 ) 2 under high pressure [86][87][88].
This analysis also tells us that the phase transition from the disordered phase to the corner-cubic phase is second-order and continuous because the RG flow suggests that when ν 4 > 0, the theory (4) flows to a conformal field theory.It remains to answer the question to which universality class the ν 4 > 0 phase transition belongs.This seemingly easy question was under debate for many years; see Ref. [59] for a review of early theoretical studies.It was not until much later that a consistent answer was obtained using three different methods: lattice Monte Carlo simulations [57], perturbative field-theory calculations [58], and the conformal bootstrap [61].It turns out that the continuous phase transition at ν 4 > 0 is in a new universality class, which is different from the O(3)-invariant Heisenberg model.The question of the universality class is closely related to the scaling dimension of the cubic anisotropy operator O = 3 i=1 (ϕ i ) 4 .Among the above-mentioned methods, the conformal bootstrap approach [89,90] provides a con-crete nonperturbative proof that ∆ O < 3 [61,62].This indicates that the O(3) conformal fixed point is unstable against cubic anisotropy, and therefore, the phase transition should be in a new universality class.
We denote the corresponding CFT the cubic fixed point.Since, in our case, the symmetric phase is a Z 2 QSL with fractionalized vison excitations, the VP-QSL transition is actually of the cubic * universality, with the same critical exponents β and ν as the cubic CFT, but it acquires a large anomalous dimension η * since the order parameter is actually made of a composite object of the fractionalized particles.Such * -type transitions between a symmetry-breaking phase and a Z 2 QSL have been shown in kagome-lattice frustrated magnetic models with (2+1)D O(2) * universality and η * ∼ 1.5 [14,32,35].

Relation to experiments on Rydberg atom arrays
Recently, Ref. [20] demonstrated that the phases of various triangular-lattice quantum dimer models can be realized using Rydberg atoms arrayed on the sites of a kagome lattice via programmable optical tweezers.In particular, this presents a direct route to the experimental realization of the fully packed quantum loop model studied in this paper in quantum simulators based on interacting Rydberg atoms.In such a setup, strong van der Waals interactions prevent the simultaneous excitation of neighboring atoms to the Rydberg state in a phenomenon known as the Rydberg blockade.This Rydberg blockade competes with the chemical potential (the laser detuning) which favors or disfavors occupation of the Rydberg states.
Such a blockaded kagome-lattice Rydberg array can be mapped to a dimer model on the triangular lattice with two dimers per site [20] in certain parameter regimes.Numerical studies [20,28] have in fact suggested that both the LN and QSL phases of the triangular-lattice QLM may arise in this setup.Therefore, it would be interesting to investigate the possible existence of the hidden-order VP phase discovered in this study in the experimental atomic system.
Importantly, this hidden order also bears significant implications for the experimental identification of Z 2 QSLs in Rydberg atom arrays.In recent works [21,22,28], a diagonal loop operator has been used to detect the Z 2 QSL phase in a model of Rydberg atoms with emergent gauge-charged Ising matter degrees of freedom [25].This operator is defined as an arbitrary loop that runs across several bond centers, i.e., Z ≡ (−1) c , where c counts the number of dimers cut by the loop.While our model has no Ising matter, such a loop operator will have a nonzero expectation value in the VP phase as well, even though it has no topological order but only a hidden translation symmetry breaking.Thus, the hidden VP phase actually reveals that the nonzero diagonal loop operator plus a disordered diagonal density is not a sufficient condition for a QSL.Discussion In this work, using the newly developed sweeping cluster QMC algorithm, supplemented with ED and symmetry analysis of the vison order parameter, we mapped out the entire phase diagram of the triangular-lattice QLM in an unbiased manner.Besides the lattice nematic and Z 2 quantum spin liquid phases, we discover a hidden vison plaquette phase via off-diagonal measurements in the dimer basis.This VP phase is invisible to diagonal measurements and obeys the local constraint of a Z 2 gauge field.Our results reveal that the LN-VP firstorder transition is triggered by the change from face-to corner-cubic anisotropy of the O(3) vison order parameter, and the VP-QSL continuous transition, driven by the condensation of visons, is of the cubic * universality class-which is very close to the (2 + 1)D O(3) one-in full consistency with recent conformal bootstrap findings on the cubic fixed point [61,62].
The VP solid phase discovered here exhibits clear translational symmetry breaking in the vison correlation functions and in the dimer hopping order in Fig. 4 but has no apparent dimer density order in our ED and QMC simulations.That is why it had been treated as a part of the QSL phase for a long time.It represents a hidden state of quantum matter and, at the same time, resolves the previous controversy regarding the phase boundary between the QSL and LN phases in Refs.[33,34] (note that Ref. [34] had observed the vanishing dimer order parameter characteristic of this phase in a parameter regime consistent with ours).
In light of recent experiments with programmable quantum simulators based on highly tunable Rydberg atom arrays, our results could also have direct experimental relevance.The experiment in Ref. [22] probes the case where the atoms are positioned on the links of the kagome lattice, connecting to the quantum dimer model on the kagome lattice [21].Our study on the triangularlattice quantum loop model connects to the case where the atoms are placed on the sites of the kagome lattice [20].Investigation of the LN, hidden-order VP, and Z 2 QSL phases-as well as their subtle phase transitions in the context of the O(3) and cubic fixed points-can inspire new experiments.In particular, for the case with Rydberg atoms on the kagome sites, an analog of the VP order is a promising possibility for the region proximate to an even QSL [25] or the so-called string phase identified in Ref. [20].Finally, the rich interplay between the VP and QSL phases is a promising avenue for related future numerical and theoretical studies of this system [25,28].

Method
Sweeping cluster algorithm.This is a quantum Monte Carlo method developed by the authors which can work well in constrained spin models [28,36,47,70,71].The key idea of the sweeping cluster algorithm is to sweep and update layer by layer along the imaginary time direction, so that the local constraints (gauge field) are recorded by update-lines.In this way, all the samplings are done in the restricted Hilbert space, i.e., the low-energy space.In this article, we can measure the information about single visons because in a strictly constrained space, the energy gap of other quasiparticles, such as spinons, becomes infinitely large and thus these quasiparticles do not exist in the restricted Hilbert space.Recently, this method has even been developed for simulating systems with soft constraints [28] or higher-order constraints [47].
These matrices generate a finite subgroup of O(3) which contain 48 elements [33].This group is isomorphic to Z 2 × S 4 , with S 4 being the permutation group of four elements.

SUPPLEMENTARY NOTE 2: ORDER PARAMETER HISTOGRAMS
Besides the order parameter histograms presented in the main text, here, we add more detailed results on the evolution from the LN to the VP solid phase.Figure S2 compiles the histograms for a L = β = 16 system, as V is varied from −1 to 1.It is clear that for cases with V < 0.1, the histograms exhibit the face-cubic distribution, and for the cases with V > 0.1, the histograms exhibit the corner-cubic distribution.
Close to the LN-to-VP transition, i.e., the panel with V = 0.1 in Fig. S2(b), the histogram develops asymptotic O(3) symmetry with very weak corner-cubic anisotropy (this is because V = 0.1 lies on the VP side of the transition).Inside the VP phase, and as V increases towards the VP-QSL continuous transition point V c = 0.598 (see below), the histograms [in panels with V = 0.3, 0.5 in Fig. S2(b)] begin to shrink in diameter, suggesting the reduction of the vison order parameter.Eventually, inside the QSL phase [panels with V = 0.7, 0.9 in Fig. S2(b)], the histogram becomes a point at (0, 0, 0), indicating that the vison order parameter completely vanishes in the topologically ordered QSL phase.

SUPPLEMENTARY NOTE 3: REAL-SPACE VISON CORRELATIONS
To investigate the symmetry-breaking pattern of the VP phase, we calculate the correlations of the operator T described in Eq. ( 9) of the main text.This operator is defined on a triangular plaquette from three related kinetic (t) terms, T i = t 1,i + t 2,i + t 3,i , [see inset of Fig. S3(a)].Here, i labels the triangle, and 1, 2, 3 denote the three kinds of kinetic (∼ t) terms (rhombi) in the Hamiltonian acting on this triangle.The correlation of the operators between the triangular plaquettes i and j is ⟨T i T j ⟩.
We calculate the real-space vison correlations, selecting the leftmost and lowest red triangle below as a reference (its correlation is identically 1), as illustrated in Fig. S3(b).Two visons living on the centers of related triangular plaquettes i and j can be connected with a open string; their correlation is ⟨v i v j ⟩ = ⟨(−1) N P ij ⟩ and N Pij is the number of dimers cut along the path P between plaquettes i and j.In the two-dimer-per-site case, ⟨v i v j ⟩ is independent on the chosen path, i.e., the visons here are gauge invariant and measurable.
We found that the T -operator and vison correlations clearly exhibit the translational symmetry breaking of the VP solid, but preserve C 6 rotational symmetry.Although the vison correlation function appears to change sign under mirror-reflection and sixfold rotations, the physical state actually preserves these symmetries, because of the two-toone correspondence between vison and dimer configurations.Individual triangles with positive (negative) values of vison correlations gather into big triangles.In addition, the central triangle of the thus-formed big triangle always has a larger absolute value, which is about three times that of the other triangles in the same big triangle.The translation periods along the s 1 and s 2 directions are both 2, which means that the unit cell is a 2 × 2 rhomboid.The sweeping cluster QMC approach employed in this work is a new method developed by us, which works well in constrained quantum lattice models [28,36,70,71].Prior to the development of sweeping cluster QMC, in order to solve the QDM or QLM types of constrained models, one had to rely on either exact diagonalization of small systems, variational approaches such as DMRG that suffer from finite-size effects on cylindrical geometries [33], projector Monte Carlo approaches (which include the Green's function [34,41,42,44,91,92] and diffusion Monte Carlo schemes [93,94]), or sampling directly in height space and discarding the unconstrained configurations [95][96][97].These projector Monte Carlo methods obey the geometric constraints but are not efficient away from the RK point [98] and only work at T = 0. Furthermore, there does not exist any cluster update schemes for the projector methods.On the contrary, the sweeping cluster algorithm, based on the world-line Monte Carlo scheme [99][100][101], is designed to sweep and update layer-by-layer along the imaginary-time direction so that the local constraints (gauge fields) are recorded by the update lines.In this way, all the samplings are performed in the restricted Hilbert space and it provides a cluster update scheme for constrained systems [71] that works at all temperatures.Proper finite-size-scaling analyses can then be carried out to explore phase transitions and critical phenomena.
To implement the constraint of having two dimer per site in our QMC simulation, we set the initial state as one of the three LN patters with the same probability [34].In the region of the LN phase, the equilibrium process will not change the orientation of the initial state; therefore, the randomly selected initial state ensures each pattern of the LN phase has an equal probability of 1  3 to emerge.In the VP and QSL regions, one can also use such initialization schemes, and we have tested that in these regions too, they produce the same QMC results as a completely random initialization.We simulate L × L triangular lattices with system sizes L = 8, 12, 16 with the inverse temperature set to β = L and 10 4 Monte Carlo samplings were used to obtain average values of the observables in all calculations.
In this study, we can measure the information about single visons because in a strictly constrained space, the energy gap of other quasiparticles such as spinons, becomes infinitely large and thus these quasiparticles do not exist in the restricted Hilbert space; the vison is thus well-defined here.
Through some simple algebra, the three anisotropy parameters can be expressed as We have also used the approximation ⟨ϕ n ⟩ 0 = ⟨ϕ n ⟩, which is valid since their difference is of linear order in the perturbation parameters ν 4 , ν 6 , or ν ′ 6 .These differences only introduce higher-order corrections in Eq. (S11).

SUPPLEMENTARY NOTE 6: DIMER DENSITY AND CORRELATION FUNCTIONS
In this section, we apply the Lanczos exact diagonalization (ED) method to a 4 × 4 (N bond = 48) lattice within the VP phase.The number of configurations in the constrained Hilbert space is 586, 695.Using the ground-state wavefunction, we analyze the configurations symmetry of eight points of cubic order parameters without statistical error as the following explanation.As shown in Fig. S2, the eight points lie in the eight octants, so each configuration of the wavefunction can be classified into a certain octant, except the ones on the boundaries between different octants.Then-akin to what is done in QMC-we classify configurations according to related octants and average over the ones in the same class to obtain the dimer density of a certain cubic order parameter point.It is worth noting that two classes of opposite octants have opposite signs for vison configurations but their dimer configurations are the same.Thus, we have to average the dimer configurations of the two opposite octants.The result [Fig.S4(a)] is similar to that obtained with QMC: there is still no dimer order to within the numerical precision of ED.We found that for each link, the real-space dimer density is 1/3.
Similarly, we average the ED vison configurations in each octant.The histogram of the order parameter and the real-space vison density obtained via ED are consistent with our QMC results.The real-space vison density is shown in Fig. S4(b).All these results strongly support the QMC conclusions that there is indeed a hidden-order vison plaquette phase between the spin liquid and nematic phases.
We also calculate the structure factor of the dimer-dimer correlation functions at the k = (0, 0) point as shown in Fig. S5.The stucture factor is defined as where ⟨D j D k ⟩ = [(⟨d 1 j d 1 k ⟩+⟨d 2 j d 2 k ⟩+⟨d 3 j d 3 k ⟩]/3 is the average dimer correlation functions of three orientation of rhombi.The dimer order parameter S(0, 0) is large in the LN phase, which preserves the transnational symmetry, and goes to zero at the first-order transition point between LN and VP phase, but it shows no sign for the continuous transition of the VP and QSL phase.we use a polynomial fit to the third order for the raw vison order parameter and Binder ratio data, and the χ 2 /d.o.f (χ 2 per degree of freedom) for each fitting curve is close to 1, where χ 2 = N i=1 (y ′ i − y i ) 2 /σ 2 i , y ′ i is the data point on the fit curve, and y i and σ i are the raw data point and the error bar, respectively.For the Binder cumulant, we observe that the crossing between the curves for system sizes L and 2L progressively shifts to the right because of the finite-size effect.Then, we set V c , β, and ν as three free parameters to carry out the sampling for data collapse both for the vison order parameter and its Binder cumulant, and find the best fitting curve in the scaled x = (V − V c )L 1/ν and x ′ = (V − V c )/V c L 1/ν with χ 2 /d.o.f ∼ 1.The three parameters (V c , β, ν) are all selected randomly within an initial region according to our practical estimation.The converged distribution of the parameters (V c , β, ν) is shown in Figs.5(b) and (e) in the main text, which compiles around 1000 independent polynomial fit results with χ 2 /d.o.f close to 1.To reduce the influence of the small system sizes, we only use L ≥ 8 for the scaling.In this fashion, we determine V c = 0.59(2), β = 0.33 (5) and ν = 0.75 (8); the obtained β and ν are consistent with the O(3) value of β = 0.3689(3) and ν = 0.7112(5) [104] within error bars.We then collapse the vison order parameter and its Binder ratio data with the obtained V c , β, and ν and the results are shown in Figs.5(c) and (f) in the main text.

1 FIG. 1 .
FIG. 1. Fully packed quantum loop model on the triangular lattice.(a) Schematic representation of the QLM; s1 and s2 are the primitive vectors.The t and V terms in the Hamiltonian (1) are depicted in the upper-right insets.The dimer configuration sketched is one of the three LN patterns, with fully packed loops along the s1 direction.(b) Phase diagram of the QLM obtained from our simulations.The first row of the left subfigure illustrates the three LN dimer configurations, corresponding to the three vison patterns shown on the second row.The triangles v1 and v2 represent two sublattices for the visons, and the red and grey colors in the triangles denote vison numbers (±1) with opposite sign.The first-order phase transition between the LN and VP states occurs at V = 0.05(5).The first row of the middle subfigure is the kinetic energy correlation pattern and the second row is the vison plaquette (VP) pattern based on QMC simulation results, respectively (see Fig.S3(a,b) in Supplementary Note 3 for values in each triangle).VP is a hidden dimer solid phase without dimer order.The right subfigures illustrate the representative dimer coverings in the Z2 QSL phase.The continuous phase transition between the VP phase and the QSL occurs at Vc = 0.59(2).

FIG. 2 .
FIG. 2.Vison O(3) order parameter |ϕ| as a function of V .The two dashed lines are guides to the eye for the position of the LN-VP first-order transition and the VP-QSL continuous transition.The error bars here represent the standard error of the mean, which is calculated as σ/ √ N , where σ is one standard deviation and N is the total number of independent samples.

FIG. 3 .
FIG. 3. Histograms of the O(3) order parameter.The histograms are plotted on the two-dimensional (ϕ1,ϕ2) plane and are obtained from the QMC data for L = 24 systems with different V .The (a,b) face-cubic anisotropies are observed inside the LN phase, while the histogram in a region of phase coexistence, plotted in (c), indicates the first-order nature of the transition.The (d, e, f) corner-cubic anisotropies are observed inside the VP phase, and (g) near the continuous phase transition point.

30 FIG. 4 .
FIG.4.Hidden order in the VP phase.(a) The definition of the operator Ti acting on triangle i.Each operator Ti includes the t-terms of the Hamiltonian (1) on the three rhombi with labels 1, 2, 3 enclosing the triangle i.(b) The structure factor at the M point, CTq=M , as a function of V ; in the VP phase, CTq=M acquires long-range order.The inset shows that the VP phase (for V = 0.1, 0.15, and 0.3, the extrapolated value at L = ∞ is 0.006(2), 0.008(2), and 0.0097(3), respectively) persists in the thermodynamic limit.The error bars here represent the standard error of the mean, which is calculated as σ/ √ N , where σ is one standard deviation and N is the total number of independent samples.

FIG. 5 .
FIG.5.The cubic fixed point.Stochastic data collapse with both the vison order parameter |ϕ| and its Binder cumulant.(a) and (b) plot the vison order parameter and its Binder cumulant U2.All data points are connected with polynomial fits to the third order and χ 2 /d.o.f (χ 2 per degree of freedom) is close to 1.The location of the crossing in (b) shifts due to finite-size effects.(c) and (d) are the distributions of β and ν with Vc obtained from the bootstrap sampling process.Each data point is determined from a fitted curve of the data collapse for |ϕ| and B2, and the the χ 2 /d.o.f of all the curves is close to 1.We obtain Vc = 0.59(2), β = 0.33(5), and ν = 0.75(8).The values of ν and β are consistent with the standard O(3) value, ν = 0.7112(5) and β = 0.3689(3).(e) and (f) show the data collapse of the vison order parameter and its Binder ratio B2.Here, the values of Vc, β, and ν are determined independently from our fitting procedures.The error bars here represent the standard error of the mean, which is calculated as σ/ √ N , where σ is one standard deviation and N is the total number of independent samples.

0 4 = 3 ( 3 − 16 FIG. 6 .
FIG.6.Evidence for the first-order phase transition.The figure shows the coefficient of the anisotropy ν4 in the effective action (4) as function of V , which is computed from the QMC histogram data of the O(3) order parameter in Fig.3.ν4 changes sign at the LN-VP transition, signifying the change from face-to corner-cubic anisotropy in the effective action of the lattice model.The error bars here represent the standard error of the mean, which is calculated as σ/ √ N , where σ is one standard deviation and N is the total number of independent samples.
FIG. S3. (a)Real-space T -operator correlations of the L = 8 system at V = 0.3.The inset shows the summation of the three t terms in the Hamiltonian.(b) Real-space vison correlations of the L = 8 system at V = 0.3; s1 and s2 are the primitive vectors.For both figures, we set the leftmost and lowest triangle as the reference.The red (grey) color conveys the positive (negative) value of the correlation.
FIG. S4. (a)Exact-diagonalization result for the real-space dimer density for the L = 4 system at V = 0.3; as previously, s1 and s2 are the primitive vectors.All bonds hold 1/3 dimer density.(b) The real-space vison density, where we have set v1(0, 0) = 1.Similar to the QMC results, four small triangles with the same sign of the vison density form bigger triangles, and the small triangles at the centers of the big ones have a larger absolute value of the density.