Programmable order by disorder effect and underlying phases through dipolar quantum simulators

In this work, we study two different quantum simulators composed of molecules with dipole-dipole interaction through various theoretical and numerical tools. Our first result provides knowledge upon the quantum order by disorder effect of the $S=1/2$ system, which is programmable in a quantum simulator composed of circular Rydberg atoms in the triangular optical lattice with a controllable diagonal anisotropy. When the numbers of up spins and down spins are equal, a set of sub-extensive degenerate ground states is present in the classical limit, composed of continuous strings whose configuration enjoys a large degree of freedom. Adopting the the real space perturbation theory, our calculation demonstrates a lifting of the degeneracy, favoring the stripe configuration. When $J$ becomes larger, we adopt the infinite projected entangled-pair state~(iPEPS) and numerically check the effect of degeneracy lifting. The iPEPS results show that even when the spin exchange coupling is strong the stripe pattern is still favored. Next, we study the dipolar bosonic model with tilted polar angle which can be realized through a quantum simulator composed of cold atomic gas with dipole-dipole interaction in an optical lattice. By placing the atoms in a triangular lattice and tilting the polar angle, the diagonal anisotropy can also be realized in the bosonic system. With our cluster mean-field theory calculation, we provide various phase diagrams with different tilted angles, showing the abundant underlying phases including the supersolid. Our proposal indicates realizable scenarios through quantum simulators in studying the quantum effect as well as extraordinary phases. We believe that our results indicated here can also become a good benchmark for the two-dimensional quantum simulators.


I. INTRODUCTION
The quantum many-body systems are often represented by Hamiltonian whose composing parts do not commute with each other.This suggests that although the properties of lowenergy states can be easily identified in a certain limit, often referred to as the classical limit, considering the whole Hamiltonian will make the estimation of ground state difficult because one needs to consider the superposition of every state in the full Hilbert space.Moreover, the extra quantum or even thermal effect may lead to further stabilization among some competing low-energy states, causing the energy level crossing or degeneracy lifting, known as the order by disorder (OBD) effect .Besides, the introduction of frustration, through the Hamiltonian [23][24][25][26] or lattice geometry [27][28][29][30][31][32], can also result in the many-fold degeneracy of ground state in the classical picture and such degeneracy enriches the quantum or thermal fluctuation.
Among all, the triangular lattice is the simplest lattice structure for inducing the geometrical frustration.Particularly, the inclusion of anisotropy could lead to the formation of long strings composed of the same spins [33].The shape of strings enjoys a great degree of freedom, leading to a manyfold degeneracy.Classically, such systems and their degeneracies can be realized, composed of self-assembled colloidal * weilintu@keio.jpparticles in a monolayer [34,35].For this soft-matter system, particles move to the opposite walls that confine them to maximize the free volume, corresponding to the up-anddown spin scenario and forming a good platform for studying the dynamics of frustration.In d dimensions, the entropy of degenerate ground states is proportional to L d−1 , with the number of sites N ∼ L d , resulting in a sub-extensive degeneracy.For the above-mentioned frustrated colloidal system, its thermal OBD effect has been studied and the authors unveiled the conclusion that the straight stripe should be favored [36].On the other hand, the quantum counterpart for such OBD effect in zero temperature is also of interest.When an extra non-commutable term, such as the spin exchange coupling, is introduced to the Ising-like Hamiltonian, the quantum fluctuation could further lower the energy and the amount is dependent on the classical configurations, resulting in the lifting of degeneracy.Starting from the identical classical system, usually the quantum and thermal fluctuations result in similar outcomes and select the same ground state for weakly frustrated magnets.However, this is not always guaranteed and many counter examples have been previously provided [13,16,21,37].Therefore, it is worthwhile to examine the quantum version of this OBD effect.However, this effect has not yet been studied to the best of our knowledge.One possible reason is because such system can hardly be prepared in real materials.However, the recent success in simulating various quantum systems with neutral Rydberg atoms [38,39] has attracted huge attention, leading to the so-called quantum simulator [40][41][42][43].It has been shown that the quantum Ising model with both longitudinal and transverse fields can be simulated [44] by utilizing the Rydberg blockade [45,46].Other spin models, such as the quantum XY model [47][48][49], XXZ model [50,51], or even the XYZ model [52], can also be generated through the microwave engineering.In the following content we will provide a possible setup to study the quantum OBD effect in our target system.We also conduct both analytical and numerical calculations and conclude that the stripe pattern is ultimately selected by the OBD effect.Our conclusion based on the theoretical analysis could become a good benchmark for the real-world device of a two-dimensional quantum simulator using Rydberg atoms in the future.
Moreover, thanks to the advance of cold-atom technique [53][54][55][56], many physical scenarios can now be simulated with emphasized quantum effect, leading to another type of quantum simulators.Especially, the dipolar quantum gases made of erbium and dysprosium recently attract wide attention for the realization of supersolid out of its Bose-Einstein condensate [57][58][59][60][61].Note that, by placing dipolar molecules in optical lattices, researchers have realized suitable quantum simulators for various Hamiltonian with electric [62][63][64][65][66] or magnetic dipoles [67].While Rydberg atoms are often used to create the electric dipole-dipole interaction to synthesize the spin exchange interaction, for the molecules carrying magnetic dipoles it can be viewed as a bosonic system probed by the Bose-Hubbard model [68][69][70].This makes the related studies for various extended Bose-Hubbard models necessary and many works have previously been done [71][72][73][74][75][76][77][78][79][80][81][82][83][84][85].In our previous work, we have studied the effect by tilting the polar angle of dipole and realized a competing scenario between different lattice structures [77].Within the original square lattice, we discovered that a three sub-lattice supersolid phase ought to show up by tilting the polar angle.Related studies have also demonstrated the potential of such techniques in fabricating artificial physical scenarios of interest.Besides the quantum OBD effect by Rydberg-atom simulator, therefore, in this work we also study the hard-core bosons with dipoledipole interaction in the triangular lattice with different tilted dipolar angles, mimicking the distortion of lattice.
In the following content, we will first study the quantum OBD effect using the real space perturbation theory (RSPT) [16,[86][87][88][89][90] and infinite projected entangled-pair state (iPEPS) [91][92][93].We first map the two-level Rydberg system to a spin-1/2 scenario and show that when the number of up spins (n ↑ ) is equal to that of down spins (n ↓ ), a subextensive degeneracy is present in the classical limit.Introducing the spin exchange coupling, energy correction can be considered with RSPT, and we focus on the leading-order correction between the stripe and kinked configurations.Our calculation shows that the energy correction favors the stripe pattern, coinciding with previous results for thermal OBD [36].As a further confirmation we apply the iPEPS calculation for a simplified toy model, showing that the stripe state is indeed more stable.Next, for the dipolar bosons in triangular optical lattice with tilted polar angle we conduct numerical simulations with the cluster mean-field theory (CMFT), which has been shown effective in capturing distinct phases in our pre-vious work [77].We unveil fruitful phase diagrams with various underlying phases, including the supersolid phase.We then conclude our work and argue some potential scopes for further studies.

A. Rydberg-atom quantum simulator
The adoption of Rydberg atoms for a quantum simulator has recently become a topic of interest.It has been shown that the quantum Ising model can be synthesized in a large scale with many interesting properties [94][95][96], as well as the potential application for quantum information [40].However, the adoption of Rydberg atoms for quantum simulation or computation is largely hindered by the limited lifetime (a few 100 µs) for low-angular-momentum (ℓ) Rydberg states [97].In contrast, when exciting the Rydberg states to a very large principal quantum number (n) and utilizing the levels with largest ℓ and magnetic quantum number (m), the lifetime of these Rydberg states can be greatly enhanced, named after the circular Rydberg states [98].Their naturally long lifetime originates from the fact that these large-ℓ states can only decay through the emission of a low-frequency photon, which is polarized parallel to the quantization axis.This feature makes the lifetime of circular Rydberg states much longer than that of regular low-ℓ Rydberg states.Since it is a single-mode decay, the lifetime can be further elongated by placing the atoms in an environment which prohibits the field of resonance of the atomic transition [99,100].Adopting a capacitor parallel to the plane where atoms reside, the lifetime again gets largely enhanced and an useful lifetime of at least 1.1 s for an atomic chain has previously been reported [101].Such measure can even be helpful for a non-zero-temperature environment and a millisecond-lived circular Rydberg state has recently been reported in room temperature [102].This makes the quantum simulation and quantum computation much more feasible [101,103,104].
From now on we denote the circular Rydberg states as |nC⟩ and such states fulfill the following condition: ℓ = |m| = n − 1.These states have a wavefunction distribution of a torus whose circular orbit has a radius r n = a 0 n 3/2 , where a 0 is the Bohr radius.The circular states can be stabilized by an external electric field (F ), which also defines the quantization axis, perpendicular to the plane of circular orbits.Through a proper measure a two-level system composed of two circular states, |n 1 C⟩ and |n 2 C⟩ where n 2 > n 1 and C means the circular state, can be generated.
The interaction between any two Rydberg atoms comes from the electric dipoles carried by the atoms themselves, leading to a manifest dipole-dipole interaction [105,106].When two atoms are far from each other the interaction hardly affects the two-atom Rydberg state and thus a second-order perturbation is adopted, giving rise to a van der Waals interaction proportional to 1/r 6 , where r represents the distance between sites [107,108].On the other hand, when the energy defect of two dipole-dipole coupled states is tuned to zero the strongest Rydberg-Rydberg interaction, proportional to 1/r 3 , appears through the Förster resonance process [109,110].According to Ref. 101 one needs to choose n 2 = n 1 ± 2 to have a flexible simulator, where both the spin-exchange and van der Waals interaction are proportional to 1/r 6 .Since the spin-exchange effect might lead to the decay of circular states to nearby degenerate states, a magnetic field (B) parallel to the electric field is cast so that the highly degenerate hydrogenic manifold of the same principal quantum number n can be lifted and the circular state is thus isolated [101].The magnetic field also prohibits the hybridization of spin up and down channel in each circular state.As a result, we obtain clean two-level Rydberg states for the purpose of quantum simulation.
The Hamiltonian for a pair of Rydberg atoms with the following basis, [|n 1 C, n 1 C⟩, |n 1 C, n 2 C⟩, |n 2 C, n 1 C⟩, |n 2 C, n 2 C⟩], can then be written as [101] where in the Pauli basis it reads with While J is independent on the magnetic and electric fields, J z and B z are flexible varying the field strength and when the intersite spacing is 5µm, h 48−48 = 2.2 GHZµm 6 , h 48−50 = 2.66 GHZµm 6 , and h 50−50 = 3.03 GHZµm 6 when F =9 V/cm and B=13 Gauss, while v 48−50 = −0.539GHZµm 6 stays independent of F and B [101].
Placing the circular Rydberg atoms in the optical lattice [111,112], a highly tunable simulator for the effective XXZ model with longitudinal field can thus be formed.In order to engineer the triangular lattice to simulate the quantum OBD effect, a better way is to trap the individual atom through the optical tweezers [113].One recent work has demonstrated such as possibility by trapping the rubidium circular Rydberg atoms into a 3 × 6 array of tweezers [114].It is reported that the trapping time can scale to over a few milliseconds while the simulation time only demands a few microseconds.A very recent paper has reported the existence of long-lived circular Rydberg states by trapping the alkaline-earth atoms through the optical tweezers [115].Their discovery also widens the potential usage of this kind of quantum simulator.Therefore, by properly manipulating the positions of optical tweezers the deformed triangular optical lattice can be generated, as shown in Fig. 1.There, the diagonal bonds (black bonds) is longer in its intersite spacing so that the sub-extensive degeneracy in the classical limit (J = 0) can be simulated.We demonstrate three different configurations in Fig. 1 for the (a) stripe, (b) kinked, and (c) mixed configurations.Mixture of stripe and kinked configurations in a larger lattice gives rise to various possible states, forming the many-fold degeneracy.Note that the diagonal spacing needs to be large enough to prohibit the off-diagonal magnetic orders.
While the basic setup has been introduced, next we will adopt the theoretical and numerical tools to study this quantum OBD effect, which can benchmark the real-world device of two-dimensional quantum simulator with Rydberg atoms.

B. Real space perturbation theory
We first introduce the RSPT and derive the formula in this section.In this work we study the effect of spin-exchange quantum fluctuation to the system of Rydberg atoms, which can be represented as ⟨ij⟩ ( Ŝ+ As a result, the target Hamiltonian represented by spin-1/2 operators can be written as ) where J z ij ≫ |J|, |B| and B represents the external Zeeman field.We assign J to be isotropic since later we will find out that only the perturbative processes concerning the diagonal spin exchange would matter.To effectively reflect the anisotropy caused by the deformed lattice, we assign that J z ij is equal to J z d when the bond is in the diagonal direction (black bonds in Fig. 1) and J z h when the bond is in the horizontal direction (green bonds in Fig. 1).When J z d < J z h and n ↑ = n ↓ , the ground state configuration enjoys a many-fold degeneracy because the same spins tend to align along the diagonal direction.However, due to the fact that there are two diagonal bonds, we can have many different degenerate configurations and the total number is equal to 2 Ly−1 , where L y is the number of rows along ŷ direction.As a result, the ground state entropy S GS ∼ L y − 1, leading to a sub-extensive degeneracy.Among these degenerate configurations there are two special patterns, the stripe and kinked states.The stripe state is composed of parallel straight lines of up and down spins (Fig. 1(a)), while for the kinked state there are three sublattices occupied by up or down spins within each vertical diamond of the triangular lattice (Fig. 1(b)).The rest degenerate states, called the mixed states, are simply the combinations of stripe and kinked states in a specific proportion (one example shown in Fig. 1(c)).
For applying RSPT, one needs to rotate the local coordinate on each site so that every spin points along the ẑ-direction.Thus, we re-write the corresponding spin-1/2 Hamiltonian in the following form with We can extend this to a more generalized Hamiltonian with longer-range hopping or even off-diagonal terms if needed, simply by modifying the form of coupling matrices.In Eq. ( 5), Note that now the S ± and S z are no longer operators.They simply represent the corresponding operation for the classical spin.Our choice of the "unit" vectors guarantees the transformation to be unitary, although its length is not necessarily the unity.This also explains why we have a pre-factor of 2 for J in Eq. ( 6).Next, we expand the equation and the result becomes After introducing the quantum fluctuation, the spin configuration does no longer align along the local êz i direction and thus we can replace S z i by S − δS i .Then, Eq. ( 7) becomes where representing the classical energy.And H unperturbed is Since Eq. ( 10) does not change the spin configuration, it will not contribute to the perturbation.At last, the perturbative Hamiltonian H p is written as where With these transverse perturbative terms, we can evaluate the energy correction with H p in n-th order with ) where E 0 is the energy for the classical configuration, |0⟩, and E ψi is the energy for intermediate state denoted by |ψ i ⟩.The summation will run over every possible intermediate processes for a given order.For any order of perturbation, its energy correction is obtained under several rules: (1) Each perturbation is evaluated within a linked cluster to ensure the extensiveness of energy correction.On each bond of such cluster H p can be acted for more than one time and the total number of execution of H p determines the order of perturbation.
(2) Any series must start and end with J (1) ij and J (3) ij , since we have rotated the local coordinate so that for every site, ij and B i only consider one spin flip and thus there must be an even number for such process in a series.Because of this constraint, one can always find a lower order series without these two terms and therefore we can ignore them in search of leading order correction.
Because of the above rules, we know that effective leadingorder perturbative series must start and end with J (2) ij .In Eq. ( 13), the denominator considers the energy difference between classical and intermediate configurations.It is given by To adopt RSPT, first we rotate the local coordinates so that S (0) i = Sê z i .Our classical configurations are composed of two different sub-lattices representing the up and down spins and then we choose where a and b denote two different sub-lattices.Without the loss of generality, we choose the planar spin orienting along x axis.In fact, in the classical limit we can simply assign ϕ a = 0 for up spin and ϕ b = π for down spin.Since we hope that our derivation can also be useful in the more general scenarios breaking the U (1) symmetry, we take the offdiagonal order into our consideration.It is also clear to see that by adopting the same coordinate of sub-lattice a or b for every site in Eq. ( 7), we can resume the original Hamiltonian, Eq. ( 4), where the translation invariance is accompanied.
To estimate the leading-order correction, we need to first look for the minimal closed cluster which differentiates the stripe and kinked patterns, and it is the vertical diamonds whose longer diagonal line is parallel to ŷ in Fig. 1.In Fig. 2(a) we demonstrate the minimal different cluster for kinked and stripe states.Then, we need to flip the spins through the intermediate process on the diamond bonds to evaluate the energy correction.Since there are four bonds on each diamond, the leading-order correction which we are looking for is of the fourth order.One of the possible intermediate processes is shown in Fig. 2(b).
The detailed equations and derivation of RSPT are shown in Appendix A and here we only demonstrate the final results.While the degeneracy stays intact when J = 0, turning on the spin-flip coupling will immediately introduce an energy correction to different patterns.For the stripe and kinked patterns, we have ϕ a = 0 and ϕ b = π for the classical configurations.Under this assumption our equations can be largely simplified.More importantly, J d(2) 3 becomes zero and thus only a few terms remain in Eqs.(A4) and (A5).The energy correction terms then become Since J z h is repulsive and thus Γ 5 = −8S 2 J z h < 0 (see Appendix A), the stripe configuration is favored and would be selected.This conclusion coincides with the previous one for the frustrated colloidal soft matter system, whose thermal OBD has been studied through the entropy estimation [36].There, the authors have shown that the straight stripe configuration possesses lower free energy, despite the fact that during the cooling process the energy barrier among different configurations will eventually lead to a mixed pattern.The effective Hamiltonian of such system can be written down as an antiferromagnetic Ising model with deformation in the triangular lattice, similar to our system considered here.At last, we emphasize that since the leading-order correction is of the fourth order, the sign of J does not affect the final conclusion.In Section III, we will consider an effective model with the ferromagnetic J.

C. iPEPS calculation
Since RSPT applies well only for small J, we would like to see if the conclusion still holds when the quantum fluctuation is strong.For this purpose, we adopt the twodimensional (2D) tensor network ansatz, the infinite projected entangled-pair state (iPEPS), to numerically calculate the simplified Eq. ( 4) where we consider the extreme case and ignore both the Ising coupling along the diagonal direction and the longitudinal field.By properly choosing the pre-designated unit-cell size (in this work 2 × 2 or 4 × 2 as shown in Fig. 1) and optimizing the d × D × D × D × D tensors, where d represents the dimension of local Hilbert space and D is the virtual bond dimension, this tensor network ansatz serves as a good variational wavefunction for the quantum many-body systems.We provide some details of iPEPS in Appendix B.
To numerically estimate the ground state ansatz of Eq. ( 18), we need to adopt numerous trials starting from different initial setups.It is because we already recognize that there are many competing states (stripe, kinked, or mixed states) so that we need to avoid the simulation being trapped in some undesirable local minima.As a result, besides the random initial tensors we also start our calculation by constructing the distinct product states, corresponding to the stripe and kinked configurations separately.Since the product state is of D = 1, we then enlarge our tensors to the assigned D, with extra tensor elements being small random numbers.Through this measure, we guide the optimization toward obtaining the stripe or kinked state.After the convergence is reached, the averaged magnetization is defined as where C means the repeating unit cell and N C is the number of lattice sites within the cell.For the stripe pattern since there are two possible choices for the iPEPS unit cell and they correspond to the wave numbers (0, π) and (π, π) after Fourier transformation, we thus define the two related orders for the stripe pattern.The first one is the diagonal magnetic moment where k 1 = (π, π) and k 2 = (0, π).r i represents the coordinate (r x i , r y i ) for each site within the unit cell.For the off-diagonal order, we define In Table I we provide the lowest energies as well as the orders for several different J.We can clearly see that for each J/J z we have m ≈ Sz , meaning that despite the different initial trials the states of lowest energy always belong to the stripe state.Although the off-diagonal orders are almost zero, ⟨ Ŝz ⟩ is highly supressed due to the quantum effect, suggesting that the initial setup for RSPT might not be well applicable.However, our simulation indicates that the conclusion does not change and the stripe state keeps being the ground state, selected by the quantum OBD effect.It is also worth mentioning that for D = 3 and J > 0.3 no matter what initial states we adopt the ground state ansatz always evolves toward the stripe state.We then conclude that our system, where the circular Rydberg atoms are placed in a deformed triangular optical lattice, favors the existence of stripe state.

III. DIPOLAR HARD-CORE BOSONS
In the previous section we have revealed the quantum OBD effect for the targeted Rydberg system and identified that the stripe state is more stable.Next, we aim to unveil the possible underlying phases for the artificial Hamiltonian of dipolar bosons in the triangular optical lattice, which can be generated with cold-atom simulators.The realization of placing magnetic cold atoms in optical lattices has been previously reported for erbium [67] and ytterbium [116], where a three-dimensional (3D) optical lattice is generated with two horizontal lasers in the x-y plane and one vertical laser reflecting the direction of gravity.By confining the atoms in a 3D lattice it would help elongate the lifetime by preventing the inelastic collision of atoms [117,118].Since the vertical tunneling can be reduced by adopting a laser with a longer wavelength, and thus a quasi-2D system can be synthesized.Ref. [67] provides a good example using the laser with wavelength (λ x , λ y , λ z ) = (532, 532, 1064) in the unit of nanometer, creating the effective 3D optical lattice with size (l x , l y , l z ) = (266, 266, 532) and a trapping potential where V i is the lattice depth and k i is the wavevector in different direction.After preparing the Bose-Einstein condensate (BEC) of the target atoms through the optical dipole trap (ODT), the BEC cloud is then adiabatically loaded to the optical lattice within a time scale of millisecond.In the end we obtain the long-lived system whose lifetime scales up to the unit of second [67,117].For our system, however, it demands a triangular optical lattice for its realization and Ref. 119   unveiled this possibility using Rb BEC.There, they adopt a three-laser trapping within the x-y plane whose wavevectors mutually share a 120 • enclosing angle.Optionally, the recent proposal of using optical tweezers to trap the cold atoms might also provide a more flexible platform in designing the quantum simulators [120,121].
Next we head to discuss the Hamiltonian of interest and its realization through cold-atom simulators.The magnetic quantum gases is composed of atoms which can be seen as magnetic dipoles even at zero magnetic field.The origin of dipole moment comes from the spin and orbital angular momentum of electrons, as well as some minor contribution from the nuclear spin.This fact results in the high susceptibility of atoms to an external Zeeman field, which we could use to control the overall orientation, expressed by the polar and azimuthal angles shown in Fig. 3, of the magnetic dipoles [68].This simplifies the general form of dipole-dipole interaction in Ref. 122 and makes it the following form where V is the interactive strength and r ij = |⃗ r i − ⃗ r j |. α ij is the included angle by the dipole moment and ⃗ r ij .This interaction contributes to the on-site repulsive interaction, inter-site interaction, and the density induced tunneling.In this work we consider the hardcore limit meaning that only one atom is allowed to occupy one site in every snapshot.This condition can be naturally fulfilled once the on-site interaction is much larger than the inter-site interaction, which can be tuned through controlling the lattice spacing [123].Once the doubly occupancy is prohibited, the density induced tunneling process also disappears.Along with the normal tunneling process, we end up obtaining the following extended Bose-Hubbard Hamiltonian with the dipolar interaction where b † i and bi represent the creation and annihilation operators of the hard-core boson, with the number operator being ni = b † i bi .V nn denotes the nn dipole-dipole interaction among bosons and we neglect its long-range tail since it decays rapidly.We emphasize that the dipole-dipole interaction in this section is different from that for the Rydberg atoms.As explained in Section II A, for the circular Rydberg atoms one makes use of the electric dipole-dipole interaction of neutral atoms to generate the spin exchange coupling.On the other hand, in this section the cold atom itself can be seen as a magnetic dipole and thus mutual interaction with the form indicated in Eq. ( 22) is cast among atoms.Moreover, the atoms act as quanta and thus tunneling effect takes place for minimizing the energy, leading to the hopping term in Eq. ( 23) that lowers the energy.
The interaction in Eq. ( 22) can be easily tilted with an external magnetic field.In this work we study the case when the dipole moments are tilted as ϕ = 0 (see Fig. 3).Thus, for the nn interaction, the two interactive terms along diagonal direction are equal and become attractive when polar angle θ is large enough.Therefore, there are two different terms for where indices h and d indicate the interaction along the horizontal and diagonal directions, as shown in Fig. 1 labeled by green and black bonds, separately.As a result, Eq. ( 23) is akin to Eq. ( 4) with negative J since there is a one-toone mapping by simply replacing b † i → Ŝ+ i , bi → Ŝ− i , and ni → Ŝz i + 1 2 [26].Therefore, by numerically studying Eq. ( 23) with different tilting angles the resulting phase diagrams cover a wide range of XXZ -like model with anisotropy in the triangular lattice.
In this work, we adopt CMFT [73,74,77] for the construction of phase diagram.Some Details of applying CMFT are indicated in Appendix C. When the strength of hopping term is small, the ground state is an ordered state with U (1) symmetry preserved.As we increase t, the solid order will dissolve and phase transits to the superfluid state through a first-order phase transition.However, in some small windows between solid and superfluid phases, two orders, namely the structural and superfluid orders, can co-exist and a supersolid phase manifests.To identify distinct phases we first introduce the order parameters adopted here.For solid and supersolid phases, they both possess structural order which is defined as where C again represents the repeating unit cell whose size is equal to 3 × 3 or 4 × 4 in our CMFT calculation, and thus N C = 9 or 16.As for the superfluid order, we calculate its condensate density defined as The above two order parameters help identify different phases in our numerical phase diagrams.
We then conduct the CMFT calculation and plot the phase diagrams for several polar angles in Fig. 4. When tilting is small (Fig. 4(a)), the diagonal solid (blue) and supersolid (purple) dominate the phase diagram before entering the superfluid phase.Its structural factor has the modulating momentum (2π/3, 2π/3), and is named after the diagonal stripe in our earlier work [77].In the supersolid phase the ñ(2π/3, 2π/3) and ρ 0 orders co-exist.This phase appears after a continuous transition from the diagonal solid.Further enlarging t, solid order disappears and the phase transits into superfluid (red) phase discontinuously.
On the other hand, the central orange lobe indicates the classically highly degenerate solid configurations (stripe, kinked, and mixed).Thanks to the diagnosis of RSPT and iPEPS in the previous section, we have learned that the true ground state must be in the stripe configuration, although the numerical energies of stripe and kinked configurations by CMFT are very close to each other due to the small t.The green phase in Fig. 4(b) and 4(c) represents a supersolid phase coming from the stripe solid and thus we call it the stripe supersolid.Note that our results indicate that the transition from stripe solid to stripe supersolid, as well as the one from stripe supersolid to superfluid are both continuous.
One can see from Fig. 4 that as we further tilt the polar angle, the central lobe grows and its corresponding supersolid phase also starts to appear (Fig. 4(b)).Finally, as the polar angle is large enough, dominant ordered phases in the phase diagram will be replaced by the stripe solid and supersolid (Fig. 4(c)).The scenario described here is just opposite to the one in our previous work, where a square lattice effectively shifts to a triangular one by tilting the polar angle [77].Here, we start from the triangular lattice where the corresponding phase possesses three sub-lattices (see, e.g., Fig. 3 in Ref. [77] for explanation), and eventually ends up with a square one with only two sub-lattices, according to the phase diagram.At last, we emphasize that including the longrange tail of dipole-dipole interaction will automatically lift the degeneracy even in the classical limit.In such scenario we expect a even more fruitful phase diagram with many different ground-state configurations but it requires a highly precise experimental tool to sort out.We will leave the further studies for this issue in the future work.

IV. CONCLUSION
In this work we study two different systems of the dipolar quantum simulator.Placing the circular Rydberg atoms in a deformed triangular optical lattice, a sub-extensive degeneracy can be realized and the spin exchange interaction, arising from the electric dipole-dipole interaction of Rydberg atoms, leads to the quantum OBD effect.The degeneracy of lowest-energy configuration in the classical limit is caused by the deformation, making the van der Waals repulsion only be seen in the horizontal direction between nearby sites.Once the spin exchange term is introduced, the overall Hamiltonian is akin to the XXZ model with anisotropic interactive potential.With RSPT and iPEPS, we predict that the stripe configuration is the true ground state, in coincidence with the thermal counterpart of this OBD effect according to the previous work [36].Next, we consider an extended Bose-Hubbard model with magnetic dipole-dipole interaction in the optical lattice, which is related to the cold-atom quantum simulator.We then provide the phase diagrams for the nn dipolar hardcore bosonic Hamiltonian with different tilting angles in the triangular lattice, using CMFT.We exploit the competing scenario proposed in our earlier work [77] and demonstrate that the effective lattice structure shifts from triangular lattice to a square one, as well as various different phases including the supersolid.
The technique of neutral atoms for quantum computing has recently attracted huge attention.By placing the Rydberg atoms with the desirable geometric controls, large-scale quantum Hamiltonian can be simulated with the quantum analog simulators [124].However, the current platforms mainly focus on the Ising model and its variants.In this work, we push one step forward and combine the spin exchange interaction as well as the van der Waal repulsion.We propose a possible experimental setup using the circular Rydberg atoms in optical lattice for realization, which benefits a further investigation by experimental groups.Our analytical and numerical studies, as well as the conclusion, might be useful benchmarks once the simulator can be put into practice, instead of checking the full phase diagram (see, e.g., Section III.B of Ref. 101).Moreover, controlling the anisotropy of the triangular lattice, we might reach a novel regime where the order disappears and a state with larger entanglement entropy manifests [44].This would benefit further studies with quantum simulators as well as numerical calculation through quantum Monte Carlo or tensor networks, and we will leave it for the future consideration.It is also important to note that our proposal unveils a possible protocol in studying the dynamics of glassy phases through quantum simulator.In Ref. [36] the authors revealed that although the straight stripe benefits from lower free energy and should be more favorable in arbitrarily low temperature region, the free-energy barrier between different configurations results in a metastable disorder state while cooling down the temperature.Similarly, if we introduce the spin exchange coupling in the same manner a non-trivial glassy dynamics [125] can likely be simulated.
Moreover, the quantum simulator made of cold atoms with magnetic dipoles plays another important role in simulating the extended Bose-Hubbard model.While each atom can be seen as a magnetic dipole, their collective alignment can be manipulated by an external field.By tilting the polar angle of dipoles, the corresponding phase diagram changes and peculiar phases also appear.This implies that through some simple action such as tilting the polar angle, many extraordinary physical scenario can be artificially realized through such quantum simulators.In sum, we believe that our both proposals indicate alternative paths in applying quantum simulators for studying the intriguing physics and demonstrate their great potential for artificial many-body systems.15) and ( 16)), we can now write down the corresponding J (1) ij and J (2) ij terms.For diagonal bonds we have where we assign J For the horizontal bonds we simply need to replace J z d in Eq. ( A1) and (A2) with J z h , but in fact it does not affect the leading-order correction since our cluster of interest is the vertical diamond (Fig. 2(a)).Moreover, all possible intermediate configurations as well as their corresponding −⟨ψ|H unperturbed |ψ⟩, labeled by different Γ and composing the denominator of energy correction, also need to be considered.For this purpose, we first list all the possible intermediate configurations during the spin flipping process, shown in Fig. 5, with each configuration's corresponding minus unperturbed energy, E 0 − E ψ , labeled by different Γ.We list their formula in Table II For simplicity, we further conclude that Along with Γ k ( 4 ) and Γ ′ k ( 4 ), we have obtained all the ingredients needed for the denominators of the energy corrections.
Finally, through considering all possible intermediate processes akin to the one in Fig. 2(b), we obtain the energy correction.For the kinked state, and for stripe state, δE A5) where δE (4) is the energy correction per site.Eqs.(A4) and (A5) will then be used for studying the degeneracy lifting.

Appendix B: Infinite projected entangled-pair state
To study our effective model when the fluctuation is strong, we have no other option but to turn to the assistance of numerical tool.Due to the frustration that hinders the usage of quantum Monte Carlo, we adopt the infinite projected entangledpair state (iPEPS) [91][92][93]126] for our purpose.There are two parts for iPEPS tensor network.One is the tensors within the repeating unit cell, where each composing unit is a rank-5 tensor with four auxiliary bonds (dimension D) and one physical bond representing the size of local Hilbert space (dimension d = 2 for S = 1/2).Since originally iPEPS was designed for the square lattice, in the triangular lattice the coordinate number for each site is equal to six.Instead of increasing the number of auxiliary bonds which would make the computation cumbersome, however, we act the local Hamiltonian also along the virtual bonds (dashed bonds in Fig. 6(a)) when calculating the energy.As a result, the overall lattice structure is effectively equivalent to the triangular lattice.
Another important part for iPEPS lies on the environment tensors.Since the thermodynamics limit can be extrapolated through the corner transfer matrix renormalization group (CTMRG) procedure [127][128][129], we need to first trace out the physical legs of each tensor with its adjoint one, forming the so-called double layered tensor where each bond is of dimension D 2 .We then adopt the CTMRG and after the fixedpoint tensors have been found, they act as the effective "bath" surrounding the bulk tensors, where we label the ones on the corners with C and ones on the edges with T .One example with the bulk size equal to 4 × 2 is provided in Fig. 6(b).
At last, we need to optimize the tensors so that the overall network can represent the ground state ansatz of the target Hamiltonian.Here we adopt the variational optimization with the usage of automatic differentiation, first introduced by Liao et al. [130].We adopt the widely applicable package, peps-torch [131], for our calculation.For readers aiming for some further information, we refer to our previous works in Refs.132 and 133.

Appendix C: Cluster mean-field theory
In this Appendix we briefly introduce the CMFT method that we use for numerically solving Eq. (23).For adopting this method, we first divide our Hamiltonian into two parts: H C within the chosen cluster, and H ∂C that contains the terms connecting the bulk to the environment on the boundary of the cluster.H C possesses the exact form of Eq. ( 23), and the mean-field decoupling only takes place in H ∂C , written as where the prime symbol indicates that this summation is between site i on the boundary of the cluster and site j connected to i outside the cluster.Our effective Hamiltonian is then written as H eff = H C + H ∂C .Although we write down the most general form in Eq. (C1), in this work we apply CMFT for the case with only the nearest-neighbor coupling.We then exactly diagonalize the effective Hamiltonian and obtain the ground state to calculate the mean-field parameters, ⟨b j ⟩ and ⟨n j ⟩, for the following iteration.After the mean-field parameters converge to consistent values, our calculation reaches its self-consistent solution.In CMFT, exact diagonalizations are performed within the chosen cluster.Therefore, it is not expected to include the long-range entanglement beyond the cluster size.To study this finite size effect and infer the phases in the thermodynamic limit, a common practice is to compute the order parameters at varying cluster sizes and extrapolate it to the infinity.
CMFT is considered to be more accurate than the regular single-site mean-field theory, for that it can capture the correct correlation within the selected cluster.Thus it is suitable for states whose correlation length is small, meaning that deep inside the ordered phase (or away from transition boundary) we can use CMFT and obtain quite accurate results.Because of the above-mentioned reason, for the boundary lines shown in Fig. 4, we have extrapolated the energies and order parameters for both the first-and second-order transition boundaries to give a better outcome.In addition, due to the randomness from the initialization process, CMFT could converge to distinct phases with the same set of parameters.In this case, one should compare the energies among them to determine the correct ground state.Accordingly, it becomes a challenge when two stable phases appear with energy differences less than or comparable to the finite size effect.
In our CMFT calculation for Eq. ( 23), it provides definitive results on most of the phases shown in Fig. 4.However, it is nontrivial to resolve the ground state phase between the stripe and kinked supersolid due to the fact that the energy difference between these two lies in the forth-order energy correction from the perturbation theory.For the 4 × 4 cluster we employed, this difference is comparable to the finite size effect coming from the boundary, as can be estimated by investigating the changes in energy upon displacements of the solid patterns.However, by checking the order parameters we can still diagnose different phases.

FIG. 1 .
FIG. 1. Configurations for (a) stripe, (b) kinked, and (c) mixed states.Red (pink) dots indicate that the Rydberg atom occupies the |n1C⟩ (|n2C⟩) state, and the black (green) bonds represent the bonds for diagonal (horizontal) interaction, written as J z d (J z h ) in Eq. (4) (also see the context).The brown dotted boxes indicate the minimum repeated unit cell for the stripe and kinked patterns, which we will later adopt for the iPEPS calculation.

FIG. 2 . 23 .
FIG. 2. (a) The four-site diamond cluster that distinguishes kinked and stripe states.Red and pink dots indicate sublattice a and b separately.For kinked state it has two configurations with three sublattice a or b while for stripe state the numbers of a and b sublattice are equal.(b) An example for the fourth-order tunneling process within the diamond cluster.The perturbative process corresponds to J (1) 12 → J (2) 14 → J (2) 34 → J (1) 23 .Black dots indicate the places where spins have been flipped. has

FIG. 3 .
FIG. 3. Schematic demonstration of dipolar interaction in the triangular lattice.Dots denote the lattice site and arrows represent the dipole polarization.Such polarization can be parametrized with polar (θ) and azimuthal (ϕ) angles.The lower right panel shows the projection of lattice from above for a better demonstration of the azimuthal angle.
L.T. acknowledge the useful discussion with Mike Zhitomirsky during the International Conference on Strongly Correlated Electron Systems 2023 (SCES 2023), Pedro Lopes from QuEra Computing Inc., Thomas Ayral during the SQAI-NCTS Workshop on Tensor Network and Quantum Embedding in Tokyo, Hyun-Yong Lee, Tsuyoshi Okubo, S. R. Ghazanfari, and Xinliang Lyu.Authors also acknowledge the Advanced Study Group (ASG) Tensor Network Approaches to Many-Body Systems organized by the Center for Theoretical Physics of Complex Systems (PCS) of the Institute for Basic Science (IBS) in Daejeon, Korea.H.-K.W. is supported by JQI-NSF-PFC (supported by NSF grant PHY-1607611).N.K. is supported by JSPS KAKENHI Grants No. JP19H01809 and No. JP23H01092.T.S. is supported by JSPS KAKENHI Grant.No. 21K03390.W.-L.T. is supported by the Center of Innovations for Sustainable Quantum AI (JST Grant Number JPMJPF2221).This work is also mainly supported by the above-mentioned JST Grant.Appendix A: Details for the real space perturbation theory In this Appendix, we work on the details for RSPT.With the choice for the local coordinate shown in the main text (Eqs.(

FIG. 5 .
FIG. 5. Possible configurations during the perturbation process for (a)kinked and (b)stripe states.For the ground state there are two configurations in kinked state and thus their evolving must be considered separately.Different configurations that give equal denominator contribution are labeled with the same Γ.

FIG. 6 .
FIG. 6.(a) Composition of the rank-5 tensors into the formation of triangular lattice.Grey solid bonds represent the auxiliary bonds with dimension D, and blue open legs are the local indices with dimension d = 2 for the S = 1/2 system.Dashed bonds are not the real legs for the tensors but we act the Hij = Ŝ+ i Ŝ− j + Ŝ− i Ŝ+ j on these bonds too for simulating the triangular lattice.(b) The double layered tensors with bond dimension D 2 in the bulk, enclosed by the blue dashed box, and the environment tensors.The dimension of bonds interconnecting C and T tensors is labeled by χ, with χ ≥ D 2 to guarantee the accuracy.