Physics of higher orbital bands in optical lattices: a review

In the tight binding regime, an optical lattice can be treated as individual harmonic oscillators, which are coupled by quantum tunnelings. On each harmonic oscillator centered at a lattice site labeled by its position $\mathbf{R}$, we have discrete energy levels with orbital wavefunctions ${{\phi}_{\alpha}}\left(\mathbf{x}-\mathbf{R}\right)$. Associated with the localized orbital wavefunctions, we can define the lattice operators ${{b}_{\alpha}}\left(\mathbf{R}\right)$. To do this, it has to be enforced that the orbital wavefunctions are orthonormal. The simple eigen wavefunctions of harmonic oscillators do not satisfy orthonormal condition, for the reason that there are overlaps between orbital wavefunctions on neighboring sites.

The procedure to construct the orthogonal basis from the localized harmonic oscillator wavefunctions is the following. We start with the harmonic oscillator wavefunctions ${{\phi}_{\alpha}}\left(\mathbf{x}-\mathbf{R}\right)$ localized on site $\mathbf{R}$. These wavefunctions are already approximately orthogonal, i.e.

$$\begin{eqnarray}&&{\int}^{}\text{d}\mathbf{x}\phi _{\alpha}^{\ast}\left(\mathbf{x}-\mathbf{R}\right){{\phi}_{{{\alpha}^{\prime}}}}\,\left(\mathbf{x}-{{\mathbf{R}}^{\prime}}\,\right)={{\delta}_{\alpha {{\alpha}^{\prime}}}}{{\delta}_{\mathbf{R}{{\mathbf{R}}^{\prime}}}}+{{\epsilon}_{\alpha \mathbf{R},{{\alpha}^{\prime}}{{\mathbf{R}}^{\prime}}}},\end{eqnarray}$$

where ${{\epsilon}_{\alpha \mathbf{R},{{\alpha}^{\prime}}{{\mathbf{R}}^{\prime}}}}$ are small numbers. By definition we know that $[\epsilon ]$ is a traceless Hermitian matrix. Then we improve this basis by introducing

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\tilde{{\phi}}\,}_{\alpha}}\left(\mathbf{x}\right)={{\phi}_{\alpha}}\left(\mathbf{x}\right)-\frac{1}{2}\underset{{{\alpha}^{\prime}}{{\mathbf{R}}^{\prime}}}{\sum}\,{{\epsilon}_{{{\alpha}^{\prime}}{{\mathbf{R}}^{\prime}},\alpha \mathbf{R}}}{{\phi}_{{{\alpha}^{\prime}}}}\left(\mathbf{x}-{{\mathbf{R}}^{\prime}}\right) \end{array}.\end{eqnarray} \tag{ 1 }$$

After that ${{\phi}_{\alpha}}\left(\mathbf{x}-\mathbf{R}\right)$ is renormalized as

$$\begin{eqnarray}&&{{\tilde{\phi}}_{\alpha}}\left(\mathbf{x}\right)\to {{\tilde{\phi}}_{\alpha}}\left(\mathbf{x}\right)/\sqrt{{\int}^{}\text{d}{{\mathbf{x}}^{\prime}}|{{{\tilde{\phi}}}_{\alpha}}\left({{\mathbf{x}}^{\prime}}\right){{|}^{2}}}.\end{eqnarray}$$

The improved wavefunctions satisfy a better approximate orthogonal condition

$$\begin{eqnarray}&&{\int}^{}\text{d}\mathbf{x}\tilde{{\phi}}\,_{\alpha}^{\ast}\,\left(\mathbf{x}-\mathbf{R}\right){{\tilde{{\phi}}\,}_{{{\alpha}^{\prime}}}}\,\left(\mathbf{x}-{{\mathbf{R}}^{\prime}}\,\right)={{\delta}_{\alpha {{\alpha}^{\prime}}}}{{\delta}_{\mathbf{R}{{\mathbf{R}}^{\prime}}}}+\mathcal{O}\left({{\epsilon}^{2}}\right).\end{eqnarray}$$

The above procedure can be iterated N times to get the orthonormal basis to the precision of $\mathcal{O}\left({{\epsilon}^{{{2}^{N}}}}\right)$.

Once we have the orthonormal basis, the tunnelings between $\mathbf{R}$ and ${{\mathbf{R}}^{\prime}}$ are calculated as

$$\begin{eqnarray}&&{{t}_{\alpha {{\alpha}^{\prime}}}}\left(\mathbf{R}-{{\mathbf{R}}^{\prime}}\,\right)={\int}^{}\text{d}\mathbf{x}\tilde{{\phi}}\,_{\alpha}^{\ast}\left(\mathbf{x}-\mathbf{R}\right)H\left(\mathbf{x}\right){{\tilde{{\phi}}\,}_{{{\alpha}^{\prime}}}}\,\left(\mathbf{x}-{{\mathbf{R}}^{\prime}}\,\right),\end{eqnarray} \tag{ 2 }$$

where $H\left(\mathbf{x}\right)$ is the Hamiltonian in the first quantization form $H=-\frac{{{\hslash}^{2}}}{2m}{{\vec{\nabla}}^{2}}+V\left(\mathbf{x}\right).$ The lattice model Hamiltonian including tunnelings is given by

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \hat{H}=\underset{\alpha {{\alpha}^{\prime}},\mathbf{R}{{\mathbf{R}}^{\prime}}}{\sum}\,{{t}_{\alpha {{\alpha}^{\prime}}}}\left(\mathbf{R}-{{\mathbf{R}}^{\prime}}\right)b_{\alpha}^{\dagger}\left(\mathbf{R}\right){{b}_{{{\alpha}^{\prime}}}}\left({{\mathbf{R}}^{\prime}}\right). \end{array}\end{eqnarray} \tag{ 3 }$$

Without truncating the basis, the Hamiltonian is exact, from which the band structure can be calculated. If we only keep the lowest harmonic wave functions, this lattice Hamiltonian gives qualitatively correct band structures for deep lattices.

The procedure described above to construct orthogonal orbital wave functions is one essential step if one uses harmonic approximation. In principle, the constructed wave functions are not the same as the maximally localized Wannier functions (Kivelson 1982, Marzari et al 2012, Uehlinger et al 2013, Ganczarek et al 2014). The procedure to calculate suchmaximally localized Wannier functions is not as straightforward and is beyond the scope of this review.

2.1.1. Multi-band Hubbard model.

Considering interacting bosonic atoms loaded on excited bands, the physics will be described by a multi-band Hubbard model

$$\begin{eqnarray}\begin{array}{*{35}{l}} H & = & \underset{\mathbf{R}{{\mathbf{R}}^{\prime}}}{\sum}\,{{t}^{(\alpha )}}\left(\mathbf{R}-{{\mathbf{R}}^{\prime}}\,\right)b_{\alpha}^{\dagger}\left(\mathbf{R}\right){{b}_{\alpha}}\left({{\mathbf{R}}^{\prime}}\,\right) \\ {} & {} & +\underset{\mathbf{R}}{\sum}\,{{U}_{{{\alpha}_{1}}{{\alpha}_{2}}{{\alpha}_{3}}{{\alpha}_{4}}}}b_{{{\alpha}_{1}}}^{\dagger}\left(\mathbf{R}\right)b_{{{\alpha}_{2}}}^{\dagger}\left(\mathbf{R}\right){{b}_{{{\alpha}_{3}}}}\left(\mathbf{R}\right){{b}_{{{\alpha}_{4}}}}\left(\mathbf{R}\right). \end{array}\end{eqnarray} \tag{ 4 }$$

With weak interaction, the coupling constants ${{U}_{{{\alpha}_{1}}{{\alpha}_{2}}{{\alpha}_{3}}{{\alpha}_{4}}}}$ can be estimated at tree-level as (Jaksch et al 1998, Liu and Wu 2006, Zhou et al 2011b, Li et al 2016, Dutta et al 2015a)

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{U}_{{{\alpha}_{1}}{{\alpha}_{2}}{{\alpha}_{3}}{{\alpha}_{4}}}}=\frac{4\pi {{a}_{s}}{{\hslash}^{2}}}{2m}{\int}^{}{{\text{d}}^{3}}\mathbf{x}\phi _{{{\alpha}_{1}}}^{\ast}\left(\mathbf{x}\right)\phi _{{{\alpha}_{2}}}^{\ast}\left(\mathbf{x}\right){{\phi}_{{{\alpha}_{3}}}}\left(\mathbf{x}\right){{\phi}_{{{\alpha}_{4}}}}\left(\mathbf{x}\right), \end{array}\end{eqnarray} \tag{ 5 }$$

where a_s is the s-wave scattering length, tunable with Feshbach resonance techniques. With fermionic atoms, we have a similar Hubbard model with interactions between hyperfine states.

In terms of field operators, the Hamiltonian of particles moving in optical lattices is

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \hat{H}={\int}^{}{{\text{d}}^{d}}\mathbf{x}{{\psi}^{\dagger}}\left(\mathbf{x}\right)\left(-\frac{{{\hslash}^{2}}}{2m}{{{\vec{\nabla}}}^{2}}+V\left(\mathbf{x}\right)\right)\psi \left(\mathbf{x}\right), \end{array}\end{eqnarray} \tag{ 6 }$$

where $\psi \left(\mathbf{x}\right)$ is a field operator. It can be either bosonic or fermionic. Statistics is irrelevant here to determine single-particle band structures. We expand the operator $\psi \left(\mathbf{x}\right)$ in the momentum basis

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \psi \left(\mathbf{x}\right)=\underset{\mathbf{K}}{\sum}\,\frac{1}{\sqrt{{{N}_{s}}}}\underset{\mathbf{k}}{\sum}\,{{{\tilde{a}}}_{\mathbf{K}}}\left(\mathbf{k}\right){{\text{e}}^{\text{i}\left(\mathbf{K}+\mathbf{k}\right)\centerdot \mathbf{x}}}, \end{array}\end{eqnarray} \tag{ 7 }$$

where $\mathbf{K}$ labels the reciprocal lattice vectors, $\mathbf{k}$ the lattice momentum, and N_s the number of lattice sites. Here and henceforth, the lattice constant is set as the length unit. Optical lattice potentials $V\left(\mathbf{x}\right)$, unlike the potentials in electronic materials, can typically be written as superpositions of just a few plane waves, i.e.

$$\begin{eqnarray}&&V\left(\mathbf{x}\right)=\underset{\mathbf{K}}{\sum}\,v\left(\mathbf{K}\right){{\text{e}}^{\text{i}\mathbf{K}\centerdot \mathbf{x}}}.\end{eqnarray}$$

For example, the potential of a square lattice created by laser is

$$\begin{eqnarray}\begin{array}{*{35}{l}} V\left(\mathbf{x}\right) & = & -{{V}_{0}}\left[{{\sin}^{2}}(kx)+{{\sin}^{2}}(ky)\right] \\ {} & = & -\frac{{{V}_{0}}}{4}\left[{{\text{e}}^{2\text{i}kx}}+{{\text{e}}^{2\text{i}ky}}+\text{c}.\text{c}.\right]+\text{const}., \end{array}\end{eqnarray}$$

where k is the wavevector of the laser beams. The Hamiltonian in momentum space reads as

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \hat{H}=\underset{\mathbf{k}}{\sum}\,\underset{{{\mathbf{K}}_{1}},{{\mathbf{K}}_{2}}}{\sum}\,{{\mathcal{H}}_{\mathbf{k}}}\left({{\mathbf{K}}_{1}},{{\mathbf{K}}_{2}}\right)\tilde{a}_{{{\mathbf{K}}_{1}}}^{\dagger}\left(\mathbf{k}\right){{{\tilde{a}}}_{{{\mathbf{K}}_{2}}}}\left(\mathbf{k}\right), \end{array}\end{eqnarray} \tag{ 8 }$$

with the matrix given by

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\mathcal{H}}_{\mathbf{k}}}\left({{\mathbf{K}}_{1}},{{\mathbf{K}}_{2}}\right)=\frac{{{\hslash}^{2}}{{\left({{\mathbf{K}}_{1}}+\mathbf{k}\right)}^{2}}}{2m}{{\delta}_{{{\mathbf{K}}_{1}}{{\mathbf{K}}_{2}}}}+v\left({{\mathbf{K}}_{1}}-{{\mathbf{K}}_{2}}\right). \end{array}\end{eqnarray} \tag{ 9 }$$

Diagonalizing this matrix, we get the band structure ${{E}_{n}}\left(\mathbf{k}\right)$ and the eigenvectors $\lambda _{\mathbf{K}}^{(n)}\left(\mathbf{k}\right)$, with n the band index. The Hamiltonian in the eigen-basis reads

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \hat{H}=\underset{n}{\sum}\,\underset{\mathbf{k}}{\sum}\,{{E}_{n}}\left(\mathbf{k}\right)b_{n}^{\dagger}\left(\mathbf{k}\right){{b}_{n}}\left(\mathbf{k}\right), \end{array}\end{eqnarray} \tag{ 10 }$$

with ${{b}_{n}}\left(\mathbf{k}\right)={\sum}_{\mathbf{K}}\lambda _{\mathbf{K}}^{(n)\ast}\left(\mathbf{k}\right){{\tilde{a}}_{\mathbf{K}}}\left(\mathbf{k}\right)$.

The Wannier basis is given by

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{b}_{n}}\left(\mathbf{R}\right)=\frac{1}{\sqrt{{{N}_{s}}}}\underset{\mathbf{k}}{\sum}\,{{b}_{n}}\left(\mathbf{k}\right){{\text{e}}^{\text{i}\mathbf{k}\centerdot \mathbf{R}}}. \end{array}\end{eqnarray} \tag{ 11 }$$

Inversely we have ${{b}_{n}}\left(\mathbf{k}\right)=\frac{1}{\sqrt{{{N}_{s}}}}{\sum}_{\mathbf{R}}{{b}_{n}}\left(\mathbf{R}\right){{\text{e}}^{-\text{i}\mathbf{k}\centerdot \mathbf{R}}}.$

The Wannier wavefunctions of the Bloch bands are given by

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{w}_{n}}\left(\mathbf{x}-\mathbf{R}\right)=\underset{\mathbf{K}}{\sum}\,\left[\frac{1}{{{N}_{s}}}\underset{\mathbf{k}}{\sum}\,\lambda _{\mathbf{K}}^{(n)}\left(\mathbf{k}\right){{\text{e}}^{\text{i}\left(\mathbf{K}+\mathbf{k}\right)\centerdot \left(\mathbf{x}-\mathbf{R}\right)}}\right]. \end{array}\end{eqnarray} \tag{ 12 }$$

The Hamiltonian can be rewritten in the Wannier basis as

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \hat{H}=\underset{\mathbf{R}{{\mathbf{R}}^{\prime}}}{\sum}\,{{t}^{(n)}}\left(\mathbf{R}-{{\mathbf{R}}^{\prime}}\,\right)b_{n}^{\dagger}\left(\mathbf{R}\right){{b}_{n}}\left({{\mathbf{R}}^{\prime}}\,\right), \end{array}\end{eqnarray} \tag{ 13 }$$

with

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{t}^{(n)}}\left(\mathbf{R}-{{\mathbf{R}}^{\prime}}\right)=\frac{1}{{{N}_{s}}}\underset{\mathbf{k}}{\sum}\,{{E}_{n}}\left(\mathbf{k}\right){{\text{e}}^{\text{i}\mathbf{k}\centerdot \left(\mathbf{R}-{{\mathbf{R}}^{\prime}}\right)}}. \end{array}\end{eqnarray} \tag{ 14 }$$

Typical values of tunnelings (tunnelings refer to tunneling matrix elements here) for s and p bands are listed in table 1.

Table 1. Tunneling amplitudes in a two dimensional square lattice with potential $V\left(\mathbf{x}\right)=-{{V}_{0}}\left[{{\sin}^{2}}(kx)+{{\sin}^{2}}(ky)\right]$.

${{V}_{0}}/{{E}_{R}}$	$4t_{nn}^{s}/{{E}_{R}}$	$4t_{nnn}^{s}/{{E}_{R}}$	$4t_{nn}^{p}/{{E}_{R}}$	$4t_{nnn}^{p}/{{E}_{R}}$
3	−0.4441	0.0449	2.0074	0.3308
5	−0.2631	0.0136	1.6912	0.2914
10	$-0.076\,73$	9.1E-4	0.9741	0.1051
20	−9.965E  −  3	1.2E-5	0.2411	5.5E-3

Note: E_R is the one photon recoil energy $\frac{{{\hslash}^{2}}{{k}^{2}}}{2m}$. $t_{nn}^{s}$ and $t_{nnn}^{s}$ are nearest neighbor and next nearest neighbor tunnelings for the lowest s band. $t_{nn}^{p}$ and $t_{nnn}^{p}$ are nearest neighbor and next nearest neighbor tunnelings in the x direction for the p_x (first excited) band.

In the definition of Wannier functions (equation (12)), there are gauge degrees of freedom $\lambda _{\mathbf{K}}^{(n)}\left(\mathbf{k}\right)\to {{\text{e}}^{\text{i}{{\theta}_{n}}\left(\mathbf{k}\right)}}\lambda _{\mathbf{K}}^{(n)}\left(\mathbf{k}\right)$. One has to make a smooth gauge choice to get localized Wannier functions (Marzari et al 2012). Wannier functions of p-bands of a square lattice are shown in figure 1.

Zoom In Zoom Out Reset image size

Here the lifetime of p-orbital condensate in a one dimensional (1D) lattice is discussed based upon Fermi's Golden rule calculation. The resulting time scale is expected to apply to two dimensional (2D) square and three dimensional (3D) cubic lattices as well (Isacsson and Girvin 2005). The p-orbital condensate wavefunction is given as

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} | \Psi \rangle =\frac{{{\left[b_{p}^{\dagger}(Q=\pi )\right]}^{N}}}{\sqrt{N}!}|\text{vac}\rangle. \end{array}\end{eqnarray} \tag{ 15 }$$

With interactions, two particles in the p-band can collide and one particle would decay to the lowest s-band and the other goes to the second excited d-band. This process is described by the following interaction term

$$\begin{eqnarray}\begin{array}{*{35}{l}} H_{\text{int}}^{spd}= & U/{{N}_{s}}\underset{{{k}_{1}}+{{k}_{2}}+{{k}_{3}}+{{k}_{4}}=0}{\sum}\, \\ {} & \times \left\{b_{s}^{\dagger}\left({{k}_{1}}\right)b_{d}^{\dagger}\left({{k}_{2}}\right){{b}_{p}}\left({{k}_{3}}\right){{b}_{p}}\left({{k}_{4}}\right)+\text{h}.\text{c}.\right\}. \end{array}\end{eqnarray}$$

The final state after the collision is

$$\begin{eqnarray}&&|{{ \Psi }_{f}};{{k}_{1}},{{k}_{2}}\rangle =b_{s}^{\dagger}\left({{k}_{1}}\right)b_{d}^{\dagger}\left({{k}_{2}}\right)\frac{{{\left[b_{p}^{\dagger}(Q)\right]}^{N-2}}}{\sqrt{(N-2)!}}|\text{vac}\rangle.\end{eqnarray}$$

The transition probability from second order perturbation theory is

$$\begin{eqnarray}&&P\left({{k}_{1}},{{k}_{2}};t\right)\approx \frac{4{{\sin}^{2}}\left(\Delta {{\epsilon}_{{{k}_{1}}{{k}_{2}}}}t/2\right)}{ \Delta \epsilon _{{{k}_{1}}{{k}_{2}}}^{2}}|\langle {{ \Psi }_{f}};{{k}_{1}},{{k}_{2}}|H_{\text{int}}^{spd}| \Psi \rangle {{|}^{2}},\end{eqnarray}$$

where $\Delta {{\epsilon}_{{{k}_{1}}{{k}_{2}}}}$ is the difference of kinetic energy between $| \Psi \rangle$ and $|{{ \Psi }_{f}};{{k}_{1}},{{k}_{2}}\rangle$. The loss rate from the p-band is obtained as

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{w}_{t}} & = & \underset{{{k}_{1}}{{k}_{2}}}{\sum}\,\frac{1}{t}P\left({{k}_{1}},{{k}_{2}};t\right){{|}_{t\to \infty}} \\ {} & \approx & \underset{{{k}_{1}}{{k}_{2}}}{\sum}\,\frac{2\pi}{\hslash}|\langle {{ \Psi }_{f}};{{k}_{2}},{{k}_{2}}|H_{\text{int}}^{spd}| \Psi \rangle {{|}^{2}}\delta \left(\Delta {{\epsilon}_{{{k}_{1}}{{k}_{2}}}}\right) \\ {} & \approx & \frac{4\pi {{U}^{2}}}{\hslash {{N}_{s}}}N(N-1){{\left[\frac{1}{\rho \left({{\epsilon}_{s}}(K)\right)}+\frac{1}{\rho \left({{\epsilon}_{d}}(-K)\right)}\right]}^{-1}}, \end{array}\end{eqnarray} \tag{ 16 }$$

where $\rho (\epsilon )$ is the density of states and K is determined by

$$\begin{eqnarray}&&{{\epsilon}_{s}}(K)+{{\epsilon}_{d}}(-K)=2{{\epsilon}_{p}}(Q),\end{eqnarray} \tag{ 17 }$$

which in general has two solutions when the band gap between s and p matches that between p and d.

The loss rate per site is

$$\begin{eqnarray}\begin{array}{*{35}{l}} w & \approx & \frac{4\pi {{(\nu U)}^{2}}}{\hslash}{{\left[\frac{1}{\rho \left({{\epsilon}_{s}}(K)\right)}+\frac{1}{\rho \left({{\epsilon}_{d}}(-K)\right)}\right]}^{-1}}, \end{array}\end{eqnarray} \tag{ 18 }$$

with ν the filling factor. The lifetime 1/w is typically short for cubic or square lattices, where the condition of equation (17) may be satisfied. It was suggested that anharmonicity (Müller et al 2007) present in the actual optical lattice potential should help suppress the decay. Nonetheless, the lifetime can be significantly improved by using double-well lattice potentials to mismatch the band gaps as first discussed in Stojanović et al (2008) and further confirmed in the experiments (Wirth et al 2011).

3.1.1. Complex ${{p}_{x}}+\text{i}{{p}_{y}}$ Bose–Einstein condensation at finite momentum.

For bosons loaded on p-orbitals of a 2D square lattice (Wirth et al 2011) in the tight-binding regime, the tunneling Hamiltonian is

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{H}_{\text{tun}}}= & \underset{\mathbf{r}}{\sum}\,\left\{{{t}_{\parallel}}\left[b_{x}^{\dagger}\left(\mathbf{r}\right){{b}_{x}}\left(\mathbf{r}+{{{\hat{a}}}_{x}}\right)+x\leftrightarrow y\right] \right. \\ {} & \left. -\,{{t}_{\bot}}\left[b_{x}^{\dagger}\left(\mathbf{r}\right){{b}_{x}}\left(\mathbf{r}+{{{\hat{a}}}_{y}}\right)+x\leftrightarrow y\right]+\text{h}.\text{c}.\right\}, \end{array}\end{eqnarray} \tag{ 19 }$$

where b_x and b_y are bosonic annihilation operators for p_x and p_y orbitals, respectively (figure 2). After a Fourier transformation, we get the energy spectra for the p_x and p_y bands. The dispersion for the p_x band is

$$\begin{eqnarray}&&{{\epsilon}_{x}}\left(\mathbf{k}\right)=2{{t}_{\parallel}}\cos \left({{k}_{x}}\right)-2{{t}_{\bot}}\cos \left({{k}_{y}}\right).\end{eqnarray}$$

Zoom In Zoom Out Reset image size

The dispersion for the p_y band is readily obtained with a lattice rotation (C₄). There are two degenerate minima—${{\mathbf{Q}}_{x}}=(\pi,0)$ and ${{\mathbf{Q}}_{y}}=(0,\pi )$ in the p-bands with the degeneracy protected by the C₄ symmetry. The ground state manifold of non-interacting p-orbital bosons is spanned by

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} |{{N}_{x}},{{N}_{y}}\rangle =\frac{{{\left[b_{x}^{\dagger}\left({{\mathbf{Q}}_{x}}\right)\right]}^{{{N}_{x}}}}{{\left[b_{y}^{\dagger}\left({{\mathbf{Q}}_{y}}\right)\right]}^{{{N}_{y}}}}}{\sqrt{{{N}_{x}}!{{N}_{y}}!}}|\text{vac}\rangle, \end{array}\end{eqnarray} \tag{ 20 }$$

which has a large degeneracy that shall be lifted by interactions.

The interaction terms of repulsive p-orbital bosons read (Isacsson and Girvin 2005, Liu and Wu 2006, Li et al 2011b)

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{H}_{\text{int}}}= & \underset{\mathbf{r}}{\sum}\,\left\{\frac{1}{2}{{U}_{1}}\left[{{n}_{x}}\left(\mathbf{r}\right){{n}_{x}}\left(\mathbf{r}\right)+{{n}_{y}}\left(\mathbf{r}\right){{n}_{y}}\left(\mathbf{r}\right)\right]\,+2{{U}_{2}}{{n}_{x}}\left(\mathbf{r}\right){{n}_{y}}\left(\mathbf{r}\right) \right. \\ {} & \left. +\,\frac{1}{2}{{U}_{3}}\left[b_{x}^{\dagger}\left(\mathbf{r}\right)b_{x}^{\dagger}\left(\mathbf{r}\right){{b}_{y}}\left(\mathbf{r}\right){{b}_{y}}\left(\mathbf{r}\right)+\text{h}.\text{c}.\right]\right\}, \end{array}\end{eqnarray} \tag{ 21 }$$

where the density operators ${{n}_{\mu}}=b_{\mu}^{\dagger}{{b}_{\mu}}$. Approximating Wannier functions by localized harmonic wavefunctions, we have

$$\begin{eqnarray}&&{{U}_{1}}=3{{U}_{2}}=3{{U}_{3}}\equiv U>0,\end{eqnarray} \tag{ 22 }$$

from which the interaction can be rewritten as

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{H}_{\text{int}}}=\frac{U}{2}\underset{\mathbf{R}}{\sum}\,\left[{{n}^{2}}\left(\mathbf{R}\right)-\frac{1}{3}L_{z}^{2}\left(\mathbf{R}\right)\right], \end{array}\end{eqnarray} \tag{ 23 }$$

with $n={\sum}_{\nu}b_{\nu}^{\dagger}{{b}_{\nu}}$ and ${{L}_{z}}=\text{i}b_{x}^{\dagger}{{b}_{y}}+\text{h}.\text{c}.$ We thus expect that the angular momentum order is 'universally' favorable in p-orbital Bose gases.

It is however worth emphasizing here that the angular momentum ordering does not rely on the strict equality in equation (22) or the interaction form in equation (23). This becomes more clear with Ginzburg–Landau or effective field theories analysis (Li et al 2012, 2016, Liu et al 2013). Detailed studies taking into account unharmonic corrections and trapping potentials also confirm that the angular momentum order indeed exists in the regimes accessible to optical lattice experiments (Collin et al 2010, Pinheiro et al 2012, Pietraszewicz et al 2013, Sowiński et al 2013).

To capture quantum/thermal fluctuations, two slowly varying bosonic fields are introduced as

$$\begin{eqnarray}&&{{\phi}_{\mu}}\left(\mathbf{x}\right)=\underset{\mathbf{k}}{\overset{ \Lambda }{\sum}}\,{{b}_{\mu}}\left({{\mathbf{Q}}_{\mu}}+\mathbf{k}\right){{\text{e}}^{\text{i}\mathbf{k}\centerdot \mathbf{x}}},\end{eqnarray}$$

where $\Lambda$ is a momentum cut off. The effective Hamiltonian of the field theory of ${{\phi}_{\mu}}\left(\mathbf{x}\right)$ is

$$\begin{eqnarray}\begin{array}{*{35}{l}} H= & {\int}^{}{{\text{d}}^{2}}\mathbf{r}\,\left[{{K}_{1}}\left({{\partial}_{x}}\phi _{x}^{\dagger}\left(\mathbf{r}\right){{\partial}_{x}}{{\phi}_{x}}\left(\mathbf{r}\right)+x\leftrightarrow y\right) \right. \\ {} & +\,\,{{K}_{2}}\left({{\partial}_{y}}\phi _{x}^{\dagger}\left(\mathbf{r}\right){{\partial}_{y}}{{\phi}_{x}}\left(\mathbf{r}\right)+x\leftrightarrow y\right)\,-\mu \left(\phi _{x}^{\dagger}{{\phi}_{x}}+x\leftrightarrow y\right) \\ {} & +\,\frac{1}{2}{{g}_{1}}\left(\phi _{x}^{\dagger}{{\phi}_{x}}\phi _{x}^{\dagger}{{\phi}_{x}}+x\leftrightarrow y\right)\,+2{{g}_{2}}\phi _{x}^{\dagger}{{\phi}_{x}}\phi _{y}^{\dagger}{{\phi}_{y}} \\ {} & \left. +\,\frac{1}{2}{{g}_{3}}\left(\phi _{x}^{\dagger}\phi _{x}^{\dagger}{{\phi}_{y}}{{\phi}_{y}}+\text{h}.\text{c}.\right)\right]. \end{array}\end{eqnarray} \tag{ 24 }$$

In a superfluid state, we have $\langle {{\phi}_{\nu}}\rangle$  =  ${{\varphi}_{\nu}}$. At mean field level, the energy of this state is

$$\begin{eqnarray}\begin{array}{*{35}{l}} E & = & \frac{1}{2}{{g}_{1}}\left(|{{\varphi}_{x}}{{|}^{4}}+\,|{{\varphi}_{y}}{{|}^{4}}\right)+2{{g}_{2}}|{{\varphi}_{x}}{{|}^{2}}|{{\varphi}_{y}}{{|}^{2}} \\ {} & {} & +\,\frac{1}{2}{{g}_{3}}\left(\varphi _{x}^{\ast 2}\varphi _{y}^{2}+\text{c}.\text{c}.\right). \end{array}\end{eqnarray}$$

From equation (22), we have ${{g}_{1}}=3{{g}_{2}}=3{{g}_{3}}>0$, and the relative phase between p_x and p_y is locked at $\pm \pi /2$, i.e. ${{\varphi}_{x}}={{\varphi}_{y}}{{\text{e}}^{\pm \text{i}\frac{\pi}{2}}}$, where the '±' sign is spontaneously chosen. The superfluid state has a staggered angular momentum order ${{(-1)}^{{{R}_{x}}+{{R}_{y}}}}\langle {{L}_{z}}\left(\mathbf{R}\right)\rangle$, which breaks time-reversal symmetry. Such a superfluid state is named transversely staggered orbital current (TSOC) superfluid. The phase configuration of this superfluid state and its momentum distribution are shown in figure 3.

Zoom In Zoom Out Reset image size

An alternative way of looking at time-reversal symmetry breaking is to project interactions into the subspace spanned by $|{{N}_{x}}{{N}_{y}}\rangle$ in equation (20). In this subspace, the interaction reads

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {} & \langle N_{x}^{\prime}N_{y}^{\prime}|{{H}_{\text{int}}}|{{N}_{x}}{{N}_{y}}\rangle \\ {} & = & \frac{{{U}_{1}}}{2{{N}_{s}}}\left(N_{x}^{2}+N_{y}^{2}\right){{\delta}_{{{N}_{x}}N_{x}^{\prime}}}{{\delta}_{{{N}_{y}}N_{y}^{\prime}}}+2\frac{{{U}_{2}}}{{{N}_{s}}}{{N}_{x}}{{N}_{y}}{{\delta}_{{{N}_{x}}N_{x}^{\prime}}}{{\delta}_{{{N}_{y}}N_{y}^{\prime}}} \\ {} & {} & +\,\frac{{{U}_{3}}}{2{{N}_{s}}}\left[\sqrt{N_{x}^{\prime}\left(N_{x}^{\prime}-1\right){{N}_{y}}\left({{N}_{y}}-1\right)}{{\delta}_{N_{x}^{\prime}{{N}_{x}}+2}}{{\delta}_{N_{y}^{\prime}{{N}_{y}}-2}}+x\leftrightarrow y\right]. \end{array}\end{eqnarray} \tag{ 25 }$$

For the two orbital components to be miscible, we need

$$\begin{eqnarray}&&2{{U}_{2}}\,-\,|{{U}_{3}}|<{{U}_{1}},\end{eqnarray} \tag{ 26 }$$

as analogous to spin miscible condition in spinor condensates (Pethick and Smith 2008).

The angular momentum correlation is given by

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{(-1)}^{{{R}_{x}}+{{R}_{y}}}}\langle {{L}_{z}}\left(\mathbf{R}\right){{L}_{z}}(0)\rangle = & -\frac{1}{N_{s}^{2}}\left[\langle {{\left(b_{x}^{\dagger}\left({{\mathbf{Q}}_{x}}\right){{b}_{y}}\left({{\mathbf{Q}}_{y}}\right)\right)}^{2}}+\text{h}.\text{c}.\rangle \right. \\ {} & \left. -\,\langle {{n}_{x}}\rangle \left(\langle {{n}_{y}}\rangle +1\right)-\langle {{n}_{y}}\rangle \left(\langle {{n}_{x}}\rangle +1\right)\right], \end{array}\end{eqnarray} \tag{ 27 }$$

with $\langle \ldots \rangle$ the ground state expectation value. With U₃  >  0, to minimize the energy in equation (25), $\langle {{\left(b_{x}^{\dagger}\left({{\mathbf{Q}}_{x}}\right){{b}_{y}}\left({{\mathbf{Q}}_{y}}\right)\right)}^{2}}+\text{h}.\text{c}.\rangle$ gets a negative value in the ground state and in the thermodynamical limit ($\langle {{n}_{x}}\rangle \gg 1$, $\langle {{n}_{y}}\rangle \gg 1$), it approaches $(-)2\langle {{n}_{x}}\rangle \langle {{n}_{y}}\rangle$. The system thus has a long range correlation in angular momentum, i.e. ${{(-1)}^{{{R}_{x}}+{{R}_{y}}}}\langle {{L}_{z}}\left(\mathbf{R}\right){{L}_{z}}(0)\rangle \overset{|\mathbf{R}|\to \infty}{{\to}}\,\text{const}.\ne 0$. The corresponding Ising order parameter is a staggered angular momentum ${{\tilde{L}}_{z}}\left(\mathbf{R}\right)\equiv {{(-1)}^{{{R}_{x}}+{{R}_{y}}}}{{L}_{z}}\left(\mathbf{R}\right)$. When U₃ is negative, $\langle {{\left(b_{x}^{\dagger}\left({{\mathbf{Q}}_{x}}\right){{b}_{y}}\left({{\mathbf{Q}}_{y}}\right)\right)}^{2}}+\text{h}.\text{c}.\rangle$ becomes positive, and the angular momentum order $\langle {{\tilde{L}}_{z}}\rangle$ vanishes and the system develops the other Ising orbital order ${{p}_{x}}\pm {{p}_{y}}$ with an order parameter ${{(-1)}^{{{R}_{x}}+{{R}_{y}}}}\langle b_{x}^{\dagger}\left(\mathbf{R}\right){{b}_{y}}\left(\mathbf{R}\right)+\text{h}.\text{c}.\rangle$. From the above analysis, the transition at U₃  =  0 is predicted to be first order (figure 4), although fluctuations may stabilize some intermediate state and the first order transition could be replaced by a sequence of double second order transitions.

Zoom In Zoom Out Reset image size

3.1.2. Symmetry based effective field theory description.

The predicted TSOC superfluid state in the p-band tight binding model is also confirmed with effective field theory (EFT) treatment (Li et al 2016), which infers that the TSOC superfluid does not necessarily require a deep lattice. In the band structure calculation for the case of lattice rotation symmetry (Wirth et al 2011), dispersion of the relevant p-band ${{E}_{p}}\left(\mathbf{k}\right)$ has two degenerate minima at ${{\mathbf{Q}}_{x}}=(\pi,0)$ and ${{\mathbf{Q}}_{y}}=(0,\pi )$, around which low energy modes can be excited due to quantum or thermal fluctuations. This leads to a two-component EFT, where the fields are introduced as

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\phi}_{\alpha}}\left(\mathbf{x}\right)={\int}^{ \Lambda }\frac{{{\text{d}}^{2}}\mathbf{q}}{{{(2\pi )}^{2}}}b\left({{\mathbf{Q}}_{\alpha}}+q\right){{\text{e}}^{\text{i}\mathbf{q}\centerdot \mathbf{x}}}, \end{array}\end{eqnarray} \tag{ 28 }$$

with $\Lambda$ a momentum cutoff and $b\left({{Q}_{\alpha}}+q\right)$ annihilation operators for the Bloch modes near the band minima. The form of EFT is determined by considering lattice rotation and reflection symmetries, under which the fields ${{\phi}_{\alpha}}$ transform as

$$\begin{eqnarray}&&\left[\begin{array}{c} {{\phi}_{x}}(x,y) \\ {{\phi}_{y}}(x,y) \end{array}\right]\to \left[\begin{array}{c} -{{\phi}_{y}}(\,y,-x) \\ {{\phi}_{x}}(\,y,-x) \end{array}\right],\end{eqnarray} \tag{ 29 }$$

and

$$\begin{eqnarray}&&\left[\begin{array}{c} {{\phi}_{x}}(x,y) \\ {{\phi}_{y}}(x,y) \end{array}\right]\to \left[\begin{array}{c} -{{\phi}_{x}}(-x,y) \\ {{\phi}_{y}}(-x,y) \end{array}\right],\end{eqnarray} \tag{ 30 }$$

respectively. The Hamiltonian density of the EFT consistent with these symmetries is

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\mathcal{H}}_{\text{eff}}}= & \phi _{x}^{\dagger}\left(\mathbf{x}\right)\left({{K}_{\parallel}}\frac{{{\partial}^{2}}}{\partial {{x}^{2}}}+{{K}_{\bot}}\frac{{{\partial}^{2}}}{\partial {{y}^{2}}}-\mu \right){{\phi}_{x}}\left(\mathbf{x}\right)+x\to y \\ {} & +\,\underset{{{\alpha}_{1}}{{\alpha}_{2}}}{\sum}\,{{g}_{{{\alpha}_{1}}{{\alpha}_{2}}}}|{{\phi}_{{{\alpha}_{1}}}}{{|}^{2}}|{{\phi}_{{{\alpha}_{2}}}}{{|}^{2}}+{{g}_{3}}\left(\phi _{x}^{\dagger 2}\phi _{y}^{2}+\text{h}.\text{c}.\right), \end{array}\end{eqnarray} \tag{ 31 }$$

with effective couplings ${{K}_{\parallel}}$, ${{K}_{\bot}}$ and g's. This form of EFT is solely symmetry based, i.e. independent of microscopic details. For weakly interacting bosons, the coupling constants in equation (31) can be calculated from microscopic models (see appendix).

In the vicinity of thermal phase transitions of the superfluid phases, classical phase fluctuations are expected to dominate the universal physics, which allows us to ignore the subdominant density fluctuations and to replace ${{\phi}_{\alpha}}$ by $\sqrt{\rho /2}{{\text{e}}^{\text{i}{{\theta}_{\alpha}}}}$ with ρ the total density. In terms of phases ${{\theta}_{\alpha}}$, the Hamiltonian density is rewritten as

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\mathcal{H}}_{\text{eff}}}=\left(-{{K}_{\parallel}}{{\left({{\partial}_{x}}{{\theta}_{x}}\right)}^{2}}-{{K}_{\bot}}{{\left({{\partial}_{y}}{{\theta}_{x}}\right)}^{2}}+x\leftrightarrow y\right)+\frac{1}{2}{{g}_{3}}{{\rho}^{2}}\cos \left(2\left({{\theta}_{x}}-{{\theta}_{y}}\right)\right). \end{array}\end{eqnarray} \tag{ 32 }$$

Bearing in mind the periodic nature of the phases ${{\theta}_{\alpha}}$, a proper lattice regularization of this EFT leads to a coupled XY model,

$$\begin{eqnarray}\begin{array}{*{35}{l}} H_{\text{phase}}^{\text{eff}}= & \underset{\mathbf{r}}{\sum}\,\left[\left\{2{{J}_{\parallel}}\cos \left({{ \Delta }_{x}}{{\theta}_{x}}\left(\mathbf{r}\right)\right)-2{{J}_{\bot}}\cos \left({{ \Delta }_{y}}{{\theta}_{x}}\left(\mathbf{r}\right)\right)\right\} \right. \\ {} & \left. +\left\{x\leftrightarrow y\right\}\right]-U\underset{\mathbf{r}}{\sum}\,{{\sin}^{2}}\left({{\theta}_{x}}\left(\mathbf{r}\right)-{{\theta}_{y}}\left(\mathbf{r}\right)\right), \end{array}\end{eqnarray} \tag{ 33 }$$

where ${{ \Delta }_{j}}{{\theta}_{\alpha}}\left(\mathbf{r}\right)={{\theta}_{\alpha}}\left(\mathbf{r}+{{\mathbf{a}}_{j}}\right)-{{\theta}_{\alpha}}\left(\mathbf{r}\right)$ with j  =  x,y. At zero temperature we have the TSOC superfluid where the phases are locked at ${{\theta}_{x}}\left(\mathbf{r}\right)={{r}_{x}}\pi +{{\theta}_{0}}$, ${{\theta}_{y}}\left(\mathbf{r}\right)={{r}_{y}}\pi +{{\theta}_{0}}+s\pi /2$, with $s=\pm$ and ${{\theta}_{0}}\in \left[0,2\pi \right)$ spontaneously chosen. At finite temperature, the coupled XY model supports two types of topological defects. The first is a vortex in the phase ${{\theta}_{0}}$, which is a point defect with logarithmic energy cost. The second is a domain wall connecting two Ising domains with different s. Upon heating the TSOC superfluid, vortex proliferation should drive a Kosterlitz–Thouless transition and the domain wall fluctuations should drive an Ising transition. Monte Carlo study finds that the Kosterlitz–Thouless transition temperature is lower than the Ising transition (Li et al 2016).

From the effective field theory analysis, the p-orbital angular momentum order, or equivalently the $\pm \frac{\pi}{2}$ phase locking, does not rely on the precise form of the interaction (equation (23)). The requirements are g₃  >  0 and two p orbitals being miscible.

3.1.3. Population of higher bands by interaction.

Here, we will focus on condensation of weakly interacting bosons in a lattice potential. With weak interaction, the condensate is well described by Gross–Pitaevskii approach where the condensate wavefunction $\phi \left(\mathbf{x}\right)$ is obtained by minimizing an energy functional

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{E}_{\text{GP}}}={\int}^{}{{\text{d}}^{d}}\mathbf{x}{{\phi}^{\ast}}\left(\mathbf{x}\right)\left(-\frac{{{{\vec{\nabla}}}^{2}}}{2m}+V\left(\mathbf{x}\right)-\mu \right)\phi \left(\mathbf{x}\right)+g|\phi \left(\mathbf{x}\right){{|}^{4}}. \end{array}\end{eqnarray} \tag{ 34 }$$

With infinitesimal interaction, the condensate wavefunction resembles the lowest band Bloch wavefunction with lattice momentum $\mathbf{k}=0$, i.e. $\phi \left(\mathbf{x}\right)\propto {\sum}_{\mathbf{R}}{{w}_{0}}\left(\mathbf{x}-\mathbf{R}\right)$. This wavefunction preserves lattice translation and time-reversal symmetries meaning $\phi \left(\mathbf{x}\right)=\phi \left(\mathbf{x}+\mathbf{a}\right)$ and $\phi ={{\phi}^{\ast}}$. From these preserved symmetries the generic form of the condensate wavefunction with weak interaction is

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \phi \left(\mathbf{x}\right)=\underset{n}{\sum}\,{{\lambda}_{n}}\underset{\mathbf{R}}{\sum}\,{{w}_{n}}\left(\mathbf{x}-\mathbf{R}\right), \end{array}\end{eqnarray} \tag{ 35 }$$

provided that there are no first order transitions. The coefficients ${{\lambda}_{n}}$ are real and the interaction induced high band condensate is at zero lattice momentum. In terms of ${{\lambda}_{n}}$, the energy E_GP reads as

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{E}_{\text{GP}}}= & \underset{n}{\sum}\,\left({{E}_{n}}\left(\mathbf{k}=0\right)-\mu \right)\lambda _{n}^{2}+{{U}_{0000}}\lambda _{0}^{4} \\ {} & +\,4\underset{n\ne 0}{\sum}\,{{U}_{000n}}\lambda _{0}^{3}{{\lambda}_{n}}+10\underset{n\ne 0,m\ne 0}{\sum}\,{{U}_{00mn}}\lambda _{0}^{2}{{\lambda}_{m}}{{\lambda}_{n}}+\mathcal{O}\left(\lambda _{n>0}^{3}\right), \end{array}\end{eqnarray} \tag{ 36 }$$

with the interactions ${{U}_{{{\alpha}_{1}}{{\alpha}_{2}}{{\alpha}_{3}}{{\alpha}_{4}}}}$ introduced in equation (5). Minimizing this energy functional leads to

$$\begin{eqnarray}\begin{array}{*{35}{l}} \lambda _{0}^{2} & \approx & \frac{\mu -{{E}_{0}}(0)}{2{{U}_{0000}}}, \end{array}\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\lambda}_{n>0}} & \approx & -\frac{2{{U}_{000n}}\lambda _{0}^{3}}{{{E}_{n}}(0)-\mu}. \end{array}\end{eqnarray} \tag{ 38 }$$

The ratio $\frac{{{\lambda}_{n>0}}}{{{\lambda}_{0}}}$ is readily given as

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \frac{{{\lambda}_{n}}}{{{\lambda}_{0}}}\approx -\frac{{{U}_{000n}}}{{{U}_{0000}}}\left[\frac{\mu -{{E}_{0}}(0)}{{{E}_{n}}(0)-\mu}\right]. \end{array}\end{eqnarray} \tag{ 39 }$$

Physically, the high band condensate is due to competition of interaction energy and lattice potential energy—the interaction favors an extended condensate; while the potential energy favors a condensate with high density at the lattice minima. The high band population due to interaction effects is found in various settings (van Oosten et al 2003, Alon et al 2005, Guo et al 2007, Kantian et al 2007, Larson et al 2009, Dutta et al 2011, Mering and Fleischhauer 2011, Sakmann et al 2011, Zhou et al 2011a, 2011b, Hofer et al 2012, Lühmann et al 2012, łącki et al 2013).

In experiments, one can make a large fraction of high band condensate with double-well lattices, where the gap between lowest two bands is typically small and the fraction can be measured with band mapping techniques (Greiner et al 2001). Experimental evidence of this phenomenon has recently been achieved (Soltan-Panahi et al 2012). Furthermore, on general argument, the double-well lattices can have the energy gap between the ground s-band and the first excited p-bands significantly smaller than that between the first excited bands and higher bands (e.g. between p and d). This mechanism suppresses the decay by energy conservation law, making the first excited bands effectively metastable (Stojanović et al 2008). This will be further discussed in section 4.

For parity-symmetric lattices, the condensate wavefunction is parity even—$\phi \left(\mathbf{x}\right)=\phi \left(-\mathbf{x}\right)$, which implies ${{\lambda}_{{{n}_{\text{odd}}}}}=0$, n_odd referring to the parity odd bands with ${{w}_{{{n}_{\text{odd}}}}}\left(\mathbf{x}\right)=-{{w}_{{{n}_{\text{odd}}}}}\left(-\mathbf{x}\right)$. At mean field level, parity odd bands do not contribute and $\langle {{b}_{{{n}_{\text{odd}}}}}\rangle =0$. However they can form pair condensate orderings—$\langle {{b}_{n}}\left(\mathbf{r}\right){{b}_{{{n}^{\prime}}}}\left(\mathbf{r}\right)\rangle$ due to Gaussian fluctuations (Zhou et al 2011a, 2011b). The mean field state is given by $|\text{M}\rangle =\exp \left({\int}^{}{{\text{d}}^{d}}\mathbf{x}\phi (x){{\psi}^{\dagger}}\left(\mathbf{x}\right)\right)|\text{vac}\rangle.$ The effective Hamiltonian of high band modes at Gaussian level reads

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{H}_{\text{eff}}}\approx & \underset{n\ne 0,\mathbf{k}}{\sum}\,\left({{E}_{n}}\left(\mathbf{k}\right)-\mu \right)b_{n}^{\dagger}\left(\mathbf{k}\right){{b}_{n}}\left(\mathbf{k}\right) \\ {} & +\,\underset{nm\ne 0,\mathbf{k}\ne 0}{\sum}\,\left[{{U}_{00nm}}\lambda _{0}^{2}b_{n}^{\dagger}\left(\mathbf{k}\right)b_{m}^{\dagger}\left(-\mathbf{k}\right)+\text{h}.\text{c}.\right]. \end{array}\end{eqnarray} \tag{ 40 }$$

From standard perturbation theory, the correction on the mean field state from high band fluctuations is

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} -\underset{nm\ne 0,\mathbf{k}}{\sum}\,\frac{{{U}_{00nm}}\lambda _{0}^{2}b_{n}^{\dagger}\left(\mathbf{k}\right)b_{m}^{\dagger}\left(-\mathbf{k}\right)|\text{M}\rangle}{{{E}_{n}}\left(\mathbf{k}\right)+{{E}_{m}}\left(-\mathbf{k}\right)-2\mu}, \end{array}\end{eqnarray} \tag{ 41 }$$

which mediates pairings

$$\begin{eqnarray}\begin{array}{*{35}{l}} \langle {{b}_{n,\mathbf{r}}}{{b}_{m,\mathbf{r}}}\rangle = & \frac{1}{{{N}_{s}}}\underset{\mathbf{k}}{\sum}\,\langle {{b}_{n}}\left(\mathbf{k}\right){{b}_{m}}\left(-\mathbf{k}\right)\rangle \\ = & -{\int}^{}\frac{{{\text{d}}^{d}}\mathbf{k}}{{{(2\pi )}^{d}}}\frac{{{U}_{00nm}}\lambda _{0}^{2}}{{{E}_{n}}\left(\mathbf{k}\right)+{{E}_{m}}\left(-\mathbf{k}\right)-2\mu}. \end{array}\end{eqnarray} \tag{ 42 }$$

In parity odd bands, bosons form pair condensate with $\langle {{b}_{{{n}_{\text{odd}}},\mathbf{r}}}\rangle =0$ and $\langle {{b}_{{{n}_{\text{odd}}},\mathbf{r}}}{{b}_{{{n}_{\text{odd}}},\mathbf{r}}}\rangle \ne 0$.

3.1.4. Three dimensional p-orbital BEC and frustrated orbital ordering.

For bosons loaded on p-bands of a three dimensional cubic lattice, the tight binding Hamiltonian is (Liu and Wu 2006)

$$\begin{eqnarray}\begin{array}{*{35}{l}} H & = & \underset{\mathbf{r}\alpha \beta}{\sum}\,\left[{{t}_{\parallel}}{{\delta}_{\alpha \beta}}-{{t}_{\bot}}\left(1-{{\delta}_{\alpha \beta}}\right)\right]\left(b_{\alpha,\mathbf{r}+{{e}_{\nu}}}^{\dagger}{{b}_{\beta \mathbf{r}}}+\text{h}.\text{c}.\right) \\ {} & {} & +\,\frac{U}{2}\underset{\mathbf{r}}{\sum}\,\left[n_{\mathbf{r}}^{2}-\frac{1}{3}\vec{L}_{\mathbf{r}}^{2}\right], \end{array}\end{eqnarray} \tag{ 43 }$$

where n and $\vec{L}$ are boson density and angular momentum operators ${{n}_{\mathbf{r}}}={\sum}_{\alpha}b_{\alpha \mathbf{r}}^{\dagger}{{b}_{\alpha \mathbf{r}}}$ and ${{L}_{\alpha \mathbf{r}}}=-\text{i}{{\epsilon}_{\alpha \beta \gamma}}b_{\beta \mathbf{r}}^{\dagger}{{b}_{\gamma \mathbf{r}}}$. Without interaction, there are three degenerate p-bands and the energy minima are at ${{\mathbf{Q}}_{x}}=(\pi,0,0)$, ${{\mathbf{Q}}_{y}}=(0,\pi,0)$ and ${{\mathbf{Q}}_{z}}=(0,0,\pi )$. The degenerate single-particle states are $|{{\mathbf{Q}}_{\alpha}}\rangle =b_{\alpha}^{\dagger}\left({{\mathbf{Q}}_{\alpha}}\right)|\text{vac}\rangle$. Thus any condensate wavefunction of a linear superposition of $|{{\mathbf{Q}}_{\alpha}}\rangle$ (Cai et al 2012b),

$$\begin{eqnarray}&&|\vec{c}\rangle =\underset{\alpha}{\sum}\,{{c}_{\alpha}}|{{\mathbf{Q}}_{\alpha}}\rangle \end{eqnarray}$$

has the same single-particle energy. Here $\vec{c}=\left({{c}_{x}},{{c}_{y}},{{c}_{z}}\right)$ is a complex vector normalized to 1, i.e. $|\vec{c}|\,=1$. This complex vector could be parametrized as (Liu and Wu 2006)

$$\begin{eqnarray}&&\left[\begin{array}{c} {{c}_{x}} \\ {{c}_{y}} \\ {{c}_{z}} \end{array}\right]={{\text{e}}^{\text{i}\varphi -\text{i}{{\text{T}}_{\alpha}}{{\theta}_{\alpha}}}}\left[\begin{array}{c} \cos (\chi ) \\ \text{i}\sin (\chi ) \\ 0 \end{array}\right],\end{eqnarray} \tag{ 44 }$$

with ${{T}_{\alpha =x,y,z}}$ the generators of SO(3) orbital rotation in the following matrix representation: ${{\left[{{T}_{\alpha}}\right]}_{\beta \gamma}}=-\text{i}{{\epsilon}_{\alpha \beta \gamma}}$.

Although the SO(3) orbital rotation is not a symmetry of the total Hamiltonian, it keeps the interaction term invariant because ${{n}_{\mathbf{r}}}$ and $\vec{L}_{\mathbf{r}}^{2}$ are both SO(3) scalars. With a condensate at the single-particle state $|\vec{c}\rangle$, the mean field interaction energy is readily given as (Liu and Wu 2006)

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{E}_{\text{int}}}=\frac{1}{2}U{{N}_{s}}n_{0}^{2}\left[1-\frac{1}{3}{{\sin}^{2}}(2\chi )\right], \end{array}\end{eqnarray} \tag{ 45 }$$

with n₀ the boson occupation number per site. The interaction energy is minimized at $\chi =\pm \frac{\pi}{4}$. Similar to the two dimensional case, the time-reversal symmetry is spontaneously broken in the p-band condensate. The ground state manifold is $U(1)\times {{Z}_{2}}\times SO(3)$ at mean field level. The ${{Z}_{2}}\times U(1)$ degeneracy remains due to the symmetries of the Hamiltonian, whereas the SO(3) degeneracy is an artifact of the mean field theory and such a degeneracy is lifted by fluctuations through an 'order by disorder' mechanism. Cai et al (2012b) carried out a variational comparison between the two superposition states

$$\begin{eqnarray}&&|\text{planar}\rangle =\frac{1}{\sqrt{2}}\left(|{{\mathbf{Q}}_{x}}\rangle +\text{i}|{{\mathbf{Q}}_{y}}\rangle \right),\end{eqnarray}$$

and

$$\begin{eqnarray}&&|\text{diag}\rangle =\frac{1}{\sqrt{3}}\left(|{{\mathbf{Q}}_{x}}\rangle +{{\text{e}}^{\text{i}2\pi /3}}|{{\mathbf{Q}}_{y}}\rangle +{{\text{e}}^{-\text{i}2\pi /3}}|{{\mathbf{Q}}_{z}}\rangle \right).\end{eqnarray}$$

It is found that the latter has lower energy under Bogoliubov approximation.

One interesting consequence is that the angular momentum ($\vec{L}$) order in the 3D p-orbital condensate state is noncollinear, which is different from the 2D case. The polarization configuration of $\langle \vec{L}\rangle$ is shown in figure 5. Such configuration exhibits non-zero chirality defined to be

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\chi}_{ijk}}={{{\vec{L}}}_{i}}\centerdot \left({{{\vec{L}}}_{j}}\times {{{\vec{L}}}_{k}}\right), \end{array}\end{eqnarray} \tag{ 46 }$$

where ijk denote nearby three sites of the four corners of a square plaquette in a clockwise direction. With thermal/quantum fluctuations, the presence of such chiral order may lead to unconventional phase transitions (Li et al 2016).

Zoom In Zoom Out Reset image size

In a relative shallow lattice, equation (43) derived under the harmonic approximation would receive significant unharmonic corrections. There Gutzwiller calculations suggest a more exotic condensate with nematic order (Collin et al 2010), which spontaneously breaks the cubic lattice symmetry.

3.1.5. Renormalization group analysis.

Fluctuation effects on p-band condensates beyond mean field theories are studied with perturbative one-loop analysis (Liu et al 2013), where the partition function takes the form

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} Z={\int}^{}D\left[\phi _{x}^{\ast},{{\phi}_{x}},\phi _{y}^{\ast},{{\phi}_{y}}\right]{{\text{e}}^{-S\left[\phi _{x}^{\ast},{{\phi}_{x}},\phi _{y}^{\ast},{{\phi}_{y}}\right]}}, \end{array}\end{eqnarray} \tag{ 47 }$$

with

$$\begin{eqnarray}\begin{array}{*{35}{l}} S= & {\int}^{}\frac{\text{d}\omega}{2\pi}\frac{{{\text{d}}^{2}}\mathbf{k}}{{{(2\pi )}^{2}}}\underset{\alpha}{\sum}\,\phi _{\alpha}^{\ast}\left(\omega,\mathbf{k}\right)\left(-i\omega +{{\epsilon}_{\mathbf{k}}}-{{r}_{\alpha}}\right){{\phi}_{\alpha}}\left(\omega,\mathbf{k}\right) \\ {} & +\,{{{\int}^{}}_{\omega,\mathbf{k}}}\left\{\underset{\alpha}{\sum}\,g_{1}^{\alpha}\phi _{\alpha}^{\ast}(4)\phi _{\alpha}^{\ast}(3){{\phi}_{\alpha}}(2){{\phi}_{\alpha}}(1) \right. \\ {} & +\,{{g}_{2}}\phi _{y}^{\ast}(4)\phi _{x}^{\ast}(3){{\phi}_{x}}(2){{\phi}_{y}}(1) \\ {} & \left. +\,{{g}_{3}}{{\left[\phi _{x}^{\ast}(4)\phi _{x}^{\ast}(3){{\phi}_{y}}(2){{\phi}_{y}}(1)+\text{h}.\text{c}.\right]}_{{}}}\underset{{}}{\overset{_{_{_{{}}}}}{{{}}}}\,\underset{{}}{\overset{{}}{{{}}}}\,\right\}. \end{array}\end{eqnarray} \tag{ 48 }$$

Here ${{{\int}^{}}_{\omega,\mathbf{k}}}={\prod}_{j=1}^{4}{\int}^{}\frac{\text{d}\omega}{2\pi}{\int}^{}\frac{{{\text{d}}^{2}}{{\mathbf{k}}_{j}}}{{{(2\pi )}^{2}}}{{(2\pi )}^{3}}\delta \left({{\mathbf{k}}_{4}}+{{\mathbf{k}}_{3}}-{{\mathbf{k}}_{2}}-{{\mathbf{k}}_{1}}\right)\delta \left({{\omega}_{4}}+\right.$ $\left. {{\omega}_{3}}-{{\omega}_{2}}-{{\omega}_{1}}\right)$ and ${{\phi}_{\alpha}}(\,j)$ denotes ${{\phi}_{\alpha}}\left({{\omega}_{j}},{{\mathbf{k}}_{j}}\right)$. For general lattices lacking of C₄ rotational symmetry, the energy potential parameters are not equal, ${{r}_{x}}\ne {{r}_{y}}$. Performing a momentum-shell renormalization group (RG) analysis, the fields are split into fast and slow parts, $\phi _{\alpha}^{>}\left(\omega,\mathbf{k}\right){{|}_{ \Lambda /s<|\mathbf{k}|< \Lambda }}$ and $\phi _{\alpha}^{<}\left(\omega,\mathbf{k}\right){{|}_{|\mathbf{k}|< \Lambda /s}}$. Following the standard Wilsonian RG procedure (integrating out fast modes and rescaling the effective action for the slow modes), the RG flow equations (or β-functions) for the potential parameters ${{r}_{\alpha}}$ to one-loop order are obtained to be

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {} & \frac{\text{d}{{{\tilde{r}}}_{x}}}{\text{d}l}=2{{{\tilde{r}}}_{x}}+4\tilde{g}_{1}^{x} \Theta \left({{{\tilde{r}}}_{x}}-1/2\right)+{{{\tilde{g}}}_{2}} \Theta \left({{{\tilde{r}}}_{y}}-1/2\right), \\ {} & {} & \frac{\text{d}{{{\tilde{r}}}_{y}}}{\text{d}l}=2{{{\tilde{r}}}_{y}}+4\tilde{g}_{1}^{y} \Theta \left({{{\tilde{r}}}_{y}}-1/2\right)+{{{\tilde{g}}}_{2}} \Theta \left({{{\tilde{r}}}_{x}}-1/2\right). \end{array}\end{eqnarray} \tag{ 49 }$$

Here the dimensionless parameters are defined as $\tilde{r}=rm/{{ \Lambda }^{2}}$ and $\tilde{g}=gm/(2\pi )$ and $\Theta (x)$ is the Heavyside step function. The RG flow equations for the quartic couplings are

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {} & \frac{\text{d}\tilde{g}_{1}^{x}}{\text{d}l}=-2\tilde{g}_{1}^{x2}-2\tilde{g}_{3}^{2}, \\ {} & {} & \frac{\text{d}\tilde{g}_{1}^{y}}{\text{d}l}=-2\tilde{g}_{1}^{y2}-2\tilde{g}_{3}^{2}, \\ {} & {} & \frac{\text{d}{{{\tilde{g}}}_{2}}}{\text{d}l}=-\tilde{g}_{2}^{2}, \\ {} & {} & \frac{\text{d}{{{\tilde{g}}}_{3}}}{\text{d}l}=-2{{{\tilde{g}}}_{3}}\left(\tilde{g}_{1}^{x}+\tilde{g}_{2}^{y}\right). \end{array}\end{eqnarray} \tag{ 50 }$$

With bare repulsive interaction, these quartic couplings are all marginally irrelevant. However they could strongly modify the RG flow of ${{r}_{\alpha}}$ before they renormalize to zero.

In the region with ${{\tilde{r}}_{x}}(0)\geqslant \frac{1}{2}$ and ${{\tilde{r}}_{y}}(0)\geqslant \frac{1}{2}$ the solutions are

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {} & {{{\tilde{r}}}_{x}}(l)={{\text{e}}^{2l}}\left[{{{\tilde{r}}}_{x}}(0)+{\int}_{0}^{l}\text{d}{{l}^{\prime}}{{\text{e}}^{-2{{l}^{\prime}}}}\left(4\tilde{g}_{1}^{x}\left({{l}^{\prime}}\right)+{{{\tilde{g}}}_{2}}\left({{l}^{\prime}}\right)\right)\right], \\ {} & {} & {{{\tilde{r}}}_{y}}(l)={{\text{e}}^{2l}}\left[{{{\tilde{r}}}_{y}}(0)+{\int}_{0}^{l}\text{d}{{l}^{\prime}}{{\text{e}}^{-2{{l}^{\prime}}}}\left(4\tilde{g}_{1}^{y}\left({{l}^{\prime}}\right)+{{{\tilde{g}}}_{2}}\left({{l}^{\prime}}\right)\right)\right]. \end{array}\end{eqnarray} \tag{ 51 }$$

In this region, ${{\tilde{r}}_{x}}$ and ${{\tilde{r}}_{y}}$ quickly run to positive infinity. In the region with ${{\tilde{r}}_{x}}<\frac{1}{2}$ and ${{\tilde{r}}_{y}}<\frac{1}{2}$, the solutions are

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {} & {{{\tilde{r}}}_{x}}(l)={{{\tilde{r}}}_{x}}(0){{\text{e}}^{2l}}, \\ {} & {} & {{{\tilde{r}}}_{y}}(l)={{{\tilde{r}}}_{y}}(0){{\text{e}}^{2l}}, \end{array}\end{eqnarray} \tag{ 52 }$$

from which the behaviors of RG flow are also fully determined by initial values of ${{\tilde{r}}_{x,y}}$. In other regions, one-loop corrections play more important roles in making the eventual values of ${{\tilde{r}}_{x,y}}$ positive or negative. Numerical studies have found interesting regions in the phase digram where ${{\tilde{r}}_{x}}(0)<0$ and ${{\tilde{r}}_{y}}(0)>0$ (or vice versa) flow to ${{\tilde{r}}_{x}}\to +\infty$ and ${{\tilde{r}}_{y}}\to +\infty$. Depending on the flow directions of ${{\tilde{r}}_{x,y}}$, four states can be identified: (1) Complex BEC (${{\tilde{r}}_{x}}\to +\infty$, ${{\tilde{r}}_{y}}\to +\infty$); (2) p_x BEC (${{\tilde{r}}_{x}}\to +\infty$, ${{\tilde{r}}_{y}}\to -\infty$); (3) p_y BEC (${{\tilde{r}}_{x}}\to -\infty$, ${{\tilde{r}}_{y}}\to +\infty$); and (4) vacuum (${{\tilde{r}}_{x}}\to -\infty$, ${{\tilde{r}}_{y}}\to -\infty$).

The RG study sketched above does not really capture the TSOC state because g₃ flows to 0, making an illusion that quantum fluctuations wash away the phase locking between p_x and p_y components. However this is not physically correct. A more careful RG study requires introducing U(1) and Z₂ order parameters to characterize the fluctuation effects in the TSOC state.

In the strongly interacting regime, bosons localize and form Mott insulator phases. Unlike the 'featureless' s-band Mott insulators, the p-band Mott insulators have orbital degrees of freedom. Details of preparation of p-band Mott states including relaxation dynamics are studied in Challis et al (2009). The orbital ordering is governed by the orbital exchange interactions which result from virtual boson tunnelings. Here we will derive the orbital super-exchange interactions and discuss the orbital frustrations on certain lattice geometries.

3.2.1. Mott states with filling factor larger than 1.

The procedure to derive super-exchange interactions is to take the local terms as the leading part and the hopping terms as perturbation. Consider the two dimensional p-band Bose gas for example. The local interaction is given by

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{H}_{U}}=\frac{U}{2}\left({{n}^{2}}-\frac{1}{3}L_{z}^{2}\right). \end{array}\end{eqnarray} \tag{ 53 }$$

It can be verified that the angular momentum operator L_z commutes with the local interaction, i.e.

$$\begin{eqnarray}&&\left[{{L}_{z}},{{H}_{u}}\right]=0.\end{eqnarray}$$

Thus the eigenstates of the local interaction can be chosen as states with definite angular momentum. For filling factor $\nu >1$, the degenerate eigenstates with lowest energy $=\frac{{{\nu}^{2}}}{3}U$ are

$$\begin{eqnarray}\begin{array}{*{35}{l}} |+\rangle & = & \frac{{{\left(b_{\uparrow}^{\dagger}\right)}^{\nu}}}{\sqrt{\nu !}}|\text{vac}\rangle, \\ |-\rangle & = & \frac{{{\left(b_{\downarrow}^{\dagger}\right)}^{\nu}}}{\sqrt{\nu !}}|\text{vac}\rangle, \end{array}\end{eqnarray}$$

where ${{b}_{\uparrow /\downarrow}}=\frac{{{b}_{x}}\pm \text{i}{{b}_{y}}}{\sqrt{2}}$. The states $|+\rangle$ and $|-\rangle$ have angular momentum $+\nu$ and $-\nu$, respectively. On a square lattice, the tunneling Hamiltonian in the transformed basis reads

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{H}_{\text{t}}}=\underset{s,{{s}^{\prime}},\mathbf{r}}{\sum}\,\left[{{T}_{s{{s}^{\prime}}}}\left(\hat{x}\right)b_{s,\mathbf{r}+\hat{x}}^{\dagger}{{b}_{{{s}^{\prime}},\mathbf{r}}}+\text{h}.\text{c}.\right]+x\leftrightarrow y, \end{array}\end{eqnarray} \tag{ 54 }$$

with the matrices

$$\begin{eqnarray}T\left(\hat{x}\right)=\left[\begin{array}{c} \frac{{{t}_{\parallel}}-{{t}_{\bot}}}{2} & \frac{{{t}_{\parallel}}+{{t}_{\bot}}}{2} \\ \frac{{{t}_{\parallel}}+{{t}_{\bot}}}{2} & \frac{{{t}_{\parallel}}-{{t}_{\bot}}}{2} \end{array}\right],\end{eqnarray} \tag{ 55 }$$

$$\begin{eqnarray}T\left(\,\hat{y}\right)=\left[\begin{array}{c} \frac{{{t}_{\parallel}}-{{t}_{\bot}}}{2} & -\frac{{{t}_{\parallel}}+{{t}_{\bot}}}{2} \\ -\frac{{{t}_{\parallel}}+{{t}_{\bot}}}{2} & \frac{{{t}_{\parallel}}-{{t}_{\bot}}}{2} \end{array}\right].\end{eqnarray} \tag{ 56 }$$

The low energy sub-space is spanned by the product states

$$\begin{eqnarray}&&|\left\{s\left(\mathbf{r}\right)\right\}\rangle \equiv {{\otimes}_{\mathbf{r}}}|s\left(\mathbf{r}\right)\rangle,\end{eqnarray}$$

where $s\left(\mathbf{r}\right)=\pm$ and $\mathbf{r}$ runs over all lattice sites. All the states in this subspace have the same energy to leading order in U and there is thus a macroscopically huge degeneracy. The corrections due to the hopping term H_t lift the degeneracy. The first order corrections vanish because H_t does not connect any states in the low energy sub-space. The second order correction is calculated by the standard perturbation theory,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \Delta E\left(|\left\{s\left(\mathbf{r}\right)\right\}\rangle \right)=\underset{m}{\sum}\,\frac{|\langle m|{{H}_{t}}|\left\{s\left(\mathbf{r}\right)\right\}\rangle {{|}^{2}}}{{{E}^{(0)}}\left(|\left\{s\left(\mathbf{r}\right)\right\}\rangle \right)-{{E}^{(0)}}(|m\rangle )}, \end{array}\end{eqnarray} \tag{ 57 }$$

where $|m\rangle$ is a higher energy state orthogonal to the product states $|\left\{s\left(\mathbf{r}\right\}\rangle \right.$, and E⁽⁰⁾ is the leading order energy.

Keeping only tunneling between nearest neighbors as in equation (54), $\Delta E\left(|\left\{s\left(\mathbf{r}\right)\right\}\rangle \right)$ simplifies to

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \Delta E\left(|\left\{s\left(\mathbf{r}\right)\right\}\rangle \right)=\underset{\langle \mathbf{r},{{\mathbf{r}}^{\prime}}\rangle}{\sum}\, \Delta E\left(|s\left(\mathbf{r}\right)s\left({{\mathbf{r}}^{\prime}}\right)\rangle \right), \end{array}\end{eqnarray} \tag{ 58 }$$

where $\mathbf{r}$ and ${{\mathbf{r}}^{\prime}}$ are adjacent sites. Calculating the energy correction on a two-site state $\Delta E\left(|s\left(\mathbf{r}\right)s\left({{\mathbf{r}}^{\prime}}\right)\rangle \right)$ is straightforward. The energy corrections are

$$\begin{eqnarray}\begin{array}{*{35}{l}} \Delta E(|++\rangle ) & = & \Delta E(|--\rangle ) \\ {} & = & \frac{3}{4}\,\left\{\frac{\nu (\nu +1)|{{t}_{\parallel}}-{{t}_{\bot}}{{|}^{2}}}{-U}+\frac{\nu |{{t}_{\parallel}}+{{t}_{\bot}}{{|}^{2}}}{-U(\nu +1)}\right\}, \\ \Delta E(|+-\rangle ) & = & \Delta E(|-+\rangle ) \\ {} & = & \frac{3}{4}\,\left\{\frac{\nu (\nu +1)|{{t}_{\parallel}}+{{t}_{\bot}}{{|}^{2}}}{-U}+\frac{\nu |{{t}_{\parallel}}-{{t}_{\bot}}{{|}^{2}}}{-U(\nu +1)}\right\}. \end{array}\end{eqnarray} \tag{ 59 }$$

Then the correction $\Delta E\left(|\left\{s\left(\mathbf{r}\right)\right\}\rangle \right)$ is given as

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \Delta E\left(|\left\{s\left(\mathbf{r}\right)\right\}\rangle \right)=\underset{\langle \mathbf{r},{{\mathbf{r}}^{\prime}}\rangle}{\sum}\,{{J}_{\nu}}s\left(\mathbf{r}\right)s\left({{\mathbf{r}}^{\prime}}\right), \end{array}\end{eqnarray} \tag{ 60 }$$

with

$$\begin{eqnarray}&&{{J}_{\nu}}=\frac{3{{\nu}^{2}}(\nu +2)}{2(\nu +1)}\frac{{{t}_{\parallel}}{{t}_{\bot}}}{U}>0.\end{eqnarray}$$

Including this correction into the Hamiltonian, we get

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \Delta \hat{H}=\underset{\langle \mathbf{r},{{\mathbf{r}}^{\prime}}\rangle}{\sum}\,{{J}_{\nu}}{{\sigma}_{y}}\left(\mathbf{r}\right){{\sigma}_{y}}\left({{\mathbf{r}}^{\prime}}\right), \end{array}\end{eqnarray} \tag{ 61 }$$

where ${{\sigma}_{y}}$ is defined to be ${{\sigma}_{y}}={{\nu}^{-1}}P{{L}_{z}}\,P$, with P a projection operator $P=\,|+\rangle \langle +|+|-\rangle \langle -|$. The orbital super-exchange makes the staggered angular momentum ordering energetically favorable.

It should be emphasized here that the energy corrections in equation (60) actually do not depend on the orientation of the link $\mathbf{r}-{{\mathbf{r}}^{\prime}}$ and that the effective Hamiltonian in equation (61) is independent of lattice geometries. Considering p-band Mott insulators on a triangle lattice, the effective orbital model is geometrically frustrated making both of ferromagnetic and antiferromagnetic correlations suppressed.

The above analysis holds in the deep lattice regime. For a relatively shallow lattice, the degeneracy in local Hilbert space could be lifted up (Collin et al 2010). Treating such effects as perturbations, based on well established results in transverse field Ising models (Sachdev 2011) the staggered angular momentum order in the Mott state is expected to be stable when the perturbations are reasonably weak as compared to the super-exchange. But we would like to emphasize that the competition of charge (atom number for neutral atoms) and spin orders in the shallow lattice regime may alter the above speculation and lead to potentially rich physics.

3.2.2. Mott state with filling factor 1.

For Mott states with filling factor $\nu =1$, the convenient basis to calculate the super-exchange interaction is the p_x, p_y basis, rather than the ${{p}_{x}}\pm \text{i}{{p}_{y}}$ basis. The generic form of interaction in equation (21) is used here. Like deriving super-exchange for filling $\nu >1$, we need to calculate the second order corrections of nearest neighbor product states—$|1,0;1,0\rangle$, $|0,1;0,1\rangle$, $|0,1;1,0\rangle$, and $|1,0;0,1\rangle$, where a notation

$$\begin{eqnarray}\begin{array}{*{35}{l}} |{{m}_{x}},{{m}_{y}};m_{x}^{\prime},m_{y}^{\prime}\rangle = & \frac{{{\left[b_{x}^{\dagger}\left(\mathbf{r}\right)\right]}^{{{m}_{x}}}}{{\left[b_{y}^{\dagger}\left(\mathbf{r}\right)\right]}^{{{m}_{y}}}}{{\left[b_{x}^{\dagger}\left({{\mathbf{r}}^{\prime}}\right)\right]}^{m_{x}^{\prime}}}{{\left[b_{y}^{\dagger}\left({{\mathbf{r}}^{\prime}}\right)\right]}^{m_{y}^{\prime}}}}{\sqrt{{{m}_{x}}!{{m}_{y}}!m_{x}^{\prime}!m_{y}^{\prime}!}}|\text{vac}\rangle \end{array}\end{eqnarray} \tag{ 62 }$$

is adopted to save writing. The zeroth order energy of these four states is U₁. The higher energy virtual states that H_t will couple to are $\frac{1}{\sqrt{2}}(|2,0;0,0\rangle \,+\,|0,2;0,0\rangle )$, $\frac{1}{\sqrt{2}}(|2,0;0,0\rangle \,-\,|0,2;0,0\rangle )$, $|1,1;0,0\rangle$, $\frac{1}{\sqrt{2}}\left(|0,0;2,0\rangle \right.$  +   $\left. |0,0;0,2\rangle \right)$, $\frac{1}{\sqrt{2}}(|0,0;2,0\rangle \,-\,|0,0;0,2\rangle )$, and $|0,0;1,1\rangle$, with corresponding energies $2{{U}_{1}}+{{U}_{3}}$, $2{{U}_{1}}-{{U}_{3}}$, ${{U}_{1}}+2{{U}_{2}}$, $2{{U}_{1}}+{{U}_{3}}$, $2{{U}_{1}}-{{U}_{3}}$ and ${{U}_{1}}+2{{U}_{2}}$. For the link with ${{\mathbf{r}}^{\prime}}=\mathbf{r}+\hat{x}$, the second order energy corrections are given by

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {} & \Delta E(|1,0;1,0\rangle )=-2t_{\parallel}^{2}\,\left\{\frac{1}{{{U}_{1}}+{{U}_{3}}}+\frac{1}{{{U}_{1}}-{{U}_{3}}}\right\}, \\ {} & {} & \Delta E(|0,1;0,1\rangle )=0, \\ {} & {} & \Delta E(|0,1;1,0\rangle )= \Delta E(|1,0;0,1\rangle )=-\frac{t_{\parallel}^{2}}{2{{U}_{2}}}. \end{array}\end{eqnarray} \tag{ 63 }$$

(Note that the transverse tunneling is neglected here, for the reason that the longitudinal tunneling is enough to lift the degeneracy and is significantly stronger than the transverse one.) Mapping p_x and p_y orbitals to the pseudo-spin 1/2 states, $\sigma =\uparrow$ and  ↓, respectively, the effective Hamiltonian on this link reads

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{H}_{\text{x}}}={{J}_{1}}{{\sigma}_{z}}\left(\mathbf{r}\right){{\sigma}_{z}}\left(\mathbf{r}+\hat{x}\right)+{{M}_{z}}\left[{{\sigma}_{z}}\left(\mathbf{r}\right)+{{\sigma}_{z}}\left(\mathbf{r}+\hat{x}\right)\right], \end{array}\end{eqnarray} \tag{ 64 }$$

with ${{J}_{1}}=-\frac{t_{\parallel}^{2}}{2}\left[{{\left({{U}_{1}}+{{U}_{3}}\right)}^{-1}}+{{\left({{U}_{1}}-{{U}_{3}}\right)}^{-1}}-{{\left(2{{U}_{2}}\right)}^{-1}}\right]$, ${{M}_{z}}=$ $-\frac{t_{\parallel}^{2}}{2}\left[{{\left({{U}_{1}}+{{U}_{3}}\right)}^{-1}}+{{\left({{U}_{1}}-{{U}_{3}}\right)}^{-1}}\right].$ Similarly, the effective Hamiltonian for the link ${{\mathbf{r}}^{\prime}}-\mathbf{r}=\hat{y}$ is obtained as

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{H}_{\text{y}}}={{J}_{1}}{{\sigma}_{z}}\left(\mathbf{r}\right){{\sigma}_{z}}\left(\mathbf{r}+\hat{y}\right)-{{M}_{z}}\left[{{\sigma}_{z}}\left(\mathbf{r}\right)+{{\sigma}_{z}}\left(\mathbf{r}+\hat{y}\right)\right]. \end{array}\end{eqnarray} \tag{ 65 }$$

Then the total effective Hamiltonian for filling $\nu =1$ on a square lattice is

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \Delta H=\underset{\mathbf{r}}{\sum}\,{{J}_{1}}\left[{{\sigma}_{z}}\left(\mathbf{r}\right){{\sigma}_{z}}\left(\mathbf{r}+\hat{x}\right)+{{\sigma}_{z}}\left(\mathbf{r}\right){{\sigma}_{z}}\left(\mathbf{r}+\hat{y}\right)\right]. \end{array}\end{eqnarray} \tag{ 66 }$$

With ${{U}_{1}}=3{{U}_{2}}=3{{U}_{3}}$, the coupling J₁ is positive and the ground state has an antiferromagnetic ordering with alternating p-orbitals (see figure 6). For a one-dimensional lattice, H_x makes p_x orbitals favorable due to the effective Zeeman term M_z (equation (64)). We mention here that including the transverse tunneling would give rise to even richer physics, e.g. an XYZ quantum Heisenberg model can emerge (Pinheiro et al 2013).

Zoom In Zoom Out Reset image size

One key difference between filling $\nu =1$ and higher fillings is that the super-exchange interaction depends on the orientation of the link ${{\mathbf{r}}^{\prime}}-\mathbf{r}$, which makes the orbital frustration on triangle/Kagome lattices even more interesting.

3.2.3. Phase diagram of p-band Bose–Hubbard model.

The phase diagram of p-band Bose–Hubbard model in a two dimensional square lattice is studied by quantum (Hébert et al 2013) and classical (Li et al 2016) Monte-Carlo simulations. For filling of two particles per site or higher, a second order quantum phase transition from antiferromagnetic Mott state to the TSOC state is found at zero temperature in the quantum Monte-Carlo study as well as in Gutzwiller approach (Larson et al 2009, Collin et al 2010, Martikainen and Larson 2012). In the weakly interacting regime at finite temperature, the fluctuations are modeled by a phase-only model studied by the classical Monte-Carlo. It is found that the TSOC state develops a two-step phase transition to the normal state, a Kosterlitz–Thouless transition followed by a higher temperature Ising transition. Sandwiched between the two transitions is a time-reversal symmetry breaking non-superfluid intermediate state. By combining the numerical results from Monte Carlo in the weak coupling regime and the analytical exact result from mapping the Mott limit of the p-band model to the orbital equivalent of the Onsager Ising model (Onsager 1944), the phase diagram for p-band Bose–Hubbard model in two dimensions is proposed (figure 7). We would like to mention here that strong correlation effects may give rise to exotic intermediate phases between p-band Mott insulator and superfluid states (Xu and Fisher 2007).

Zoom In Zoom Out Reset image size

3.3.1. Nested Fermi surface—FFLO state.

Searches for superconducting Fulde–Ferrell–Larkin–Ovchinnikov (FFLO) phases with spatially varying order parameters have been attracting tremendous interest in both atomic gases and electronic materials. It appears that the parameter window for this novel state to occur in two or three dimensions is quite narrow for conventional settings (Radzihovsky and Sheehy 2010). (In an optical lattice, the spin imbalance window to reach the FFLO state is larger but only near half filling (Loh and Trivedi 2010).) In one dimension the parameter regime is considerably larger and much progress has been made to find FFLO phases (Yang 2001, Feiguin and Heidrich-Meisner 2007, Guan et al 2007, Hu et al 2007, Orso 2007, Parish et al 2007, Casula et al 2008, Kakashvili and Bolech 2009, Liao et al 2010). Nonetheless, long range order is prohibited due to fluctuations in one dimension, which makes it challenging to observe FFLO states even in one dimension.

Due to intrinsic anisotropy of p-orbital wavefunctions, FFLO states in p-orbital fermion systems are found to occur in a wide window in two dimensions or even in three dimensions (Cai et al 2011), not restricted to half filling. Here, we shall focus on a two dimensional square lattice, where the tunneling Hamiltonian is (Cai et al 2011)

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{H}_{0}} & = & {{t}_{\parallel}}\underset{\mathbf{r},\sigma}{\sum}\,\left[c_{x,\sigma,\mathbf{r}}^{\dagger}{{c}_{x,\sigma,\mathbf{r}+{{{\hat{a}}}_{x}}}}+c_{y,\sigma,\mathbf{r}}^{\dagger}{{c}_{y,\sigma,\mathbf{r}+{{{\hat{a}}}_{y}}}}+\text{h}.\text{c}.\right] \\ {} & {} & -\,\mu \underset{\mathbf{r},\sigma}{\sum}\,{{n}_{\sigma}}\left(\mathbf{r}\right)-\frac{h}{2}\left({{n}_{\uparrow}}\left(\mathbf{r}\right)-{{n}_{\downarrow}}\left(\mathbf{r}\right)\right), \end{array}\end{eqnarray} \tag{ 67 }$$

with ${{c}_{x,\sigma}}$ the fermionic annihilation operators for p_x(p_y) orbital with pseudo-spin $\sigma =\uparrow /\downarrow$ and ${{n}_{\sigma}}={\sum}_{\alpha}c_{\alpha,\sigma}^{\dagger}{{c}_{\alpha,\sigma}}$ the density operators for spin σ. The transverse tunneling ${{t}_{\bot}}$ ($\ll {{t}_{\parallel}}$) is neglected for simplicity and in presence of spin-imbalance, this leads to perfect nesting of p-orbital Fermi surfaces (see figure 8).

Zoom In Zoom Out Reset image size

In atomic gases the pseudo-spin components are hyperfine states. The interactions between them maybe engineered by s-wave Feshbach Resonance. With the Feshbach Resonance, a matured technique in experiments, the induced interactions are given as (Zhang et al 2010a)

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{H}_{\text{int}}} & = & U\underset{\mathbf{r}}{\sum}\,\left[{{n}_{x\uparrow,\mathbf{r}}}{{n}_{x\downarrow,\mathbf{r}}}+{{n}_{y\uparrow,\mathbf{r}}}{{n}_{y\downarrow,\mathbf{r}}}\right] \\ {} & {} & -\underset{\mathbf{r}}{\sum}\,J\left[{{{\vec{S}}}_{x,\mathbf{r}}}\centerdot {{{\vec{S}}}_{y,\mathbf{r}}}-\frac{1}{4}{{n}_{x,\mathbf{r}}}{{n}_{y,\mathbf{r}}}\right] \\ {} & {} & +\underset{\mathbf{r}}{\sum}\,{{V}_{xy}}\left[c_{x\uparrow,\mathbf{r}}^{\dagger}c_{x\downarrow,\mathbf{r}}^{\dagger}{{c}_{y\downarrow,\mathbf{r}}}{{c}_{y\uparrow,\mathbf{r}}}+\text{h}.\text{c}.\right]. \end{array}\end{eqnarray} \tag{ 68 }$$

At tree level, J, V_xy and U are related: $J=\frac{2U}{3}$, ${{V}_{xy}}=\frac{U}{3}$. With attractive interaction, U  <  0, the induced Cooper parings are ${{ \Delta }_{x,\mathbf{r}}}=\langle {{c}_{x\uparrow,\mathbf{r}}}{{c}_{x\downarrow,\mathbf{r}}}\rangle$, ${{ \Delta }_{y,\mathbf{r}}}=\langle {{c}_{y\uparrow,\mathbf{r}}}{{c}_{y\downarrow,\mathbf{r}}}\rangle$, ${{ \Delta }_{xy,\mathbf{r}}}=\langle {{c}_{x\uparrow,\mathbf{r}}}{{c}_{y\downarrow,\mathbf{r}}}\rangle$, and ${{ \Delta }_{yx,\mathbf{r}}}=\langle {{c}_{y\uparrow,\mathbf{r}}}{{c}_{x\downarrow,\mathbf{r}}}\rangle$. The BCS mean field Hamiltonian is

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{H}_{\text{M}}}= & {{H}_{0}}-U\underset{\mathbf{r}}{\sum}\,\left[{{ \Delta }_{x,\mathbf{r}}}c_{x\uparrow,\mathbf{r}}^{\dagger}c_{x\downarrow,\mathbf{r}}^{\dagger}+{{ \Delta }_{y,\mathbf{r}}}c_{y\uparrow,\mathbf{r}}^{\dagger}c_{y\downarrow,\mathbf{r}}^{\dagger}+\text{h}.\text{c}.\right] \\ {} & +\,\frac{1}{2}J\underset{\mathbf{r}}{\sum}\,\left\{{{ \Delta }_{xy,\mathbf{r}}}c_{x\downarrow,\mathbf{r}}^{\dagger}c_{y\uparrow,\mathbf{r}}^{\dagger}+{{ \Delta }_{yx,\mathbf{r}}}c_{x\uparrow,\mathbf{r}}^{\dagger}c_{y\downarrow,\mathbf{r}}^{\dagger} \right. \\ {} & \left. +\,{{ \Delta }_{yx,\mathbf{r}}}c_{y\uparrow,\mathbf{r}}^{\dagger}c_{x\downarrow,\mathbf{r}}^{\dagger}+{{ \Delta }_{xy,\mathbf{r}}}c_{x\uparrow,\mathbf{r}}^{\dagger}c_{y\downarrow,\mathbf{r}}^{\dagger}+\text{h}.\text{c}.\right\} \\ {} & -\,{{V}_{xy}}\underset{\mathbf{r}}{\sum}\,\left[{{ \Delta }_{y,\mathbf{r}}}c_{x\uparrow,\mathbf{r}}^{\dagger}c_{x\downarrow,\mathbf{r}}^{\dagger}+{{ \Delta }_{x,\mathbf{r}}}c_{y\uparrow,\mathbf{r}}^{\dagger}c_{y\downarrow,\mathbf{r}}^{\dagger}+\text{h}.\text{c}.\right]. \end{array}\end{eqnarray} \tag{ 69 }$$

The last term V_xy  <  0 locks the phase difference between ${{ \Delta }_{x}}$ and ${{ \Delta }_{y}}$ at 0, which plays a central role in making the superconducting coherence two-dimensional. Solving the gap equation numerically yields

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {} & {{ \Delta }_{xy}},{{ \Delta }_{yx}}\approx 0, \end{array}\end{eqnarray} \tag{ 70 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {} & {{ \Delta }_{x,\mathbf{r}}}={{ \Delta }_{y,\mathbf{r}}}=| \Delta |\cos \left(\mathbf{Q}\centerdot \mathbf{r}\right). \end{array}\end{eqnarray} \tag{ 71 }$$

The mean field phase diagram has been mapped out in Cai et al (2011). The parameter regime supporting FFLO states is considerably large in spin imbalanced p-band fermions.

3.3.2. Nested Fermi surface—density stripes.

Besides the superconducting stripes in FFLO states, p-orbital fermions also naturally support density stripe orders, again from Fermi surface nesting (figure 9). Stripe and checkerboard orders are found in a system of spinless fermions loaded up to p-orbital bands (Zhang et al 2012), where the Hamiltonian takes a simple form

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} H=\underset{\mathbf{r}\alpha \beta}{\sum}\,\left({{t}_{\alpha \beta}}c_{\alpha,\mathbf{r}+{{{\hat{a}}}_{\beta}}}^{\dagger}{{c}_{\alpha,\mathbf{r}}}+\text{h}.\text{c}.\right)-\,\mu \underset{\mathbf{r}\alpha}{\sum}\,{{n}_{\alpha \mathbf{r}}}+g\underset{\mathbf{r}}{\sum}\,{{n}_{x,\mathbf{r}}}{{n}_{y,\mathbf{r}}}, \end{array}\end{eqnarray} \tag{ 72 }$$

with ${{t}_{\alpha \beta}}=\left[{{t}_{\parallel}}{{\delta}_{\alpha \beta}}-{{t}_{\bot}}\left(1-{{\delta}_{\alpha \beta}}\right)\right]$. Fermi surface nesting of this system is pictorially illustrated in figure 9. Fermi surfaces of p_x and p_y bands are approximately perpendicular to each other, which greatly suppresses the Cooper instability. The reason is that for the spinless case, the onsite interaction can lead to only cross-orbital Cooper pairing that is antisymmetric in px and py orbitals. The nearly orthogonal geometry of the two Fermi surfaces now makes it impossible to condense such Cooper pairs at a single center of mass momentum. In contrast, each p_x (p_y) particle-hole pair in the density channel composed of one particle and one hole within the p_x (p_y) band benefits from the fermi surface nesting. To simultaneously condense particle-hole pairs in each orbital band, the wave-vector

$$\begin{eqnarray}&&{{\mathbf{Q}}_{1,2}}=\left(2{{k}_{\text{F}}},\pm 2{{k}_{\text{F}}}\right)\end{eqnarray}$$

is most favorable (see figure 9).

Zoom In Zoom Out Reset image size

To characterize the Fermi surface nesting effect observed in figure 9, one can look at the density–density correlations, which can be calculated by the field theory with partition function $Z=\text{Tr}{{\text{e}}^{-\beta H}}={\int}^{}D\left(\psi _{\alpha}^{\ast}\left(\mathbf{r},\tau \right){{\psi}_{\alpha}}\left(\mathbf{r},\tau \right)\right)\exp \left(-{{S}_{\text{F}}}\right)$, and the action

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{S}_{\text{F}}} & = & {\int}^{}\text{d}\tau \underset{\mathbf{r},\alpha}{\sum}\,\psi _{\alpha}^{\ast}\left(\mathbf{r},\tau \right)\left({{\partial}_{\tau}}-\mu \right){{\psi}_{\alpha}}\left(\mathbf{r},\tau \right) \\ {} & {} & +\underset{\mathbf{r}\alpha \beta}{\sum}\,{{t}_{\alpha \beta}}\left(\psi _{\alpha}^{\ast}\left(\mathbf{r}+{{{\hat{a}}}_{\beta}},\tau \right){{\psi}_{\alpha}}\left(\mathbf{r},\tau \right)+\text{h}.\text{c}.\right) \\ {} & {} & +\,g\underset{\mathbf{r}}{\sum}\,\psi _{x}^{\ast}\left(\mathbf{r},\tau \right)\psi _{y}^{\ast}\left(\mathbf{r},\tau \right){{\psi}_{y}}\left(\mathbf{r},\tau \right){{\psi}_{x}}\left(\mathbf{r},\tau \right). \end{array}\end{eqnarray} \tag{ 73 }$$

The density fields are defined as ${{\rho}_{\alpha}}\left(\mathbf{r},\tau \right)={{\psi}^{\ast}}\left(\mathbf{r},\tau \right)\psi \left(\mathbf{r},\tau \right)$, and density–density correlations are given by

$$\begin{eqnarray}&&{{ \Pi }_{\alpha \beta}}(q)=\frac{T}{{{N}_{s}}}\langle {{\rho}_{\alpha,q}}{{\rho}_{\beta,q}}\rangle,\end{eqnarray}$$

with $q=\left(\mathbf{q},\text{i}\omega \right)$,

$$\begin{eqnarray}&&{{\rho}_{\alpha q}}=\underset{\mathbf{r}}{\sum}\,{\int}^{}\text{d}\tau \rho \left(\mathbf{r},\tau \right){{\text{e}}^{\text{i}\mathbf{q}\centerdot \mathbf{r}-\text{i}\omega \tau}}=\underset{k}{\sum}\,\tilde{\psi}_{\alpha}^{\dagger}(k+q){{\tilde{\psi}}_{\alpha}}(k),\end{eqnarray}$$

and ${{\tilde{\psi}}_{\alpha}}(k)=\frac{1}{\sqrt{\beta {{N}_{s}}}}{\int}^{}\text{d}\tau {\sum}_{\mathbf{r}}{{\psi}_{\alpha}}\left(\mathbf{r},\tau \right){{\text{e}}^{-\text{i}\vec{k}\centerdot \mathbf{r}+\text{i}\omega \tau}}$. It is useful to split into two channels—number density (${{\rho}_{+}}$) and orbital density (${{\rho}_{-}}$): ${{\rho}_{\pm }}(q)={{\rho}_{x,q}}\pm {{\rho}_{y,q}}$. The correlations in these two separate channels are defined as

$$\begin{eqnarray}&&{{ \Pi }_{\pm }}=\frac{T}{{{N}_{s}}}\langle {{\rho}_{\pm }}(q){{\rho}_{\pm }}(-q)\rangle.\end{eqnarray}$$

Summing up ring diagrams under random phase approximation (RPA), the correlations are obtained as (Zhang et al 2012)

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{ \Pi }_{\pm }}\left(\mathbf{Q},0\right)=\frac{2{{\chi}^{0}}}{1\pm g{{\chi}^{0}}}, \end{array}\end{eqnarray} \tag{ 74 }$$

with ${{\chi}^{0}}$ given by

$$\begin{eqnarray}&&{{\chi}^{0}}=D\left({{E}_{\text{F}}}\right)\ln \left(\frac{{{\omega}_{D}}}{T}\right),\end{eqnarray}$$

in the limit of ${{t}_{\bot}}\to 0$. Here $D\left({{E}_{\text{F}}}\right)$ is the density of states near the Fermi surface and ${{\omega}_{D}}$ is some energy cutoff in the field theory. Due to the logarithmic divergence in ${{\chi}^{0}}$, any arbitrarily weak attractive (repulsive) interaction g  <  0 (g  >  0) induce divergence of ${{ \Pi }_{+}}$ (${{ \Pi }_{-}}$) at sufficiently low temperature. The divergence of ${{ \Pi }_{+}}$ and ${{ \Pi }_{-}}$ indicates long range ordering of charge density wave (CDW) and orbital density wave (ODW), respectively. Transitions to these density waves are studied with mean field theory, where the Hamiltonian is approximated by

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{H}_{\text{MF}}}= & \underset{\mathbf{r}\alpha \beta}{\sum}\,{{t}_{\alpha \beta}}\left(c_{\alpha,\mathbf{r}+{{{\hat{a}}}_{\beta}}}^{\dagger}+\text{h}.\text{c}.\right)-\mu \underset{\mathbf{r}\alpha}{\sum}\,{{n}_{\alpha \mathbf{r}}} \\ {} & +\,g\underset{\mathbf{r}}{\sum}\,\left({{n}_{x,\mathbf{r}}}{{M}_{y,\mathbf{r}}}+{{n}_{y,\mathbf{r}}}{{M}_{x,\mathbf{r}}}-{{M}_{x,\mathbf{r}}}{{M}_{y,\mathbf{r}}}\right), \end{array}\end{eqnarray} \tag{ 75 }$$

with ${{M}_{\alpha,\mathbf{r}}}=\langle {{n}_{\alpha,\mathbf{r}}}\rangle$. Self-consistent mean field calculations confirm that repulsive and attractive interactions favor CDW and ODW, respectively. The density patterns of these density waves are shown in figure 9.

The order parameter for the charge density wave phase is introduced by

$$\begin{eqnarray}&&\rho \left(\mathbf{r}\right)=\left[{{\phi}_{1}}{{\text{e}}^{\text{i}{{\mathbf{Q}}_{1}}\centerdot \mathbf{r}}}+{{\phi}_{2}}{{\text{e}}^{\text{i}{{\mathbf{Q}}_{2}}\centerdot \mathbf{r}}}+\text{c}.\text{c}.\right]+\text{const}.,\end{eqnarray}$$

where ${{\phi}_{1}}$ and ${{\phi}_{2}}$ are complex valued fields slowly varying in space. The phenomenological free energy reads

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} F={\int}^{}{{\text{d}}^{2}}\mathbf{r}\underset{j=1,2}{\sum}\,\left(K|\vec{\nabla}{{\phi}_{j}}{{|}^{2}}+\,r|{{\phi}_{j}}{{|}^{2}}+\,u|{{\phi}_{j}}{{|}^{4}}+\,v|{{\phi}_{1}}{{|}^{2}}|{{\phi}_{2}}{{|}^{2}}\right). \end{array}\end{eqnarray} \tag{ 76 }$$

Here, incommensurate filling is assumed and the theory has an emergent $U(1)\times U(1)$ symmetry; otherwise the terms such as $\left(\phi _{j}^{p}+\text{c}.\text{c}.\right)$ are allowed and the theory has lower symmetry.

The effective couplings, K, r, u, and v in equation (76) have been connected to microscopic parameters by field theory calculations (Zhang et al 2012). At low temperature, we have (r  <  0, u  >  2v), and the system is in a striped CDW phase with the wavevector ${{\mathbf{Q}}_{1}}$ or ${{\mathbf{Q}}_{2}}$ spontaneously chosen. Assuming ${{\mathbf{Q}}_{1}}$ is spontaneously chosen there is an algebraic long range order in ${{\phi}_{1}}$, i.e. $\langle \phi _{1}^{\ast}\left(\mathbf{r}\right){{\phi}_{1}}\left({{\mathbf{r}}^{\prime}}\,\right)\rangle \propto \frac{1}{|\mathbf{r}-{{\mathbf{r}}^{\prime}}{{|}^{\gamma}}}.$ At higher temperature it is found that the striped CDW phase first melts to a nematic phase through an Ising transition and then to normal through a Kosterlitz–Thouless transition.

3.3.3. Strongly correlated orbital models.

At half filling p-orbital fermions described by the Hamiltonian in equation (72) exhibits a Mott transition with strong repulsive interaction, which is studied in Wu (2008b) and Zhao and Liu (2008). In the fermionic Mott state, like in the bosonic case, fermions are localized on each lattice site. As a result, the low energy physics is described by an effective model of super-exchange interactions.

Considering a link $\langle \mathbf{r},{{\mathbf{r}}^{\prime}}\,\rangle$ the super-exchange interactions are determined by energy corrections on the states $|1,0;1,0\rangle$, $|0,1;0,1\rangle$, $|0,1;1,0\rangle$ and $|1,0;0,1\rangle$, where a notation is taken from equation (62) with the bosonic operators ${{b}_{\alpha}}$ replaced by fermionic ones ${{c}_{\alpha}}$. Suppose this link is in the x direction, the tunneling is then $H_{\text{t}}^{x}=-{{t}_{\parallel}}c_{x,\mathbf{r}}^{\dagger}{{c}_{x,{{\mathbf{r}}^{\prime}}}}+\text{h}.\text{c}.$, with the transverse tunneling neglected. From standard perturbation theory, the energy corrections due to virtual fermion fluctuations are

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {} & \Delta E(|0,1;0,1\rangle )= \Delta E(|1,0;1,0\rangle )=0, \end{array}\end{eqnarray} \tag{ 77 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {} & \Delta E(|1,0;0,1\rangle )= \Delta E(|0,1;1,0\rangle )=-\frac{t_{\parallel}^{2}}{U}. \end{array}\end{eqnarray} \tag{ 78 }$$

Mapping p_x (p_y) to pseudo-spin $\uparrow$ (↓) states, the super-exchange interactions are given in a compact form as

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} h_{\text{eff}}^{x}={{J}_{z}}{{\sigma}_{z}}\left(\mathbf{r}\right){{\sigma}_{z}}\left({{\mathbf{r}}^{\prime}}\,\right)+\text{const}., \end{array}\end{eqnarray} \tag{ 79 }$$

with ${{J}_{z}}=\frac{t_{\parallel}^{2}}{2U}$. Rotating p-orbitals by an angle θ, we have the following transformation

$$\begin{eqnarray}&&\left(\begin{array}{c} {{c}_{x}} \\ {{c}_{y}} \end{array}\right)\to \mathcal{U}\,(\theta )\left(\begin{array}{c} {{c}_{x}} \\ {{c}_{y}} \end{array}\right),\end{eqnarray} \tag{ 80 }$$

with

$$\begin{eqnarray}\mathcal{U}\,(\theta )=\left(\begin{array}{c} \cos \theta & \sin \theta \\ -\sin \theta & \cos \theta \end{array}\right).\end{eqnarray}$$

For a link oriented at an angle θ with respect the x axis, the super-exchange interaction reads as

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} h_{\text{eff}}^{\theta}={{J}_{z}}{{{\tilde{\sigma}}}_{z}}\left(\mathbf{r}\right){{{\tilde{\sigma}}}_{z}}\left({{\mathbf{r}}^{\prime}}\right), \end{array}\end{eqnarray} \tag{ 81 }$$

with ${{\tilde{\sigma}}_{z}}={{\mathcal{U}}^{\dagger}}(\theta ){{\sigma}_{z}}\mathcal{U}(\theta )=\sin (2\theta ){{\sigma}_{x}}+\cos (2\theta ){{\sigma}_{z}}$.

On a square lattice, the Hamiltonian describing the orbital order is

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} H_{\text{eff}}^{\text{sq}}={{J}_{z}}\underset{\mathbf{r}}{\sum}\,{{\sigma}_{z}}\left(\mathbf{r}\right){{\sigma}_{z}}\left(\mathbf{r}+{{\mathbf{a}}_{x}}\right)+{{\sigma}_{z}}\left(\mathbf{r}\right){{\sigma}_{z}}\left(\mathbf{r}+{{\mathbf{a}}_{y}}\right), \end{array}\end{eqnarray} \tag{ 82 }$$

which has the same form as the p-band Mott insulator of bosons. This Hamiltonian supports an alternating ${{p}_{x}}/{{p}_{y}}$ order as shown in figure 6. On a honeycomb lattice, the super-exchange Hamiltonian is

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} H_{\text{eff}}^{\text{hc}}={{J}_{z}}\underset{\mathbf{r},j}{\sum}\,{{T}_{j}}\left(\mathbf{r}\right){{T}_{j}}\left(\mathbf{r}+{{\mathbf{e}}_{j}}\right), \end{array}\end{eqnarray} \tag{ 83 }$$

with

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{T}_{1}} & = & -\frac{\sqrt{3}}{2}{{\sigma}_{x}}+\frac{1}{2}{{\sigma}_{z}}, \\ {{T}_{2}} & = & \frac{\sqrt{3}}{2}{{\sigma}_{x}}+\frac{1}{2}{{\sigma}_{z}}, \\ {{T}_{3}} & = & -{{\sigma}_{z}}. \end{array}\end{eqnarray} \tag{ 84 }$$

Here the summation ${\sum}_{\mathbf{r}}$ includes one set of the 'A' sublattices (figure 10). This model, dubbed quantum ${{120}^{\circ}}$ model (Wu et al 2007, Zhao and Liu 2008, Wu 2008b), is geometrically frustrated. The complication of this model originates precisely from the spatial nature of orbitals, which makes orbital degrees of freedom drastically different from real spins.

Zoom In Zoom Out Reset image size

In three dimensions on a diamond lattice, the p-orbital exchange interaction leads to an exact orbital Coulomb phase characterized by ice rules and emergent gauge structures (Chern and Wu 2011).

3.3.4. Anti-Ferromagnetic phases of spinor p-orbital fermions.

As motivated by understanding the role of magnetism in high temperature superconductors, studies of antiferromagnetic transitions in s-band fermions attracted tremendous interest, but the transition temperature is still out of reach for current cooling techniques. One way to improve the transition temperature could be provided by considering p-band fermions. The antiferromagnetic transition of spin-1/2 fermions loaded in p-bands of a 3D cubic lattice are studied in Wu and Zhai (2008), where half filling (three fermions per lattice site) is assumed. The Hamiltonian describing such a system is $H={{H}_{0}}+{{H}_{\text{int}}}$, with

$$\begin{eqnarray}&&{{H}_{0}}=\underset{\mathbf{r}\alpha \beta}{\sum}\,\left({{t}_{\alpha \beta}}c_{\alpha,\mathbf{r}+{{{\hat{a}}}_{\beta}}}^{\dagger}{{c}_{\alpha,\mathbf{r}}}+\text{h}.\text{c}.\right),\end{eqnarray}$$

and

$$\begin{eqnarray}&&{{H}_{\text{int}}}=\underset{\mathbf{r}}{\sum}\,\left\{U\underset{\alpha}{\sum}\,{{n}_{\mathbf{r},\alpha,\uparrow}}{{n}_{\mathbf{r},\alpha,\downarrow}}+W\underset{\alpha \ne \beta}{\sum}\,\left({{n}_{\mathbf{r},\alpha,\uparrow}}{{n}_{\mathbf{r},\beta,\downarrow}}\right. \right. \\ &&\left. \left. +c_{\mathbf{r},\alpha,\uparrow}^{\dagger}c_{\mathbf{r},\beta,\downarrow}^{\dagger}{{c}_{\mathbf{r},\alpha,\downarrow}}{{c}_{\mathbf{r},\beta,\uparrow}}+c_{\mathbf{r},\alpha,\uparrow}^{\dagger}c_{\mathbf{r},\alpha,\downarrow}^{\dagger}{{c}_{\mathbf{r},\beta,\downarrow}}{{c}_{\mathbf{r},\beta,\uparrow}}\right)\right\}.\end{eqnarray}$$

With strong repulsion, we can project to the low-energy subspace determined by H_int; projecting out high energy subspace will contribute to super-exchange interactions. The low-energy states are the four degenerate components of total spin-3/2: $|\uparrow \uparrow \uparrow \rangle$, $\frac{1}{\sqrt{3}}(|\downarrow \uparrow \uparrow \rangle +|\uparrow \downarrow \uparrow \rangle +|\uparrow \uparrow \downarrow \rangle ),$ $\frac{1}{\sqrt{3}}(|\downarrow \downarrow \uparrow \rangle +|\uparrow \downarrow \downarrow \rangle +|\downarrow \uparrow \downarrow \rangle ),$ and $|\downarrow \downarrow \downarrow \rangle,$ where a notation $|{{s}_{1}}{{s}_{2}}{{s}_{3}}\rangle =c_{x,{{s}_{1}}}^{\dagger}c_{y,{{s}_{2}}}^{\dagger}c_{z,{{s}_{3}}}^{\dagger}|\text{vac}\rangle$ is used. This is manifestation of the Hund's rule. It is quite involved to perform the quantum mechanical second order perturbation theory here. A more elegant way is to take the Brillouin–Wigner approximation where the super-exchange interaction is given by

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{H}_{J}}=-{{P}_{G}}{{H}_{0}}{{P}_{e}}H_{\text{int}}^{-1}{{P}_{e}}{{H}_{0}}{{P}_{G}}. \end{array}\end{eqnarray} \tag{ 85 }$$

Here P_e and P_G mean projections onto excited and low-energy subspaces, respectively.

The calculations are greatly simplified by two observations (Wu and Zhai 2008). Firstly, the hopping processes only take place within the same orbital. Secondly, all terms in H₀ acting on the low-energy subspace create eigenstates of H_int with the same excitation energy U  +  2W. Then the super-exchange Hamiltonian H_J on a link $\langle \mathbf{r},{{\mathbf{r}}^{\prime}}=\mathbf{r}+\hat{x}\rangle$ is given by

$$\begin{eqnarray}&&\underset{\alpha}{\sum}\,\frac{2|{{t}_{\alpha,x}}{{|}^{2}}}{U+2W}{{P}_{G}}\left(\underset{s{{s}^{\prime}}}{\sum}\,c_{\mathbf{r},\alpha,s}^{\dagger}{{c}_{\mathbf{r},\alpha,{{s}^{\prime}}}}c_{{{\mathbf{r}}^{\prime}},\alpha,{{s}^{\prime}}}^{\dagger}{{c}_{{{\mathbf{r}}^{\prime}},\alpha,s}}\right){{P}_{G}}.\end{eqnarray}$$

By symmetrizing $c_{\mathbf{r},\alpha,s}^{\dagger}{{c}_{\mathbf{r},\alpha,{{s}^{\prime}}}}$, one can show

$$\begin{eqnarray}&&{{P}_{G}}c_{\mathbf{r},\alpha,s}^{\dagger}{{c}_{\mathbf{r},\alpha,{{s}^{\prime}}}}{{P}_{G}}=\frac{1}{3}{{P}_{G}}\underset{\beta}{\sum}\,\left(c_{\mathbf{r},\beta,s}^{\dagger}{{c}_{\mathbf{r},\beta,{{s}^{\prime}}}}\right){{P}_{G}}.\end{eqnarray}$$

Then the super-exchange Hamiltonian is obtained to be

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{H}_{J}}=J\underset{\langle \mathbf{r},{{\mathbf{r}}^{\prime}}\rangle}{\sum}\,{{{\vec{S}}}_{\mathbf{r}}}\centerdot {{{\vec{S}}}_{{{\mathbf{r}}^{\prime}}}}, \end{array}\end{eqnarray} \tag{ 86 }$$

with ${{\vec{S}}_{\mathbf{r}}}=\frac{1}{2}{{P}_{G}}{\sum}_{\alpha s{{s}^{\prime}}}c_{\mathbf{r},\alpha,s}^{\dagger}{{\vec{\sigma}}_{s{{s}^{\prime}}}}{{c}_{\mathbf{r},\alpha,{{s}^{\prime}}}}{{P}_{G}}$, and the effective coupling

$$\begin{eqnarray}&&J=4\left(t_{\parallel}^{2}+2t_{\bot}^{2}\right)/(9U+18W)>0.\end{eqnarray} \tag{ 87 }$$

The effective description is an isotropic spin-3/2 Heisenberg model. We remind the reader that this is the model for half filling, with the full Hilbert space of each site being spanned by three p-orbitals and two spins. Hund's rule reduces the low energy subspace to the total spin-3/2 space. The ground state of the system thus has an antiferromagnetic long range order. This antiferromagnetic order is destroyed by thermal fluctuations when the temperature is above Néel temperature  ∼J.

In optical lattice experiments considerable efforts have been made to create topological bands. Neutral atoms loaded in such bands would experience effective magnetic fields due to non-trivial Berry curvatures. These experimental developments are motivated by consideration of novel quantum many-body states such as quantum Hall states and topological insulators/superconductors. While previous experiments largely focused on manipulating different hyperfine states of atoms with synthetic gauge fields, recent theoretical studies (Liu et al 2010, 2014, Liu et al 2016, Sun et al 2011, 2012b, Li et al 2013, Dutta et al 2014a, 2014b, Yin et al 2015) point to alternate ways to achieve topological bands by considering high orbital states in the optical lattices of non-standard geometry.

3.4.1. Topological sp-orbital ladder.

A one dimensional ladder composed of two chains of s and p-orbitals is shown in figure 11. The two orbitals are level in energy, and other lower orbitals are energetically separated with a large gap, and thus can be neglected when considering the sp ladder. The Hamiltonian describing this orbital ladder system is given by

$$\begin{eqnarray}{{H}_{0}}=\underset{j}{\sum}\,C_{j}^{\dagger}\left[\begin{array}{c} -{{t}_{s}} & -{{t}_{sp}} \\ {{t}_{sp}} & {{t}_{p}} \end{array}\right]{{C}_{j+1}}+\text{h}.\text{c}.-\underset{j}{\sum}\,\mu C_{j}^{\dagger}{{C}_{j}},\end{eqnarray} \tag{ 88 }$$

where $C_{j}^{\dagger}=\left[a_{s}^{\dagger}(\,j),a_{{{p}_{x}}}^{\dagger}(\,j)\right]$, with $a_{s}^{\dagger}(\,j)$ and $a_{{{p}_{x}}}^{\dagger}(\,j)$ being fermion creation operators for the s- and p_x-orbitals on the A and B chain respectively. The relative sign of the hopping amplitudes is fixed by parity symmetry of the s and p_x orbital wave functions. As depicted in figure 11, the hopping pattern plays a central role in producing a topological phase. With a proper global gauge choice, t_s, t_p and t_sp are all positive. Focusing on half filling with chemical potential $\mu =0$, the Hamiltonian is particle-hole symmetric under transformation ${{C}_{j}}\to {{(-1)}^{j}}C_{j}^{\dagger}$. Heuristically, topologically non-trivial band structure of the sp-orbital ladder may be speculated by rewriting the staggered quantum tunneling as

$$\begin{eqnarray}&&{{t}_{sp}}\underset{j}{\sum}\,\left[C_{j}^{\dagger}\left(-\text{i}{{\sigma}_{y}}\right){{C}_{j+1}}+\text{h}.\text{c}.\right]\end{eqnarray}$$

resembling the spin–orbit interactions when the s and p orbitals are mapped to pseudo-spin (1/2) states. The physics of the sp orbital ladder is also connected to the more familiar frustrated ladder with magnetic π-flux, but the sp-ladder appears much easier to realize in optical lattice experiments.

Zoom In Zoom Out Reset image size

In the momentum space, the Hamiltonian takes a suggestive form

$$\begin{eqnarray}&&\mathcal{H}(k)={{h}_{0}}(k){\mathbb 1}+\vec{h}(k)\centerdot \vec{\sigma},\end{eqnarray} \tag{ 89 }$$

where ${{h}_{0}}(k)=\left({{t}_{p}}-{{t}_{s}}\right)\cos (k)$, h_x  =  0, ${{h}_{y}}(k)=2{{t}_{sp}}\sin (k)$ and ${{h}_{z}}(k)=-\left({{t}_{p}}+{{t}_{s}}\right)\cos (k)$. Here, ${\mathbb 1}$ is the unit matrix, and ${{\sigma}_{x}}$, ${{\sigma}_{y}}$ and ${{\sigma}_{z}}$ are Pauli matrices in the two-dimensional orbital space. The energy spectrum consists of two branches,

$$\begin{eqnarray}&&{{E}_{\pm }}(k)={{h}_{0}}(k)\pm \sqrt{h_{y}^{2}(k)+h_{z}^{2}(k)}.\end{eqnarray}$$

An interesting limit is that when ${{t}_{s}}={{t}_{p}}={{t}_{sp}}$, the two bands are both completely flat. As the momentum k is varied from $-\pi$ through 0 to $+\pi$, crossing the entire Brillouin zone, the direction of the vector $\vec{h}(k)$ winds an angle of $2\pi$. In the notation of Wen (2012), the sp-orbital ladder belongs to the symmetry group $G_{++}^{-+}(U,T,C)$, as it has both particle-hole and time-reversal symmetries in addition to the usual charge U(1) symmetry. At half filling, it is characterized by an integer topological invariant, in this case the winding number 1.

A manifestation of the nontrivial band topology is existence of edge states. It is easiest to show the edge states in the flat band limit, ${{t}_{s}}={{t}_{p}}={{t}_{sp}}\equiv t$, by introducing auxiliary operators, ${{\phi}_{\pm }}(\,j)=\left[{{a}_{p}}(\,j)\pm {{a}_{s}}(\,j)\right]/\sqrt{2}$. Then the Hamiltonian only contains coupling between ${{\phi}_{+}}$ and ${{\phi}_{-}}$ of nearest neighbors,

$$\begin{eqnarray}&&{{H}_{0}}\to 2t\underset{j}{\sum}\,\phi _{-}^{\dagger}(\,j){{\phi}_{+}}(\,j+1)+\text{h}.\text{c}.\end{eqnarray}$$

One sees immediately that the edge operators ${{\phi}_{+}}(1)$ and ${{\phi}_{-}}\left({{N}_{s}}\right)$ are isolated from the bulk, i.e. decoupled from the rest of the system. These modes describe two edge states at zero energy. Away from the flat band limit, the wavefunctions of the edge states analytically constructed in Li et al (2013) are found not to confine strictly at the ends, but instead decay exponentially with a characteristic length scale

$$\begin{eqnarray}&&\xi =2/\log \left(|\left(\sqrt{{{t}_{s}}{{t}_{p}}}+{{t}_{sp}}\right)/\left(\sqrt{{{t}_{s}}{{t}_{p}}}-{{t}_{sp}}\right)|\right).\end{eqnarray}$$

Here, recall that the implicit length unit is the lattice constant along the ladder leg direction. For ${{t}_{sp}}=\sqrt{{{t}_{s}}{{t}_{p}}}$, which includes the flat band limit, the decay length ξ vanishes and we have sharply confined edge states.

A topological phase transition to a trivial insulator state can be driven by inducing a coupling between s and p orbitals, $\Delta H={{ \Delta }_{y}}{\sum}_{j}C_{j}^{\dagger}{{\sigma}_{y}}{{C}_{j}}$, which can be engineered by rotating the atoms locally on each site (Gemelke et al 2010). For the coupling strength ${{ \Delta }_{y}}$ greater than some critical value $\Delta _{y}^{c}$, Berry phase vanishes and the system becomes a trivial band insulator, and the zero energy edge states disappear. Such a phase transition can be detected by measuring the density correlation between two ends in experiments.

Regarding practical experimental realizations, careful treatments of band structures and Wannier functions are required as the details of tight binding models could receive significant corrections beyond harmonic approximations (equation (2)) (Ganczarek et al 2014). One controllable way to couple s and p orbitals is to use a one dimensional shaking lattice (Sowiński 2012, łącki and Zakrzewski 2013, Zhang and Zhou 2014, Dutta et al 2015b, Przysiężna et al 2015, Sträter and Eckardt 2015, Zhang et al 2015), which has recently been realized in experiments (Fort et al 2011, Parker et al 2013, Khamehchi et al 2016, Niu et al 2015, Weinberg et al 2015). The other way to systematically control the sp-orbital coupling is to consider a noncentrosymmetric lattice where the coupling can be turned on and off by manipulating inversion symmetries (Liu et al 2016).

3.4.2. Topological semimetal from mixing p and d orbitals.

We now turn to two dimensions and study how degeneracy of higher orbital bands may give rise to topological phases (Liu et al 2010, Sun et al 2012b). Consider a double-well optical lattice of the configuration shown in figure 12 (Sun et al 2012b). By the space group symmetry (D₄) of the lattice, the two p-orbital states (p_x and p_y) are degenerate with the lowest d-orbital (i.e. ${{d}_{{{x}^{2}}-{{y}^{2}}}}$) at high symmetry points in the momentum space. The lattice configuration is found to exhibit degenerate p and d orbitals. Considering a square lattice with three orbitals on each site (p_x, p_y and ${{d}_{{{x}^{2}}-{{y}^{2}}}}$), the Hamiltonian of the tight binding model takes the following form (Sun et al 2012b)

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{H}_{0}}= & \delta \underset{\mathbf{r}}{\sum}\,d_{\mathbf{r}}^{\dagger}{{d}_{\mathbf{r}}}-{{t}_{dd}}\underset{\mathbf{r}}{\sum}\,\left(d_{\mathbf{r}}^{\dagger}{{d}_{\mathbf{r}+{{\mathbf{a}}_{x}}}}+d_{\mathbf{r}}^{\dagger}{{d}_{\mathbf{r}+{{\mathbf{a}}_{y}}}}+\text{h}.\text{c}.\right) \\ {} & +\,{{t}_{\parallel}}\underset{\mathbf{r}}{\sum}\,\left(\,p_{x,\mathbf{r}}^{\dagger}{{p}_{x,\mathbf{r}+{{\mathbf{a}}_{x}}}}+p_{y,\mathbf{r}}^{\dagger}{{p}_{y,\mathbf{r}+{{\mathbf{a}}_{y}}}}+\text{h}.\text{c}.\right) \\ {} & -\,{{t}_{\bot}}\underset{\mathbf{r}}{\sum}\,\left(\,p_{x,\mathbf{r}}^{\dagger}{{p}_{x,\mathbf{r}+{{\mathbf{a}}_{y}}}}+p_{y,\mathbf{r}}^{\dagger}{{p}_{y,\mathbf{r}+{{\mathbf{a}}_{x}}}}+\text{h}.\text{c}.\right) \\ {} & +\,{{t}_{pd}}\underset{\mathbf{r}}{\sum}\,\left(d_{\mathbf{r}}^{\dagger}{{p}_{x,\mathbf{r}+{{\mathbf{a}}_{x}}}}-p_{x,\mathbf{r}}^{\dagger}{{d}_{\mathbf{r}+{{\mathbf{a}}_{x}}}}+d_{\mathbf{r}}^{\dagger}{{p}_{y,\mathbf{r}+{{\mathbf{a}}_{y}}}}-p_{y,\mathbf{r}}^{\dagger}{{d}_{\mathbf{r}+{{\mathbf{a}}_{y}}}}+\text{h}.\text{c}.\right), \end{array}\end{eqnarray} \tag{ 90 }$$

where ${{\mathbf{a}}_{x}}$ (${{\mathbf{a}}_{y}}$) is the lattice vector in x (y) direction, and ${{p}_{x,\mathbf{r}}}$, ${{p}_{y,\mathbf{r}}}$ and ${{d}_{\mathbf{r}}}$ are fermionic annihilation operators for p_x, p_y and ${{d}_{{{x}^{2}}-{{y}^{2}}}}$ orbitals at site $\mathbf{r}$. The amplitudes of tunneling between these orbitals at nearby sites are t_dd, ${{t}_{\parallel}}$, ${{t}_{\bot}}$, and t_pd. With a proper gauge choice, these tunneling amplitudes are all positive. Here a point group D₄ and time-reversal symmetries have been assumed. This tight binding Hamiltonian can be realized by a double-well optical lattice potential

$$\begin{eqnarray}\begin{array}{*{35}{l}} V(x,y) & = & -{{V}_{1}}\left[\cos (kx)+\cos (ky)\right] \\ {} & {} & +{{V}_{2}}\left[\cos (kx+ky)+\cos (kx-ky)\right]. \end{array}\end{eqnarray} \tag{ 91 }$$

Zoom In Zoom Out Reset image size

A typical configuration and the experimental protocol to realize it are shown in figure 12 (Sun et al 2012b). By the point group symmetry (D₄) of the lattice, the two p-orbital states (p_x and p_y) are degenerate at high symmetry points in the momentum space. By dialing the relative strength of V₁ and V₂, the two p-orbitals may be tuned to degeneracy with the lowest d-orbital (i.e. ${{d}_{{{x}^{2}}-{{y}^{2}}}}$). That corresponds to the control of the value of the band gap δ in equation (90). Band structure calculation has confirmed that the relevant physics is captured by the tight binding model (Sun et al 2012b).

In the momentum space, the tight binding Hamiltonian becomes,

$$\begin{eqnarray}&&{{H}_{0}}=\underset{\mathbf{k}}{\sum}\,\left(d_{\mathbf{k}}^{\dagger},p_{x,\mathbf{k}}^{\dagger},p_{y,\mathbf{k}}^{\dagger}\right)\mathcal{H}\left(\mathbf{k}\right)\left(\begin{array}{c} {{d}_{\mathbf{k}}} \\ {{p}_{x,\mathbf{k}}} \\ {{p}_{y,\mathbf{k}}} \end{array}\right),\end{eqnarray} \tag{ 92 }$$

with $\mathcal{H}\left(\mathbf{k}\right)$ given by

$$\begin{eqnarray}\left(\begin{array}{c} -2{{t}_{dd}}\left(\cos {{k}_{x}}+\cos {{k}_{y}}\right)+\delta & 2\text{i}{{t}_{pd}}\sin {{k}_{x}} & 2\text{i}{{t}_{pd}}\sin {{k}_{y}} \\ -2\text{i}{{t}_{pd}}\sin {{k}_{x}} & 2{{t}_{\parallel}}\cos {{k}_{x}}-2{{t}_{\bot}}\cos {{k}_{y}} & 0 \\ -2\text{i}{{t}_{pd}}\sin {{k}_{x}} & 0 & 2{{t}_{\parallel}}\cos {{k}_{y}}-2{{t}_{\bot}}\cos {{k}_{x}} \end{array}\right).\end{eqnarray}$$

Depending on the value of the energy difference δ, there are two types of band structures for this model. For $\delta >4{{t}_{dd}}+2{{t}_{\parallel}}-2{{t}_{\bot}}$, ${{d}_{{{x}^{2}}-{{y}^{2}}}}$ orbitals are weakly hybridized with p-orbitals; for $0<\delta <4{{t}_{dd}}+2{{t}_{\parallel}}-2{{t}_{\bot}}$, the orbitals are strongly hybridized. For the latter case, a band touching point between the top and middle bands shows up at $\left({{k}_{x}}=0,{{k}_{y}}=0\right)$ ($\Gamma$ point). This band touching point has non-trivial topological property, which is characterized by the Berry flux defined as the contour integral of the Berry connection in the momentum space,

$$\begin{eqnarray}&&{{\gamma}_{n}}={\oint}_{\mathcal{C}}\text{d}\mathbf{k}\centerdot {{\mathbf{A}}_{n}}\left(\mathbf{k}\right),\end{eqnarray}$$

with n the band index, $\mathcal{C}$ a close contour enclosing the band-touching point, and the Berry connection ${{\mathbf{A}}_{n}}\left(\mathbf{k}\right)=\text{i}\langle {{u}_{\mathbf{k}}}|{{\partial}_{\mathbf{k}}}|{{u}_{\mathbf{k}}}\rangle$, where $|{{u}_{\mathbf{k}}}\rangle$ is the eigenstate of the Hamiltonian $\mathcal{H}\left(\mathbf{k}\right)$. The Berry flux ${{\gamma}_{n}}$ is quantized to an integer multiplied by $2\pi$, and only two cases ${{\gamma}_{n}}=0$ or π are distinguishable without any symmetry requirement due to the gauge choice in $|{{u}_{n}}\left(\mathbf{k}\right)\rangle$. However, with space-inversion symmetry, we can restrict $I|{{u}_{n}}\left(\mathbf{k}\right)\rangle =\,|{{u}_{n}}\left(-\mathbf{k}\right)\rangle$, with I the space-inversion operator. The Berry flux then becomes well defined up to $\text{mod}~4\pi$ (Sun et al 2012b). For the band touching point considered here, ${{\gamma}_{n}}$ is $2\pi$, and this band touching is topologically protected (in presence of symmetry). Filling fermions up to such a touching point gives rise to a topological semimetal.

A more illuminating way to show the topological protection is to construct an effective two band Hamiltonian in the vicinity of $\Gamma$ point. Near this point, the ${{d}_{{{x}^{2}}-{{y}^{2}}}}$ orbital band is far below in energy and can thus be eliminated. With standard perturbation theory, the effective Hamiltonian is given to second order as (Sun et al 2012b)

$$\begin{eqnarray}{{\mathcal{H}}_{\text{eff}}}=\left(\begin{array}{c} {{\mathcal{H}}_{22}} & {{\mathcal{H}}_{23}} \\ {{\mathcal{H}}_{32}} & {{\mathcal{H}}_{33}} \end{array}\right)-\frac{1}{{{H}_{11}}-\mu}\left(\begin{array}{c} {{\mathcal{H}}_{21}}{{\mathcal{H}}_{12}} & {{\mathcal{H}}_{21}}{{\mathcal{H}}_{13}} \\ {{\mathcal{H}}_{31}}{{\mathcal{H}}_{12}} & {{\mathcal{H}}_{31}}{{\mathcal{H}}_{13}} \end{array}\right),\end{eqnarray}$$

with μ the chemical potential of the topological semimetal. Further expanding momentum around 0, the effective Hamiltonian takes the following form

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\mathcal{H}}_{\text{eff}}}=\frac{{{t}_{1}}+{{t}_{2}}}{2}\left(k_{x}^{2}+k_{y}^{2}\right){\mathbb 1}+ & 2{{t}_{3}}{{k}_{x}}{{k}_{y}}{{\sigma}_{x}}+\frac{{{t}_{1}}-{{t}_{2}}}{2}\left(k_{x}^{2}-k_{y}^{2}\right){{\sigma}_{z}}, \end{array}\end{eqnarray} \tag{ 93 }$$

where ${{t}_{1}}={{t}_{\parallel}}+\frac{4t_{pd}^{2}}{2{{t}_{\parallel}}-2{{t}_{\bot}}+4{{t}_{dd}}-\delta}$, ${{t}_{2}}=-{{t}_{\bot}}$, and ${{t}_{3}}=\frac{2t_{pd}^{2}}{2{{t}_{\parallel}}-2{{t}_{\bot}}+4{{t}_{dd}}-\delta}$. The absence of ${{\sigma}_{y}}$ component is protected by time-reversal and space-inversion symmetries. The energy gap near $\Gamma$ point is $2|\vec{h}|$, with $\vec{h}$ a planar vector $\vec{h}\left(\mathbf{k}\right)=\left(2{{t}_{3}}\,{{k}_{x}}\,{{k}_{y}},\frac{{{t}_{1}}-{{t}_{2}}}{2}\left(k_{x}^{2}-k_{y}^{2}\right)\right)$. The vector $\vec{h}$ forms a vortex configuration with winding number 2 in the momentum space. At the vortex core (the $\Gamma$ point $\mathbf{k}=(0,0)$), it is guaranteed that $\vec{h}=0$, which means the degeneracy (or band touching) point is topologically protected.

A question naturally arising is whether the required time-reversal and space inversion symmetries can spontaneously break at low temperature. Renormalization group analysis (Sun et al 2009, 2012b) points to the spontaneous symmetry breaking of time-reversal, and a state with angular momentum order $\langle \text{i}p_{x}^{\dagger}{{p}_{y}}+\text{h}.\text{c}.\rangle$ is stabilized at low temperature, if the interaction is repulsive. Taking this into the effective Hamiltonian, a gap opens at $\Gamma$ point. As a result, the topological semimetal gives way to an insulator state at low temperature. This insulator is topologically non-trivial with finite Chern number. If the bare interaction is attractive, the renormalization equation shows that it flows to the fixed point of zero (usually called a marginally irrelevant term). In other words, the topological semimetal phase is stable against any attractive interaction in the perturbative renormalization group sense.

3.4.3. Nearly flatbands with nontrivial topology.

In the system described by equation (92) at low temperature, the developed angular momentum order generates an additional coupling between two p orbitals,

$$\begin{eqnarray} &&\Delta H=\underset{\mathbf{k}}{\sum}\,\text{i} \Delta p_{x,\mathbf{k}}^{\dagger}{{p}_{y,\mathbf{k}}}+\text{h}.\text{c}.,\end{eqnarray}$$

which breaks time-reversal symmetry and thus allows the Chern number to be non-trivial. With the parameter choice $\delta =-4{{t}_{dd}}+2{{t}_{\parallel}}+ \Delta -2{{t}_{\parallel}} \Delta /\left(4{{t}_{\parallel}}+ \Delta \right)$ and ${{t}_{\bot}}={{t}_{\parallel}} \Delta /\left(4{{t}_{\parallel}}+ \Delta \right)$, the energies of the top band at $\Gamma$ and M points are equal (Sun et al 2011). Varying $\Delta$ with ${{t}_{dd}}={{t}_{pd}}={{t}_{\parallel}}=t$ fixed, they found that the ratio of the bandwidth/band gap is minimized ($\approx 1/20$) at $\Delta /t=2.8$ for the top band. The top and bottom bands carry opposite Chern numbers  ±1, and are thus topologically non-trivial while the middle band is topologically trivial. Such nearly flatbands with nontrivial topology mimic the Landau levels of 2D electron gas in strong magnetic fields. The flatness is crucial to reach fractional topological states in lattice models. Further numerical investigations have shown that fractional quantum Hall states are supported when the flatbands are filled at certain fractional filings (Sheng et al 2011, Wang et al 2011).

Earlier discussion of the lifetime of p-orbital BEC based on Fermi's golden rule calculation largely relies on single-particle picture, and may underestimate many-body effects. Numerical studies based upon Gross–Pitaevskii approach (equation (34)) (Martikainen 2011, Xu et al 2013) indeed find interesting phenomena beyond the scope of Fermi's golden rule treatment. To study the TSOC superfluid in continuous space, a variational condensate wavefunction is taken,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\psi}_{0}}\left(\mathbf{x}\right)=\underset{\mathbf{K}}{\sum}\,\left[{{X}_{\mathbf{K}}}{{\text{e}}^{\text{i}\left({{\mathbf{Q}}_{x}}+\mathbf{K}\right)\centerdot \mathbf{x}}}+\text{i}{{Y}_{\mathbf{K}}}{{\text{e}}^{\text{i}\left({{\mathbf{Q}}_{y}}+\mathbf{K}\right)\centerdot \mathbf{x}}}\right], \end{array}\end{eqnarray} \tag{ 94 }$$

with $\mathbf{K}$ the reciprocal lattice vectors, ${{X}_{\mathbf{K}}}$ and ${{Y}_{\mathbf{K}}}$ the variational parameters, and ${{\mathbf{Q}}_{x}}$ (${{\mathbf{Q}}_{y}}$) the minima for the p_x (p_y) band. This wavefunction is superposed of two Bloch functions, and it breaks lattice translation symmetry, as required to describe the TSOC superfluid. The key features of TSOC superfluid which are time-reversal symmetry breaking and staggered orbital current, as predicted based on tight binding models, are confirmed in the numerical calculations for continuous space (Xu et al 2013).

The stabilities of the TSOC superfluid state are investigated within the time-dependent GP equation,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \text{i}{{\partial}_{\tau}}\psi \left(\mathbf{x},\tau \right)=\left[-\frac{{{{\vec{\nabla}}}^{2}}}{2M}+V\left(\mathbf{x}\right)+g|\psi {{|}^{2}}\right]\psi \left(\mathbf{x},\tau \right). \end{array}\end{eqnarray} \tag{ 95 }$$

Rewriting $\psi \left(\mathbf{r},\tau \right)$ into condensate and fluctuation parts,

$$\begin{eqnarray}&&\psi \left(\mathbf{x},\tau \right)={{\psi}_{0}}\left(\mathbf{x},\tau \right)+{{u}_{\mathbf{q}}}\left(\mathbf{x},\tau \right){{\text{e}}^{\text{i}\mathbf{q}\centerdot \mathbf{x}}}+v_{\mathbf{q}}^{\ast}\left(\mathbf{x},\tau \right){{\text{e}}^{-\text{i}\mathbf{q}\centerdot \mathbf{x}}}\end{eqnarray}$$

the time-dependent GP equation determines the dynamics of fluctuations (Xu et al 2013)

$$\begin{eqnarray}&&\text{i}{{\partial}_{\tau}}\left(\begin{array}{c} {{u}_{\mathbf{q}}}\left(\mathbf{x},\tau \right) \\ {{v}_{\mathbf{q}}}\left(\mathbf{x},\tau \right) \end{array}\right)={{\sigma}_{z}}{{\mathcal{K}}_{\mathbf{q}}}\left(\begin{array}{c} {{u}_{\mathbf{q}}}\left(\mathbf{x},\tau \right) \\ {{v}_{\mathbf{q}}}\left(\mathbf{x},\tau \right) \end{array}\right),\end{eqnarray} \tag{ 96 }$$

with

$$\begin{eqnarray}{{\mathcal{K}}_{\mathbf{q}}}=\left(\begin{array}{c} \mathcal{L}\left(\mathbf{q}\right) & g\psi _{0}^{2} \\ g\psi _{0}^{\ast 2} & \mathcal{L}\left(\mathbf{q}\right) \end{array}\right),\\ &&\mathcal{L}\left(\mathbf{q}\right)=-\frac{{{\left(\vec{\nabla}+\text{i}\mathbf{q}\right)}^{2}}}{2M}+V\left(\mathbf{x}\right)+2g|{{\psi}_{0}}{{|}^{2}}.\end{eqnarray} \tag{ 97 }$$

Note that the vector $\mathbf{q}$ is the lattice momentum after doubling periods to make $\mathbf{q}$ a good quantum number, and that ${{u}_{\mathbf{q}}}$ and ${{v}_{\mathbf{q}}}$ are periodic—${{u}_{\mathbf{q}}}\left(\mathbf{x}+2{{\mathbf{a}}_{x}}\right)={{u}_{\mathbf{q}}}\left(\mathbf{x}+2{{\mathbf{a}}_{y}}\right)={{u}_{\mathbf{q}}}\left(\mathbf{x}\right)$, ${{v}_{\mathbf{q}}}\left(\mathbf{x}+2{{\mathbf{a}}_{x}}\right)={{v}_{\mathbf{q}}}\left(\mathbf{x}+2{{\mathbf{a}}_{y}}\right)={{v}_{\mathbf{q}}}\left(\mathbf{x}\right)$. The eigenvalues of ${{\sigma}_{z}}{{\mathcal{K}}_{\mathbf{q}}}$ determine the Bogoliubov spectra, which are studied for square and checkerboard lattices. The fluctuations would grow in time if the eigenvalues are imaginary, leading to dynamical instability. This instability is cross checked by simulating real time dynamics in the continuous space where the optical lattice is treated exactly by a periodic potential (Xu et al 2013), beyond the standard tight-binding model approximation.

For a square lattice, the TSOC superfluid state is found to be dynamically unstable unless the interaction strength is extremely weak. In presence of dynamical instability, the lifetime of the TSOC superfluid state in a simple square lattice could be tens of milliseconds, rendering that such a state is experimentally unreachable for the simple square lattice. This conclusion is fully consistent with the early experimental finding of a relatively fast decay of the p-orbital atoms in a quasi-1D lattice system (Müller et al 2007). In contrast, for the checkerboard lattice as used in experiments (Ölschläger et al 2011, Wirth et al 2011), when the lattice is not too shallow and the interaction is not too strong, the TSOC superfluid state is shown numerically to be dynamically stable. This is consistent with the long lifetime as observed in experiments. Similar improvement with superlattices is also found in one dimension (Martikainen 2011). When the interaction is stronger than some critical value, the TSOC superfluid is no longer dynamically stable even for the checkerboard lattice. Based on the dynamical stability, a phase diagram is predicted in Xu et al (2013), which is consistent with experimental observations.

Another way to understand the dynamical instability is to look at the energy cost for fluctuations ${{u}_{\mathbf{q}}}$, ${{v}_{\mathbf{q}}}$, which takes the following form Wu and Niu (2001),

$$\begin{eqnarray}&&\delta {{E}_{\mathbf{q}}}={\int}^{}{{\text{d}}^{2}}\mathbf{x}\left(u_{\mathbf{q}}^{\ast}\left(\mathbf{x}\right),v_{\mathbf{q}}^{\ast}\left(\mathbf{x}\right)\right){{\mathcal{K}}_{\mathbf{q}}}\left(\begin{array}{c} {{u}_{\mathbf{q}}}\left(\mathbf{x}\right) \\ {{v}_{\mathbf{q}}}\left(\mathbf{x}\right) \end{array}\right).\end{eqnarray} \tag{ 98 }$$

The fact that the eigenvalues of ${{\sigma}_{z}}{{\mathcal{K}}_{\mathbf{q}}}$ are imaginary implies the matrix ${{\mathcal{K}}_{\mathbf{q}}}$ is not positive definite (although the reverse may not be true), which means that the variational ansatz in equation (94) is not a stable saddle point of the GP energy functional. This in principle indicates tendency of forming some crystalline ordering (Li et al 2011a).

The other type of instability is Landau instability for the reason that there are always Bogoliubov modes causing the free energy to be negative for p-orbital BEC, which means the state is a local saddle point that can decay into the lowest s-band. However this instability is less important than the dynamical instability within the lifetime of experiments. The time scale for Landau instability to destroy the p-orbital BEC is estimated to be 500 ms while it is found to be around 10 ms in numerical simulations for dynamical instability. Although the p-orbital BEC is not strictly a metastable state due to Landau instability, it is fairly stable within the experimentally relevant time-scale. In the checkerboard lattice experiment (Wirth et al 2011) where each lattice site actually represents an elongated tube in the third direction, the dynamical phenomena are even richer. For example, a collision process with two atoms decaying into the lowest band is allowed as the energy could be released to the kinetic motion in the third direction (Paul and Tiesinga 2013).

The dynamical instability of excited band condensate in a double-well lattice has also been studied in detail, and the loop structure in Bogoliubov spectra is found to be correlated with the dynamical instability (Hui et al 2012).

Coherent bosonic cold atoms were observed in the higher bands of an optical lattice in the pioneering experiments of accelerating lattices (Browaeys et al 2005) and of cross-band Raman transitions (Müller et al 2007).

In the experiment of Müller et al (2007), the sample is prepared by first loading a Bose–Einstein condensate of ⁸⁷Rb atoms into a deep symmetrically simple cubic 3D optical lattice formed by three far detuned laser standing waves. For this deep lattice, it can be treated as an array of 3D harmonic oscillators with discrete vibrational levels, which can be labeled as $|{{m}_{x}}{{m}_{y}}{{m}_{z}}\rangle$ with m_j the vibrational quantum number along the j axis. Population transfer in these orbital levels can be controlled using a stimulated two-photon Raman process with propagating laser beams along the x axis (see figure 13), which provides an inter-orbital coupling

$$\begin{eqnarray}&&{{ \Omega }_{\text{eff}}}|m_{x}^{\prime}{{m}_{y}}{{m}_{z}}\rangle \langle {{m}_{x}},{{m}_{y}}{{m}_{z}}|,\end{eqnarray}$$

with ${{ \Omega }_{\text{eff}}}$ the effective Rabi frequency. The experiment restricts the Raman coupling to the lowest Bloch bands and demonstrates orbital transition from the $|0\,0\,0\rangle$ state (s-orbital) to $|100\rangle$ (p_x-orbital). Rabi oscillations between the two orbitals have been observed. A maximal transfer efficiency of nearly $80 \%$ is achieved.

Zoom In Zoom Out Reset image size

The decay of atoms into the lowest orbital due to collisional events has also been measured. The lifetime was found to be 10–100 times longer than the tunneling scale. Emergence of coherence compatible with a Bose–Einstein condensation to a nonzero momentum state has been seen; yet the experimental system was anisotropic and the predicted ${{p}_{x}}+\text{i}{{p}_{y}}$-wave condensate was not studied for the absence of p_x and p_y orbital symmetry.

After the early observation of higher band population (Müller et al 2007, Johnson et al 2009), long-lived Bose–Einstein condensate in the high-bands was not achieved until the groundbreaking experiment (Wirth et al 2011). In this experiment, a square optical lattice, composed of two classes (A and B) of (tube-shaped) lattice sites is used (see figure 14). Formed by two standing waves oriented along the x and y axes with polarization along the z axis, the lattice potential is

$$\begin{eqnarray}\begin{array}{*{35}{l}} V(x,y) & = & -\frac{{{V}_{0}}}{4}|\eta \left[\left(\hat{z}\cos (\alpha )+\hat{y}\sin (\alpha )\right){{\text{e}}^{\text{i}kx}}+\epsilon \hat{z}{{\text{e}}^{-\text{i}kx}}\right] \\ {} & {} & +\,{{\text{e}}^{\text{i}\theta}}\hat{z}\left({{\text{e}}^{\text{i}ky}}+\epsilon {{\text{e}}^{-\text{i}ky}}\right){{|}^{2}}, \end{array}\end{eqnarray} \tag{ 99 }$$

where $\eta \approx 0.95$ accounts for a small difference in the powers directed to interferometer branches, $\epsilon \approx 0.81$ accounts for the imperfect retro-reflections, and the angle α permits tunability of anisotropy in the x-y plane. An isotropic p-band with degenerate band minima arises when $\cos (\alpha )\approx \epsilon$ (or $\alpha ={{\alpha}_{\text{iso}}}\approx \pi /5$). The controllability of the phase difference θ allows to adjust the relative depth of potentials at A and B sites, which is crucial in this experiment to populate higher bands. For $\theta <\pi /2$ the A sites are shallower than the B sites and vice versa.

Zoom In Zoom Out Reset image size

Initially a Bose–Einstein condensate of rubidium (⁸⁷Rb) atoms is prepared and the lattice potential is adiabatically turned on with $\theta =0.38\pi$ such that B-sites are much lower than A. A lowest band lattice Bose–Einstein condensate is thus created with most of atoms confined in B sites. Then θ is rapidly increased to a final value ${{\theta}_{f}}>\pi /2$ such that the s-orbitals in the B sites are level with the p-band of the lattice in energy. In doing so, atoms are efficiently transfered to the p-band. Since this preparation procedure is abrupt, the prepared state is not immediately a condensate state but rather an incoherent state, in which the atomic distribution in the Brillouin zone is fairly uniform. Surprisingly, after some holding time around 10 ms, sharp peaks arise at p-band minima and the p-band Bose–Einstein condensate spontaneously emerges. In theory the emergence of phase coherence is beyond the scope of Gross–Pitaevskii approach, and can be studied by constructing a quantum rotor model (Sau et al 2012), where the dynamics is well captured by the truncated Wigner approximation (Polkovnikov et al 2002).

The p-band condensate is not a ground state of the system but a metastable state; decaying into the lowest band is unavoidable. In this checkerboard lattice, the band gap between p-band and the lowest band is largely mismatched with the gap between p-band and the higher band, and Fermi's golden rule calculation (see section 2.3) predicts a significant improvement of stability. In the experiment, the lifetime of the p-band condensate could reach 100 ms or longer.

From the measurements of momentum distribution, the experimental evidence of p-band condensate is conclusive. However there is no direct evidence for the orbital ordering in the TSOC state as predicted in theory. As a step further, a phase diagram is mapped out with varying α (controlling the anisotropy) and the phase diagram is quantitatively consistent with theoretical predictions (Ölschläger et al 2013). The remarkable consistency of experimental observations with theories strongly suggests the p-band condensate be a TSOC state. Yet, direct evidence of the orbital order requires further experimental investigation.

Population of even higher bands, say f-bands, is also achieved in this checkerboard lattice (Ölschläger et al 2011) thanks to the tunability of relative depth between two sublattices. Similar procedure was implemented as in preparing the p-band condensate. The resulting f-band condensate also has a complex nature. The condensate wavefunction locally resembles the superposition ${{\psi}_{[3,0]}}\pm \text{i}{{\psi}_{[0,3]}}$ of eigenfunctions ${{\psi}_{[n,m]}}$ of a 2D harmonic oscillator with n and m oscillator quanta in x and y directions, which has a spatial $\left(2{{x}^{3}}-3x\right)\pm \text{i}\left(2{{y}^{3}}-3y\right)$ dependence locally. The complex f-band condensate emerges from the same mechanism as the TSOC state of the p-band, namely maximizing the local angular momentum.

Besides the way of loading atoms into the excited bands demonstrated in the checkerboard lattice, there are other possibilities, for example by Bloch oscillation techniques (Larson and Martikainen 2011, Tarruell et al 2012) or by vibrating lattices (Sowiński 2012, łącki and Zakrzewski 2013).

Observations of higher bands in optical lattices are achieved in the early experiments manipulating double-well lattices, which were largely motivated by implementing coherent control of quantum degrees of freedom (Sebby-Strabley et al 2006, Anderlini et al 2007, Cheinet et al 2008, Lundblad et al 2008, Trotzky et al 2008).

Here we use the experiment (Sebby-Strabley et al 2006) to demonstrate how the higher bands are populated in double-well lattices and what consequent observables are achieved. This double-well lattice is a two dimensional lattice formed by superimposing two lattices with orthogonal polarizations. Having a laser setup as shown in figure 15(a), the electric field generated by the four laser beams is $\text{Re}\left[\vec{E}(x,y)\right]{{\text{e}}^{\text{i}\omega t}}$, with

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \vec{E}(x,y)=E\left({{\text{e}}^{\text{i}kx}}+{{\text{e}}^{\text{i}(2\theta +2\phi -kx)}}\right){{{\hat{e}}}_{1}}+E\left({{\text{e}}^{\text{i}(\theta -ky)}}+{{\text{e}}^{\text{i}(\theta +2\phi +ky)}}\right){{{\hat{e}}}_{2}}, \end{array}\end{eqnarray} \tag{ 100 }$$

where $k=2\pi /\lambda$ (λ is the wavelength of the laser light), $\theta =k{{d}_{1}}+\delta \theta$, and $\phi =k{{d}_{2}}+\delta \phi$ (the extra phase shifts $\delta \theta$ and $\delta \phi$ are polarization dependent and can be controlled in experiments). We have neglected several imperfections such as imperfect alignment and reflections for simplicity here. In experiments these imperfections could cause technical challenges. For light polarizations being all in plane such that ${{\hat{e}}_{1}}=\hat{y}$, ${{\hat{e}}_{2}}=\hat{x}$, we have a laser intensity field

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{I}_{xy}}(x,y)/{{I}_{xy,0}}=2\cos \left(2kx-2{{\theta}_{xy}}-2{{\phi}_{xy}}\right)+2\cos \left(2ky+2{{\phi}_{xy}}\right)+4, \end{array}\end{eqnarray} \tag{ 101 }$$

with subscripts in $\theta$ and ϕ specifying the polarization dependence. For the out-of-plane case, ${{\hat{e}}_{1}}={{\hat{e}}_{2}}=\hat{z}$, the laser intensity field is

$$\begin{eqnarray}&&\begin{array}{*{20}{l}}{{I_z}(x,y)/{I_{z,0}} = 16{{\left[ {\cos (\frac{k}{2}(x + y) - \frac{{{\theta _z}}}{2}} \right]}^2}\,{{\left[ {\cos (\frac{k}{2}(x - y) - \frac{{{\theta _z}}}{2} - {\phi _z}} \right]}^2}.}\end{array} \end{eqnarray} \tag{ 102 }$$

Zoom In Zoom Out Reset image size

The laser field creates an optical potential $V(x,y)\propto$ $\left({{I}_{xy}}(x,y)+{{I}_{z}}(x,y)\right)$. With in-plane and out-plane polarized laser beams combined, a double-well lattice can be created (figure 15(b)).

Ground state can be achieved by adiabatically loading atoms into the lattice. For the double-well lattice, different from simple Bravais lattices, the band gap could be very small compared with the energy scale $\hslash T_{\text{load}}^{-1}$, with T_load the loading time. Then the Landau–Zenner transitions across the lowest and first excited bands can be significant. The population of the first excited band causes the oscillations in the momentum distribution measured in time-of-flight, which are observed in experiments.

The relation between the observed oscillations in the momentum distribution and the population of the excited band can be quantified by constructing a two-band model,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} H=\underset{\mathbf{r},{{\mathbf{r}}^{\prime}}}{\sum}\,\phi _{\mathbf{r}}^{\dagger}{{T}_{\mathbf{r}{{\mathbf{r}}^{\prime}}}}{{\phi}_{{{\mathbf{r}}^{\prime}}}}, \end{array}\end{eqnarray} \tag{ 103 }$$

with ${{\phi}_{\mathbf{r}}}={{\left[{{\phi}_{A,\mathbf{r}}},{{\phi}_{\text{B},\mathbf{r}}}\right]}^{T}}$ where ${{\phi}_{A}}$ and ${{\phi}_{\text{B}}}$ are annihilation operators for the localized orbitals, ${{w}_{A}}\left(\mathbf{x}-\mathbf{r}\right)$ and ${{w}_{\text{B}}}\left(\mathbf{x}-\mathbf{r}\right)$, in the two sub-wells at site $\mathbf{r}$ in the double-well lattice. In momentum space, the Hamiltonian then reads $H={\sum}_{\mathbf{k}}{{\phi}^{\dagger}}\left(\mathbf{k}\right)\mathcal{H}\left(\mathbf{k}\right)\phi \left(\mathbf{k}\right),$ with $\phi \left(\mathbf{k}\right)$ Fourier transform of ${{\phi}_{\mathbf{r}}}$. After loading bosonic atoms into the lattice, the condensate is a superposition of the ground state and excited state at lattice momentum $\mathbf{k}=0$,

$$\begin{eqnarray}&&|\psi \rangle ={{\psi}_{g}}|g\rangle +{{\psi}_{e}}|e\rangle.\end{eqnarray}$$

Writing $\mathcal{H}\left(\mathbf{0}\right)$ as

$$\begin{eqnarray}&&\mathcal{H}\left(\mathbf{0}\right)={{h}_{0}}{\mathbb 1}+{{h}_{x}}{{\sigma}_{x}}+{{h}_{y}}{{\sigma}_{y}},\end{eqnarray}$$

the dynamics of the state $|\psi \rangle$ is given as $|\psi (t)\rangle ={{\psi}_{g}}{{\text{e}}^{\text{i} \Delta t/2}}|g\rangle +$ ${{\psi}_{e}}{{\text{e}}^{-\text{i} \Delta t/2}}|e\rangle$, with $\Delta =2\sqrt{h_{x}^{2}+h_{z}^{2}}$. In terms of ${{\phi}_{A,B}}\left(\mathbf{k}\right)$ basis, we have

$$\begin{eqnarray}\begin{array}{*{35}{l}} |\psi (t)\rangle = & \left\{\left[{{\psi}_{e}}{{\text{e}}^{-\text{i} \Delta t/2}}\cos (\gamma /2)-{{\psi}_{g}}{{\text{e}}^{\text{i} \Delta t/2}}\sin (\gamma /2)\right]\phi _{A}^{\dagger}(0) \right. \\ {} & \left. +\,\left[{{\psi}_{e}}{{\text{e}}^{-\text{i} \Delta t/2}}\sin (\gamma /2)+{{\psi}_{g}}{{\text{e}}^{\text{i} \Delta t/2}}\cos (\gamma /2)\right]\phi _{\text{B}}^{\dagger}(0)\right\}|0\rangle, \end{array}\end{eqnarray}$$

with γ the polar angle of the vector $\left({{h}_{z}},{{h}_{x}}\right)$. The momentum distribution is then given as

$$\begin{eqnarray}\begin{array}{*{35}{l}} n\left(\mathbf{k}\right)= & \text{const}.+2\text{Re}\left[\psi _{e}^{\ast}{{\psi}_{g}}{{\text{e}}^{\text{i} \Delta t}}\left(\tilde{{w}}\,_{A}^{\ast}\left(\mathbf{k}\right)\cos (\theta /2)+\tilde{{w}}\,_{\text{B}}^{\ast}\left(\mathbf{k}\right)\sin (\theta /2)\right) \right. \\ {} & \left. \times \left({{\tilde{{w}}\,}_{\text{B}}}\left(\mathbf{k}\right)\cos (\theta /2)+{{\tilde{{w}}\,}_{A}}\left(\mathbf{k}\right)\sin (\theta /2)\right)\right] \end{array}\end{eqnarray} \tag{ 104 }$$

where ${{\tilde{w}}_{A,B}}\left(\mathbf{k}\right)$ is the Fourier transform of ${\sum}_{\mathbf{r}}{{w}_{A,B}}\left(\mathbf{x}-\mathbf{r}\right)$. The population fraction of the excited band could thus be extracted from the dynamical evolution of momentum distribution.

Although the above discussions were restricted to the setup in the experiment (Sebby-Strabley et al 2006), the coherent oscillation in time-of-flight is a generic phenomenon when a superposed state of ground and excited bands is prepared. And indeed similar oscillations are observed in other double-well lattices as well (Anderlini et al 2007, Müller et al 2007, Trotzky et al 2008).

Early theoretical studies of p-band condensates focus on the case with the point group D₄ symmetry. For the lattice potential realized in the experiment of Hamburg (equation (99)), the point group symmetry is maintained only for the ideal case $\epsilon =1$ and $\alpha =0$, where the potential reduces to $V=-{{V}_{0}}\left({{\eta}^{2}}{{\cos}^{2}}kx+{{\cos}^{2}}ky+2\eta \cos \theta \cos kx\cos ky\right)$. For the realistic situation with $\epsilon <1$, the D₄ symmetry is thus broken and only reflection symmetry with respect to the x-axis is preserved. The asymmetry could be partially compensated by setting $\alpha ={{\alpha}_{\text{iso}}}$, for which the potential reads $V=-{{V}_{0}}\epsilon \left[{{\eta}^{2}}\epsilon {{\cos}^{2}}kx+{{\cos}^{2}}ky\right]-{{V}_{0}}\epsilon \eta \cos kx\left[\cos (ky+\theta )+\right.$ $\left. {{\epsilon}^{2}}\cos (ky-\theta )\right]$. The consequences of asymmetry are studied in detail in Cai and Wu (2011) and Shchesnovich (2012).

The band structure is calculated by plane-wave expansion (Cai and Wu 2011). The reciprocal lattice lattice vectors are defined as ${{\mathbf{G}}_{m,n}}=m{{\mathbf{b}}_{1}}+n{{\mathbf{b}}_{2}}$, with ${{\mathbf{b}}_{1,2}}=(\pm \pi /a,\pi /a)$ (a the lattice constant). Taking the single-particle Hamiltonian ${{H}_{0}}=-{{\hslash}^{2}}{{\vec{\nabla}}^{2}}/(2M)+V\left(\mathbf{x}\right)$, the diagonal matrix elements are $\langle \mathbf{k}+{{\mathbf{G}}_{m,n}}|{{H}_{0}}|\mathbf{k}+{{\mathbf{G}}_{m,n}}\rangle =$ ${{E}_{r}}\left\{{{\left[a{{k}_{x}}/\pi +(m-n)\right]}^{2}}+{{\left[a{{k}_{y}}/\pi +(m+n)\right]}^{2}}\right\},$ with E_r the single-photon recoil energy, and the off-diagonal matrix elements are

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {} & \langle \mathbf{k}|{{H}_{0}}|\mathbf{k}+{{\mathbf{G}}_{\pm 1,0}}\rangle =-\frac{{{V}_{0}}}{4}\eta \epsilon \left(\cos \alpha {{\text{e}}^{\mp \text{i}\theta}}+{{\text{e}}^{\pm \text{i}\theta}}\right), \\ {} & {} & \langle \mathbf{k}|{{H}_{0}}|\mathbf{k}+{{\mathbf{G}}_{0,\pm 1}}\rangle =-\frac{{{V}_{0}}}{4}\eta \left(\cos (\alpha ){{\text{e}}^{\pm \text{i}\theta}}+{{\epsilon}^{2}}{{\text{e}}^{\mp \text{i}\theta}}\right), \\ {} & {} & \langle \mathbf{k}|{{H}_{0}}|\mathbf{k}+{{\mathbf{G}}_{\pm 1,\mp 1}}\rangle =-\frac{{{V}_{0}}}{4}\epsilon {{\eta}^{2}}\cos \alpha, \\ {} & {} & \langle \mathbf{k}|{{H}_{0}}|\mathbf{k}+{{\mathbf{G}}_{\pm 1,\pm 1}}\rangle =-\frac{{{V}_{0}}}{4}\epsilon \cos \alpha. \end{array}\end{eqnarray} \tag{ 105 }$$

There are four time-reversal invariant points in the Brillouin zone, O  =  (0,0), ${{X}_{\pm }}=\left(\pm \frac{\pi}{2a},\frac{\pi}{2a}\right)$, and $M=\left(\frac{\pi}{a},\frac{\pi}{a}\right)$, at which the Bloch functions are real valued. The band spectra are symmetric at these points, and consequently ${{\partial}_{\mathbf{k}}}\varepsilon \left(\mathbf{k}\right)=0$, which means that they are saddle points in the band structure. For the choice $\alpha ={{\alpha}_{\text{iso}}}$, the second band has double degenerate minima at X₊ and X₋. For $\alpha <{{\alpha}_{\text{iso}}}$ ($\alpha >{{\alpha}_{\text{iso}}}$), X₊ (X₋) becomes the unique band minimum.

To investigate the interaction effects, the Gross–Pitaevskii equation

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \left\{-\frac{{{\hslash}^{2}}{{{\vec{\nabla}}}^{2}}}{2M}+{{V}_{\text{eff}}}\left(\mathbf{x}\right)\right\} \Psi \left(\mathbf{x}\right)=E \Psi \left(\mathbf{x}\right), \end{array}\end{eqnarray} \tag{ 106 }$$

with ${{V}_{\text{eff}}}\left(\mathbf{x}\right)=V\left(\mathbf{x}\right)+g\rho | \Psi \left(\mathbf{x}\right){{|}^{2}}$, is solved self-consistently by assuming the condensate wavefunction is a superposition of Bloch functions at ${{X}_{\pm }}$,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \Psi \left(\mathbf{x}\right)=\cos (\delta ){{\psi}_{{{X}_{+}}}}\left(\mathbf{x}\right)+\sin (\delta ){{\text{e}}^{\text{i}\phi}}{{\psi}_{{{X}_{-}}}}\left(\mathbf{x}\right). \end{array}\end{eqnarray} \tag{ 107 }$$

The Bloch functions ${{\psi}_{{{X}_{\pm }}}}$ have nodal lines in space, while the variational condensate wavefunction could avoid nodal lines by having complex values (with $\delta \ne 0$ or $\pi /2$, and $\phi \ne 0$). The complex solution is spatially more uniform and thus more favorable by interactions, but at the same time costs more kinetic energy when $\alpha \ne {{\alpha}_{\text{iso}}}$.

The competition between interactions and anisotropy leads to an interesting phase diagram containing two real and one complex states of Bose–Einstein condensation. The Gross–Pitaevskii approach finds second order transitions at zero temperature (Cai and Wu 2011). The phase transitions can be understood within a Ginzburg–Landau theory,

$$\begin{eqnarray}\begin{array}{*{35}{l}} F= & -{{r}_{1}}|{{\psi}_{+}}{{|}^{2}}-{{r}_{2}}|{{\psi}_{-}}{{|}^{2}}+\,{{g}_{1}}|{{\psi}_{+}}{{|}^{4}}+\,{{g}_{2}}|{{\psi}_{-}}{{|}^{4}} \\ {} & +\,{{g}_{3}}|{{\psi}_{+}}{{|}^{2}}|{{\psi}_{-}}{{|}^{2}}+\,{{g}_{4}}\left(\psi _{+}^{\ast 2}\psi _{-}^{2}+\text{c}.\text{c}.\right), \end{array}\end{eqnarray} \tag{ 108 }$$

with ${{\psi}_{\pm }}$ describes the condensate component at ${{X}_{\pm }}$. The Umklapp term g₄  >  0 favors the complex state. Assuming ${{r}_{1}},{{r}_{2}}$, and ${{g}_{3}}-2{{g}_{4}}$  >0, the complex state occurs in the regime

$$\begin{eqnarray}&&\frac{{{g}_{3}}-2{{g}_{4}}}{2{{g}_{2}}}<\frac{{{r}_{1}}}{{{r}_{2}}}<\frac{2{{g}_{1}}}{{{g}_{3}}-2{{g}_{4}}}.\end{eqnarray} \tag{ 109 }$$

The predicted phase diagram is confirmed in the experiment (Ölschläger et al 2013).

Direct measurement of orbital ordering, namely the staggered angular momentum, was thought to be an experimental challenge, which motivates a theoretical proposal of using quench dynamics (Li et al 2016). The key idea could be understood by drawing an analogy between the two orbital states at each site $\left(\,{{p}_{x}},{{p}_{y}}\right)$, and a pseudospin-1/2 degrees of freedom ($\uparrow$, ↓). In this analogy, the ${{p}_{x}}\pm \text{i}{{p}_{y}}$ state corresponds to a pseudospin pointing along the y direction in spin space. Applying a 'magnetic field' along the x direction to this pseudospin should then induce Larmor precession, leading to periodic oscillations of the z-magnetization, corresponding to the population imbalance between two p-orbitals, $\Delta N=N\left(\,{{p}_{x}}\right)-N\left(\,{{p}_{y}}\right)$. Here we consider a square lattice. We can take a certain initial state and then quickly turn on a strong 'magnetic field'

$$\begin{eqnarray}&&{{H}_{\text{mag}}}=\underset{\mathbf{r}}{\sum}\,{{(-1)}^{{{r}_{x}}+{{r}_{y}}}}\lambda \left(\mathbf{r}\right)\left[b_{x}^{\dagger}\left(\mathbf{r}\right){{b}_{y}}\left(\mathbf{r}\right)+\text{h}.\text{c}.\right]\end{eqnarray} \tag{ 110 }$$

at time $\tau =0$. For simplicity, the 'magnetic field' is assumed to be strong enough to completely dominate the short-time dynamics. If initially a staggered superposition ${{p}_{x}}\pm {{\text{e}}^{\text{i}\theta}}{{p}_{y}}$ is prepared, all local Larmor precessions add up to produce a macroscopic oscillation in the orbital imbalance $\Delta N$. This imbalance evolves within a Heisenberg picture as

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \frac{\text{d} \Delta N\left(\mathbf{r},\tau \right)}{\text{d}\tau}=-\text{i}\left[\Delta N\left(\mathbf{r},\tau \right),{{H}_{\text{mag}}}\right]=-2\lambda \left(\mathbf{r}\right)L_{z}^{\text{stag}}\left(\mathbf{r},\tau \right), \end{array}\end{eqnarray} \tag{ 111 }$$

with $L_{z}^{\text{stag}}$ the staggered angular momentum operator, whose time evolution is described by

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \frac{\text{d}L_{z}^{\text{stag}}}{\text{d}\tau}=2\lambda \left(\mathbf{r}\right) \Delta N\left(\mathbf{r},\tau \right). \end{array}\end{eqnarray} \tag{ 112 }$$

This leads to oscillations in $\langle \Delta N\left(\mathbf{r},\tau \right)\rangle$,

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {} & \langle \Delta N\left(\mathbf{r},\tau \right)\rangle \\ {} & = & \langle \Delta N\left(\mathbf{r},0\right)\rangle \cos \left(2\lambda \left(\mathbf{r}\right)\tau \right)-\langle L_{z}^{\text{stag}}\left(\mathbf{r},0\right)\rangle \sin \left(2\lambda \left(\mathbf{r}\right)\tau \right) \\ {} & \equiv & A\left(\mathbf{r}\right)\cos \left(2\lambda \left(\mathbf{r}\right)\tau +\phi \left(\mathbf{r}\right)\right), \end{array}\end{eqnarray} \tag{ 113 }$$

where $\langle \Delta N\left(\mathbf{r},0\right)\rangle$ and $\langle L_{z}^{\text{stag}}\left(\mathbf{r},0\right)\rangle$ denote the orbital imbalance and staggered angular momentum for the initial state. The trigonometric form of this time-dependent equation thus defines the quantities $A\left(\mathbf{r}\right)$ and $\phi \left(\mathbf{r}\right)$, ready to compare with the experimental measurement of $\Delta N$.

Neglecting spatial inhomogeneity in $\lambda \left(\mathbf{r}\right)$ and $\phi \left(\mathbf{r}\right)$, we can set $\lambda \left(\mathbf{r}\right)=\lambda$ and $\phi \left(\mathbf{r}\right)=\phi$, and extract the initial angular momentum order from the amplitude A and the phase shift ϕ in the dynamics of the spatially averaged orbital imbalance $\overline{\langle \Delta N(\tau )\rangle}=1/{{N}_{s}}{\sum}_{\mathbf{r}} \Delta N\left(\mathbf{r},\tau \right)$. The coefficient λ can be read off from the oscillation period ${{\tau}_{Q}}\equiv \pi /\lambda$. The orbital population imbalance can be measured directly in time-of-flight experiments.

For a C₄ symmetric initial state with non-zero staggered angular momentum, but no orbital imbalance, $\overline{\langle \Delta N(\tau )\rangle}$ is expected to oscillate with a non-zero amplitude and phase shift $\phi =\pm \pi /2$ whose sign will fluctuate from realization to realization. By contrast $\overline{\langle \Delta N(\tau )\rangle}=0$ should be observed for a completely disordered state. The amplitude of the signal is thus a direct measure of the staggered angular momentum order parameter.

In the case that C₄ symmetry is explicitly broken as achieved in the recent experiment, a state with an initial orbital imbalance but no angular momentum order would exhibit oscillations with a finite amplitude but no phase shift, i.e. $\phi =0$. In contrast, for a state with angular momentum order, the spontaneous time-reversal symmetry breaking yields a finite phase shift $\phi \ne 0$, which would vanish in a singular fashion as we tune from the angular momentum ordered to disordered regime through a second order phase transition.

The required coupling H_mag can be engineered by adding a quench potential ${{V}_{\text{mag}}}\left(\mathbf{x}\right)$ modulated in the (1,1) direction with respect to the original lattice potential. The add-on potential generates a coupling between p_x and p_y orbitals

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \epsilon \left(\mathbf{r}\right)\approx \frac{\hslash}{4m{{\omega}_{0}}}\frac{{{\partial}^{2}}{{V}_{\text{mag}}}\left(\mathbf{r}+l\left[\frac{{{\mathbf{a}}_{x}}+{{\mathbf{a}}_{y}}}{\sqrt{2}}\right]\right)}{{{a}^{2}}\partial {{l}^{2}}}{{|}_{l\to 0}}, \end{array}\end{eqnarray} \tag{ 114 }$$

where ${{\omega}_{0}}$ is the harmonic oscillator frequency of the lattice wells hosting the p-orbitals and $a=\,|{{\mathbf{a}}_{x}}|\,=\,|{{\mathbf{a}}_{y}}|$ is the lattice constant. The above estimate for the coupling strength is valid in the tight binding regime when the quench potential is weak as compared with the original optical lattice. Without loss of generality, one may consider an add-on optical potential of the form

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{V}_{\text{mag}}}\left(\mathbf{x}\right)=- \Gamma {{\cos}^{2}}\left(\frac{2\nu +1}{4}\left({{\mathbf{K}}_{x}}+{{\mathbf{K}}_{y}}\right)\centerdot \mathbf{x}\right), \end{array}\end{eqnarray} \tag{ 115 }$$

with some integer $\nu \geqslant 0$, a positive amplitude $\Gamma$, and $\left({{\mathbf{K}}_{x}}\right.$, $\left. {{\mathbf{K}}_{y}}\right)$ denoting the primitive vectors of the reciprocal lattice. This potential leads to a p_x/p_y coupling

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \epsilon \left(\mathbf{r}\right)=\frac{{{E}_{\text{r}}}}{\hslash {{\omega}_{0}}}\frac{ \Gamma }{4}{{(2\nu +1)}^{2}}{{(-1)}^{{{r}_{x}}+{{r}_{y}}}}, \end{array}\end{eqnarray} \tag{ 116 }$$

with E_r the photon recoil energy with wave number $1/2|{{\mathbf{K}}_{x}}+{{\mathbf{K}}_{y}}|$. The staggering factor in the engineered coupling is crucial to probe the staggered angular momentum order.

This quench proposal brings other interesting possibilities in addition to providing a method to probe orbital order. For instance, one can simulate spin dynamics in solid state materials by studying orbital dynamics of p-band bosons. One advantage about orbital dynamics is that engineering artificial effective magnetic fields is intrinsically easier due to the the spatial nature of orbital degrees of freedom than engineering real staggered magnetic fields.

There is another proposed scheme to measure the inter-orbital phase coherence in ${{p}_{x}}\pm \text{i}{{p}_{y}}$ superfluid by Raman transition (Cai et al 2012a). In the ${{p}_{x}}\pm \text{i}{{p}_{y}}$ superfluid, condensation takes place at X₊ and X₋ and the condensate state is $| \Psi \rangle \propto {{\left(b_{{{X}_{+}}}^{\dagger}+{{\text{e}}^{\text{i}\theta}}b_{{{X}_{-}}}^{\dagger}\right)}^{N}}|0\rangle$ in general. The idea is to transform the phase coherence to number difference in momentum space. With a Raman operation, bosons in the original condensate can be transfered to a state with

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {} & b_{{{X}_{+}}}^{\prime}=\frac{1}{\sqrt{2}}\left({{b}_{{{X}_{+}}}}-\text{i}{{\text{e}}^{\text{i}\phi}}{{b}_{{{X}_{-}}}}\right), \\ {} & {} & b_{{{X}_{-}}}^{\prime}=\frac{1}{\sqrt{2}}\left({{b}_{{{X}_{-}}}}-\text{i}{{\text{e}}^{-\text{i}\phi}}{{b}_{{{X}_{+}}}}\right). \end{array}\end{eqnarray} \tag{ 117 }$$

With $\phi =0$, the phase coherence in the ${{p}_{x}}\pm \text{i}{{p}_{y}}$ state is then transformed as

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \langle \text{i}b_{{{X}_{+}}}^{\dagger}{{b}_{{{X}_{-}}}}+\text{h}.\text{c}.\rangle =\langle b_{{{X}_{-}}}^{\prime \dagger}b_{{{X}_{-}}}^{\prime}-b_{{{X}_{+}}}^{\prime \dagger}b_{{{X}_{+}}}^{\prime}\rangle \equiv \delta {{n}^{\prime}}, \end{array}\end{eqnarray} \tag{ 118 }$$

which can be extracted in time-of-flight experiment.

The required Raman transition can be implemented by two traveling-wave laser beams along different directions with corresponding wave vector ${{\mathbf{k}}_{1,2}}$ and frequency ${{\omega}_{1,2}}$ (Duan 2006). These laser beams induce an effective Raman Rabi frequency with a spatially varying phase $\Omega \left(\mathbf{x},t\right)={{ \Omega }_{0}}{{\text{e}}^{\text{i}\left(\delta \mathbf{k}\centerdot \mathbf{x}-\delta \omega t+\phi \right)}}$, where $\delta \mathbf{k}={{\mathbf{k}}_{1}}-{{\mathbf{k}}_{2}}$, $\delta \omega ={{\omega}_{1}}-{{\omega}_{2}}$, and ϕ is the relative phase between the two laser beams (see figure 16(a)). The effective Hamiltonian for the Raman process is described by

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{H}_{R}}={\int}^{}\text{d}\mathbf{x} \Omega \left(\mathbf{x},t\right){{\phi}^{\dagger}}\left(\mathbf{x}\right)\phi \left(\mathbf{x}\right)+\text{h}.\text{c}., \end{array}\end{eqnarray} \tag{ 119 }$$

where $\phi \left(\mathbf{x}\right)$ is the boson annihilation operator in continuous space. The generated spatially dependent potential couples the two condensate components at the two momentum points (in the Hamburg experiment (Wirth et al 2011) ${{X}_{\pm }}=(\pm \pi /2,\pi /2)$, requiring $\delta \mathbf{k}={{X}_{+}}-{{X}_{-}}=(\pi,0)$).

Zoom In Zoom Out Reset image size

To avoid complications of interband transitions (with band gap $\Delta$) and dynamics caused by tunnelings (t), an optimal choice for the Raman coupling strength is $t\ll \hslash {{ \Omega }_{0}}\ll \Delta$. For the experimental situation, the Raman coupling strength should be chosen to be ${{ \Omega }_{0}}\approx 2\pi \times 0.5$ kHz. Thus the required duration of the Raman pulse is around 1 ms. To get efficient Raman operation, the frequency $\delta \omega$ should match the energy difference between the initial and final states which is around a few Hz. Therefore the phase accumulation $\delta \omega t$ within the duration of Raman pulse is negligible. With this approximation the Raman coupling is simplified to be

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{H}_{R}}\approx \underset{\mathbf{k}}{\sum}\,{{\text{e}}^{\text{i}\phi}}\lambda \left(\mathbf{k}\right)b_{\mathbf{k}+\delta \mathbf{k}}^{\dagger}{{b}_{\mathbf{k}}}+\text{h}.\text{c}. \end{array}\end{eqnarray} \tag{ 120 }$$

Here $\lambda \left(\mathbf{k}\right)$ is the $\mathbf{k}$ dependent effective coupling, which can be calculated from the Bloch functions. For the Hamburgh experiment, it is estimated that $\lambda \left({{X}_{\pm }}\right)\approx 0.98{{ \Omega }_{0}}\equiv \lambda$. Choosing the duration of the Raman pulse to be $\lambda \delta t=\pi /4$, the required state transfer in equation (117) is achieved. The resultant density difference is

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \delta {{n}^{\prime}}=\langle \text{i}{{\text{e}}^{\text{i}\phi}}b_{{{X}_{+}}}^{\dagger}{{b}_{{{X}_{-}}}}+\text{h}.\text{c}.\rangle. \end{array}\end{eqnarray} \tag{ 121 }$$

For the ${{p}_{x}}\pm \text{i}{{p}_{y}}$ superfluid, the density difference would be $\delta {{n}^{\prime}}\propto \cos (\phi )$. With $\phi =0$, $\delta {{n}^{\prime}}=\langle \text{i}b_{{{X}_{+}}}^{\dagger}{{b}_{{{X}_{-}}}}+\text{h}.\text{c}.\rangle$ represents the order parameter of the complex orbital ordering (figure 16).

In a recent experiment (Kock et al 2015), that generalizes the idea of Young's double slits, an interference measurement has been implemented to detect the inter-orbital phase coherence in the ${{p}_{x}}+\text{i}{{p}_{y}}$ superfluid. In this experiment, two independent copies of the lattice condensates are prepared with the experimental setup as illustrated in figure 17. The condensates are simultaneously prepared in the second band in two spatially separated regions of the lattice. After the state preparation, all potentials are switched off. The zeroth-order Bragg peaks observed in the xy-plane carry interference patterns in the z direction due to overlapping contributions from the condensates originally separate in space. In the simplified picture approximating the two condensates by two point sources, the wave length of the density grating in the interference is ${{\lambda}_{z}}=\frac{2\pi \hslash {{t}_{\text{TOF}}}}{m{{d}_{z}}}$, with t_TOF the time of ballistic expansion, d_z the spatial separation of the two condensates. This estimate is quantitatively consistent with experimental results.

Zoom In Zoom Out Reset image size

In the ballistic expansion, the Bragg peaks (labeled by 1, 2, 3 and 4 in figure 17) yield the Fourier components of the condensate wavefunction, and we can associate a phase for each component, ${{\theta}_{j=1,2,3,4}}$. Since the spatially separate condensates are decoupled, they carry different phases, ${{\theta}_{j}}$ and $\theta _{j}^{\prime}$. From the relative phase $\Delta {{\theta}_{j}}={{\theta}_{j}}-\theta _{j}^{\prime}$, we can introduce $\Delta {{\theta}_{i,j}}= \Delta {{\theta}_{i}}- \Delta {{\theta}_{j}}$, which directly determines the correlation among the interference patterns in the Bragg peaks. If $\Delta {{\theta}_{i,j}}=0$ (π), the density patterns of the ith and jth peak are positively (negatively) correlated. The interference patterns obtained in experiments yield that $\Delta {{\theta}_{1,3}}= \Delta {{\theta}_{2,4}}=0$, over 420 independent realizations, and that $\Delta {{\theta}_{1,2}}= \Delta {{\theta}_{1,4}}= \Delta {{\theta}_{2,3}}= \Delta {{\theta}_{3,4}}$, and their value spontaneously chooses 0 or π. The interference measurement unambiguously tell that the phase of different momentum components is indeed correlated. To the best of our knowledge, the experiment (Kock et al 2015) appears to be the first phase sensitive measurement which poses an important constraint on the nature of p-orbital Bose–Einstein condensates. It is desirable that future experiments can directly probe the phase lock between the condensate components at two band minima, corresponding to the X₊ and X₋ points in the paper (Kock et al 2015).

The crystal structure of the atomic ions in solids provide confining potential for electrons due to strong Coulomb force. Electrons in solids are usually nearly localized on atomic ions and the resulting orbital wavefunctions (or the shape of the electron cloud) are determined by the strong confining potential. This orbital degree of freedom is of great importance in correlated materials such as transition metal oxides (Tokura and Nagaosa 2000). Many intriguing phenomena such as metal-insulator transitions and colossal magnetoresistance can be attributed (or partially attributed) to the interplay of d-orbitals with charge and spin degrees of freedom.

Considering a transition-metal oxide material with perovskite crystal structure, d-orbital electrons localized on the transition-metal atom are surrounded by six oxygen ions O²⁻, which give rise to crystal field and consequent energy splitting of the d-orbitals. Orbital wavefunctions pointing towards the negative-charged oxygen ions (the e_g orbitals, ${{d}_{{{x}^{2}}-{{y}^{2}}}}$ and ${{d}_{3{{z}^{2}}-{{r}^{2}}}}$) have higher energy compared with those pointing in other orientations (the t_2d orbitals, d_xy, d_yz and d_xz) due to Coulomb repulsion (see figure 18). The spatial nature of orbital makes it intrinsically attached to the crystal fields, even in the absence of the relativistic spin–orbital interaction, and this intrinsic coupling of orbital degree of freedom to crystal fields and the resultant crystal symmetry make it distinct from real spins. When orbitals are modeled as pseudo-spins, the model Hamiltonian is in general lack of SU(2) symmetry. Consider a typical Mott insulator LaMnO₃ as an example. A neutral Mn atom has an electron configuration $3{{d}^{5}}4{{s}^{2}}$. Losing three electrons, Mn³⁺ in this material has four electrons in those five d-orbitals. From Hund's rule, the spins are aligned ferromagnetically, and there are thus two possibilities for e_g orbitals with either ${{d}_{{{x}^{2}}-{{y}^{2}}}}$ or ${{d}_{3{{z}^{2}}-{{r}^{2}}}}$ being occupied. This represents the orbital degree of freedom in this Mott insulator, which can be modeled as pseudo-spins T_x,y,z. The model Hamiltonian is

$$\begin{eqnarray}&&H=\underset{\mathbf{r}{{\mathbf{r}}^{\prime}}}{\sum}\,J_{\mathbf{r}{{\mathbf{r}}^{\prime}}}^{\alpha \beta}{{T}_{\alpha}}\left(\mathbf{r}\right){{T}_{\beta}}\left({{\mathbf{r}}^{\prime}}\right),\end{eqnarray}$$

which is typically not SU(2) symmetric. With a long range orbital order, spin magnetism would be strongly affected by so called Jahn-Teller effect (Jahn and Teller 1937).

Zoom In Zoom Out Reset image size

Most p-orbital solid state materials, for example the semiconducting silicon and graphene, are actually weakly correlated. However, recent studies in one oxide heterostructure LAO/STO have found that correlated physics such as ferromagnetism emerges from the effective p-orbitals, where p_x and p_y are mimicked by d_xz and d_yz orbitals (the degeneracy with d_xy orbital is broken due to lack of out-of-plane inversion symmetry at the interface). In d-orbital systems, correlated physics usually emerges due to large Hubbard U interaction because of the tight confinement of these orbitals. The emergence of correlated physics in p-orbital systems on the other hand could be attributed to a different origin, which is the quasi-one dimensionality (Chen and Balents 2013, Li et al 2014). In one dimension at low filling, the magnetic susceptibility diverges as ${{\chi}_{1d}}\sim 1/{{\rho}^{2}}$, where ρ is the occupation number per site. Even for infinitesimal interaction U, there is a strong interaction effect: the ratio of the interacting to free fermion susceptibility diverges, ${{\chi}_{1d}}/{{\chi}_{ff}}\to \infty$ for $\rho \to 0$. A general result for the free energy (per site) versus magnetization at low density is obtained to be

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} F=2\rho \,{{J}_{\text{eff}}}{{F}_{1}}\left[\frac{M}{\rho},\frac{{{k}_{\text{B}}}T}{{{J}_{\text{eff}}}}\right]-{{J}_{\text{H}}}{{M}^{2}}, \end{array}\end{eqnarray} \tag{ 122 }$$

where M is the magnetization (per site), J_H is the Hund's rule coupling, J_eff is the effective antiferromagnetic coupling, and F₁[m,t] is the free energy per site of the one-dimensional antiferromagnetic chain, with reduced magnetization m and temperature t (this is known from thermodynamic Bethe ansatz). The effective coupling J_eff is reasonably conjectured to scale as ${{J}_{\text{eff}}}\propto {{\rho}^{3}}$ (Chen and Balents 2013). From the free energy, the Hund's energy is dominant and favors a ferromagnetic state with sufficiently low density for arbitrarily weak Hund's coupling J_H. A rigorous work (Li et al 2014) studies the higher filling regime (but assumes no double occupancy), where a ferromagnetic ground state for p-orbital fermions is proved based on transitivity and non-positivity of the many-body Hamiltonian. Further studies are required to find out the boundary of the ferromagnetism in p-orbital fermions.

In material science, design of materials for applications is an important subject. Recent developments involve engineering heterostructures with hybrid materials. For example oxide heterostructures such as LaAlO₃/SrTiO₃ and GdTiO₃/SrTiO₃ have been created and extensively studied. While the properties of many materials can be calculated within the density functional theory (DFT), this approach fails for ones of strong correlation for which d-orbital electrons typically play an important role. At the same time, these strongly correlated materials could have fascinating properties including important applications. High T_c superconductivity belongs to this class. Lots of efforts have been made in searching for materials with higher T_c, but there is no real improvement in the last two decades. Lack of reliable tools in predicting T_c leaves the design of high T_c superconductivity essentially to empirical trials, which are costly in both time and materials. Developing new tools to simulate strongly correlated materials by incorporating correlation effects in DFT has triggered tremendous interest but appears to be very challenging.

To address the challenge of simulating correlated d-orbital electrons in classical computers, one alternative way is to create synthetic orbital matter with optical lattices and take it as a quantum orbital simulator. With this optical-lattice-based quantum orbital simulator, the ultimate procedure for material design would be—(1) conceive a particular design of materials; (2) determine the orbital configuration of the imagined material by quantum chemistry; and (3) apply cold atoms in optical lattices to simulate the properties. In such a way, we could explore the imagined quantum materials for desired properties, bypassing the often tedious chemical process of really fabricating them from electronic compounds. This would significantly speed up the material design and should help improve key quantities of great interest, for instance, the value of critical temperature T_c of superconductivity in future. Although the optical lattice experiment is still at a very early stage, with future developments, synthetic orbital matter in optical lattices could be extremely helpful to the design of real materials.

Finally, we would like to point out that orbital degrees of freedom are found to play an important role for a vast majority of intriguing electronic quantum materials that condensed matter physicists have found since 1970s. Magnetic materials of spin only are an important class of systems that have been studied with great progress and remain to pose new challenges, such as frustrated magnets possibly showing spin liquid phases. In fact, the spin-only systems represent a small fraction of the world of real materials. Furthermore, past theoretical studies predicted exotic phenomena for model systems that have no spin but only orbital degrees of freedom. Such hypothetical models, which previously might have seemed too special and excessive, now become readily realizable with optical lattices. On this regard, using higher orbital bands of the optical lattice appears to open up a new front to explore orbital physics, both for understanding the electronic systems and for exploring artificial quantum orbital-only models that have no prior analogue in solids.

Coherent dynamics across different bands has been observed in many experiments (Jona-Lasinio et al 2003, Sebby-Strabley et al 2006, Anderlini et al 2007, Cheinet et al 2008, Trotzky et al 2008, Zhai et al 2013, Hu et al 2015). In particular the recent experiments (Zhai et al 2013, Hu et al 2015) have demonstrated fast coherent controllability of orbital degrees of freedom. These experimental developments open up possibilities of studying many-body dynamics of high orbital, where the observed Rabi-like oscillations between different bands can be affected by interaction. One particular example would be orbital Josephson effect, which has been studied for double-well potentials (García-March et al 2011, Garcia-March et al 2012, Gillet et al 2014, Garcia-March and Carr 2015). This effect has also been seen in numerical simulations of a dynamical procedure, proposed to detect the $p+\text{i}p$ BEC (Cai et al 2012a, Li et al 2016).

The orbital Josephson effect is expected to be generic for various experimental setups for high orbital atoms. Here we consider the specific setup proposed to probe the complex order (see section 4.5). Assuming all atoms condense, the dynamics is then approximately captured by a two-mode Hamiltonian,

$$\begin{eqnarray}\begin{array}{*{35}{l}} H= & \lambda b_{{{\mathbf{K}}_{1}}}^{\dagger}{{b}_{{{\mathbf{K}}_{2}}}}+\text{h}.\text{c}.+{{g}_{1}}\left(b_{{{\mathbf{K}}_{1}}}^{\dagger}b_{{{\mathbf{K}}_{1}}}^{\dagger}{{b}_{{{\mathbf{K}}_{1}}}}{{b}_{{{\mathbf{K}}_{1}}}}+{{\mathbf{K}}_{1}}\to {{\mathbf{K}}_{2}}\right) \\ {} & +\,{{g}_{2}}b_{{{\mathbf{K}}_{1}}}^{\dagger}{{b}_{{{\mathbf{K}}_{1}}}}b_{{{\mathbf{K}}_{2}}}^{\dagger}{{b}_{{{\mathbf{K}}_{2}}}}+\,{{g}_{3}}\left(b_{{{\mathbf{K}}_{1}}}^{\dagger}b_{{{\mathbf{K}}_{1}}}^{\dagger}{{b}_{{{\mathbf{K}}_{2}}}}{{b}_{{{\mathbf{K}}_{2}}}}+\text{h}.\text{c}.\right), \end{array}\end{eqnarray} \tag{ 123 }$$

where ${{b}_{{{\mathbf{K}}_{1,2}}}}$ are the two condensed modes and the last term g₃ is a Umklapp process. Following the treatment of Josephson effect developed for double-well Bose–Einstein condensates (Smerzi et al 1997, Zapata et al 1998), the dynamical state could be approximated by

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} | \Psi (t)\rangle =\frac{1}{\sqrt{N!}}{{\left({{\psi}_{1}}(t)b_{{{\mathbf{K}}_{1}}}^{\dagger}+{{\psi}_{2}}(t)b_{{{\mathbf{K}}_{2}}}^{\dagger}\right)}^{N}}|0\rangle. \end{array}\end{eqnarray} \tag{ 124 }$$

The corresponding time-dependent Gross–Pitaevskii equation is (Cai et al 2012a)

$$\begin{eqnarray}\begin{array}{*{35}{l}} \text{i}{{\partial}_{t}}{{\psi}_{1}}(t) & = & \lambda {{\psi}_{2}}(t)+\left(2{{g}_{1}}|{{\psi}_{1}}{{|}^{2}}+\,{{g}_{2}}|{{\psi}_{2}}{{|}^{2}}\right){{\psi}_{1}}+2{{g}_{3}}\psi _{1}^{\ast}\psi _{2}^{2}, \\ \text{i}{{\partial}_{t}}{{\psi}_{2}}(t) & = & \lambda {{\psi}_{1}}(t)+\left(2{{g}_{1}}|{{\psi}_{2}}{{|}^{2}}+\,{{g}_{2}}|{{\psi}_{1}}{{|}^{2}}\right){{\psi}_{2}}+2{{g}_{3}}\psi _{2}^{\ast}\psi _{1}^{2}. \end{array}\end{eqnarray} \tag{ 125 }$$

To make the dynamics more physical, one can rewrite the wavefunctions ${{\psi}_{j}}(t)$ in terms of densities and phases as

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {} & {{\psi}_{1}}\to \sqrt{{{\rho}_{1}}}{{\text{e}}^{\text{i}{{\theta}_{1}}}}, \\ {} & {} & {{\psi}_{2}}\to \sqrt{{{\rho}_{2}}}{{\text{e}}^{\text{i}{{\theta}_{2}}}}. \end{array}\end{eqnarray}$$

The equation of motion is most easily derived by constructing the Lagrangian, which takes the form,

$$\begin{eqnarray}\begin{array}{*{35}{l}} L= & -{{\rho}_{1}}{{\partial}_{t}}{{\theta}_{1}}-{{\rho}_{2}}{{\partial}_{t}}{{\theta}_{2}}-\left\{2\lambda \sqrt{{{\rho}_{1}}{{\rho}_{2}}}\cos \left({{\theta}_{2}}-{{\theta}_{1}}\right)+2{{g}_{3}}{{\rho}_{1}}{{\rho}_{2}}\cos \left(2\left({{\theta}_{2}}-{{\theta}_{1}}\right)\right) \right. \\ {} & \left. +\,{{g}_{1}}\left(\rho _{1}^{2}+\rho _{2}^{2}\right)+{{g}_{2}}{{\rho}_{1}}{{\rho}_{2}}\right\}. \end{array}\end{eqnarray} \tag{ 126 }$$

From Euler–Lagrangian equations,

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {} & {{\partial}_{t}}{{\rho}_{j}}=-\frac{\partial L}{\partial {{\theta}_{j}}}, \\ {} & {} & \frac{\partial L}{\partial {{\rho}_{j}}}=0, \end{array}\end{eqnarray} \tag{ 127 }$$

one gets

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\partial}_{t}}{{\rho}_{1}}=-{{\partial}_{t}}{{\rho}_{2}}=2\lambda \sqrt{{{\rho}_{1}}{{\rho}_{2}}}\sin \left({{\theta}_{2}}-{{\theta}_{1}}\right)+4{{g}_{3}}{{\rho}_{1}}{{\rho}_{2}}\sin \left(2\left({{\theta}_{2}}-{{\theta}_{1}}\right)\right), \\ {{\partial}_{t}}{{\theta}_{1}}=-\lambda \sqrt{\frac{{{\rho}_{2}}}{{{\rho}_{1}}}}\cos \left({{\theta}_{2}}-{{\theta}_{1}}\right)-2{{g}_{3}}{{\rho}_{2}}\cos \left(2\left({{\theta}_{2}}-{{\theta}_{1}}\right)\right)-2{{g}_{1}}{{\rho}_{1}}-{{g}_{2}}{{\rho}_{2}}, \\ {{\partial}_{t}}{{\theta}_{2}}=-\lambda \sqrt{\frac{{{\rho}_{1}}}{{{\rho}_{2}}}}\cos \left({{\theta}_{1}}-{{\theta}_{2}}\right)-2{{g}_{3}}{{\rho}_{1}}\cos \left(2\left({{\theta}_{1}}-{{\theta}_{2}}\right)\right)-2{{g}_{1}}{{\rho}_{2}}-{{g}_{2}}{{\rho}_{1}}. \end{array}\end{eqnarray}$$

To make a direct connection to Josephson effects, the number imbalance and phase difference are defined to be $z={{\rho}_{1}}-{{\rho}_{2}}$ and $\phi ={{\theta}_{1}}-{{\theta}_{2}}$, whose dynamical evolution is governed by

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {} & {{\partial}_{t}}z=2\left(\lambda \sqrt{1-{{z}^{2}}}\sin (\phi )+{{g}_{3}}\left(1-{{z}^{2}}\right)\sin (2\phi )\right), \\ {} & {} & {{\partial}_{t}}\phi =\left(2{{g}_{1}}-{{g}_{2}}\right)z-\frac{2\lambda z}{1-{{z}^{2}}}\cos (\phi )-2{{g}_{3}}z\cos (2\phi ). \end{array}\end{eqnarray} \tag{ 128 }$$

Compared with Josephson effects in double-well Bose–Einstein condensates (Smerzi et al 1997, Zapata et al 1998), the key difference is that here we have $\sin (2\phi )$ and $\cos (2\phi )$) terms which are generated by the Umklapp process g₃. In the ground state, these terms give rise to the spontaneous time-reversal symmetry breaking.

In the noninteracting limit, ${{g}_{1,2,3}}\to 0$, the Rabi-like oscillation with frequency $2\lambda$ is easily recovered. In the linear regime, $|z|\ll 1$, the dynamics in z and ϕ is simplified to

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {} & {{\partial}_{t}}z\approx \left(2\lambda +4{{g}_{3}}\right)\delta \phi, \\ {} & {} & {{\partial}_{t}}\delta \phi \approx \left(2{{g}_{1}}-{{g}_{2}}-2\lambda -2{{g}_{3}}\right)z, \end{array}\end{eqnarray}$$

assuming $\phi \ll 2\pi$. This gives rise to oscillatory dynamics with a frequency

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\omega}_{\text{real}}}=\sqrt{\left(2\lambda +4{{g}_{3}}\right)\left(2\lambda -2{{g}_{1}}+{{g}_{2}}+2{{g}_{3}}\right)}, \end{array}\end{eqnarray} \tag{ 129 }$$

which is the Josephson frequency for a real superposition state ${{p}_{x}}+{{p}_{y}}$. For the complex superposition ${{p}_{x}}\pm \text{i}{{p}_{y}}$, expressing ϕ in terms of fluctuation field $\delta \phi$, $\phi \to \frac{\pi}{2}+\delta \phi$ ($\delta \phi \ll 2\pi$), the linear dynamics is

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {} & {{\partial}_{t}}z\approx -2{{g}_{3}}\left(\delta \phi -\frac{\lambda}{{{g}_{3}}}\right), \\ {} & {} & {{\partial}_{t}}\delta \phi \approx \left(2{{g}_{1}}-{{g}_{2}}+2{{g}_{3}}\right)z, \end{array}\end{eqnarray}$$

which predicts a Josephson frequency

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\omega}_{\text{complex}}}=\sqrt{2{{g}_{3}}\left(2{{g}_{1}}-{{g}_{2}}+2{{g}_{3}}\right)}, \end{array}\end{eqnarray} \tag{ 130 }$$

with $\lambda /{{g}_{3}}$ assumed to be small. In the Josephson effects, the frequency is different from that in non-interacting Rabi oscillations. This frequency difference is also seen in the numerical simulations based on Gross–Pitaevskii equations (Cai et al 2012a) and Gutzwiller methods (Li et al 2016).

The nonlinear effects of dynamics in equation (128) are expected to be more interesting, because of the $\sin (2\phi )$ term, than the usual Josephson physics of double-wells. For example the analogy of self-trapping effect in double-wells would certainly exist in this orbital setting, and very likely would lead to new possibilities beyond the standard double-well Josephson effect. Details of such orbital Josephson effects call for further theoretical and experimental investigations.

Orbital degree of freedom can certainly be mapped to pseudo-spins. In doing so, spin–orbit couplings of certain types usually arise naturally due to the spatial nature of orbitals (Liu et al 2010, 2014, Sun et al 2012a, 2012b, Li et al 2013, Belemuk et al 2014, Zhou et al 2015). The tunneling Hamiltonian of orbital models mixes different orbitals. In particular mixing of different parities could lead to non-trivial effective spin orbit couplings and consequent topological properties. Mixing of s and p-orbitals in a ladder system (Li et al 2013) closely mimics the one dimensional spin orbital coupling recently engineered in cold gases by Raman transitions (Lin et al 2011, Galitski and Spielman 2013). Such sp orbital mixing is recently achieved in a shaken lattice experiment of ¹³³Cs Bose–Einstein condensates (Parker et al 2013) and a similar band structure with double minima like the spin–orbit coupled case is indeed obtained. Mixing of p and d-orbitals gives rise to the phases of topological semimetal and topological insulator (Sun et al 2012b). One recent work shows that mixing of p-orbitals in spin imbalanced fermions leads to topological superconductivity with novel features (Liu et al 2014). We note however that some of these novel predictions are made for fermionic species of atoms, whereas the high band experiments have been explored only for bosons so far as to this time. Further experimental developments are expected.

With strong repulsion, particles could form Mott states with the charge degrees of freedom frozen. The orbital ordering in Mott states is then described by super-exchange interactions of orbitals, which typically depend on the orientation of links. This orientation dependent orbital super-exchange gives rise to novel pseudo-spin models such as quantum 120⁰ model (Zhao and Liu 2008, Wu 2008b) (see equation (83)).

With spin–orbit couplings, many interesting quantum phases such as skyrmions and topological states have been investigated. The connection of orbital physics to spin–orbit coupling suggests possibilities of novel orbital states. One reason to study spin–orbit coupled physics in orbital systems (with atoms loaded into higher bands) is that there appears no additional heating in this system, in contrast with the heating challenge faced by engineered spin–orbital couplings by the advanced Raman laser technique. In this regard, orbital physics provides an alternative platform to investigate spin–orbit coupled phenomena, which is a direction worth future exploration.

In recent optical lattice experiments (Aidelsburger et al 2013, Miyake et al 2013, Parker et al 2013, Struck et al 2013, Jotzu et al 2014, Niu et al 2015, Weinberg et al 2015), periodically driven systems have been developed with a motivation to create exotic atomic phases. In such systems time reversal symmetry is explicitly broken. With the driving frequency matching band gaps, energetically separated orbital bands can be efficiently coupled.

Here we use one example to demonstrate the key idea of using lattice shaking to induce/control orbital couplings. Consider a one dimensional shaking lattice as implemented in experiments (Parker et al 2013). The time-dependent optical potential of this lattice reads

$$\begin{eqnarray}&&V(x,t)={{V}_{0}}\cos \left[k\left(x-{{x}_{0}}(t)\right)\right],\end{eqnarray} \tag{ 131 }$$

with x₀(t) a periodic function, ${{x}_{0}}(t)={{X}_{0}}\sin (2\pi t/T)$. Taking X₀  =  0, we have a static lattice potential where s and p orbital bands are decoupled and well separated by an energy gap. With weak driving, we have $V(x,t)\approx {{V}_{0}}\left[\cos (kx)+k{{x}_{0}}(t)\sin (kx)\right]$. The time-dependent term introduces an effective coupling between s and p orbitals, approximately given by

$$\begin{eqnarray}&&{{\lambda}_{sp}}=k{{V}_{0}}{{x}_{0}}(t){\int}^{}\text{d}x\sin (kx)w_{s}^{\ast}(x){{w}_{p}}(x),\end{eqnarray} \tag{ 132 }$$

with ${{w}_{\nu}}(x)$ the orbital wavefunction. With frequency $2\pi /T$ matching the band gap, the system is approximately described by a static two-band model with s and p orbitals coupled, under a rotating wave approximation.

It appears natural to engineer orbital couplings by lattice modulation/shaking techniques. But the problem is that heating effects are fundamentally unavoidable in periodically driven quantum systems. Since periodic driving breaks time translational symmetry, energy is no longer a conserved quantity. It follows that driven systems (assuming ergodicity) at long time would necessarily be described by infinite temperature ensemble. Nonetheless, there could be long lifetime transient states that manifest interesting topological features. This requires more careful treatment of quantum dynamics than just solving for the ground states of effective static Hamiltonians. One way out is to combine with dissipation. Driven-dissipative orbital models may exhibit steady quantum many-body states with interesting topological properties. This is worth future exploration.

For bosons, firstly, it remains open how to experimentally reach the Mott insulator phases of the p-band and study the p-band superfluid-Mott insulator transition. The current experiments at Hamburg are performed with a two-dimensional checkerboard lattice and a relatively shallow harmonic trap in the third dimension. Introducing an additional optical lattice potential in the third dimension is required to access the Mott regime. Unfortunately that would also increase the on-site interaction between p-orbital bosons, which leads to faster decay (Hemmerich 2014).

Secondly, it is intriguing to find out what type of new topological defects, other than vortices, may possibly occur in the staggered ${{p}_{x}}\pm \text{i}{{p}_{y}}$-orbital Bose–Einstein condensate. The state breaks not only U(1) but also other interesting symmetries that are usually not broken in other conventional Bose condensates, including for example, time-reversal, lattice translational and rotational symmetries. On the general ground of broken symmetries, new classification of topological defects is expected but remains unknown.

For fermions, the stability of the p and higher orbital bands is protected by Fermi statistics, if the experimental system is prepared with the lowest ground band being completely filled, as opposed to the method of band population inversion (Müller et al 2007, Wirth et al 2011, Ölschläger et al 2012, Ölschläger et al 2013, Kock et al 2015). Nevertheless, this approach would require a higher density of fermions, which in turn requires a higher efficiency of cooling fermions down to degeneracy. The recent breakthrough in the Rice experiment of fermions on lattice (Hart et al 2015) is promising for studying the higher orbital bands.