Emergent Orbital Skyrmion Lattice in a Triangular Atom Array

Multi-orbital optical lattices have been attracting rapidly growing research interests in the last several years, providing fascinating opportunities for orbital-based quantum simulations. Here, we consider bosonic atoms loaded in the degenerate $p$-orbital bands of a two-dimensional triangular optical lattice. This system is described by a multi-orbital Bose-Hubbard model. We find the confined atoms in this system develop spontaneous orbital polarization, which forms a chiral Skyrmion lattice pattern in a large regime of the phase diagram. This is in contrast to its spin analogue which largely requires spin-orbit couplings. The emergence of the Skyrmion lattice is confirmed in both bosonic dynamical mean-field theory (BDMFT) and exact diagonalization (ED) calculations. By analyzing the quantum tunneling induced orbital-exchange interaction in the strong interaction limit, we find the Skyrmion lattice state arises due to the interplay of $p$-orbital symmetry and the geometric frustration of the triangular lattice. We provide experimental consequences of the orbital Skyrmion state, that can be readily tested in cold atom experiments. Our study implies orbital-based quantum simulations could bring exotic scenarios unexpected from their spin analogue.

These recent developments open up unprecedented opportunities for orbital-based quantum simulations, by which the realizable quantum many-body states and phenomena could strongly deviate from the spin analogue due to the fundamental difference in their symmetry [16].
In this letter, we study interacting p-orbital bosons in a two-dimensional triangular optical lattice, which corresponds to the experimental setups of hexagonal lattices [2,4] with a large sublattice potential imbalance [17].This system is described by a multi-orbital Bose-Hubbard model [16,18].The quantum many-body phases of this system are investigated in the strong interaction regime by BDMFT and ED calculations.We find an orbital Skyrmion lattice state emerges in the Mott insulating regime.The Skyrmion lattice state has a composite chirality that spontaneously breaks the timereversal symmetry.In the orbital setting, two different types of Skyrmion lattices occur simultaneously in the degenerate quantum many-body ground states, which is in sharp contrast to the Dzyaloshinskii-Moriya scenario systems [19][20][21][22][23]. Near the Mott-superfluid transition, we find the chiral Skyrmion lattice melts to quantum states with stripe and ferro-orbital orders.
Model and Hamiltonian.The system of interacting spinless bosonic atoms, loaded into the p-orbital bands of a two-dimensional (2D) triangular optical lattice, can be described by a multi-orbital Bose-Hubbard model in the tight-binding limit [16,18] where t and t ⊥ denote the hopping amplitudes between two nearest-neighboring p-orbitals along the parallel and the perpendicular directions, respectively.The lattice annihilation operators p m,r ≡ (p x,r e x + p y,r e y ) • e m with the unit vectors e 1 = e x and e 2,3 = ± Here, e m is shown in Fig. 1(a), and p x,r (p y,r ) denotes the annihilation operator for the p x (p y ) orbital degree of freedom at site r.µ is the chemical potential, n r = ν p † ν,r p ν,r with ν = x, y, and U , U 1 and U 2 denote the interaction strengths with U − 2U 1 = U 2 as a result of symmetry analysis for the triangular lattice.In the deep lattice limit, the harmonic approximation of the Wannier function implies U = 3U 1 = 3U 2 , and consequently the interactions (U and U 1,2 ) take a simplified form as ), with the orbital angular momentum L z,r ≡ −i(p † x,r p y,r − p † y,r p x,r ) [12].Weak interaction limit.To understand the many-body phenomena, we first discuss the physics in the weakly interacting superfluid regime with t ,⊥ U and U 1,2 , where the bosons are expected to condense.For triangular lattices, the Brillouin zone forms the shape of a regular hexagon with the edge length 4π/3, where the lattice constant is set to be the unit of length.The singleparticle spectrum of the noninteracting p-band bosonic system is shown in Fig. 1(b),(c), where we plot the dispersion of the lowest p-orbital band of a 2D triangular lattice.For t ⊥ = 0, the system supports three degenerate minima located at M 0 = (0, 2π/ √ 3) and M ± = (±π, π/ √ 3) [Fig.1(b)].For t = 0, band minima move to the center of Brillouin zone [Fig.1(c)].Due to the competition between hopping t and t ⊥ , it is expected to develop rich orbital orders with the appearance of orbital angular momentum L z,r = 0. We find that the system supports two different types of condensates in the weakly interacting regime.The ground state in the limit of t ⊥ t is a stripe superfluid phase [S SF phase in Fig. 1(d)], which is described by In the regime of t t ⊥ , the system condenses at the center of Brillouin zone and demonstrates a ferro-orbital order [H SF phase in Fig. 1(d)], with Here, |0 is the vacuum state, N denotes the total number of particles, σ r = ±1 is the sign of the staggered orbital angular momentum, β r = 0 (σ r = 1) or π/2 (σ r = −1) in the stripe direction, and β r = π in the homogeneous direction.Note here that α = π/6 in the regime of t t ⊥ and t ,⊥ /U 1 [24] and the bosons condense at two of the three degenerate minima as observed in the experiments [17], which is consistent with our numerical simulations.
Orbital Skyrmion lattice state at strong interaction.Considering the tunability of Hubbard parameters experimentally [11], we extend our study to the strongly interacting regime of the spinless p-orbital bosons in the triangular lattice, described by Eq. ( 1).To analyze quantum ground states of the many-body system, a bosonic version of dynamical mean-field theory is implemented.We remark here that BDMFT is an extension of fermionic dynamical mean-field theory, and suitable to treat strongly correlated systems for the full range of couplings from Mott insulator to superfluid.To accommodate long-range orders that spontaneously break lattice-translational symmetry, we generalize a real-space BDMFT [25] for our system of spinless p-orbital bosons in the triangular lattice.The technical details can be found in Supplementary Materials [26].
At strong interactions, bosons form a Mott insulating state.Through BDMFT calculation, we find this Mott state develops spontaneous orbital polarization forming a Skyrmion lattice that breaks timereversal and lattice-translational symmetries.Interestingly, the BDMFT calculation reveals two different types of Skyrmion textures, i.e., Skyrmion and anti-Skyrmion lattices as shown in Fig. 1(e), where an orbital polarization vector S r = S x r , S y r , S z r is defined, with ).The two Skyrmion lattice phases are connected by the T × I symmetry, with T and I being time-reversal and space-reflection (p x,r → −p x,r ) symmetries, respectively.We remark here that the emergent orbital Skyrmion texture observed here is solely induced by onsite interactions, and the underlying physics is the interplay of p-orbital symmetry and geometric frustration of the triangular lattice, captured by an effective orbital-exchange model (this model will be discussed later).To show the robustness of Skyrmion texture against quantum fluctuations, we map out the full t − t ⊥ phase diagram for filling n r = 1 and interactions U = 3U 1 = 3U 2 in the framework of BDMFT, as shown in Fig. 1(d).We find that the orbital Skyrmion lattice is robust against quantum fluctuations and explores a wide regime in Mott phases.Only for sufficiently large asymmetry between the two hopping amplitudes, two other Mott phases develop instead, including stripe-orbital phase (S MI ) breaking time-reversal symmetry, and ferro-orbital phase (H MI ) respecting time-reversal symmetry (realspace orbital textures shown in Fig. S1 [26]).All of these orbital-ordered phases are found to persist up to the superfluid transition.After the Mott-superfluid transition, the system demonstrates two superfluid phases, where one is a stripe-orbital phase (S SF ) breaking time-reversal, lattice-translational and rotational symmetries [24], and the other a ferro-orbital phase (H SF ) breaking timereversal symmetry, consistent with Eq. ( 2) and (3).Note here that the phase diagram is symmetric upon orbital interchange in the low-hopping regime, which is also manifested in the effective orbital-exchange model [Eq.( 4)].
To quantify phase boundaries in Fig. 1(d), we intro-duce superfluid order φ ν ≡ r | p ν,r |/N lat , stripe order Θ stripe ≡ S z + − S z − , and scalar spin chirality χ = S r • (S r+e1 × S r+e2 ) , where N lat is the number of lattice sites, and S z + ( S z − ) denotes the z-component of orbital polarization per site on the stripe with positive (negative) value.We clearly observe nonzero values of scalar spin chirality χ = 0 in the Skyrmion phase, as shown in Fig. 2(a).For larger hopping amplitudes, we find a Mott phase transition from the Skyrmion to the S MI phase, indicated by the absence of χ = 0 and the appearance of Θ stripe = 0.The corresponding contour plots of static spin structure factor F z (k) = i,j e ik•(ri−rj ) S z ri S z rj are shown in the inset of Fig. 2(a) for different Mott phases, where i and j denote the lattice sites.Increasing hopping amplitudes further, atoms delocalize with the coexistence of stripe order Θ stripe = 0 and superfluid order φ x,y = 0.In addition, we observe a superfluid phase transition from a ferro-to a stripe-orbital order, as shown in Fig. 2(b).We remark here that the phase transitions are found to be discontinues within BDMFT.
Orbital exchange and effective model construction.To explain the underlying mechanism of the orbital textures in the deep Mott regime with t ,⊥ U and unit filling, we construct an effective orbital-exchange model for Eq. ( 1).In our case, the orbital-exchange interactions arise from the virtual hopping processes induced by t ,⊥ , and are obtained from the perturbative expansion of tunneling processes up to third order (thirdorder expansion will be justified later).By introducing the projection operator P to describe the Hilbert space of the singly occupied Mott state, the effective Hamiltonian reads H eff P|ψ = EP|ψ , where −1 QH t P with Q = 1 − P, and H t being the hopping part of Eq. (1).Due to E ∼ t 2 /U , we obtain QHQ − E ≈ QHQ.
Generally, the orbital polarization operator S r changes with bond orientations [24,[27][28][29].It is convenient to introduce the rotation direction e θ for orbital polarization operator, i.e., e x θ = cos (2θ) e x + sin (2θ) e y , e y θ = −sin (2θ) e x + cos (2θ) e y and e z θ = e z , for a bond directing at angle θ with the x axis.With the definition above, we finally obtain an anisotropic orbital-exchange model for the triangular lattice system with interactions U = 3U 1 = 3U 2 [26], where {u, v, w} = {x, y, z}, θ m is the angle with the x axis for the bond e m , and J and J denote the orbitalexchange terms from second-O(t 2 ,⊥ /U ) and third-order O(t 3 ,⊥ /U 2 ) tunneling processes, respectively.In the absence of third-order interactions, the effective model reduces to an XY Z model with orbital-exchange pa- rameters J x = −3(t 2 + t 2 ⊥ )/2U , J y = 3t t ⊥ /U , and J z = 9t t ⊥ /U [16].Generally, J x dominates in the regime of t ⊥ t or t t ⊥ , where in-plane ferro-orbital order develops, and J z dominates the remains (t ≈ t ⊥ ), where the system favors out-of-plane Ising-orbital order for bipartite lattices.For triangular lattices, however, the exchange coupling J z (t ≈ t ⊥ ) results in Ising-type frustration [30][31][32][33][34] forming novel quantum phases [33,35].In addition, the orbital-exchange interactions in Eq. ( 4) are strongly anisotropic, as a result of the anisotropic porbital hopping, leading to unique properties in the triangular lattice, as shown below.
We numerically solve the frustrated orbital-exchange model by ED with Quspin python package [36,37].Here, we mainly consider lattices with periodic boundary conditions [26].Phase diagrams of the orbitalexchange model are shown in Fig. 3(a),(b).To distinguish different Mott-insulating phases, both fidelity metric g (Fig. S4 [26]) [38][39][40] and static spin structure factor F z (k) are utilized.Within ED, we find that the orbital Skyrmion phase is described by the effective orbitalexchange model with leading-order O(t 2 ,⊥ )/U 2 tunneling processes, as shown in Fig. 3(a).Considering the absence of Dzyaloshinsky-Moriya interactions [19][20][21][22][23] in the effective model, the mechanism for generating orbital Skyrmion texture is a result of the interplay of hoppinginduced anisotropic orbital-exchange interactions and geometric frustration of the triangular lattice.The leadingorder orbital-exchange model, however, only favors two Mott phases, i.e., Skyrmion and H MI [Fig. 3 ter including subleading-order O(t 3 ,⊥ /U 2 ) tunneling processes, ED resolves three Mott phases, i.e., Skyrmion, S MI and H MI [Fig.3(b)], whose conclusion is consistent with the prediction of BDMFT for the extended Bose-Hubbard model [Fig.1(d)].The corresponding static spin structure factors F z (k) for each phase are shown in Fig. 3(c-e) (consistent with the inset of Fig. 2(a), obtained from BDMFT).We remark here that subleadingorder orbital-exchange interactions are absent in the previous studies [24,[27][28][29] and should be included in the effective model to obtain the complete Mott phases.
To obtain more insights of the Skyrmion phases, we investigate ground-state degeneracy, scalar spin chirality χ, and real-space orbital textures.As shown in Fig. 4(a), we find that the Skyrmion phase is actually a gapped phase with a ground-state degeneracy and displays a finite scalar spin chirality order χ = 0 [26], where both the degeneracy and χ disappear in other two Mott phases.We remark here that one only expects an approximate degeneracy in simulations, due to finite-size effects.The ground-state degeneracy indicates the possibility of different types of orbital Skyrmion textures.As shown in Fig. 4(b),(c), two different types of real-space orbital textures are resolved within ED [26].This prediction is consistent with BDMFT results [Fig.1(e)].In addition, we also calculate orbital correlations between different lattice sites for the Skyrmion phase, and observe long-range correlations for the three components of the orbital polarization [26].
Experimental detection.One key feature of the orbital Skyrmion lattice phase is the momentum structure shown in Fig. 2 and 3. Its spin analogue has been revealed by neutron scattering experiments [22].The momentum structure of the orbital Skrymion lattice state can be probed by combining inter-orbital transition tech-niques [3,41] and Bragg spectroscopy [42].With the inter-orbital transition techniques [41], which has been demonstrated in experiments [3], the orbital texture can be converted to density modulations, which maintain the same crystal structure.The periodic density modulations can then be probed by the standard Bragg spectroscopy in cold atom experiments [42].

Conclusion.
We study cold atoms loaded in the porbital band of a triangular optical lattice, and find a chiral orbital Skyrmion lattice phase in a large part of the phase diagram.This quantum state emerges due to natural anisotropic orbital-exchange interaction for p-orbital bosons, unlike the conventional Dzyaloshinskii-Moriya scenario.In this multi-orbital setting, the Skyrmion and anti-Skyrmion lattice states are exactly degenerate due to time-reversal symmetry, in contrast to the widely studied Skyrmion lattice states in spin systems.The exotic orbital polarization texture of the orbital Skyrmion state can be probed by Bragg spectroscopy, a technique accessible to most cold atom experiments.
For the results presented in this work, we consider the system with lattice sites up to N lat = 576 and periodic boundary conditions.The maximum occupation number of the orbital for each normal bath is four to guarantee convergence in our simulations.In the calculations, random initial values are utilized for different lattice sites to break lattice-translational symmetry.We mainly focus on the unit filling case with n px,i + n py,i = 1 in the Mottinsulating regime.The many-body phase diagram is shown in Fig. 1 in the main text.In the BDMFT calculation, we determine the phase boundaries by superfluid order φ x,y , stripe order Θ stripe , and scalar spin chirality χ, as shown in Fig. 2 in main text.Actually, the Mott-insulating phases can also be distinguished by real-space orbital texture S r .As shown in Fig. S1 We also study the robustness of Skyrmion texture against chemical potential.Our calculated filling-dependent phase diagrams are presented in Fig. S2, as a function of chemical potential µ and hopping amplitudes (a) t = t ⊥ , and (b) t = 5t ⊥ , respectively.We observe that there are three many-body quantum phases in the parameter regime studied here, including the Skyrmion, S MI , and S SF phases.The orbital Skyrmion texture is robust and explores a large region of the phase diagrams in the lower hopping regime, indicating large opportunities for experimentally observing the many-body quantum phase.Our studies indicate that the Mott phase with Skyrmion texture is a general long-range order, stabilized in a large parameter regime.

EFFECTIVE ORBITAL-EXCHANGE MODEL
Since the Mott-insulating phase is a state with suppressed number fluctuations, it is convenient to divide the Hilbert space with projection operator P. For the Mott-insulator state with unit filling in the strong coupling limit |t ,⊥ | U , Following Eq. (S29), we can obtain a third-order effective Hamiltonian along the bond direction e 1 where ∆ denotes a triangular lattice with three sites (i, j, k) = (i, i + e 1 , i + e 2 ).The parameters for the three-site , and , respectively, and the parameters for the two-site terms are , and . The final effective Hamiltonian is written as where the J xy term is the z -component of the cross product.We remake here that the cross-product term plays tiny role in our simulations, which does not influence the orbital structures but only shifts the phase boundary slightly.
The orbital Skyrmion texture is actually a result of anisotropic orbital-exchange interactions.

5 FIG. 1 .
FIG. 1. (Color online) (a) The geometry of the twodimensional triangular lattice.(b),(c) The energy dispersion of the lowest p-orbital band of a two-dimensional triangular lattice, with (b) t ⊥ = 0, and (c) t = 0. (d) Many-body phase diagram in terms of hopping amplitudes t and t ⊥ for filling nr = 1, obtained from bosonic dynamical mean-field theory, where the system favors three Mott phases with Skyrmion, stripe (SMI), ferro-orbital (HMI) textures, and two superfluid phases with stripe (SSF) and ferro-orbital (HSF) angular momentum.(e) Real-space orbital textures for Skyrmion (upper) and anti-Skyrmion (down) lattices, where the arrows represent the projection in the xy-plane of orbital polarization texture Sr , and the color denotes the z component.The interactions are U = 3U1 = 3U2.
FIG. 2. (Color online) Phase transitions of ultracold bosonic gases in p-orbital bands of a 2D triangular lattice for different hopping amplitudes (a) t = 5t ⊥ and (b) t ⊥ = 0.055U , obtained via bosonic dynamical mean-field theory [26].Indicated by the dashed lines, the system demonstrates (a) a Mott transition from the Skyrmion to the stripe phase (SMI), followed by the appearance of a stripe superfluid phase (SSF), and (b) a superfluid transition from the ferro-(HSF) to the stripe-orbital phase (SSF) upon increasing hopping t .Inset: Contour plots of spin structure factor F z (k) for Mott phases.The interactions U = 3U1 = 3U2, and filling nr = 1.

FIG. 3 .
FIG. 3. (Color online) Phase diagrams of (a) second-and (b) third-order anisotropic orbital-exchange models as a function of hopping amplitudes t and t ⊥ , obtained from exact diagonalizations for a lattice of N lat = 24 sites [26].The effective model with third-order exchange interactions favors three Mott-insulating phases with Skyrmion, stripe (SMI), and ferro-orbital (HMI) textures, consistent with BDMFT results.(c-e) Contour plots of spin structure factor F z (k) for different Mott phases.The interactions are U = 3U1 = 3U2.

5 FIG. 4 .
FIG. 4. (Color online) (a) Low energy spectra Ei and scalar spin chirality χ as a function of t ⊥ /U , obtained from exact diagonalizations, where the red dots denote the ground-state degeneracy and the green squares are for the excited states.ED calculations of real-space (b) Skyrmion and (c) anti-Skyrmion orbital textures [26], where the arrows represent the projection in the xy-plane of orbital polarization texture Sr , and the color denotes the z component.Here, t /U ≡ 0.01, U = 3U1 = 3U2, and the lattice size N lat = 3 × 6.

5 FIG
FIG. S1. (Color online) Real-space distributions of orbital polarization texture Sr for different Mott phases.The Skyrmion phase is shown in the first column.The SMI and HMI phases are given by the second and third columns, respectively.Here, the length of the arrows represents the amplitude of the local vectors in the xy-plane, and the color denotes the z-component.
FIG. S2. (Color online)Filling-dependent phase diagrams of ultracold bosonic gases in p-orbital bands of a two-dimensional triangular lattice for different hopping amplitudes (a) t = t ⊥ and (b) t = 5t ⊥ , obtained via bosonic dynamical meanfield theory.The system supports two Mott-insulating phases with different types of orbital orders, including Skyrmion and stripe (SMI) orbital textures, and the superfluid phase with stripe orbital angular momentum (SSF).The interactions are U = 3U1 = 3U2.
r,m (J x + J x ) S r • e x θm S r+em • e x θm + J y + J y S r • e y θm S r+em • e y θm + (J z + J z ) S r • e z θm S r+em • e z θm + J xy S r • e x θm S r+em • e y θm − S r • e y θm S r+em • e x θm + r,u,v,w J uvw S r • e u θ1 S r+e1 • e v θ1 S r+e2 • e w θ1 , (S31) FIG. S3. (Color online) Clusters used in the ED calculations.a1 = a/2, √ 3a/2 and a2 = (a, 0) are the primitive vectors of the triangular lattice.The clusters 12, 16, 24, and 24b contain the M points in the reciprocal space, and the lattices 12 and 24b contain the K momentum points.We mainly utilize the cluster 24b in our ED calculation to contain both M and K points.
FIG. S5. (Color online) Real-space orbital correlations of the Skyrmion phase.The red dots indicate a positive correlation sites 0 and j, and the blue dots denote a negative correlation.The lattice structure is 24b, t /U = 0.01, and t ⊥ /U = 0.026.
, and p m,r ≡ (p x,r e x + p y,r e y • e m with e 1 = e y and e 2,3 = − √ 3 2 e x ± 1 2 e y for hopping t ⊥ .