Bosonic Weyl excitations induced by $p$-orbital interactions in a cubic optical lattice

Weyl fermions, predicted by H. Weyl in 1929 Weyl (1929), are massless spin-1/2 particles deriving from the Weyl equation in relativistic quantum field theory. They still haven’t been directly observed as elementary particles. However, the existence of Weyl quasiparticles has been proved in real solid-state materials Wan et al. (2011); Burkov and Balents (2011); Weng et al. (2015); Huang et al. (2015a); Xu et al. (2015a); Lv et al. (2015); Huang et al. (2015b). Distinct from traditional topological insulators Haldane (1988); Hasan and Kane (2010); Qi and Zhang (2011); Bansil et al. (2016), the topological signature of Weyl materials is embodied in three-dimensional (3D) linear dispersion relation near Weyl points, which are paired band crossings and can not be removed unless the Weyl points with opposite charges are merged Turner et al. (2013); Wehling et al. (2014); Armitage et al. (2018).

Ultracold atomic gases provide a highly controllable platform to realize desired quantum states of matter Lewenstein et al. (2007); Bloch et al. (2012); Gross and Bloch (2017); Schäfer et al. (2020). Recently, a great number of significant progresses have been made for simulating Weyl points in ultracold atoms Jiang (2012); Dubček et al. (2015); Zhang et al. (2015); He et al. (2016); Li et al. (2016a); Kong et al. (2017); Zheng et al. (2019); Lu et al. (2020); Wang et al. (2021a); Xu et al. (2014); Liu et al. (2015); Xu et al. (2015b); Lang et al. (2017); Li et al. (2021); Wu et al. (2017), especially an ideal Weyl semimetal band is realized in a quantum gas with 3D spin-orbit coupling Wang et al. (2021a). Most existing researches concentrate on Weyl points induced by spin-orbit couplings and artificial gauge fields in single-particle bands. A few studies focus on Weyl Bogoliubov excitations for fermionic gases Xu et al. (2014); Liu et al. (2015); Xu et al. (2015b); Lang et al. (2017); Li et al. (2021) and bosonic gases Wu et al. (2017). However, the relation between interactions and the emergence of Weyl points has remained elusive.

On the other hand, the investigation of interaction-induced topological phases for bosonic excitations is attracting increasing attention. In two-dimensional (2D) optical lattices, motivated by theoretical proposals Xu et al. (2016); Di Liberto et al. (2016); Luo et al. (2018); Huang et al. (2021), interaction-induced atomic chiral superfluid with topological excitations in a Bose gas was realized in a recent experiment Wang et al. (2021b). Even with these successful instances, the relevant researches are still lacking for the 3D topological phases. Similar to Weyl magnons Li et al. (2016b); Mook et al. (2016); Su et al. (2017); Li et al. (2017); Su and Wang (2017); Owerre (2018); Jian and Nie (2018), interacting bosons in optical lattices also provide a new paradigmatic type of bosonic Weyl quasiparticle excitations.

In this paper, we propose a scheme to realize bosonic Weyl excitations induced by $p$-orbital interactions in a cubic optical lattice. The emergence of Weyl points needs either spatial inversion symmetry breaking or time-reversal symmetry breaking. The special repulsive interaction of $p$ orbitals can induce spontaneous time-reversal symmetry breaking for the bosonic superfluid Liu and Wu (2006). We thus apply this mechanism to realize desired Weyl points for Bogoliubov excitations of the time-reversal symmetry breaking superfluid. The rest of the paper is organized as follows. In Sec. II, we introduce a cubic lattice by a layer-by-layer approach with a special design for $p$ orbitals. In Sec. III, we find the interaction will induce a $p_{x}\pm ip_{y}$ bosonic chiral superfluid with time-reversal symmetry breaking after loading bosons into $p$ orbitals. Therefore Weyl points, which are absent in the single-particle spectra, emerge in the Bogoliubov excitation spectra. In Sec. IV, we discuss the effective two-band models near Weyl points which are obtained by the Krein-unitary perturbation theory deriving from the Schrieffer-Wolff transformation. The edge states of the bosonic excitation modes are also discussed. In Sec. V, we propose a realistic scheme to create the cubic optical lattice in ultracold atom experiment. Distinct from the majority of previous proposals, our scheme does not rely on external Raman lasers or atomic internal states. This paper is concluded in Sec. VI.

We construct a 3D cubic lattice via stacking a 2D square lattice along $z$ direction. One unit cell of the lattice contains two inequivalent lattice sites denoted by A and B, respectively. This is shown in Fig. 1(a). Three $p$-orbitals with $p_{x,y}$ orbitals located at A sites and $p_{z}$ orbitals located at B sites are isolated from remaining orbitals by tuning the lattice parameters. We thus consider loading bosonic atoms into these three $p$-orbitals. The single-particle Hamiltonian can be written as

Here, $\hat{p}_{x,y}$ and $\hat{p}_{z}$ are bosonic annihilation operators for $p_{x,y}$ orbitals at A sites and $p_{z}$ orbitals at B sites, respectively. Tight-binding parameters $t_{x}$, $t_{z}$, $t_{xz}$, $\delta$ are all positive. $\mathbf{e}_{x,y,z}$ are vectors of length $a_{0}$ along $x,y,z$ directions, where $a_{0}$ is the distance between the nearest-neighbor sites. $m=\pm 1$, $n=\pm 1$, and $\alpha=x,y$. Tunnelings between $p_{x,y}$ and $p_{z}$ orbitals in the same $xy$ plane are prohibited due to odd parity of $p$-orbitals. Applying the Fourier transformation, the single-particle Hamiltonian in quasimomentum space can be obtained as $\hat{H}_{0}=\sum_{\mathbf{k}}\hat{\Psi}_{\mathbf{k}}^{\dagger}\mathcal{H}_{0}(\mathbf{k})\hat{\Psi}_{\mathbf{k}}$, where $\hat{\Psi}_{\mathbf{k}}=(\hat{p}_{x,\mathbf{k}},\ \hat{p}_{y,\mathbf{k}},\ \hat{p}_{z,\mathbf{k}})^{\mathrm{T}}$ and $\mathcal{H}_{0}(\mathbf{k})$ is given by

where we have defined $\tilde{k}_{x,y,z}\equiv a_{0}k_{x,y,z}$.

The single-particle energy spectra are obtained by diagonalizing $\mathcal{H}_{0}(\mathbf{k})$. Since the system preserves time-reversal and spatial inversion symmetries, there is no Weyl point for the single-particle spectra. In the following parts, we consider to lift the $p_{z}$ orbitals comparing to $p_{x}$ and $p_{y}$ orbitals, which leads to a positive and large value of $\delta$. In this case, the energy minima are degenerate and located at the plane of $k_{z}=0$. This degeneracy can be broken if we consider more hopping terms, e.g. including next-to-next nearest neighbor tunnelings. For a real optical lattice system as discussed later, the energy minima are located at $M$ point and double degenerate (see Appendix C). We thus consider a Bose-Einstein condensate with atoms condensing at $M$ point.

The interaction among atoms gives rise to the interaction-part Hamiltonian

Here, $\hat{n}_{xy,\mathbf{r}}=\hat{p}_{x,\mathbf{r}}^{\dagger}\hat{p}_{x,\mathbf{r}}+\hat{p}_{y,\mathbf{r}}^{\dagger}\hat{p}_{y,\mathbf{r}}$ and $\hat{L}_{z,\mathbf{r}}=-i(\hat{p}_{x,\mathbf{r}}^{\dagger}\hat{p}_{y,\mathbf{r}}-\hat{p}_{y,\mathbf{r}}^{\dagger}\hat{p}_{x,\mathbf{r}})$. $U_{z}$ and $U_{x}$ are positive for the repulsive interaction. Obviously, the interaction between atoms favors a ground state with $\langle\hat{L}_{z,\mathbf{r}}\rangle\neq 0$, leading to a time-reversal symmetry breaking condensate.

We use the mean-field theory to determine the ground state. Since the atoms are assumed to condense at $M$ point, we obtain the ground-state energy functional as

Here, we define $\bm{\varphi}\equiv(1/\sqrt{N_{0}})(\langle\hat{p}_{x,\mathbf{M}}\rangle,\langle\hat{p}_{y,\mathbf{M}}\rangle,\langle\hat{p}_{z,\mathbf{M}}\rangle)^{\mathrm{T}}=(\varphi_{1},\varphi_{2},\varphi_{3})^{\mathrm{T}}$, and $N_{0}$ is the particle number of condensate. $\rho$ is the atom number per unit cell. Via minimization of the energy functional, we obtain two degenerate ground states with $\bm{\varphi}=\bm{\varphi}_{\pm}$, where $\bm{\varphi}_{\pm}\equiv(1/\sqrt{2})(1,\pm i,0)^{\mathrm{T}}$. The corresponding order parameters are given by

Each ground state breaks the time-reversal symmetry and in the meantime preserves the inversion symmetry. One ground state is shown in Fig. 2(a) for the case with weak interaction among atoms.

We then investigate the Bogoliubov excitations for the time-reversal symmetry breaking ground state with $\bm{\varphi}=\bm{\varphi}_{+}$. The Bogoliubov-de Gennes (BdG) Hamiltonian can be written as $\hat{H}_{\mathrm{BdG}}=(1/2)\sum_{\mathbf{k}\neq\mathbf{M}}\hat{\bm{\phi}}_{\mathbf{k}}^{\dagger}\mathcal{H}_{\mathrm{BdG}}(\mathbf{k})\hat{\bm{\phi}}_{\mathbf{k}}$, where $\hat{\bm{\phi}}_{\mathbf{k}}=(\hat{p}_{x,\mathbf{k}},\hat{p}_{y,\mathbf{k}},\hat{p}_{z,\mathbf{k}},\hat{p}_{x,2\mathbf{M}-\mathbf{k}}^{\dagger}\hat{p}_{y,2\mathbf{M}-\mathbf{k}}^{\dagger},\hat{p}_{z,2\mathbf{M}-\mathbf{k}}^{\dagger})$ and $\mathcal{H}_{\mathrm{BdG}}(\mathbf{k})$ takes the form of

Note that we use “$-\mathbf{k}$” rather than “$2\mathbf{M}-\mathbf{k}$” because the quasimomentum $2\mathbf{M}$ is equivalent to zero. The matrix $h_{a}(\mathbf{k})$ is given by $h_{a}(\mathbf{k})=\mathcal{H}_{0}(\mathbf{k})-\mu 1_{3\times 3}+\rho\mathrm{diag}(4U_{x}/3,4U_{x}/3,0)$. Here $\mu=E_{0}(\mathbf{M})+2U_{x}\rho/3$ with $E_{0}(\mathbf{M})=-2t_{x}$ being the single-particle energy of ground state, $1_{3\times 3}$ is the identity matrix. The matrix $h_{x}(\mathbf{k})$ is given by

To satisfy the bosonic commutation relation, the bosonic BdG Hamiltonian $\mathcal{H}_{\mathrm{BdG}}(\mathbf{k})$ is diagonalized by a paraunitary matrix $T_{\mathbf{k}}$ rather than a unitary matrix as $T_{\mathbf{k}}^{\dagger}\mathcal{H}_{\mathrm{BdG}}(\mathbf{k})T_{\mathbf{k}}=\mathrm{diag}(\mathcal{E}_{\mathbf{k}},\mathcal{E}_{-\mathbf{k}})$. The paraunitary matrix $T_{\mathbf{k}}$ satisfies the relations of $T_{\mathbf{k}}^{\dagger}\tau_{z}T_{\mathbf{k}}=\tau_{z}$ and $T_{\mathbf{k}}\tau_{z}T_{\mathbf{k}}^{\dagger}=\tau_{z}$, where $\tau_{z}=\sigma_{z}\otimes 1_{3\times 3}$ with $\sigma_{z}$ being the Pauli matrix. The excitation spectra $\mathcal{E}_{\mathbf{k}}$ along the line with $(k_{x},k_{y})=(\pi/a_{0},0)$ are depicted with the gray solid lines in the right panel of Fig. 2(b). We note that the degeneracy among $p_{x,y}$ bands (marked by the blue solid line and the yellow dashed line in the left panel of Fig. 2(b)) is lifted by the atomic interaction for the Bogoliubov excitations. The corresponding excitation spectra are given by

where $\mathcal{E}_{i}$ denotes the excitation energy for the $i$-th excitation band. These feature also appears for excitation spectra along the line with $(k_{x},k_{y})=(0,0)$.

Therefore, numerous band crossing points emerge for the Bogoliubov excitations of the time-reversal symmetry breaking superfluid. As we confirmed in the following section, these crossing points are indeed Weyl points. In total, we find that there are eight Weyl points. Half of them are located at the line with $(k_{x},k_{y})=(0,0)$ and the other half are located at the line with $(k_{x},k_{y})=(\pi/a_{0},0)$. On each line, two Weyl points arise at the band crossing point between the first and second excitation bands (characterized by $\mathcal{E}_{1}$ and $\mathcal{E}_{2}$) and another two Weyl points are located at band crossing points between the second and third excitation bands (characterized by $\mathcal{E}_{2}$ and $\mathcal{E}_{3}$). The topological charges of the Weyl points can be obtained by numerical calculations. We show the positions and charges of the Weyl points in Fig. 2(c).

A Hamiltonian to describe the dispersion for a Weyl point can be written as $\mathcal{H}_{W}(\mathbf{k})=\sum_{i,j}k_{i}v_{ij}\sigma_{j}$, where $j=x,y,z$. The topological charge of a Weyl point is determined by $c=\text{sgn}[\det(v_{ij})]$, which equals to $\pm 1$. We apply the perturbation transformation on the bosonic BdG Hamiltonian $\mathcal{H}_{\mathrm{BdG}}(\mathbf{k})$ to obtain an effective low-energy Hamiltonian to describe Weyl points. More details can be found in Appendix A and the related literature Zhou et al. (2020); Wan et al. (2021); Massarelli et al. (2022).

We first define a non-Hermitian matrix $H\equiv\tau_{z}H_{\mathrm{BdG}}(\mathbf{k})$, which is Krein Hermitian with respect to $\tau_{z}$. $H$ can be decomposed into a diagonal matrix $H^{0}$ and a perturbation part $H^{\prime}$, where $H^{\prime}$ can be further decomposed into a block-diagonal matrix $H^{1}$ and a block off-diagonal matrix $H^{2}$. We can block-diagonalize $H$ by the Schrieffer-Wolff transformation

where $e^{W}$ is Krein-unitary. Up to second order, $\tilde{H}$ has the form

where $W$ is determined by $[H^{0},W]+H^{2}=0$.

The Weyl points emerge at the linear band crossings and are mainly induced by the hybridization among $p_{+}$ and $p_{z}$ bands or $p_{-}$ and $p_{z}$ bands in the excitation spectra, where $p_{\pm}=p_{x}\pm ip_{y}$. Therefore, we should transform the bases from $p_{x,y}$ orbitals to $p_{\pm}$ orbitals, and then only focus on the subspace of $p_{+,z}$ or $p_{-,z}$ during the perturbation process. The perturbation method is not available for the entire Brillouin zone due to the divergence, but it does not influence the series expansion near the Weyl points. We apply the series expansion up to first order obtain the final result which satisfies the form of the Weyl Hamiltonian $\mathcal{H}_{W}(\mathbf{k})$. We thus can derive topological charges for Weyl points. In the right panel of Fig. 2(b), the energy spectra for the effective model with only $p_{+,z}$ orbitals and the series expansion near a Weyl point in the effective model are depicted with the blue dashed lines and the red dashed lines, respectively.

The topological charge of a Weyl point is also associated with the Chern number of a 2D surface surrounding the Weyl point. Based on the discussion in Ref. Shindou et al. (2013), the Chern number of $n$-th band can be obtained by the integral of Berry curvature as $c_{n}=1/(2\pi)\int d\mathbf{k}\mathbf{B}_{n}(\mathbf{k})$. The Berry curvature is defined as $\mathbf{B}_{n}(\mathbf{k})=\nabla_{\mathbf{k}}\times\mathbf{A}_{n}(\mathbf{k})$ with $\mathbf{A}_{n}(\mathbf{k})=i\langle T^{n}_{\mathbf{k}}|\tau_{z}\nabla_{\mathbf{k}}|T^{n}_{\mathbf{k}}\rangle$, where $|T^{n}_{\mathbf{k}}\rangle$ is the $n$-th column eigenvector of $T_{\mathbf{k}}$. Numerically, we consider all six surfaces of a tiny cube enclosing a Weyl point and calculate their Berry curvatures by the algorithm described in Ref. Fukui et al. (2005). The Chern number of all surfaces is the topological charge of the Weyl point. The resulting Chern numbers exactly match with the charges obtained from the series expansion. We show the charge of each Weyl point in Fig. 2(c) with $c=+1$ and $c=-1$ being distinguished by red and blue colors, respectively.

Due to bulk-edge correspondence, there are also topological protected edge states analogous to the Fermi arcs in electronic systems for the open system. We thus set an open boundary along $(\hat{x}-\hat{y})/\sqrt{2}$ direction and a periodic boundary along $(\hat{x}+\hat{y})/\sqrt{2}$ and $\hat{z}$ directions, where $k_{\parallel}$ denotes the quasimomentum component along $(\hat{x}+\hat{y})/\sqrt{2}$. One unit cell of the semi-infinite systems contains 30 unit cells of the lattice. In Fig. 3(a) and (b), we show the energy spectra of the semi-infinite system. The density distributions of typical edge states on the lattice surfaces illustrate the bulk-boundary correspondence as shown in Fig. 3(c). More details can be found in Appendix B.

A 2D bipartite square optical lattice has been realized in the experiment Wirth et al. (2011). This optical lattice comprises a deep potential well and a shallow potential well in one unit cell. The lattice potential in $xy$ plane is given by

Here, $k_{L}=2\pi/\lambda$ with $\lambda=1064$ nm being the wavelength of laser beams. This bipartite square lattice is the 2D counterpart of the 3D cubic optical lattice we desired.

We further consider to add extra beams propagating along $x\pm z$ and $y\pm z$ directions to create a 3D cubic lattice with lattice potential $V(\mathbf{r})=V_{\mathrm{2D}}(x,y)+V^{\prime}(\mathbf{r})$, where

As illustrated in Fig. 4, the 3D optical lattice comprises two potential wells in one unit cell and the relative energy offset between the two wells is adjustable. We label the deep well as A and the shallow well as B. The horizontal and vertical slices of the lattice potential are shown in Fig. 4(b) and (c). The unit cell is a cuboid and the nearest distance between A and B is $a_{0}=\pi/k_{L}$. By choosing suitable potential parameters, $p_{x,y}$ orbitals in A site and $p_{z}$ orbital in B site are coupled and far away from other orbitals as shown in Fig. 4(d).

We also perform plane-wave expansion to calculation the energy spectra for the 3D lattice potential and confirm that the tight-binding model of Eq. (1) indeed well describe the band structure. A more practical tight-binding model with more hoppings and the corresponding Bogoliubov excitations have also been discussed. These can be found in Appendix C.

One may load the atoms into the $p$ orbital bands by applying the band swapping method. The interaction strength can be adjusted by the Feshbach resonance technique. The Weyl excitations of the bosonic superfluid may be detected by the momentum-resolved Bragg spectroscopy Ernst et al. (2010). A density wave can form along the edge owing to an interference of the background condensate and the edge mode Furukawa and Ueda (2015). The Hall responses of the bosonic topological Bogoliubov excitations have also been discussed in Ref. Huang et al. (2022).

We investigate bosonic Weyl excitations in a 3D cubic optical lattice with Weyl points emerging in the Bogoliubov excitation spectra. The bosonic Weyl excitations are induced by the interaction of $p$ orbitals. Furthermore, an possible experimental scheme based on the existing experiment of a 2D bipartite optical square lattice is also proposed. This work complements the search of Weyl systems in orbital optical lattices and motivates future experimental studies with realistic ultracold atomic platforms.