XY\!Z Quantum Heisenberg Models with p-Orbital bosons

# Xyz Quantum Heisenberg Models with p-Orbital bosons

Fernanda Pinheiro Department of Physics, Stockholm University, Se-106 91 Stockholm, Sweden NORDITA, KTH Royal Institute of Technology and Stockholm University, Se-106 91 Stockholm, Sweden    Georg M. Bruun Department of Physics and Astronomy, University of Aarhus, DK-8000 Aarhus C, Denmark    Jani-Petri Martikainen COMP Center of Excellence, Department of Applied Physics, Aalto University, Fi-00076, Aalto, Finland    Jonas Larson Department of Physics, Stockholm University, Se-106 91 Stockholm, Sweden
July 22, 2019
###### Abstract

We demonstrate how the spin- quantum Heisenberg model can be realized with bosonic atoms loaded in the band of an optical lattice in the Mott regime. The combination of Bose statistics and the symmetry of the -orbital wave functions leads to a non-integrable Heisenberg model with anti-ferromagnetic couplings. Moreover, the sign and relative strength of the couplings characterizing the model are shown to be experimentally tunable. We display the rich phase diagram in the one dimensional case, and discuss finite size effects relevant for trapped systems. Finally, experimental issues related to preparation, manipulation, detection, and imperfections are considered.

###### pacs:
03.75.Lm, 67.85.Hj, 05.30.Rt

Introduction.– Powerful tools developed recently to unravel the physics of many-body quantum systems offer an exciting new platform for understanding quantum magnetism. It is now possible to engineer different systems in the lab that mimic the physics of theoretically challenging spin models maciek (), thereby performing “quantum simulations” qs (). Along these lines, systems of trapped ions and of polar molecules are promising candidates. Trapped ions, for example, have already been employed to simulate both small trap1 () and large trap2 () numbers of spins. In these setups, however, sustaining control over the parameters becomes very difficult as the system size increases. Furthermore, due to trapping potentials realizations are limited to chains with up to 25 spins. It is also very difficult to construct paradigmatic spin models with short range interactions using systems of trapped ions. Similar limitations appear when using polar molecules, where the effective spin interactions polar1 (); polar2 () are obtained from the intrinsic dipole-dipole interactions. Due to the character of the dipolar interaction, these systems give rise to emergent models that are inherently long range and the resulting couplings usually feature spatial anisotropies.

Short range spin models can instead be realized with cold atoms in optical lattices maciek (). A bosonic system in a tilted lattice has recently been used to simulate the phase transition in a 1D Ising model spin1 (). Fermionic atoms were employed to study dynamical properties of quantum magnetism for spin systems spin4 (); spin2 (). This idea, first introduced in Ref. spin3 (), has also been applied to other configurations, and simulation of different types of spin models have been proposed spinmodels (). However, due to the character of the atomic -wave scattering among the different Zeeman levels, such mappings usually yield effective spin models supporting continuous symmetries like the model. But as the main goal of a quantum simulator is to realize systems that cannot be tackled via analytical and/or numerical approaches, it is important to explore alternative scenarios that yield low symmetry spin models with anisotropic couplings and external fields.

In this paper we propose such a scenario by demonstrating that bosonic atoms in the first excited band ( band) of a two-dimensional (2D) optical lattice can realize the spin- quantum Heisenberg model in an external field. Systems of cold atoms in excited bands feature an additional orbital degree of freedom girvin () that gives rise to novel physical properties maciekliu (), which include supersolids ss () and other types of novel phases novphas (), unconventional condensation nonzero (), and frustration frust (). Also a condensate with a complex order parameter was recently observed experimentally TMuller (); complex (). The dynamics of bosons in the band include anisotropic tunneling and orbital changing interactions, where two atoms in one orbital state scatter into two atoms in a different orbital state. This is the key mechanism leading to the anisotropy of the effective spin model obtained here: These processes reduce the continuous symmetry characteristic of the model, which would effectively describe fermions in the band ref23 (), into a set of discrete symmetries characteristic of the model. In addition, due to the anomalous -band dispersions the couplings of the resulting spin model can favor for anti-ferromagnetic order even in the bosonic case.

We also demonstrate how further control of both the strength and sign of the couplings is obtained by external driving. This means that one can realize a whole class of anisotropic models with ferromagnetic and/or anti-ferromagnetic correlations. To illustrate the rich physics that can be explored with this system we discuss the phase diagram of the chain in an external field. This case exhibits ferromagnetic as well as anti-ferromagnetic phases, a magnetized/polarized phase, a spin-flop and a floating phase eran (). Finite size effects relevant for the trapped case are examined via exact diagonalization. This reveals the appearance of a devil’s staircase manifested in the form of spin density waves. Finally, we discuss how to experimentally probe and manipulate the spin degrees of freedom.

-orbital Bose system.– We consider bosonic atoms of mass in a 2D optical lattice of the form . Assuming that all atoms are in the first excited bands, the tight-binding Hamiltonian is

 ^H=−∑ij,αtαij^a†i,α^aj,α+∑i,α[Uαα2^ni,α(^ni,α−1)+Epα^ni,α] +∑i,α≠α′(Uαα′^ni,α^ni,α′+Uαα′2^a†i,α^a†i,α^ai,α′^ai,α′). (1)

Here creates a bosonic particle in the orbital at site , , and the sum is over nearest neighbors . The tunneling matrix elements are given by where is the Wannier function of orbital at site . Note that is anisotropic. For instance, a boson in the -orbital has a much larger tunneling rate in the -direction than in the -direction. The coupling constants are given by , with the onsite interaction strength determined by the scattering length. The last term in (1) is the orbital changing term describing the flipping of a pair of atoms from the state to the state . Note that this term is absent in the case of fermionic atoms.

Effective spin Hamiltonian.– We are interested in the physics of the Mott insulator phase with unit filling in the strongly repulsive limit . Projecting onto the Mott space of singly occupied sites with the operator , the Schrödinger equation becomes with . Here and  auerbach (). Since , we can take .

The space of doubly occupied states of a given site is three-dimensional and spanned by , , and . In this space, it is straightforward to find from (1), and subsequent inversion yields

 ^H−1Q=⎛⎜ ⎜⎝Uyy/U2−Uxy/U20−Uxy/U2Uxx/U20001/2Uxy⎞⎟ ⎟⎠ (2)

with . In particular, the off-diagonal terms in derive from the orbital changing term. Using (2) we can now calculate all possible matrix elements of in the Mott space,

 ^HMott=−∑ij,α⎛⎝2|tαij|2U¯α¯αU2^ni,α^nj,α+|tαij|22Uxy^ni,α^nj,¯α −2txijtyjiUxyU2^a†i,α^ai,¯α^a†j,α^aj,¯α+txijtyji2Uxy^a†i,α^ai,¯α^a†j,¯α^aj,α) (3)

where , and . By further employing the Schwinger angular momentum representation, , and , together with the constraint , we can (ignoring irrelevant constants) map (3) onto a spin-1/2 model in an external field coms ()

 ^HXYZ=∑⟨ij⟩Jij[(1+γ)^Sxi^Sxj+(1−γ)^Syi^Syj] +∑⟨ij⟩Δij^Szi^Szj+h∑i^Szi. (4)

Here, means summing over each nearest neighbor pair only once. The couplings are given by , , and . The magnetic field is , where is the onsite energy of the orbital .

Equation (4) is a main result of this paper. It demonstrates how -orbital bosons in a 2D optical lattice can realize the quantum spin- Heisenberg model. Several interesting facts should be noted. First, due to the symmetry of the -orbitals girvin () and therefore . Furthermore, since we have anti-ferromagnetic instead of the usual ferromagnetic couplings for bosons. Also, we obtain the model when . The presence of can be traced to the orbital changing term in Eq. (1), which reduces the continuous symmetry of and to a set of symmetries. The symmetries reflect the ‘parity’ conservation in the original bosonic picture which classifies the many-body states according to total even or odd number of atoms in the and orbitals. Since the orbital changing term is absent for fermions, the model with anisotropic coupling is a peculiar feature of bosons in the band. We emphasize that the above derivation makes no assumptions regarding the geometry of the 2D lattice - i.e. it can be square, hexagonal etc.

1D phase diagram.– To illustrate the rich physics of the model, we now focus on the case of a 1D lattice where where quantum fluctuations are especially pronounced. Note that by increasing both the lattice amplitude and spacing in the direction keeping , one can exponentially suppress tunneling in the direction to obtain a 1D model, while the and orbitals are still quasi-degenerate liu1D (). In the 1D setting, we will drop the ”direction” subscript on the coupling constants.

For 1D, the importance of the orbital changing term can be further illuminated, by employing the Jordan-Wigner transformation for fermionic operators . The result is the fermionic Hamiltonian

 ^HK/J=∑n[(^c†n^cn+1+^c†n+1^cn)+γ(^c†n^c†n+1+^cn+1^cn)+ ΔJ(^c†n^cn+12)(^c†n+1^cn+1−12)+hJ(^c†n^cn−12)]. (5)

We see that leads to a pairing term that typically opens a gap in the energy spectrum. Incidentally the limit of in Eq. (5) is a realization of the Kitaev chain kitaev ().

The schematic phase diagram is illustrated in Fig. 1 (a). At zero field, the model is integrable Takhtadzhan (). For large positive values of the system is anti-ferromagnetic (AFM) in the direction. Small values of are characterized by Néel ordering in the direction and the system is in the so-called spin-flop phase (SF). The line for large negative values of is characterized by a ferromagnetic phase (FM) in the direction, and for all the cases, the limit of large external field displays a magnetized phase (PP), where the spins align along the orientation of the field in the direction. These three phases also characterize the phase diagram of the model in a longitudinal field mikeska (). However, for non-zero anisotropy , a gapless floating phase (FP) emerges between the SF and the AFM phases which is characterized by power-law decay of the correlations bak (); eran (); com (). The transition from the AFM to the FP is of the commensurate-incommensurate (C-IC) type whereas the transition between the FP and SF phases is of the Berezinsky-Kosterlitz-Thouless (BKT) type. For there is a first order transition at between the two polarized phases. Finally, there is an Ising transition between the PP and the SF phases.

The experimental realization of the Heisenberg model will inevitably involve finite size effects due to the harmonic trapping potential. Within the local density approximation, the trap renormalizes the couplings so that they become spatially dependent fep1 (), but this effect can be negligible if the orbitals are small compared to the length scale of the trap. In the regime of strong repulsion, the main effect of the trap is instead that it gives rise to “wedding cake” structures with Mott regions of integer filling. This effect was observed in the lowest band Bose-Hubbard model maciek (), and predicted theoretically to occur for anti-ferromagnetic systems AndersenBruun (). To examine finite size effects, we have performed exact diagonalization in a chain with 18 spins with open boundary conditions. Figure 1 (b) displays the resulting finite size ’phase diagram’. The colors correspond to different values of the total magnetization of the ground state. While the PP phase and the AMF phase are both clearly visible, the numerical results reveal a step like structure of the magnetization in between the two phases. We attribute these steps in to a devil’s staircase structure of spin-density-waves (SDW). As we see from Fig. 1 (b), it is only possible to give a numerical result for the PP-SF Ising transition. In particular, the C-IC and BKT transitions are overshadowed by the transitions between SDW. In the thermodynamic limit the staircase becomes complete and the changes in become smooth. One then recovers the phase diagram of Fig. 1 (a). These transitions, between different SDW, are more pronounced for moderate systems sizes. For a typical experimental system with 50 sites, for example, we estimate 15 different SDW between the AFM and PP phases.

Measurements and manipulations.– While time-of-flight measurements can reveal some of the phases complex (), single-site addressing techniques single_site () will be much more powerful when extracting correlation functions. To address single orbital states or even perform spin rotations, one may borrow techniques developed for trapped ions leshouches (). Making use of the symmetries of the and orbitals, stimulated Raman transitions can drive both sideband and carrier transitions for the chosen orbitals in the Lamb-Dicke regime. These transitions can be made so short that the system is essentially frozen during the operation. Driving sideband transitions in this way, spin rotations may be implemented. For example, a spin rotation around is achieved by driving the red-sidebands for both orbitals coms (). As a result, the two orbitals are coupled to the orbital in a configuration and in the large detuned case an adiabatic elimination of the band gives an effective coupling between the and orbitals bruce (). This scheme, thus, realizes an effective spin Hamiltonian with the effective Rabi frequencies and the detuning. Alternatively, Stark-shifting one of the orbitals results in a rotation around . Since the spin operators do not commute, any rotation can be realized from these two operations. Performing fluorescence on single orbital states by driving the carrier transition acts as measuring . This combined with the above mentioned rotations makes it possible to measure the spin at any site in any direction coms (); leshouches ().

Tuning of couplings.– For a square optical lattice, we have . Moreover, in the harmonic approximation , from which it follows that and . This gives ferromagnetic couplings for the component of neighboring spins, while the interactions between and between the components have anti-ferromagnetic couplings. We now show how the relative strength and sign of the different couplings can be controlled by squeezing one of the orbital states. Such squeezing can be accomplished by again driving the carrier transition of either of the two orbitals, dispersively with a spatially dependent field coms (). The shape of the drive can be chosen such that the resulting Stark shift is weaker in the center of the sites, resulting in a narrowing of the orbital. To be specific, assume that the ratio of the harmonic length scales of the and orbitals in the direction is tuned. A straightforward calculation using harmonic oscillator functions yields and . The coupling constants now depend on as and . The inset in Fig. 3 displays the three coupling parameters as a function of for . We see that the relative size and even the sign of the couplings can be tuned by varying . In particular, while always has AFM couplings, they can be made both FM or AFM for and . In the main part of Fig. 3, we sketch the different accessible models as a function of and . This clearly demonstrates that one can realize a whole class of spin chains by using this method.

Experimental realization.– In Ref. TMuller (), the experimental realization of -orbital bosons in an effective 1D optical lattice with a life-time of several milliseconds was reported. With an average number of approximately two atoms per site, the atoms could tunnel hundreds of times in the band before decaying. Since the main decay mechanism stems from atom collisions TomaszSpin (); girvin (), an increase of up to a factor of 5 in the lifetime is expected when there is only one atom per site TMuller (). Typical values of the couplings can be estimated from the overlap integrals of neighboring Wannier functions. Considering Rb atoms, m and , and , we obtain and the characteristic tunneling time s. This corresponds to a few dozens of times smaller than the expected lifetimes TMuller (), which should allow for experimental explorations of our results since relaxation typically occurs on a scale less than ten tuneling times GreinerPaper (). In addition, as pointed out in coms (), it is possible to increase the lifetimes even further with the use of external driving.

A major experimental challenge is to achieve a unit filling of the band. This could be achieved by having an excess number of atoms in the band and then adiabatically opening up the trap such that the unit filling is reached. A minority of sites will still be populated, however, by immobile -orbital atoms. Since the interaction energy between - and -orbital atoms is higher than between two -orbital atoms, processes involving -orbital atoms will be suppressed. The presence of atoms in the band corresponds therefore to introducing static disorder in the system coms (). This may affect correlations dis (), but the qualitative physics will remain unchanged for concentrations close to a unit filling. A more detailed study of this interesting effect is beyond the scope of the present work.

As a final remark we note that the spin correlations discussed here will emerge at temperatures  spin3 (). In addition, we estimate the required entropy Kollath () by equating the critical temperature to the gap between the ground and first excited states in the anti-ferromagnetic phase. Using the energy spectrum obtained from exact diagonalization, yields the entropy per particle . Experimentally one has in fact already achieved  GreinerTalk (), which indicates that our results are within experimental reach.

Conclusions.– We showed that the Mott regime of unit filling of bosonic atoms in the first excited bands of a 2D optical lattice realizes the spin- quantum Heisenberg model. We then illustrated the rich physics of this model by examining the phase diagram of the 1D case. Finite size effects relevant to the trapped systems were discussed in detail. We proposed a method to control the strength and relative size of the spin couplings thereby demonstrating how one can realize a whole class of models. We finally discussed experimental issues related to the realization of this model. We end by noting, that recent experiments reported a 99% loading fidelity of bosons into the -band d_band (), which indeed opens possibilities to probe rich physics beyond spin- chains.

Acknowledgments.– We thank Alexander Altland, Alessandro De Martino, Henrik Johannesson, Stephen Powell, Eran Sela, and Tomasz Sowiński for helpful discussions. We acknowledge financial support from the Swedish research council (VR). GMB acknowledges financial support from NORDITA.

## References

• (1) M. Lewenstein et al., Adv. Phys. 56, 243 (2007).
• (2) R. P. Feynman, Int. J. Theor. Phys. 21, 467 (1982).
• (3) A. Friedenauer et al., Nature Phys. 4, 757 (2008); K. Kim et al., Nature 465, 590 (2010); R. Islam et al., Nature Commun. 2, 377 (2011); R. Islam et al., arXiv:1210.0142; P. Richerme et al., arXiv:1303.6983.
• (4) J. W. Britton et al., Nature 484, 489 (2012).
• (5) A. Micheli, G. K. Brennen, and P. Zoller, Nature Phys. 2, 341 (2006).
• (6) B. Yan et al., arXiv:1305.5598.
• (7) J. Simon et al., Nature 472, 307 (2011).
• (8) T. Fukuhara et al., arXiv:1305.6559.
• (9) J. S. Krauser et al., Nature Phys. 8, 813 (2012); D. Greif et al., arXiv:1212.2634
• (10) L.-M. Duan, E. Demler, and M. D. Lukin, Phys. Rev. Lett. 91, 090402 (2003).
• (11) E. Altman, W. Hofstetter, E. Demler, M. Lukin, New J. Phys. 5 113, (2003), J. RadicÄ, A. DiCiolo, K. Sun, and V. Galitski, Phys. Rev. Lett. 109, 085303, (2012).
• (12) A. Isacsson and S. M. Girvin, Phys. Rev. A 72, 053604 (2005).
• (13) M. Lewenstein and W. V. Liu, Nature Phys. 7, 101 (2011).
• (14) V. W. Scarola and S. DasSarma, Phys. Rev. Lett. 95, 033003 (2005).
• (15) C. Xu and M. P. A. Fisher, Phys. Rev. B 75, 104428 (2007); J. Larson, A. Collin, and J.-P. Martikainen, Phys. Rev. A 79, 033603 (2009); A. Collin, J. Larson, and J. -P. Martikainen, Phys. Rev. A 81, 023605 (2010).
• (16) W. V. Liu and C. Wu, Phys. Rev. A 74, 013607 (2006); C. Wu, Mod. Phys. Lett. B 23, 1 (2009).
• (17) Z. Cai, Y. Wang and C. Wu, Phys. Rev. B 86, 060517(R) (2012).
• (18) T. Müller, S. Fölling, A. Widera, and I. Bloch, Phys. Rev. Lett. 99, 200405 (2007).
• (19) G. Wirth et al., Nature phys. 7, 147 (2011); P. Soltan-Panahi et al., Nature Phys. 8, 71 (2012).
• (20) C. Wu, Phys. Rev. Lett. 100, 200406 (2008); E. Zhao and W. V. Liu, Phys. Rev. Lett. 100, 160403 (2008).
• (21) E. Sela, A. Altland, and A. Rosch, Phys. Rev. B 84, 085114 (2011).
• (22) A. Auerbach, Interacting Electrons and Quantum Magnetism, (Springer, New York, 1998); F. Essler et al., The One-Dimensional Hubbard Model, (Cambridge University Press, Cambridge, 2005).
• (23) For more detailed calculations see the suplementary material.
• (24) X. Li, Z. Zhang, and W. V. Liu, Phys. Rev. Lett. 108, 175302 (2012).
• (25) A. Y. Kitaev Usp. Fiz. Nauk. (Suppl.) 171 10, 2001.
• (26) R. J. Baxter, Exactly Solvable Models in Statistical Mechanics (Academic Press, London, 1982)
• (27) H. Mikeska and H.-J. Kolezhuk, in Quantum Magnetism, U. Schollwöck, J. Richter, D. J. J. Farnell, R. F. Bishop (eds.), (Springer Verlag, Berlin, 2004).
• (28) P. Bak, Rep. Prog. Phys. 45 587 (1982).
• (29) In terms of bosonization (E. Miranda, Brazilian J. Phys. 33, 3 (2003)) and renormalization group arguments, the FP is characterized by irrelevant umklapp terms and accordingly described by the Luttinger liquid theory. Upon entering the XY phase these terms are no longer irrelevant and the phase becomes gapped. eran ().
• (30) F. Pinheiro, J.-P. Martikainen, and J. Larson, Phys. Rev. A 85, 033638 (2012).
• (31) B. M. Andersen and G. M. Bruun, Phys. Rev. A 76, 041602 (2007).
• (32) J. Pietraszewicz et. al., arXiv:1303.5232v2 (2013).
• (33) W. S. Bakr et al., Nature 462, 74 (2009); J. F. Sherson et al., Nature 467, 68 (2010).
• (34) Chapter 5 and 6, Quantum Entanglement and Information Processing, Edited by D. Esteve, J.-M. Raimond, and J. Dalibard, (Eslevier, Amsterdam, 2004).
• (35) B. W. Shore, Manipulating Quantum Structures Using Laser Pulses, (Cambridge University Press, Cambridge, 2011).
• (36) M. Greiner et al, Nature, 415, 39-44 (2002).
• (37) C. A. Doty and D. S. Fisher, Phys. Rev. B 45, 2167 (1992).
• (38) L. Pollet et al., New J. Phys., 10 065001 (2008).
• (39) Markus Greiner, private communication.
• (40) Zhai, Yueyang, et. al., arXiv:1306.3313 (2013).

SUPPLEMENTARY MATERIAL

## Appendix A Derivation of the effective spin model

We are interested in the strong coupling regime where the system is deep in the Mott insulator phase with a unit filling of the lattice sites. A natural way of analyzing this limit involves the use of projection operators that divide the Hilbert space of the associated eigenvalue problem in orthogonal subspaces according to site occupations. We define the and operators that project, respectively, into the subspace of states with a unit occupation and into the perpendicular subspace. They decompose the eigenvalue equation , with its associated energy, in the form

 (^Q^Ht^P+^Q^Ht^Q+^Q^HU^P+^Q^HU^Q)|Ψ⟩=E^Q|Ψ⟩(^P^Ht^P+^P^Ht^Q+^P^HU^P+^P^HU^Q)|Ψ⟩=E^P|Ψ⟩, (6)

where is the interaction part of the Hamiltonian. Since , , , and all vanish, it follows that

 ^Q|Ψ⟩=−1^Q^H^Q−E^Q^Ht^P|Ψ⟩. (7)

By further substitution of Eq. (7) in the eigenvalue equation, we are left with the Hamiltonian which describes the one particle Mott phase of -orbital bosons

 ^HMott=−^P^Ht^Q1^Q^HU^Q−E^Q^Ht^P. (8)

So far this result is exact. It explicitly shows the role of the tunneling in the system, namely of coupling the subspace of states where the sites have unitary occupation with the states that have one site doubly occupied. First, a particle tunnels, say, from the site to , where it interacts with another particle according to what is described by . After interaction, one of the particles is brought back to the site , and the final state is again characterized by lattice sites with a unit filling.

Equation (8) is the starting point in the derivation of the effective Hamiltonian describing the Mott phase of -orbital bosons. The procedure is developed here for an effective 1D system with dynamics along the -axis, but generalization to the 2D lattice is straightforward. Realization of the 1D configuration relies on the adjustment of the lattice parameters, that should contain potential wells much deeper in the than in the axis, but in such a way that the quasi degeneracy between the different orbital states is still maintained. This means that , and furthermore, due to the strong coupling regime condition, we also have that , .

Under these assumptions, the operator in Eq. (8) can be expanded to second order in () in analogy to the customary procedure used for the Hubbard model at half filling Korepin (). In the tight-binding regime considered here, it is enough to consider the 2-site problem. The basis spanning the subspace of states with unit filling is

 HP={|x,x⟩|x,y⟩,|y,x⟩,|y,y⟩},

where represents the state with an -orbital atom in site and a -orbital atom in site . The relevant states for the doubly occupied sites is

 HQ={|0,2x⟩,|0,2y⟩,|0,xy⟩},

which span the intermediate states of the projection operation. We notice, however, that due to the possibility of transferring population between the different orbital states, the projection of the Hamiltonian in the subspace is not diagonal in this basis of intermediate states. This is a peculiarity of the present model and derives entirely from the orbital changing collisions. As a consequence, we compute , with by calculating the projected Hamiltonian in the subspace and taking its corresponding inverse. Since , it is justified to ignore and to consider . Explicitly,

 ^HQ=⎛⎜⎝UxxUxy0UxyUyy0002Uxy⎞⎟⎠

giving

 ^H−1Q=⎛⎜ ⎜⎝Uyy/U2−Uxy/U20−Uxy/U2Uxx/U2000Uxy/2⎞⎟ ⎟⎠,

where .

We determine the final form of the effective Hamiltonian by computing the relevant matrix elements of (8). To this end, we consider in detail all the different cases where the resulting action of the operator of Eq. (8) in the states of the subspace yield non vanishing contribution.

From states of the type

 ^a†α,i^aα,j^H−1Q^a†α,j^aα,i|αi,αj⟩=^a†α,i^aα,j^H−1Q√2|0,2αj⟩
 =√2^a†α,i^aα,j(UββU2|0,2αj⟩−UαβU2|0,2βj⟩)=2UββU2|αi,αj⟩

the effective Hamiltonian acquires a term of the form

 −∑⟨i,j⟩∑α2|tαij|2UββU2^nα,i^nα,j.

In these and the following expressions, it is understood that . In the same way, from the states of the type ,

 ^a†α,i^aα,j^H−1Q^a†α,j^aα,i|αi,βj⟩=^a†α,i^aα,j^H−1Q|0,αjβj⟩
 =12Uxy^a†α,i^aα,j|0,αjβj⟩=12Uxy|αi,βj⟩,

corresponding to the operator

 −∑⟨i,j⟩∑α|tαij|22Uxy^nα,i^nβ,j.

From the states of the type and the following process

 ^a†α,i^aα,j^H−1Q^a†β,j^aβ,i|βi,αj⟩=^a†α,i^aα,j^H−1Q|0,αjβj⟩
 =12Uxy^a†α,i^aα,j|0,αjβj⟩=12Uxy|αi,βj⟩,

the Hamiltonian gains a contribution as

 −∑⟨i,j⟩∑αtαjitβij2Uxy^a†α,i^aβ,i^a†β,j^aα,j

Finally, we consider the states of the type ,

 ^a†α,i^aα,j^H−1Q^a†β,j^aβ,i|βi,βj⟩=^a†α,i^aα,j^H−1Q√2|0,2βj⟩
 =√2^a†α,i^aα,j(UααU2|0,2βj⟩−UxyU2|0,2αj⟩)=−2UxyU2|αi,αj⟩,

that contribute to the effective Hamiltonian with a term that changes the orbital states of the atoms in both sites

 ∑⟨i,j⟩∑α,α≠β2tαjitβijUxyU2^a†α,i^aβ,i^a†α,j^aβ,j.

The resulting expression for the effective Hamiltonian corresponds thus to

 ^HMott=−∑⟨i,j⟩∑α(2|tα|2UββU2^nα,i^nα,j+|tα|22Uxy^nα,i^nβ,j−2txtyUxyU2^a†α,i^aβ,i^a†α,j^aβ,j+txty2Uxy^a†α,i^aβ,i^a†β,j^aα,j). (9)

We now use the orbital states to define the Schwinger spin operators

 ^Sz=12(^a†x^ax−^a†y^ay)^S+=^Sx+i^Sy=^a†x^ay^S−=^Sx−i^Sy=^a†y^ax, (10)

and together with the constraint of unit occupation of the lattice sites in the Mott phase, i.e. , we rewrite Eq. (9) as

 ^HMott=−∑⟨i,j⟩(Jzz^Szi^Szj+Jxx^Sxi^Sxj+Jyy^Syi^Syj)−∑iJz^Szi.

Thus, within the strong coupling regime, the physics of the Mott insulator phase is equivalent to the spin-1/2 Heisenberg model in an external field. In terms of the lattice parameters, the expressions for the various couplings follow

 Jxx=2txtyUxy(1−4U2xyU2)
 Jyy=2txtyUxy(1+4U2xyU2)
 Jzz=4|tx|2UyyU2+4|ty|2UxxU2−|tx|2Uxy−|ty|2Uxy
 Jz=4|tx|2UyyU2−4|ty|2UxxU2+(Eosx−Eosy)

In terms of Eq. (3) of the main text, we can identify , , and .

## Appendix B Single site addressing of orbital states

Single site addressing for the present setup implies selective detection/manipulation of the two orbitals. Since the spin is encoded in external spatial degrees of freedom rather than internal atomic electronic states, methods like those described in Refs. ion1 () would not work. To control the spatial state of the atoms at single sites we may instead apply methods borrowed from trapped ion physics leshouches (). Similar methods were already employed in the experiment TMuller () in order to load bosons from the band to the band. Müller et al. of Ref. TMuller () did not, however, consider single site addressing and more importantly they did not discuss control of the orbital degree of freedom.

Two internal atomic electronic states, e.g. an and an state for Rb atoms, are Raman coupled with two lasers. This transition is described by the matrix element where and are the laser amplitudes and wave vectors, respectively, and the detuning of the transitions relative to the ancilla electronic state. The spatial dependence of the lasers will induce couplings between vibrational states of the atom, i.e. different bands. The time duration for a -pulse, for example, can be made very short by making the effective Rabi frequency large. In particular, this time can be made short on any other time scale in the system and one can approximately consider the system dynamics frozen during the applied pulse. Indeed, the same assumption applies to any single site addressing in optical lattices. Furthermore, by driving resonant two-photon transitions we do not need to worry about accidental degeneracies between other undesired states.

Deep in the Mott insulator phase, as considered in this work, we can approximate single sites with two dimensional harmonic oscillators with frequencies . The Lamb-Dicke parameters leshouches (); jonas () become , and within the Lamb-Dicke regime () we can neglect multi-phonon transitions. Thus, in one dimension we have three possible transitions: () Carrier transitions - with no change in the vibrational state, () red sideband transitions - which lower the vibrational state with one quantum, and () blue sideband transitions - which raise the vibrational state with one quantum. The various possibilities are demonstrated in Fig. 3.

Since the different transitions are not degenerate, it is possible to select single transitions by carefully choosing the frequencies of the lasers. Moreover, considering for example , i.e. no component in the direction, it is possible to only address the -orbital. Thus, we have a method to singly address the different orbitals. Full control is achieved when every unitary , where and is an effective rotation angle, can be realized. To start with the simplest example, implementation of , we first note that since we are considering the case with a single atom per site such that it is enough to realize the operation of . This is nothing but a phase shift of one of the orbitals. This is most easily done by driving the carrier transition off-resonantly for one of the two orbitals. Since the driving is largely detuned it only results in a Stark shift of the orbital.

The operation is preferably achieved by simultaneously driving off-resonantly the red sidebands of the two orbitals. The -band will never get populated due to the large detuning while instead the transition between the two orbitals can be made resonant. More precisely, for the three involved states (with the last entry in the ket-vector being the -orbital) the resulting coupling Hamiltonian in the rotating wave approximation has the form a -coupled system bruce ()

 ^HV=⎡⎢⎣00Ω100Ω2Ω1Ω2δ⎤⎥⎦, (11)

where and have been taken real and for now spatially independent. For we adiabatically eliminate the state to obtain the desired Hamiltonian generating the rotation , namely

 ^Hx=[0ΩΩ0]=Ω^Sx. (12)

Note that if the Raman transition between the two orbitals is not resonant, such an action performs a combination of an - and -rotation. To perform -rotations, one could either adjust the phases of the lasers or simply note that . With this at hand, any manipulation of single site spins can be performed. To measure the spin state in a given direction one should combine the rotations with single site resolved fluorescence (i.e. measuring haroche (). More precisely, since the drive laser can couple to the two orbitals individually, one orbital will be transparent to the laser while the other one will show fluorescence. In other words, one measures on a single site. Other directions of the spin are measured in the same way, but after the correct rotation has been implemented to it. Furthermore, with the help of coincident detection it is possible to also extract correlators  blochnew (). Since there is a single atom at every site, the “parity problem” ion1 () of these techniques deriving from photon induced atom-atom collisions is avoided and thereby loss of atoms will not limit our measurement procedure. This summarizes how preparation, manipulations, and detection of single site spins can be performed.

Finally we note that the methods discussed above can be used in a broader context. For example, there is a transition between two -orbital atoms (one - and one -orbital atom) and one - and the -orbital atom joana (). This transition is resonant for any parameters and and could in principle cause rapid decay of the -band state, or even Rabi-type oscillations between the bands. We note, however, that in the experiment of Ref. TMuller () the collisional decay mechanism was surprisingly small despite this resonant transition. Nevertheless, one could suppress this resonant transition to increase the life-time even further with the technique described above: By driving the red sideband for the two -orbital states dispersively, the and bands will be repelled and thereby this breaks the resonance condition for .

## Appendix C External parameter control

The ideas of the previous section can also be utilized to change the system parameters. The simplest example is the application of which implements a shift in the external field . Apart from the external field, it is also desirable to control the coupling in the component of the spin, , and especially to tune it from ferromagnetic into anti-ferromagnetic.

Using the fact that we have

 Δ≈−|tx|2(4UyyUxxUyy−U2xy−1Uxy). (13)

This is most easily estimated in the harmonic approximation. Introducing the widths of the orbital wave functions for the spatial directions , in this limit

 Uxx=Uyy=3Uxy≡u0σxσyσz, (14)

where is an effective interaction strength (proportional to the -wave scattering length). We notice that even though the use of lattice Wannier functions yields a different ratio between and from what is obtained in the harmonic limit tomasz (), it does not affect the qualitative picture of the results discussed here. Using (14) in the expression for we find

 Δ=−|tx|23σxσyσz2u0