Interaction-induced chiral p_x \pm i p_y superfluid order of bosons in an optical lattice

The study of superconductivity with unconventional order is complicated in condensed matter systems by their extensive complexity. Optical lattices with their exceptional precision and control allow one to emulate superfluidity avoiding many of the complications of condensed matter. A promising approach to realize unconventional superfluid order is to employ orbital degrees of freedom in higher Bloch bands. In recent work, indications were found that bosons condensed in the second band of an optical chequerboard lattice might exhibit p_x \pm i p_y order. Here we present experiments, which provide strong evidence for the emergence of p_x \pm i p_y order driven by the interaction in the local p-orbitals. We compare our observations with a multi-band Hubbard model and find excellent quantitative agreement.


INTRODUCTION
Understanding the role of unconventional order for superconductivity is a fundamental task in low temperature physics. Prominent examples in condensed matter are transition metal oxides [1]. Studying the order parameters in these systems is complicated by their vast complexity. A widely debated example is the chiral p x + i p y order possibly formed in strontium ruthenates [2], which has recently attracted much interest because its topological nature may give rise to Majorana fermions [3]. Optical lattices in their lowest band have proven to be a useful experimental arena to emulate superfluidity with exceptional precision and control [4,5]. However, since under most general circumstances bosonic ground state wavefunctions are necessarily positive definite [6,7] and hence topologically trivial, the realization of unconventional superfluid order in optical lattices with bosons is not straight forward. Possible approaches, presently receiving great attention, are based upon amending the scalar light-shift potentials of conventional optical lattices by static abelian and even non-abelian artificial gauge fields [8][9][10][11] or upon use of dynamical lattice potentials [12][13][14]. A method, more closely geared to electronic matter, is to employ atoms in metastable higher bands which provide orbital degrees of freedom [15,16]. Recently, we have shown that upon suitable control of band relaxation bosons can be condensed in the second band of a bipartite square lattice [17]. Momentum spectra were observed consistent with chiral p x ± i p y superfluid order characterized by a spontaneously formed pattern of staggered local angular momenta, which breaks time-reversal symmetry. In the present work, by means of the following line of arguments we present clear evidence that p x ± i p y is in fact formed: We show that for the lowest energy state of a bosonic Hubbard-model accounting for the second, third and fourth bands, the repulsive interaction in the p-orbitals stabilizes p x ± i p y order in analogy to Hund's second rule in multi-electronic atoms. We calculate a characteristic phase diagram with respect to a change of the interaction in the local p-orbitals and an adjustable distortion of the lattice, which tunes the energy minima of the second band. Experimental observations are presented, which show excellent quantitative agreement with the theoretical predictions of this phase diagram. We finally discuss excited state scenarios, which are compatible with the previously observed momentum spectra, but inconsistent with the experimental signatures reported in this work.

LATTICE POTENTIAL
We produce a two-dimensional optical potential comprising deep and shallow wells (A and B in Fig. 1 (a)) arranged as the black and white fields of a chequerboard with an average well depth V 0 and an adjustable relative potential energy offset ∆V [17][18][19]. In the xy-plane the optical potential is given by | η e ikx + x e −ikx + e iβ e iky + y e −iky | 2 . Adjustment of β permits controlled tuning of ∆V ≡ V 0 η(1 + x )(1 + y ) cos(β). A weak harmonic potential along the z-direction provides elongated tubular lattice sites. If η = x = y = 1, the lattice potential possesses perfect C 4 rotation symmetry. In our experiment, we are constrained to fixed parameter values η = 1.03, and x = 0.93 and hence C 4 symmetry is weakly broken. In contrast to Ref. [17], careful power and polarization management of the lattice beams permits controlled adjustment of arbitrary values of y around unity. The second Bloch band, shown in Fig. 1 (b), provides two inequivalent local minima at the edge of the 1st Brillouin zone (denoted by X + and X − ), which are energetically degenerate if the lattice displays C 4 rotation symmetry. By adjusting the lattice distortion parameter y (see Appendix), we can continuously tune their energy difference ∆E ≡ E(X − ) − E(X + ) ∼ 1 − y . In a tight-binding picture, the quantum states corresponding to X ± may be approximated by Bloch states |ψ ± with real-valued Bloch functions ψ ± (r) composed of local p-orbitals (p x , p y ) in the deep wells and local s-orbitals in the shallow wells. Denoting their occupations per unit cell as n p and n s with n 0 ≡ n s + n p , the relative occupations ν p ≡ n p /n 0 and ν s ≡ n s /n 0 only depend on the spatial shape of ψ ± (r) (but not on the local chemical potential) and hence can be tuned via adjustment of ∆V (see Appendix).

MULTI-BAND HUBBARD MODEL
In order to predict the nature of the expected quantum phases for different values of n p and ∆E, we employ a bosonic multi-band Hubbard-HamiltonianĤ (see Appendix) accounting for the three tight-binding bands associated with the three local orbitals p x , p y and s. We include nearest-neighbour (NN) and next-nearest-neighbour (NNN) tunneling processes, on-site interactions for the p-and s-orbitals, and a term accounting for a potential energy difference between p-and s-orbitals, required for modeling the tuning of ∆V . Furthermore, an extra term is included, which introduces a quadrupolar anisotropy of the tunneling between p-orbitals with respect to the (e x + µ e y )-directions (with e x,y shown in Fig. 1 (a) and µ ∈ {−1, 1}). This term with the real amplitude J ,a permits to model the tuning of ∆E according to the relation ∆E ≈ 8ν p J ,a , which is obtained by evaluating the single-particle spectrum ofĤ in momentum space at X ± . Comparing the lowest two tight-binding bands ofĤ with the second and third bands of a full two-dimensional band calculation for the experimentally realized lattice potential lets us determine the tunneling parameters. The largest tunneling amplitude J sp arises for NN-tunneling between s-and p-orbitals (see Appendix).
For ∆E = 0, any linear combination |θ, φ ≡ sin(θ) |ψ + + cos(θ) e iφ |ψ − minimizes the single-particle energy. For repulsive collisions, minimization of the interaction energy in the p-orbitals requires maximal angular momentum [15,16] and hence φ = ±π/2 and θ = π/4 because the local superposition states p x ± ip y are maximally delocalized such that the atoms can optimally avoid each other. For ∆E = 0, depending on the sign of ∆E, either |ψ + or |ψ − minimize the single-particle energy. When ∆E is increased beyond some critical value, the gain of single particle energy exceeds the cost of interaction energy required for eliminating angular momentum, and hence also the total energy is minimized by one of the states |ψ ± . A mean field analysis for the general case ∆E = 0 shows that the ground state ofĤ can be in fact approximated as |θ, ±π/2 = sin(θ) |ψ + ± i cos(θ) |ψ − with the mixing angle θ plotted in the phase diagram in the centre of Fig. 2 versus J ,a and n p U p (U p is the on-site collision energy per particle in the p x -and p y -orbitals). Note that sin 2 (θ) and cos 2 (θ) quantify the relative populations of the condensation points X + and X − , respectively. This phase diagram comprises three regions, separated by second-order phase boundaries highlighted by white dashed lines. In regions (I) and (III) one finds θ = 0 and θ = π/2 respectively (hence, only one of the condensation points X ± in Fig. 1 (b) is occupied), while in region (II) the simple relation holds [20]. The mixing of the two condensates is governed by a competition between the gain of single-particle energy per unit cell n 0 ∆E introduced by the lattice distortion and the interaction energy per unit cell n 2 p U p gained by maximizing angular momentum in the p-orbitals. The pictograms on the upper edge of the phase diagram illustrate the orbital and local phase ordering predicted in the three regions. In regions (I) and (III) the order parameters are real and their local phases indicated by the colors and numbers are arranged in order to maximize tunneling, while interaction energy does not play a role. In region (II) an inherently complex-valued order arises with orbital currents and plaquette currents highlighted by circular and straight arrows, respectively. Interaction energy sets the relative local phases between p x -and p y -orbitals at the same site to be π/2, while the phase relations between orbitals at neighboring sites are arranged to maximize tunneling. The unit cell of the order parameter comprises four unit cells of the lattice potential and time-reversal is equivalent to a shift by one unit cell of the lattice potential. Hence, time-reversal symmetry is broken. Our theoretical considerations are consistent with a numerical analysis based upon the Gross-Pitaevskii equation [21] and a renormalization group analysis [22].

EXPERIMENT
Our experimental procedure begins with a Bose-Einstein condensate of rubidium atoms ( 87 Rb) loaded into the lowest band. By means described in Ref. [17] the atoms are excited to the second band. Their temperature remains sufficiently low, such that after typically 10 ms they condense in the energy minima of the band. In our experiments we wait for 80 ms to ensure complete thermalisation. Below the phase diagram in Fig. 2, we show momentum spectra observed in the different areas (I), (II) and (III). In regions (I) and (III) the observed Bragg-resonances unambiguously prove the realization of the predicted standing-wave order in configuration space sketched in the corresponding pictograms. In region (II), where both points X + and X − are populated, the momentum spectra appear as superpositions of those in regions (I) and (III). This clearly confirms that the underlying quantum state is composed of two states |ψ ± corresponding to condensates at X ± with p x±y order. However, a more precise determination of its nature requires additional information. In the following, we show that the state prepared in our experiment closely follows the phase diagram in Fig. 2 derived for the ground state |θ, ±π/2 ofĤ, which we refer to as scenario A.
In Fig. 3 we compare the observed populations of the condensation points X ± as a function of ν p and the lattice distortion with the predictions of scenario A. The experimental procedure yielding the data points (filled red disks) in Figs.3 (a)-(f) is as follows: For different settings of the lattice distortion parameter y , ∆V is varied adiabatically and the normalized mean occupation difference between both condensation points ν dif ≡ ( n + − n − )/( n + + n − ) is recorded and plotted versus ν p . Here, n ± are the populations observed in X ± in a single measurement and . . . denotes the average over multiple experimental realizations. For the small interaction energies realized in our system -with U p /J sp on the order of a few percent -and temperatures well below the condensation temperature, our mean field zero temperature analysis appears well adapted to model the observations if the finite size of the lattice with spatially varying local values of n p and U p is accounted for. In fact, upon applying a local density approximation, the observations are remarkably well reproduced by the theoretical predictions for scenario A. Using θ according to Eq. 2 with the local values of ∆E, n p and U p , we calculate the corresponding local value of ν dif,th ≡ sin 2 (θ) − cos 2 (θ) at each plaquette in the lattice and subsequently apply an average over all plaquettes to obtain ν dif,th , which is plotted as the solid green lines in Figs.3 (a)-(f). The grey areas reflect the largest uncertainty in this calculation given by the calibration of the lattice distortion parameter. A discussion of the technical details is given in the Appendix. According to the observations in Fig. 3, n + and n − approach each other if the distortion is reduced or if ν p is increased, in excellent quantitative agreement with the predictions for region (II) of the phase diagram in Fig. 2, where the value of J ,a and hence ∆E required to produce imbalanced condensate fractions grows with increasing ν p . This striking behavior is hence identified as a signature of p y ± ip y -order described by scenario A. The physical mechanism may be regarded as an analogue of Hund's second rule for repulsively interacting bosons: maximal local angular momentum is favorable in order to minimize interaction energy.
Aside from the ground state scenario A, the distribution of the Bragg-resonances in region (II) is also compatible with three possible excited-state scenarios: the coexistence of spatially separated real phases |ψ ± (scenario B), a real coherent superposition |θ, (1 ∓ 1)π/2 ≡ sin(θ) |ψ + ± cos(θ) |ψ − with p x+y ± p x−y order (scenario C), and an incoherent mixture described by the density operator ρ θ ≡ sin 2 (θ) |ψ The compelling quantitative agreement of our observations with the predictions for the ground state scenario A is underlined by the fact that the excited state scenarios B, C, and D lead to predictions in explicit contrast to the observations. The phase separation scenario B exhibits spatially separate p x+y and p x−y -orbitals corresponding to the two condensates |ψ ± . Hence, interaction energy does not favor equal populations of the condensation points and thus ν dif,th cannot acquire any dependence upon ν p , in sharp contrast to the observations. The absence of phase separation in our experiment is also supported by the fact that in our Bragg spectra we never observe different coherence areas for the Bragg peaks corresponding to the two condensation points. This would be expected, if the measured condensate populations n + and n − would be attributed to spatially separated |ψ ± -domains of different size. The irrelevance of scenario B in our experiments is not surprising: Coexisting spatially separated domains of real phases |ψ + and |ψ − as in scenario B are energetically more costly than either of the pure phases |ψ ± . A change of the domain sizes merely requires a local redistribution of particles between the p x+y and p x−y -orbitals. Particle transport over many lattice sites is not necessary. Hence, equilibration should occur within a few tunneling times leading to the formation of either |ψ + or |ψ − . The real superposition state of scenario C exhibits spatially separate p x and p y -orbitals. For equal populations of the condensation points every second local s-orbital remains unoccupied due to destructive interference, which is energetically strongly unfavorable. In fact, minimization of the energy of |θ, (1 ∓ 1)π/2 with respect to θ yields either of the values θ = 0 or θ = π/2 only depending on the sign of ∆E and hence ν dif,th can only take the values ±1 (see Sec. X of Appendix). Finally, although the incoherent mixture ρ θ in scenario D provides orthogonal p x±y -orbitals at the same lattice site, the interaction energy gained by equally distributing the atoms among these orbitals is smaller than for scenario A because of the indeterminate phase relation between the two states |ψ ± . We have minimized the energy of ρ θ with respect to θ in order to determine ν dif,th for this scenario. The result (see Appendix, Eq.A3) is plotted as the dashed blue lines in Fig. 3, which obviously disagrees with the observations. In fact, significantly larger values of ν p would be required to equilibrate the condensate fractions as compared to scenario A. Scenario D also appears implausible because the two incoherently superimposed condensates can exchange particles through binary collisions in the shared local s-orbitals of the shallow wells (see pictograms of regions (I) and (III)). This should rapidly degrade the coherence in each condensate, which is not observed. Further insight is gained by analyzing the fluctuations of ν ± ≡ n ± / n + + n − . Histograms recorded for ∆E ≈ 0 show fluctuations of ν + − ν − , well described by Gaussians, centered at ν + = ν − (see Fig. 4 (a)). This observation appears incompatible with the phase separation scenario B, for which, similarly as in Ref. [14], a double peak structure should be expected instead, since equal populations n ± of the two involved condensation points are not preferred. In Fig. 4 sum −∆ν 2 dif )/4 significantly deviates from zero, i.e., the fluctuations of ν ± are strongly correlated. While the comparably small fluctuations of the total condensate fraction ∆ν sum show basically no dependence on ν p , ∆ν dif is sizable and notably decreases as ν p is increased, showing that the interaction in the p-orbitals tend to lock the populations condensed in X + and X − . This observation again clearly rules out scenarios B and C, which do not prefer equal values of n ± as ν p is increased. It clearly reflects the physics of scenario A captured in the phase diagram in Fig. 2. As ν p is decreased, the critical point is approached along the vertical J ,a = 0 -line where all three phases adjoin. The larger ν p the higher the interaction energy to be paid for deviations of ν + − ν − from zero yielding suppression of fluctuations ∆ν dif .

CONCLUSION
Optical lattices with p x ± i p y order open exciting perspectives for future research. Matter wave interference techniques [23] could be used to further study the mutual coherence of the condensates at X + and X − . Imaging of the atoms with single-site resolution as demonstrated in Refs. [24,25] might allow one to directly observe the local angular momentum of the wave function and hence to explore the spontaneous symmetry breaking process. Proceeding to deeper potential wells, one may access the strongly correlated regime where a rich phase diagram of Mott insulators with distinct orbital ordering is expected [26][27][28]. One may also explore topologically protected features in higher bands [29] and, if fermions are used, simulate forms of topological matter [30] resembling those discussed in the context of electronic systems. APPENDIX A I. Lattice potential. Using an interferometric lattice set-up [18], we produce a two-dimensional optical potential comprising deep and shallow wells (A and B in Fig. 1 (a) of the main text) arranged as the black and white fields of a chequerboard with an average well depth V 0 and an adjustable relative potential energy offset [17]. In the xy-plane the optical potential is given by Adjustment of β with a precision exceeding π/300 permits controlled tuning of ∆V ≡ V 0 η(1 + x )(1 + y ) cos(β). A weak harmonic potential (with 40 Hz vibrational frequency) along the z-direction provides elongated tubular lattice sites. If η = x = y = 1, the lattice potential possesses perfect C 4 rotation symmetry. In our experiment, due to unavoidable imperfections of the lattice set-up, we are constrained to fixed parameter values η = 1.03, and x = 0.93 and hence C 4 symmetry is weakly broken. In contrast to Ref. [17], the optical set-up permits controlled adjustment of arbitrary values of y within an interval including y = 1. This is accomplished as follows: the optical standing wave along the y-axis is obtained by a retro-reflected laser beam. The linear polarization of the incoming beam can be rotated with a retardation plate. After retro-reflection the polarization is rotated to precisely match with the z-direction, which exclusively contributes to the lattice potential.
A II. Band structure: tight-binding picture and full-band calculation. One may understand the single-particle bands within a simplified tight-binding (TB) picture. As sketched in Fig. 5 (a), each local vibrational orbital of the two types of wells gives rise to a Bloch band. We operate in the regime where the 2nd, 3rd and 4th bands are significantly closer to each other than their separation from the 1st band, which correspond to the local s-ground states in the deep wells. We may write ∆V = V sp + ∆V sp with a potential energy offset V sp depending on the adjusted value of V 0 , and ∆V sp denoting the potential energy difference between the local s-states in the shallow wells and the local p x and p y orbitals in the deep wells (cf. Fig. 5 (a)). For V 0 = 6 E rec used in our experiment, V sp = 5.7 E rec . If ∆V sp < 0, the second band arises primarily from the local s-states in the shallow wells (upper graph in (a)), while for ∆V sp > 0 it corresponds to a superposition of the degenerate p x -and p y -orbitals in the deep wells (lower graph in a). A two-dimensional band calculation for the potential in Eq. (3) yields the true 2nd, 3rd and 4th Bloch bands. In Fig. 5 (b) and (c), these bands are plotted along a roundtrip in the 1st Brillouin zone connecting the points Γ, X + , M, X − , Γ (indicated in (d)). The average well depth is V 0 = 6 E rec and ∆V = 5.7 E rec , which corresponds to ∆V sp ≈ 0 in the TB picture in (a). In (b) the idealized C 4 symmetric case is shown, compared in (c) to the experimentally implemented case with y = 1. The band degeneracies arising at the Γ-and M -points in (b) are lifted in (c) due to the lattice distortion. In either case the 2nd band possesses two inequivalent band minima at the X − -and X + -points with energies E(X ± ), which are degenerate even in the case (c) despite the broken C 4 -symmetry.
Only if y is tuned away from unity, this degeneracy is lifted. Experimentally, tuning of ∆E ≡ E(X − ) − E(X + ) is accomplished by tuning of y with ∆E( y ) ∼ 1− y quantified by a band calculation for the potential in Eq. (3). In the vicinity of the X − -and X + -points, where the condensed atoms reside, the bands in (b) and (c) are approximately equal.
A III. Bosonic multi-band Hubbard model. We employ the multi-band Hubbard-Hamiltonian with the summation indices µ, ν ∈ {−1, 1}, σ ∈ {x, y}, d ν ≡ e x + νe y , and e x,y as shown in Fig. 1 (a) of the main text. This Hamiltonian accounts for all possible tunneling processes between nearest-(NN) and next-nearest neighbors (NNN), as is indicated in Fig. 5 (f): NN-tunneling between s-orbitals and p-orbitals (J sp ), NNN-tunneling between s-orbitals (J ss ), NNN-tunneling between p x -orbitals or p y -orbitals (J ), NNN-tunneling between p x -and p y -orbitals (J ⊥ ), a term accounting for a potential energy difference between p-and s-orbitals (∆V sp ) and the on-site interactions for p-orbitals (U p ) and s-orbitals (U s ), respectively, withn s,R ≡ŝ † Rŝ R ,n p,R ≡p † x,Rp x,R +p † y,Rp y,R andL p,R ≡ i(p † x,Rp y,R −p † y,Rp x,R ). Furthermore, an extra term scaling with J ,a is included, which introduces a quadrupolar anisotropy of the tunneling between p-orbitals with respect to the d ν -directions. This term permits to model the tuning of the energy difference ∆E between the minima of the second band. Diagonalizing the kinetic part ofĤ in momentum space at X ± yields the relation ∆E ≈ 8ν p J ,a with the relative occupation ν p of the p-orbitals given as ν p ≈ 1 2 [1 + ∆V sp /(32J 2 sp + ∆V 2 sp ) 1/2 ].
A IV. Determination of hopping parameters and phase diagram. The Hubbard Hamiltonian of Eq. (4) yields three TB bands associated with the three local orbitals involved. Setting J ,a = 0 we match the lowest two TB bands with the full bands in Fig. 5 (b). This lets us determine the tunneling parameters as J ss ≈ 0, J sp = 0.12 E rec , J = 0.07 J sp , J ⊥ = 0.15 J sp , and ∆V sp = 0.3 J sp . The resulting TB bands are plotted in Fig. 5 (e) showing good agreement with the full bands. The smaller bandwidth of the 3rd TB band as compared to the corresponding 4th full band indicates the proximity of higher bands (p-orbitals in the shallow wells) in the true lattice potential, not accounted for in the TB-description. A non-zero value of J ,a introduces an imbalance for NNN-tunneling between p-orbitals along the x + y and x − y-directions, which acts to lift the degeneracy of the two band minima of the lowest TB band. The phase diagram in Eq. (1) of the main text is derived as follows: the Hamiltonian in Eq. (4) including J ,a is rewritten, using a composition of the lattice via plaquettes with four sites. We diagonalize the kinetic energy in momentum space and assume that the system is condensed at the two inequivalent energy minima of the lowest TB band. Applying a mean-field approximation, we calculate the interaction and minimize the total energy with respect to the relative phase between the two possible condensation states and with respect to their relative weights.
A V. Momentum spectra and band populations. Momentum spectra are obtained by rapidly (< 1 µs) extinguishing the lattice potential, permitting a free expansion of the atomic sample during 30 ms, and subsequently recording an absorption image. Band populations are measured as follows: the population of the n-th band is transferred into the n-th Brillouin zone by adiabatically terminating the lattice potential in 400 µs, followed by a ballistic expansion of 30 ms. An absorption image of the atomic density distribution is recorded and the populations in the different Brillouin zones are counted.
A VI. Tuning of ν p , effect of band relaxation. The required tuning of the relative population of the p-orbitals ν p is accomplished via adjustment of ∆V . Since this tuning is essential, we have studied it in some detail. The value of ν p corresponding to some β and y (and hence ∆V and ∆E) is obtained by integrating |ψ ± (r)| 2 over those fields of the chequerboard lattice comprising the deep wells with the Bloch functions ψ ± (r) at the condensation points X ± derived from a band calculation employing the potential of Eq. (3). ν p scales monotonously with ∆V (as shown by the solid black trace of Fig. 6 (a)) and it is practically independent of the lattice distortion parameter y and hence of ∆E. This behavior is consistent with the approximate analytic expression derived from the Hubbard model in Sec. A III.
Changing ∆V also has a significant impact on the time-scale of band relaxation, as is shown by the data points in Fig. 6 (a). The band lifetime is found to be maximal (≈ 230 ms) if most atoms reside in the local s-orbitals (ν p ≈ 0). This is expected, since for these atoms there is no local state with lower energy available, which could give rise to relaxation. The initial preparation of the condensate is carried out in this configuration. Following a holding time of 80 ms to reach complete equilibrium, a subsequent adiabatic increase of ν p during 3 ms is directly visible in the momentum spectra: the higher order Bragg-peaks become more populated due to the increased contribution of the local p-orbitals, which due to their nodal structure comprise higher momenta. This is shown in Fig. 6 (b) and (c) with the corresponding calculation in Fig. 6 (d) yielding good agreement. Increasing values of ν p are accompanied by faster band decay. Since a single tunneling process between adjacent wells is sufficient for adjusting ν p , adiabaticity only requires tuning times of a few tunneling times of about 1 ms.
The lifetime of the atoms in the second band ( Fig. 6 (a) red disks) is measured as follows: after forming the condensate at ∆V = 3E rec , ∆V is tuned in 3 ms (sufficiently slow to permit tunneling) to the desired value. After a variable duration the population of the second band is determined (see V.). The wings of the decaying populations are fitted by exponentials with the 1/e-times plotted as the red disks in Fig. 6 (a). The lifetimes of the condensed fraction (blue squares in Fig. 6 (a)) are obtained by an analog procedure, however counting the number of atoms in the lowest order Bragg peaks of a momentum spectrum. the linear polarization of the incoming beam in the y-branch of the lattice potential is rotated such that n + ≈ n − , which corresponds to ∆E( y = 1) = 0. Arbitrary values of ∆E are adjusted by rotating the polarization away from this position by precisely quantified amounts and determining the corresponding values of ∆E via a band calculation for the potential in Eq. (3). The estimated error in the determination of y is ∆ y ≈ 2.5 · 10 −3 .
A VIII. Determination of n 2 p U p , n 0 ∆E and ν dif,th . We account for the isotropic harmonic potential with approximately 40 Hz trap frequency superimposed upon the two dimensional lattice of Eq. (3) with η = x = y = 1 and tubular lattice sites extending along the z-direction. The total number of condensed atoms, determined by fitting Gaussians to all visible Bragg peaks of a momentum spectrum and counting the atoms, is N ≈ 1.5 · 10 4 before band relaxation sets in. Assuming a Thomas-Fermi density distribution along the weakly confined z-direction, we calculate the local number of atoms per unit cell n 0,R at the Bravais lattice site R, and hence the local number of particles in the s-orbitals and p-orbitals, n s,R = ν s n 0,R and n p,R = ν p n 0,R , respectively. In the center of the lattice n 0,R ≈ 50. In order to determine the local interaction energies per particle in the p-orbitals and s-orbitals at site R, U p,R and U s,R , the local wells at the A-and B-sites are approximated to 6th order and the corresponding wave-functions of the local p-orbitals and s-orbitals are calculated. With this input, we determine the local energies n 2 p,R U p,R and n 0,R ∆E. By applying Eq. 2 of the main text, the corresponding local value of θ and hence ν dif,th = sin 2 (θ) − cos 2 (θ) is obtained. Averaging over all positions in the lattice yields n 2 p U p , n 0 ∆E, used in Fig. 7 (a) and ν dif,th plotted in Fig. 3 of the main text.
A IX. Measurement of ν dif and comparison to theory ν dif,th . In order to obtain Fig. 3 (a)-(f) of the main text, in Eq. (3) y is adjusted to the six different values 1 − y ∈ {0.019, 0.014, 0.009, 0.004, 0.001, −0.001} and β is varied for each setting. Momentum spectra are recorded and the number of atoms n ± in the two zero order Bragg peaks corresponding to X ± are counted. A data point in Fig. 3 (a)-(f) represents an average over 8 measurements. We thus obtain ν dif ≡ ( n + − n − )/( n + + n − ) for different combinations of the parameters β and y . From β one obtains ∆V (see Sec.A I) and thus ν p (see Sec.A VI). The value ∆E is determined from a band calculation for the potential in Eq. (3) (see Sec.A VII). For each data point, n 0 and n p = ν p n 0 are derived at the time of the measurement, thus accounting for band relaxation loss. With this the energies n 2 p U p , n 0 ∆E and the occupation difference ν dif,th are obtained (see Sec.A VIII). The data points in each of the graphs (a)-(f) in  Fig. 3 (b), the local value of ν dif,th (R) is plotted across the lattice. The lattice is parametrized by the Bravais vectors R = µ d+ + ν d−.
identified here in Fig. 7 (a). For the case of the red dashed path corresponding to Fig. 3 (b), the local value of ν dif,th is plotted across the lattice in Fig. 7 (b). This confirms that, particularly for the data points in the center of the path, finite size effects are in fact substantial and have to be accounted for.
A X. Mean field predictions of ν dif,th for real superposition state and incoherent mixture. We have also minimized the total energy for the real coherent superposition sin(θ) |ψ + ± cos(θ) |ψ − (scenario C) and for the incoherent mixture of spatially superimposed condensates sin 2 (θ) |ψ + ψ + | + cos 2 (θ) |ψ − ψ − | (scenario D) in order to calculate the local mixing angle θ as a function of ∆E and ν p and hence to determine ν dif,th = sin 2 (θ) − cos 2 (θ). For scenario C we obtain the simple result that ν dif,th is constrained to the values ±1 depending on the sign of ∆E with no dependence upon ν p . For scenario D, values of ν dif,th deviating from ±1 require ν p > 2/3 and satisfy