Phases of d-orbital bosons in optical lattices

We explore the properties of bosonic atoms loaded into the d bands of an isotropic square optical lattice. Following the recent experimental success reported in [Y. Zhai et al., Phys. Rev. A 87, 063638 (2013)], in which populating d bands with a 99% fidelity was demonstrated, we present a theoretical study of the possible phases that can appear in this system. Using the Gutzwiller ansatz for the three d band orbitals we map the boundaries of the Mott insulating phases. For not too large occupation, two of the orbitals are predominantly occupied, while the third, of a slightly higher energy, remains almost unpopulated. In this regime, in the superfluid phase we find the formation of a vortex lattice, where the vortices come in vortex/anti-vortex pairs with two pairs locked to every site. Due to the orientation of the vortices time-reversal symmetry is spontaneously broken. This state also breaks a discrete Z2-symmetry. We further derive an effective spin-1/2 model that describe the relevant physics of the lowest Mott-phase with unit filling. We argue that the corresponding two dimensional phase diagram should be rich with several different phases. We also explain how to generate antisymmetric spin interactions that can give rise to novel effects like spin canting.


I. INTRODUCTION
Interest in systems of cold atoms in optical lattices has greatly increased during the last decade [1] partly because of their versatility as simulators of quantum systems. More precisely, the flexibility, control, and cleanness of these systems have led to numerous realisations of paradigm many-body models of condensed matter physics. The field was greatly stimulated with experimental explorations of the Mott insulator-superfluid transition of the Bose-Hubbard (BH) model [2] following the theoretical proposal of Ref. [3]. Today the list of achievements is long, but to mention just a few: singlesite resolved detection [4], simulation of magnetism and spin models [5], realising synthetic magnetic fields [6] and topological states [7], demonstration of a fermion Mott insulator [8], as well as various dynamical studies like equilibration and Lieb-Robinson spread of correlations [9]. Common to all the above examples is that the essential physics appear on the lowest energy band, the s band. In particular, for spinless particles on the s band, the onsite atomic states are not degenerate. We know, however, that many interesting phenomena in condensed matter theories have their origin in the degeneracies of onsite states/orbitals. For example, understanding 'orbital physics' [10] is essential for giving a full description of superfluid properties of 3 He [11] and the metalinsulating transitions in metal oxides [12]. Spurred by this, in recent years the first steps towards controlled studies of orbital physics with cold atoms have been taken [13].
Degenerate orbital atomic states appear naturally on excited bands of optical lattices, and it was predicted that bosonic atoms loaded into the first excited, the p bands, of a two dimensional (2D) square optical lattice give rise to interesting physics beyond that found on the lowest band [14]. In particular, the superfluid order parameter at zero temperature is complex valued and time-reversal symmetry is spontaneously broken [14][15][16]. Also the insulating phases provide rich physics with the possibilities to study exotic quantum magnetic phases [17,18]. In different lattice geometries, the system may display 'stripped phases' [19], and in 3D the state of bosons occupying an isotropic cubic lattice become frustrated both in the insulating as well as in the superfluid phase [20,21]. Fermionic systems in the p band are also known to feature rich physics, with a plethora of different phases [22].
Experimentally, populating the p bands with bosonic atoms was first realised by accelerating the lattice potential in such a way that the system traversed through a Landau-Zener transition non-adiabatically [23]. However, a more thorough experimental analysis of the coherence and life-time of p-band bosons was performed by Müller et al. [24] using a two-laser Raman loading technique from an insulating s-band state. In particular, it was found that the atoms on the p band relax very rapidly (on the order of a few tunneling times), and that the life-time is relatively long in comparison to the tunneling time. This enables coherence to be established from the initial Mott insulator and the possibility to detect dynamical processes. More recently, a 2D superlattice was used in order to prepare bosonic atoms in hybridised states composed of s and p orbitals [25]. In this experiment, and in a follow up one [26], clear evidences of a complex order parameter were presented. The same loading method was further benchmarked by loading atoms arXiv:1501.03514v1 [cond-mat.quant-gas] 14 Jan 2015 into f -orbital states [27], i.e. the third set of excited bands.
A completely different approach to initialise bosonic atoms in excited bands was proposed and demonstrated by Zhai et al. [28]. There the atoms start in a zero momentum eigenstate (held in a broad harmonic trap) and then a sequence of 'on/off pulses' are applied to the atoms, where the 'on pulse' consists of a standing wave which couple the zero momentum states to higher momentum states by photon recoils. With this method the atoms were loaded into the d band with a fidelity as high as 99%, which is much better than what was reported for the p bands (around 80%). This technique is, however, more easily applied for atom transfer between bands of even or odd parities, e.g. s → d or p → f and so on, and is therefore not suitable for preparing p-band atoms starting from the ground state. Surprisingly, the properties of d-band bosons remain unexplored. Considering the recent experiment [28], it is therefore of importance to work out the relevant physics of such systems.
Here we analyse the zero temperature phase diagram of d-band bosons in an isotropic 2D square lattice. On the d bands we have three orbitals which are generally denoted d x 2 , d y 2 , and d xy . However, due to the anharmonicity of the single site potential well, degeneracy is broken and the energy of the d xy orbital is somewhat higher than the other two orbitals which stay degenerate. The Mott insulator boundaries are determined within the Gutzwiller method. At low atomic densities, the d xy -orbitals are only weakly populated and we can therefor reduce the number of orbitals to two, i.e. we can describe the physics with an effective two-orbital BH model. The properties of the superfluid phase are analysed analytically by minimising the onsite energies. Like for the p-band superfluid phase, the order parameter is complex but due to the particular shape of the d orbitals a new type of vortex lattice is formed with two vortex/anti-vortex pairs fixed to every site. At higher densities, the general structure of the vortex lattice persists but now also the d xy -orbitals becomes populated.
One difference with respect to p-band BH models is the lack of a characteristic Z 2 -parity symmetry for bosonic atoms in the d band. The lack of this symmetry stems from density-assisted orbital-changing collision terms. Nevertheless, there is another Z 2 -symmetry that is spontaneously broken in the superfluid phase. In the Mott insulating phase, we employ perturbation theory to derive an effective spin model, where it is seen that the new interaction terms appear in the form of an external field. In particular, the first insulating phase (the Mott with unit filling) can be described by an XYZ Heisenberg model in an external field. The model is, of course, still supporting the Z 2 -symmetry and we expect a rich phase diagram within the insulating phase with the possibility of ferromagnetic phases with broken symmetry. It should be noted that the phase diagram of the 2D XYZ model in an external field is not fully known. We also present a brief discussion on how the so called Dzyaloshinsky-Moriya (DM) interactions (or antisymmetric spin exchange) may appear when the lattice is no longer isotropic but the onsite degeneracy is kept. While not studied in detail, we point out that the insulating phases with higher fillings are also of great interest. Such phases will represent exotic magnetic states of higher spins or quasi-spins of SU (N ) models.

II. d-BAND BOSONS
Throughout we will assume that all atoms reside on the d bands, i.e. the loading is perfect and that any timescales are short compared to the decay times, driven by intra-atomic collisions, of these meta-stable atoms. The second condition implies that relaxation occurs on a time shorter than a few µs which is the experimentally measured life-time of d-band atoms [28]. Relaxation typically happens during few tunneling times [24], which allows coherence on these excited bands to form before substantial dissipation/decoherence sets in. This has in particular been experimentally demonstrated for bosons loaded in the p [24,25] and the f [27] bands.
Before discussing the physics on the d bands, let us first recapitulate the main results for the system of bosons in the p bands of 2D square lattices. On the isotropic lattice, p-band bosons support two degenerate orbitals and a very interesting phase diagram. As for the regular sband BH model [3,29], the interplay between repulsive interaction and kinetic hopping leads to Mott insulating phases with integer fillings n 0 = 0, 1, 2, ... and condensed superfluid phases [14,16]. The orbital degree of freedom implies, however, that both the insulating and the superfluid phases are non-trivial. In the condensate, the superfluid order parameter ψ(x, y) is a complex function representing a condensate with a vortex pinned to every site. The properties of the tunneling of p-band atoms in 2D further implies that vortices in neighboring sites have angular momentum quantised in opposite directions (i.e. winding numbers are ±1) [14][15][16]. Thus, there is a 'checkerboard' lattice of vortices and the two possible configurations underline the discrete Z 2 -parity symmetry of the model, as well as the spontaneous breaking of time-reversal symmetry. The orbital degree of freedom also makes the insulating phases very rich. The n 0 = 1 Mott phase with one atom per site can be mapped onto an effective Heisenberg XYZ model [18], whose phase diagram in 1D is qualitatively known [30] but not known in 2D. In the spin language of the Heisenberg XYZ model, the breaking of the Z 2 -parity symmetry represents for example an Ising-like transition.
The goal of the present paper is to demonstrate that also the experimentally relevant d-band bosonic model hosts very interesting physics with some notable differences from the physics in the p band. A first guess would be that on the d bands, one again encounters a checkerboard of quantised vortices but with winding numbers ±2 instead. This would be the direct generalisation of the re-sults appearing on the p bands. To construct such states all three orbitals need to be populated, but as we will see in the following only two of them are actually populated due to broken degeneracy. As a result we find a different type of vortex lattice. However, before discussing the possible phases on the d band, we give the corresponding many-body Hamiltonian in the next Section, and outline how we will analyse it within the Gutzwiller mean-field approach.
A. Single particle properties Figure 1: The three d-bands; E x 2 (qx, qy), E y 2 (qx, qy), and Exy(qx, qy) for the potential amplitude V = 20Er. The bandwidth gives an estimate of the tunneling strength and the curvature the sign of the tunneling coefficient The single particle Hamiltonian for an atom with mass m in an isotropic 2D square optical lattice is [1] ∂ 2 x + ∂ 2 y +Ṽ sin 2 (kx) + sin 2 (ky) . (1) For convenience we introduce dimensionless parameters by letting the recoil energy E r = 2 k 2 /2m set the energy scale and the inverse wave number k −1 the characteristic length scale. Thus, in dimensionless variables, The spectrum E ν (q x , q y ), consisting of bands separated by forbidden gaps, is characterised by three quantum numbers, a discrete band index ν and two continuous quasi momenta q x and q y . In scaled dimensions, the first Brillouin zone is defined by q x , q y ∈ [−1, +1). Increasing the potential amplitude V implies that the band gap increases while the widths of the energy bands shrink. The 'flatness' of the bands determines the mobility of the atom. Consequently, as V increases the mobility is reduced and the bands become flatter. In the limit of very deep lattice each site mimics a two-dimensional harmonic oscillator and hence the isotropic 2D case is characterised by a single s band, two degenerate p bands, three degenerate d bands and so on for the higher excited bands. The similarity with an isotropic 2D harmonic oscillator can also be seen from the quantum numbers of the onsite orbitals; if |n x , n y represents the eigenstate of the oscillator with quantum numbers n x and n y then the ground-state is the 'vacuum' |0, 0 (s-orbital state), the first excited states are |0, 1 and |1, 0 (p-orbital states), and the second excited states are |0, 2 , |2, 0 , and |1, 1 (d-orbital states). For finite, but large V , the separations between the different set of bands (s, p, d, ...) may remain much larger than the band widths (at least for the lower bands). However, due to the anharmonicity of the potential wells the bands are no longer equidistant from each other, and even the three d bands may split. As an example, in Fig. 1 we display the three d bands (E x 2 ,E y 2 , and E xy ) for a potential of amplitude V = 20. As we can see, the d xy band has a higher energy than the other two. The eigenstates of the single particle HamiltonianĤ sp are the Bloch states φ νq (x) with ν = (n x , n y ) the band index, q = (q x , q y ) the quasimomentum and x = (x, y). Taking the modified (i.e. restrict the integration to the first Brillouin zone) Fourier transform of these, one obtains the Wannier functions w νi (x) where i = (i x , i y ) with i x , i y ∈ N labelling the site. Contrary to the delocalised Bloch functions, the Wannier functions are localised around the site i and can therefore be normalised in the usual way. In addition, although they are not the eigenfunctions of the single-particle problem, Wannier functions are orthogonal and therefore provide an alternative basis for describing such systems. In the infinitely deep lattice limit, they become simply the harmonic oscillator x|n x , n y eigenfunctions localised at each site. For large but finite V , the harmonic oscillator eigenfunctions are not exact, but they still preserve the general structure of w νi (x). That is, for the s band, where ν = (0, 0), w si (x) is approximately Gaussian. For the p bands, of ν = (1, 0) and (0, 1), the Wannier functions w pxi (x) and w pyi (x) are approximately Gaussian in one direction and have a node in the transverse direction (the subscripts x and y tell the direction of the node). In the d bands, ν = (2, 0), (1, 1), and (0, 2), the three Wannier functions are shown in Fig. 2. These are the three atomic d-orbital states, and the notation used is borrowed from atomic/chemical physics: The d x 2 -orbital state has two nodes in the x direction and no node in the y direction. The opposite is true for the d y 2 orbital, while the d xy or-bital has a single node in both directions. Furthermore, since the square lattice is separable, Wannier functions in 2D can be constructed as a product of 1D Wannier functions. Explicit expressions for the d-band Wannier functions are: where we have introduced the subscripts s, p, d to label the Wannier functions of the different bands in 1D.
Naturally, the shape of the Wannier functions play a crucial role in the dynamics of the system: On a singleparticle level, since the d α orbital (α = x 2 , y 2 ) is much wider in the direction of the two nodes, tunneling processes parallel to this direction happen with considerably larger amplitude than in the perpendicular one. Thus, a d x 2 -orbital atom is more mobile in the x than in the y direction. The d xy orbital, however, is equally mobile in both directions. This direction-dependent propagation results in an anisotropy in the problem, that is also present in the p-band system [14][15][16] (although absent on the s-band). This is also reflected as an anisotropy in the single-particle spectrum. As illustrated in Fig. 1, E x 2 has a much larger curvature (inverse effective mass) in the x direction than in the y direction. Moreover, the sign of the curvature is also of importance since it determines whether the effective mass is positive (particle-like) or negative (hole-like), and as we will see below, it also sets the sign of the tunneling coefficient in the extended BH model. Furthermore, on a many-body level the amplitude of scattering processes is determined by overlaps of products of Wannier functions, and therefore their shape strongly influences the possible interaction processes as well as their strengths.

B. Many-body Hamiltonian
The Hamiltonian of the full many-body problem describing contact interactions between the neutral atoms is given bŷ where the coupling constant U 0 = 4π 2 a/m, with m the atomic mass and a the s-wave scattering length. The field operatorsΨ(x) andΨ † (x) annihilate and create a particle at the position x and obey the bosonic commutation In order to constrain the atoms to the d bands (i.e., impose the single-band approximation), we expand the field operators in the corresponding Wannier functions, with the sum running over all the lattice sites, and the bosonic creation/annihilation operators obey d βi ,d † βj = δ ij , where β = {x 2 , y 2 , xy} and all remaining commutators vanish. Following the usual prescription, Ψ(x) and its Hermitian conjugateΨ † (x) are inserted into Eq. (3). We then impose the tight-binding approximation such that the tunneling is restricted to nearestneighbour sites and interactions to only dominant onsite terms. This yields the final form of the many-body Hamiltonian describing bosons in the d band. For later convenience, we separate it into two different parts, one containing all the processes that involve the d xy orbital, and one part that only include contributions of the d x 2 and d y 2 orbitals:Ĥ Furthermore, the two separate parts can be split up intô The free part of the Hamiltonian accounts for the onsite energies E β of different bands together with the chemical potential µ and is given by with the particle number operatorsn βi =d † βid βi . The kinetic energy contributions arê with i,j σ the sum over nearest neighbours in the direction σ = {x, y}, i,j the sum over nearest neighbours in all directions and α = {x 2 , y 2 } the label of the orbital state. The tunneling coefficients satisfy where t α and t α ⊥ are the tunneling coefficients in the directions parallel and perpendicular to the two nodes of the orbital state. Furthermore t α , t α ⊥ > 0 and t p < 0. Explicit expression of the overlap integrals are given in Appendix A.
The interacting part of the Hamiltonian (6) is characterised by various processes. It contains the densitydensity interactionŝ (10) the orbital changing interactionŝ and finally, the density assisted orbital changing collisionŝ Other interaction terms vanish due to the symmetries of the orbitals. Like for the tunnelling coefficients, the various coupling constants can again be found in the Appendix A. As pointed out above, the system of bosons in the p band supports a Z 2 -symmetry which can be spontaneously broken to give a vortex-checkerboard condensed phase or a spin-flop phase in the superfluid or insulating phase respectively. This symmetry arises in that case because p-orbital bosons in different orbital states scatter in pairs and therefore the number of particles in each orbital state is preserved modulo 2. On the d band, this conserved quantity is broken by the density-assisted interactions (12). There exists, however, another non-trivial discrete Z 2 -symmetry in the many-body Hamiltonian (5). With the notation i = (i x , i y ), we let i = (i y , i x ), and the Hamiltonian is invariant under the transformation This swaps both the d x 2 and d y 2 orbitals as well as the lattice indices i x and i y . Note that the interchange of indices is necessary due to the tunneling anisotropy of the d x 2 -and d y 2 -orbital atoms (t d σ = t s σ ). In the next section we will see that this symmetry can indeed be spontaneously broken to give interesting phases. We should also mention that in addition there are two more Z 2 -symmetries. Since the Hamiltonian only consists of quadratic or quartic terms, it is invariant under a signchange of all operators. This parity symmetry can be further decided into two discrete symmetries, one which change the signs of the d x 2 and d y 2 orbitals, leaving the d xy orbital unaltered, and one with the reverse operation. These last two symmetries are trivial in the case when the d xy orbital is unpopulated as in Subsecs. III B 1 and III C.

C. Mean-field approaches
In the simplest mean-field analysis, the ground-state |Ψ 0 of the many-body system can be expressed as the direct product [32,33] with respect to the (complex) amplitudes ψ dαi and ψ dxyi and the particle conservation constrain where N tot is the total number of particles. Note, in particular, that since the HamiltonianĤ BH is normally ordered [34], the energy expectation values are simply obtained by replacing operatorsd αi ,d xyi ,d † αi , andd † xyi by the corresponding coherent state amplitudes ψ dαi , ψ dxyi , ψ * dαi , and ψ * dxyi respectively. This crude mean-field method is capable of giving a qualitative correct picture deep in the superfluid phase with a large or moderate number of particles per site, as will be utilised in Subsection III B. However, once the interaction is increased relative to the tunneling, the system becomes insulating and the above approach fails to predict such a phase: the method calls for improvement.
As a next step to advance on this coherent-state ansatz, (16) we consider intra-site particle fluctuations, but still without intersite correlations where every single site state is with |n x 2 , n y 2 , n xy i a Fock state with n x 2 d x 2 -orbital atoms in site i and so on. Note that normalisation implies nx,ny,nxy c (i) nxnynxy 2 = 1 for every site i. The corresponding semi-classical energy functional is similar to the above, E Gutz c (i) n x 2 n y 2 nxy = Gutz Ψ|Ĥ BH |Ψ Gutz , and the amplitudes c (i) n x 2 n y 2 nxy are obtained again via self-consistently minimisation of the energy functional. This procedure, called the Gutzwiller mean-field method, is particular for its ability of identifying the existence of Mott insulating phases [3,35]. However, although the factorisation (17) becomes a more accurate approximation in higher dimensions, it already provides a good estimate of the insulating-superfluid phase transition in 2D [36].

III. PHASE DIAGRAM
In this section we characterise different phases that can be expected for bosonic atoms in the d bands. To quantitatively determine the phase boundaries of the insulating phases we use the Gutzwiller method discussed in the previous section. For the p-band system [14], the properties of the superfluid phase can be understood from analysing the single site problem on a mean-field level and then considering how the presence non-zero tunneling terms establish global long range order. Applying the same recipe for the d-band system fails for low atomic densities; it suggests that all three orbitals are populated while from the Gutzwiller analysis it follows that almost all population is distributed in the d x 2 and d y 2 orbitals. For this reason, when identifying the superfluid state in cases where the number of atoms per site n 0 5 with the coherent state mean-field approach, it is reasonable to constrain the analysis to account for only two of the orbitals. At higher densities however, the analysis must include all three orbitals.
Since mean-field methods are less reliable for probing the physics in Mott insulator phases, this regime will be studied here in terms of an effective spin model derived to describe the limit where the magnitude of any tunneling coefficient |t| |U |, where U is the typical interaction strength.

A. Superfluid-insulator boundaries
The boundaries separating the Mott insulators from the condensed superfluid phase are determined, as mentioned above, by the Gutzwiller ansatz wave function (17). The condensate order parameters for the three orbitals are given by ψ d β = Gutz Ψ|d d β |Ψ Gutz . Within the Gutzwiller method, the insulating phase is described by vanishing order parameters, while it is non-zero in the condensed superfluid phase.
Before discussing the insulating-superfluid transition we make a brief detour on other possible phases absent in our case due to the imposed approximations. When the lattice amplitude is not too large, the larger width of the higher band Wannier functions may lead to nonnegligible contributions from interaction between neighbouring sites. The effective model then includes nearest neighbour interaction terms, which if large enough (compared to the onsite interaction terms), may lead to supersolid phases [37]. The existence of such a phase have been demonstrated for bosonic atoms in the p band [31] and occurs when the strength of interactions on neighbouring sites is as large as one fourth of the onsite interaction strength. For 30 ≤ V ≤ 70 as considered in this work, however, the strength of interactions on neighbouring sites is always less than 1% and such terms can be safely neglected. Nevertheless, it could be interesting to analyse the possibility of using the d rather than the p bands for realising such novel phases. As one would be required to relax the tight-binding, and most likely also the single-band approximation [38], this study goes beyond the scope of the present work.
Let us now turn to the Gutzwiller predictions for the transitions. The self-consistent minimisation of the Gutzwiller energy functional E Gutz c (i) n x 2 n y 2 nxy is based on the use of the Nelder-Mead simplex algorithm [39]. Since we are interested in the phases for the system with different densities, the chemical potential µ (independent of the orbital flavour) is added to the Hamiltonian, see Eq. (7). For the numerics we truncate the number of Fock states in the sum (18) to four for every orbital, i.e. up to three atoms for each orbital state and nine atoms in total per site. This gives an onsite Hilbert space dimension of 64, and we can accurately capture the first four Mott lobes. Going to higher insulating states would require Fock states of larger particle number, and thereby a rapid slow down of the numerics. Furthermore, we scale Eq. (6) with the interaction strength U .
The results of the numerical Gutzwiller calculation are presented in Fig. 3, showing the superfluid order parameter in the µt-plane. In this range of the system's parameters, the occupation of the d xy orbital n xy ≈ 0 everywhere. Throughout, proper Wannier functions are used in the computation of the overlap integrals (A1) presented in the Appendix A, and we consider the potential amplitudes 30 ≤ V ≤ 70. These large potential amplitudes (up to 70 recoil energies) are required to ensure The corresponding order parameter for the dxy-orbital is approximately zero except for t/U < 10 −3 where it is, however, at least one order of magnitude smaller than ψ). As is seen, the chemical potential is varied such that the first four Mott insulators are illustrated (dark blue regions representing a vanishing superfluid order parameter). For this plot, the relative strengths between the interaction terms has been taken as {Uxy, Up, Upx, Upxy, Uxxy}/U = {0.17, 0.9, 0.3, 0.04, −0.03} corresponding to a lattice with amplitude V = 40. With this lattice depth, E x 2 = E y 2 = 0 and Exy = 1.6U , and the relative tunnelings {t α ⊥ , |t p |}/t α = {0.002, 0.007}. the validity of the single-band and tight-binding approximations in the d band, that accommodates very broad Wannier functions. The order parameters for the d x 2 and d y 2 orbital are identical and for simplicity we call it ψ = | d x 2 i | = | d y 2 i |. Note that the effective Gutzwiller tunneling t used in the figure is twice the sum of the sand d-band tunnelings. More precisely, out of the four nearest-neighbours in the 2D system, tunneling processes in the d band occur with (1D) d-band tunneling strength to two of the neighbouring sites, and with (1D) s-band tunneling strength on the remaining two neighbouring sites. The general structure of the Mott lobes are similar to those found for the regular one orbital BH model using mean-field methods. That is, the extent of the higher lobes fall off as roughly n −1 0 for the filling factor n 0 [36]. For higher densities n 0 , the typical interaction energy becomes considerably larger than the energy gap separating the E xy from the E x 2 and E y 2 bands, and therefore the occupation of the d xy orbital is expected to increase. This regime can be explored numerically but requires a higher truncation in the number of particles, which makes computations much more costly. To circumvent this issue, we studied the system at fixed chemical potential µ/U and for various values of t/U and indeed, the results show that with the the same parameters as in

B. Onsite superfluid states
The Gutzwiller mean-field method of the previous Subsection provided us with an estimate of the Mott lobes, but the orbital dependence within the superfluid phase was left implicit. We will approach this problem next using of the coherent state mean-field ansatz for clarity.

Moderate atomic densities
Generalising the results from the p bands [14,16] we would expect a doubly quantised vortex on every site for the d-band system. Such a state takes the form where N s is the number of atoms at the site and is fixed a priori. However, according to the previous Subsection, the occurrence of a state as that of Eq. (19) is precluded for lower atomic densities by the negligible occupation of the d xy orbital. Moreover, from only w d x 2 (x) and w d y 2 (x) it is not possible to construct a doubly quantised vortex state, which leads to the conclusion that the superfluid order parameter is a state of a different type.
Following the results of the previous Subsection and focusing first on the low-density case, we approximate ψ dxy = 0, and rewrite the single-site order parameter as where 0 ≤ θ, ≤ π and 0 ≤ ϕ ≤ 2π. Even though this is valid for systems with onsite particle number limited to 10, we believe this gives a good picture of the superfluid state. Indeed, as we show in the next Subsection, a similar state is found when computation includes the third orbital state.
which after minimisation yields the fixed point (θ 0 , ϕ 0 ) θ 0 = π/4, ϕ 0 = arccos (−U xxy /U xy ) , where ∂ θ E cl [θ 0 , ϕ 0 ] = ∂ ϕ E cl [θ 0 , ϕ 0 ] = 0. This solution, is depicted in Fig. 4 (a) and (b), for the parameters corresponding to a lattice with V = 40. The first plot (a) shows the atomic onsite density |Ψ vor (x)| 2 and the second (b) the phase ϕ of the order parameter. As is evident, the onsite condensate hosts two vortex/anti-vortex pairs. Each vortex has a phase winding of ±2π around each singular point. As argued above, while the naïve harmonic approximation on the d bands suggests a single multiply quantised vortex in each site, they do not occur here due to broken degeneracy of the orbitals. It has also been known that doubly quantised vortices should be energetically unstable in a harmonically trapped Bose-Einstein condensate [40]. Throughout, proper Wannier functions are used in the computation of the overlap integrals (A1) presented in the Appendix A, and we consider the potential amplitudes 30 ≤ V ≤ 70. These large potential amplitudes (up to 70 recoil energies) are required to ensure the validity of the single-band and tight-binding approximations in the d band, that accommodates very broad Wannier functions. Typical ratios between the interaction strengths for V = 40 are presented in the caption of Fig. 3.
The condensates at neighbouring sites have the same ordering of the vortex pairs since t α , t α ⊥ > 0, which supports a "0" or "2π" phase locking between consecutive sites. Note that this is always the solution in case of occupation of only the d x 2 and d y 2 orbitals, and that it does not depends on U . In addition, the configurations defined by the vortices orientations spontaneously break the Z 2 -symmetry defined in Eq. (15) as well as time-reversal symmetry. For the solution given in Fig. 4 (a) and (b), and starting from the vortex in the upper left corner, the phase winds clockwise/counterclockwise/clockwise/counter-clockwise.
The alternative configuration with broken Z 2 -symmetry has instead counter-clockwise/clockwise/counterclockwise/clockwise winding. Now the natural question is what are the properties of the phase with restored Z 2 -symmetry? To answer this, we first notice that according to (22), the shape of the two-vortex pairs depend on the ratio R ≡ U xxy /U xy . When R → 1 the vortex/anti-vortex approach each other and annihilates at R = 1. For R > 1, the order parameter reads instead The atomic density and the order parameter phase are shown in Fig. 4 (c) and (d). The density vanishes identically around a circle, where also the phase makes a πjump from ϕ 0 = 0 to ϕ 0 = π. Such a behaviour describes a dark soliton, immobile and confined to the lattice site.
If the ratio R could be externally controlled, then the system could be driven through a phase transition from a 'soliton'-superfluid into a symmetry broken 'vortex'superfluid. In addition, since the derivative of the classical energy is discontinuous at the critical point R c = 1, it suggests that the transition is of the first order. This Here we assumed that the atomic density is rather low such that we can neglect any population of the dxy orbital. In the upper two plots we consider a lattice with an amplitude V = 40. It is seen in (b) that the phase of the order parameter winds 4π at the four points where the density vanish. This reflects the presence of four vortices -two vortex/anti-vortex pairs. In the lower two plots we use the same parameters but put Uxxy = −2Uxy in order to reach the regime where the state qualitatively changes. Here a dark soliton is separating the central peak from the surrounding circle.
is not so surprising, since the states separated by the critical point are topologically different. How to control R experimentally may be non-trivial, but nevertheless, the present analysis sheds some light on the underlying physics of d-band bosons.
We end this Subsection by noticing that on the p band in 3D the orbitals are also three-fold degenerate. In this case an artificial instability appears whenever the harmonic approximation is considered, i.e. by replacing the Wannier functions by the corresponding harmonic oscillator eigenfunctions, which implies that it can become favourable to populate the orbitals unequally [16,41]. On the d band, the depopulation of the d xy orbital does not derive from such a spurious effect, but instead, it appears together with the self-consistent minimisation of the full (lattice) energy functional, that includes both the interacting and tunneling processes simultaneously. This means, furthermore, that the combined (independent) knowledge of onsite properties with additional considerations accounting for effects of the tunneling -as one usually proceeds for studying the system in the p band, is not enough to determine the physics in the d band.  1) where all three orbitals are considerably populated. The difference between the upper and lower plots is the sign of the order parameter ψ dxy (the two states have equal energy), which reflects the states in two neighbouring sites in the lattice. Despite the fact that the dxy orbital is largely populated in this case, the general structure of the superfluid state shown in Fig. 4 (a) and (b) survives, i.e. the onsite order parameter hosts two vortex/anti-vortex pairs. However, the state gets distorted with two of the vortices at the very edge of the distribution. As for Fig. 4, the potential amplitude V = 40.

High atomic densities
In this Subsection we consider a larger particle number n 0 1 such that all three orbitals are populated. Relaxing the assumption ψ dxy = 0 adds a new variable to be accounted for in the minimisation of the energy functional. By using the parametrisation of the three-orbital order parameter in spherical coordinates, the normalization constraint is automatically taken care of and by further choosing the overall phase to be zero we are left with four angle parameters; E cl [θ, φ, ϕ, ϑ]. Since in the general case we do not find analytical solutions for the minimisation (fix point) the system is studied numerically (again utilising the Nelder-Mead simplex algorithm). When occupation of the d xy -orbital increases, it is in principle possible to form doubly excited vortices in each of the lattice sites. However, as argued above, it is questionable whether such state could be energetically favourable. Another aspect arising with all three orbitals populated is that the other two Z 2 -symmetries discussed after Eq. (15) become relevant as will be demonstrated. These two symmetries are reflected in the fact that the energy E cl is invariant under either (ψ d x 2 , ψ d y 2 , ψ dxy ) ↔ (−ψ d x 2 , −ψ d y 2 , ψ dxy ) and/or (ψ d x 2 , ψ d y 2 , ψ dxy ) ↔ (ψ d x 2 , ψ d y 2 , −ψ dxy ) (we note that these two are equivalent via a gauge transformation of the overall phase).
The numerically obtained ground state is displayed in Fig. 5 (a) and (b). Here, the population of the three orbitals d x 2 , d y 2 , and d xy are 0.36, 0.33, and 0.31 respectively, and hence none of the orbital states dominates the others. The effect of the non-zero population of the d xy orbital is to squeeze and rotate (in this example counterclockwise) the atomic distribution without destroying the two vortex/anti-vortex pairs; the larger population in the d xy orbital the more rotated and squeezed is the atomic distribution. Thus, as is clear from the deformations in Fig. 5, vortices of (for example) positive winding are closer to the centre of the site where the density is higher. This suggests that as interactions become stronger and the different orbitals become effectively more degenerate, the squeezing is also increased and finally two of the vortices appear were the density is vanishingly small. In this case the onsite states start to approach doubly quantised vortex states. This is an interesting observation, since as we have pointed out earlier such a vortex is typically not stable. However, one should bare in mind that here we have two singly excited vortices coming infinitely close to each other. Furthermore, this continuous deformation does not imply a violation of conserved angular momentum; the state only appears as a doubly excited state, and in the full lattice model the doubly excited vortices would alternate between neighbouring sites such that the total angular momentum vanishes.
To understand the full superfluid state that extends over the whole lattice we recall that the tunneling amplitude t p on the p band comes with a negative sign. This implies that ψ dxy changes sign between neighbouring sites while ψ d x 2 and ψ d y 2 do not. This is illustrated in Fig. 5 (c) and (d). The reversed sign of ψ dxy has the effect of rotating the distribution in the opposite direction (i.e. here clockwise), but the winding numbers are unaffected. Now recall that the Z 2 -symmetry of Eq. (15) corresponds to reversal of the phase winding. Taking all into account we have that in the symmetry broken phase the onsite rotation direction alternates between neighbouring sites (breaking of the parity symmetry (ψ d x 2 , ψ d y 2 , ψ dxy ) ↔ (ψ d x 2 , ψ d y 2 , −ψ dxy )) and furthermore the phase winding is determined from breaking of the symmetry (15). In addition, we have numerically checked that in the same way as in the previous subsection, the state in the phase with restored symmetry is characterised by a dark soliton.
C. Insulating phases for filling n0 = 1 Deep inside the Mott insulating phases the relevant physics of the many-body system can be studied in terms of an effective (pseudo) spin model derived from the perturbative treatment of the tunneling processes relative to the interaction [42]. By mapping the problem onto a spin system, the effective dimension of the Hilbert space is greatly reduced and the relevant physics becomes more transparent and easily analysed. Moreover, even though the properties of many spin models are still unknown, known results of previous studies can be applied to the present case.
Since the different Mott lobes are characterised by a definite integer filling n 0 , the derivation of the effective spin model can be most naturally handled with the use of operators that project the eigenvalue problem onto orthogonal particle number subspaces of the Hilbert space. As a consequence, the result of the perturbative treatment is not general to any insulating state, and thus we restrict the study here to the first Mott lobe (n 0 = 1) which is most relevant for experimental studies. We proceed by defining theP andQ operators,P 2 =P and Q 2 =Q, that project, respectively, onto the space of states with unit filling, and the states with at least one doubly occupied site. For the lowest Mott insulator we may again impose the assumption that the d xy -orbital state is negligibly populated, and that the only relevant degrees of freedom considered stem from the d x 2 and d y 2 orbitals.
The eigenvalue problem can then be written as (for a more detailed derivation see the Appendix and Refs. [18,21] whereĤ dkin is given by Eq. (8) and the interaction partĤ U =Ĥ dden +Ĥ dc +Ĥ do , cf. Eqs. (10)- (12). Due to the assumptions of the Mott phase, the resolvent 1/ QĤ UQ − E can be expanded such that the effective Hamiltonian contains contributions to second order in t/U . In addition, the tight-binding approximation allows proceeding with the analysis as a two-site problem. Therefore, we define the basis spanning theP subspace by where |α, β =d † αid † βj |0 corresponds to the state with a d α -orbital atom in the site i and a d β -orbital atom in the neighbouring site j, α, β = {x, y}. In the same way, the relevant states in the basis of the subspaceQ of doubly occupied sites is given by with |0, 2α = 2 −1/2d † αjd † αj |0 and |0, αβ =d † αjd † β,j |0 . Due to the possibility of transferring population between the different orbital states viaĤ dc andĤ do , we notice, however, that the projection of the Hamiltonian onto the H Q subspace is non-diagonal in the basis of intermediate states of the perturbation theory. A practical way of dealing with this situation is to notice that the contributions of the different processes can be obtained from computation of (Ĥ Q − E), whereĤ Q =QĤ UQ , with a subsequent inversion. In addition, since the energy E ∼ t 2 /U we have (Ĥ Q − E) −1 ≈Ĥ −1 Q . More explicitly, ordering the basis of H Q according to (28), and with We determine the final form of the effective Hamiltonian by computing the relevant matrix elements of Eq. (26). This yields (see Appendix B for more details of the derivation) By further employing the Schwinger angular momentum representation [43] together with the constraint of unit fillingn i =n x 2 i + n y 2 i = 1, Eq. (32) can be mapped into a spin-1/2 XYZ model with DM interactions [44,45] and an external field (34) Here J = t x σ t y σ U 2 /Λ, with U 2 = 4(U xx U yy − U 2 xy ) and Λ given by Eq. (31), the anisotropy parameter γ = −4U 2 xy /U 2 , However, in the isotropic lattice as considered in this work, the DM interactions vanish due to the Z 2symmetry discussed in Eq. (15), which in the spin lan- This symmetry is broken in anisotropic lattices. In fact, there t x = t y and t x ⊥ = t y ⊥ , and in such case, since δ = 0, the DM terms reappear in the effective spin model.
Several interesting facts should be noted. First, that in the same way as for the system of bosons in the p band [18], d-orbital bosons provide an alternative realisation of the XYZ Heisenberg model, albeit here with the presence of an external field even in the isotropic case. This is a consequence of the fact that the density assisted processes, Eq. (12), break the Z 2 -parity symmetry associated to orbital changing interactions of Eq. (11) which preserves the number of orbital atoms modulo two. This also explains the presence of the S x i term in the Hamiltonian and it readily follows that in the absence of this field the parity symmetry is restored. Second, that in the same way as for the external field, the DM interaction results from the density-assisted processes. Third, that in the limit of vanishing density-assisted interactions, the effective spin Hamiltonian becomes identical to the corresponding one of the p-orbital system, of Ref. [18]. That this indeed should be the case can be understood from the fact that the anisotropy parameter γ is a consequence of the orbital changing processes. However, a difference between the two models is that while the t α t β ⊥ < 0 for the p-orbital system, t α t β ⊥ > 0 in the d band, and therefore the effective spin model favours primarily ferromagnetic alignment at neighbouring sites in all the pseudo-spin components.
To the best of our knowledge, the phase diagram of this model is not fully known in 2D even for the simplest case of the isotropic lattice, where δ = 0. However, the underlying Z 2 -symmetry suggests the possibility of a rich phase diagram. We note that for Γ/J ∆/J (|γ| < 1) the system is characterised by a highly magnetised state in the x-spin component, while in the opposite limit it is ferromagnetic in the z-component of the spin (noting that ∆ > 0). The later is a symmetry broken phase with magnetisation in either positive or negative z-direction. For J ∼ γ, ∆ another possible symmetry broken ferromagnetic phase could appear with the spin in the xy-plane or possibly a gapless floating phase which exists in the 1D XYZ model. In 1D the qualitative features of the XYZ phase diagram are known [30]. This has been studied in terms of p-band bosons [18], but following the same procedure of that study to reduce the dimensionality to an effective 1D model one would inevitably generate the DM interaction terms for the system in the d band. This, however, gives interesting possibilities since this model is known to host an interesting phase diagram with ferromagnetic and Luttinger liquid phases [46].

IV. CONCLUDING REMARKS
Motivated by a recent experiment [28], we have studied the zero temperature properties of bosonic atoms on the d bands of an isotropic square optical lattice. Due to the particular shapes of the onsite d orbitals, the phases characterising the ground-state of this model are different from those analysed in the past for systems on the p band [14][15][16]. This is not surprising if we think of the orbitals in the harmonic approximation where the p orbitals carry angular momentum components l = ±1 and the d orbitals l = 0, ±2. The |l| = 1 angular momentum of the p orbital is reflected in the singly excited vortex at each site in the superfluid phase. The higher angular momentum components of the d orbitals implies that the onsite state can take different forms than a state with single vortices. The direct generalisation of the superfluid phase on the d band would be a checkerboard lattice of a single |l| = 2 vortex located at every site. However, we found instead that two vortex/anti-vortex pairs appear on each site. The common wisdom is that multiple excited vortices are energetically unstable to form several singly excited vortices [40], and one could argue that this explains why we should not find doubly excited vortices in this system. This simplified argument cannot be true here since the total angular momentum should be preserved, which is not the case for a state carrying two vortex/anti-vortex pairs where the total angular momentum vanishes.
The inter-site phase locking of the different intra-site order parameters is determined by the signs of the different tunneling coefficients. For the d band this implied that the corresponding single site vortex states are ordered in the same way in all the sites. This is in contrast with the checkerboard type (or anti-ferromagnetic) arrangement of vortices in the superfluid phase of the system in the p band. This ferromagnetic phase is a symmetry broken phase with respect to the direction of phase windings of the vortices. The spontaneous breaking of the symmetry (15) is also accompanied by breaking of time-reversal symmetry, which otherwise would imply the possibility of switching between the two degenerate vortex states of different orientation. In the symmetry preserved superfluid state, a dark ring soliton is formed on every site. This phase may, however, not occur naturally in experiments since the strength of the density assisted interaction coefficient must become comparable to the coefficient of the orbital changing interaction. This said, it does not necessarily mean that this phase cannot be monitored but then by external driving like suggested in Ref. [18].
Like for the superfluid phase, the orbital structure also makes the physics of the insulating phases very intriguing. We focused on the lowest filling n 0 = 1, where the Mott phase was shown to feature very interesting properties. With a perturbative treatment, this system was mapped onto an XYZ Heisenberg spin model in an external field. The phase diagram of this model is not known in 2D, but by considering limiting cases we argued that the model should posses a very rich phase diagram. We further showed the appearance of DM interactions in the case where the lattice is no longer isotropic. In particular, by fine tuning the lattice is is possible to restore the degeneracy between the d x 2 and d y 2 orbitals and still break the lattice isotropy. In the harmonic approximation this accounts to fulfilling the condition V x k 2 x = V y k 2 y [17] (in unscaled variables, where V µ and k µ are the potential amplitudes and wave numbers in the directions µ = x, y). Thus, it is possible to study effects arising from the DM interactions, of which the most famous effect is perhaps that of spin canting where an anti-ferromagnetic state builds a finite magnetisation or where the magnetisation in a ferromagnetic state gets quenched [45]. The presence of these terms imply the breaking of the Z 2 -symmetry (15) and as a result the transitions should turn into first order or of Berezinskii-Kosterlitz-Thouless type [47].
In 3D, the n 0 = 1 Mott insulator could be effectively described by an SU (3) pseudo-spin model, as was recently shown for the system in the p band [21]. The density assisted orbital chaining collision terms (12) present on the d band would, however, give some additional terms to the effective model in comparison to that emerging on the p band. Now, what could we expect from the physics on other insulating phases? On the n 0 = 2 insulating phase, the projected (two atoms/site) single site Hilbert space is spanned by the three states {|xx , |xy , |yy }. In the Schwinger spin representation (33) with the constraintn i = 2 one would obtain an effective spin-1 XYZ model. The integer spin implies the possibility of a gapless Haldane phase to appear [48]. Going up into the Mott lobes with larger filling, the effective spin models would be characterised by higher spin, but as long as we remain in 2D, the (pseudo) spins in the two-orbital case would be generators of the SU (2) group. If occupation of thed xy orbital becomes non-negligible, then the pseudospins are the generators of the SU (3) group even in the 2D lattice, with a corresponding spin model with DM interactions in all the components.
Since this work presents the first theoretical study of the physics in the d band, we have limited the analysis to the isotropic 2D lattice. Naturally, other lattice geometries and in other dimensions the physics is expected to change as well as if one considers fermionic atoms instead of bosonic ones. Another aspect left out here is the influence of a trapping potential. On the p band it was recently demonstrated that the presence of a harmonic trap is a most simple way to directly detect outcomes of the anisotropic tunneling on the excited bands [49]. Similar effects as those found for the confined p-band system would also appear on the d band, i.e. time-of-flight detection of the freely expanding atomic cloud or singlesite addressing [50] would reveal information about the intrinsic anisotropic tunneling on the d bands.