Engineering non-binary Rydberg interactions via phonons in an optical lattice

Coupling electronic and vibrational degrees of freedom of Rydberg atoms held in optical tweezer arrays offers a flexible mechanism for creating and controlling atom-atom interactions. We find that the state-dependent coupling between Rydberg atoms and local oscillator modes gives rise to two- and three-body interactions which are controllable through the strength of the local confinement. This approach even permits the cancellation of two-body terms such that three-body interactions become dominant. We analyze the structure of these interactions on two-dimensional bipartite lattice geometries and explore the impact of three-body interactions on system ground state on a square lattice. Focusing specifically on a system of $ ^{87} $Rb atoms, we show that the effects of the multi-body interactions can be maximized via a tailored dressed potential within a trapping frequency range of the order of a few hundred kHz and for temperatures corresponding to a $>90\% $ occupation of the atomic vibrational ground state. These parameters, as well as the multi-body induced time scales, are compatible with state-of-the-art arrays of optical tweezers. Our work shows a highly versatile handle for engineering multi-body interactions of quantum many-body systems in most recent manifestations on Rydberg lattice quantum simulators.

Introduction.-In the past years Rydberg atoms [1][2][3] held in optical tweezer arrays have emerged as a new platform for the implementation of quantum simulators and, potentially, also quantum computers [4][5][6][7][8][9][10].One- [6], two- [11] and threedimensional [12] arrays containing hundreds of qubits are in principle achievable and the wide tunability of Rydberg atoms grants high flexibility for the implementation of a whole host of quantum many-body spin models.The physical dynamics of these quantum simulators takes place in the electronic degrees of freedom which mimic a (fictitious) spin particle.Effective magnetic fields and interactions are achieved via light-shifts effectuated by external laser fields and the electrostatic dipolar interaction between Rydberg states.Additional tuning with electric [13] and magnetic fields [14] permits the realization of exotic interactions, allowing for the study of ring-exchange Hamiltonians [15][16][17][18], frustrated-spin models [19][20][21] or crystallization phenomena [22][23][24].Within this context, in the last decade systems with tunable two-and three-body interactions [25][26][27][28][29] have attracted a lot of attention since the latter are responsible for the emergence of many exotic quantum states of matter, ranging from topological phases [30,31] to spin liquids [32,33].
In this work we put forward a new mechanism for engineering non-binary interactions in Rydberg tweezer arrays [6,9,[34][35][36][37][38][39][40][41][42][43].Here, each atom is held in place by a strong local harmonic potential.The simultaneous excitation of neighboring atoms to the Rydberg state gives rise to a mechanical force that couples the electronic degrees of freedom to the local phonon modes.We show that this coupling gives rise to effective spin-spin interactions between excited atoms.Similar mechanisms in which effective inter-particle interactions arise as a consequence of the coupling with an extra degree of freedom have been extensively studied in condensed matter systems.Here, well-known examples include the electron-electron interaction mediated by lattice phonons in metals [44] and the indirect spin-spin couplings [45] due to the Ruderman-Kittel-Kasuya-Yosida [46], superexchange [47], and Dzyaloshinskii-Moriya mechanisms [48].In these cases, integrating out the extra degree of freedom typically results in two-body effective interactions between the remaining degrees of freedom.Crucially, in our system, since spins and phonons are coupled via pairs of Rydberg atoms, not only two-body but also three-body effective interactions arise.We analyze in depth the interplay between the various effective couplings in the case of two-dimensional (2D) bipartite lattice geometries, demonstrating that regimes dominated by three-body interactions can be achieved.Our results show that the multibody interactions arising from the electron-phonon coupling are highly tunable and can drive non-trivial phase transitions in the ground state of a Rydberg spin system.By tuning the local harmonic potentials, we show that checkerboard, striped, and clustered phases occur as well as signatures of frustration phenomena.Our work is directly relevant for recent developments on the domain of quantum simulation with Rydberg tweezer arrays where it highlights a so far unanticipated mechanism for experimentally realizing exotic interactions.2D model.-Weconsider a 2D lattice of N Rydberg atoms in the x − y plane, whose sites are labeled by k = (k x , k y ).The electronic degree of freedom is modeled as effective twolevel system (with |↓ and |↑ denoting the ground state and the Rydberg excited state, respectively) [3,40].The two levels are coupled by a laser with Rabi frequency Ω and detuning ∆ [see Fig. 1(a)].Each of the atoms, with mass m, is trapped in a strong three-dimensional harmonic potential, characterized by trapping frequencies ω µ along the directions µ = x, y, z.The atomic motion inside the confining potential can then be described in terms of the bosonic operators b k,µ .The FIG. 1. Setup.(a) Each atom is modeled as a two-level system with ground state |g and excited Rydberg state |r .The two levels are coupled by a laser with Rabi frequency Ω and detuning ∆.The atom is trapped inside a tight harmonic optical tweezer (grey) and, at low temperature, it occupies the ground state of the associated phonon degree of freedom (red).For simplicity, we assume that Rydberg and ground state experience the same trapping potential.(b) Energy diagram of a two atom system arranged along the x−axis.When both atoms are excited to the Rydberg state, |rr , they experience, in addition to the electronic dipolar interaction V, the potential change δV arising as a consequence of the coupling between spin and phonon degrees of freedom and consisting of both two-and three-body contributions; see text for details.This also results in a state-dependent displacement δx 1,2 of the atoms from their equilibrium position x 0 1,2 , separated by the lattice spacing a.
Hamiltonian describing the single-particle dynamics is Here, where the prime in the sum implies that terms with equal indices are excluded.Note that in Eq. ( 1) we have assumed the same trapping frequencies ω µ for atoms in the ground and in the Rydberg state.This "magic" condition can be realized in Rydberg tweezer arrays through bottle beam traps [49].Furthermore, a small frequency mismatch between the two states does not affect our central results, as discussed in the Supplemental Material (SM) [50].At low temperature each atom oscillates around the minimum of its local potential, r 0 k , and its position can thus be written as r k = r 0 k + δr k , with δr k;µ = µ (b † k;µ + b k;µ ) being the atomic displacements from equilibrium.Here, = ( x , y , z ) is the vector of the characteristic lengths associated with the harmonic trapping potentials in the three spatial directions with µ = (2mω µ ) −1 .As a consequence, the two-body interaction depends on the displacements: Clearly, this implies that a coupling between electronic and vibrational degrees of freedom emerges.This situation is reminiscent of a mechanism for creating long-range spin models in arrays of trapped ions [51][52][53].In that case, the interplay of long-range Coulomb repulsion between the ions and laser induced spin-dependent forces results in an effective long-range spin-spin interaction and allows to simulate a rich variety of quantum systems.However, in contrast to the ions, Eq. ( 2) implies that in our setup the potential V(r k , r m ) couples electronic and vibrational degrees of freedom only when two atoms are excited, which is the origin of many-body spin interaction terms.
To demonstrate this, we focus on the strong confinement regime, in which the displacements δr k are much smaller than inter-atomic distances.Indeed, this represents the typical situation in Rydberg quantum simulators [6,9,38,39].By expanding the potential in Eq. ( 2) in a Taylor series to the first order in δr, the atom-atom interaction Hamiltonian acquires the form where Finally, since the spin-phonon coupling in Eq. ( 3) is linear in the bosonic operators, we can apply a polaron transformation, U (see SM [50]), to decouple spin and phonon dynamics.We obtain [50,51,54] with Here, we have introduced the coefficients V k;p,q = µ ( ω µ ) −1 W k,p;µ W k,q;µ and V k;m ≡ V k;m,m .Equations (6a) and (6b) show that, as consequence of the spin-phonon coupling, an effective atom-atom interaction emerges.The latter consists of an extra two-body [Eq.(6a)] and a novel threebody term [Eq.(6b)], whose strengths are both ∝ V k;p,q .Importantly, the coefficients V k;p,q depend on the trapping frequencies ω µ and are therefore tunable via the harmonic confinement.
The term H res in Eq. ( 5) describes a residual spin-phonon coupling, which is negligible in the limit |W k,m;µ | ω µ [50][51][52]54].In this regime the phonon dynamics decouples from the spins.The approximation further improves at temperatures low enough to ensure a 90% population in the vibrational ground state.Such temperatures can be experimentally achieved in state-of-the-art optical tweezers via Raman sideband cooling [34,55].Details on the validity of this spin-phonon decoupling approximation are provided in next section and in SM [50].
Microwave dressed Rydberg states.-Thestrength of the phonon-mediated effective interactions in Eqs.(6a) and (6b) is directly connected to the strength of the dipolar ones: This is because the coefficients W k,m;µ are proportional to the gradient of V(r k , r m ).Typical dipolar interactions exhibit a power-law behavior ∝ |r k − r m | −α (e.g., α = 6 for a van der Waals potential).In which case, one generally finds that V(r k , r m ) V k;m .This means that, in common situations, phonon-mediated interactions only represent a small correction.However, the interaction potential between excited atoms can be tailored via microwave (MW) dressing of two different Rydberg states [50,56,57], allowing to make the effective interactions dominant.In Fig. 2(a) we show one possible realization of such potential, obtained via MW dressing of the atomic levels |65S and |75P of 87 Rb atoms arranged on a square lattice.Here, a and a NNN = √ 2 are the distances at equilibrium between nearest neighbors (NNs) and next-nearest neighbors (NNNs), respectively.By properly choosing the MW field parameters (see SM for details [50]), the potential can be parameterized, to a good degree of approximation, as [50].MW dressing allows to control the values of the constants C 1,2 and c 1,2 in Eq. ( 7) independently and, in turn, to tune the strength of the dipolar potential (as well as its gradient) at NN and NNN distances, denoted by V 1 and V 2 , respectively.
Phonon-mediated interactions.-Forthe case shown in Fig. 2(a), we have V 1 / ≈ 0, V 2 / ≈ 2π × 0.3 MHz, and −1 dV/dR| R=a = 2π × 1.45 MHz.In this way we can thus achieve regimes dominated by the phonon-mediated interactions, whose strength along the µ direction is described by the parameter In this case, Eqs.(6a) and (6b) become where, explicitly, The symbols k, m and k, m denote the sum over NNs and NNNs, respectively, while k, p, q implies that the sum is restricted to sites satisfying Note that, due to the presence of the factors R0 k,m;µ , the terms ∝ V 3,µ strongly depend on the lattice geometry and, as we will show for the case of bipartite lattices, they give rise to anisotropic contributions in atom-atom interactions even if original dipolar forces are isotropic.
The strength of the phonon-mediated interactions can be tailored by tuning the trapping frequencies ω µ [see Eq. ( 8)], which are typically of the order of hundreds kHz [6,9,38,40].In particular, Eq. (9a) implies that it is possible to make the overall two-body term vanish and maximize the effects of three-body interactions.Recalling Eq. ( 5), in order to decouple the electronic and vibration degrees of freedom and to focus only on the spin dynamics we have to require |W k,m;µ | ω µ .On the other hand, to access regimes governed by the effective two-and three-body interactions, one should also consider 4) and (8), the above conditions translate into the following bounds on ω µ , In Fig. 2(b), we show typical values of the effective interaction strength V 3 for a square lattice geometry.The gray region denotes the regime where |W k,m;µ | ≤ ω µ , while along the red dashed curve As discussed in more details in SM [50], the leftmost condition in Eq. (10), which holds for any value of Ω, can be relaxed in the strong (effective) interaction regime, where V 3 /Ω |W k,m;µ |/( ω µ ), while the spin-phonon decoupling becomes exact in the classical limit (i.e., with vanishing Rabi frequency Ω).Thus, as can be seen in Fig. 2(b), the regime with V 2 ∼ V 3 can be accessed experimentally and corresponds to trapping frequencies and coupling strengths achievable in Rydberg atom tweezer arrays [6,40,42].Finally, we note that the time scales associated with the effective interaction dynamics, τ 3 = /V 3 , are < 50 µs in a wide region of the parameter space [i.e., the non-hatched area in Fig. 2  and are thus significantly shorter than the lifetime of the Rydberg states used in tailoring the MW-dressed potential of Fig. 2(a), which are of the order of hundreds of µs [50].

Phase diagram for a bipartite lattice.
-We now focus on a system of atoms arranged on a bipartite lattice and investigate the effects of the interplay between (two-body) dipolar and effective (two-and three-body) interactions on its phase diagram.The simplest case of a square lattice is shown in Fig. 3(a).The different contributions to atom-atom interaction are listed in Fig. 3(a,c).Importantly, the lattice-dependent structure of H 3B in Eq. (6b) implies that effective two-body interactions are attractive while, on the contrary, three-body terms have a repulsive character.This feature is quite general and, e.g., in Ising spin models on non-bipartite lattices (triangular, kagome) it could be employed to implement frustrated interactions [15,[19][20][21].The study of such phenomena will constitute the focus of future investigations.Due to the competition between two-and three-body interactions, we expect that different phases emerge.To map out the phase diagram we consider the classical limit (i.e., with vanishing Rabi frequency Ω) and determine its ground state through a classical Metropolis algorithm [58,59] by employing an annealing scheme [60].
Results are displayed in Fig. 4(a,b,c).Here, we show the behavior in the V 2 − V 3 plane (with V 1 > 0 and V 3, x = V 3,y ) of the average value of the Rydberg excitation density, n , of the density of dimers n dim , and of the density of trimers n trim [see Fig. 3(c,d)].Beyond the trivial states with all excited and all de-excited atoms, four further phases emerge, see Fig. stripe, dominated by NNN two-body (attractive) interaction ∝ V 2 , (3) frustrated striped phase with one missing line [here, the trimers occurring in (2) are melted due to the three-body repulsive contribution ∝ V 3 ], and (4) four-excitation clustered phase, dominated by attractive two-body interactions ∝ V 3 .Concerning this latter, we note that the transition is not as sharp as the other ones.Indeed, as can be seen from the last panel of Fig. 4(d), the lattice is not entirely covered by fourparticle clusters.This may suggest either that (4) is a liquid phase or that it represents a critical region.A full covering can be obtained for V 2 > 0, where attractive NNN interactions contribute to enhance the energy gain in forming clusters.
Interestingly, effective interactions due to spin-phonon coupling give rise to finite-size frustration phenomena even in a square lattice in the presence of isotropic dipolar interactions.This is manifest in the emergence of the different striped phases (2) and (3): see Fig. 4, which displays the case of a lattice with an even number of sites.On the contrary, if an odd number of sites is considered only a single regular striped phase emerges in this region of the phase diagram.However, a frustrated phase forms inside phase (1) (see SM [50]).
In non-square lattices, the geometrical factors characterizing phonon-mediated interactions [see Eq. ( 9)] give rise to anisotropic two-and three-body contributions even if the original dipolar interactions between atoms are isotropic.This can be seen in Fig. 3(d), where the various interaction contributions arising in a honeycomb lattice are displayed.Here, though the phase diagram is similar to the one shown in Fig. 4, non-trivial and anisotropic system configurations emerge [50].
The various phases shown in Fig. 4 can be probed in state-ofthe-art Rydberg simulators consisting of 2D defect-free arrays of optical tweezers [37].Indeed, as shown in SM [50], a significant part of the phase diagram in Fig. 4 can be mapped out by employing trapping frequencies ω µ ranging from a few tens to a few hundreds kHz, while the required dipolar interaction couplings are of the order of few MHz.The desired many-body states can be prepared by real-time control of Rabi frequency and detuning via a generalization of the rapid adiabatic passage protocol proposed in Refs.[22,23] and demonstrated in Ref. [24].The latter is perfectly compatible with the time scales associated with the effective interactions and, in turn, with the lifetime of the Rydberg states we considered.
Conclusions.-We have shown that electron-phonon interactions in Rydberg lattice quantum simulators permit the engineering of tunable multi-body interactions.We have illustrated the underlying mechanism in bipartite lattices, discussing in particular the case of an isotropic square lattice, where we studied the phase diagram in the classical limit.Going beyond this limit and considering the impact of quantum fluctuations (Ω > 0) will be possible in Rydberg quantum simulator experiments.Many future directions of this work can be envisioned: In particular, we expect that, as a consequence of the lattice-dependent structure of the induced interactions, peculiar two-and three-body terms would arise in non-bipartite lattices (e.g., triangular, kagome), allowing for the investigation of frustrated magnetism in spin models with non-trivial multi-body interactions.Furthermore, the mechanism leading from the spin-phonon coupling to effective many-body interactions can be generalized to different kinds of bare atom-atom potentials (e.g., exchange interactions, oscillating potentials) and may allow for engineering effective interactions with different structure and/or even n−body (with n > 3) contributions.
The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement No. 335266 ("ESCQUMA"), the EPSRC Grant No. EP/M014266/1, the EPSRC Grant No. EP/R04340X/1 via the QuantERA project "ERyQSenS", and the Deutsche Forschungsgemeinschaft (DFG) within the SPP 1929 Giant interactions in Rydberg Systems (GiRyd).The simulations used resources provided by the University of Nottingham High-Performance Computing Service.

Supplemental Material for "Engineering non-binary Rydberg interactions via phonons in an optical lattice"
In this Supplemental Material we provide additional details on the microwave dressing scheme employed to tailor the atom-atom potential in the main text.We then discuss in more details the approximations made in the main text, namely the small displacement expansion leading to Eq. ( 4) and the spin-phonon decoupling in system dynamics, and we investigate their validity.We comment about the values of ω µ corresponding to the phase diagram in Fig. 4 of the main text.We then discuss finite-size frustration phenomena in the phase diagram of a system arranged on a square lattice with an odd number of lattice sites.Finally, we inspect anisotropic effects in the case of a honeycomb lattice, showing its phase diagram and some typical configurations.

MICROWAVE DRESSING OF RYDBERG ATOMS
To engineer an atom-atom interaction potential leading to dominant phonon-mediated interactions, we consider two Rydberg atoms coupled by a microwave (MW) field [56,57].The two relevant levels on which we focus are denoted by |r and |p .The Hamiltonian is where Ω MW is the MW Rabi frequency and ∆ MW is the detuning.Here, C r , C p and C 3 are the respective dispersion coefficients of two atoms in |r and |p state, and of the dipolar interaction.n (j) σ = |σ (j) σ (j) | is the density operator associated to the state |σ (j) .
By considering different Rydberg states and MW parameters it is possible to tailor the strength of the nearest neighbor (NN) and next-nearest neighbor (NNN) interactions, denoted by V 1 and V 2 respectively, and to control the slopes of the potential (forces) at these points.As stated in the main text [see Fig. 2(a) and Eq. ( 10)], we are interested in a case with small values of V 1 > 0 and V 2 < 0 but with a large magnitude of the gradient of the interaction potential at NN distances, denoted by a.As an example, we consider two 87 Rb atoms and focus on the MW dressing of Rydberg levels |r = |65S and |p = |75P [see Fig. S1(a)].Their lifetime is 325.03 µs in |r and 917.87 µs in |p state [61].The corresponding interaction strengths are C r / = 2π × 0.37 × 10 3 GHz µm 6 , C p / = −2π × 1.84 × 10 3 GHz µm 6 , and C 3 / = 2π × 1.52 MHz µm 3 .For a MW Rabi frequency Ω MW / = 2π × 230 MHz and detuning ∆ MW / = −2π × 170 MHz, the higher energy potential V h (R) at the top of panel (b), which is also shown in panel (c) and whose derivative is reported in panel (d), meets the above criteria if we consider a lattice spacing a ≈ 5.30µm.In particular, note that the slope at NNN distance, √ 2a, is an order of magnitude smaller than the slope at a. Therefore, contributions to the phonon-mediated interactions (which are proportional to the gradient of the interaction potential; see main text) arising from distances R > a can be neglected.

DETAILS ON THE APPROXIMATIONS AND VALIDITY CHECKS
We provide here additional details about the derivation of Eq. ( 3) and the polaron transformation employed in the main text to decouple the spins and phonon dynamics in Eq. ( 5).In particular, in the strong confinement regime, the Rydberg interaction potential V(r k , r m ) can be expanded in a Taylor series to the first order in δr as By substituting Eq. into the interaction Hamiltonian H int [Eq.( 2) of the main text] we obtain Eq. ( 3) of the main text, , and the prime in the sum meaning k m.
Since in Eq. (S3) the coupling between spin and phonon degrees of freedom is linear in the phonon operator b k;µ , their dynamics can be decoupled by applying to the full system Hamiltonian H = H sp + H int the canonical (polaron) transformation defined by [51,54] We obtain with H spin the Hamiltonian governing the (decoupled) spin dynamics, H ph = k,µ ω µ b † k;µ b k;µ the phonon Hamiltonian, and H 2B and H 3B defined in Eq. ( 6) of the main text.The last term is a residual spin-phonon coupling and is given by Effects of this residual coupling are therefore negligible if |W k,m;µ | ω µ .Furthermore, two additional interesting regimes can be identified.In the classical limit Ω → 0 (considered in the last section of the main text), the canonical transformation results in a complete spin-phonon decoupling, On the other hand, we note that, in the strong effective two-and three-body interactions regime, the error introduced by H res can be neglected if V 3,µ /Ω |W k,m;µ |/ ω µ .The latter condition leads to We briefly comment here on the situation in which the trapping frequencies corresponding to the ground state, ω µ , and to the Rydberg state, ω µ = ω µ +δ µ are different.Here, δ µ is the trapping frequency mismatch.In this case, the phonon part of the Hamiltonian can be written as [64] and, thus, an extra spin-phonon coupling contribution appears.Nevertheless, a renormalized version of the canonical transformation in Eq. (S4), with results in a transformed Hamiltonian in which the spin terms contained in H spin have exactly the same structure as in Eq. (S5), though with renormalized coefficients.On the other hand, the residual spin-phonon coupling contribution H res of Eq. (S6) now contains the extra term k,µ ω µ δ µ n k b † k;µ b k;µ .The latter limits the validity of our approximation up to times t (ω µ δ µ ) −1 , which for standard experimental parameters are larger than the typical interaction induced time scales [see Fig. S4(b)].Moreover, ponderomotive bottle beam traps have recently been proposed for the realization of Rydberg tweezer arrays satisfying the ground-Rydberg "magic" condition [49], allowing to minimize the trapping frequency mismatch δ µ .
To check the validity of the approximations we made in the main text (i.e, the small displacement expansion in Eq. ( 3) and the decoupling between phonon and spin dynamics in mapping out the phase diagram), we consider a minimal system of three atoms aligned along the x−axis with an intermediate value of Ω, corresponding the worse case scenario.See Fig. S2(a).For the sake of simplicity, we only consider one longitudinal phonon mode per atom.We focus on the spin dynamics starting from a product state in which the phonons are initialized in a thermal state at inverse temperature β, ρ ph = e −βH ph /Tr e −βH ph , and with the spins prepared in the pure state ρ spin = |ψ 0 ψ 0 |, with |ψ 0 = s z 1 , s z 2 , s z 3 and s z i the projection of the i−th spin on the z−axis.We then evaluate numerically [62,63] the fidelity F (t), defined as the overlap between the spin state at time t, ψ approx (t) , obtained in the small displacement and spin-phonon decoupling limit (i.e., with time-evolution governed by H spin = H L + H 2B + H 3B ) and the full system density matrix, evolved with the full Hamiltonian H of Eq. ( 2) in the main text: As can be seen from Fig. S2(b,c), F (t) remains close to 1 over a significant time window, especially when small temperatures are considered [panel (b)].Here, the latter have been chosen in order to correspond to an average phonon occupation number n ph = 0.5 [panel (b)] and n ph = 0.05 [panel (c)], which can be achieved in optical tweezer setups via Raman sideband cooling [34,55].In Fig. S2(d-i) and Fig. S3 we compare the time evolution of the expectation values of the population of the Rydberg state of atom 1, n 1 , and of various projectors on states with σ z 1 , σ z 2 , σ z 3 spin components along the z−direction, e.g., P ↑↑↑ = n 1 n 2 n 3 , for the dynamics governed by H and H spin for two different initial spin configurations.In both cases, the agreement between the approximate and full dynamics is excellent.

TRAPPING FREQUENCIES RANGE, EFFECTIVE INTERACTIONS TIME SCALES AND CLASSICAL GROUND STATE PREPARATION
In Fig. S4(a) we reproduce the V 2 − V 3 plane of Fig. 4 in the main text and show the values of the trapping frequencies ω µ corresponding to the values of V 3 used in drawing the system phase diagram.As can be seen, in the majority of the V 2 − V 3 plane, the values of ω µ are both within the approximation validity regimes shown in Fig. 2   values.Figure S4(b) displays the time scale associated with the effective two-and three-body interactions, i.e. 1/τ 3 = V 3 / , which are considerably lower than the lifetime of the Rydberg levels considered in engineering the MW dressed potential in Section .
To illustrate how the classical ground states of the system shown in Fig. 4 of the main text can be prepared, we focus on a 3 × 3 spin square lattice described by H spin .Its phase diagram, obtained through exact diagonalization and corresponding to the same parameters used in Fig. 4, is shown in Fig. S4(c).Note that, as a consequence of the reduced system size, the only non-trivial phase emerging here is the striped one, [phase (2) in Fig. 4].As an example, the ground state corresponding to the parameters identified by the red cross in Fig. S4(c) can be prepared by following the rapid adiabatic passage procedure described in Ref. [24]: By setting a large value of the detuning ∆ in and a vanishing value of the Rabi frequency Ω, the system is initialized at first in a state with all the atoms in their ground state (i.e., with n = 0).Then, Ω is slowly increased up to Ω in , inducing a coupling between states belonging to the various manifolds with different Rydberg excitation density.Subsequently, the detuning is decreased from ∆ in to ∆ fin and, finally, Ω is reduced to its final value Ω fin = 0.As can be seen from the lowest, thick curve in panels (d) and (e), this procedure leads to a degenerate ground state with the desired value of n and well-seperated from high-lying many-body eigenstates.Here, the ground state degeneracy is associated with all the possible arrangements of the excited atoms.As shown in Ref. [24], an overall laser sweep duration of a few µs is sufficient to ensure a good degree of adiabaticity and, therefore, this protocol is perfectly compatible with the typical time scales of our system [see, e.g., Fig. S4(b)].

PHASE DIAGRAM FOR A SQUARE LATTICE WITH ODD NUMBER OF SITES
In the main text we have investigated the phase diagram of a system of atoms arranged on a square lattice with an even number of sites (see Fig. 4 in the main text).In that case, different phases are present and, in particular, finite-size frustration phenomena emerge [see, e.g., phases (2) and (3) in Fig. 4 in the main text].Here, we inspect what happens if a square lattice with an odd number of sites is considered.The phase diagram and some typical system configurations are shown in Fig. S5.The two (frustrated) phases ( 2) and (3) of the even number of sites case are now merged [see phase (3) in Fig. S5(a)] in a single phase showing regular strips of dimers.However, finite-size frustration phenomena emerge in the small V 3 region of the phase diagram: see phases (1) and (2) in Fig. S5, which consists of stripes of single excitations with a (a) dimer or (b) a missing strip, respectively.We note that the four-excitation clustered phase occurs in this case as well and exhibits a behavior similar to the one discussed in the main text.

PHASE DIAGRAM FOR A HONEYCOMB LATTICE
In this section we finally investigate the phase diagram for a honeycomb lattice.In this case, as stated in the main text, the geometrical factors present in Eq. ( 11) can be exploited to realize non-trivial and anisotropic interactions [see Fig.  (a,b,c) show the phase diagram in the V 2 − V 3,y plane (with V 1 > 0) of a system with L 2 = 64 unit cells (and 2L 2 = 128 sites) with fixed V 3, x / = 2π × 1 MHz.Even if the phase diagram looks similar to the one obtained for a square lattice and discussed in Fig. 4 of the main text, anisotropic phenomena can be seen in typical system configurations [see Fig. S6(d)].In particular, as can be seen from Fig. S6(f), with the previous choice of parameters (e.g., with V 3,x / = 2π × 1 MHz) the three-body interaction along the x−direction changes from repulsive to attractive when V 3,y = 3V 3, x .Above this threshold we thus expect the formation of "chains" of excitations along x to be energetically preferred.Indeed, the emergence of such structures is clearly visible in the typical configurations associated to phases (2) and (3) in Fig. S6(d).

and σ µ k
are the Rydberg number operator and Pauli matrices acting on the atom at site k and position r k = (x k , y k , z k ), respectively.Any two atoms at lattice positions r k and r m , if excited to the Rydberg state, interact through the two-body potential V(r k , r m ), which depends on the inter-particle distance |r k − r m | [1, 3, 7].The overall Hamiltonian is therefore H = H sp + H int with

FIG. 3 .
FIG. 3. Interaction terms in bipartite lattices.(a) Square and (b) honeycomb lattice with NNs (orange dots) and NNNs (blue squares)interacting through dipolar interactions (red, solid and blue, dashed lines, respectively).As a consequence of the phonon-mediated effective multi-body interaction terms arise.Two-body (green, left) and three-body (purple, right) contributions along the horizontal (solid) and vertical (dotted) direction are shown, with their corresponding sign, in panels (c) and (d) for the square and honeycomb geometries, respectively.Note that, in the latter case, horizontal (solid) terms contribute to both x and y directions, resulting in anisotropic interactions.
FIG. 4. Phase diagram (square lattice) as a function of V 2 / and V 3 / .In panel (a) we show the average density of Rydberg excitations n as a function of V 2 / and V 3 / (units 2π × MHz), with V 3, x = V 3,y , for a square lattice.Note that V 3,µ ∼ ω −2 µ and, therefore, it can be controlled by the confinement strength.In panels (b) and (c) are reported the density of dimers n dim and trimers n trim , respectively.Dashed red lines represent guides for the eye to distinguish between the different system phases.In panel (d) we show typical configurations in the different regions of the phase diagram.Dark (red) spots correspond to excited atoms.See text for details.In all panels, L = 10, V 1 / = 2π × 0.2 MHz and ∆/ = 2π × 1 MHz.
FIG. S1.MW dressed potential.(a) Each atom is modeled as a three level system, with ground state |g and Rydberg states |r and |p .The ground state |g and the Rydberg state |r are laser coupled Rabi frequency Ω and detuning ∆.A MW field couples the two Rydberg states |r and |p with Rabi frequency Ω MW and detuning ∆ MW .(b) Interaction potentials (units 2π ×MHz) as a function of the inter-atom distance R (units µm) for a pair of 87 Rb Rydberg atoms with |r = |65S , |p = |75P , Ω MW / = 2π × 230 MHz, and ∆ MW / = −2π × 170 MHz.Dashed lines denote the interaction potential corresponding to Ω MW = 0. Here, we have Rydberg states |r (black) and |p (blue), and dipolar interaction (red).After turning on the MW coupling, the potentials are modified by the mixing of different states (solid lines).We will focus on the repulsive potential at the top of the panel, V h (R) (blue curve).(c) Plot of the highest potential V h (R) in (a) in the presence of MW coupling.With respect to (b), the potential is shifted by V/ = 2π × 320.4 MHz in order to make it asymptotically vanishing.For a two-dimensional square lattice with a = 5.3 µm, values at NNs and NNNs areV 1 = V h (a) = 2π × 0.01 MHz and V 2 = V h ( √ 2a) = −2π × 0.15MHz, respectively.(d) The negative slope (force) of the potential shown in (c).The slope is −2π × 1.45 MHz/µm at R = a and 2π × 0.10 MHz /µm at R = √ 2a.
(a)  and within experimentally achievable

5 Pσ z 2 σ z 3 (
FIG. S2.Full vs approximated dynamics.(a) One dimensional three-atom setup used to verify the approximations made in the main text.For simplicity, we assume one longitudinal phonon mode per atom.(b,c) Fidelity F (t), quantifying the overlap between ψ approx (t) and the exact density matrix (ρ spin ⊗ ρ ph )(t) for the system shown in panel (a).Different curves correspond to different initial atomic configurations (see legend) while phonon modes are initialized in a thermal state with average phonon occupation (b) n ph = 0.5 and (c) n ph = 0.05.Parameters are:C 1 /( a 6 ) −1 = 2π × 2.5 MHz, C 2 /( a 6 NNN ) −1 = −0.1C 1 /( a 6 ) −1 , ω = 2π × 0.3 MHz, Ω/ = 2π × 1 MHz, ∆/ = 2π × 1MHz, and a = 5µm.(d) Time evolution of the Rydberg population of atom 1 and (e-i) of the projectors P σ z 1 σ z 2 σ z 3(t), with the phonon modes initialized in a thermal state with average phonon occupation number n ph = 0.05 and the spins in the ↓↓↓ configuration.Blue, solid lines correspond to the approximate dynamics governed by H spin , while the yellow dashed curves to the one generated by the full Hamiltonian H. Same parameters as in panels (b,c).

25 P
FIG. S3.Full vs approximated dynamics -spin initialized in an excited state.(a-f) Same as FigS2(d-i) with spins initialized in the ↓↑↑ configuration.
FIG. S4.Frequency range, effective interaction time scales, and ground state preparation.(a) Density plot in the V 2 − V 3 plane of the trapping frequencies (units 2π × MHz) corresponding to the values of the effective interaction V 3 / on the y−axis.(b) Density plot of the time scale τ 3 = /V 3 associated with the effective interaction V 3 / in the C 1 /( a 6 ) − ω plane (units 2π × MHz).As in Fig. 2(b) in main text, the gray area denotes the regime in which |W k,m;µ | ≤ ω µ (with equality corresponding to the red solid line).Along the red dashed curve V 3 = V 2 .Here, C 2 /( a 6 NNN ) = −0.1C 1 /( a 6 ).(c) Exact phase diagram for a 3 × 3 square lattice spin system for V 1 / = 2π × 0.2 MHz, V 2 / = −2π × 0.7 MHz, V 3 / = 2π × 1 MHz, Ω/ = 2π × 0.01 MHz, and ∆/ = 2π × 1 MHz [corresponding to the red cross in panel (c)].(d, e) Many-body energy spectrum of H spin for a 3 × 3 square lattice as a function of (d) ∆/ and (e) Ω/ (units 2π × MHz).Here, ∆ in / = 2π × 10 MHz, ∆ fin / = 2π × 1 MHz, Ω in / = 2π × 1 MHz, and Ω fin = 0.The colormap indicates the average excitation density n in the various eigenstates.The lowest, thick line denotes the ground state.For a large value of the detuning ∆ the spectrum separates in different manifolds with fixed number of excited atoms (see labels).In the limit of vanishing Ω the ground state becomes 6−fold degenerate, where each one of the 6 eigenstates corresponds to a specific arrangement of the N ex = 5 excited atoms on the 3 × 3 lattice.

FIG. S5 .
FIG.S5.Phase diagram for a square lattice with an odd number of lattice sites as a function of V 2 and V 3 .In panel (a) we show the average density of Rydberg excitations n as a function of V 2 / and V 3, x / = V 3,y / (units 2π × MHz) for a square lattice.In panels (b) and (c) are reported the density of dimers n dim and trimers n trim [see Fig.3(b)in main text for details], respectively.Dashed red lines represent guides for the eye to distinguish between the different system phases.In panel (d) we show typical configurations in the different regions of the phase diagram.Dark (red) spots correspond to excited atoms.In all panels, L = 9, V 1 / = 2π × 0.2 MHz and ∆/ = 2π × 1 MHz.