Trimer quantum spin liquid in a honeycomb array of Rydberg atoms

Quantum spin liquids are elusive but paradigmatic examples of strongly correlated quantum states that are characterized by long-range quantum entanglement. Recently, the direct signatures of a gapped topological $\mathbb{Z}_2$ spin liquid have been observed in a system of Rydberg atoms arrayed on the ruby lattice. Here, we illustrate the concrete realization of a fundamentally different class of spin liquids in a honeycomb array of Rydberg atoms. Exploring the quantum phase diagram of this system using both density-matrix renormalization group and exact diagonalization simulations, several density-wave-ordered phases are characterized and their origins explained. More interestingly, in the regime where third-nearest-neighbor atoms lie within the Rydberg blockade radius, we find a novel ground state -- with an emergent $\mathrm{U}(1)\times \mathrm{U}(1)$ local symmetry -- formed from superpositions of classical {\it trimer} configurations on the dual triangular lattice. The fidelity of this trimer spin liquid state can be enhanced via dynamical preparation, which we explain by a Rydberg-blockade-based projection mechanism associated with the smooth turnoff of the laser drive. Finally, we discuss the robustness of the trimer spin liquid phase under realistic experimental parameters and demonstrate that our proposal can be readily implemented in current Rydberg atom quantum simulators.

The zoo of quantum spin liquids [1,20,21], however, has many other species including, for instance, states where the invariant gauge group is U(1) instead of Z 2 .Such a U(1) quantum spin liquid is particularly interesting from the perspective of fundamental physics as it hosts an emergent gapless excitation, termed a photon in analogy to conventional electromagnetism [22].Unfortunately, experimental efforts to realize such a phase in rare-earth pyrochlore materials [23] are often complicated by competing microscopic interactions.Standard parton mean-field constructions [17,24] also do not guarantee the stability of the gapless state under monopole fluctuations [25].These varied considerations highlight the importance of discovering robust candidates for studying the rich physics of U(1) gauge theories.
In this work, motivated by the recent experimental progress in trapping neutral atoms, we explore the possibility of finding such spin liquids in programmable Rydberg platforms.We show that, on the honeycomb lattice, strong van der Waals interactions between Rydberg atoms lead to an effective "blockade" constraint [30][31][32][33] that can be mapped to a "trimer" constraint on an underlying triangular lattice [26].Related blockade-induced dimer constraints have played a central role in the recent proposals of gapped dimer-RVB Z 2 spin liquids of Rydberg atoms on both the kagome [34][35][36] and the ruby [37] lattice.On the honeycomb lattice, however, the possibility for a novel trimer quantum spin liquid (TQSL) state arises owing to the correspondence with trimer-as opposed to dimer-coverings.While such an RVB state of trimers was recently classified as a gapless liquid with an emergent U(1) × U(1) local symmetry [28], a microscopic model supporting such a phase and its corresponding physical realizations has yet to be found.Here, we bridge this gap and present conclusive evidence of the realization of TQSLs with honeycomb Rydberg arrays for a range of experimentally relevant parameters.
In order to explore the nature of the ground states, we turn to density-matrix renormalization group (DMRG) FIG. 1. Trimer quantum spin liquid (TQSL) on a honeycomb lattice of Rydberg atoms.(a) In the regime where three nearest-neighbor Rydberg atoms are within the blockade radius (k = 3 shell, shown in orange), the blockade-obeying configurations map exactly to the trimer coverings of the triangular lattice with vertices at the centers of the hexagons [26].
The number of trimer coverings on the triangular lattice grows exponentially with the system size [27].(b) The TQSL is an equal superposition of exponentially many trimer coverings (one covering shown with filled triangles), and is characterized by a U(1) × U(1) local symmetry due to the tripartite nature of triangular lattice with respect to trimers [28].For the tripartition and the trimer configuration shown, we assign two sets of electric fields directed from A to B and from B to C sublattices (arrows).The two U(1) degrees of freedom can be related to two conservation laws, as the independent fluxes are equal to the charges NA − NB and NB − NC enclosed by a closed loop.(c) Quantum phase diagram of Rydberg atoms on a 32 × 4 honeycomb lattice retaining three strongest interactions, as obtained by DMRG [29].The boundaries of the three ordered phases (Néel, columnar, and brick) are mapped out by entanglement entropy, energy susceptibility, and fidelity susceptibility peaks (full lines).In addition, a region with a large entanglement entropy is distinguished by fidelity and energy susceptibility (dashed lines) measurements in the regime where third-nearest neighbors are blockaded.The properties of this unordered phase agree with the expected properties of the TQSL state on a finite cylinder.[38,39] simulations and exact diagonalization (ED) of the so-called "PXP model" on finite-size clusters [40].We find clear signatures of the TQSL phase-including, for example, a high fidelity overlap with the perfect TQSL state-and demonstrate its robustness to real experimental conditions that include long-ranged interactions, relaxed boundary conditions, and experimental state preparation protocols.In addition, we find that the experimental protocol leads to an enhancement of TQSL fidelities compared to the ground state.Our understanding of the underlying mechanism thereof leads us to conjecture a universal fidelity enhancement for any state that is a superposition of configurations with the maximum allowed Rydberg excitations (subject to blockade constraints).Finally, we discuss the experimental signatures of the TQSL states and demonstrate that our proposal can be implemented and studied in today's Rydberg atom quantum simulators.
Trimer model mapping.-Weconsider a system of Rydberg atoms arrayed on a honeycomb lattice with the distance between nearest-neighbour sites being a.This can be achieved experimentally by placing the neutral atoms in optical tweezers and arranging them using spatial light modulators, with currently attainable system sizes in excess of 200 atoms [13].The atoms are driven between the ground (|g ) and highly excited Rydberg states (|r ) by a coherent laser drive with Rabi frequency Ω and detuning ∆, leading to the Hamiltonian where is the reduced Planck constant, i denotes the lattice sites, n i ≡ |r i r| counts the occupation of the excited states, and V ij are the van der Waals interactions between atoms in Rydberg states.The van der Waals interactions fall off with the distance between atoms, R, as V (R) = C 6 /R 6 and are central to the phenomenon of Rydberg blockade that we utilize.More precisely, of the neighboring atoms lying within a distance R b (the blockade radius), defined by V (R b ) ≡ Ω, only one can be excited to the Rydberg state, leading to the blockade mechanism.
Choosing a blockade radius such that the k = 3 nearest neighbors of an atom on the honeycomb lattice are blockaded, the Rydberg blockade constraint becomes identical to a trimer constraint.The trimer in question is a covering of three edges forming a triangle within the triangular lattice built from the vertices placed at the centers of the honeycomb lattice.The constraint enforces that no edge or vertex is shared between trimers [26].This trimer mapping is illustrated in Fig. 1(a), where the blockaded neighbors of a central honeycomb atom are shown to match the corresponding forbidden triangular trimers.Thus, the space of maximally filled blockaded configurations-also referred to as the maximum independent set (MIS) subspace [41,42]-can be matched to trimer coverings of the dual triangular lattice.The number of covering configurations is found to scale exponentially with the total system size (area) [27], despite the stiffness of individual trimer moves.
While the classical trimer model is interesting in its own right as it was numerically found to avoid ordering for any density [26], we are interested in a trimer model with quantum fluctuations that can be realized in Rydberg quantum simulators.The possibility for novel physics in the quantum model arises by considering a trimer state that is an equal superposition of all of the exponentially many trimer covering configurations (TC): where D MIS denotes the dimension of the trimer covering subspace (i.e., the number of different trimer coverings), and the sum extends over all the trimer coverings TC i .Such a state was recently classified by Giudice et al. [28] as a gapless U(1) × U(1) spin liquid.To see the emergence of the two local U(1) symmetries, one can tripartition the triangular lattice such that the trimers cover one site of each sublattice (A, B, C) and then assign electric fields on A-B and B-C trimer bonds, as shown in Fig. 1(b) and argued in Ref. 28.The total A-B (B-C) flux through a closed loop is then given by the difference between the number of A and B (B and C) sites enclosed by the loop, showing the presence of two independent U(1) symmetries (see Supplementary Information [43] for more details).The gaplessness of the TQSL state then follows from Polyakov's results [44] and was also confirmed numerically by Ref. 28 (see Fig. 8(d) therein).
While the abovementioned trimer mapping provides a general starting point, the nature of the quantum ground state of the microscopic model has to be carefully investigated to establish the existence of a quantum spin liquid phase, because other possibilities, including trivial disordered or valence bond solid (VBS) states, cannot be excluded a priori.Here, we do so by demonstrating that certain ground states of the Hamiltonian (1) are adiabatically connected to the perfect TQSL state.
DMRG phase diagram.-Wefirst explore the quantum phase diagram of the Rydberg Hamiltonian on the honeycomb lattice using DMRG on long cylinders of finite sizes [29].We use bond dimensions of up to 1800 and retain the three strongest Rydberg interactions in Eq. ( 1), resulting in good convergence.The additional details of the numerical calculations are presented in the Methods section.We use the following quantities to map the boundaries of the different phases: where S vN is the von Neumann entanglement entropy (ρ 0 being the reduced density matrix for half of the system), χ is the energy susceptibility (E 0 being the ground state energy), and F is the fidelity susceptibility (|Ψ being the ground-state wavefunction).The phase diagram is presented in Fig. 1(c).Among the main features of the phase diagram are the presence of three ordered phases in three different blockade regimes.Intriguingly, an additional unordered region with a large entanglement entropy is apparent, clearly separated from the trivial disordered phase by a phase transition.The Rydberg density profile in all of the regions is presented in Fig. 2, while the static structure factors for the different phases are given in the Supplementary Information [43].The first ordered phase appears in the nearest-neighbor (k = 1) blockade regime, and it corresponds to the Néel phase with a staggered order within the honeycomb unit cell.The Néel state hosts a domain wall in the middle of the 32 × 4 honeycomb cluster employed here due to two different domains being preferred by the open boundaries on each end.
The next-nearest-neighbor blockade (k = 2) leads to the stabilization of the columnar phase, characterized by the pattern of Rydberg excitations in Fig. 2. We note, however, that there exists an extensive classical "string" degeneracy in this k = 2 blockade regime (see Supplementary Information [43]).From that classical manifold of states, the columnar pattern is stabilized by both third-neighbor Rydberg interactions and quantum fluctuations.Both of these effects prefer maximal distance between further-neighboring Rydberg excitations (see Ref. 43), a condition satisfied by the columnar state.Incidentally, thermal fluctuations in the classical hardcore boson model on the honeycomb lattice lead to the same ordered phase in this regime [26].
The final ordered phase, which we call the "brick" phase [45], appears in the third-neighbor blockade regime (k = 3), with the Rydberg excitations patterned on the so-called brick lattice [45] (see Fig. 2).Similar to the Néel state, open boundaries prefer different brick domains, leading to a domain wall in the middle of the cluster considered.The appearance of the brick ground state observed in our simulations is a consequence of quantum fluctuations breaking the exponential degeneracy of valid trimer coverings for the k = 3 blockaded classical model via an order-by-disorder mechanism [46,47].The stabilization of the brick phase, in particular, can be understood by noting that the quantum fluctuations prefer maximally flippable configurations.The brick phase indeed satisfies the maximum flippability condition, as also observed in classical simulations of the hard trimer model relevant to twisted bilayer graphene [45].Additional classical interactions on top of the hard trimer model were shown therein to favor maximal flippability and, in turn, also the brick state.
At intermediate detunings, however, an additional region emerges in the k = 3 blockade regime at fillings close to 1/6.This region shows no order and has a high entanglement entropy throughout.On changing the detuning at a fixed blockade radius, the fidelity and energy susceptibilities manifest a clear peak, as presented in Fig. 3, presumably stemming from a nonadiabatic change in the wavefunction compared to the trivial disordered phase.To probe the intrinsic (bulk) nature of the transition, we also calculate the energy susceptibility difference between 32 × 4 and 24 × 4 clusters (χ b ), thus subtracting out the effects of four boundary columns of atoms at each end of the system.The energy susceptibility peak persists after such boundary subtraction with sizable magnitude, pointing to a putative transition in the bulk (inset of Fig. 3).We label this highly entangled region as TQSL in the phase diagram of Fig. 1 due to its separation from the trivial disordered phase and its appearance in the k = 3 regime with the density of ≈ 1/6 expected for a TQSL state.In addition, the scaling of the gap with the system size and the enhanced susceptibility towards boundaryinduced density oscillations (see Supplementary Information [43]) are seemingly suggestive of a gapless state in the thermodynamic limit, broadly consistent with the expectations for a TQSL.We note, however, that the state is generically observed as gapped for finite clusters (since FIG. 3. Transition into the TQSL regime.On going from the trivial disordered phase to larger ∆/Ω, clear peaks are visible in the fidelity (evaluated from Eq. ( 3) for δ = 0.1 on a 32 × 4 cluster) and energy susceptibilities that are preserved upon boundary subtraction (between 32×4 and 24×4 clusters; inset).The second peak at a larger detuning arises from the transition into the (columnar) ordered phase.
the discrete, allowed momenta need not coincide with the gapless point in momentum space) in the absence of twisted boundary conditions or flux insertion [48][49][50].In the remainder of this paper, we focus on this intriguing TQSL regime identified by our DMRG simulations and demonstrate the existence of a true spin liquid ground state.
TQSL in the PXP model.-Inorder to analyze the existence, characterization, and experimental feasibility of preparation of the TQSL state, we perform large-scale exact diagonalization calculations [40].For these simulations, we employ a hard-constraint approximation to the full Rydberg Hamiltonian, known as the PXP model [51], to enable us to reach large system sizes.The essence of this approximation is to eliminate states violating the Rydberg blockade from the Hilbert space.This is achieved by making the first k Rydberg interactions infinite and discarding the longer-range interactions while projecting the Rabi-oscillation term into the subspace of allowed configurations, leading to the Hamiltonian: where P is the projector onto the blockade subspace and ).This allows us to explore the k = 3 blockade regime with ED [40,52] on clusters of up to 60 sites with periodic boundary conditions in both directions, thus also better simulating the bulk of a large system.The PXP approximation has been effectively employed to understand a variety of phenomena in Rydberg systems, including quantum scars [11,53,54], emergent lattice gauge theories [55], and gapped spin liquids on the ruby lattice [15,37].
We map out the ground-state phase diagram of the FIG. 4. TQSL in the PXP model.A robust TQSL is detected in the PXP model [40,52] with up to third-nearest neighbors blockaded.(a) Overlap of the ground state with the perfect TQSL state reveals three distinct regions in the 10 × 6 shape 1 cluster with periodic boundary conditions [43], namely, a trivial phase at small ∆/Ω, a trimer RVB region with high TQSL fidelity in the middle, and a VBS state with decreasing fidelity at high ∆/Ω.This is further confirmed by considering the fidelity susceptibility [evaluated from Eq. ( 3) with δ = 0.05], wherein two clear peaks stemming from the trivial-TQSL and TQSL-VBS transitions appear.(b) The TQSL fidelity in the spin-liquid region is preserved upon increasing the effective system size, i.e., the number of trimer coverings (equivalent to the MIS degeneracy, DMIS), irrespective of cluster size or aspect ratio.The inset shows the different cluster shapes employed in simulations.The cluster from (a) is emphasized with a green triangle.(c) The spin-liquid nature of the state with high TQSL overlap is seen by considering the ground-state superposition structure for the typical 6 × 10 shape 3 cluster at ∆/Ω = 3.0, circled in (b).All MIS configurations have dominant almost-equal weights and equal phases (see Supplementary Information [43]), while non-MIS configurations have an order-of-magnitude smaller weights.
PXP model as a function of the tuning parameter, ∆/Ω.The phase boundaries are determined by considering the overlap of the ground state with the perfect TQSL state in Eq. ( 2) (equal superposition of all MIS states for the given cluster), | Ψ|TQSL |, the fidelity and energy susceptibilities, as well as changes in the low-lying energy spectrum of the model (see Supplementary Information [43]).To explore the robustness of our predictions, we consider three different honeycomb cluster shapes with several system sizes and aspect ratios for each [43].For all the shapes and sizes probed, we find a sizable TQSL region.The TQSL phase is identified by a high fidelity with respect to the perfect TQSL state as well as by explicitly checking that the ground state is predominantly near-equal weight and phase superposition of all trimer configurations [43].An example for a particular cluster is shown in Fig. 4(a), where the TQSL overlap and the fidelity susceptibility are plotted as a function of the detuning.Three regions, separated by fidelity susceptibility peaks, are observed: a trivial disordered phase for small detunings, a trimer RVB phase with a high TQSL overlap in the intermediate regime, and a phase with decreasing fidelity at high detunings pointing to the formation of a VBS state.In the majority of the other clusters explored, the VBS state is completely absent up to ∆/Ω = 5.This behavior, though driven by quantum fluctuations, is reminiscent of the effect of thermal fluctuations in the classical PXP-equivalent model at finite temperatures, which also lacks order for k = 3 at any density [26].Furthermore, increasing the effective system size, as measured by the number of classically degenerate trimer coverings for a given cluster, leads to no drop in the TQSL fidelity in the spin liquid region, as shown in Fig. 4(b).Lastly, examination of the structure of the ground states in the RVB phase reveals that they are predominantly equalweight and equal-phase superposition of all trimer coverings (see Fig. 4(c) and Supplementary Information [43]), cementing the TQSL nature of the state.Dynamical preparation of TQSL states.-Wenow explore the feasibility of preparing the TQSL state with an experimentally relevant quasi-adiabatic protocol, illustrated in Fig. 5(a).The protocol, of total time T , consists of starting from an initial state where all atoms are in |g and increasing Ω to a desired value at a fixed large negative detuning in the first segment of duration 0.1T .This is followed by increasing ∆ to its desired final value, and then finally, an Ω off ramp of length 0.1T at a fixed detuning.The pulses are then smoothed with a Gaussian kernel to eliminate short timescale effects.
The results for the TQSL overlap of the state at the end of the ramp as a function of the total time T are showcased in Fig. 5(b).The obtained fidelities are, in most cases [43], several orders of magnitude above the groundstate fidelity.It also appears that the fidelities can approach arbitrarily close to one with increasing T .The prepared fidelity depends only weakly on detuning within the TQSL phase, while it drops in the trivial and VBS phases [43].This remarkable fidelity enhancement points to an important role that the quasi-adiabatic preparation protocol might play in the preparation of spin liquid states in general and is consistent with the recent results reported in simulations of the Z 2 spin liquid on the ruby lattice [15,56].
Here, we are able to assign the origin of the fidelity enhancement to the off-ramp part of the pulse.As shown in Fig. 5(c), the fidelity reached during the constant-Ω part of the pulse is of the order of the ground-state fidelity.In the off-ramp part (last 0.1T of time), a significant en- showing that the gain compared to the ground state stems mostly from the off-ramp part of the protocol, to the left of the gray dashed lines denoting T down .The off ramp effectively acts as a projector to the MIS subspace [43] in a mechanism that is expected to be a general feature of Rydberg systems.
hancement is seen.This is in agreement with the experimental observations on the ruby lattice [57] and is shown to be valid for both honeycomb and ruby lattices [43].We explain this off-ramp fidelity enhancement by a projection mechanism.First, we note that the time-dependent Hamiltonian during the off ramp can be thought of as a sequence of Hamiltonians with ever-increasing values of the ratio ∆/Ω.This leads to greater penalties for state admixtures with less than the maximum allowed number of Rydberg excitations (non-MIS configurations), thus leading to an effective projection to the MIS subspace.We test this hypothesis by comparing the TQSL overlap of the final prepared state to that of the state obtained by projecting |Ψ(0.9T ) on to the MIS subspace.We find that while the fidelity of the projection is always higher than that of the prepared one, the two approach each other for long total preparation times (Fig. 5(b) and [43]).While such a mechanism appears to be connected to the PXP model's details, we show in the next section that it applies more generally to Rydberg systems in a slightly modified form.
Robustness of the TQSL state.-Given the theoretically expected RVB nature of the TQSL state, it is important to probe its robustness to perturbations stemming from the "tails" of the van der Waals interactions, which break the classical degeneracy of the trimer configurations.These (experimentally relevant) interaction tails discarded in the pure PXP model of Eq. ( 4) are always present in real atomic systems, so their inclusion is necessary in a realistic model.Note that the PXP approximation of the hard blockade was already relaxed completely in the DMRG studies above.Proceeding further, we now add the interaction tails up to R = 3a to the PXP Hamiltonian of Eq. ( 4).The strength of the interaction tails is controlled by the dimensionless parameter R b /a that we explore in the realistic range between 2 and √ 7 for the k = 3 regime.
The resulting ground-state fidelities for several values of R b /a are shown for a typical case in Figs.6(a), with additional clusters presented in the Supplementary Information [43].Without tails, this cluster does not show a VBS phase up to ∆/Ω = 5.0, thus presenting an extended TQSL region with the wavefunction character shown in Fig. 4(c).Taking the tails into account, a sizable spinliquid region, manifested as a high TQSL fidelity plateau, survives to large R b /a at intermediate detunings.The size of the region and the maximum TQSL fidelity reduces upon increasing the strength of the interaction's tails.In addition, this region is now followed by the VBS plateau with fidelities independent of the interaction tails' strengths.
Next, we explore the quasi-adiabatic preparation protocol of Fig. 5(a) for the TQSL state in the presence of long-ranged interaction tails.The optimal TQSL overlap during the preparation protocol is shown in Fig. 6(b) for parameters corresponding to peaks of the ground-state fidelity in the TQSL phase [solid symbols in Fig. 6(a)].We observe that the preparation fidelity-despite being lower than that for the pure PXP case and falling off with increasing R b /a-still significantly outperforms the ground-state fidelity.Unlike the pure honeycomb PXP case, there now exists some optimal value of the total preparation time that is linear in the system size [43], similar to the case of the pure PXP model on the ruby lattice [56].
In order to gain further insight into the fidelity enhancement observed with dynamical state preparation, we also consider the fidelity as a function of the preparation time, as showcased in Fig. 6(c).We find that in the first 0.9T segment of the protocol (denoted by gray dashed lines), the TQSL overlap from this dynamical preparation is similar to the one obtained from the ground state [see Fig. 6(a)]; however, the Ω off-ramp (the last 0.1T ) leads to a significant fidelity enhancement over FIG. 6.Effect of interaction tails.To probe the robustness of the TQSL, we add the interaction tails up to R = 3a to the PXP Hamiltonian.(a) The ground-state overlap with the perfect TQSL state as a function of ∆/Ω for different interaction strengths (R b /a) shows that the spin liquid survives the long-ranged interactions on the typical cluster with the large MIS degeneracy, though with a smaller parameter region for the TSQL phase and lower fidelities.(b) The adiabatic preparation protocol at parameters corresponding to the peak ground-state fidelity in the TQSL region, denoted by filled symbols in (a), still shows maximum preparation fidelities significantly above the ground state ones.In contrast to the system without tails, an optimal total preparation time exists, showing the significance of semi-adiabatic effects.(c) The TQSL overlap during the adiabatic sweep continues to show a sizable gain during the off ramp, with an additional nonmonotonic feature.The behavior during the off ramp in the presence of interaction tails can still be explained by the universal off-ramp projection mechanism [43].
the ground state.Compared to the case without interaction tails, an additional upturn in 1 − | Ψ(t)|TQSL | is consistently observed at the end of the ramp, leading us to consider the optimal protocol fidelity at intermediate times instead of the final prepared fidelity in Fig. 6(b).This upturn can also be explained within the off-ramp projection mechanism.The semi-adiabatic protocol suppresses the destabilizing effect of the interaction tails for intermediate total preparation times, leading to a state before the off ramp (at 0.9T ) that has equal MIS weights and phases, but still large admixtures of non-MIS configurations.Then, at the start of the off ramp, an effective projection to the k = 3 MIS subspace takes place.However, once Ω drops to a value such that the next shell becomes effectively blockaded, i.e., V k+1 (R b /a) = Ω(t * ), the projection to the MIS subspace for k = 4 is the effective description of the off-ramp Hamiltonian evolution.Therefore, if one wishes to optimize for k = 3 ground states, the off ramp should be sharply cut off before t * .This picture is independent of the PXP-type approximations and generalizes well to the preparation of entangled quantum states or to the optimization algorithms [42,58] arising from blockade physics in Rydberg atom simulators.
Discussion.-Thus far, we have shown how a highly entangled TQSL phase can emerge in a honeycomb lattice of Rydberg atoms and presented evidence for its existence and robustness on finite-sized clusters.Exploring the experimentally accessible preparation protocols, we also uncovered an off-ramp fidelity enhancement mechanism potentially relevant to a broad range of quantum state preparation tasks in Rydberg platforms.We now turn to the question of the experimental characterization of this novel TQSL state.
The spin liquid state that we report here is directly accessible in current-generation Rydberg atom simulators that can realize the relevant lattice geometry, achieve the necessary parameter regimes, and employ the quasi-adiabatic preparation protocol.From our simulations, the parameter range to search for the TQSL phase in experiments corresponds to R b /a ≈ 2.0-2.4 and ∆ max /Ω max ≈ 1.0-4.0.We consider 87 Rb atoms and laser coupling a hyperfine ground state to a 70S 1/2 Rydberg state; a realistic choice of Ω max = 4.0 × 2π MHz and C 6 = 8.6×10 5 ×2π MHz µm 6 leads to a honeycomb lattice constant of a ≈ 3.2-3.9µm, therefore easily accommodating L > 200 atoms in a 100 × 100 µm array.The typical preparation times of ∼ 3-5 µs correspond to T /L ∼ 0.1, which, though currently less than optimal, should still lead to sizable TQSL fidelities and are comparable to the ones employed in the ruby-lattice experiments preparing a Z 2 spin liquid phase [15].Further improvements to experimental coherence times by increasing the laser power and moving to larger intermediate-state detunings can extend T /L to the optimal preparation time in experiments.
The main experimental signature to characterize such a state would be obtained by testing its trimer character and the associated U(1) × U(1) symmetry upon sampling the wavefunction in the experimentally accessible Z-basis.For each snapshot of the TQSL state, a corresponding trimer structure with fluxes can be assigned as presented in Fig. 1.One can then test for the U(1) × U(1) conservation law by evaluating the flux enclosed in a closed loop and comparing it to the theoretical expectations based on the occupation of different sublattices within the loop (see Supplementary Information [43]).Additionally, the resonance between different trimer configurations can be explored using X-basis measurements.These measurements would require a global basis rota-tion applied in the experimental protocol for X-loop measurements, akin to the measurements performed for the ruby lattice [15,37].The relevant X-loop operators to be probed are those that flip between different trimer configurations; these can be found by considering the difference between two irregular breathing honeycomb lattices [27,45,59,60] describing the underlying trimer configurations, as sketched in Fig. 8 of the Supplementary Information [43].Such X-loop operators are expected to decay exponentially with the perimeter of the loop only in the TQSL phase, akin to the ruby-lattice X-loop operators.
Other useful probes of the TQSL phase may be more indirect.For instance, the TQSL should be featureless in the bulk, thus setting it apart from proximate ordered phases in measurements of the static structure factor.However, it should also be distinguishable from the trivial disordered state via energy susceptibility measurements that exhibit a peak at the trivial-TQSL transition [43]; note that the energy susceptibility can be experimentally extracted from the total density in the Z-basis, as χ = ∂ n /∂∆.
Our study opens up several new research directions.Theoretically, it would be important to understand whether a U(1) × U(1) spin liquid state can be stabilized in (2+1) dimensions by coupling the gauge field to gapless fermionic matter and to analytically demonstrate the irrelevance of perturbations (in the renormalization-group sense) about the spin-liquid fixed point [25].Numerically, a question for future work would be to establish the TQSL as a stable phase of matter in the thermodynamic limit, perhaps using methods such as infinite DMRG [61] or quantum Monte Carlo, which can provides new insights beyond exact diagonalization calculations that inevitably suffers from finite-size effects.Experimentally, the preparation and characterization of the TQSL state would pave a new path towards exploring its physics, including the robustness of the state, its gauge field dynamics, and the fractionalization of excitations.For example, in order to probe the excitations of the TQSL, Rydberg spectroscopy consisting of time-dependent correlator measurements can be employed using locally addressable Rydberg arrays [62][63][64].Furthermore, it would be interesting to investigate the nature of the quantum critical points leading out of the TQSL phase as well as the associated nonequilibrium quantum many-body dynamics and potential dynamical phase transitions.Finally, the generic understanding of fidelity enhancements developed in this work can be useful for not only adapting the preparation protocol to obtaining strongly correlated states, but also solving hard optimization problems on Rydberg atom simulators.

METHODS
Exact diagonalization.-Exactdiagonalization simulations of the PXP Hamiltonian were performed using the Bloqade [40] and Generic Tensor Networks [52] packages with periodic boundary conditions on a torus.The Hamiltonians in the blockaded subspace were generated by Bloqade's routines and diagonalized using the Lanczos scheme.The perfect TQSL state was generated by finding all MIS configurations and calculating the MIS degeneracy using Generic Tensor Networks.Three distinct shapes were explored with a variety of aspect ratios and system sizes of up to 60 sites.These cluster shapes are shown in Fig. 4(b), which presents an example of a N x = 4, N y = 4 shape 1 cluster, a N x = 6, N y = 2 shape 2 cluster, and a N x = 4, N y = 3 shape 3 cluster.DMRG.-TheDMRG calculations were performed using the ITensor package [29].The geometry studied here consisted of a long cylinder (with open boundary conditions along the shorter edge and periodic along the longer) shown in Fig. 7 for a N x = 12, N y = 4 system, with the phase diagram of Fig. 1 constructed for a 32 × 4 system.Phase diagrams for 24 × 4, 26 × 4, and 32 × 3 cylinders, as well as for a 32 × 4 system with the addition of interaction tails up to fifth-nearest neighbors (R = 3a) were also fully reconstructed, with the phase boundaries qualitatively matching the ones for the 32 × 4 system, including the TQSL region.Other aspect ratios were employed for probing the properties of the TQSL regime, with N x = 16-32 and N y = 3-6.The sweep protocol employed was based on typically O(200) sweeps with bond dimensions in the range 400-1800, depending on the transverse system size (N y ).Initially, many sweeps were performed at a relatively small bond dimension (< 100), with the bond dimension being progressively increased in later stages.A gradually decreasing noise term was added until the final stages of the sweep to prevent the DMRG from being stuck in local minima.The protocol achieved good convergence with typical discarded weights below 10 −10 and the relative change in the ground state energy after the final increase in the bond dimension being below 10 −6 .Supplementary Information: "Trimer quantum spin liquid in a honeycomb array of Rydberg atoms" In this Supplementary Information, we discuss in more detail the experimental detection of the TQSL state and present additional results from both DMRG and ED calculations supporting the TQSL character of the state and its robustness on a wide range of geometries.In Sec.I, we provide details of the trimer mapping relevant for the experimental detection of the TQSL; Sec.II and Sec.III present additional DMRG and ED results, respectively; Sec.IV includes additional data on the dynamical preparation of the TQSL state; finally, Sec.V shows the numerical results for the second-neighbor-blockaded PXP model for comparison.

I. TRIMER MODEL MAPPING AND EXPERIMENTAL SIGNATURES OF THE TQSL
We start by presenting the details of the trimer model mapping relevant to the experimental detection of the TQSL phase.This actually turns out to be closely related to the question of classifying the degenerate trimer configurations.It is known [27,45,59] that all the trimer configurations can be represented by a corresponding irregular honeycomb lattice that tiles the plane, as shown in Fig. 8(a-c).The exponential degeneracy of trimers is then directly related to the "breathing" degree of freedom for each of the irregular honeycombs: the freedom to contract and expand while keeping the aspect ratio fixed.Such a move is local, as it only involves changes in the immediate vicinity of the breathing honeycomb.This has an important consequence for the construction of experimental probes of the TQSL state.In particular, a local operator comprised of spin flips ( i σ x i ) that corresponds to a particular breathing move will, in general, permute within the classes of states in the MIS subspace.Thus, measuring the expectation value of different breathing operators, explicitly constructed from irregular honeycomb mappings, can distinguish between the TQSL and ordered trimer phases.The breathing operators are analogous to the X-loops for the ruby lattice [15].They also naturally require measurements in the X-basis, necessitating a global basis rotation performed after a quench to noninteracting atom regimes or employing hyperfine mapping and single-qubit rotations [64].
The X-loop breathing operators are expected to decay rapidly with increasing loop and system size due to the structure of the trimer subspace and monomer fluctuations.The fluctuation-induced decay is expected to scale exponentially with the loop perimeter only in the TQSL phase, in complete analogy to the ruby-lattice X-loops [15,37].The perimeter-law scaling is showcased for a breathing move in Fig. 8(d), where any difference between two loop-connected configurations is limited to the yellow-shaded region of two triangular layers around the loop's perimeter.The additional prefactor for the X-loop's intensity is related to the trimer subspace's fractionalization into sectors described by topologically distinct irregular honeycomb lattices [59].This prefactor is, in the worst case, proportional to the inverse of the number of topological classes, scaling as 1/L 2 for system size L [59].The constant nature of the prefactor for a given system size and shape still allows for probing the perimeter-law scaling of the X-loops, and thus, the TQSL phase.
A second set of operators can be used to experimentally distinguish between the trivial and TQSL states, allowing for a definite experimental detection of the TQSL phase in conjunction with the breathing operators.These operators check the two U(1) conservation laws [28] by evaluating the electric flux and conserved charges for an arbitrary closed loop, as illustrated in Fig. 8(d).For a given snapshot of the honeycomb lattice, a related triangular lattice can be defined, tripartitioned, and covered with trimers.Electric fields related to the associated U(1) gauge degrees of freedom can then be assigned on the trimers, as described in the main text, following which, their flux (Φ) can be evaluated and compared to the enclosed charge.For A-B electric fields, Z AB = Φ AB − (N A − N B ) equals 0 for the perfect trimer covering, and likewise for B-C.Thus, an expectation value close to 0 of these Z-basis-accessible operators on arbitrary closed loops distinguishes between trivial and trimer states.These operators are direct analogs of the closed Z-loops [15] employed for detection of the Z 2 topological phase on the ruby lattice.
As an example, we numerically evaluate the Z-loop operators across the phase diagram of the 10×6 shape 1 cluster, showcased earlier in Fig. 4(a) of the main text.In particular, Z AB and Z BC are calculated for the Z-loop presented in Fig. 8(d) and shown in Fig. 8(f).The expectation values for both of the loop operators are high in the trivial phase and approach the expected value of zero in the TQSL phase.The second transition between the TQSL and the VBS present on this cluster is not distinguishable from Z-loop measurements alone.[27,45,59] one can start from an ordered trimer state and note that flipping a line of trimers still leads to a valid trimer covering.(b) The same is true for a confluence of three flipped lines ending with a filled (shaded) or empty (unshaded) trimer.(c) All the trimer configurations reachable from the original ordered state can then be described by any irregular honeycomb tiling of the plane.The number of such configurations grows exponentially with the system size.This can be seen by noting that each of the irregular honeycombs has a "breathing" degree of freedom-it can contract and expand while keeping its aspect ratio.Such a move is local (shown in red), changing only the configuration in the vicinity of the breathing irregular honeycomb, making it a potentially useful experimental probe of the resonance between different trimer configurations.(d) This experimental probe can be explicitly constructed for a given breathing move as it is equivalent to a product of σx operators needed to expand the loop.An example of a simple breathing move is shown with configurations before and after the move differing only within the yellow shell around the perimeter of the loop.The moves needed to expand the honeycomb, in this case, lead to the many-body X-loop operator that is a product of 24 individual honeycomb σx operators.(e) The main experimental probe that definitively distinguishes between trivial and TQSL phases entails checking the two U(1) conservation laws [28]

A. Structure factors of the phases
The static structure factor quantifies the correlations between Rydberg excitations in the various phases and presents a direct experimental classification tool for the density-wave-ordered ones.The structure factors were calculated as the Fourier transform of density-density correlations: where r i denote the positions of the atoms.The results are presented in Fig. 9 for representative points in the different phase regions with the same parameters as in Fig. 2 of the main text.For all the structure factors shown, one has to take into account the effect of finite-size clusters with open boundary conditions as well as the choice of one of the degenerate ordered configurations made by DMRG due to boundary conditions that explicitly break the point group symmetries of the lattice.In the experimental scenario, with a larger system, different ordering patterns may occur in different snapshots, averaging out to a symmetric structure factor.From Fig. 9, we see that the disordered phase, as expected, exhibits only a trivial peak at the Γ point of the first Brillouin zone defined by the reciprocal lattice vectors, b 1 = 2π(−1, −1/ √ 3) and b 2 = 2π(0, 2/ √ 3).This is similar to the structure factor in the TQSL regime, which too shows a complete absence of sharp ordering peaks.On top of the expected Γ-point peak, broad features can also be seen in the structure factor, stemming from the short-range correlations induced by the strong k = 3 trimer constraint.
In contrast, the structure factor of the Néel phase displays strong Bragg peaks in the second Brillouin zone because the Néel ordering occurs within the unit cell of the honeycomb lattice.The doubling of the peaks in the second Brillouin zone is a consequence of the two domains present in our numerically obtained Néel state.The columnar and the brick phases, on the other hand, break the translational symmetry of the lattice and can be characterized by considering the structure factor peaks inside the first Brillouin zone.For the columnar phase, the main ordering peaks occur at ±(b 1 /2 + b 2 ) with two additional sets of secondary peaks, the first appearing at ±(3b 1 /4 + b 2 /2) and ±(−b 1 /4 + b 2 ), while the second is at ±b 1 /2.All of these peaks are expected from the real-space structure of the columnar phase.Finally, the main ordering peaks in the brick phase are at ±(2b 1 /3 + b 2 /2) and ±2b 1 /3, with secondary peaks at ±(b 1 /3 + b 2 ), ±(b 1 /3 + b 2 /2), and ±b 2 /2.The brick-phase peaks are all consistent with the real-space ordering observed and show doubling due to the presence of two distinct domains.

B. Gap and density fluctuations in the TQSL regime
In this section, we present some additional properties of the TQSL region observed with DMRG that hint at the possibly gapless character of the state in the thermodynamic limit.Figure 10(a) shows the entanglement entropy, which attains large values in the TQSL region, with the only phase transition discernible being that into the ordered (columnar) phase.In Fig. 10(b), we show the energy gap, ε, as a function of the detuning for several system sizes across the trivial and TQSL regions.The gap is generically expected to be nonzero in a finite system like the one studied.In general, a gapless spin liquid on finite clusters without specially engineered boundary conditions and flux insertion (inaccessible in this case) becomes gapped [48][49][50].Still, a qualitative distinction is apparent between the trivial phase, where the gap is independent of the transverse size, and the TQSL region, where the gap sharply drops with increasing transverse size, alluding to a potentially gapless state in the limit of infinite transverse size.Going from trivial disordered phase to the TQSL region, one detects a strong decay of the gap as a function of ∆/Ω on the one hand, as well as a change in the behavior with increasing transverse system size.While in the trivial phase, the gap is approximately constant for different system sizes, in the TQSL region, the gap generally drops with increasing system size, which might hint at a gapless state in the thermodynamic limit.Inset: The second derivative of the gap (Ny = 5) also displays a peak akin to the energy and fidelity susceptibilities.(c) Similarly, on a plot of the column density fluctuations, a strong susceptibility towards fluctuations is seen on going from the trivial phase (∆/Ω < 1.9) to larger ∆/Ω.Although these fluctuations decay faster than the expected Friedel oscillations for a truly gapless state, the result is broadly consistent with a gapless state in the thermodynamic limit.
True extrapolation to the thermodynamic limit is presently not viable with this dataset of small transverse sizes and increasing the transverse size is, unfortunately, exponentially costly with respect to the DMRG bond dimension.Interestingly, the second derivative of the gap with respect to the detuning also exhibits a peak consistent with those in the energy and fidelity susceptibilities, as shown in the inset of Fig. 10(b).
Another interesting feature of the ground state in the TQSL regime is seen in the column-averaged density oscillations presented in Fig. 10(c).The system, for a given N x , hosts N x /2 columns, and the relative density oscillations across the columns compared to the average density are shown.Going from the trivial phase at ∆/Ω = 1.5 to ∆/Ω = 2.5, one notices a stark increase in the density oscillations' magnitude in both the bulk and at the boundaries.These density oscillations are induced by the open boundaries [66,67] and in a true gapless state, they correspond to the Friedel oscillations expected to decay with a power law.In the small transverse-size systems we employ, an exponential (rather than power-law) decay is observed in the TQSL regime.However, as the gap is reduced, the oscillations are enhanced in the bulk.The possibility of an increase in susceptibility towards boundary-induced oscillations being a signature of the TQSL state is compounded by the observation of the same feature in the cylindrical clusters studied with ED and presented in Sec.III C.

III. EXACT DIAGONALIZATION OF THE k = 3 PXP MODEL
A. Dependence on system size and shape A crucial part of establishing the presence of the TQSL state in ED simulations is to probe different system sizes, aspect ratios, and shapes.The most important takeaway from such an analysis was presented in Fig. 4(b), where it was shown that the large TQSL overlap persists upon increasing the effective system size parameter-the MIS degeneracy-for a given cluster.Our ED results point to a significantly better ground-state spin liquid fidelity and fidelity retention with increasing MIS degeneracy, than observed in the ED simulations of the PXP model on rubylattice clusters for the Z 2 spin liquid [56].Here, in Fig. 11, we show the TQSL overlap as a function of the detuning for a range of clusters with the largest MIS degeneracies.
The typical result valid for the clusters with MIS degeneracies above O(100) as well as for most of the smaller MIS degeneracy clusters is the presence of two phases only: trivial at small detunings and TQSL at large ones, as exemplified by the high TQSL fidelity for the largest detunings probed.The only deviation to this trend occurs for clusters with a large system size (10 × 6) but comparatively small (< 100) MIS degeneracy.This destabilizes the TQSL state at large detunings and leads to the appearance of the VBS phase.This is seen clearly for the 10 × 6 shape 1 cluster, where the TQSL region that survives is evident in the overlap.On the 10 × 6 shape 3 cluster, the TQSL region is still present, as detected by the changing superposition structure of the ground state (see Sec. III B), but only in a diminished detuning range (1.3-1.9 ∆/Ω), making it hard to establish two transitions in the overlap and the fidelity susceptibility.FIG.11.Ground-state fidelity for various system shapes and sizes.The ground-state TQSL overlap is plotted as a function of ∆/Ω for five system sizes with the largest MIS degeneracies for each of the probed shapes.One sees that in all but two of the clusters, a TQSL state is reached with high fidelity at ∆/Ω 1.5, with the VBS phase absent up to the largest detunings probed.The two remaining 10 × 6 clusters have relatively small MIS degeneracies (60 and 90, respectively), leading to a less robust TQSL state, and show the presence of the VBS state at larger ∆/Ω.

B. Structure of the ground states and low-energy spectrum
We now establish additional signatures showing the presence of the TQSL phase in the ED simulations and provide details of the ground-state wavefunctions obtained.In Fig 12(a), the lowest 100 eigenenergies for the 10 × 6 shape 1 cluster are shown.The spectrum displays two transition points, as evidenced by the qualitative changes in energy splittings appearing at the same positions as peaks in the fidelity susceptibility in Fig. 4(a).The first qualitative change can be related to the MIS manifold becoming dominant for the low-energy eigenstate.This is seen in the spectrum as the appearance of a large energy splitting ∝ ∆ between two groups of states corresponding to predominantly |MIS| and |MIS| − 1 superpositions.The energy splittings within the groupings are determined by Ω-induced hoppings between different states with the same number of Rydberg excitations and decay as a high power of Ω/∆.The second transition can be connected to level crossings and energy splittings within the MIS manifold itself, signaling the onset of the VBS order.The trivial-TQSL transition can also be detected in the (experimentally accessible) energy susceptibility, as shown in Fig. 12(b) for two representative clusters, or equivalently, the average density as a function of detuning.The second (TQSL-VBS) transition, however, is hard to detect from the density, as both states have an average density very close to 1/6 as expected for trimer coverings.In Fig. 13, we dwell on the details of the ground states in the TQSL regions of two representative clusters, by showing the weights and phases of the most significant states in the superposition.Since these states lie within the TQSL phase, the highest-weight states belong to the MIS subspace, with all the MIS configurations included in the superposition.Furthermore, the perfect TQSL state should have equal weights and equal phases for all of the MIS configurations.This is indeed very close to the actual ground state realized in the typical cluster of Fig. 13(a) where all the MIS states participate in the ground state with equal phases and with weights within 2 × 10 −4 of the relative spread.The magnitude of the weights themselves is ∼ 0.86 of the perfect TQSL value of 1/ √ 486, with the remaining weights being distributed mostly among states in the |MIS|−1 subspace with an order-of-magnitude smaller individual configuration weights (see Fig. 4(c) in the main text).This is the main cause of the decrease in fidelity relative to the perfect TQSL state and shows the significance of trimer-monomer fluctuations that, in turn, determine the effective physics within the MIS subspace.We emphasize that this remarkable similarity to the perfect TQSL state is shown here for the cluster with the largest MIS degeneracy, and it corresponds to the typical case of clusters with high TQSL fidelities from Fig. 11.The near worst-case scenario exemplified by the cluster of Fig. 13(b) still shows a similar structure with MIS admixtures being dominant with equal phases.The difference arises due to higher weight inequality, with half of the weights being close to the perfect TQSL value and the other half at ∼ 0.75 thereof.In contrast, deep in the VBS region of the 10 × 6 shape 1 cluster [Fig.13(c)], one observes a ground state that is predominantly a superposition of only half of the MIS configurations, with the contribution from the rest of the states in the MIS subspace strongly suppressed such that the next highest weights states arise from the |MIS| − 1 subspace.The |MIS| − 1 configurations are also accompanied by a relative phase of π with respect to the MIS configurations in the wavefunction.

C. Effect of open boundaries
In order to probe the robustness of the TQSL state found in ED numerics as well as to compare more directly between the ED and DMRG simulations, we perform exact diagonalization calculations on clusters with cylindrical boundary conditions.For this purpose, we take the cluster shapes and aspect ratios previously considered with periodic boundaries and relax the open boundaries at the shorter edge of the system.The results are presented in Fig. 14 for several representative clusters.
Inspecting the TQSL overlap as a function of the detuning in Fig. 14(a) shows that the cylindrical TQSL fidelity remains large for the typical clusters, signifying the preservation of the TQSL phase.The fidelity drops somewhat in the less-favorable cluster (10 × 6), which has a high boundary-to-bulk atom ratio (1/3).Even in this case, however, the phase region at ∆/Ω 1.5 can be described as a TQSL from the analysis of the ground-state structure.In Fig. 14(b), we show the density profile in the TQSL phase for two clusters.The similarity to the TQSL excitation density profiles found previously with DMRG is apparent, with the overall density in the bulk approaching 1/6 and the boundary inducing density fluctuations.This is further supported by Fig. 14(c), where we show the column-averaged relative density oscillations.These exhibit a strong enhancement across the the trivial-TQSL transition, in qualitative similarity to the ones observed using DMRG [see Sec.II B and Fig. 10(c)].

D. Effect of interaction tails
The second aspect of robustness particularly relevant to experimental Rydberg systems is the effect of the interaction tails initially discarded in the simulation.While we have already summarized the major results in the main text, here, we present additional details with different system sizes and shapes in mind.The results are plotted in Fig. 15 for several clusters previously studied without tails in Fig. 11.A typical cluster shows the presence of the TQSL region in the intermediate detuning parameter range, with decreasing fidelity as the interaction strength increases.As noted in the main text, at large ∆/Ω, a VBS phase is always present with the same fidelity plateau reached for different values of the interaction tails' strengths.The worst-case scenario of the 10 × 6 shape 1 cluster still shows a TQSL region at intermediate detunings, as exemplified by fidelity peaks and confirmed by considering the wavefunctions' structures, although in a much diminished detuning range and with lower fidelities, vanishing at high R b /a completely.

IV. DYNAMICAL PREPARATION OF TQSL STATES
A. Dependence on detuning, system size, and cluster shape In this section, we provide additional data on the dynamical preparation of TQSL states for different clusters.We first show, in Fig. 16, the maximum TQSL overlap reached during dynamical preparation with fixed total time as a function of the detuning for several system sizes with and without interaction tails.This effectively allows us to extract the dynamical phase diagram.The results without tails are qualitatively similar to the ground-state fidelity from Fig. 15, with a typical 6×10 shape 3 cluster showing a large TQSL region and no VBS phase, while the worst-case 10 × 6 shape 1 cluster shows a dip in the fidelity at large detunings pointing to the onset of the VBS.The prepared fidelities in the TQSL region are significantly above the ground-state ones.In the system with tails, a sizable region of high TQSL fidelity is observed at intermediate detunings, with the eventual onset of the VBS now visible for all clusters.Perhaps most remarkably, the 10 × 6 shape 1 cluster with tails clearly reveals the presence of a sizable TQSL phase with the dynamical preparation scheme, despite the diminished ground-state TQSL regime.
The data for the final prepared TQSL fidelity is shown in Fig. 17(a,b) as a function of the total preparation time FIG.16.Dynamically prepared TQSL fidelity.The maximum TQSL overlap reached during dynamical state preparation with fixed total time, as a function of the detuning for several system sizes with and without interaction tails.While the dynamical phase diagram obtained from the TQSL fidelity is reminiscent of the ground-state ones, the prepared fidelities are significantly higher, and the presence of a sizable TQSL region is evident even in the 10 × 6 shape 1 cluster.
for different clusters with no interaction tails.The typical cluster in the TQSL region displays behavior akin to the 6 × 8 shape 3 cluster presented here and to the cluster from Fig. 5 of the main text.The fidelity appears to be able to approach arbitrarily close to unity with increasing total preparation time.A similar behaviour is observed on the 10 × 6 shape 1 cluster shown in Fig. 17(b).In addition, we also examine two clusters with a qualitatively different behavior.The 10 × 6 and 6 × 6 shape 3 clusters exhibit a limited maximum fidelity and correspondingly, the presence of an optimal total preparation time.The optimal prepared overlap is still in excess of ∼ 0.95 for these cases.The difference is possibly related to the smaller MIS degeneracies for these clusters, which is seen to destabilize the TQSL ground state for the 10 × 6 cases (see Sec. III A).An upper-bounded fidelity that reaches an optimal value at intermediate preparation times is also the typical behavior observed for the k = 3 PXP model without tails on the ruby lattice [56], as well as for the honeycomb lattice with long-ranged interactions explored here.Figure 17(c,d) presents the optimal TQSL overlap reached during the adiabatic preparation sweep as a function of the total preparation time for three additional clusters with different interaction tail strengths.The detuning is selected to correspond to the maxima of the TQSL ground-state fidelities from Fig. 11, which does not generically correspond to the optimal detunings for the dynamical preparation seen in Fig. 16.The obtained results are in agreement with the ones described in Fig. 6 of the main text.The main features observed are: (1) fidelities significantly above the ground-state ones, (2) the existence of an optimal total preparation time, and (3) a decrease in the optimal fidelity with increasing interaction strengths.Moreover, fluctuations in the optimal prepared fidelity are observed for short total preparation times.These stem from the short time-scale effects (see Fig. 18) induced by the preparation protocol that has steeper ramps for smaller T and are more pronounced in smaller clusters.

B. Universal off-ramp projection
Following our previous discussion, we now show more examples of the off-ramp fidelity-enhancement mechanism in action, including the adiabatic preparation of the dimer RVB (dRVB) state on the ruby lattice, as well as supporting evidence for the projection mechanism of the off ramp.With this in mind, we highlight three distinct cases where the fidelity enhancement mechanism applies in Fig. 18: the PXP model on the honeycomb lattice without tails, the PXP model on the honeycomb lattice with interaction tails, and the no-tail PXP model on the ruby lattice.In all of these cases, the same preparation protocol as in Fig. 5 of the main text was used, with ∆/Ω in the main part of the ramp corresponding to the TQSL/dRVB region.The clusters employed for the ruby lattice are the same as the 48-and 36-site clusters used in Ref. 56.In all three scenarios, after the first 0.9T of preparation time (denoted by gray dashed lines), the fidelity obtained with this dynamical preparation scheme is similar to that of the ground state (see Figs. 11,15 and [28]).In the last 0.1T time segment corresponding to the off ramp, a large gain in the fidelity is observed.In the cases without tails, this can yield several orders of magnitude of fidelity enhancement.We note that the off-ramp part of the pulse is responsible for the high RVB fidelities reported in the recent ruby-lattice simulations [56] of the Z 2 spin liquid.More significantly, in the experimental systems with tails, the off ramp still leads to practically relevant fidelity enhancements over the ground state, regardless of the cluster shape and size.An additional upturn appears at the end of the ramp, understanding which in terms of the projection mechanism of the In all the cases presented, significant fidelity gain is obtained in the off-ramp part of the protocol, which can be explained by the projection to the MIS subspace during the off ramp.In the case with interaction tails, after Ω drops below the appropriate value for the next blockade shell, projection to the smaller MIS subspace starts to occur, leading to a final upturn in infidelity.
is not a resonating state, eliminating the possibility for spin-liquid physics in this regime.We contrast this with a recent study of similar string degeneracy on the kagome lattice [34], where the classical string configurations could be mapped to a dimer model on the medial triangular lattice, potentially giving rise to a Z 2 spin liquid; however, no such mapping exists for the honeycomb strings.Nonetheless, a disordered quantum string state arising from highly constrained classical dynamics might be an interesting platform for exploring the recent proposal of emergent glassy dynamics [68] in a nonequilibrium setting.As discussed in the main text, quantum fluctuations break the string degeneracy and ultimately favor a state with ordered antiparallel strings, the columnar state.The choice of antiparallel strings is not surprising, given that the distance between excited atoms is larger for antiparallel than for parallel strings, leading to a greater energy gain from fluctuating configurations (note that thermal fluctuations make the same choice of ordered state in the classical model [26]).The argument above for maximizing the distance between neighboring excitations also explains why the interaction tails prefer antiparallel strings, leading to the same columnar state.
The columnar order sets in only at ∆/Ω > 3.0, as evidenced by the DMRG phase diagram.Here, we explore whether there is a possibility to access the potentially novel string MIS state in the intermediate detuning regime using the k = 2 PXP model.The effective system size as measured by the MIS degeneracy is significantly smaller than in the k = 3 case, a consequence of the degeneracy scaling with the linear system size, making the ED argument for such a state weaker.Still, as shown in Fig. 19(b) for the PXP model without tails, a high string MIS fidelity is observed in the intermediate detuning regime.However, including the interaction tails up to R = 3a leads to a large decrease in fidelity and the loss of the string MIS state, making it unlikely to be realizable in an experiment.The fragility of the string MIS state lies in contrast to the observed remarkable robustness of the TQSL state to interaction tails.

FIG. 2 .
FIG. 2. Phases on the honeycomb lattice.Rydberg excitation density profiles at representative points of the DMRG phase diagram from Fig. 1.While the columnar phase is defectless for this cluster, the Néel and brick phases host a domain wall in the middle of the lattice due to their incommensurability with the open boundary.The two ordered domains on each end of the system are indicated by red boxes, while the primitive unit cells for the ordered phases are delineated in blue.In contrast to the ordered phases, in the TQSL region, a state with no density-wave order and bulk density close to the expected value of 1/6 is observed, on top of which, density oscillations spread inwards from the boundaries.

FIG. 5 .
FIG. 5. Adiabatic preparation of the TQSL.(a) We test the feasibility of preparing the TQSL state with an experimentally relevant adiabatic preparation protocol.The ramp-down time (T down = 0.9T ) is denoted by a dashed gray line.(b) This results in prepared state fidelities that are several orders-of-magnitude better than the ground-state ones, as seen by plotting the prepared TQSL overlap as a function of the total preparation time (in units of 2π/Ω), rescaled by the system size (L = NxNy), for parameters representative of the TQSL phase.The blue dashed line shows the projection of |Ψ(0.9T ) onto the MIS subspace.(c) TQSL overlap during the quasi-adiabatic sweep presented for the filled data points in (b),showing that the gain compared to the ground state stems mostly from the off-ramp part of the protocol, to the left of the gray dashed lines denoting T down .The off ramp effectively acts as a projector to the MIS subspace[43] in a mechanism that is expected to be a general feature of Rydberg systems.

FIG. 7 .
FIG. 7. Cluster shape used in DMRG simulations.DMRG was performed on cylindrical boundary conditions with small transverse cylinder sizes.

FIG. 8 .
FIG. 8. Trimer model mapping and experimental probes of the TQSL.(a) To classify the trimer degeneracy on the triangular lattice[27,45,59] one can start from an ordered trimer state and note that flipping a line of trimers still leads to a valid trimer covering.(b) The same is true for a confluence of three flipped lines ending with a filled (shaded) or empty (unshaded) trimer.(c) All the trimer configurations reachable from the original ordered state can then be described by any irregular honeycomb tiling of the plane.The number of such configurations grows exponentially with the system size.This can be seen by noting that each of the irregular honeycombs has a "breathing" degree of freedom-it can contract and expand while keeping its aspect ratio.Such a move is local (shown in red), changing only the configuration in the vicinity of the breathing irregular honeycomb, making it a potentially useful experimental probe of the resonance between different trimer configurations.(d) This experimental probe can be explicitly constructed for a given breathing move as it is equivalent to a product of σx operators needed to expand the loop.An example of a simple breathing move is shown with configurations before and after the move differing only within the yellow shell around the perimeter of the loop.The moves needed to expand the honeycomb, in this case, lead to the many-body X-loop operator that is a product of 24 individual honeycomb σx operators.(e) The main experimental probe that definitively distinguishes between trivial and TQSL phases entails checking the two U(1) conservation laws[28] relating the electric field (arrows) flux with enclosed charges (number of different sublattice sites) along an arbitrary closed loop (example shown in blue).(f) Two Z-loop operators evaluated for the loop shown in blue in (e) for the 10 × 6 shape 1 cluster from Fig. 4(a) of the main text.The two U(1) conservation laws are violated in the trivial phase, but are approximately satisfied in both the TQSL and VBS phases.
FIG. 8. Trimer model mapping and experimental probes of the TQSL.(a) To classify the trimer degeneracy on the triangular lattice[27,45,59] one can start from an ordered trimer state and note that flipping a line of trimers still leads to a valid trimer covering.(b) The same is true for a confluence of three flipped lines ending with a filled (shaded) or empty (unshaded) trimer.(c) All the trimer configurations reachable from the original ordered state can then be described by any irregular honeycomb tiling of the plane.The number of such configurations grows exponentially with the system size.This can be seen by noting that each of the irregular honeycombs has a "breathing" degree of freedom-it can contract and expand while keeping its aspect ratio.Such a move is local (shown in red), changing only the configuration in the vicinity of the breathing irregular honeycomb, making it a potentially useful experimental probe of the resonance between different trimer configurations.(d) This experimental probe can be explicitly constructed for a given breathing move as it is equivalent to a product of σx operators needed to expand the loop.An example of a simple breathing move is shown with configurations before and after the move differing only within the yellow shell around the perimeter of the loop.The moves needed to expand the honeycomb, in this case, lead to the many-body X-loop operator that is a product of 24 individual honeycomb σx operators.(e) The main experimental probe that definitively distinguishes between trivial and TQSL phases entails checking the two U(1) conservation laws[28] relating the electric field (arrows) flux with enclosed charges (number of different sublattice sites) along an arbitrary closed loop (example shown in blue).(f) Two Z-loop operators evaluated for the loop shown in blue in (e) for the 10 × 6 shape 1 cluster from Fig. 4(a) of the main text.The two U(1) conservation laws are violated in the trivial phase, but are approximately satisfied in both the TQSL and VBS phases.

FIG. 9 .
FIG. 9. Structure factors of the phases.Static structure factors for different regions of the phase diagram, as obtained by DMRG, are shown.The first Brillouin zone is denoted by a white hexagon.(a) The structure factor for the trivial disordered phase displays only a trivial peak at the Γ point.(b) The same is true in the TQSL regime (indicating the absence of order) along with additional broad features stemming from the trimer constraint being fulfilled.(c) For the Néel state-where the order is present within the honeycomb unit cell and does not break the translational symmetry of the lattice-sharp Bragg peaks are visible in the second Brillouin zone.In contrast, the columnar (d) and brick (e) phases break translational symmetry and present additional nontrivial ordering peaks.

FIG. 10 .
FIG. 10.Properties of the TQSL region in DMRG.(a) In contrast to the ordered phases, the TQSL region does not show a peak in von Neumann entanglement entropy and is characterized by a large entanglement entropy.(b)Going from trivial disordered phase to the TQSL region, one detects a strong decay of the gap as a function of ∆/Ω on the one hand, as well as a change in the behavior with increasing transverse system size.While in the trivial phase, the gap is approximately constant for different system sizes, in the TQSL region, the gap generally drops with increasing system size, which might hint at a gapless state in the thermodynamic limit.Inset: The second derivative of the gap (Ny = 5) also displays a peak akin to the energy and fidelity susceptibilities.(c) Similarly, on a plot of the column density fluctuations, a strong susceptibility towards fluctuations is seen on going from the trivial phase (∆/Ω < 1.9) to larger ∆/Ω.Although these fluctuations decay faster than the expected Friedel oscillations for a truly gapless state, the result is broadly consistent with a gapless state in the thermodynamic limit.

FIG. 12 .
FIG. 12. Exact-diagonalization signatures of transitions into the TQSL.Besides the signatures in TQSL overlap and fidelity susceptibility, transitions into the TQSL can be followed by (a) evolution of the low-lying energy spectrum where the character of the spectrum exhibits changes at both trivial-TQSL and TQSL-VBS transitions, and (b) the density or the energy susceptibility, which only detect the first transition, similar to Ref. 37.

FIG. 13 .
FIG. 13.Structure of the TQSL wavefunction.The weights and phases of different configurations in the ground-state wavefunction on two representative clusters.(a) An almost perfect equal-superposition of all MIS states describes the 6 × 10 shape 3 ground state in the TQSL region.The remaining weight is distributed among configurations with monomers with at least an order-of-magnitude smaller coefficients (see Fig. 4(c) in the main text).This is representative of all the shapes showing the perfect TQSL state in Fig. 11.(b) For the 10 × 6 shape 1 cluster in the TQSL region, the ground state is still well described by a superposition of all the MIS states with equal phases, but with a more pronounced weight asymmetry.(c) In contrast, deep in the VBS region of the 10 × 6 shape 1 cluster, the ground state is predominantly a superposition of only 29 out of 60 MIS states, with the next-highest contribution coming from the |MIS| − 1 subspace.

FIG. 14 .
FIG. 14.Effect of open boundaries.(a) The TQSL overlap for several representative clusters with cylindrical boundary conditions drops somewhat compared to the system on a torus, but the TQSL nature of the state is still preserved.(b) The density profile of the clusters in the TQSL state shows the primary effect of the open boundaries, namely, inducing the density fluctuations.(c) Columnar density fluctuations are enhanced in the TQSL region compared to the trivial state, leading to similar qualitative behavior as seen in the DMRG calculations [Fig.10(c)] with the equivalent boundary conditions.

FIG. 15 .
FIG. 15.Effect of interaction tails.The TQSL overlap as a function of ∆/Ω for several representative clusters with varying interactions, ranging from no tails (R b /a = 0) to a tail strength equivalent to 2 < R b /a < √ 7 in the k = 3 regime.The clusters that hosted a perfect TQSL state in the case without tails, (a), (b), and (d), show a sizable TQSL regime in the intermediate ∆/Ω region for all interaction strengths probed.The relatively less-robust cluster (c) now shows a significantly diminished TQSL region that completely vanishes at larger tail strengths.

FIG. 18 .
FIG. 18. Universality of the off-ramp projection mechanism.Spin liquid overlap during adiabatic preparation for several clusters and total protocol times in the case of (a) the honeycomb lattice without tails (TQSL fidelity), (b) the honeycomb lattice with tails (TQSL fidelity), and (c) the ruby lattice without tails for clusters explored in Ref. 56 (dimer RVB fidelity).In all the cases presented, significant fidelity gain is obtained in the off-ramp part of the protocol, which can be explained by the projection to the MIS subspace during the off ramp.In the case with interaction tails, after Ω drops below the appropriate value for the next blockade shell, projection to the smaller MIS subspace starts to occur, leading to a final upturn in infidelity.

FIG. 19 .
FIG.19.String degeneracy and the ground state in the second-neighbor blockade regime.(a) In the k = 2 regime, there is a classical sliding degeneracy of strings, arising because every two neighboring strings can be in parallel (top) or antiparallel (bottom) configurations.This degeneracy scales exponentially with the linear system size in contrast to the total system size scaling of the trimer degeneracy.(b) Overlap of the k = 2 PXP ground state with an equal superposition of string states ("string MIS") as a function of the detuning.While the string MIS fidelity is quite high for the system without tails, the state is not robust to the interaction tails.