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).

0.0.0.1 Model Hamiltonian and phase diagram.–

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.

0.0.0.2 Classification and variational wave functions from parton construction.–

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.

0.0.0.3 Overlaps.–

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.

0.0.0.4 Excitation spectra.–

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.

0.0.0.5 Chiral currents.–

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.

0.0.0.6 Conclusions.–

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].

in 1,...,

References↩︎