Introduction

In quantum many-body systems, topologically ordered phases1,2,3,4 and their fractionalized excitations2,5 have attracted much attention because they hold the promise of providing solutions to many mysteries of highly entangled quantum matter as encountered in unconventional superconductors6,7,8, quantum moiré materials9,10, frustrated quantum magnets11,12,13,14,15,16,17,18,19, and—more recently—programmable quantum simulators20,21,22,23,24,25,26,27,28,29,30. Concomitantly, the quantum phase transitions between topologically ordered phases and conventional symmetry-breaking states have also been actively investigated with a variety of theoretical and numerical approaches14,17,25,28,31,32,33,34,35,36.

One of the most widely studied models to realize such phase transitions, in either numerical simulations or cold-atom experiments, is the quantum dimer model (QDM)6,8,37. In QDMs, \({{\mathbb{Z}}}_{2}\) topologically ordered phases originate from the famous Rokhsar-Kivelson (RK) point8,38,39,40 and acquire extended phase space on nonbipartite triangular and kagome lattices36,39,41,42,43. However, away from the RK point, the overall structures of the phase diagrams on various 2D geometries—either bipartite or nonbipartite—such as square and triangular lattices, are still under intensive investigation24,33,34,36,44,45,46,47. This is due to the lack of controlled methodology to solve these strongly correlated lattice models with local constraints, and discoveries that either confirm or defy earlier nonrigorous arguments.

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 constraints33,34,52,53. The QLM is also related to the resonant-valence-loop phase in frustrated spin-1 models54,55, in analogy to the relation between the one-dimer-per-site QDM and spin-\(\frac{1}{2}\) models. The triangular-lattice QLM has been investigated by field-theoretical analyses and various numerical approaches33,34, but, so far, no consistent phase diagram has been obtained.

At the RK point, the model realizes a \({{\mathbb{Z}}}_{2}\) quantum spin liquid (QSL) phase, which has even \({{\mathbb{Z}}}_{2}\) topological order and hosts fractionalized quasiparticles called visons. Tuning away from the RK point, the model also 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 studies33,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 \({{\mathbb{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) calculations33, or variational quantum Monte Carlo (QMC) simulations with intermediate system sizes of 12 × 1234. 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 point56,57,58,59,60. However, this assumption has been disproven by recent conformal bootstrap analyses61,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 study52 on the honeycomb-lattice 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 triangular-lattice QLM requires an additional and consistent understanding that unifies the aforementioned recent developments33,34,52,61, and provides the correct mechanism of the vison-condensation transition from the \({{\mathbb{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 phase63 and the \({{\mathbb{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 \({{\mathbb{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 points56,57,58,59,60 in (2 + 1)D61. 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 arrays20,21,28,64,65,66,67,68,69.

Results

The constrained model

The Hamiltonian of the QLM on a triangular lattice is defined as

(1)

where α represents the three possible orientations of all plaquettes of the triangular lattice, as shown in Fig. 1a for an example of the LN state. The local constraint of our model requires two dimers to touch every triangular-lattice site in any configuration, thereby forming the fully packed quantum loops34. 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 \({{\mathbb{Z}}}_{2}\) QSL solution39. We set t = 1 as the unit of energy in our simulations.

Fig. 1: Fully packed quantum loop model on the triangular lattice.
figure 1

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. S3a, 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 \({{\mathbb{Z}}}_{2}\) QSL phase. The continuous phase transition between the VP phase and the QSL occurs at Vc = 0.59(2).

Full size image

We employ the recently developed powerful sweeping cluster quantum Monte Carlo method28,36,47,70,71 to solve this model. Our simulations are performed on the triangular lattice with periodic boundary conditions and system sizes N = 3L2 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 modes33,72,73,74, which is given by

$${\phi }_{j}=\sum\limits_{{{{{{{{\bf{r}}}}}}}}}({v}_{1,{{{{{{{\bf{r}}}}}}}}},{v}_{2,{{{{{{{\bf{r}}}}}}}}})\cdot {{{{{{{{\bf{u}}}}}}}}}_{j}{e}^{i{{{{{{{{\bf{M}}}}}}}}}_{j}\cdot {{{{{{{\bf{r}}}}}}}}},\quad j=1,2,3,$$
(2)

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 vivj  =  \({(-1)}^{{N}_{{P}_{ij}}}\), with \({N}_{{P}_{ij}}\) being the number of dimers cut along the path P between triangular plaquettes i and j36. The label (v1,r, v2,r) denotes visons living at the centers of the two sublattices of the triangular plaquette at r [see the left subfigure in Fig. 1b]. In our simulations, we choose a reference vison on a certain triangular plaquette, e.g., by setting v1,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 (v1,r, v2,r) for all r. The three momenta M1  =  \((\pi /\sqrt{3},\pi /3)\), M2  =  \((\pi /\sqrt{3},-\pi /3)\), and M3 = (0, 2π/3) correspond to the low-energy modes of the vison dispersion and are invariant under the projective symmetry transformations33,74 as we show in Supplementary Note 1. The associated uj are (1, 1)T for M1 and M2, and (1, −1)T for M3. 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 Supplementary Note 1.

According to the analysis of the effective critical theory33,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 \({\mathbb{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 \({\mathbb{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 \(| {{{{{{{\boldsymbol{\phi }}}}}}}}| =\sqrt{{\phi }_{1}^{2}+{\phi }_{2}^{2}+{\phi }_{3}^{2}}\), plotted in Fig. 2, as the order parameter to roughly detect the solid phases and their transitions to the \({{\mathbb{Z}}}_{2}\) QSL phase. As illustrated in Fig. 2, the order parameter ∣ϕ∣ clearly vanishes via a two-step 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). We thus ask: why does the nature of the order seem so different when viewed in terms of dimers and visons?

Fig. 2: Vison O(3) order parameter ∣ϕ∣ as a function of V.
figure 2

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 \(\sigma /\sqrt{N}\), where σ is one standard deviation and N is the total number of independent samples.

Full size image

To further reveal the nature of the symmetry-breaking 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. 3c, 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.

Fig. 3: Histograms of the O(3) order parameter.
figure 3

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–f corner-cubic anisotropies are observed inside the VP phase, and g near the continuous phase transition point.

Full size image

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 decade33,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. S3b 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 triangles75. 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 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. 3d–f is indeed reflected in other correlation functions. Firstly, this can be seen from the vison correlations 〈vivj〉, 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:

$${T}_{i}={t}_{1,i}+{t}_{2,i}+{t}_{3,i},$$
(3)

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. 4a]. This is a natural choice of a composite order parameter because it preserves the threefold rotational symmetry. The correlation function of CT  ≡  〈TiTj〉 shows a structure-factor peak at the M point of momentum space in the VP phase, as shown in Fig. 4b. The corresponding real-space CT correlation function—plotted in Fig. S3a and sketched in the middle subfigure of Fig. 1b—clearly reveals the pattern of translational symmetry breaking and the associated 2 × 2 unit cell.

Fig. 4: Hidden order in the VP phase.
figure 4

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 \(\sigma /\sqrt{N}\), where σ is one standard deviation and N is the total number of independent samples.

Full size image

We can also try to identify CTq=M with operators of the low-energy effective action (4). The operator should break translational symmetry while preserving the sixfold rotational symmetry. Since CTq=M is constructed from dimers in (3), it should also be gauge invariant. The natural candidate is therefore X = ϕ1ϕ2 + ϕ2ϕ3 − ϕ3ϕ1, 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 CTq=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.4257,61,62, which is close to the critical exponent for the (2+1)D O(2)* transition observed before in a frustrated kagome magnet14,32.

The cubic* phase transition

From the histograms in Fig. 3e–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 point56,57,58,59,60. However, recent conformal bootstrap analyses61,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 above—belongs 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.0161). For example, the most recent classical Monte Carlo simulations give − 0.00001 < νO(3) − νcubic < 0.00007 and ηO(3) − ηcubic = − 0.00061(10)76, confirming an earlier prediction of a perturbative calculation77. 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.

Fig. 5: The cubic fixed point.
figure 5

Stochastic data collapse with both the vison order parameter ∣ϕ∣ and its Binder cumulant. a, 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, 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, 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 \(\sigma /\sqrt{N}\), where σ is one standard deviation and N is the total number of independent samples.

Full size image

Therefore, in Fig. 5, we perform stochastic data collapse with both the vison order parameter ϕ and its Binder cumulant78,79 as shown in Fig. 5a, 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 Fig. 5b, e, gives independent estimates of Vc = 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 Fig. 5e, 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 transitions33, with the O(3) order parameter ϕ in (2), can be cast into the action S = ∫dtd2x L with the Lagrangian

$$L=\sum\limits_{i=1}^{3}{({\partial }_{\mu }{\phi }_{i})}^{2}+r\sum\limits_{i=1}^{3}{\phi }_{i}^{2}+\mu {\left(\sum\limits_{i = 1}^{3}{\phi }_{i}{\phi }_{i}\right)}^{2}+{\nu }_{4}\sum\limits_{i=1}^{3}{\phi }_{i}^{4}+\ldots ,$$
(4)

coupled to a \({{\mathbb{Z}}}_{2}\) gauge theory. The … denote higher-order terms, whose explicit forms can be found in the SI. The \({{\mathbb{Z}}}_{2}\) gauge theory does not change the dynamics of the system; its effect is to gauge out operators that are odd under the \({{\mathbb{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 three-dimensional cubic symmetry. This cubic model has a long history since its first appearance in describing the structural phase transition of perovskites56,59 (for a review, see ref. 60). More recently, the experimentally observed structural phase transitions of single-layer transition metal dichalcogenides80,81,82, such as MoS2, have also been described by the effective action above83.

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

$$\langle {\phi }_{1}\rangle =\pm v,\quad \langle {\phi }_{2}\rangle =\pm v,\quad \langle {\phi }_{3}\rangle =\pm v.$$
(5)

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 Fig. 3d–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

$$\langle {\phi }_{1}\rangle =\pm v^{\prime} ,\quad \langle {\phi }_{2}\rangle =\langle {\phi }_{3}\rangle =0,$$
(6)

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 Fig. 3a, b), which has also been seen in the honeycomb-lattice transverse-field Ising model52 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.

Fig. 6: Evidence for the first-order phase transition.
figure 6

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 \(\sigma /\sqrt{N}\), where σ is one standard deviation and N is the total number of independent samples.

Full size image

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 as52 (see derivations in Supplementary Note 5)

$$\begin{array}{rc}{\nu }_{4}&=-\frac{1}{vol}\frac{15\sqrt{\pi }}{\langle {\phi }^{4}\rangle }\left(\langle {Y}_{4}^{0}\rangle +\frac{\langle {\phi }^{2}{Y}_{4}^{0}\rangle \langle {\phi }^{4}\rangle \langle {\phi }^{6}\rangle -\langle {Y}_{4}^{0}\rangle {\langle {\phi }^{6}\rangle }^{2}}{{\langle {\phi }^{6}\rangle }^{2}-\langle {\phi }^{4}\rangle \langle {\phi }^{8}\rangle }\right),\end{array}$$
(7)

where \({Y}_{4}^{0}\) and \({Y}_{6}^{0}\) are two spherical harmonics given by \({Y}_{4}^{0}=3(3-30{\cos }^{2}\theta +35{\cos }^{4}\theta )/(16\sqrt{\pi })\), \({Y}_{6}^{0}=\sqrt{13}(-5+105{\cos }^{2}\theta -315{\cos }^{4}\theta +231{\cos }^{6}\theta )/(32\sqrt{\pi })\), and θ is the polar angle of the order parameter ϕ = (ϕ1, ϕ2, ϕ3) which is defined in (2) for every component. The volume factor is vol = L2β. Notice the leading term in the above equation essentially measures the angular dependence of the expectation of 〈ϕ(θ, ψ)〉, projecting onto \({Y}_{4}^{0}\), 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 Fig. 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 lattices84,85; the latter is in fact related to ongoing experimental efforts in understanding the thermodynamic data of the Shastry-Sutherland material SrCu2(BO3)2 under high pressure86,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 simulations57, perturbative field-theory calculations58, and the conformal bootstrap61. 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 \({{{{{{{\mathcal{O}}}}}}}}=\mathop{\sum}_{i = 1}^{3}{({\phi }_{i})}^{4}\). Among the above-mentioned methods, the conformal bootstrap approach89,90 provides a concrete nonperturbative proof that \({{{\Delta }}}_{{{{{{{{\mathcal{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 \({{\mathbb{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 \({{\mathbb{Z}}}_{2}\) QSL have been shown in kagome-lattice frustrated magnetic models with (2+1)D O(2)* universality and η* ~ 1.514,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 site20 in certain parameter regimes. Numerical studies20,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 \({{\mathbb{Z}}}_{2}\) QSLs in Rydberg atom arrays. In recent works21,22,28, a diagonal loop operator has been used to detect the \({{\mathbb{Z}}}_{2}\) QSL phase in a model of Rydberg atoms with emergent gauge-charged Ising matter degrees of freedom25. 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 \({{\mathbb{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 \({{\mathbb{Z}}}_{2}\) gauge field. Our results reveal that the LN–VP first-order 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 point61,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 lattice21. Our study on the triangular-lattice quantum loop model connects to the case where the atoms are placed on the sites of the kagome lattice20. Investigation of the LN, hidden-order VP, and \({{\mathbb{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 QSL25 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 system25,28.

Method

Sweeping cluster algorithm

This is a quantum Monte Carlo method developed by the authors which can work well in constrained spin models28,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 constraints28 or higher-order constraints47.