Exotic quantum clusters and non-equilibrium dynamics of Rydberg excitations in one-dimensional optical lattices

In this mini-review, we report results from M. Mattioli, et al. [Phys. Rev. Lett. 111, 165302 (2013)], M. Dalmonte, et al. [Phys. Rev. B 92, 045106 (2015)] and M. Mattioli, et al. [New J. Phys. 17, 113039 (2015)], where it is shown that Rydberg atoms trapped in one-dimensional optical lattices are a useful tool to investigate the equilibrium phase diagram and the non-equilibrium dynamics of extended Hubbard models and Kinetically Constrained Models, respectively. Atoms weakly-dressed to an high-lying Rydberg state, which interact with a constant potential extended over several lattice sites, can be in an exotic quantum liquid state, the cluster Luttinger liquid phase [42, 43]. Furthermore, we show how a many-body model of interacting three-level atoms in the V-shaped configuration, where one of the level is a Rydberg state, might relax to equilibrium according to the same rules, so-called kinetic constraints, which are known to reproduce the characteristic dynamical arrest and separation of timescales of real glass-forming materials [62].


Introduction
Alkali Rydberg-excited atoms [1], whose relatively simple electronic structure resembles the one of hydrogen atoms, have exaggerated properties that depend on their principal quantum number n [2,3].
As an example, the van der Waals interaction between two Rydberg atoms is known to scale with the interatomic distance r as ∼ n 11 /r 6 [2]. While interactions between ground-state atoms are mainly weak and contact-like [4] and those between ions are non-adjustable even though large [5], Rydberg atoms possess both strong and tuneable interactions over distances of several μm. The tunability is twofold: on one hand, the sign of the interaction can be changed from positive (repulsive) to negative (attractive), and from isotropic (e.g. for s states) to anisotropic (for higher angular momentum states); on the other hand, in presence of external fields or close a e-mail: marco.mattioli@uibk.ac.at

2994
The European Physical Journal Special Topics to Förster resonances, it becomes dipolar-like, ∼ n 4 /r 3 [2], thereby enhancing its long-range character with respect to the van der Waals case.
The advance in atomic molecular and optical physics, that has nowadays given the possibility to trap and laser-excite Rydberg atoms in experiments, has sparkled intense theoretical and experimental activity with the aim to simulate condensed matter strongly-correlated systems. Indeed, the combination of the high degree of control typical of cold and ultracold atom experiments, together with the aforementioned strong, tuneable and long-range interactions, elect Rydberg atoms as ideal candidates to study spin models and novel many-body phases either at equilibrium or at non-equilibrium .
At present, the investigation of static (equilibrium) properties of spin models focuses on two main branches: In the so-called frozen gas regime [32], where kinetic energies and spontaneous emission rates from Rydberg states (∼ kHz) are much smaller than interaction energies (∼ GHz), and atomic motion is negligible within experimental timescales, an ensemble of Rydberg atoms can be mapped to a fully coherent many-body two-level spin Hamiltonian, whose two states are encoded in the ground-and Rydberg states. When, instead, kinetic energies are comparable to interactions (both being larger than spontaneous emission rates from Rydberg states) and atomic motion is relevant, Hubbard-like models can be investigated. For Rydberg atoms, the requirement is to go beyond the frozen-gas regime: in [33], it has been proposed to do so by weakly-admixing the Rydberg state to the ground-state using a far-off-resonant laser in order to obtain an interaction strength which is considerably smaller than the bare Rydberg interaction and comparable to average atomic kinetic energies (see introduction in Sect. 2). Furthermore, the resulting effective potential between dressed ground-state atoms, for distances smaller than a critical cut-off radius, remarkably flattens off and becomes constant. The latter feature is peculiar and does not occur typically in nature; nevertheless it gives the possibility to explore exotic models and states of matter which have been only theoretically investigated in literature. As an example, recent works [34][35][36][37] based on the dressing proposal, suggest that the ground-state of Rydberg-dressed atoms can be in a supersolid phase, an elusive quantum state of matter theoretically predicted almost 50 years ago [38][39][40] and characterised by the unusual coexistence of both crystalline and superfluid orders [41]. In one-dimensional (1d) systems, instead, as discussed in detail in Sect. 2 and [42,43], it has been proven that these so-called soft-shoulder potentials lead to the breakdown of the standard Luttinger liquid theory [44], known to capture the physics of a plethora of different models in one-dimension [45][46][47][48][49][50][51][52][53][54].
As anticipated above, also non-equilibrium aspects of spin models can be studied with Rydberg atoms. For example, recently it has been shown that dynamical constraints typical of certain idealized models, so-called Kinetically Constrained Models [55] (KCMs), appear naturally in the transient relaxation of ensembles of Rydberg atoms towards equilibrium. In KCMs, these kinetic constraints, which allow or forbid transitions among different configurational spin states depending on the state of neighboring spins, determine the dramatic increase of relaxation timescales and the emergence of dynamical heterogeneities [56,57] predicted and observed in real structural glasses [58][59][60][61]. In [62], whose main results are summarized in Sect. 3, dynamical constraints are encoded in the interplay between the Rydberg-blockade [63][64][65][66] and the laser detuning from the Rydberg state, in that the presence of a Rydberg excitation at a certain position facilitates (i.e. promotes) the creation of a further neighbouring Rydberg excitation if the laser detuning from the Rydberg state exactly compensates for the interaction energy shift between the two Rydberg excitations. Other proposals on how to encode similar kinetic constraints with Rydberg atoms can be found in [67][68][69][70][71][72]. We finally note that, recently, dynamical properties of a manybody Rydberg system exhibiting manifest kinetic constraints have been measured in [73], confirming the potential of cold Rydberg experiments to study phenomena, like glassy dynamics, typical of soft-condensed matter physics.
The remainder of this mini-review paper is organized as follows. In Sect. 2 we study the equilibrium properties of clusters of Rydberg excitations in 1d optical lattices. First we determine the frustrated underlying cluster structure in the classical limit (Sect. 2.1). In the quantum regime, we then study the strongly-interacting regime in perturbation theory (Sect. 2.2) and finally we establish an effective cluster Luttinger liquid theory (Sect. 2.3) whose predicitons, which differ from the standard Luttinger liquid scenario ones, are corroborated by numerical simulations presented in Sect. 2.4. In Sect. 3, we investigate the non-equilibrium dynamics of Rydberg excitations in 1d optical lattices. After introducing the model (Sect. 3.1), we review properties of a single three-level atom in the V-configuration (Sect. 3.2) and we then introduce (Sect. 3.3) and study numerically (Sect. 3.4) the dynamical behaviour of the many-body model, which we find to display signatures of glassy dynamics, including dynamical heterogeneities and separation of relaxation timescales.

Cluster Luttinger liquids of Rydberg-dressed atoms in optical lattices
In [42,43], we studied the phase diagram of an extended Bose-Hubbard model, described by the following Hamiltonian Here, b † i (b i ) are bosonic/fermionic creation (annihilation) operators at the site i (L being the total number of sites), n i = b † i b i is the number operator, J is the tunnelling rate and V is the soft-core off-site interaction potential which extends up to the critical cut-off radius R c .
These particular features of the interaction potential can be experimentally realised with Rydberg-dressed cold gases. In the weak-dressing regime, Ω |Δ|, atoms in their ground-states are off-resonantly coupled with Rabi frequency Ω and red detuning Δ < 0 to an high-lying Rydberg state. The resulting effective potential as a function of the relative distance x between pair of atoms reads [33] V where the cut-off radius is the Condon radius R c = [C 6 /(2|Δ|)] 1/6 and C 6 is the van der Waals coefficient of the addressed Rydberg state. At large distances x R c , V (x) reduces to the usual repulsive van-der-Waals interaction between Rydberg atoms ∝ x −6 , suppressed by a factor [Ω/(2Δ)] 4 , since only a small fraction [Ω/(2Δ)] 2 of the Rydberg state is admixed to the bare ground-state. However, for x < R c , a double Rydberg excitation is prevented by the dipole blockade and V (x) saturates to a universal constant value Ω 4 /8Δ 3 .

Classical cluster exchange model
The two relevant length scales in the classical limit [J = 0 in Eq. (1)] are R c and the average distance between particles r = 1/ρ 0 , with ρ 0 = N/L the particle density and 2996 The European Physical Journal Special Topics N the total number of particles in the system. This leads to three possible regimes: (i) liquid r > R c , (ii) crystal r = R c and (iii) cluster liquid r < R c . The groundstate, in the regime of our interest r < R c , is made of two different kind of building blocks, say A and B, characterized by a different number of particles and a fixed number R c of holes (i.e. lattice site where there is no atom). Furthermore, a given ρ 0 corresponds to a certain ratio of the number of blocks A and B, N A and N B , respectively. The classical ground-state consists of all possible permutations of the configurations (arrangements) between blocks A and B with the minimal energy. For open boundary conditions, it is possible to compute the degeneracy of the groundstate, which we define as the total number of minimal energy configurations where M = N A + N B = (L − N )/R c is the total number of clusters in the chain. Results exchange model rely on the fundamental assumption that all clusters are objects without an internal structure. In the following, we discuss the physical interpretation behind the aforementioned statement and its regime of validity. In Sect. 2.1.1, we show, that clusters are stable due to the presence of large energy barriers, which prevent hopping of particles to other clusters. In Sect. 2.1.2, we probe the energy landscape of the system and show that a series of local minima appear which are typical of clustered systems studied in soft-condensed matter physics [74][75][76][77][78][79].

Tunnelling energy barriers
As an example, let us consider the instructive case with R c = 4 and ρ 0 = 1/4, where an A-block is formed by two particles (followed by R c = 4 empty sites), and a B-block is formed by a single particle (followed by R c = 4 empty sites). As explained above, all configurations which can be written as a sequence of A-blocks and B-blocks, with the constraint N B = 2N A , are ground-states. We now identify, for later convenience, a cluster-particle as one of the two particle of the A-block. In the lowest-energy configuration, the cluster-particle sits precisely next the other particle, i.e. there are no holes in between them. Thus, in the ground-state, the displacement Δs (in units of lattice spacing a) of the cluster-particle from its ground-state position is, by definition, equal to 0. Due to the presence of finite-range interactions in our model, any Δs = 0 of the cluster-particle costs some energy. In our example, Δs = 4 corresponds to the exchange between an A-block with a neighboring B-block. We define the energy cost associated with this exchange process as tunneling energy barrier (where tunneling is here inspired by the corresponding quantum mechanical exchange process). Fig. 1(c) compares the tunnelling energy barrier of the box potential with the Rydberg-dressed potential from Eq. (2). We note that the barrier of the former is considerably larger than the one of the latter, suggesting that stronger interactions would be needed in the quantum regime to observe the new quantum cluster state in a realistic experiment. We quantify this feature by computing explicitly the barrier for both potentials. In the box case, the only energy penalty paid in one of the possible ground-state configurations is the pair interaction of the two particles in the A-block, E 0,Δs=0 = V . If instead Δs = 2, E 0,Δs=2 = 2V , as the cluster-particle interacts with two neighbors. The resulting energy barrier is thus ΔE 0 = E 0,Δs=2 −E 0,Δs=0 = V . In the case of the Rydberg-dressed potential, we have to consider all terms of pairwise interactions given by V (x) from Eq. (2). For the case of Fig. 1(b), the tunnelling The physical move associated to an A-with-B block exchange is the tunnelling of one cluster-particle (red sphere) from a doubly-occupied A-block (green) to a neighbouring singly-occupied B-block (red). Δs is defined as the displacement of a cluster-particle with respect to its ground-state original position, Δs = 0. Δs = 4 is again a ground-state configuration and corresponds to the exchange process completed. (b) For R c = 4, the ground-state and the maximally excited (Δs = 2) configurations are shown, including the energy penalties paid due to the presence of all other particles. Other excited states are Δs = 1 and 3 (not shown). (c) At Δs = 2, the tunnelling energy barrier of the box potential (black) is roughly a factor of 2 larger than the Rydberg-dressed barrier (red). Image taken from [42]. (This figure is subject to copyright protection and is not covered by a Creative Commons license.)

energy barrier is then ΔE
This is about a factor of two smaller than for the case of the box potential.
We conclude by remarking that the presence of these tunnelling energy barriers strongly suppresses fluctuations of the cluster-particle position for finite (but large) ratios of V /J in the quantum regime J = 0. This justifies to extend the assumption of no internal correlations between cluster particles of our classical exchange model to be valid also for the low-energy field theory we formulated in Sect. 2.3.

Energy landscapes of test particles
We conclude the discussion of classical energies by considering the energy landscape F (x), defined as the energy experienced by a test particle added at position x, while keeping all other particles in the original position. Such test particle will then 'sit' in correspondence of the local minima of F (x). For example, in Fig. 2, the test particle in the Rydberg-dressed potential will prefer to be next to singly-occupied B-blocks, precisely at positions x = 7, 10, 23, 26. For the box case, the flat shape of the interaction potential increases the total number of local minima of F (x), which are located at positions x = 5, 7, 8, 9, 10, 12, 21, 23, 24, 25, 26, 28. The complex structure of the classical ground-states results in the emergence of an exotic quantum liquid once quantum fluctuations are introduced. In Sect. 2.2, we will first present a strong-coupling approach in the V J limit, and in Sect. 2.3 a modified bosonization treatment which embodies the classical cluster constraints studied in Sect. 2.1. The combination of the two approaches allows us to gain a for the box potential (black), the discrete (red) and free-space (red, dashed) Rydberg-dressed potentials. Noticeably, all landscapes have similar barrier heights, i.e. h ryd 0.919 h box . However, the many-step-like structure of the Rydberg-dressed potential renders the cluster structures less stable with respect to the box case. We note that the physical interpretation of these barrier heights is different with respect to the tunnelling energy barrier of Fig. 1: here, F (x)/V denotes the energy cost of introducing a single test particle on top of a classical ground-state configuration, while ΔE 0/V in Fig. 1 describes the energy cost of moving a particle within the cluster ground-state manifold. Image taken from [42]. (This figure is subject to copyright protection and is not covered by a Creative Commons license.) qualitative analytical understanding of the cluster Luttinger liquid phase, which we quantitatively investigate in Sect. 2.4 with numerical simulations using the so-called Density-Matrix Renormalization Group algorithm (DMRG) [80,81].

Strong-coupling approach to cluster manifolds: the XX model
After having identified the manifold of clustered ground-states in Sect. 2.1, it is now possible to derive an effective Hamiltonian describing the strong-coupling limit J V of our extended Hubbard model. In order to do so, we define the projector G on the classical ground-state manifold. In second order perturbation theory, the Hamiltonian in Eq.
where G = 1 − F, 1 is the identity matrix, E 0 is the classical ground-state energy within the cluster manifold, H V is the the classical energy functional [interaction term in Eq. (1)] and H J is the perturbation Hamiltonian [hopping term in Eq. (1)]. First, we define effective spin-1/2 operatorsS j as follows: Ordering the cluster configuration with an index j ∈ [1, M], we associate to the position of each A-type cluster a spin-up, and for each B-type, a spin-down. This is a one-to-one mapping of the Hilbert space defined by F to the Hilbert space of a spin-1/2 chain with M sites. The effective Hamiltonian can be then cast in a compact form as a spin chain. In the illustrative example of R c = 2 and ρ 0 = 2/5, the second order effective Hamiltonian becomes whereS + j (S − j ) is the usual raising (lowering) operators associated toS j . The strongcoupling limit of Eq. (1) then reduces to a system of hard-core bosons (or spin-1/2) hopping in an artificial lattice generated by the underlying cluster structure. The cluster Luttinger liquid, which we will investigated in Sect. 2.3, interpreted as a Luttinger Liquid of clusters of particles [42,43], is adiabatically connected to the strong-coupling Hamiltonian (5).

Phenomenological cluster Luttinger liquid theory: beyond the standard Luttinger liquid theory
Here we present the modified Luttinger liquid theory [44], which we call cluster Luttinger liquid Theory, as it is based on the underlying classical cluster structure. Even though our model is a discrete one, let us start with the most general form of an interacting 1d Hamiltonian in the continuum where is the bosonic quantum field at position x, the first integral represents the kinetic term and the second (double) integral is the potential contribution, with V (x − y) the (translational invariant) pairwise interaction between two particles at position x and y. We will verify that the discrete nature of our model does not affect significantly the predictions of the theory we developed.
In the density-phase representation where we have introduced the phase field θ(x), the following commutation relation holds Let us now introduce a new field, Π(x), describing the oscillation of ρ(x) on top of the mean density ρ 0 , such that The oscillation field and the phase field are conjugate of each other since Up to know this is the standard hydrodynamic approach [44]. In order to go beyond it, as in the traditional Haldane's treatment, we introduce the auxiliary field ϕ(x) such that The auxiliary field is assumed to increase by 2π every time it encounters a particle while scanning the system in all its length L, i.e.
[see Fig. 3(a)]. If we assume an homogeneous distribution of particles, it is possible to rewrite the density operator as

3000
The European Physical Journal Special Topics = Fig. 3. Configurations allowed for the field ϕ(x). (a) In the Luttinger liquid scenario, ϕ(x) describes all possible particle configurations and takes integer values at the position of each particle (red circles). The integrated particle density (thin orange line) follows the profile of ϕ(x). (b) In the cluster Luttinger liquid scenario, ϕ cl (x) is constrained by the cluster structure, in which a certain number of particles are packed next to each other (blue circles). The integrated particle density (thick orange line) does not follow the behaviour of ϕ cl (x), but jumps whenever a two-particle cluster is encountered. Image taken from [43]. (This figure is subject to copyright protection and is not covered by a Creative Commons license.) the quantity in brackets indicating that particles are located at points x where ϕ(x ) = 2πn. Finally, using the Poisson summation formula, (with A a constant) and the fact that x > 0, it is possible to show that where, without loss of generality, we have redefined the summation index as p = −k/2 and Let us now consider a more complicated shape of the initial particle distribution which reflects the underlying classical cluster structure: Here f (x m ), which satisfies M m=1 f (x m ) = N , can be equal to either 1 or 2. Due to our choice of the density of particles, this means that we can have clusters of either 1 or 2 particles. In the center equality, ∈cl represents the sum over all two-particle clusters. The key point is that the first sum indicates that there are M = (L − N )/R c clusters, each one containing at least one particle, while the second sum takes into account the possibility of having clusters formed by two particles. We now introduce a new field, ϕ cl (x), schematically depicted in Fig. 3(b), which accounts for quantum fluctuations on top of the classical cluster density such that in analogy to ϕ(x) in the standard Haldane scenario. Using we can rewrite the ansatz for the density field as Next we need to Fourier transform Eq. (20). In order to do so we expand its sum as where, in analogy to Eq. (17), ∈cl is over all two-particle clusters. The first sum can be easily Fourier transformed exploiting again the Poisson summation formula Now we need to Fourier transform ∈cl , which gives as a result just renormalized coefficients in the previous expression, thus not affecting its functional form. We notice that they are nevertheless expected to be very small for M N , that is when the number of two-particle clusters is not too large.
In analogy to Eq. (16), we rescale the field as where σ = M/N and then we can express the density as Here the prefactor σ has been introduced in order to compensate for the delta function renormalized coefficients a p . The single particle bosonic operators ψ(x) [see Eq. (7)] can also be computed applying an analogous procedure to get where θ cl (x) is the conjugate field of ϕ cl (x), i.e.
in analogy to Eq. (10). This can be easily verified considering the role of the factor β in Eq. (25): if we approximate the density operators with its non-oscillating part, we need, from the analogous of Eq. (8), the following commutation to be fulfilled

3002
The European Physical Journal Special Topics Equation (27) is satisfied if β = σ, which then automatically ensures that Eq. (26) holds.
The effective low-energy Hamiltonian associated to the fields ϕ cl (x) and θ cl (x) [or, equivalently, ϕ cl (x) and θ cl (x)] is a compactified boson theory ( = 1) where v cl is the (cluster) sound velocity while K cl is the (cluster) Luttinger parameter. The Hamiltonian in Eq. (28) is gapless and conformal, with central charge c = 1. All correlation functions of the microscopic operators can be evaluated using conventional techniques. This theory is very similar to a Luttinger liquid, with the exception that the mapping between microscopic and low-energy degrees of freedom presents remarkable differences, which are visible in the correlation functions we are going to study in Sect. 2.4.

Numerical results
Here we investigate the difference between a standard Luttinger liquid and a cluster Luttinger liquid comparing the analytical results of our (modified) cluster Luttinger liquid theory 2.3 with DMRG numerical simulations [80,81]. The modified functional form of the density operator [Eq. (24)] allows us to make predictions on the scaling of the density-density correlation functions with α 1 , α 2 and γ 1 non-universal coefficients. Equation (29) displays a different spatial modulation with respect to the standard Luttinger liquid correlation function where K is the so-called Luttinger parameter and α is a constant. In fact, in Eq. (30), the periodicity of the oscillation is determined by the particle density 2πρ 0 , while in the cluster Luttinger liquid it is determined by the cluster density, 2πρ 0 σ. This feature is reflected also in the static structure factor S(q) = 1 + ρ 0 L 0 dx e iqx ρ(x)ρ(0) , i.e. the Fourier transform of the density-density correlation function. In the extended Hubbard model of Eq. (1), assuming periodic boundary conditions (PBCs), S(q) can be written as where both indices , have to be summed over all lattice sites L. The factor 1/L ensures that S(q) is finite also in the themordinamic limit. Peaks in S(q) are not located at momenta associated with the particle density, as in standard Luttinger liquids. For example, the lowest-momentum peak of the cluster Luttinger liquid structure factor, as it can be verified in Fig. 4(a), is located at The bosonic field operator we calculated in Eq. (25) allows to predict the behaviour of the single-particle Green's function, ψ † (x)ψ(0) , which also signals the departure In both panels plots in the range [π, 2π] are not shown because, due to PBCs, they are symmetric with respect to [0, π]. Image adapted from [42]. (This figure is subject to copyright protection and is not covered by a Creative Commons license.) from the standard Luttinger liquid scenario. In particular, its Fourier transform, the momentum distribution m(q) = L 0 dx e iqx ψ † (x)ψ(0) or, in our model, exhibits a peak at the same critical momentum q of S(q) (see Fig. 4(b) and [42]). We note that other signatures of the breakdown of the standard Luttinger Liquid theory can also be found using level spectroscopy techniques [82].

Glassy dynamics of Rydberg excitations in optical lattices
After dealing with static equilibrium properties of 1d clusters of Rydberg excitations, here we study the dynamical properties of the relaxation of Rydberg excitations in 1d optical lattices that 'hop' from one site to another according to the same kinetic constraints of KCMs, the idealized lattice models known to capture dynamical slowing down and space-time heterogeneities of structural glasses.

The model
We schematize each atom, trapped at its site in a largely spaced 1d optical lattice [24,[83][84][85][86][87][88], with lattice unit a, as a V-shaped three-level system. Atoms are resonantly driven with Rabi frequency Ω e from the electronic ground-state |g to a low-lying excited state |e , which is short-lived (i.e. its associated spontaneous emission rate γ e is large), and far-off-resonantly driven to a metastable Rydberg state |r with detuning Δ r and Rabi frequency Ω r (see Fig. 5).
Tuning Ω e allows the exploration of the crossover from classical to coherent dynamics of Rydberg excitations trough the 1d chain. More precisely, for large Ω e , the relaxation of Rydberg excitations resembles the non-equilibrium dynamics of a classical KCM, the so-called One-Spin Facilitated Model [55,69] (1-SFM). As already mentioned, kinetic constraints in our model are a consequence of the interplay between Rydberg blockade [63][64][65][66] and the laser detuning Δ r . In contrast, for vanishing Ω e , coherences of the |g ↔ |r transition dominate frequencies Ω e (Ωr) between |g -|e (|g -|r ) and γe is the decay rate due to spontaneous emission from |e . The anti-blockade condition V = Δ r , i.e. when the laser detuning Δr with respect to the |g -|r atomic transition is equal to the interactions between neighbouring Rydberg-excited atoms V , ensures that whenever an atom is in |r (defect, blue sphere), an atom next to it, initially in |g (facilitated atom, red sphere), will also be favoured to be excited to |r . An atom in |g far from a defect (non-facilitated atom, black sphere) is blocked, because its |g -|r transition is far-off-resonant with large Δ r . Image adapted from [62]. (This figure is subject to copyright protection and is not covered by a Creative Commons license.) the dynamics within the lifetime of the addressed Rydberg state, given by the inverse of its spontaneous emission rate γ −1 r [24]. Let us now make the following assumptions, which will be valid later on, unless otherwise stated: the Rydberg state is metastable, that is γ r = 0; interactions between Rydberg-excited atoms are only nearest-neighbor. Since the real potential between two Rydberg excitations is of the van der Waals type [2], V (x) = V /x 6 , the approximation translates to V (1) = V and V (x 2) = 0 (x here is in units of a); -Ω r Ω 2 e /γ e , for reasons that will be clear later on.
The Hamiltonian of the system described above, in the rotating frame (with = 1) and after performing the rotating wave approximation, reads

A single-atom in the V-configuration
Before studying the many-body model, we review here the properties of a single three-level atom in the V-configuration, which has been already extensively studied in literature, both theoretically [89][90][91][92][93][94][95] and experimentally [96,97]. For a single atom the Hamiltonian in Eq. (34), reduces to In the Born-Markov approximation, the master equation describing the dissipative evolution of the reduced system density matrix (t) is in the Lindblad form [98] ( where c = |g e| is the quantum jump operator, which (ideally) projects instantaneously the atomic state into the ground-state. Let us imagine to measure photons emitted from |e with a perfectly efficient photodetector. The theory of continuous measurement [95] ensures that, under continuous monitoring of the strong transition |g ↔ |e with our photodetector and during time intervals where no photons are detected, the system wave function |ψ(t) evolves according to the Schrödinger equation i d dt |ψ(t) = H eff |ψ(t) , with an effective non-Hermitian Hamiltonian and H given by Eq. (35). The general solution for the system wave function is where λ n and |u n are the eigenvalues and associated eigenvectors of H eff , respectively [99]. The coefficients c n can be determined from the initial condition, e.g. |ψ(0) = |g . In the limit Ω r Ω e , γ e , Δ r , λ n and |u n can be computed in perturbation theory.

Brillouin-Wigner perturbation theory
In the following, we summarize the steps for the calculation of λ n and |u n of H eff in Eq. (37) within a Brillouin-Wigner perturbation theory approach. We compute explicitly λ 3 and |u 3 (the other eigenvalues and eigenvectors can be in analogy straightforwardly computed). In the {|g , |e , |r } basis, H eff has the following matrix form Assuming Ω r Ω e , γ e , Δ r , one can split Eq. (39) in  which diagonalizes H 0 , such that where U † is the Hermitian conjugate of U , E m = i 4 −γ e + γ 2 e − 4Ω 2 e and E p = − i 4 γ e + γ 2 e − 4Ω 2 e . Applying the same transformation on H 1 ,H 1 = UH 1 U † , we getH where Ω c = Ω r cosh { 1 2 arctanh 2Ωe γe } and Ω s = Ω r sinh { 1 2 arctanh 2Ωe γe }. The unitary transformation U defines a change of basis from {|g , |e , |r } to, say, We now define two operators P and Q which project onto |r and its complementary subspace, respectively. Since U does not affect the |r subspace [see Eq.
while Q can be easily obtained from the orthogonality property of projectors, i.e. P + Q = 1, where 1 is the identity operator. The first-(E 1 ) and second-order (E 2 ) corrections to the eigenvalue E 0 = [PH 0 P] 3,3 = [PH 0 P] 3,3 = −Δ r of the transformed unperturbed HamiltonianH 0 are and Here, [A] i,j indicates the ith row and jth column of the matrix A and R 0 = (E 0 1 − QH 0 Q) −1 is the so-called resolvent. Up to second order in Ω r , λ 3 can be approximated as We can now calculate the first-order correction to |r originated by the other two eigenstates ofH 0 , |− and |+ , both mixed to |r by the transformed perturbation HamiltonianH 1 . For example, the eigenvector |u 3 associated to λ 3 , is: As mentioned above, the coefficients c n can be calculated from the initial condition. For example, setting |ψ(0) = |g , we get

Quantum trajectories of a V-shaped three-level atom
The eigenvalues λ n are complex numbers with negative imaginary parts Im{λ n }. As a consequence, |ψ(t) does not conserve the norm during the evolution. For Ω r Ω 2 e /γ e , one of the eigenvalues of H eff , exactly λ 3 in Eq. (49), has a much less negative imaginary part with respect to the other two, |Im{λ 3 }| |Im{λ 1 }|, |Im{λ 2 }| [99]. This determines long tails in the delay function which is the conditional probability that no photon has been detected until time t, given a photon count has occurred at t = 0. Note that D 0 (t) is exactly the norm of |ψ(t) , and its monotonous decay is dominated, in the long time limit, by Im{λ 3 }. Thus Im{λ 3 } is responsible of the existence of a small but finite probability of a long time window with no photon counts. When such a long interval of no photon counts occurs, the atom is prepared in the Rydberg state as |ψ(t) is dominated by |u 3 [Eq. (50)]. The theory of continuous measurement teaches us that the master equation (36) can be obtained simulating several stochastically different system realizations (or trajectories) with wavefunction |ψ α (t) and then averaging over them. The index α labels each of the N traj independent trajectories, which can be reproduced according to the following algorithm [91]: 1. Prepare the system in a certain (normalized) state |ψ(t 0 ) α at the initial time t 0 . 2. Compute the probability of the occurrence of a quantum jump at t 0 from p α (t 0 ) = γ e ψ (t 0 )|c † c|ψ(t 0 ) α dt 1, where dt is an ideally infinitesimal time step. The corresponding probability of no quantum jump is 1 − p α (t 0 ) ≈ 1. 3. According to these probabilities, randomly choose a no jump/jump event. 4. In case of no jump, propagate |ψ(t 0 ) α with H eff [Eq. (39)] to obtain the normalized wavefunction |ψ(t 0 + dt) α = e −iH eff dt |ψ(t 0 ) α /||...|| = (1 − iH eff dt)|ψ(t 0 ) α /||...||. Otherwise, project the wave function into the ground-state, |ψ(t 0 + dt) α = √ γ e c|ψ(t 0 ) α /||...|| = |g . 5. Repeat the procedure above for each time t ∈ [t 0 , t max ] until the desired final trajectory time t max is reached.
A typical experimental sequence of photon detections, that is of quantum jumps in the trajectory language, looks like an alternation of bright and dark periods. During a bright period, labelled by 0, the atomic population is rapidly cycled between |g and |e and the time between two successive photon detections is of the order of the lifetime of |e , γ −1 e . This means that the photodetector will count a large number of photons in a bright period. In contrast, during a dark period, labelled with 1, the atomic population is shelved in |r and no photons are detected. An approximation of the exact time evolution of (t) from Eq. (36) can be obtained simulating N traj different trajectories and finally averaging over them, In the following, the observable we are interested in is the Rydberg population, that is, the expectation value of the Rydberg projector |r r| in |ψ(t) α , Averaging r α (t) over N traj leads to In analogy to Eq. (53), lim Ntraj→∞ R(t) = rr (t), where rr (t) ≡ r| (t)|r is the Rydberg population whose evolution is governed by the single-atom master Eq. (36). As expected, Fig. 6 shows that R(t) converges to rr (t) upon increasing N traj .

Calculation of single-atom rates
In the following we illustrate three independent methods for the estimation of the average duration of bright and dark periods T 0 and T 1 , respectively, or of their inverse, the transition probabilities (rates) from 0 to 1, Γ 0→1 ≡ T −1 0 , and from 1 to 0, Γ 1→0 ≡ T −1 1 . In [93], these rates are found to be together with rr (0) = 0, and together with rr (0) = 1. The rates in Eqs. (56) and (57) can thus be evaluated fitting the slope of rr (t) within the initial times of the evolution and the proper initial condition (e.g., from Fig. 6, it is possible to estimate Γ 1→0 since the initial condition is rr (0) = 1). An analytical derivation of the transition rates is also possible. In the long time limit the delay function D 0 (t) ∼ p e −t/T1 , where p is the (small) probability to be in a dark period [94,95]. Comparing it with Eq. (52) one gets and where the results are obtained substituting λ 3 from Eq. (49) and c 3 from Eq. (51). In order to compute Γ 0→1 (Δ r ), we note that, during a bright period, |r is almost never populated and the effective emission rate will be that of the (resonantly driven) two-level system |g ↔ |e with decay γ e , i.e. γ e Ω 2 e /(γ 2 e + 2Ω 2 e ). Then, Γ 0→1 (Δ r ) is simply this two-level emission rate multiplied by p [99], If Δ r = 0 the two rates simplify to Finally, the third independent method to compute the rates relies on the direct numerical calculation of the average over trajectories of the time-averaged length of  In all three methods, parameters are: Ω e = γe and Ωr = 0.03γe. Within the QT method, simulations of N traj = 6000 trajectories up to times γet = 10 6 have been performed.

Many atoms in the V-configuration
We now switch to the many-body model. In analogy to the single-atom case [Eq. (36)], the many-body master equatioṅ with c i = |g e| i the quantum jump operator at site i, describes the driven-dissipative dynamics with the coherent evolution given by the Hamiltonian in Eq. (34). Furthermore, let us assume: (i) Ω r Ω 2 e /γ e , as in Sect. 3.2, for well defined quantum jumps to occur, (ii) Δ r Ω e , Ω r , γ e and (iii) V = Δ r , which we call the anti-blockade condition. The many-body transition rates, analogous to the single-particle rates in Eqs. (58) and (60), can be determined as follows. In a 1d system, due to our assumption V (x 2) = 0, every process at a certain site i can be influenced at most by its left and right neighbors, i ± 1. The many-body rates can be calculated from Eqs. (58) and (60) (58) and (60), Δ r appears only in even powers, both off-resonant rates = 0, 2 are equally suppressed with respect to the resonant rate = 1, as highlighted in Table 2.
In absence of conditions (i), (ii) and (iii), the strong transition |g ↔ |e might significantly increase the probability Γ 0(1)→1(0) (±Δ r ) of (originally) suppressed processes of excitation (deexcitation) to (from) |r to take place. Since rates differ by several order of magnitudes (see Table 2), the dynamics shows a large separation of timescales (cf. Sect. 3.4).
Allowed ( = 1) and suppressed ( = 0, 2) processes are analogous to those who characterize Spin Facilitated Models, as discussed in the following.

Spin facilitated models
Let us now introduce here the basics of the dynamical behaviour of the 1-SFM and the analogy between its many-body transition rates and the rates of our model. Spin Facilitated Models (SFMs) are a class of KCMs that exhibit dynamical heterogeneities and slowing down of relaxation timescales [55]. Even though they can be defined in a variety of lattice geometries and dimensions, from now on we will limit ourselves to 1d. Let us now consider the energy function of non-interacting spins in an external field, introduced by Fredrickson and Andersen in [100,101], where n i = 0, 1 is a classical spin variable. In our setup, n i = 0 (1) corresponds to the ith atom being in a bright (dark) period. In the SFM language, n i = 1 states are called defects, in that their number is smaller than n i = 0 states. A system described by Eq. (63) does not feature a finite temperature phase transition for a field strength K = 0. The equilibrium concentration of defects, reads where β = 1/(k B T ), k B is the Boltzmann's constant and T is the equilibrium temperature of the system. Furthermore, for every K > 0, lim T →0 d eq = 0, in line with the intuition that at low T the concentration of defects should be small. From now on, without loss of generality, we set k B = 1 and we fix the temperature scale setting K = 1.
Despite the simple ground-state of Eq. (63), the non-equilibrium dynamics towards it can be non-trivial. In general, in 1d, the evolution of a classical system of spins is described by a rate equation for the probabilities p(n, t) of being in the many-body spin state specified by the vector n = (n 1 , ..., n i , ..., n N ) at time t: Here, Γ(m → n) is the rate of the transition from m to n = m. The equilibrium probabilities π(m, t) satisfy the detailed balance Γ(m → n) π(m, t) = Γ(n → m) π(n, t) with respect to Eq. (63). Under the assumption of local spin-flips (i.e., at every timestep only one spin state at a certain site of the chain can change) and Glauber dynamics [102], the transition rate from (n 1 , ..., n i , ..., n N ) to (n 1 , ..., 1 − n i , ..., n N ) at site i, Γ ni→1−ni , is equal to 1/(1 + e βΔE ), with ΔE = 1 − 2n i the energy difference between the final and initial spin configurations. Since, for local spin flips, ΔE = ±1 are the only allowed outcomes, from Eq. (64) the single-particle rates can be written as Thus, the asymmetry between the ratio of the rates of creation over destruction of a defect fixes the temperature T at which the system of non-interacting spins should thermalize in the long-time limit. (Note that Γ 0→1 → Γ 1→0 for T → ∞.) Up to now, all transitions between different configurations are allowed. The key idea behind SFMs (and KCMs in general) is to introduce restrictions in these transitions, which are exactly the kinetic constraint we were talking about in Sect. 1. The general rule can be stated as: a spin flip is allowed only if there are enough defects in the neighbourhood that can facilitate such process. For the so-called 1d 1-SFM, the dynamical rule requires that a spin can flip provided at least one of its two neighbors is a defect. This facilitation condition transforms Eq. (66) into The association between rates in Eq. (68) and the many-body rates of our model is straightforward: n i−1 = n i+1 = 0 corresponds to = 0 for which no spin flip at site i is allowed, while either {n i−1 = 1, n i+1 = 0} or {n i−1 = 0, n i+1 = 1} corresponds to = 1 and i is facilitated to flip its spin. The only difference arises in the case of {n i−1 = n i+1 = 1}: while rates in the 1-SFM are doubled with respect to the case {n i−1 = 1, n i+1 = 0}, for = 2, rates are suppressed as much as in the case of = 0. This will lead to different steady states which we will discuss in Sect. 3.4.

Numerical results
In the limit of strong driving to the intermediate state Ω e ∼ γ e , the dynamics of the system is dominated by dissipation. As a method to study the system in the presence of both coherent driving and dissipation, we use quantum trajectory simulations (cf. 3.2). Figure 7 depicts a typical trajectory, labelled by α, of the Rydberg population r α (t) from Eq. (54). While in the non-interacting case, Δ r = V = 0 [panel (a)], the dynamics of Rydberg excitations is unconstrained and atoms behave independently from each other, as expected, in the strongly-interacting case, Δ r = V = 0 [panel (b)], clear spatial correlations emerge and the dynamics of Rydberg excitations resembles the ones of defects in the 1d 1-SFM [55]. The dynamics of defects in the 1d 1-SFM has been shown to be mainly diffusive [103] and it can be explained introducing three different species of atoms, which we label, in analogy to the SFM language, as defects, that is atoms in |r ; -facilitated atoms, that is atoms in either |g or |e with one defect as nearestneighbor; -non-facilitated atoms, that is atoms in either |g or |e with no defect as nearest-neighbor.
Let us illustrate the diffusion process of defects in more detail. A typical trajectory is depicted in Fig. 7(b) and marked by a blue ellipses. A defect can facilitate the excitation of a neighboring atom. In this example, the atom to the right is excited. For a short time, both are in the excited state but soon the original atom is deexcited (soon meaning that the probability of de-excitation is much larger than the probability of excitation): therefore, the defect has effectively moved to the right by one site. Another possible move is highlighted in Fig. 7(b) with a red ellipses. In this second example, the left neighbor is excited but later on de-excited. Here, the defect did effectively not move. The diffusion process can also be defined more quantitatively in some limits. For example, in the regime T 1 (or, equivalently, d eq 1), Γ 0→1 (0)/Γ 1→0 (0) ∼ d eq ∼ e −β and the diffusion constant is d eq /2 [55]. However, in the intermediate regime that characterizes our work, the diffusion constant does not have such a simple analytical expression, since Γ 0→1 (0)/Γ 1→0 (0) = 0.33 is not much smaller than 1.

Separation of timescales of defects, facilitated and non-facilitated atoms
In the following, we study the relaxation of the three different species individually. Numerical simulations of N traj stochastically different trajectories are averaged in order to approximate the dynamics of Rydberg excitations, according to the continuous measurement theory. We plot R i (t) from Eq. (55) for each atom i as a function of time in Fig. 8. Facilitated atoms (red lines) show a dynamical behaviour analogous to the one of resonantly-driven single-atoms, due to the anti-blockade condition which acts as a kinetic constraint in the 1-SFM [cf. Eq. (68)]. In fact, the inset of Fig. 8 shows that, in proximity of t = 0, the profiles of facilitated atoms are very well approximated by the resonant (Δ r = 0) single-atom Rydberg population rr (t) (see 3.2), the associated transition rate being Γ 0→1 (0).
Contrarily, the dynamics of defects (blue lines) and non-facilitated atoms (black lines) is slowed down with respect to the one of facilitated atoms. At t = 0, the weak transition |g ↔ |r of a defect is far-off-resonant and Γ 1→0 (Δ r ) is orders of magnitude suppressed with respect to Γ 1→0 (0) [cf. Table 2]. A defect can in fact relax only after one of its neighboring facilitated atoms is excited to the Rydberg state. Similarly, at t = 0, the weak transition |g ↔ |r of a non-facilitated atom is far-off-resonant and Γ 0→1 (Δ r ) is orders of magnitude suppressed with respect to Γ 0→1 (0). The inset of Fig. 8 shows that profiles of non-facilitated atoms agree with the far-off resonant (Δ r = 0) single-atom Rydberg population rr (t): non-facilitated atoms, in order to be excited to the Rydberg state, have to wait for one of their two neighbors to be Rydberg-excited too. Note that, for larger N , we predict the existence of non-facilitated atoms with profiles of R i (t) even slower than the ones in Fig. 8. This can be regarded as a form of hierarchical relaxation [71], in that the more distant a non-facilitated atom is with respect to a defect at t = 0, the more slowed-down the evolution of R i (t) will be. Hierarchical relaxation is a feature typical of classical SFMs and more in general a signature of glassy dynamics predicted in soft-condensed matter systems [104]. Off-resonant rr (t) with rr (0) = 1 well approximates also the short time dynamics of defects (not shown). Image taken from [62]. (This figure is subject to copyright protection and is not covered by a Creative Commons license.) Moreover, from Fig. 8 it is possible to infer that the steady-state value r of the average concentration of defects, does not coincide with the expected (thermal) concentration of Rydberg excitations in the 1-SFM d eq = γ 2 e 2(γ 2 e + Ω 2 e ) .
Our system, thus, does not thermalize in the long-time limit, but it seems to converge to a different steady-state. The reason has to be found in the presence of avoided merging processes, highlighted in Fig. 7(b) with a yellow ellipses. Two Rydberg excitations, separated by one lattice site in the ground-state, cannot merge, in that the atom in the middle cannot be Rydberg-excited since it feels a detuning 2V . This destroys the facilitating property of the anti-blockade condition V = Δ r . Let us finally mention that, in theory, one could think of workaround the problem of avoided merging processes and reproduce the exact dynamics of the 1d 1-SFM introducing a third laser which couples weakly |g ↔ |r with the same small Rabi frequency Ω r and with large detuning 2Δ r from |r (in addition to the original laser with Rabi frequency Ω r and detuning Δ r ). In this configuration, processes of merging of Rydberg excitations would be allowed and we expect the steady-state to be also the thermal state of the 1-SFM, i.e. r = d eq .

Conclusions
We have studied the equilibrium and non-equilibrium dynamics of clusters of Rydberg excitations in one-dimension. We have demonstrated the stabilization of cluster quantum liquids with qualitatively new features with respect to the standard Luttinger liquid scenario. These features, which result from frustration and cluster formation in the corresponding classical ground-state structure, can be readily observed in stateof-the-art experiment with Rydberg-dressed atoms in optical lattices.
In addition, we have proposed a model system of interacting V-shaped three-level atoms, where one of the level is a Rydberg state, and studied the resulting glasslike dynamics of Rydberg excitations. We found that the anti-blockade regime, that is where the nearest-neighbour interaction among Rydberg-exicted atoms equals the single atom laser detuning from the Rydberg state, implements a form of kinetic constraint, causing the dynamical arrest and increase of relaxation timescales typical of glass-forming systems.