Quantum Spin Ice and dimer models with Rydberg atoms

Quantum spin ice represents a paradigmatic example on how the physics of frustrated magnets is related to gauge theories. In the present work we address the problem of approximately realizing quantum spin ice in two dimensions with cold atoms in optical lattices. The relevant interactions are obtained by weakly admixing van der Waals interactions between laser admixed Rydberg states to the atomic ground state atoms, exploiting the strong angular dependence of interactions between Rydberg p-states together with the possibility of designing step-like potentials. This allows us to implement Abelian gauge theories in a series of geometries, which could be demonstrated within state of the art atomic Rydberg experiments. We numerically analyze the family of resulting microscopic Hamiltonians and find that they exhibit both classical and quantum order by disorder, the latter yielding a quantum plaquette valence bond solid. We also present strategies to implement Abelian gauge theories using both s- and p-Rydberg states in exotic geometries, e.g. on a 4-8 lattice.


I. INTRODUCTION
The ice model has been fundamental in furthering our understanding of collective phenomena in condensed matter and statistical physics: in 1935 Pauling provided an explanation of the 'zero-point entropy' of water ice [1] as measured by Giauque and Stout [2], while Lieb demonstrated with his exact solution of the ice model in two dimension [3] that there exist phase transitions with critical exponents different from those of Onsager's solution of the Ising model. The experimental discovery [4] of a classical spin version of the ice model [5] has in turn generated much interest in the magnetism community [6][7][8][9].
More recently, quantum ice models [10][11][12][13][14] have attracted a great deal of attention in the context of phases exhibiting exotic types of orders, such as resonating valence bond liquids [15,16] or quantum Coulomb phases [12,[17][18][19]. They form part of a broader family of models, which also includes quantum dimer models or other quantum vertex models [20], in which locally a hard constraint is imposed, such as the ice rules defined below. Such a constraint can then endow the configuration space with additional structure -most prominently, an emergent gauge field which can be the basis of the appropriate effective description at low energies [20][21][22]. This is an important phenomenon as it is perhaps the simplest way of obtaining gauge fields as effective degrees of freedom in condensed matter physics. More broadly, this is part of a longrunning search for magnetic materials hosting quantum spin liquids [7].
In the present work we address the problem of physically realizing models of quantum spin ice and quantum dimer models in two spatial dimensions with ultracold atoms in optical lattices. Our proposal builds on the recent experimen- * alexander.glaetzle@uibk.ac.at tal advances, and opportunities in engineering many-body interactions with laser excited Rydberg states [23][24][25][26][27][28][29][30][31][32][33][34][35][36][37]. In particular, this will allow us to develop a Rydberg toolbox for the complex interactions required in 2D quantum ice models [8]. Our investigation also fits into the broader quest for the realization of synthetic gauge fields with cold atoms. While much effort is being devoted to the generation of static gauge fields [38], e.g. on optical flux lattices, here we follow the strategy of generating a dynamical gauge field [12, 17-19, 21, 22, 39-42] emerging upon imposition of the ice rule.
While in condensed matter systems the interactions underlying ice and spin ice arise naturally in a 3D context [see Fig. 1(a)], the 2D quantum ice on a square lattice requires a certain degree of fine tuning of the relevant interactions [see Fig. 1(b)]. In 3D spin ice materials, for example, the ions of magnetic rare-earth atoms reside on a pyrochlore lattice, representing a network of corner-sharing tetrahedra. Magnetic interactions in combination with crystal fields give rise to a low energy manifold of states on each tetrahedron consisting of six configurations, in which two spins point inward and two spins point outward [compare Fig. 1(c)]. In a similar way, in water ice each O 2− atom in a tetrahedrally coordinated framework has two protons attached to it, giving rise to a manifold of energetically degenerate configurations. 2D models of ice and spin ice can be understood as projection of the pyrochlore on a square lattice (see Fig. 1), where again the low energy configurations of spins residing on the links obey the "ice rule" two-in and two-out at each vertex. While these 2D ice models play a fundamental role in our theoretical understanding of frustrated materials, a physical realization requires a precise adjustment of the underlying interactions -different local configurations which are symmetry distinct need to be at least approximately degenerate; the required fine-tuning however needs to be delicately directionally dependent, as in the pure ice model, the ratio of some interactions between different pairs of equidistant spins vanishes, see Fig. 1(b). Things are not all hopeless, however, as there exist a number of settings in which partial In spin ice materials the magnetic moments (yellow arrows) of rare-earth ions are located on the corners of a pyrochlore lattice, which is a network of corner-sharing tetrahedra. They behave as almost perfect Ising spins and point along the line from the corner to the centre of the tetrahedron, either inward or outward. Due to the different Ising-axes of the spins this results in an effectively antiferromagnetic interaction which is frustrated. (b) Projecting the 3D pyrochlore lattice onto a 2D square lattice yields a checkerboard lattice where tetrahedrons are mapped onto crossed-plaquettes (lightblue). Interactions between two spins located on or lattice sites have to be (i) step-like as a function of the distance, (ii) anisotropic and (iii) require a bipartite labelling of the lattice sites. (c) Degenerate ground state configurations of spins on a crossed-plaquette. They obey the ice-rules, which enforce two spins pointing inward and two spins pointing outward at each vertex.
progress has been made to realizing such models. In two dimensions, artificial structures using nanomagnetic [43] or colloidal arrays [44] have been proposed, including strategies for tuning the interactions appropriately [45]. The present proposal with Rydberg atoms is unique, however, as it combines both the possibilities of engineering the complex interactions using Rydberg interactions with the accessibility of the quantum regime in cold atom experiments. Alkali atoms prepared in their electronic ground state can be excited by laser light to Rydberg states, i.e. states of high principal quantum number n [46][47][48][49]. These Rydberg atoms interact strongly via the van der Waals interaction exhibiting the remarkable scaling V VdW ∼ n 11 , and which exceed typical ground state interactions of cold atoms by several orders of magnitude. In an atomic ensemble the large level shifts associated with these interactions implies that only a single atom can be excited to the Rydberg state, while multiple excitations are suppressed within a blockade radius determined by the van der Waals interactions and laser parameters [50,51]. This blockade mechanism results in novel collective and strongly correlated many-particle phenomena such as the formation of superatoms, and Rydberg quantum crystals [46][47][48][49]. In present experiments the emphasis is on isotropic van der Waals interactions [23][24][25][26][27][28][29][30][31][32][33][34][35][36], which, for example, can be obtained by exciting Rydberg s-states using a two-photon excitation scheme. In contrast, we will be interested below in excitations of Rydberg p-states, where the van der Waals interactions can be highly anisotropic, and we will discuss below in detail the controllability of this anisotropy, shape and range of these interactions via atomic and laser parameters for the case of Rb atoms. In the present context this will provide us with the tools to engineer the required complex interaction patterns for 2D quantum spin ice and dimer models. The setup we will discuss will consist of cold atoms in optical lattices, where the strong Rydberg interactions are weakly admixed to the atomic ground states [52][53][54][55][56], thus effectively dressing ground state atoms to obtain the complex interactions in atomic Hubbard models required for the realization of 2D quantum spin ice and dimer models.
We then proceed to numerically analyze the family of Hamiltonian realizable with this toolbox. We verify that it contains two phases exhibiting distinct types of order by disorder [57,58]. Most remarkably, quantum order by disorder -due to the presence of quantum dynamics in the ice model -realizes the plaquette valence bond solid as an unusual non-Néel phase of a frustrated magnet. This terminates when classical degeneracy lifting takes over.
The latter phase is conventional in that it is diagnosed by a conventional spin order parameter, which would manifest itself in a Bragg peak in the structure factor. By contrast, the valence bond solid would be diagnosed by higher order 'string' correlators, and is thus fundamentally different from some other instances of quantum order by disorder, where finally the mechanism, but not quite so much the order parameter, is exotic [58,59]. Notably, precisely such order parameters have become accessible to experimental measurements recently [60], so that our proposal not only covers the setting for realizing quantum order by disorder, but also the means for detecting it.
In addition, at finite temperature, we find that the classically ordered state, even in the absence of quantum dynamics, melts into a classical version of a Coulomb phase, namely a Coulomb gas in which thermally activated plaquettes violating the ice rule play the role of positive and negative charges. As these interact via an entropic two-dimensional (logarithmic) Coulomb law, this phase is only marginally confined [61,62].
While throughout the paper we will be mostly interested in quantum ice models, which require the development of advanced interaction-pattern design, we will also discuss a second strategy to implement constrained dynamics, and in particular quantum dimer models, with Rydberg atoms. It relies on combining simple interaction patterns, such as the ones generated by s-states, with complex lattice structures, which can be realized either via proper laser combination or by the recently developed optical lattice design with digital micromirror devices [63]. These models extend the class of dynamical gauge fields in AMO systems to non hyper-cubic ge-ometries. Overall, the ability to synthetically design Abelian dynamical gauge fields with discrete variables also establishes interesting connections with high energy physics, where these theories are usually refereed as quantum link models [64,65]. Within this context, the key developments in engineering pure gauge theories can be combined with other schemes, where dynamical matter is included, which have already been proposed in the context of cold atom gases [65].
The the paper is structured as follows. In Sec. II, we briefly provide the background on quantum ice models needed for digesting the remainder of the material. Sec. III outlines the implementation of a family of model Hamiltonians approximating the quantum ice model using atoms in optical lattices weakly admixed with a Rydberg p-state. Our model Hamiltonians are then analyzed for their phase diagram in Sec. IV. Sec. V presents strategies to implement simpler Abelian gauge theories using both s-and p-Rydberg states in exotic geometries, e.g. a 4-8 lattice. The paper closes with a summary and outlook in Sec. VI.

II. THE QUANTUM ICE MODEL
This section provides a brief overview over the statistical mechanics of the ice model, the emergence of a gauge field, and the challenges in realizing such a model experimentally.
A. The configuration space: ice rules and emergent gauge fields The ice model on the square lattice, also known the sixvertex model or simply square ice [3], has Ising degrees of freedom residing on the links of a square lattice. They can either be thought of as 'fluxes', Ŝ z i , pointing in either of the directions along the bond i, see layer (i) of Fig. 2(a); or equivalently can be mapped onto spins S z i which point up or down depending on the direction of the flux [see layer (ii) of Fig. 2(a)].
Only configurations satisfying the ice rule are permitted, which stipulate that the spins on each vertex add up to zero -there are 4 2 = 6 ways of arranging this, see panel (i) of Fig. 2(b). The number of configurations satisfying the ice rule grows exponentially with the size of the system -for a lattice of N spins, there are (4/3) 3N 4 ice states [3]. The origin of the emergent gauge field is transparent in flux language, where it implies that the lattice divergence of the flux field vanishes: defining the x(y) component of a twodimensional vector flux b to be the flux along the corresponding links emanating from a vertex in the positive x(y) direction, one has (1) Note that a gauge field a has appeared naturally as a consequence of enforcing the ice rule, just as it does in magnetostatics, where Maxwell's law for the magnetic field ∇ · B = 0 leads to the introduction of the familiar vector potential A.
The spins of the 2D ice model on a checkerboard lattice can be interpreted as (i) "fluxes", {Ŝ z i }, pointing either inward or outward of a specific crossed plaquette. This requires a bipartite labelling of the plaquettes (light blue and light magenta plaquettes) since an outward pointing flux vector corresponds to an inward pointing flux vector for the neighboring plaquette. (ii) They can be interpreted as spins, {S z i } aligned perpendicular to the plane, pointing either up (red arrows) or down (black arrows). In the right inset we identify a spin pointing up, S z i = + 1 2 , with a flux vector pointing from the magenta to the blue plaquette,Ŝ z i = + 1 2 and vice versa. (iii) Spins, {S z i }, can be mapped onto hard-core bosons, n i ∈ {0, 1}. Here, e.g. a particle (red circle), n i = 1, corresponds to a spin pointing upward, S z i = + 1 2 , while an empty lattice site (white circle), n i = 0, corresponds to a spin pointing downward, S z i = − 1 2 . (b) The six icerule states correspond to vertex configurations with two hard-core bosons and two empty lattice sites.
In the present example in two dimensions, where the flux is a two component vector b, and the scalar constraint ∇ · b = 0 fixes one degree of freedom, a only has one physical degree of freedom left -it can be thought of as a scalar, usually referred to as a height: a = hz 'in the z-direction' [61,62]. Defects in this height field -forbidden in the six vertex model but allowed when violating the ice rule comes only with a finite energy penalty -are then known as charges or monopoles, which carry a gauge charge with respect to the emergent gauge field.
Having enforced the ice rule, the natural degree of freedom is thence an emergent gauge field a -it is in this way that gauge fields quite generically emerge in condensed matter physics, with a constraint arising either from the need to satisfy a dominant term in the Hamiltonian, or a microscopic relation on the local Hilbert space [20]. B. Realization, and fine-tuning in d = 2 The ice rule on a given vertex involves four spins, but it can be enforced via a pairwise interaction: if all four spins on a vertex interact antiferromagnetically and equally -described by the Hamiltonian where + denotes a crossed plaquette in 2D or a tetrahedron in 3D -the resulting ground states are those which obey the ice rules, see panel (i) of Fig. 2(b). In three dimensions, equality of the pairwise interactions can be symmetry-generated -by placing the spins on the corners of a tetrahedron, any antiferromagnetic interaction depending only on the distance between the spins will yield the ice rule. By contrast, in two dimensions, a tetrahedron becomes a square with interactions also across the diagonal (Fig. 1), which are no longer symmetry equivalent to those along the edges [66].
In particular, interactions, V i j (r), between two spins i and j located on the bonds of a checkerboard lattice separated by a distance r, have to fulfill three demeaning properties [see Interactions have to be strongly anisotropic.
This is illustrated in panel (ii) of Fig. 1(b). Particles which belong to the same vertex interact strongly (red arrow), while particles which do not belong to the same vertex do not interact (gray arrow). Thus, for particles in panel (ii) one needs an interaction which satisfied V (ϑ = 0) = 0 (gray arrow) and V (ϑ = π/2) =Ṽ 0 (red arrow), where the angle ϑ is defined in the inset.

2.
Step-like potentials: All four particles which belong to the same vertex (enclosed by light blue squares) interact with the same strengthṼ 0 , independent of their distance, either a or a √ 2, where a is the lattice spacing. Obviously, an interaction of the form 1/|r| α would not suffice. It is therefore necessary to have step-like potentials which fulfill V i j (|r| < r c ) =Ṽ 0 0 and 3. Bipartite-lattice structure: Furthermore, panel (iii) of Fig. 1(b), shows that the desired interaction properties cannot be satisfied by a homogeneous interaction pattern, but require a bipartite structure [squares and circles in Fig. 1(b)] where the angular dependence on the interaction depends on the lattice bipartition. For example, in the last paragraph we enforced that V (ϑ = 0) = 0 [see panel (ii)] but the opposite is true for particles, see panel (iii). Here, V (ϑ = π/2) = 0 but V (ϑ = 0) =Ṽ 0 . On top of that, mixed interactions between and particles on the 45 degree lines should obey V (ϑ = ±π/4) =Ṽ 0 in order to ensure that all six possible interactions at a specific vertex are the same, see panel (i).
It is these three countervailing requirements that we manage to satisfy approximately by using Rydberg dressed atoms to engineer an appropriate quantum Hamiltonian (Sec. III).
C. Adding quantum dynamics, and quantum order by disorder While the properties of the two-dimensional ice model were broadly understood a long time ago, the question what a quantum version of it would look like was not posed until much later [10]. Unlike in, say, a transverse field Ising model, where the simplest quantum dynamics consists of reversing a single spin, the ice model does not permit such single-site configuration changes, as these would lead to a violation of the ice rule.
The smallest cluster which may flip consists of a closed flux loop around a plaquette, denoted by (see Fig. 3), and this will be the second ingredient that our work will implement (Sec. III). What this amounts to in the language of gauge theory is the addition of a field conjugate to the height/gauge field a = hz -or in more familiar parlance of electromagnetism, the appearance of an (emergent) magnetic field alongside an (emergent) electric one [22].
In two dimensions, adapting a celebrated result by Polyakov (which does not apply straightforwardly as it is based on Lorentz invariance which does not hold a priori for our emergent field), it is known that the (emergent) electromagnetism is confining. As a consequence, the emergent excitations cannot spread freely over the system, being bounded by an effective string tension due to the gauge fields. Concretely, one finds a phenomenon known as order by disorder [58]. The quantum dynamics mixes the degenerate ice states into a superposition to form the quantum ground state. Even though the quantum dynamics induces fluctuations ('disorder'), the resulting ground state exhibits long-range order. This order takes the form known as a plaquette valence bond solid (Fig. 3) [67], which breaks translational symmetry. Such valence bond solids occur frequently in the theory of quantum magnets, but they are not commonly realized in experiment.
In Sec. IV, we show that the model Hamiltonian we provide a recipe for does exhibit this kind of order-by-disorder plaquette phase, and we discuss how to detect this kind of exotic order. In addition, we find that for weak quantum dynamics, a different, classical type of symmetry-breaking occurs (see in Fig. 3(c)). This happens because different ice states are only approximately degenerate for our engineered Hamiltonian, and the residual energy differences are sufficient to select a particular ordered configuration. D. Relation between quantum ice, Bose-Hubbard models and dimer models As a starting point for our implementation, we will consider a hard-core extended Bose-Hubbard Hamiltonian on a 2D checkerboard lattice: Here, b † i (b i ) is an operator that creates (annihilates) a hardcore boson on site i which obey an on-site contained b 2 i = b †2 i = 0. The rate J h is the nearest neighbor (NN) hopping amplitude and V describes a repulsion between all atoms sitting close to the same vertex. The operator n i = b † i b i counts the number of bosons at site i and can be either zero or one, n i ∈ {0, 1}. The summation runs over nearest neighbors only. The hard-core boson model can be mapped to a spin-1/2 model using the transformation [68] Expanding the last term gives the two-body interaction proportional to J i j S z i S z j and an additional magnetic field term proportional to J i j S z i which is constant after fixing an initial number of particles. This will fix the gauge sector in the gauge theory description [20].
In order to implement the constrained model of Eq. (2) we demand (i) anisotropic and (ii) step-like interactions between (iii) two species of particles, as discussed in Sec II B, that is V i j has to fulfill i jṼ i j n i n j =Ṽ 0 withṼ 0 a constant interaction between all particles belonging to the same vertex denoted by +. Under this assumptions, in the limitṼ 0 J h the Bose-Hubbard Hamiltonian of Eq. (4) maps onto the spin ice Hamiltonian of Eq. (2). The specific form of V i j ensures that all six interactions between particles which belong to the same vertex are equal and interactions between particles which do not belong to the same vertex vanish. In the case of total half-filling of the initial bosons, N = L/2, one has i S z i = 0: this fixes the effective dynamics on the aforementioned ice manifold of interest. In the case different initial fillings are considered, one has access to different quantum dynamics: a notable case is the N = L/4 case, which defines a constrained dynamics on a manifold where a single boson sits close to each vertex [10]. The effective description is then the same of hard-code quantum dimer models on a square lattice.

III. QUANTUM ICE WITH RYDBERG-DRESSED ATOMS: EXPLOITING P-STATES
We now turn to the realization of the extended 2D Bose Hubbard Hamiltonian of Eq. (4) with cold atoms in optical lattices. The key challenge is the implementation of the in-teractionsṼ i j with constraints represented in Eq. (6). We will show below that this can be achieved via the very anisotropic Rydberg interactions involving laser excited p-states of Rubidium atoms.
A. Single-particle Hamiltonian on a bi-partite lattice In our setup we consider Rubidium 87 Rb atoms prepared in an internal ground state, which we choose as |g ≡ |F = 2, m F = 2 . The atoms are trapped in a 2D square optical lattice in the xz-plane created by two pairs of counter propagating laser beams of wavelength λ and wave vector k = 2π/λ, and strongly confined in the y-direction by an additional laser. Note that by tilting the laser beams by an angle α we can adjust the lattice spacing in the xz-plane to any value a = λ/[2 sin(α/2)] ≥ λ/2 [69] (see below). Quantum tunneling allows the atoms to hop between different lattice sites, thus realizing the kinetic energy term with hopping amplitude J h of the single band Hubbard model (4). Furthermore, we work in the hardcore boson limit, i.e. U J h , in which multiple occupancy in a single site is energetically prohibited.
As already discussed in the context of Fig. 1(b) we want to distinguish between and sites in the 2D lattice. This bipartite labeling of the optical lattice is essential for realizing the complex interaction patternṼ ,Ṽ andṼ discussed in Sec. II B, which underlies the second term of the extended Bose Hubbard Hamiltonian (4). In our scheme, we assume that atoms on lattice sites are excited by laser light to the Rydberg 2 P 3/2 -state |r = |n 2 P 3/2 , m = 3/2 z , whereas atoms at sites are excited to |r = |n 2 P 3/2 , m = 3/2 x . Here the subscripts x and z indicate the different local quantization axes for the and sites. Note that the Rydberg states of interest are the stretched states of the fine structure manifold, i.e.
(ii) (i) Figure 4. We consider 87 Rb atoms loaded in a square optical lattice with lattice spacing a and lattice sites labeled alternating as and . Additional AC Stark lasers (magenta and light blue arrows) with wave vectors k 1 = ±2π/λ ACẑ and k 2 = ±2π/λ ACx , respectively, form two pairs of standing waves each with periodicity b = λ AC /2, which are rotated by 45 degrees with respect to the initial lattice. In order to create local quantization axes alongẑ orx for or lattice sites, respectively, we require that that atoms located a lattice site only feel the intensity maxima of the light blue laser with k 1 ∼ z, while lattice site only feel the intensity maxima of the magenta laser with k 2 ∼x. This can be achieved by adjusting the initial trapping lattice by tilting the corresponding trapping lasers with an angle α such that a = λ/[2 sin(α/2)] ≥ λ/2 [69] in order to fulfill b = √ 2a. The two AC Stark lasers have a polarization σ + and resonantly couple the n 2 P 3/2 manifold to a lower lying n D 3/2 manifold thereby inducing an AC Stark shift on each Zeeman m-level in the n 2 P 3/2 manifold except for the maximum stretched |n 2 P 3/2 , 3/2 z,x states. This locally isolates the |n 2 P 3/2 , 3/2 z and |n 2 P 3/2 , 3/2 x state at lattice sites and , respectively, in energy by at least E AC (see left and right panels). A global Rydberg laser (dark blue arrow) with detuning ∆ r E AC propagating along the y-direction then selectively admixes the states |n 2 P 3/2 , 3/2 z and |n 2 P 3/2 , 3/2 x at lattice sites and , respectively, to the ground state |g . states with maximum m = 3/2 value for the given angular momentum j = 3/2. We will show in the next section that the van der Waals interactions between these polarized Rydberg p-states realize naturally the complex interaction pattern required by Eq. (6) for quantum spin ice. By weakly admixing these Rydberg states to the atomic ground state with a laser [see Sec. III C], the ground state atoms will inherit these interaction patterns, thus realizing the interaction term in the extended Bose Hubbard Model of Eq. (4), including the constraints enforced by the interactions satisfying Eq. (6).
It is essential in our scheme that we energetically isolate the stretched states |n 2 P 3/2 , m = 3/2 x,z from the other m-states in the given fine structure manifold. This is necessary to protect these states from mixing with other Zeeman m-levels. Such unwanted couplings can be induced by van der Waals interactions (see Sec. III B below), or via the light polarization of the Rydberg laser. This energetic protection requires an (effective) local magnetic field, which for the and sites points in the x and z-direction, respectively. Strong local fields with spatial resolution on the scale given by the lattice spacing can be obtained via AC Stark shifts, combining the m-dependence of atomic AC tensor polarizabilities with spatially varying polarization gradients. Fig. 4 outlines a scheme, where we superimpose two pairs of counter propagating laser beams of wavelength λ AC (light blue and magenta arrows). They create a standing wave pattern (light blue and magenta gradients), such that lattice sites only see the intensity maxima of the standing wave propagating along the z-direction (light blue laser), while lattice sites see the intensity maxima of the standing wave propagating along the x-direction (magenta laser).
The AC Stark lasers have polarization σ + and resonantly couple the nP 3/2 manifold to a lower lying n D 3/2 manifold (see magenta and blue arrows in the left and right panels of Fig. 4, respectively). This induces an AC-Stark shift on each Zeeman m-level in the n 2 P 3/2 manifold. The Rabi frequency is proportional to Ω AC ∼ (m − 3 2 )(m + 5 2 ) nP 3/2 ||r||n D 3/2 E with E the electric field strength of the AC Stark lasers. In this configuration the stretched states |n 2 P 3/2 , m = 3/2 z,x of interest are not affected by the AC Stark lasers. The minimum shift (as a function of m 3/2) is denoted E AC , which has to obey E AC V off and E AC ∆ r in order to suppress mixing between different m-states due to van der Waals interactions and the excitation laser. Here, V off is the largest off-diagonal van der Waals matrix element in the n 2 P 3/2 manifold (see App. C).
The AC Stark lasers will create an additional trapping potential, V AC (r i )|g g| i , for ground state atoms with minima not commensurate with the initial trapping lattice. In order to not distort the desired lattice structure this additional potential must not be larger than the initial lattice potential, see App. A.
It is then possible to dress the ground state atoms with either the |r = |n 2 P 3/2 , m = 3/2 z or the |r = |n 2 P 3/2 , m = 3/2 x Rydberg state by a single, global laser with Rabi frequency Ω r and detuning ∆ r propagating in the direction perpendicular to the plane, i.e. k r ∼ y (dark blue arrow in Fig. 4). In the local x-and z-basis this laser will couple to all four |m z,x levels with different weights (see App. B). Since the states |m 3/2 z,x are energetically separated by at least E AC from the |m = 3/2  6 , between a pair of 87 Rb atoms in the |3/2, 3/2 z ≡ |n 2 P 3/2 , 3/2 z ⊗ |n 2 P 3/2 , 3/2 z state (solid lines) and in the |1/2, 1/2 z ≡ |n 2 P 3/2 , 1/2 z ⊗ |n 2 P 3/2 , 1/2 z state (dashed lines). We plot the angular part of the rescaled interaction energy, A (n) m 1 ,m 2 as a function of the angle ϑ for various values of the principal quantum number n in atomic units. Here, A (n) 3/2,3/2 (solid lines) corresponds to the angular part of the interaction, V , between two atoms excited to the |r Rydberg state and shows a characteristic ∼ sin 4 ϑ shape due to the dominant S 1/2 -channel. Residual interactions at ϑ = 0 are very small and arise from channels coupling to virtual D-states. (see Tab. I). (b) The angular characteristic of the interaction between two atoms both in the "stretched" Rydberg state, |r = |n 2 P 3/2 , 3/2 z = |np, 1 z | 1 2 1 2 z , can be qualitatively understood from its angular part |np, 1 z and the dominating S -channel: atom A prepared in a |np, 1 z state can make a virtual transition to a lower-lying |ns, 0 z state (red arrow in the lower left panel) while emitting a photon. If this photon propagates along the z-direction it has polarization σ + and cannot be absorbed by atom B. Therefore, atom A and atom B will not interact, i.e. V (ϑ = 0) = 0. If this photon propagates alone the x-direction it is linear polarized with a polarization vector along the y-direction. In the frame of atom C this photon will drive both σ + -and σ − -transitions and thus can be absorbed. Hence, atom A can interact with atom C, i.e. V (ϑ = π/2) 0. state a laser with detuning ∆ r E AC and wave vector k ∼ y will selectively admix the states |3/2 z and |3/2 x at lattice sites and , respectively, to the ground state |g with an effective Rabi frequency Ω r = Ω r /(2 √ 2). The single particle Hamiltonian describing the laser dressing in a frame rotating with the laser frequency for an atom i then becomes where α i ∈ { , } depends on the lattice site of the i-th atom.
In the weakly-dressing regime, ∆ r Ω r , the new dressed ground states are | i ≡ |g i + Ω r /(2∆ r )|r i or | i ≡ |g i + Ω r /(2∆ r )|r i if α i = or , respectively. Thus, each ground state atom gets a small admixture of one of the Rydberg states, depending on the sublattice. Due to the weak admixture of the Rydberg states, the dressed ground states get a comparatively small decay rateΓ = (Ω r /2∆ r ) 2 Γ, where Γ is the decay rate of the bare Rydberg state, which has to be much smaller than the relevant system energy scales discussed below.

B. Interactions between p-states
Below we will consider the van der Waals interactions, V , V and V , between pairs of atoms prepared in the bare Rydberg states |r = |n 2 P 3/2 , m = 3/2 z and |r = |n 2 P 3/2 , m = 3/2 x . For Rubidium atoms excited to Rydberg p-states, these van der Waals forces are strongly anisotropic [70][71][72][73]. Fig. 5(a) shows the angular part of the van der Waals interaction, V , for different n-values, which is in very good approximation proportional to while the actual strength depends on the principal quantum number n and scales as n 11 away from the Förster resonance at n = 38. Similarly, one finds for the interaction between which can be obtained by rotating the coordinate system by π/2. Mixed interactions such as are shown in Fig. 16 (App. D) and have two asymmetric maxima at ϑ = ±π/4. The Rydberg states |r and |r therefore realize the desired angular interaction properties as discussed in Sec. II B. Together with the possibility of creating softcore potentials (see the following subsection), the anisotropy of these interactions leads naturally to the desired interaction pattern illustrated in Fig. 1(b) and demanded by Eq. (6). These interactions underly our realization of the Bose Hubbard Hamiltonian (4). We now detail the physical mechanism which generates these anisotropic interactions, and describe how to derive the aforementioned results. Van der Waals interactions between two atoms i and j prepared in a given Rydberg state arise from the exchange of virtual photons: atom i in a Rydberg state |r i can for example virtually undergo a dipole allowed transition to a lower-lying electronic state |α while emitting a photon. If this virtual photon reaches atom j during its lifetime, it can excite the second atom to an electronic state |β . This then leads to correlated oscillations of instantaneously induced dipoles in both atoms which give rise to the non-retarded van der  Table I. Two atoms both in the |r = |n 2 P 3/2 , 3/2 z = |np, 1 z | 1 2 1 2 z state can couple to six channels. Each channel ν has a characteristic angular dependency (D ν ) 3 2 3 2 which contributes with weight C (ν) 6 . The total interaction can be obtained by summing over all channels, i.e.
2 /r 6 . It turns out that two atoms in the |r Rydberg state dominantly couple to the S 1/2 + S 1/2 channel with a characteristic angular dependence ∼ sin 4 ϑ.
In the case of Rydberg p-states, the angular distribution of these emission and absorption processes of virtual photons in combination with the angular momentum structure of the atomic orbitals leads to nontrivial anisotropic van der Waals interactions. We now focus on the van der Waals interaction V , between both atoms in the |r = |n 2 P 3/2 , m = 3/2 z Rydberg state with quantization axis along the z-direction. Mixed interactions, V , and interactions between both atoms in the |r = |n 2 P 3/2 , m = 3/2 x Rydberg state, V , will be derived in App. C. The latter can simply be obtained by rotating the xz-plane by 90 degrees, i.e. V (r, ϑ) = V (r, ϑ − π/2), while in order to calculate mixed interactions one has to calculate off-diagonal matrix elements in the n 2 P 3/2 manifold.
In general, van der Waals interactions arise as a second order process from dipole-dipole interaction, where V dd couples the initial Rydberg states |r i , r j to virtual intermediate states |α, β and back. Here, d (i) is the dipole operator of the i-th atom and r = rn = (r, ϑ, ϕ) is the relative distance between atom i and atom j with n a unit vector and (r, ϑ, ϕ) the spherical coordinates. It is convenient to rewrite the latter expression in a spherical basis [75] m 1 ,m 2 ;M the Clebsch-Gordan coefficients, and Y m l the spherical harmonics.
Due to the dipole selection rules, states in the n 2 P 3/2 manifold can only couple to states in a n S 1/2 , n D 5/2 or n D 3/2 manifold. It turns out that for 87 Rb the dominating channel is P 3/2 + P 3/2 −→ S 1/2 + S 1/2 , which can be explicitly seen from Table I for various n levels. In order to simplify the following discussion we will first focus on this channel and neglect all other channels including D 3/2 and D 5/2 states which lead to small imperfections discussed in App. C.
For a single atom, the |n 2 P 3/2 , 3 2 state is a stretched-state which reads in the uncoupled basis |n 2 P 3/2 , 3 2 ≡ |np, 1 ⊗ | 1 2 , 1 2 . Thus, it can be factorized into an angular and a spin degree of freedom. Since the dipole-dipole interaction,V dd (r), does not couple spin degrees of freedom, the angular dependence of the van der Waals interaction is determined solely by the angular part of the wave function, which is |np, 1 . Figure 5(b) illustrates the interaction between two atom initially prepared in a |np, 1 state as a function of ϑ for ϑ = 0 (atom A and B) and ϑ = π/2 (atom A and C). We first consider atom A in the lower left corner of Fig. 5(b). Initially prepared in a |np, 1 state it can make a virtual transition to a lower-lying |ns, 0 state [red arrow in the lower left panel of Fig. 5(b)] while emitting a photon. The corresponding angular distribution of the spontaneously emitted photon has the same characteristic as light emitted by a classical dipole tracing out a circular trajectory in the x-y plane [75]. In general, it is elliptically polarized with cylindrical symmetry, but in particular there are two specific directions: (i) light emitted along the z direction (ϑ = 0) is circularly polarized, rotating in the same way as the dipole. Thus, a photon emitted in the z direction has polarization σ + and carries one unit of angular momentum such that the total angular momentum of the combined system atom-photon is conserved. A second atom, labeled as atom B in Fig. 5(b), located on the z axis (ϑ = 0) cannot absorb this photon [see red arrow in the upper left panel of Fig. 5(b)], since only a |n s, 0 state is available. The same result can be derived from Eq. (11), which for ϑ = 0 simplifies tô and couples only states with initial magnetic quantum number m 1 , m 2 to states with m 1 ± 1, m 2 ∓ 1, such that the total angular momentum M = m 1 + m 2 is conserved. Therefore, the dipole-dipole matrix element vanishes, np1, np1|V (AB) dd (ϑ = 0)|ns0, n s0 = 0, and hence atoms A and B do not interact.
(ii) light emitted into the x-y plane (ϑ = π/2) is linearly polarized, with a polarization vector lying in the x-y plane and perpendicular to the emission direction. A third atom, labeled as atom C in Fig. 5(b), located on the x axis is able to absorb this linearly polarized photon emitted by atom A, which in the frame of atom C corresponds to a superposition of σ + and σ − polarized light [see red arrows in the right panel of Fig. 5(b)]. For ϑ = π/2, Eq. (11) contains a sum over cos[ π 2 (µ + ν)], and the only non-vanishing combinations for µ = −1 are ν = ±1. Thus, the dipole matrix element np1, np1|V (AC) dd (ϑ = π/2)|ns0, n s0 is non zero and atom A and C will interact.
In general, for the dominant channel P 3/2 + P 3/2 −→ S 1/2 + S 1/2 only the term d (1) −1 d (2) −1 with µ = ν = −1 in Eq. (11) can contribute to the dipole-dipole matrix element and thus the van der Waals interaction between both atoms in a n 2 P 3/2 , m = 3/2 states becomes V (n) 2 ) 2 ∼ sin 4 ϑ for this channel. Residual interactions at ϑ = 0 and π come from couplings to D 3/2,5/2 channels which are small, see Tab. I. Thus, using Rubidium n 2 P 3/2 states with m j = 3/2 allows an almost perfect realization of an anisotropic interaction with vanishing interaction along one axis and large interaction along a perpendicular axis. Note that interactions between two atoms in a |n 2 P 3/2 , m = 3/2 are negative (attractive) for n > 38 and positive (repulsive) for n < 38, where a Förster resonance at 38P 3/2 + 38P 3/2 −→ 38S 1/2 + 39S 1/2 changes the sign of the interaction [70][71][72][73]. Figure 5(a) shows the result of the the full calculation of the van der Waals interactions between n 2 P 3/2 , 3/2 states including all channels and summing over n and n levels between n ± 10. The full calculation agrees well with the simplified picture discussed above and illustrated in Fig. 5(b) since the dominating channel is the one coupling to S 1/2 states.

C. Soft-core potentials
In the previous section we showed how to engineer the anisotropic part of the interactions required by Eq. (6). We now discuss how to create soft-core potentials by weakly admixing the Rydberg state p-states to the atomic ground state. Following the s-state case [55], this leads to an effective interaction between dressed ground state atoms,Ṽ i j , with a softcore shape and an anisotropic plateau radius. This guarantees that interactions between atoms sitting on a square lattice at different distances a and √ 2a experience the same interaction potential [76], as required by Eq. (6) and illustrated in Fig. 1(b).
The single atom configuration we have in mind was introduced in Sec. III A and is governed by the Hamiltonian of Eq. (7). Pairwise interactions [see panel (a) of Fig. 6] between N atoms both excited to the Rydberg states |r i |r j are described by where V i j (r i j ) = A i j (ϑ i j , ϕ i j )/r 6 i j is the van-der-Waals interaction potential between the Rydberg states of atom i and atom j discussed in the previous section and (r i j , ϑ i j , ϕ i j ) are the spherical coordinates of the relative vector. In the dressing limit, Ω r ∆ r , atoms initially in their electronic ground states |g 1 . . . |g N are off-resonantly coupled to the Rydberg states |r 1 . . . |r N . As a consequence, the new dressed ground states inherit a tunable fraction of the Rydberg interaction. The effective interaction potential between N atoms in their dressed ground states, |g 1 . . . |g N , can be obtained by diagonalizing the Hamiltonian for a fixed relative position and zero kinetic energy. The total Hamiltonian H has block structure where H n governs the dynamics in the subspace with n-Rydberg excitations present (see App. E), while the Ω n matrices describe the coupling between adjacent sectors n and n − 1 due to the laser. Only subspaces H n≥2 contain the interaction potentials V i j since we assume that ground and Rydberg states do not significantly interact. Here, H 0 described the dynamics within the ground state manifold. In the following we will use Brillouin-Wigner perturbation theory (see App. E) based on the small parameter Ω r /(2∆ r ) 1 in order to derive the effective potential between dressed ground state atoms. For red detunings, ∆ r < 0, and repulsive Rydberg interactions V i j > 0, one finds up to fourth order in Ω r /(2∆ r ) for the position dependent energy shift of the dressed ground states a sum of binary interactions of the form with being the anisotropic Condon radius.
In the case of anisotropic Rydberg interactions, the Condon radius depends on the angular pattern of the van der Waals interaction A(ϑ, ϕ) which can be tuned by choosing a particular Rydberg state. Additionally, the Condon radius can be scaled by changing the detuning ∆ r of the dressing laser. Figure 6(b) shows typical examples of the dressed ground state potential,Ṽ i j , for different Condon radii, r c . For large distances, r r c , the dressed ground state potential is proportional to the Rydberg interaction,Ṽ i j = Ω 4 r /(2∆ r ) 4 V i j ∼ 1/r 6 i j , reduced by a factor [Ω r /(2∆ r )] 4 arising from the small probability to excite the atomic ground state to the subspace of two atoms in the Rydberg state, governed by H 2 . However, for small distances, r < r c , when two atoms are within the Condon radius, the excitation to the Rydberg states becomes ineffective due to the large total detuning |∆ r | + V i j (Dipole blockade), and the effective ground state interaction,Ṽ i j ≈Ṽ 0 [1 − (r/r c ) 6 ] for r < r c , saturates at a constant valuẽ which is independent of the strength or form of the Rydberg-Rydberg interactions [52][53][54][55][56]. The resulting specific step-like form ofṼ i j is shown in Fig. 6(b). The presence of a plateau at short distances, r < r c and a rapid decrease of the potential at r ∼ r c , whereṼ i j ∼ 1/r 6 i j , allows to engineer approximately equal interactions between atoms within r c independent of their specific distance. At the same time, long-range interactions, such as next-nearest-neighbor (NNN), are substantially suppressed.
In the following we want to combine the effective steplike ground state potentials of Eq. (16) with anisotropic interactions discussed in the previous section. Panel (c) of Fig. 6 shows a 2D contour plot ofṼ i j /Ṽ 0 in the x-z plane for 87 Rb ground state atoms dressed with the Rydberg state |r = |32 2 P 3/2 , m j = 3/2 z . The dashed yellow equipotential lines have the form of a figure-eight which follows from the V i j ∼ sin 4 ϑ dependency of the Rydberg interaction, discussed in Sec. III B. Residual interactions along the z-axis (ϑ = 0) comes from channels coupling to D-states (see Tab. I) and depend on the principal quantum number n as shown in Fig. 7.
Once defined on the top of an underlying lattice, the combination of anisotropic interactions between Rydberg p-states and step-like potentials via dressing techniques allows to design interaction potentials as the ones shown in Fig. 6(c). Here, atoms (black circles) regularly arranged on a square lattice with lattice spacing a interact with J, J and J along the ± 45-degree lines, the x-axis and the z-axis, respectively. It is possible to tune the interaction strength J, J and J over a large range by e.g. changing the detuning, ∆ r , or the principal quantum number, n, of the Rydberg state. In particular one can realize an interaction pattern where atoms sitting at different distances a and √ 2a interact with equal strength, that is J ≈ J , while J J, thus realizing a frustrated J − J model. Note that the interaction symmetry in this case is triangular on top of an square lattice.

D. Explicit numbers and discussion of imperfections
As an explicit example we consider the 29 2 P 3/2 Rydberg manifold of 87 Rb. We resonantly couple the 29 2 P 3/2 manifold to the lower-lying 7 2 D 3/2 manifold, as illustrated in Fig. 4 (n = 29 and n = 7), with a laser of wavelength λ AC = 3.296 µm. This results in a lattice spacing a = λ AC /(2 √ 2) = 1.16 µm which can be adjusted by tilting the trapping lasers by an angle α = 39 degrees (see Sec. III A).
Adjusting the detuning ∆ r of the Rydberg laser allows to tune the length scale and the imperfections in Eq. (16). These are (i) small long-range interactions between nearest-neighbor lattice sites and (ii) deviations form the constraint model of Eq. (6). Here, for example we use ∆ r = 2π × 400 kHz which yields the following interaction pattern between particles labeled in Fig. 6 By varying the Rabi frequency of the Rydberg laser Ω r = 2π × (80, 120, 160) kHz one obtains = Ω r /2∆ r = (0.10, 0, 15, 0.20) which gives rise to an effective ground state interactionṼ 0 = Ω 4 r /8∆ 3 r = 2π × (80, 410, 1290) Hz. This is much larger than the effective decay rate from the dressed ground stateΓ = 2 Γ = 2π × (33, 75, 133) Hz, and larger than a corresponding tunneling rate between the minima. Here, Γ = 2π × 3.3 kHz is the decay rate form the Rydberg states.
There is an ample choice in the parameter regimes available as a function of the n-level. Away from the Förster resonance at n = 38, it is possible to engineer infra-red lattices which allow for comparable timescales between the interactions induced by the dressing, and the tunneling matrix elements of the atoms on the original square lattice. Going higher in n, closer to the Förster resonance, allows faster timescales and slower decays: however, in this case the infra-red laser has a strong influence on the underlying lattice, excluding the possibility of using conventional single particle tunneling to induce quantum fluctuations. On the other hand, one can profit here from the richness of the Rydberg manifolds involved, realizing the hopping matrix element as a spin-exchange coupling between different atoms sitting at different potential minima [77]. In both cases above, the interaction pattern will depend on the specific targeted n, as discussed in Sec. III B.
As the qualitative (and in many respects quantitative, as indicated in Table I and Fig. 7) shape of the interactions will be very similar in the interval of interest n = 25 − 37, we will focus in the following on a single case sample to underpin the stability of the many-body effects we are interested in.

IV. NUMERICAL RESULTS
In this section, we consider the properties of the approximate realization of the quantum spin ice model Hamiltonian proposed above. We demonstrate that, as a function of the strength of the quantum dynamics, the ground state has two regimes reflecting two distinct forms of ordering (Sec. IV B). One, stabilized via a quantum order by disorder mechanism, generates the above mentioned plaquette phase for sufficiently strong quantum dynamics. As it is weakened, there is a transition into a phase with classical ordering, which is stabilized by the long-range parts of the dipolar couplings and which breaks translational symmetry in a different way. In addition, we show that even without quantum dynamics, there is an interesting thermal phase transition to an approximate realization of a (classical) Coulomb phase, with only a very small density of defects (plaquettes violating the ice rule) of around 5% (Sec. IV C 2). We discuss signatures of these items in various quantities, in particular proposing a simple correlation function in which the quantum plaquette order will be visible, and which should be accessible in cold atom experiments via in-situ parity measurements [33,78,79].

A. General definitions and conventions
We begin with some general definitions and technical details of our numerical study. We consider both the unconstrained spin-1/2 model H from by Eq. (5), as well as the projected model H 2 inside the spin ice manifold: Here i j denote nearest-neighbor (NN) sites on the 2D checkerboard lattice and (i jkl) label the four sites around empty square plaquettes. In our Exact Diagonalizations (ED) we have considered finite-size clusters with periodic boundary conditions and N=16, 32, 36, 64 and 72 sites, see details in App. F. To treat these clusters with ED, we exploit translational symmetry, point group operations (the model has C 2v symmetry), as well as spin inversion (S z → −S z ) inside the total magnetization sector S z = 0. Consequently, the eigenstates are labeled by linear momentum k, the irreducible representations of the point group of k, and the parity under spin inversion. We note that, whereas the quantum phase is quite robust, the classical phase is considerably less so, reflecting the many nearly-degenerate classical ice states. We illustrate this in App. G by imposing a variable cut-off on the long-range aspect of the dipolar couplings J i j : by neglecting terms weaker than a cutoff J c , we find a set of states with different classical orders, which settle down into the correct ground state without truncation for J c no larger than 0.001.
B. The two zero-temperature phases: Low-energy spectroscopy and ground state diagnostics Figure 8 shows the low energy spectra of H as a function of J ⊥ for N=32 (a) and N=36 (b), and that of H 2 as a function of t for N=64 (c). All spectra correspond to the total magnetization sector S z = 0 and a cutoff value of J c = 0.001. In all spectra, there is a manifold of low-lying states that is well separated from higher-energy excitations. Provided they become degenerate in the thermodynamic limit, these states are the finite-size fingerprints of the spontaneously symmetry broken phases [80][81][82][83][84]: their multiplicities and symmetry content reveal the nature of the ground state. The structure of the low-lying energy states show consistently two qualitatively different phases. One, which is adiabatically connected to the classical limit J ⊥ = 0, and the other which is stabilized for large enough J ⊥ or t.
We begin with the classical phase, focusing on the N=32 (a) and N=64 (c) results first. Here we find four low-lying states which become exactly degenerate as J ⊥ → 0. We find translational symmetry breaking with ordering wavevector Q = (− π 2 , π 2 ), as illustrated in Fig. 9(d). The nature of this phase is revealed by the spin-spin correlation profiles of Figs. 9(a-b), with alternating up-down spins along one of the two diagonal directions of the lattice. The vanishing of cor-relations on every second diagonal line arises due to the existence of two states compatible with the non vanishing correlations on the other diagonals. For a finite cluster, these appear with equal weight and thus average out, while in the thermodynamic limit, symmetry breaking selects either one of the two spontaneously. Finally, the N=36 cluster cannot accommodate the Q = (− π 2 , π 2 ) phase (see App. G), which is why the low-lying sector of Fig. 8(b) has a different structure (and, in fact, higher ground state energy per site, see Table II).
Turning to the quantum phase, the N=32 and 64-site spectra give the onset of this phase around J ⊥ 0.23 and t 0.28, respectively. Beyond this point, the spin structure factor (not shown) is completely structureless, indicative of the absence of magnetic (classical) ordering. Since the imperfections in the present spin-ice model are expected to become irrelevant for large enough J ⊥ , this phase must be the plaquette phase of the pure spin-ice model [11,67] and the pure Heisenberg model [85]. The standard diagnostic for this phase is the dimer-dimer (or energy-energy) correlations, and indeed the correlation profiles of Fig. 9(c) show a strong Q = (π, π) response within one sublattice of empty plaquettes. This is consistent with the structure of the low-lying spectra which show two low-lying states with momenta k = 0 and (π, π) which come almost on top of each other for N=64, see Fig. 8(c). Note that for N=32, there is a third low-lying state (with k = 0) which is however not related to the physics at the thermodynamic limit but it is specific to the special topology of this cluster [86].
Further information about the two phases is given in Fig. 10, which shows the GS expectation values of the longitudinal and transverse NN spin-spin correlations for all symmetry-inequivalent bonds, as well as the square magneti- zation of crossed plaquettes. The former describe how the energy is distributed over the bonds and over the different directions in spin space, while the latter is a measure of the admixture from states outside the spin ice manifold. First, the NN correlations show that the spins fluctuate mostly along the zaxis for small J ⊥ , as expected. More importantly, most of the energy comes from antiferromagnetic bonds along one of the two diagonal directions (bonds labeled 's 1 c 2 ' in the inset of Fig. 10), which is a clear signature of the presence of strongly asymmetric spin-spin correlations in this regime. This asymmetry, which is inherited by the point group symmetry (C 2v ) of the model, is more directly revealed in the spin structure factor discussed above. Second, the qualitative change in the behavior of the NN correlations around J ⊥ ∼ 0.2, reflects the presence of the phase transition in this region. Finally, the square of the total magnetization per crossed plaquette reveals that the spin-ice manifold remains well protected up to relatively high J ⊥ .
C. Further insights in the classical limit J ⊥ = 0

Momentum space minimization
The nature of the classical phase and the role of the dipolar couplings can be understood in more detail by a closer examination of the limit J ⊥ = 0 using a classical minimization treatment in momentum space [87][88][89][90]. The checkerboard lattice has a square Bravais lattice with two sites per unit cell. In the following, sites are labeled as i → (R, α), where R gives the position of the unit cell, and α = 1-2. For J ⊥ = 0, we can replace S z i → 1 2 σ i , where σ i = ±1. The total energy then reads Using σ R,α = 1 √ N uc k e ik·R σ k,α , where N uc = N/2 is the number of unit cells, and J Rα,R α = J R−R ,αα (from translational   invariance), yields where the 2 × 2 interaction matrix Λ(k) is given by Let us denote by λ 1,2 (k) and v 1,2 (k) the eigenvalues and the corresponding (normalized) eigenvectors of Λ(k), with λ 1 (k) ≤ λ 2 (k). Minimizing λ 1 (k) over the entire BZ of the model provides a lower bound for the energy [87][88][89][90]. The corresponding eigenvector is a faithful ground state provided it satisfies the spin length constraint at all sites. The minimization can be done both for the infinite lattice and for the finite lattices studied by ED by simply scanning through the allowed momenta of each cluster. The latter are discussed in App. G and are useful for clarifying the various finite-size effects in our ED data. Here we focus on the infinite lattice case. Figure 11 shows the momentum dependence of the low-energy branch λ 1 (k) for cutoff values J c =0.3, 0.1, 0.01 and 0.001. For J c = 0.3, which amounts to keeping only the dominant, NN couplings (i.e. three couplings per site), the minimum sits at Q = (π, π) and corresponds to the wellknown Néel phase with AFM correlations along both the horizontal and the vertical directions of the lattice. This phase is stabilized by the imbalance in the NN imperfections, which favors the first two vertex configurations in Fig. 2(b). However, further-neighbor interactions destabilize the Néel phase and lead to a different minimum. For J c = 0.1, which amounts to keeping seven interactions per site, the minimum of λ 1 (k) now sits at the two M points of the BZ, Q = (π, 0) and (0, π), which correspond to a stripy AFM alignment of the spins in the horizontal or the vertical direction of the lattice.
Lowering J c further shifts the minimum to two incommensurate (IC) positions, ±Q IC , which are extremely close to the commensurate ±(−π/2, π/2) points. For example, for J c = 0.01 (12 interactions/site), J c = 0.001 (31 interactions/site) and J c = 10 −6 (299 couplings/site), the minima sit respectively at Q IC = 0.473(−π, π), 0.457(−π, π) and 0.462(−π, π). At the same time, the corresponding eigenvector v 1 (Q IC ) cannot be used to construct a state satisfying the spin length constraint at all sites of the system simultaneously. This means that the present method cannot deliver the true ground state of the system and that λ 1 (Q IC ) serves only as a lower energy bound. Physically, the system may accommodate the tendency for incommensurate correlations by forming long-wavelength modulations of the local (−π/2, π/2) order parameter, in analogy e.g. to the anisotropic Ising model with competing interactions (the so-called ANNNI model) [91][92][93][94]. We should remark however that the energy landscape around the IC minimum is very flat and its distance from (− π 2 , π 2 ) is very small, so in principle such discommensurations (if any) should appear at much longer distances than the ones considered in our finite-lattice calculations, and indeed the length scales over which cold atom realizations are uniform (in both density and interaction patterns) on account of the parabolic confining potential. To confirm this point we have performed classical Monte Carlo (CMC) simulations on N = 2 × L × L-site clusters with periodic boundary conditions, see details in App. H. All results up to L = 36 give consistently the (− π 2 , π 2 ) state without any sign of domain-wall discommensurations, implying that at least for these distances the system locks-in to the closest commensurate Q = (−π/2, π/2) phase.

Thermal phase transition into a classical Coulomb phase
Given the finite energy gap above the commensurate Q = (− π 2 , π 2 ) state at J ⊥ = 0 (see Fig. 8), one expects that this phase survives against thermal fluctuations up to a finite temperature T C . To confirm this picture, and to find the numerical value of T C , we have performed classical MC simulations at finite temperatures. The first two panels of Fig. 12 show the T-dependence of the structure factor S(Q) at Q = (−π/2, π/2), and the specific heat per site for systems up to N=2x28x28 sites. The results demonstrate clearly the thermal phase transition, with T C 0.185. The third panel shows the Tdependence of the three different types of crossed plaquette configurations: the ice-rule 2in-2out states, and the defected 3up-1down (or 3down-1up) states and 4-in (or 4-out) states. The defects are almost entirely of the 3up-1down type, but their density remains very small up to the transition temperature (about 5%). So the classical phase gives way to a Coulomb gas [61,62], an approximate realization of a classical Coulomb phase, with only a very small density of defects. This phase is marginally confined on account of the logarithmic nature of the interactions between the defects; given the non-vanishing defect density above T C , their correlations are expected to exhibit a screened (Debye) form [95]. As we have shown in the spin ice example, weakly Rydberg-dressed atoms in optical lattices provide a perfect platform to investigate quantum magnetism in AMO settings. From the one hand, the starting building block being XXZ models allows to use the zz component to impose constraints, and to employ the exchange terms to generate the dynamics in perturbation theory. This is generally the best way to get effective many-body interactions from simple two-or one-body ones: the desired terms emerge at the lowest non-vanishing order in the expansion. From the other hand, the fact that interactions decay like van der Waals beyond a set radius makes imperfections intrinsically local in any dimension, circumventing additional dimensional effects which can emerge with power law interactions with a slower decay.
The possibility of engineering the spin ice dynamics discussed above resulted from the combination of two features: (i) a simple lattice structure, a square lattice, on which the spins are arranged on, and (ii) a complicated interaction pattern resulting from the condition of equal zz interaction strength around each vertex -requiring both anisotropic and plateau-like interactions. In AMO settings, a second strategy can be pursued, where simple interaction patterns are combined with exotic lattice geometries to induce emergent gauge invariance. This way, the complication of realizing fine-tuned pattern is transferred to a complicated lattice geometry, which might be realized provided the correspondent light-pattern is realizable. Below, we briefly summarize the key features of this second strategy, and discuss a specific application thereof.
The first ingredient is simple interactions, which are in general repulsive and have an isotropic plateau-like structure. Within the context of Rydberg atoms, they are usually present in two cases. The first one are ground states atoms coupled via a two-photon transition to an s-state, whose interactions are isotropic in the full 3D space up to corrections of few percent [46][47][48][49][70][71][72][73]. The corresponding frozen regimes have already been accessed in a series of experiments [23][24][25][26][27][28][29][30][31][32][33][34][35][36][37]96]. A second way of realizing isotropic interaction is to consider atoms trapped on a 2D plane, and dressed to a p-state polar-ized along the direction orthogonal to the plane itself. This way, differently with respect to the case discussed above, the in-plane interaction is basically isotropic -virtual exchange of photons within the Rydberg manifold is always allowed.
The second key ingredient is a set of complicated lattices where the aforementioned interactions, once considered as a sharp plateau ones, are sufficient to define a classical limit where there is a set of degenerate classical ground states, increasing extensively with system size. Given the shape of the interactions, one can identify the possible lattices as follows. First, we define as b the largest distance between sites belonging to the same simplices [8] (or gauge cell), that is, the unit cell where the Gauss law of interest is defined (in the spin ice case, these are the squares with crosses). Secondly, we define as c the smallest distance between sites which do not belong to the same gauge cell. It is then clear that, in the case where b < c, a plateau interaction of range b < r bc < c can generates the desired constraint in each gauge cell. In case this is not true (like, e.g., in the square ice case), additional features are needed such as angular dependence.
Some examples of the lattices which satisfy the previous property are illustrated in Fig. 13, together with the corresponding gauge cells. The list includes several 2D lattices which have already been realized in AMO settings, such as the Kagome lattice with triangular gauge cells (but not with hexagonal gauge cells), the Ruby lattice and the Honeycomb lattice. In general, various lattices can be constructed satisfying the previous properties taking triangular cells as gauge cells. In the 3D case, a simple example is the pyrochlore lattice, already discussed in Ref. [39] in the context of polar gases. Interestingly, in this former case, dipolar interactions behave in a very similar manner to simple plateau-like ones due to their symmetry content [97].
Once taken at the proper filling factor for the underlying Bose-Hubbard Hamiltonian, all of those lattices generate in perturbation theory quantum dimer or quantum loop models [8]. These are naturally described by emergent gauge theories: however, the gauge symmetry itself is not always straightforwardly determined given the gauge symmetry of the microscopic constituents. While in our cases the latter is always of U(1) type (the number of bosons in each gauge cell is preserved along the dynamics), the underlying gauge symmetries can be discrete. This is the case, e.g., of the Kagome Bose-Hubbard model with the hexagonal cell as gauge cell, where the proper low-energy theory is a Z 2 gauge theory which can undergo deconfinement, and thus stabilize a spinliquid phase.
The procedure to derive the proper dimer model dynamics given a lattice and an Ising constraint is outlined in Ref. [8]. Below, we illustrate a simple example of how complicated lattices can meet simple interactions to let a quantum dimer model emerge by focusing on the concrete example of the 4-8 lattice.
B. Emergent quantum dimer dynamics on a 4-8 lattice from an XXZ model The 4-8 lattice (also known as CAVO lattice) [98] represents a useful example to illustrate how the combination of a complicated lattice with simple Ising interactions can lead to intriguing quantum dynamics. The lattice structure for the underlying bosons we start from is the squagome lattice [99,100], illustrated in Fig. 14: once the triangles are identified as the gauge cells, it is easy to see that b = a, c = √ 2a, so that a plateau interaction of range 1 < r/a < √ 2 can indeed enforce constraints on the gauge cells. Since each site is shared by 3 gauge cells, a filling fraction of n = 1/3 atoms per site, combined with the plateau interactions, will generate a degenerate manifold H 4−8 of classical ground states where for each triangle a single site is occupied [see Fig.14(b)].
When formulated in spin language with S z j = n j − 1/2, the Hamiltonian has (trivially) a set of U(1)-like conserved charges at each triangle, that is: Once an additional, small term inducing quantum fluctuations is introduced: tunneling between the different classical degenerate minima becomes possible within perturbation theory, while still preserving the set of conserved charges in Eq. (21). Notice that we used two different matrix elements for particle tunneling around the squares (J) and around the octagons (J). Two kind of moves are allowed. At second order, two particles sitting along a diagonal of a square plaquette can resonantly flip to (a) (b) Figure 13. Lattices possessing the properties discussed in Sec. V A: Panel (a) shows a kagome lattice with triangles as gauge cells (shaded area) and panel (b) a honeycomb lattice with hexagons as gauge cells. In both cases, the maximal intraplaquette euclidean distance (yellow dashed circle) is smaller than the minimal interplaquette distance (gray arrow). The radius of the needed plateau-like interaction is described by the yellow circles.
(a)  23) and (24). Panel (f) shows the optical lattice pattern of Eq. (25) as described in the text. Darker areas correspond to deeper potentials. sit on the other diagonal: where we have numbered the sites of the square plaquette in clockwise order. The next non-vanishing contribution takes place at fourth order, where particles sitting at the edges of each octagonal plaquette can re-arrange via an extended ringexchange: where we have numbered the sites of the octagonal plaquette in clockwise order. The two terms are illustrated in Fig. 14(d)-(e). We now reformulate the problem in terms of dimer models, which allows to set up a proper description in terms of effective degrees of freedom. In order to do that, we follow the procedure exemplified in Ref. [8,[101][102][103], and illustrated in Fig. 14(b)-(c): we define a new lattice, the so called simplex lattice, whose vertices are the middle points of each gauge cell, and whose bonds connect vertices of gauge cells which share a single site: each bond sits on a vertex of the original lattice. Then, we introduce dimer variables on the bonds as follows: i) if a bond sits on a site which is occupied by a boson, we draw a dimer, ii) if not, we leave the bond empty. This way, the Gauss law of Eq. (21) is easily reformulated as a conservation law of a single dimer at each vertex. The lattice on the top of which the quantum dimer model is defined is then a 4-8 lattice: as it is bipartite, the corresponding low-energy theory is a U(1) gauge theory, which can then display different confined phases as a function of the two kinetic energy terms for the dimers H and H . This setup might constitute then a perfect setting for the investigation of the competition between different RVB solid orders and the transitions between them. The corresponding periodic structure can be either realized using digital-micromirror-devices (DMD) [63], or by using an optical potential of the form where V 1 (x, y) = cos(πx) 2 + cos(πy) 2 − 2 cos(0.55) cos(πx) cos(πy) is a 2D lattice created by two 1D standing waves with phase difference φ = 0.55 and anti-parallel polarizations e 1 ·e 2 = −1.
The second 2D lattice is created by lasers with three times the frequency and orthogonal polarization, Both lattices are rotated by 45 degrees, respectively. The full lattice structure is illustrated in Fig. 14(f), and realizes the squagome lattice potential of interest.

VI. CONCLUSIONS AND OUTLOOK
In summary, we have shown how dynamical gauge fields emerging from frustration can be ideally realized in cold atom systems by employing optical lattices combined with Rydberg interactions, allowing to probe gauge theory phenomena in a variety of models. In particular, we analyzed in detail the case of quantum square ice, a paradigmatic example of frustrated statistical mechanics, both at the few-and at the many-body level.
From the atomic physics side, the key element of our implementation is represented by the tunable interaction pattern generated by Rydberg p-states. Here, prominent atomic physics features can be exploited in order to generate (repulsive) anisotropic interactions that allow to enforce the complex gauge constraints of square ice models. The possibility of generating such anisotropic interaction patterns enriches the cold atom Hubbard toolbox of another potential feature, which can find different applications in many-body physics even beyond engineering complicated and fine tuned lattice constraints. A straightforward extension could be the realization of U(1) discrete gauge theories on a cubic lattice, where the Polyakov argument is then circumvented and thus a stable deconfined phase, a gapless spin liquid, can be stabilized. Moreover, the non-trivial interactions within the Rydberg manifold can also be used to construct directly the spin model of interest, thus ensuring larger energy scales with respect to real tunneling dynamics.
From the many-body side, we have provided numerical evidence that typical imperfections generated by the Rydberg interactions still allow the observation of a non-trivial state of matter, a plaquette valence bond crystal. Moreover, we have shown how a cold atom suited detection technique can be identified, by performing parity measurements along the plaquettes, which directly identifies the spontaneous symmetry breaking of a discrete lattice symmetry. An additional point is that, even in the absence of quantum dynamics, the engineered interactions stabilise a magnetically ordered state with a large unit cell at low temperature, which gives way to a classical Coulomb gas, a marginally confining two-dimensional Coulomb phase with a small but nonzero density of charges in the form of thermally activated plaquettes violating the ice rule.
The advanced interaction engineering available within the Rydberg toolbox paves the way toward the realization of different constrained dynamics even beyond quantum ice. In the last section, we discussed some 2D examples of quantum dimer and quantum loop models (quantum link models) that can be realized by combining isotropic plateau-like interactions with exotic lattices. This different route, which replaces the complex interaction pattern of square ice with a complex optical lattice setting, represents an interesting, complementary approach for the microscopic realization of gauge theories, which can benefit from the recent developments of in situ imaging and digital-mirror-device optical lattice techniques.
Different directions can be pursued further following the lines discussed here. A first, interesting extension would be to understand whether different kinds of anisotropic interactions can play a significant role in engineered Ising constraints in cold atom systems. In particular, anisotropic interactions between Rydberg d-states of 87 Rb atoms have been recently demonstrated in Ref. [37]: as their angular dependence dif-fers from the one discussed here, it can constitute yet another tool in order to realize complicated, fine-tuned interaction patterns. Secondly, the present proposal, which generates pure gauge theories, can be combined in a modular way with previous ones [104] in such a way that either fermionic or bosonic matter can be included into the dynamics. Microscopically, this would require an additional fermionic (bosonic) species to be trapped onto a lattice whose minima sits at the centre of each vertex. The combination of gauge fields and dynamical matter fields will then allow to investigate scenarios such as 2D quantum electrodynamics in a quantum link formulation in the fermionic case, or the Fradkin-Shenker scenario of Higgs physics in the bosonic one. Finally, cold atom realizations can also provide a suitable platform for the investigation of dynamical effects in quantum dimer models and gauge theories in general; it'd be interesting to see whether simple observables and experimental procedures can be implemented, to described complex many-body phenomena such as string dynamics [105] in the presence of static charges [40], or the dynamical properties of thermally activated monopoles on top of a vacuum state. The AC Stark lasers introduced in Sec. III A will create an additional trapping potential, V AC (r i )|g g| i , for ground state atoms with minima not commensurate with the initial trapping lattice. In order to not distort the desired lattice structure this additional potential must not be larger than the initial trapping potential. The dominant effect comes from a second order Stark effect by off-resonantelly coupling the 5S state to the first excited state 5P and is given by V AC = Ω 2 5s5p /(2∆ 5s5p ) with Rabi frequency Ω 5s5p = 2d 5s5p E/ and detuning ∆ 5s5p = 2πc(λ −1 5s5p − λ −1 AC ). Here, d 5s5p = 5S |d|5P is the transition dipole matrix element and λ 5s5p the transition wavelength. Fig. 15(a) shows the desired trapping lattice created by two counter propagating laser beams which form a ground state potential V trap (z, x) = cos 2 kz + cos 2 kx. The dashed black lines indicate the 0.9 level lines of the AC-Stark potential The maxima are localized at the and lattice sites, respectively, as required in Sec. III A. Fig. 15(b) shows the total potential, V tot = V trap + αV AC , in the case of (a) (b) Figure 15. Contour plots of the total trapping potential, V tot (x, z).  Fig. 15(a) and (b)]. In the case of equal strength of the trapping lattice and additional lattice created by the AC-Stark lasers the potential minima are still located at the same position, but slightly elongated. Note that the potential barrier between neighboring lattice sites is about 1/2 smaller than without the additional AC Stark laser. This will lead to higher tunneling rates compared to the case without the AC Stark laser.

Appendix B: Global Rydberg laser excitation
In the following we show that it is possible to weakly admix the locally polarized Rydberg states of Sect. III A to the electronic ground state |g using a single laser with a wave vector k ∼ y and polarization σ + (see Fig. 4) In the local x-and z-basis this laser will couple to all four m j -levels with different weights, i.e.
where we used the irreducible representation of a rotation in the j = 3/2 subspace, D[R(α, β, γ)] = e −iαJ z e −iβJ y e −iγJ z , with matrix elements Since the states |m 3/2 z,x are energetically separated by at least E AC from the |m = 3/2 state a laser with detuning ∆ R E AC and wave vector k ∼ y will selectively admix the states |3/2 z and |3/2 z at lattice sites and , respectively, to the ground state |g with an effective Rabi frequency Ω R /(2 √ 2).

Appendix C: Van der Waals interactions
In this appendix we briefly summarize the technical details in order to calculate the angular dependent van der Waals interactions of Sec. III B. Due to the odd parity of the electric dipole operators d (i) µ and d ( j) ν , the dipole-dipole interaction, V dd , of Eq. (11) can only couple states with initial angular (total) momentum ( j) to states with new angular (total) momentum ± 1 ( j or j ± 1). Therefore, the number of possible "channels" n jm 1 +n jm 2 −→ n j m +n j m for which the matrix element n jm 1 ; n jm 2 |V (i j) dd |n j m ; n j m is non-zero are limited. While there is no selection rule for possible final principal quantum numbers n and n which solely determine the overall strength of the matrix element, the dipole-dipole matrix element is only non-zero if the magnetic quantum numbers and the spherical component of the dipole operator fulfill m 1 + µ = m and m 2 + ν = m . If the energy difference δ αβ = E(α) + E(β) − 2E(n j), between the initial states n j and the intermediate states α ≡ n α α j α m α and β ≡ n β β j β m β of the atoms is larger than the dipole-dipole matrix element connecting those states the dominant interaction is of van der Waals type which arises from V dd in second order Here,V vdW is an operator acting in the degenerate manifold of magnetic sublevels withP i j = |n jm i , n jm j n jm i , n jm j | a projector into the n j-manifold andQ α,β = |α, β α, β| a projector on a specific state in the complementary space. The sum is over all two-atom energy levels, where the indices α ≡ n α α j α m α and β ≡ n β β j β m β denote a full set of quantum numbers that specify the states. Due to the electric dipole selection rules discussed above this sum can be split up into channels denoted by ν = ( α , j α ; β , j β ). Eq. (C1) can be written asV vdW = ν C (ν) 6 D ν (ϑ, ϕ)/r 6 , where C (ν) 6 contains the radial part of the matrix elements which accounts for the overall strength of the interaction and is independent of the magnetic quantum numbers. Here, R j i = drr 2 ψ n i , i , j i (r) * r ψ n j , j , j j (r) is the radial integral calculated with radial wave functions ψ n j , j , j j (r) obtained using the model potential from [106]. The matrix on the other hand is a matrix in the subspace of magnetic quantum numbers which contains the relative angles between the two atoms (s = 1/2) As an example we show the D 1 matrix for the first channel P 3/2 + P 3/2 −→ S 1/2 + S 1/2 in the subspace of states | 3 2 , 3 2 , | 3 2 , 1 2 , | 3 2 , − 1 2 and | 3 2 , − 3 2 where the first atom is fixed in the m = 3/2 state. In general one has to diagonalize the operatorV vdW in the degenerate Zeeman subspace in order to obtain the new eigenenergies and eigenstates in the presence of interactions. If an external electric or magnetic field separates an initial two atom state |m 1 , m 2 from all other Zeeman sublevel such that the energy difference is larger than the vdW coupling matrix elements then it is possible to simply take expectation values of V (n) m 1 ,m 2 (r) = m 1 , m 2 |V vdW |m 1 , m 2 in order to obtain the interaction potential of two atoms initially in the |m 1 , m 2 state. The angular dependence of the van der Waals interaction between two Rydberg atoms in a | or in a | state, V (r, ϑ) = V (r, ϑ − π/2) ∼ sin 4 ϑ/r 6 are the same up to a rotation by 90 degrees and show the typical anisotropic behavior discussed in Sec. III B [see solid lines in Fig. 5(a)] On the other hand, the angular dependence of the mixed interactions between two Rydberg atoms in a | state, shown in Fig. 16(a), exhibits two asymmetric maxima at ϑ = ±π/4. The asymmetry arises from off-diagonal matrix elements, e.g.  Fig. 16 shows a contour plot of the mixed interaction,Ṽ /Ṽ 0 of Eq. (16), between the dressed ground state atoms | in the middle and the surrounding | atoms. Interactions with the neighboring | atoms (red solid arrows) are strong, ∼Ṽ 0 , while interactions with nextnearest-neighbor | atoms (red dotted arrows) are strongly suppressed due to the plateau structure of the potential. In this appendix we review Brillouin-Wigner perturbation theory for N atoms in order to obtain the effective ground state potentials of Sec. III C. We are interested in finding the new ground state |G of the Hamiltonian (14) which is adiabatically connected to |G = N i=1 |g i for Ω → 0.
For the following analysis it is convenient to split the dynamics in subspaces containing 0, 1, 2 etc. Rydberg excitations H 0 = {|g 1 , . . . , g N }, H 1 = {|g 1 , . . . , g i−1 , r i , g i+1 , . . . , g N }, H 2 = {|g 1 , . . . , g i−1 , r i , g i+1 , . . . , g j−1 , r j , g j+1 , . . . , g N }, We denote the matrices of size (dim H n ) 2 on the diagonal with H n . They describe the dynamics in the subspace H n with a fixed number of Rydberg states n, while the Ω n matrices of size (dim H n−1 × dim H n ) describe the coupling between adjacent sectors n and n − 1 due to the laser. Note, that only subspaces H n≥2 contain the interaction potentials V i j since we assume that ground and Rydberg states do not significantly interact. E.g. we find H 0 = 0 and H 1 = −∆1 N , where 1 N is the N × N identity matrix. Finally, we split the Hamiltonian into two parts, H 0 and H 1 , accounting for the diagonal terms and the laser coupling, respectively.
In the following we use Brillouin-Wigner perturbation theory to find the energy of the ground state of H up to fourth order in H 1 /H 0 . Therefore, we define the projector P = |G G| and its complement Q = 1 − P for which PH 0 Q = 0. Splitting the Hamiltonian in the corresponding subspaces yields and For the zero order term we find E (0) G = PH 0 P = H 0 = 0. The first order contribution is also zero, since E (1) G = PH 1 P = 0. In order to calculate higher order terms we first need to calculate the resolvent operator in the Q space The second order contribution to the ground state energy is E (2) G = PH 1 QRQH 1 P. Using PH 1 Q = (Ω 1 , 0, 0) , and Ω 1 = Ω1 T , where 1 is a n × 1 vector containing only 1's, we get E (2) G = Ω 1 (−H −1 1 )Ω † 1 = N Ω 2 ∆ . The second order contribution accounts for the single-particle light shift due to the laser, E (2) G /N. The third order contribution to the ground state energy vanishes E (3) G = PH 1 QR QH 1 Q − E (1) G RQH 1 P = 0. For the fourth order contribution we obtain Using the matrix representation in the subspaces yields The first line of the latter equation contains (for the first time) the interaction potential V i j which is present in all terms H n≥2 . The last line can be simplified to −N 2 Ω 4 ∆ 3 , while the explicit form of the first one depends on the number of particles.