March 22, 2023

We investigate the nature of quantum phases arising in chiral interacting Hamiltonians recently realized in Rydberg atom arrays. We classify all possible fermionic chiral spin liquids with \(\mathrm{U}(1)\) global symmetry using parton construction on the honeycomb lattice. The resulting classification includes six distinct classes of gapped quantum spin liquids: the corresponding variational wave functions obtained from two of these classes accurately describe the Rydberg many-body ground state at \(1/2\) and \(1/4\) particle density. Complementing this analysis with tensor network simulations, we conclude that both particle filling sectors host a spin liquid with the same topological order of a \(\nu=1/2\) fractional quantum Hall effect. At density \(1/2\), our results clarify the phase diagram of the model, while at density \(1/4\), they provide an explicit construction of the ground state wave function with almost unit overlap with the microscopic one. These findings pave the way to the use of parton wave functions to guide the discovery of quantum spin liquids in chiral Rydberg models.

*Introduction.–* There is presently considerable interest in studying strongly correlated phases of matter in synthetic quantum systems based on Rydberg atom arrays [1], [2]. Stimulated by early experiments realizing symmetry-protected topological phases in one dimension [3], these platforms are now able to realize frustrated Hamiltonian dynamics in two dimensions [4]–[6], thus providing unparalleled opportunities to realize quantum spin
liquids (QSLs) – elusive, exotic states of matter which have captivated the attention of physicists for decades [7]–[10]. One route to access QSLs is based on the realization of frustrated Ising models in the so-called frozen gas regime [2], [11], [12]:
several theoretical works have proposed different realistic scenarios for both gapped and gapless phases of matter [13]–[16], with pioneering experiments already reporting evidence for deconfinement [6]. These models resemble closely situations investigated in the context of quantum dimer models [17],
providing direct link between gauge theories and experimental settings [18]–[20].

Over the last two years, a new route has been paved in a very different experimental regime, where the dynamics solely takes place within the Rydberg manifold. The resulting Hamiltonians naturally feature various forms of chiral multi-body interactions [21]–[24], which have already been experimentally demonstrated [25]. These classes of dynamics differ fundamentally from traditional Ising- and Heisenberg-type frustrated magnets and, while very promising since they display chiral terms, it is presently not even clear what classes of quantum spin liquids these can stabilize, and in which parameter regimes they might be observed.

In this work, we provide a general framework to describe chiral spin liquids (CSLs) in Rydberg atom honeycomb arrays. This framework is based on a systematic CSLs classification [26]–[28] using a fermionic spinon construction [29], [30] that yields Gutzwiller-projected parton wave function *ansätze* for the many-body ground state. The resulting classification differs substantially from
those of Heisenberg-type regimes: it rules out the possibility of gapless Dirac spectrum, while predicting several, distinct topological phases.

Combining variational wave functions obtained from the classification with exact diagonalization (ED) methods, we demonstrate that the former capture the intermediate liquid regime of the chiral Rydberg model at both 1/2 [24] and 1/4 density [23], which - surprisingly - encode the same form of topological order: a two-fold ground state degeneracy and a fractionalized Chern number of \(C=1/2\) per state. These two CSLs represent distinct phases, as they are characterized by different projective representations of the lattice symmetries in the underlying fermionic spinon space. Remarkably, the model at 1/2 density corresponds to integer filling of the single-particle band. The CSL found in our model thus represents an interaction-driven topological phase generated from a trivial band insulator [31], setting an open quest recently put forward in Ref. [24]. We then corroborate the topological character of this phase by analyzing the pattern of currents at the edges of a cylinder using the density-matrix-renormalization-group (DMRG) [32]–[35], which shows substantial counter-propagating nearest-neighbor edge currents, and offers a simple mean to experimentally detect such phase. In the 1/4-density case, our results allow to frame the recent observation of CSL [23] within a rigorous classification, as well as providing a genuine understanding of the system wave function (which, even for systems up to 40 spins, reaches overlaps exceeding 0.96).

We consider a system with atoms arranged on a honeycomb lattice. In addition to the atomic ground state \(|0\rangle\), each atom can be found in two Rydberg states, \(| + \rangle\) \(=\) \(a^\dagger | 0 \rangle\) and \(| - \rangle\) \(=\) \(b^\dagger | 0 \rangle\) (differing, e.g., by their electronic angular momentum). Atomic motion is irrelevant on typical experimental timescales: below, we focus on the dynamics of the Rydberg excitations. The model Hamiltonian is [21], [22]:

\[\begin{align} \nonumber H_0 = &\sum_{i \neq j} \begin{pmatrix} a_i \\ b_i \end{pmatrix}^\dagger \begin{pmatrix} -t^a_{ij} && w_{ij} e^{-i 2\phi_{ij}} \\ w_{ij} e^{i 2\phi_{ij}} && -t^b_{ij} \end{pmatrix} \begin{pmatrix} a_j \\ b_j \end{pmatrix} \\ &+ \frac{\mu}{2} \sum_i (n^a_i- n^b_i). \label{eq58fullModel} \end{align}\tag{1}\]

The first term represents hopping processes (of excitations) between different sites, with the real hopping conserving the internal state, and the complex hopping resulting in a change of the internal state. The amplitudes for real and complex hoppings are given by \(t^a_{ij}\), \(t^b_{ij}\), and \(w_{ij} e^{\pm i 2\phi_{ij}}\) respectively, with \(\phi_{ij}\) being the polar angle between the two Rydberg atoms on the sites \(i\) and \(j\). All the amplitudes scale as \(1/r_{ij}^3\). The second term represents the energy difference between the two internal states, with \(n^a_i\) and \(n^b_i\) being the particle number operators for the \(| + \rangle\) and \(| - \rangle\) states, respectively.

Here, we focus on the regime \(\mu \gg t^a_{ij},t^a_{ij},w\), in which case the internal state \(| + \rangle\) can be adiabatically eliminated. We further make an approximation by considering only nearest-neighbor (NN) interactions in Eq. 1, with NN hopping amplitudes \(t^a, t^b,\) and \(w e^{\pm i 2\phi_{ij}}\). At leading order, the effective model is given by [25] \[\begin{align} \nonumber H = &-J \sum_{\langle ij \rangle} b^\dagger_j b_i - 2 g J \sum_{\langle \langle ij \rangle \rangle} b^\dagger_j b_i e^{\pm 2\pi i/3} (1-n_{ij}) + \text{h.c.} \\ &+ 4 g J \sum_{\langle ij \rangle} n_i n_j, \label{eq58model} \end{align}\tag{2}\] where \(J\) \(=\) \(t^b\) and \(g\) \(=\) \(w^2/(2\mu)\). The complex phases \(e^{\pm 2\pi i/3}\) in the next-nearest-neighbor (NNN) hopping are illustrated in Fig. 1a. The NNN hoppings explicitly break time-reversal and reflection symmetries, but preserves their combination. The Hamiltonian has U(1) symmetry related to particle-number conservation. Note that, in the language of spins, the U(1) symmetry corresponds to spin-rotation symmetry around the \(z\)-axis. Hereafter, we set the energy scale to \(J\) \(=\) \(1\).

The phase diagram of the model at 1/2-density has been studied in Ref. [24], where three different phases were found for \(g\) \(\geq\) \(0\). For \(g\) \(\lesssim\) \(0.4\), the phase is a Bose-Einstein condensate (BEC) [36], while for \(g\) \(\gtrsim\) \(0.9\) the phase exhibits spiral or \(120^\circ\) spin order. The intermediate phase between \(0.4\) \(\lesssim\) \(g\) \(\lesssim\) \(0.9\), shows no clear order, and it is believed to be a candidate for a spin liquid. However, its true nature remains unclear, also due to hard-to-interpret spectral properties.

At 1/4-density, Ref. [23] investigated the full model in Eq. 1 and provided clear numerical evidence for a fractional Chern insulating state: that included ground state degeneracy and Chern number compatible with a \(\nu\) \(=\) \(1/2\) bosonic Fractional Quantum Hall (FQH) state. Building on such numerical understanding, we will show below how that reflects into a very clear ansatz for the system wave function, informed by our classification.

In order to construct a spin liquid wave function, a method based on fermionic representation of spins have been introduced in [29], [30]. The main idea is to fractionalize the spin \(S=1/2\) operators into fermionic spinon operators as \(S^a\) \(=\) \(1/2 f_i^\dagger \sigma^a_{ij} f_j\), where \(\sigma^a\) are Pauli matrices, with the constraint of one spinon per site. It is convenient to introduce a two-component spinor \(\Psi\) \(=\) \((f_{\uparrow} \;f_{\downarrow}^\dagger)^T\). Directly rewriting the spins in terms of spinons gives rise to quartic spinon interactions, which after performing mean-field approximation, leads to a quadratic spinon Hamiltonian \[\label{eq58mf95Ham} H_{MF} = \sum_{ij} \Psi^\dagger_i u_{ij} \Psi_j + \text{h.c.},\tag{3}\] where \(u_{ij}\) are the mean-field amplitudes. The spinon interactions include hopping and pairing terms. The mean-field Hamiltonian is invariant under global spin rotation around the \(z\)-axis [37], [38]. The matrices \(u_{ij}\) can be written as \(u_{ij}\) \(=\) \(u_{ij}^\mu \sigma^\mu\), where \((\sigma^\mu)\) \(=\) \((i\tau^0,\tau^a)\), \(u_{ij}^\mu\) are complex parameters and \(\tau^a\) are Pauli matrices. Real \(u_{ij}^\mu\) correspond to singlet terms, while imaginary \(u_{ij}^\mu\) correspond to triplet terms [39]. Different mean-field ansätze are described by different (gauge-inequivalent) \(\{u_{ij}\}\) on the links of the lattice. Finally, a physical spin state is obtained by applying Gutzwiller projection \(|\psi \rangle\) \(=\) \(P_G |\psi_{MF} \rangle\), with \(P_G\) \(=\) \(\prod_i n_i (2-n_i)\), to the mean-field ground state \(|\psi_{MF} \rangle\).

A method to systematically classify all possible spin liquids within this parton construction has been introduced by Wen [26], [27], [40], based on projective symmetry groups (PSG). It has subsequently been extended to classify spin liquid phases in the absence of time-reversal (i.e., CSL) [28] and SU(2) spin-rotation [37], [38] symmetries. Here, we are interested in a CSL which breaks time-reversal and reflection symmetries but preserve their combination, and which preserves U(1) spin-rotation symmetry. Such chiral mean-field states are stable beyond mean-field treatment, as the mean-field gauge fluctuations are gapped out by the Chern-Simons mechanism [41].

The PSG classification of symmetric spin liquids on the honeycomb lattice has been worked out in [42], where 160 different algebraic PSG’s are found. In the absence of time-reversal symmetry, we do not need to specify the SU(2) representation of the time-reversal operation [28]. Thus, we find that the number of different classes of algebraic PSG is reduced to 24 (see [43]). Each PSG is characterized by the representations of reflection, \(g_\sigma(A,B)\), and \(\pi /3\) rotation, \(g_R(A,B)\), for each sublattice \(A\) and \(B\) (for more technical details, see [43]).

Given the model that we study, we focus on those ansätze that have nonzero mean-field amplitudes on the NN and NNN links. This leaves 6 distinct ansätze, which are listed in the Fig. 1b. The last 2 columns indicate the symmetry-allowed terms in the mean-field Hamiltonian on the NN and NNN links [44]. Their amplitudes are taken as variational parameters in the following section.

Note that if the ansätze are restricted to NN interactions, the mean-field states are gapless with Dirac spectrum (in particular, ansatz no. 1 corresponds to the SU(2) algebraic spin liquid (ASL) state discussed in [45], or equivalently the u-RVB state discussed in [42]). Thus, the resulting states after Gutzwiller projection describe a Dirac spin liquid (DSL). However, this DSL ansatz submanifold preserves time-reversal, which is explicitly broken by our Hamiltonian. This excludes the possibility of a DSL being stabilized in chiral systems such as our model.

To determine whether the intermediate phase of the model in Eq. 2 is described by one of the ansätze above, we optimize the variational parameters by maximizing the overlap of the exact ground state of the Hamiltonian with the wave function ansatz, for each of the 6 ansätze. The optimization is performed on a 16-site cluster at \(g\) \(=\) \(0.7\). We find that the ansatz with the largest overlap with the ground state at 1/2-density is ansatz no. 1. We have also checked that the optimal parameters do not differ much from the optimal parameters on the smaller clusters.

The corresponding ansatz with the optimal parameters describes a gapped phase, with two completely filled bands separated by a finite energy gap to the lowest unfilled band. The Gutzwiller projection of such an ansatz has been shown to yield a topological CSL [41], [46], [47], which is a lattice analogue of \(\nu\) \(=\) \(1/2\) bosonic FQH Laughlin state [48]. The topological nature of such CSL is manifested by the two-fold topological degeneracy of states on a torus. These degenerate states can be constructed by threading fluxes along the non-contractible loops on a torus [49], which can be implemented by twisting the boundary conditions of the spinons, \(\Psi\rightarrow e^{i\theta/2}\Psi\). Although there are four states that can be constructed with \(\theta_x,\theta_y\) \(\in\) \(\{0,\pi \}\), they only span a 2-dimensional space, resulting in two topological states. We verify this numerically by computing the overlap matrix for the four states, defined as \(O_{ij}\) \(=\) \(\langle \psi_j | \psi_i \rangle\). We find that the rank of the overlap matrix is 2, within a numerical accuracy of \(10^{-2}\). The 2 independent states are then constructed from the eigenvectors of the overlap matrix with non-zero eigenvalues. Furthermore, we computed the many-body Chern number by numerically integrating the Berry curvature and obtain fractionalized \(C\) \(=\) \(1/2\) per state [43]. The Chern number, along with the two-fold degeneracy, are consistent with the properties of \(\nu\) \(=\) \(1/2\) bosonic Laughlin state.

Next, we compute the overlap of the ground state of the Hamiltonian with the two topological states [50], [51], \[\mathcal{O}^{ED}_{GW}=\sqrt{|\langle \psi_{ED} | \psi_{GW}^{1} \rangle|^2 + |\langle \psi_{ED} | \psi_{GW}^{2} \rangle|^2 }.\] where \(|\psi_{ED} \rangle\) is the ground state obtained with ED and \(|\psi_{GW}^{1}\rangle\) and \(|\psi_{GW}^{2}\rangle\) are the two topological states. The results are shown in Fig. 2 for different system sizes. It can be seen that the overlap remains large in the middle phase with increasing system size, indicating that the ground state is strongly related to the topological states. Interestingly, on a 32-site cluster, we find that the ground state in the intermediate phase is not in the rotation-neutral sector. Specifically, if we take the operator \(R_{\pi/3}\) which generates a \(\pi/3\) rotation around the center of a honeycomb plaquette, then \(R_{\pi/3} | \psi_{ED} \rangle\) \(=\) \(e^{-2\pi i/3} | \psi_{ED} \rangle\). Remarkably, one of the two topological states is in the same rotation sector as the ground state of the Hamiltonian, i.e., the eigenvalue of \(R_{\pi/3}\) is \(e^{-2\pi i/3}\). This nontrivial observation strongly supports the hypothesis that the intermediate phase is described by the wave function ansatz.

To compare with the 1/2-density case, we performed the same parameter optimization procedure for the 1/4-density case. In [23], it was shown that a CSL emerges at 1/4-density for the full model in Eq. 1. For the effective model in Eq. 2, we found that the CSL phase emerges in a narrow range around \(g\) \(=\) \(0.2\). For 1/4-density, a gapped phase can be obtained within the parton construction when the mean-field ansatz has a doubled unit-cell. We obtain large overlaps with the ansatz no. 4 at small-size clusters in the CSL phase. Fig. 2b shows the overlaps for 1/4-density with the optimized parameter for different system sizes. We found that the overlap remains huge in the CSL phase, reaching 0.96 at the largest system size we considered, \(L\) \(=\) \(40\).

In Fig. 3, we present the excitation spectrum in the momentum sector \(k\) \(=\) \((0,0)\) at \(g\) \(=\) \(0.1\) and \(g\) \(=\) \(0.7\), along with the overlaps \(\mathcal{O}^{ED}_{GW}\) for each eigenstate. Note that, since the two topological states \(|\psi_{GW}^{1}\rangle\) and \(|\psi_{GW}^{2}\rangle\) lie in the momentum sector \(k\) \(=\) \((0,0)\), only the eigenstates in this sector can have non-zero overlap. In agreement with [24], we observe no approximate two-fold degeneracy in the ED spectra, which would have been expected in a CSL. Nevertheless, this can be attributed to finite-size effects, which may significantly modify the spectra on small-size clusters. It is therefore possible that one of the low-lying states corresponds to another topological ground state, which becomes degenerate with the true ground state in the thermodynamic limit. To test this hypothesis, it is useful to examine the overlaps of the low-lying levels. If an eigenstate describes the topological ground state of the CSL, it would have a sizable overlap with the wave function ansatz. Indeed, at \(g\) \(=\) \(0.7\), we observe that the overlap is highest for the ground state, and that there is a low-lying state with a modest overlap. In contrast, at \(g\) \(=\) \(0.1\), the overlaps do not exhibit any clear pattern for each system size.

In [24], a DSL was proposed as one of possible scenarios, based on the observation that the gap to the first excited state varies significantly with twist angle when imposing twisted boundary conditions. In light of this, we analyze the excitation gaps as a function of twist angles \(\theta_{1,2}\) along the lattice vectors \(\vec{a}_{1,2}\). In Fig. 4, we show the excitation spectra obtained with ED on a 24-site cluster in the intermediate phase. We observe that while the first gap may become very small at some isolated points in twist space, the second gap and charge gap remain wide open. This contrasts with the expected behavior of a DSL, where all gaps would exhibit a vanishing behavior with respect to twist angles [52]–[54]. In addition, the transfer matrix spectrum does not show any signatures of Dirac cones [43]. This is consistent with the known instability of DSL against time-reversal symmetry breaking perturbations.

Quantum Hall states can be identified through the pattern of the currents in the system. With a finite gap in the bulk and gapless edge excitations on the boundary, it is expected that the currents are large at the edges and vanish in the bulk. Furthermore, the current can be readily measured in experiments, making it a convenient tool for diagnosing the phase in experimental setups. The NN and NNN currents can be derived using the continuity equation, resulting in \[\begin{align} & \mathcal{J}_{\text{nn}}=i\langle b^\dagger_j b_i - b_j b^\dagger_i \rangle, \nonumber \\ & \mathcal{J}_{\text{nnn}}=2 g i e^{\pm 2\pi/3 i}\langle (1-n_{i j}) (b^\dagger_j b_i - b_j b^\dagger_i) \rangle, \label{eq58current} \end{align}\tag{4}\] where \(i\) and \(j\) are nearest-neighbors or next to nearest-neighbors respectively.

We perform DMRG simulations on a finite cylinder to compute the edge currents as a function of \(g\), as shown in Fig. 1c [55]. It can be seen that the NN edge currents, computed across two rungs in one of the edges, in the intermediate phase are significantly larger compared to those in the neighboring phases. The transition points are in good agreement with those found in [24]. Furthermore, we show the full (NN and NNN) current profile in the intermediate phase (\(g\) \(=\) \(0.74\)) in Fig. 1d. It is clear that, in the CSL phase, large NN currents are only observed at the edges, while they vanish in the bulk (see [43] for the same in other phases). This result further confirms the FQH nature of the intermediate phase. Fig. 1c also exhibits signatures of an another intermediate phase with non-zero edge currents for \(0.25 \lesssim g \lesssim 0.4\). While we discuss more in detail this extra phase in the supplemental material [43], we leave its full characterization for future work.

In this work, we systematically classify CSLs on the honeycomb lattice using the PSG analysis. We show that the CSL wave functions constructed from the Gutzwiller-projected parton wave functions are able to capture the intermediate disordered phase in chiral Rydberg atom arrays. In particular, our results resolve the previously unclear nature of the intermediate phase found in Ref. [24]. Our work provides a general framework which can be utilized to search for CSLs in other lattice models. Given the fast experimantal progress in the field, it would be interesting to extend our approach to other lattices, which are immediately available in tweezer arrays [2], [56].

We thank F. Becca, M. Fleischhauer, and S. Ohler for insightful discussions. M.D. and P.S.T. are grateful to Y. Iqbal for explanations and illuminating clarifications on the parton construction. The work of M.D. and P.S.T. was partly supported by the ERC under grant number 758329 (AGEnTh), by the MIUR Programme FARE (MEPH), and by the European Union’s Horizon 2020 research and innovation programme under grant agreement No 817482 (Pasquans). P.S.T. acknowledges support from the Simons Foundation through Award 284558FY19 to the ICTP. G.G. acknowledges support from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC2111 – 390814868 and from the ERC grant QSIMCORR, ERC-2018-COG, No. 771891. T.C. acknowledges the support of PL-Grid Infrastructure for providing high-performance computing facility for a part of the numerical simulations reported here.

in 1,...,

[1]

C. Gross and I. Bloch, Quantum simulations with ultracold atoms in optical lattices, https://doi.org/10.1126/science.aal3837.

[2]

A. Browaeys and T. Lahaye, Many-body physics with individually controlled Rydberg atoms, https://doi.org/10.1038/s41567-019-0733-z.

[3]

S. de Léséleuc, V. Lienhard, P. Scholl, D. Barredo, S. Weber, N. Lang, H. P. Büchler, T. Lahaye, and A. Browaeys, Observation of a symmetry-protected
topological phase of interacting bosons with Rydberg atoms, https://doi.org/10.1126/science.aav9105.

[4]

P. Scholl, M. Schuler, H. J. Williams, A. A. Eberharter, D. Barredo, K.-N. Schymik, V. Lienhard, L.-P. Henry, T. C. Lang, T. Lahaye, A. M. Läuchli, and A. Browaeys, Quantum
simulation of 2D antiferromagnets with hundreds of Rydberg atoms, https://doi.org/10.1038/s41586-021-03585-1.

[5]

S. Ebadi, T. T. Wang, H. Levine, A. Keesling, G. Semeghini, A. Omran, D. Bluvstein, R. Samajdar, H. Pichler, W. W. Ho, S. Choi, S. Sachdev, M. Greiner, V. Vuletić, and M. D.
Lukin, Quantum phases of matter on a 256-atom programmable quantum simulator, https://doi.org/10.1038/s41586-021-03582-4.

[6]

G. Semeghini, H. Levine, A. Keesling, S. Ebadi, T. T. Wang, D. Bluvstein, R. Verresen, H. Pichler, M. Kalinowski, R. Samajdar, A. Omran, S. Sachdev, A. Vishwanath, M. Greiner,
V. Vuletić, and M. D. Lukin, Probing topological spin liquids on a programmable quantum simulator, https://doi.org/10.1126/science.abi8794.

[7]

L. Savary and L. Balents, Quantum spin liquids: a review, https://doi.org/10.1088/0034-4885/80/1/016502.

[8]

X.-G. Wen, Colloquium: Zoo of quantum-topological phases of matter, https://doi.org/10.1103/RevModPhys.89.041004.

[9]

C. Lacroix, P. Mendels, and F. Mila, eds., https://doi.org/10.1007/978-3-642-10589-0(Springer Berlin Heidelberg, 2011).

[10]

R. Moessner and J. E. Moore, https://doi.org/10.1017/9781316226308(Cambridge University Press, 2021).

[11]

M. D. Lukin, M. Fleischhauer, R. Cote, L. M. Duan, D. Jaksch, J. I. Cirac, and P. Zoller, Dipole Blockade and Quantum Information Processing in Mesoscopic Atomic Ensembles,
https://doi.org/10.1103/PhysRevLett.87.037901.

[12]

P. Schauss, Quantum simulation of transverse Ising models with Rydberg atoms, https://doi.org/10.1088/2058-9565/aa9c59.

[13]

R. Samajdar, W. W. Ho, H. Pichler, M. D. Lukin, and S. Sachdev, Quantum phases of Rydberg atoms on a kagome lattice, https://doi.org/10.1073/pnas.2015785118.

[14]

R. Verresen, M. D. Lukin, and A. Vishwanath, Prediction of Toric Code Topological Order from Rydberg Blockade, https://doi.org/10.1103/PhysRevX.11.031005.

[15]

G. Giudici, M. D. Lukin, and H. Pichler, Dynamical Preparation of Quantum Spin Liquids in Rydberg Atom Arrays, https://doi.org/10.1103/PhysRevLett.129.090401.

[16]

G. Giudice, F. M. Surace, H. Pichler, and G. Giudici, Trimer states with \(\mathbb{Z}_{3}\) topological order in Rydberg atom arrays,
https://doi.org/10.1103/PhysRevB.106.195155.

[17]

R. Moessner and S. L. Sondhi, Ising models of quantum frustration, https://doi.org/10.1103/PhysRevB.63.224401.

[18]

A. W. Glaetzle, M. Dalmonte, R. Nath, I. Rousochatzakis, R. Moessner, and P. Zoller, Quantum Spin-Ice and Dimer Models with Rydberg Atoms,
https://doi.org/10.1103/PhysRevX.4.041037.

[19]

P. S. Tarabunga, F. M. Surace, R. Andreoni, A. Angelone, and M. Dalmonte, Gauge-Theoretic Origin of Rydberg Quantum Spin Liquids,
https://doi.org/10.1103/PhysRevLett.129.195301.

[20]

R. Samajdar, D. G. Joshi, Y. Teng, and S. Sachdev, Emergent \(\mathbb{Z}_{2}\) Gauge Theories and Topological Excitations in Rydberg Atom
Arrays, https://doi.org/10.1103/PhysRevLett.130.043601.

[21]

D. Peter, N. Y. Yao, N. Lang, S. D. Huber, M. D. Lukin, and H. P. Büchler, Topological bands with a Chern number \(C=2\) by dipolar exchange
interactions, https://doi.org/10.1103/PhysRevA.91.053617.

[22]

S. Weber, S. de Léséleuc, V. Lienhard, D. Barredo, T. Lahaye, A. Browaeys, and H. P. Büchler, Topologically protected edge states in small Rydberg
systems, https://doi.org/10.1088/2058-9565/aaca47.

[23]

S. Weber, R. Bai, N. Makki, J. Mögerle, T. Lahaye, A. Browaeys, M. Daghofer, N. Lang, and H. P. Büchler, Experimentally Accessible Scheme for a Fractional Chern Insulator in
Rydberg Atoms, https://doi.org/10.1103/PRXQuantum.3.030302.

[24]

S. Ohler, M. Kiefer-Emmanouilidis, and M. Fleischhauer, Quantum spin liquids of rydberg excitations in a honeycomb lattice induced by density-dependent peierls phases,
https://doi.org/10.1103/PhysRevResearch.5.013157.

[25]

V. Lienhard, P. Scholl, S. Weber, D. Barredo, S. de Léséleuc, R. Bai, N. Lang, M. Fleischhauer, H. P. Büchler, T. Lahaye, and A. Browaeys, Realization of a Density-Dependent
Peierls Phase in a Synthetic, Spin-Orbit Coupled Rydberg System, https://doi.org/10.1103/PhysRevX.10.021031.

[26]

X.-G. Wen, Quantum orders and symmetric spin liquids, https://doi.org/10.1103/PhysRevB.65.165113.

[27]

X.-G. Wen, Quantum order: a quantum entanglement of many particles, https://doi.org/10.1016/s0375-9601(02)00808-3.

[28]

S. Bieri, C. Lhuillier, and L. Messio, Projective symmetry group classification of chiral spin liquids, https://doi.org/10.1103/PhysRevB.93.094437.

[29]

G. Baskaran and P. W. Anderson, Gauge theory of high-temperature superconductors and strongly correlated Fermi systems, https://doi.org/10.1103/PhysRevB.37.580.

[30]

G. Baskaran, Z. Zou, and P. Anderson, The resonating valence bond state and high-Tc superconductivity A mean field theory,
https://doi.org/10.1016/0038-1098(87)90642-9.

[31]

In a different context, another example of such occurrence on a Kagome lattice is reported in Ref. .

[32]

S. R. White, Density matrix formulation for quantum renormalization groups, https://doi.org/10.1103/PhysRevLett.69.2863.

[33]

S. R. White, Density-matrix algorithms for quantum renormalization groups, https://doi.org/10.1103/PhysRevB.48.10345.

[34]

U. Schollwöck, The density-matrix renormalization group in the age of matrix product states, https://doi.org/10.1016/j.aop.2010.09.012.

[35]

R. Orús, A practical introduction to tensor networks: Matrix product states and projected entangled pair states,
https://doi.org/10.1016/j.aop.2014.06.013.

[36]

As we report in the supplemental material , an extra intermediate phase occurs for \(0.25 \lesssim g \lesssim
0.4\).

[37]

T. Dodds, S. Bhattacharjee, and Y. B. Kim, Quantum spin liquids in the absence of spin-rotation symmetry: Application to herbertsmithite,
https://doi.org/10.1103/PhysRevB.88.224413.

[38]

J. Reuther, S.-P. Lee, and J. Alicea, Classification of spin liquids on the square lattice with strong spin-orbit coupling, https://doi.org/10.1103/PhysRevB.90.174417.

[39]

The “singlet” and “triplet” terminology is derived from the discussion of SU(2) spin-rotation symmetric Hamiltonians, which have been extensively studied in the
literature. If the state does not spontaneously break the spin-rotation symmetry, only singlet terms are present. However, in the presence of spin-rotation symmetry-breaking perturbations, both terms may be present, and more generally the mean-field
Hamiltonian has to be written in terms of a four-component spinor \(\Psi=(f_{\uparrow} f_{\downarrow}^\dagger f_{\downarrow} -f_{\uparrow}^\dagger)^T\). If the spin-rotation is only broken down to U(1), it is still possible
to use a two-component spinor representation as in Eq. \(\eqref{eq:mf_Ham}\), where the singlet and triplet terms can be present without mixing.

[40]

X.-G. Wen, https://doi.org/10.1093/acprof:oso/9780199227259.001.0001(Oxford University Press, 2007).

[41]

X. G. Wen, F. Wilczek, and A. Zee, Chiral spin states and superconductivity, https://doi.org/10.1103/PhysRevB.39.11413.

[42]

Y.-M. Lu and Y. Ran, \({\mathbb{Z}}_{2}\) spin liquid and chiral antiferromagnetic phase in the Hubbard model on a honeycomb lattice,
https://doi.org/10.1103/PhysRevB.84.024420.

[43]

See supplemental material.

[44]

We include all terms that are allowed by symmetry for each link, which in principle can be present in a mean-field state. However, if the mean-field Hamiltonian is restricted in the
range of interactions (such as up to NNN interactions in our case), some of the terms can be removed by a gauge transformation.

[45]

M. Hermele, SU(2) gauge theory of the Hubbard model and application to the honeycomb lattice, https://doi.org/10.1103/PhysRevB.76.035125.

[46]

Y. Zhang, T. Grover, and A. Vishwanath, Topological entanglement entropy of \({\mathbb{Z}}_{2}\) spin liquids and lattice Laughlin states,
https://doi.org/10.1103/PhysRevB.84.075128.

[47]

Y. Zhang and A. Vishwanath, Establishing non-Abelian topological order in Gutzwiller-projected Chern insulators via entanglement entropy and modular \(\mathcal{S}\)-matrix, https://doi.org/10.1103/PhysRevB.87.161113.

[48]

R. B. Laughlin, Anomalous Quantum Hall Effect: An Incompressible Quantum Fluid with Fractionally Charged Excitations, https://doi.org/10.1103/PhysRevLett.50.1395.

[49]

J.-W. Mei and X.-G. Wen, Modular matrices from universal wave-function overlaps in Gutzwiller-projected parton wave functions,
https://doi.org/10.1103/PhysRevB.91.125123.

[50]

A. Wietek, A. Sterdyniak, and A. M. Läuchli, Nature of chiral spin liquids on the kagome lattice, https://doi.org/10.1103/PhysRevB.92.125122.

[51]

A. Wietek and A. M. Läuchli, Chiral spin liquid and quantum criticality in extended \(S=\frac{1}{2}\) Heisenberg models on the triangular
lattice, https://doi.org/10.1103/PhysRevB.95.035141.

[52]

Y.-C. He, M. P. Zaletel, M. Oshikawa, and F. Pollmann, Signatures of Dirac Cones in a DMRG Study of the Kagome Heisenberg Model,
https://doi.org/10.1103/PhysRevX.7.031020.

[53]

S. Hu, W. Zhu, S. Eggert, and Y.-C. He, Dirac Spin Liquid on the Spin-\(1/2\) Triangular Heisenberg Antiferromagnet,
https://doi.org/10.1103/PhysRevLett.123.207203.

[54]

F. Ferrari, A. Parola, and F. Becca, Gapless spin liquids in disguise, https://doi.org/10.1103/PhysRevB.103.195140.

[55]

We consider two different cylinder geometries for our calculations, namely the geometries I and II, as demarcated in .

[56]

D. Barredo, V. Lienhard, S. de Léséleuc, T. Lahaye, and A. Browaeys, Synthetic three-dimensional atomic structures assembled atom by atom,
https://doi.org/10.1038/s41586-018-0450-2.