Spin-Induced Orbital Frustration in a Hexagonal Optical Lattice

Complex lattices provide a versatile ground for fascinating quantum many-body physics. Here, we propose an exotic mechanics for generating orbital frustration in hexagonal lattices. We study two-component (pseudospin-$1/2$) Bose gases in $p$-orbital bands of two-dimensional hexagonal lattices, and find that the system exhibits previously untouched orbital frustration as a result of the interplay of spin and orbital degrees of freedom, in contrast to normal Ising-type orbital ordering of spinless $p$-orbital band bosons in two-dimensional hexagonal lattices. Based on the classification by symmetry analysis, we find the interplay of orbital frustration and strong interaction leads to exotic Mott and superfluid phases with spin-orbital intertwined orders, in spite of the complete absence of spin-orbital interaction in the Hamiltonian. Our study implies many-body correlations in a multi-orbital setting could induce rich spin-orbital intertwined physics in complex lattice structures.

Introduction. Frustration, caused by lattice geometry or long-range interactions [1], often demonstrates new emergent structures for strongly correlated quantum matter, where exotic phases appear, such as spin liquid phases [2][3][4], and topological states [5]. In these complex magnetic systems, a spin cannot find an orientation which simultaneously favors all the spin-spin interactions with its neighbors, as a result of the spin frustration of anti-ferromagnetic spin states [6]. In addition to spin, other fundamental physical quantities, such as orbital degree of freedom, can show frustration in a complex lattice as well, which provides an opportunity to investigate new orbital physics [7][8][9][10][11][12][13][14]. More interestingly, the combining of spin and orbital degrees of freedom in many-body systems [15][16][17][18][19][20][21][22] potentially supports unconventional frustrations and exotic phenomena in complex lattices.
p-orbital ultracold atomic gases provide a new platform for investigating the interplay of spin and orbital degrees of freedom, such as spin-orbital coupling in a square lattice [23,24], where the key element is onsite interactions for building many-body correlations. However it remains unknown how the exotic mechanism of interaction driven spin-orbit coupling carries over to other lattices, such as hexagonal lattices, where one has successfully loaded ultracold 87 Rb atoms into sp-orbital bands of a hexagonal lattice and observed Potts-nematic superfluid [25]. Generally, the orbital degree of freedom develops Ising-type order for ultracold spinless bosons in p-orbital band of hexagonal lattices [12], whereas the bosons demonstrate ferromagnetic order for spin degree of freedom in manybody ground states [26,27], indicating no frustrations for spin or orbital degree of freedom in hexagonal lattices.
We propose a novel mechanics for generating orbital frustrations in spinful bosonic systems in hexagonal lattices. Motivated by recent achievements of p-(a) -1.5 π π π π π π -π -π -π -π -π -π A B e 2 e 1 e 3 e x e y k y k x P ↑,r P ↓,r P r (Color online) (a) The geometry of the twodimensional hexagonal lattice with lattice vector em (left), and the first Brillouin zone (right). (b-d) The lowest p-orbital band of the noninteracting spinful bosonic system in a twodimensional hexagonal lattices, with (b) t ⊥ = 0, (c) t = t ⊥ , and (d) t = 0. (e-f) Cartoons of real-space orbital polarization Pσ,r for strongly interacting many-body phases in porbital bands of the two-dimensional hexagonal lattice, where (e) spinless bosons demonstrate out-of-plane Ising-type orbital order, and (f) spinful cases inplane orbital textures, indicating spin-induced orbital frustrations.
orbital physics of hexagonal lattices [25], we study twocomponent bosonic gases loaded into the p-orbital bands of a two-dimensional (2D) hexagonal lattice. Induced by the interplay of spin and orbital degrees of freedom, we observe orbital frustrations in 2D hexagonal lat-arXiv:2101.05367v2 [cond-mat.quant-gas] 24 Aug 2021 tices, namely inplane circulating orbital textures, in contrast to out-of-plane Ising-type orbital ordering for the spinless case. The emergent frustrations induce exotic spin-orbital intertwined orders in both strongly interacting Mott-insulating and superfluid states, which spontaneously break time-reversal symmetry, in contrast to the time-reversal-even order established for the square lattice [23]. This fascinating order can be probed with timeof-flight measurements in optical-lattice experiments.
Model and Hamiltonian. We consider bosonic atoms prepared in two hyperfine ground states which are loaded into the p-orbital bands of a 2D hexagonal optical lattice. The corresponding annihilation operators for the bosonic particles are denoted as p νσ,r , with r the position of lattice sites, ν = x, y labeling the p x and p y orbital degrees of freedom, and σ =↑, ↓ the two pseudospin (hyperfine) states. For a compact notation, we introduce two-component spinors Φ ν,r = [p ν↑,r , p ν↓,r ] T . Under the single-mode approximation, the 2D hexagonal lattice can be described by a generalized Bose-Hubbard model, Here, the hopping amplitudes between two nearestneighboring p-orbitals along the parallel and the perpendicular directions are denoted as t and t ⊥ , respectively. The unit vectors e 1,2 = ± √ 3 2 e x + 1 2 e y and e 3 = −e y are illustrated in Fig. 1(a). The lattice vectors are introduced as a 1 = e 1 − e 2 , and a 2 = e 1 − e 3 . The position of A and B sublattices locates at r = l 1 a 1 + l 2 a 2 + e 2 , and r = l 1 a 1 + l 2 a 2 − e 1 , respectively, with l 1 and l 2 being integer numbers. The lattice annihilation operators Φ m,r ≡ (p x↑,r e x + p y↑,r e y ) · e m , (p x↓,r e x + p y↓,r e y ) · e m T for hopping t , and Φ m,r ≡ (p x↑,r e x + p y↑,r e y · e m , (p x↓,r e x + p y↓,r e y · e m T with e 1,2 = − 1 2 e x ± √ 3 2 e y and e 3 = e x for hopping t ⊥ . The chemical potential µ controls the number density, and U 0,2 denotes the interaction strengths [23]. We introduce the occupation number operator, n r = ν Φ † ν,r Φ ν,r , the angular momentum, L z,r = i[Φ † x,r Φ y,r − Φ † y,r Φ x,r ], the spin moment, S r = ν Φ † ν,r σΦ ν,r , and a spin angular-momentum coupled operator, LS z, For convenience, we also introduce U ↑ = U ↓ = U 0 + U 2 , and U ↑↓ = U 0 − U 2 , which correspond to intra-and inter-spin interactions, respectively, and are controllable in experiments via Feshbach resonances [28,29]. We first discuss the single-particle spectrum of the noninteracting p-band bosonic system, described by Eq. (S15) with U 0 = U 2 = 0. For hexagonal lattices, each unit cell consists of two nonequivalent sites, and each of which contains two p-orbitals, producing four bands in the system for each pseudospin σ [30,31]. In Fig. 1 we plot the dispersion of the lowest p-orbital band of the noninteracting spinful bosonic system in a 2D hexagonal lattice, with (b) t ⊥ = 0, (c) t = t ⊥ , and (d) t = 0, respectively. In the limiting case with t ⊥ = 0 or t = 0, the system supports a flat band over the entire Brillouin zone in the hexagonal optical lattice. The band flatness dramatically enhances interaction effects, which has been shown to support various exotic phases [32][33][34]. In the regime of t ≈ t ⊥ , however, the lowest p-orbital band is no longer rigorously flat but develops a finite width [ Fig. 1(c)]. In this dispersive band, bosons tend to form condensation. For the multi-orbital nature of this system, it is expected to develop rich spin and orbital orders in the strongly interacting limit with interaction strength comparable to the band width, which has been achieved in experiments [25].
Effective orbital-exchange model. We focus on the physics in the strongly interacting regime with t ,⊥ U ↑,↓ and U ↑↓ , and study the influence of virtual hopping processes on the orbital ordering. To quantify the orbital order, we introduce an orbital polarization vector P σ,r = P x σ,r , P y σ,r , P z σ,r , withP x σ,r ≡ 1 2 (p † xσ,r p xσ,r − p † yσ,r p yσ,r ),P y σ,r ≡ 1 2 (p † xσ,r p yσ,r + p † yσ,r p xσ,r ), andP z σ,r ≡ 1 2i (p † xσ,r p yσ,r − p † yσ,r p xσ,r ). For comparison, we first discuss the deep Mott physics of spinless p-oribtal band bosons in optical lattices with filling n = 1, which is given by an effective orbitalexchange model where ij denotes the nearest-neighbor sites i and j. J x = −3(t 2 + t 2 ⊥ )/2U , J y = 3t t ⊥ /U , and J z = 9t t ⊥ /U with U being the onsite interactions between spinless bosons [12]. In the regime of t ≈ t ⊥ , the J z term dominates, and the Ising exchange favors antiparallel configuration of nearby orbitals. This out-of-plane Ising-type order P z i survives in the hexagonal lattice [35]. For the spinful bosonic gases, we focus on the spinmixed regime with U ↑,↓ ≥ U ↑↓ > 0, and the lowest Mott lobe with filling n ↑ + n ↓ = 2. In the zero-hopping limit, the single-site ground state |ψ = 1 √ 2 |p x↑ , p y↓ − |p x↓ , p y↑ with P ↑,↓ = 0, indicating the absence of orbital polarization in the deep Mott-insulating regime. Away from the atomic limit but still in the strongly interacting regime, one can apply second-order perturbation theory and generate an effective Hamiltonian. To emphasis the interaction effects, we relax the superexchange processes for the spin-↑ component by freezing the spin-↓ one. In the regime of t ≈ t ⊥ , the effective  [38]. There are three quantum many-body phases with different spin-orbital intertwined orders, including the spin-orbital intertwined Mott-insulating phase (sMI) with P 2 = 0 by breaking time-reversal symmetry, the spinorbital intertwined superfluid phase (sSF) with both P 2 = 0 and ρc = 0 (definition in the main text) by breaking both time-reversal, Uc(1) gauge, and Us(1) spin-rotational symmetries, and the unordered Mott-insulating phase (uMI) in the absence of symmetry breaking. (c) The uMI-sMI-sSF phase transitions along µ/U ↑↓ = 0.45 [38].

orbital-exchange model reads
where . For the hexagonal lattice, this effective model is transferred to H eff = m,r∈A J xP↑,rP↑,r+e m withP ↑,r denoting the operator along the bond e m [36,37], indicating the inplane orbital frustration induced by the spin degree of freedom. In other words, the system may support a phase transition from an unordered Mott insulator to an orbital ordering phase determined by frustrated orbital exchange, which is totally distinct from the spinless case in the hexagonal lattice.
Frustration induced symmetry-breaking orders. In order to determine the quantum ground-state orders in the spinful p-orbital lattice (Eq. (S15)), we perform a bosonic dynamical mean-field theory (BDMFT) calculation recently adapted to honeycomb lattices [39,40] and p-orbital bands [23]. To accommodate stripe-like orders that spontaneously break lattice-translational symmetry, we apply the real-space BDMFT [41] to our system of spinful p-orbital bosons in the hexagonal lattice. The technical details are described in supplementary materials [35].
As expected for strongly interacting bosons in a lattice, we find both Mott insulating and superfluid phases in our system. Unlike the s-orbital spinful Bose-Hubbard model [26], the Mott phase of the two-component p-orbital bosons exhibits exotic orbital orders in the intermediate hopping regime. The orbital polarization vector P σ,r has spatial modulations with characteristic momenta Q ± = ± 4π 9 (e 1 − e 2 ) (P z σ,r = 0), forming a inplane circulating texture structure in the hexagonal lattice Fig. 1(f) and 3(c) , which is totally distinct from the spinless case [35] Fig. 1(e) , indicating orbital frustration induced by spin-dependent coupling. On the A-and B-sublattices, we have P σ,r ∼ s σ [cos θ A,r , sin θ A,r ] and ∼ s σ [cos θ B,r , sin θ B,r ], respectively, with s ↑,↓ = ±, and θ A,r = Q + · r, θ B,r = Q − · r + π, or the other way around. The sublattice symmetry is spontaneously broken. The orbital frustration also dominates in the strongly correlated superfluid phase, as shown in Fig. 2

and 3.
Orbital frustration stabilizes novel long-range order in our spinful p-orbital lattice, where there are 16 potential spin-orbital intertwined orders, and can be classified according to their symmetry transformations [35]. Three different types of spin-orbit intertwined orders could potentially emerge for the hexagonal lattice, including A 1 -Odd, A 2 -Even, and E-Odd, which are labeled according to their transformations under the C 3v and T symmetries [35]. In the Mott regime, we find the E-Odd order, with local observables d z where the phase is referred to as the spin-orbital intertwined Mott-insulator (sMI) below. This order forms a two dimensional representation (E) of the lattice rotation C 3v symmetry group, and breaks the time-reversal symmetry, in sharp contrast to spin angular-momentum coupled order as established for spinful p-orbital bosons in the square lattice [23]. To accommodate the realspace distribution, a new order parameter is defined as represents a counter-clockwise rotation matrix by an angle (. . .), and N lat is the total number of lattice sites. Its square, P 2 , is a C 6 rotation invariant, which characterizes the strength of the spin-orbital intertwined order, as shown in Fig. 2(c) and 3(d).
In the superfluid regime, besides the d z x 2 −y 2 ,r and d z xy,r orders, we also find that d , as a consequence of spontaneous symmetry breaking of U c (1) and U s (1), where the many-body phase is denoted as the spin-orbital intertwined superfluid phase (sSF). The local orders in the sSF phase are still E-Odd orders, which correspond to the E representation of the lattice C 3v symmetry group, and remain time-reversal odd [35], in distinction to the emergent ordering as found in the square lattices [8,12,23,42]. In addition, we find the bosons in the sSF phase con- dense into a state that is a linear combination of singleparticle states at two inequivalent momenta Q + and Q − , as shown in Fig. 4(a)(b), indicating lattice-translationalsymmetry breaking. The corresponding superfluid order φ νσ,r ≡ p νσ,r takes a form of Here, ρ c ≡ νσ,r |φ νσ,r | 2 /N lat represents the averaged number of condensed particles per lattice site. We have λ + x = 1/2, λ + y = i/2, and λ − ν = [λ + ν ] * , to accommodate the local p x ± ip y character of the Bloch function at the Q ± -momentum points, and γ is a real positive number to be consistent with the orbital order shown in Fig. 2  and 3.
Many-body phase diagrams. We summarize our zerotemperature phase diagrams of spinful p-orbital bosons of 2D hexagonal lattices in Fig. 2, and 3, obtained by bosonic dynamical mean-field theory. Here, averaged particle number of the condensate ρ c , and spin-orbital intertwined order P 2 , are utilized to identify the quantum phase transitions. Without loss of generality, we mainly focus on the case with interactions U ↑ = U ↓ = U ↑↓ , which are good approximations for 87 Rb atoms.
Our calculated filling-dependent phase diagrams are presented in Fig. 2, as a function of chemical potential and hopping amplitudes (a) t ≡ t = t ⊥ , and (b) t ≡ t = 5t ⊥ , respectively. In the lower hopping regime, two Mott-insulating phases are found: an unordered Mott insulator (uMI) in the absence of any symmetry breaking, and a sMI phase with time-reversal T and SU (2) spin-0 6000 0 6000 -π π -π π π π -π -π n ↑,k n ↓,k rotational symmetries being spontaneously broken (the U s (1) spin-rotational symmetry develops as a subgroup instead). In the larger hopping but still strongly interacting regime, the atoms are delocalized and a sSF phase is energetically favorable, which breaks time-reversal T , U s (1) spin-rotational, and U c (1) gauge symmetries. In Fig. 2(c), the order parameters as a function of hopping amplitudes are presented to characterize the uMI-sMI-sSF transitions for µ/U ↑↓ = 0.45. We remark here that the phase transitions can be continues or discontinues, depending on the parameters of the model, which is beyond the scope of the present work. Normally, the hopping amplitudes of p-orbital band bosons can be tuned separated in the perpendicular and parallel directions, respectively. The hopping-dependent phase diagrams are shown in Fig. 3(a)(b) for different interactions with fixed filling n = 2. We observe that there are also three many-body quantum phases for the parameter regime studied here, including the uMI, sMI, and sSF phases, where the spin-orbital intertwined phases (sMI and sSF) are stable and explore a large part of the phase diagrams, indicating large opportunities for experimentally observing these many-body quantum phases. Interestingly, we observe that pseudospin ordering S pν ,r ≡ 1 2 p † νσ,r F σσ p νσ ,r of the spin-orbital intertwined phase is colinear with Z 3 order in the hexagonal lattice, which coincides with the frustration-induced orbital textures P σ,r via breaking lattice-translational symmetry, as shown in Fig. 3(c). Here, F σσ is the spin matrices for pseudospin-1/2. We also study the influence of the inter-spin interaction on the phase diagrams. As shown in Fig. 3(b), we observe the Mott-insulating phase expands for U ↑ = U ↓ = 2U ↑↓ .
Experimental proposal and detection. For the experimental detection, we expect the superfluid phase having the spin-orbital intertwined order is accessible to the present experiments [25,[43][44][45][46]. The symmetry-breaking pattern in the momentum distribution (Fig. 4) can be readily tested. Each atomic pseudospin component predominantly condense at Q + and Q − , respectively, and the symmetry between them is spontaneously brokenthe symmetry breaking pattern is opposite for the two components. The spatial dependence of the spin-orbital intertwined order in the superfluid phase is determined by the phase coherence between the two momenta, which can be probed by Raman-assisted interference [12].
Conclusion. We study two pseudospin-component bosons loaded into p-orbital bands of a two-dimensional hexagonal lattice. We find exotic circulating orbital textures with the two-pseudospin components rotating in opposite directions, which are stabilized by orbital frustration in hexagonal lattices. Accommodating this spontaneous orbital ordering, we find a spin-orbital intertwined order in both superfluid and Mott insulating phases, which corresponds to the two-dimensional E representation of the lattice C 3v group, and is timereversal odd. Our work opens up wide opportunities for spontaneous spin-orbital coupled physics in hexagonal optical lattices. For example, the interplay of spin degrees of freedom with the renormalization stabilized Potts nematicity in an sp 2 -orbital hybridized hexagonal lattice [25] is expected to support even more exotic spontaneous spin-orbital coupled phenomena.

BOSONIC DYNAMICAL MEAN-FIELD THEORY
The long-range order of many-body ground-state states of bosonic gases in p-orbital bands loaded into a 2D hexagonal optical lattice, described by Eq. (1) in the main text, can be obtained from diagonal Green's functions G µν,ij (τ, τ ) = p µ,ri (τ )p † ν,rj (τ ) , and off-diagonal Green's functions G µν,ij (τ, τ ) = p µ,ri (τ )p ν,rj (τ ) . For the essential four-component bosonic system studied here, they are 32 independent Green's functions. In practice, we can combine different Green's functions to obtain the long-range orders listed in Table I.
To obtain these Green's functions, we utilize a bosonic version of dynamical mean-field theory (BDMFT) on the two-dimensional situation. BDMFT has been developed to provide a non-perturbative description of zero-and finitetemperature properties of the Bose-Hubbard model [S23, S41, S47-S53], whose reliability of this approach has been compared against the quantum Monte-Carlo simulations [S54]. Actually, BDMFT has already been implemented in the honeycomb-lattice structure to study the bosonic Haldane model [S39] and Kane-Mele-Hubbard model [S40]. Recently, a four-component bosonic dynamical mean-field theory has been developed to study the multi-species bosonic system in p-orbital bands of a two-dimensional square lattice [S23]. In this paper, we modify our method to study the bosonic ultracold gases loaded into the frustrated hexagonal lattices.
As in fermionic dynamical mean-field theory [S55], the main idea of the BDMFT approach is to map the quantum lattice problem with many degrees of freedom onto a single site. To obtain a proper local model that is approximate to the original lattice Hamiltonian, we perform an effective integration over all off-site degrees of freedom and keep only terms of subleading order. Finally, we find that the local Hamiltonian is given by a bosonic Anderson impurity Hamiltonian where : . . . : denotes normal ordering, t σ is the nearest-neighbor hopping amplitude for the four components denoted as σ, the chemical potential and interaction terms are directly inherited from the Hubbard Hamiltonian. The bath of condensed bosons is represented by the Gutzwiller term with superfluid order parameters φ σ for the component σ.
The bath of normal bosons is described by a finite number of orbitals with creation operatorsâ † l and energies l , where these orbitals are coupled to the impurity via normal-hopping amplitudes V σ,l and anomalous-hopping amplitudes I. Local symmetry-breaking orders for two-component bosons loaded into p-orbital bands of a hexagonal lattice. The bilinear operators Φ † ν σzΦ ν keep the Us(1) spin-rotational symmetry preserved, whereas the operators Φ † ν σx,yΦ ν break this symmetry. Note here that we omit the lattice position r in the operator Φν,r, to shorten the donation. Operators The anomalous hopping terms are needed to generate the off-diagonal elements of the hybridization function. Note here that we useb σ to denote the four-component bosonic annihilation operator, to shorten the notation of the function.
To solve the Anderson Hamiltonian, we apply exact diagonalization as a solver to solve the local problem [S55, S56]. After diagonalization, we obtain the eigenstates, eigenenergies, and local Green's functions in the Lehmannrepresentation where ω n denotes Matsubara frequency, the superfluid order parameter φ σ = b σ 0 , and the notation . . . 0 means that the expectation value is calculated in the cavity system [S48]. Then, the local self energy for each site can obtained via Dyson equation: where G A,σσ (iω n ) denotes the noninteracting Weiss Green's function of the Anderson impurity site. After obtaining the self energy for all the sites, we can employ the Dyson equation in real-space representation to compute the interacting lattice Green's function: where the noninteracting lattice Green's function G −1 0 (iω n ) = (µ+iω n )1−t, with the matrix of hopping t determined by lattice structures. Note here that the site-dependence of the Green's functions and self energies are shown by boldface quantities that denote a matrix form with site-indexed elements.
The self-consistency loop is solved as follows: starting from an initial choice for the Anderson parameters and the superfluid order parameters, the Anderson Hamiltonian is constructed in the Fock basis and diagonalized to obtain the eigenstates and eigenenergies. The eigenstates and energies allow us to calculate the superfluid order parameter, the impurity Green's functions and self-energies, and then obtain the lattice Green's functions via Eq. (S5). Subsequently, new Anderson parameters are obtained, by comparing the new Green's functions with the old ones. With these new Anderson parameters, the procedure is iterated until convergence is reached.

SYMMETRY ANALYSIS
For the complexity of our spinful p-orbital lattice, we classify the potential long-range orders according to their symmetry transformations. The system has a U c (1) phase symmetry associated with particle-number conservation in ν,r Φ † ν,r Φ ν,r , and a U s (1) spin-rotational symmetry associated with spin conservation in ν,r Φ † ν,r σ z Φ ν,r . Apart from these continuous symmetries, it also has discrete lattice rotation C 3v , and time-reversal T symmetries. Different symmetry-breaking orders in this system would give rise to local quadratic observables of different symmetry properties listed in Table I. The time-reversal invariant spin angular-momentum intertwined order LS z,r = i[Φ † x,r σ z Φ y,r − Φ † y,r σ z Φ x,r ] has been found for spinor bosons loaded into p-orbital bands of a square lattice [S23]. For the hexagonal lattice, three different types of spin-orbit intertwined orders could potentially emerge, including A 1 -Odd, A 2 -Even, and E-Odd, which are labeled according to their transformations under the C 3v and T symmetries, as shown in Table I.

EFFECT ORBITAL-EXCHANGE MODEL OF SPINFUL BOSONIC GASES IN p-ORBITAL BANDS OF A HEXAGONAL LATTICE
In this part, we derive the effective orbital-exchange model for the spinful bosonic gases in p-orbital bands of a hexagonal lattice. We consider the lowest Mott lobe with filling n = n ↑ + n ↓ = 2 trapped in a single site. In this atomic limit with t = t ⊥ = 0, the system is described by We construct the basis as |p νσ , p ν σ , where ν = x, y and σ =↑, ↓, under which the above Hamiltonian can be expressed as After diagonalizable the matrix, we obtain the eigenstates and the eigenenergies From these eigenenergies, we find that the ground state is |ψ 9 in the regime of U ↑,↓ ≥ U ↑↓ > 0. With the above eigenvectors, we next calculate the expectation values of the operators for the ground state, where the orbital polarization operators are given by (S11) To emphasis the interaction effects, we relax the superexchange processes for the spin-↑ component by freezing the spin-↓ one with n ↓ = 1. In the low-hopping regime, one can obtain an effective orbital-exchange Hamiltonian from second-order perturbation theory where α and β label the ground states, γ denotes the excited states, and E are the eigenenergies of the zeroth order part H 0 , with V being the hopping terms treated as perturbation. Finally, we obtain an effective orbital-exchange Hamiltonian for the spin-↑ component where ij denotes the nearest-neighbor sites i and j of a lattice, J x = − (t 2 +t 2 ⊥ ) To generalize to the hexagonal lattice, we can rotate p-orbitals at an angle of θ with the x axis. Along this direction, the operators transform as b xσ = cosθ b xσ + sinθ b yσ and b yσ = −sinθ b xσ + cosθ b yσ . Accordingly, the orbital polarization operators P x σ = cos2θ P xσ + sin2θ P yσ , P y σ = −sin2θ P xσ + cos2θ P yσ , and P z = P z . In the regime of t ≈ t ⊥ , the J x term dominates in the effective model (S13), which is given by for the hexagonal lattice, with P ↑,r denoting the operator along the bond e θ directing at angle θ with the x axis.

PHASE DIAGRAM OF SPINLESS BOSONIC GASES IN p-ORBITAL BANDS OF A HEXAGONAL LATTICE
We consider spinless bosonic atoms prepared in the hyperfine ground state which is loaded into the p-orbital bands of a two-dimensional (2D) hexagonal optical lattice. The corresponding annihilation operators for the bosonic particles are denoted as p νσ,r , with r the position of lattice sites, ν = x, y labeling the p x and p y orbital degrees of freedom. Under the single-mode approximation, the 2D hexagonal lattice can be described by a generalized Bose-Hubbard model, Here, the hopping amplitudes between two nearest-neighboring p-orbitals along the parallel and the perpendicular directions are denoted as t and t ⊥ , respectively. The unit vectors e 1,2 = ± √ 3 2 e x + 1 2 e y and e 3 = −e y are illustrated in Fig. 1(a) in the main text. The lattice vectors are introduced as a 1 = e 1 − e 2 , and a 2 = e 1 − e 3 . The position of A and B sublattices locates at r = l 1 a 1 + l 2 a 2 + e 2 , and r = l 1 a 1 + l 2 a 2 − e 1 , respectively, with l 1 and l 2 being integer numbers. The lattice annihilation operators p m,r ≡ (p x,r e x +p y,r e y )·e m for hopping t , and p m,r ≡ (p x,r e x +p y,r e y ·e m with e 1,2 = − 1 2 e x ± √ 3 2 e y and e 3 = e x for hopping t ⊥ . µ denotes the the chemical potential, and U denotes the interaction strengths. We introduce the occupation number operator, n r = ν p † ν,r p ν,r , the angular momentum, L z,r = i[p † x,r p y,r − p † y,r p x,r ]. As shown in the effective spin model, the system can develop two types of long-range order in the Mott-insulating phases. We confirm this conclusion via numerical simulations, based on bosonic dynamical mean-field theory, and observe that there are three different quantum phases in the spinless bosonic gases in the p-orbital band of hexagonal lattices, including two Mott-insulating phases and one superfluid phase, as shown in Fig. S1. In the low-hopping regime t ≈ t ⊥ U , the J z term dominates and the system favors a Mott phase (L M I ) with staggered orbital angular momentum P z r = 0, and in the regime t t ⊥ or t ⊥ t , the J x term dominates and a new Mott phase appears with the absence of staggered orbital angular momentum. In the larger hopping regime, the atoms delocalize and the superfluid phase (L SF ) with staggered orbital angular momentum develops.

PHASE DIAGRAM OF HETERONUCLEAR BOSE-BOSE MIXTURES IN HEXAGONAL LATTICES
Up to now, theoretical calculations were mainly performed for the symmetric parameters with t and t ⊥ being identical for the two pseudospin components, which is the most relevant case for the present experiments. However, Bose-Bose mixtures may consist of different atomic species, such as 87 Rb and 7 Li, which generally do not have these symmetric properties. Without loss of generality, we still focus on U ↑ = U ↓ = U ↑↓ and t σ = t ⊥σ , where t σ and t ⊥σ are the hopping amplitudes along the parallel and perpendicular directions for pseudospin σ, respectively.
As shown in Fig. S2, a phase diagram is demonstrated as a function of t ↑ ≡ t ↑ = t ⊥↑ and t ↓ = t ↓ = t ⊥↓ , where four phases appears, including unordered Mott phases (uMI), spin-orbital intertwined Mott phases (sMI), spin-orbital smectic superfluid phase (sSF), and a new phase (sMI+sSF) when the hopping amplitude for one component is small and the other one relatively large. In this new phase (sMI+sSF), we observe that one species is the Mott-insulating state, and the other one the superfluid phase. In addition, we observe exotic orbital textures P σ,r for this sMI+sSF phase, with the Mott-insulating species rotating clockwise and the superfluid one anti-clockwise, or vice versa, as shown in the inset of Fig. S2, indicating the spin-orbital intertwined order of the many-body ground state emerging.