Orbital ordering of ultracold alkaline-earth atoms in optical lattices

We report on a dynamical mean-field theoretical analysis of emerging low-temperature phases in multicomponent gases of fermionic alkaline-earth(-like) atoms in state-dependent optical lattices. Using the example of $^{173}$Yb atoms, we show that a two-orbital mixture with two nuclear spin components is a promising candidate for studies of not only magnetic but also staggered orbital ordering peculiar to certain solid-state materials. We calculate and study the phase diagram of the full Hamiltonian with parameters similar to existing experiments and reveal an antiferroorbital phase. This long-range-ordered phase is inherently stable, and we analyze the change of local and global observables across the corresponding transition lines, paving the way for experimental observations. Furthermore, we suggest a realistic extension of the system to include and probe a Jahn-Teller source field playing one of the key roles in real crystals.


I. INTRODUCTION
In solid-state materials, electrons can occupy different orbital states, which usually determine their directional mobility. Besides spin and charge, this orbital degree of freedom plays an important role in interacting electron systems and lies at the heart of intriguing many-body phenomena such as colossal magnetoresistance, heavy fermions, and the Kondo effect [1][2][3].
In particular, orbital ordering is one of the key phenomena in materials with multiorbital structure such as transitionmetal oxides. Similar to the ordered pattern of spins in the ground state of an antiferromagnet, electrons from different -orbital manifolds can spatially arrange in these materials [4]. While great advances have been made in both the experimental observation and the theoretical description of orbital ordering [5][6][7][8], the microscopic origin of the processes is still not fully understood. Numerical calculations for these systems are challenging due to the simultaneous presence of electron-electron as well as electron-phonon interactions, which both can effect orbital ordering [9,10]. Therefore, quantum simulations of the corresponding model Hamiltonians could shed light onto the competing mechanisms and the nature of orbitally-ordered phases.
Ultracold atoms have become a versatile and clean platform for quantum simulations of solid-state systems in the last decade [11]. In particular, ultracold fermionic atoms in optical lattices have allowed to study the single-orbital Fermi-Hubbard model [12][13][14][15], which is believed to describe certain high-temperature superconductors [16][17][18]. More recently, a two-orbital Fermi-Hubbard system has been realized with ultracold alkaline-earth(-like) atoms (AEAs) in a statedependent optical lattice (SDL) [19]. For atoms of this kind, the availability of the long-lived metastable 3 P 0 electronic * a sotnikov@kipt.kharkov.ua state (denoted as ) in addition to the 1 S 0 ground state (denoted as ) allows populating two orbital states of the lattice with distinct kinetic and interaction properties, which makes these atoms attractive candidates for the study of orbital phenomena [20][21][22].
In this work, we report on the possibility to approach and simulate orbital ordering with AEAs in SDLs. Our study is based on dynamical mean-field theory (DMFT) applied to a two-orbital Fermi-Hubbard model with parameters closely related to existing experimental implementations with 173 Yb atoms [19]. Nevertheless, our results are also applicable to other fermionic AEAs due to the similarity of relevant interaction properties [20]. We calculate the phase diagram for realistic orbital fillings and, in addition to multiple magneticallyordered phases, we find a particular stable orbitally-ordered phase. The ordering instability in this system is driven by the different intra-and interorbital on-site interaction of atoms in the and state, corresponding to electron-electron interactions in a solid state system. Transitions to this long-rangeordered phase result in noticeable changes of experimentallyaccessible observables, which we determine for the fraction of doubly-occupied lattice sites, the orbital density distribution in a harmonic trap, and nearest-neighbor correlations. We also show how the influence of electron-phonon interactions in the form of the Jahn-Teller effect (JTE) can be probed with a suitable superlattice potential or by adjusting the density of atoms appropriately.

II. SYSTEM, MODEL, AND METHOD
We consider a two-orbital mixture of AEAs in the lowestenergy band of a square optical lattice [see Fig. 1(a)], which is state-dependent, i.e., the lattice depth differs for the two orbital states and . In addition to the orbital degree of freedom, atoms can occupy one of two nuclear spin states (denoted by ↓ and ↑), which are equally populated. Within the tight-binding approximation, this system can be described by (1) illustrating the interaction parameters. The lattice bonds are shown as gray lines and the blue (yellow) circles refer to ( ) atoms in different spin states as indicated by the arrows. (b) Illustration of the state-dependent lattice potential (solid lines) which introduces distinct hopping amplitudes > for and atoms (blue and yellow circles). the following two-orbital Hubbard model [20,21]: where ⟨ , ⟩ denotes the set of nearest-neighbor lattice sites, † ( ) is the fermionic creation (annihilation) operator of an atom in the orbital ∈ { , } with spin ∈ {↑, ↓}, and = † is the corresponding density operator. Among the relevant Hubbard parameters are the orbital hopping amplitudes , the intraorbital interactions as well as the interorbital direct interaction and the exchange interaction ex . We obtain these parameters from the experimentally determined -wave scattering lengths and a bandstructure calculation. In our theoretical approach, the chemical potentials and determine the densities and of atoms in the corresponding orbital state. We restrict our study to the regimes of low lattice fillings, = ( + ) ≤ 2, to avoid three-body losses [23].
In the following, we focus on parameters for a realistic 173 Yb system with a state-dependent optical lattice, which localizes atoms, i.e., > as depicted in Fig. 1(b). This ensures lossy collisions between atoms are strongly suppressed [19,20,24]. We consider a combination of a square state-dependent lattice and a sufficiently deep perpendicular state-independent lattice (see Appendix A), which ensures the system is quasi-two-dimensional (quasi-2D) and can be described by the Hamiltonian in Eq. (1). We note that our theoretical calculation can be directly extended to a threedimensional system [25], but the 2D geometry has advantages for experimental realizations, especially for quantum gas microscopy [11]. In general, and atoms experience different lattice potentials due to distinct polarizabilities ( ) ≠ ( ) of the respective states at a given wavelength . At fixed depth of the SDL, the ratio of orbital mobility, ∕ , can be tuned by adjusting this wavelength appropriately. The associated polarizability ratio = ( )∕ ( ) increases monotonously between certain atomic transition wavelengths, in particular, from 1 at the so-called magic wavelength ( = 759 nm) to 3.3 at = 670 nm [19,26]. We consider this experimentally accessible regime of to tune the orbital-dependent mobility, which enhances or suppresses ordered phases.
The on-site interaction energies of 173 Yb atoms in the proposed setup demonstrate the hierarchy < ∼ ex < due to the relatively large scattering length of the orbitally-symmetric state [27], which contributes both to and ex [20,28]. This contrasts with the well-known Slater-Kanamori parametrization of the Coulomb interaction in correlated electron systems [29], where the largest quantity is the Coulomb parameter = = , and the Hund's coupling is bounded from above, ex ≤ ∕3, which ensures repulsive interactions for all spin and orbital components. In the case of cubic symmetry in the orbital space, the direct interaction amplitude is given as = − 2 ex , such that a different hierarchy is observed, ex ≤ < . Nevertheless, the relatively small and amplitudes do not have any strong implications on magnetic phases that can be approached with cold 173 Yb atoms [30]. Moreover, as we show below, the difference between and ( − ex ) as well as , makes the system more susceptible to an orbital ordering instability, which also appears in transition metal compounds [4].
For studying low-temperature phases, we employ DMFT, which is approximative and only exact in the limit of an infinite-dimensional system [31]. Nevertheless, it has become a well-accepted method successfully applied to stronglycorrelated electron systems and has also found applications in the description of ultracold atoms in optical lattices [32][33][34]. The calculations are performed with an exactdiagonalization impurity solver [35], which preserves the SU(2) spin-rotational symmetry of the two-orbital Hubbard model in Eq. (1) [36]. To measure observables in the symmetry-broken phases, we perform doubling of the unit cell, i.e., allow two different DMFT solutions on neighboring lattice sites.

III. RESULTS
For the central case in our DMFT analysis, we choose a fixed set of Hubbard parameters, = 0.26 , = 6.8 , = 17 , = 32 , and ex = 23 , corresponding to a typical quasi-2D SDL with the polarizability ratio = 2.1 (see Appendix B). The exact choice of parameters is not crucial but motivated by the experimental accessibility and the signatures of orbital ordering considered in our study. We briefly note that the interaction energies and ex are comparable to the band gap of the SDL and need to be renormalized accordingly. Renormalizing these parameters is particularly non-trivial in our regime of quasi-2D geometry and mixed confinement. Instead, we perform the renormalization on the basis of an approximative scheme and verify independently that our main results are not sensitive to the precise magnitude of and ex (see Appendix B). The low-temperature phase diagram derived from the DMFT calculation is shown in Fig. 2(a). For certain orbital fillings ( , ) and at sufficiently low temperature ∕ ≤ 1 ( = 1 is used below), we distinguish ferromagnetic (FM), antiferromagnetic (AFM), and antiferroorbital (AFO) longrange-ordered phases as illustrated in Fig. 2 Since magnetically-ordered phases are not the main focus of the current study, we only briefly comment on the key observations. According to the diagram shown in Fig. 2(a), the -AFM and AFM instabilities appear along diagonal lines in the -plane, where the total density is ≈ 1 or ≈ 2, respectively. These are the Mott-insulating regimes with one exception at = 2, where the system becomes an insulator with vanishing local magnetic moments. At ≈ 1 and variable , the Hamiltonian in Eq. (1) can be mapped to the double-exchange model (FM Kondo-lattice model) with both FM and AFM terms [37]. We note that away from the polarizability ratio = 1, FM ordering becomes stabilized with increasing due to the larger exchange interaction ex and the stronger localization of atoms. In contrast, the AFM phases involving the orbital become suppressed due to the strong localization of atoms and an increase of the local interaction amplitudes.
The AFO phase is characterized by the alternating occupation of neighboring lattice sites with atoms in different orbitals. In analogy to the Néel order of spins, the lattice can be viewed as a set of two sublattices: one is dominantly occupied by pairs of atoms, while the other is dominantly occupied by single atoms in the state. This configuration is similar to orbital ordering in solid-state materials, where sublattices are formed by electrons occupying different orbital angularmomentum states of the lattice ions [4]. In contrast to real crystals, our proposed implementation does not introduce directional or interorbital hopping, which makes it more feasible for direct experimental realizations.
A peculiar feature of the system under study is that the AFO phase is stabilized in a wide region around = 1 and = 0.5 as can be seen in Fig. 2(a). In most regions, this phase is accompanied by charge order, i.e., the periodic modulation of the total density . The transition to the longrange-ordered state is purely driven by an interplay between the direct and superexchange interaction amplitudes. This can be intuitively understood in the strong coupling limit, In this case, the dominant superexchange amplitude ∼ 2 ∕ ef f reduces the total energy in an arrangement of pairs of atoms next to atoms on neighboring lattice sites, which yields the antiferrorbital order illustrated in Fig. 2(b). We note that the system is actually in a slightly different regime with intermediate coupling ( ≲ ef f ), for which there is no exact analytical formula of the corresponding amplitude, but the intuitive picture remains valid. Furthermore, it is worth mentioning that the AFO phase is not limited to the chosen 2D geometry. Supporting calculations performed for both a three-dimensional [38] and a onedimensional system (matrix product state algorithm [39]), reveal qualitatively similar correlations in the orbital domain at comparable densities.
For our chosen set of Hubbard parameters, the AFO phase is enhanced compared to most magnetically-ordered phases and the DMFT analysis yields a transition point at the critical temperature ∕ = 0.31. We also obtain characteristic values for the entropy per particle required for orbital ordering, which are calculated from the Maxwell relation for the local density of atoms (see Appendix C). We estimate the maximal entropy per particle in the bulk for the AFO phase to be ≈ 0.8, which is related to approximately a tenth of the Fermi temperature in a harmonic trap under the assumption that loading into the optical lattice is adiabatic [40][41][42]. Degenerate Fermi gases of 173 Yb atoms in the state have been reported in a similar entropy and temperature regime [23], but the preparation of two-orbital mixtures is more challenging, since atoms need to be partly excited into the state for the desired orbital population (see Appendix A).
In Fig. 3, we explore the influence of the polarizability ratio on the critical temperature and entropy of the AFO phase. While the orbitally-ordered phase vanishes completely in the vicinity of = 1, it is stabilized with increasing and the critical temperature only changes negligibly for > 2. Similarly, the critical entropy increases at small , reaches its maximal value at moderate coupling ( ≈ 2), and then slowly decreases due to a stronger suppression of particle-number fluctuations at larger values of the interaction strengths. This critical behavior is analogous to the one observed in the proximity of AFM phases [36] and motivates the choice of = 2.1 for our study.
In particular, the AFO phase can be probed experimen- tally since it covers a sizeable fraction of the -phase diagram, as can be seen in Fig. 2(a). Due to its relative stability against particle-number fluctuations (thermally-induced metal-insulator crossover region at = 1.5) and almost equidistant separation from other insulating regimes ( = 1, = 2, and = 1), it should only require relatively coarse tuning of the respective densities.
First, we analyze how the AFO ordering can be detected by measuring the fraction of lattice sites occupied by pairs of atoms,  = ⟨ ↑ ↓ ⟩. This observable can be probed experimentally by measuring the atom number upon removal of atoms on doubly-occupied lattice sites with a resonant photoassociation pulse, a well-established measurement technique, which has been successfully applied to ultracold 173 Yb atoms in optical lattices [43]. As shown in Fig. 4(a), we first keep the atomic densities fixed, = 1.0 and = 0.5, and vary the temperature. In addition, we also study the  dependencies at fixed temperature but variable or to quantify the sensitivity on the density in each orbital [see Fig. 4(b)-(c)]. The temperature and density dependencies of  clearly indicate the enhancement of local pairing of atoms in the AFO phase. Below the critical temperature, this effect increases but approaches a saturated regime for temperatures ∕ ≤ 0.2. The site-averaged  signal in the AFO phase at ∕ = 0.2 differs from the one obtained in an artificially-restricted normal phase at the same temperature by ≈ 8%. Additional calculations with variable Hubbard parameters show that this value can be slightly increased by a reduction of the lattice depth in any of the three spatial directions.
In comparison to the global observable  , the AFO phase shows much stronger signatures in local quantities such as the in-trap density distribution, which can be measured directly with high-resolution in-situ imaging. Ultracold atoms trapped in an attractive optical lattice potential usually experience harmonic confinement due to the curvature of the Gaussian laser beams. The resulting smooth change of the chemical potential can then lead to the coexistence of multiple phases in a single trap, such as the well-known shell structure consisting of spatially-alternating Mott-insulating and metallic regions. We explore the density profiles in a harmonic trap across the thermally-induced AFO phase transition in Fig. 5(a)-(c). Below the critical temperature, in parallel with the development of AFO correlations across the trap, we observe the formation of a Mott-insulating plateau at = 1.5, which is clearly visible in Fig. 5(a). Interestingly, the phase separation of and atoms, as well as another Mott-insulating plateau at = 1, can already be observed above the critical temperature due to a decoupling from the superexchange energy scales [see Fig. 5(c)]. This could allow detecting this signature as a precursor of the AFO phase in an experiment, even above the actual transition point ( ∕ ≲ 0.8).
Next, we focus on density correlations between individual lattice sites, which could be directly probed with single-site resolved imaging of and atoms [11]. In Fig. 5(d)-(e), we show temperature dependencies of the densities and on two neighboring lattice sites across the AFO phase transition. In general, these show strong signals from pair formation and redistribution of atoms in different orbitals on a checkerboard-like pattern in the lattice [see Fig. 2(b)]. Already slightly below , the density ( ) reaches 1.5 (0.25) on the first site and 0.5 (0.75) on the neighboring site. We also expect the build-up of spatial correlations beyond nearestneighbor sites, whose amplitudes cannot be accurately calculated within DMFT but could be probed in the experiment. With recent advances in quantum gas microscopy of AEAs [44,45], local density correlations could provide a direct detection method of the AFO phase and its properties. Moreover, local control in these experiments could allow to precisely engineer and study excitations in the AFO regime [4].
Finally, let us briefly comment on a potential extension of our proposal to host not only the superexchange-driven mechanism for orbital ordering but also include a source field comparable to the Jahn-Teller effect in transition metal oxides [4]. For the analog of the antiferro-distortive JTE (as in K 2 CuF 4 ) [46], we consider adding a superlattice structure for the orbital. For the proposed 173 Yb system, an additional optical lattice with the wavelength ≈ 1380 nm and | | ≫ 1 would produce a suitable superlattice potential, which acts predominantly on atoms and lowers their energy at every second site of the original SDL. The DMFT analysis for the corresponding staggered potential confirms that the AFO phase can be substantially extended to higher temperatures with a transformation of the second-order transition point to a crossover regime due to explicit symmetry breaking by the superlattice. In contrast, the analogue of the ferro-distortive JTE (as in La 2 CuO 4 ) [46] could be probed without additional potentials. Its destructive impact on staggered orbital ordering can be analyzed by varying / along a line of constant total density, in particular, = 1.5, shown as dotted line in Fig. 2(a).

IV. SUMMARY AND OUTLOOK
We show that AEAs in SDLs are promising candidates for the experimental observations of orbital ordering phenomena and potentially could improve the understanding of related mechanisms in solid-state materials. In particular, by means of changing the lattice depth and polarizability ratio between different orbital states, a capability to enhance or suppress the superexchange contributions to the AFO ordering instability is demonstrated. At the same time, in a well-controlled and independent manner, contributions analogous to the JTE in crystals could be explored by adjusting the orbital densities or by introducing a superlattice potential. The rich structure of the phase diagram revealed in this study also makes AEAs in SDLs suitable for studies of open questions on the critical behavior and excitations in transition-metal oxides hosting orbitally-ordered as well as various magnetic and superconducting phases [3][4][5][6]47].
Our analysis oriented towards experimental implementations with 173 Yb atoms reveals that the SDL substantially increases the difference between the intraorbital interactions, ( − ) ≳ . Therefore, the AFO instability crucially depends on the energy gap to the closest interorbital excitation [( − ex − ) for ex > 0]. This small gap gives the largest contribution to the corresponding AFO superexchange amplitude, which depends less on and the energy of the other interorbital excitation [( + ex − ) for ex > 0]. Therefore, similar calculations and experiments could be realized with related species, such as 87 Sr or 171 Yb. While the former and 173 Yb have comparable ordering of the interaction parameters [48], the latter features antiferromagnetic exchange interaction ex < 0 [49] and almost vanishing | | ≪ [50], which could provide an interesting extension of the phase diagram discussed in our study.
At higher spin symmetry, the AFO phases may demonstrate unconventional space modulations involving more than two sublattices. These could naturally be studied with 173 Yb when the large SU( ≤ 6) symmetry in the and orbital is utilized [20]. Another related effect concerns the potential magnetic order of the SU(2)-symmetric mixture in the AFO phase at very low temperatures, which requires a comprehensive analysis of potential sublattice structures and remains an interesting task for future theoretical research. density in a 2D SDL at low enough temperatures. We first focus on the optical lattice implementation and briefly outline possible state preparation techniques.
For 173 Yb, a monochromatic SDL at a wavelength of 670 nm (polarizability ratio = 3.3) has been implemented in one dimension [19] and can be realized similarly in 2D. For our choice of = 2.1, theoretical calculations of the polarizability [26] yield a wavelength of 690 nm, which is accessible with commercial laser systems. We note that the precise value of this wavelength has only negligible influence on the results discussed in the main text. For the SDL, we consider a fixed lattice depth of , = 5 SDL rec ( atoms) to ensure strong suppression of next-nearest-neighbor tunneling and the validity of the tight-binding approximation. For the strong confinement along , we consider a deep magic-wavelength ( = 759 nm) lattice, = 18 rec , such that the system is in the quasi-2D regime. Here, SDL rec = ℎ × 2.4 kHz and rec = ℎ × 2.0 kHz refer to the recoil energy from a photon of the SDL or magicwavelength lattice, respectively.
The two-orbital mixture can be prepared in the optical lattice by optically exciting part of the atoms with an appropriate laser pulse [19]. Besides the orbital degree of freedom, 173 Yb atoms feature six nuclear spin states in and , with ∈ {−5∕2, −3∕2, … , +5∕2}. Due to SU( )-symmetric collisions, a stable subset of these states can be prepared and used in the experiment [20,28]. For the realization of the Hamiltonian in Eq. (1), we only consider two spin states, = −5∕2 and +5∕2 (denoted by ↓ and ↑), as discussed in the main text. For the state preparation, we suggest utilizing two additional spin states, = −3∕2 and +3∕2 (denoted by ↙ and ↗). Optical pumping on the intercombination line allows to prepare atoms in an imbalanced mixture of all four spin states such that ( ↙ + ↗ )∕( ↓ + ↑ ) equals the desired ratio of ∕ . Subsequent transfer of the ancillary states (↙ and ↗) into the orbital with circularly-polarized light yields the desired densities and with spin states ↓ and ↑ (see Ref. [51] for a similar technique).

Appendix B: Hubbard parameters
We calculate the Hubbard parameters from the numerical solution of the band structure of a separable three-dimensional optical lattice in the tight-binding approximation and with the corresponding lattice depths discussed in Appendix A. Since and atoms experience different lattice depths in the ( , ) plane and along , we use independent band structures for each orbital and spatial direction. We list all relevant parameters in Table I for a range of polarizability ratios considered in the main text.
The on-site interaction strength ′ is typically calculated from the corresponding s-wave scattering length ′ , Here, is the atomic mass and ( ) is the Wannier function of the corresponding orbital ∈ { , } derived from the band-structure calculation. In the limit of large scattering lengths comparable to the lattice spacing, ∼ lat , contributions from higher bands of the optical lattice become sizeable. Nevertheless, such a system can still be described within the lowest-band approximation by absorbing these contributions into renormalized Hubbard parameters [52,53]. In our case, lat ≈ 6 × 10 3 0 is the (smallest) lattice constant with 0 the Bohr radius. For the intraorbital scattering lengths = 199 0 [50], = 306 0 [28], and the interorbital singlet scattering length − = 220 0 [27,28], the corrections are small and neglected. However, the large orbitally-symmetric scattering length + ≈ 2 × 10 3 0 [27] leads to a significant correction of the corresponding amplitude + , which would otherwise exceed the band gap.
The system discussed in the main text features anisotropic and mixed confinement due to the SDL and the quasi-2D geometry, which prevents us from directly applying existing results for the renormalization of + [52]. Instead, we use the geometric mean of both orbitals as the effective lattice depth and approximate each lattice site with a harmonic oscillator potential [19,54]. In addition, we apply first-order perturbation theory to account for the anharmonic cosine potential of the optical lattice. Finally, we assume spatial separability of the problem and calculate two independent solutions for the ( , ) as well direction, which we combine into the single interaction amplitude = 2∕3 , 1∕3 . When applied to an isotropic system with comparable lattice depths, our results reasonably agree with Ref. [52]. We find an on-site interaction energy + in excess of the band gap of atoms along and by up to 60% ( = 3.3), which suggests that our approximate approach fails to correctly predict the renormalized Hubbard parameter. Although the effective + in the experiment will be different, we verify that the phases discussed in the main text are robust against variation of this parameter on a similar scale.
The large scattering length + also causes an increased relevance of non-Hubbard terms in the Hamiltonian, specifically, direct off-site interactions and density-assisted tunneling [53]. While we expect the former to be negligible in our regime, the latter could become comparable to for the orbitallysymmetric interaction channel. We cannot directly incorporate this term into our DMFT calculation, but the main effect will be a renormalization of the hopping amplitudes for sites occupied simultaneously by and atoms. In principle, these excitations should mainly occur virtually in the AFO phase at ≤ 1.0 and ≤ 0.5. At higher densities, we assume the effects can be absorbed into a modified + , respectively and ex , which again should not alter the phase diagram significantly.

Appendix C: DMFT calculation
In the DMFT analysis, we employ an exact diagonalization solver for the Anderson impurity problem with up to four bath orbitals per each spin and orbital component. The DMFT self-consistency conditions for two sublattices are applied in the analysis of the AFO and AFM phases, while the normal and FM phases are analyzed within the single-site lattice projection [30].
We obtain the inhomogeneous distributions in the harmonic trap and entropy dependencies within the local density approximation. The entropy is calculated by numerical integration of the Maxwell relation, = ∫ d ( ∕ ) on the interval from the vacuum state, ( 0 , 0 ) = 0, to the chemical potential values and , which yield the desired densities of atoms, = 1 and = 0.5, in particular.
For the phase diagram in Fig. 2, we fit the parameters ,