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

TL;DR

This paper proposes a novel scheme to realize interaction-induced Weyl points in bosonic superfluids within a 3D optical lattice, revealing topologically protected excitations analogous to electronic Weyl semimetals.

Contribution

It introduces a method to generate Weyl points for bosonic Bogoliubov excitations via $p$-orbital interactions in a cubic optical lattice, a first in bosonic systems.

Findings

Interaction-induced Weyl points in bosonic superfluids

Topologically protected edge modes in the superfluid

Time-reversal symmetry breaking due to $p$-orbital interactions

Abstract

Weyl points exist in a fascinating topological state of matter with linear band crossings analogous to magnetic monopoles. Tremendous efforts have been devoted to investigate fermionic topological matters with Weyl points in the single-particle band dispersion. It remains elusive for realizing interaction-induced Weyl points, especially for bosons. Motivated by recent experimental progress in ultracold atoms, we propose a scheme to create Weyl points for Bogoliubov excitations of a bosonic superfluid in a three-dimensional cubic optical lattice. The unique design of the lattice leads to interaction-induced time-reversal symmetry breaking for a $p$-orbital superfluid, which in turn induces Weyl Bogoliubov excitations. Analogous to Weyl semimetals of electronic systems, the superfluid also support topologically protected edge modes due to the bulk-boundary correspondence.

Equations

$${\hat{H}}_0=$$

$$-\sum _{\mathbf{r}\in A,\pm }{t}_x\left({\hat{p}}_{x,\mathbf{r}}^{†}{\hat{p}}_{x,\mathbf{r}\pm {\mathbf{e}}_z}+{\hat{p}}_{y,\mathbf{r}}^{†}{\hat{p}}_{y,\mathbf{r}\pm {\mathbf{e}}_z}\right)$$

$$+\sum _{\mathbf{r}\in A,m,n,\alpha }\left(mn{t}_{xz}{\hat{p}}_{\alpha ,\mathbf{r}}^{†}{\hat{p}}_{z,\mathbf{r}+m{\mathbf{e}}_{\alpha }+n{\mathbf{e}}_z}+\text{h.c.}\right).$$

$$\left(\begin{array}{ccc}-2t_x\cos ({\tilde{k}}_z)& & -4t_{xz}\sin ({\tilde{k}}_z)\sin ({\tilde{k}}_x)\\ & -2t_x\cos ({\tilde{k}}_z)& -4t_{xz}\sin ({\tilde{k}}_z)\sin ({\tilde{k}}_y)\\ \text{h.c.}& & \delta +2t_z\cos ({\tilde{k}}_z)\end{array}\right),$$

$${\hat{H}}_I=$$

$$+\frac{U_x}{2}\sum _{\mathbf{r}\in A}\left[{\hat{n}}_{xy,\mathbf{r}}\left({\hat{n}}_{xy,\mathbf{r}}-\frac{2}{3}\right)-\frac{{\hat{L}}_{z,\mathbf{r}}^2}{3}\right].$$

$$E_{\mathrm{MF}}(\mathbf{\phi })=$$

$$+\frac{U_x\rho }{6}(4|{\phi }_1{|}^2|{\phi }_2{|}^2+{\phi }_1^{\ast 2}{\phi }_2^2+{\phi }_2^{\ast 2}{\phi }_1^2).$$

$$⟨\hat{\Psi }(\mathbf{r})⟩=e^{i\mathbf{M}\cdot \mathbf{r}}\frac{1}{\sqrt{2}}(1,\pm i,0{)}^{\mathrm{T}}.$$

$$\mathcal{H}_{\mathrm{BdG}}(\mathbf{k})=\left(\begin{array}[]{cc}h_{a}(\mathbf{k})&h_{x}(\mathbf{k})\\ h_{x}^{*}(-\mathbf{k})&h_{a}^{*}(-\mathbf{k})\end{array}\right).$$

$$h_{x}(\mathbf{k})=\frac{U_{x}\rho}{3}\left(\begin{array}[]{ccc}1&i&0\\ i&-1&0\\ 0&0&0\end{array}\right).$$

$${\mathcal{E}}_1(k_z)$$

$${\mathcal{E}}_2(k_z)$$

$${\mathcal{E}}_3(k_z)$$

$$\tilde{H}$$

$$=H+[H,W]+\frac{1}{2}[[H,W],W]+\cdots ,$$

$$\tilde{H}=H^0+H^1+\frac{1}{2}[H^2,W],$$

$$V_{2D}(x,y)=-V_0{\left|\cos (k_{L}x)+\exp (i\theta )\cos (k_{L}y)\right|}^2.$$

$$V^{\prime }(\mathbf{r})=-V^{\prime }{\cos }^2[k_L(x+z)]-V^{\prime }{\cos }^2[k_L(x-z)]$$

$$-V^{\prime }{\cos }^2[k_L(y+z)]-V^{\prime }{\cos }^2[k_L(y-z)],$$

$$T_{\mathbf{k}}^{†}{\tau }_z{T}_{\mathbf{k}}=T_{\mathbf{k}}{\tau }_z{T}_{\mathbf{k}}^{†}={\tau }_z$$

$$T_{\mathbf{k}}^{†}{H}_{\mathrm{BdG}}(\mathbf{k})T_{\mathbf{k}}={\tau }_z(T_{\mathbf{k}}^{-1}({\tau }_z{H}_{\mathrm{BdG}}(\mathbf{k}))T_{\mathbf{k}}).$$

$$H=H^0+H^{\prime }=H^0+H^1+H^2,$$

$$\tilde{H}$$

$$=H+[H,W]+\frac{1}{2}[[H,W],W]+\cdots ,$$

$$[H^0,W]+H^2=0,$$

$$W_{ml}=-\frac{H_{ml}^{\prime }}{E_m^0-E_l^0},$$

$$\tilde{H}=H^0+H^1+\frac{1}{2}[H^2,W]+\cdots .$$

$${\tilde{H}}_{m{m}^{\prime }}=$$

$$+\frac{1}{2}\sum _l{H}_{ml}^{\prime }{H}_{l{m}^{\prime }}^{\prime }\left(\frac{1}{E_m^0-E_l^0}+\frac{1}{E_{m^{\prime }}^0-E_l^0}\right).$$

$${\hat{\mathbf{\varphi }}}_{\mathbf{k}}^{†}{\mathcal{H}}_{\mathrm{BdG}}(\mathbf{k}){\hat{\mathbf{\varphi }}}_{\mathbf{k}}={\hat{\mathbf{\varphi }}}_{\mathbf{k}}^{\prime †}{U}^{†}{\mathcal{H}}_{\mathrm{BdG}}(\mathbf{k})U{\hat{\mathbf{\varphi }}}_{\mathbf{k}}^{\prime },$$

$$U=\left(\begin{array}[]{ccc}1&1\\ i&-i\\ &&\sqrt{2}\end{array}\right)\oplus\left(\begin{array}[]{ccc}1&1\\ -i&i\\ &&\sqrt{2}\end{array}\right)/\sqrt{2}.$$

$${\hat{H}}_{\mathrm{BdG}}$$

$$\displaystyle=\frac{1}{2}\sum_{\mathbf{k}\neq\mathbf{M}}\left(\hat{\bm{\gamma}}_{\mathbf{k}}^{\dagger},\hat{\bm{\gamma}}_{2\mathbf{M}-\mathbf{k}}^{\top}\right)T_{\mathbf{k}}^{\dagger}\mathcal{H}_{\mathrm{BdG}}(\mathbf{k})T_{\mathbf{k}}\left(\begin{array}[]{c}\hat{\bm{\gamma}}_{\mathbf{k}}\\ \hat{\bm{\gamma}}_{2\mathbf{M}-\mathbf{k}}^{\dagger\top}\end{array}\right).$$

$$T_{\mathbf{k}}=\left(\begin{array}[]{cc}\mathbf{u}_{\mathbf{k}}&\mathbf{v}_{2\mathbf{M}-\mathbf{k}}^{*}\\ \mathbf{v}_{\mathbf{k}}&\mathbf{u}_{2\mathbf{M}-\mathbf{k}}^{*}\end{array}\right)$$

$$⟨\Omega |{\hat{\gamma }}_{j,\mathbf{k}}\hat{N}{\hat{\gamma }}_{j,\mathbf{k}}^{†}|\Omega ⟩-⟨\Omega |\hat{N}|\Omega ⟩$$

Full text

I Introduction

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.

II Theoretical model

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

[TABLE]

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

[TABLE]

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

[TABLE]

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.

III Bogoliubov excitations

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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).

IV Topology and edge states

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

[TABLE]

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

[TABLE]

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.

V Experimental implementation

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

[TABLE]

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

[TABLE]

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).

VI Conclusion

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.

Acknowledgements.

We acknowledge the helpful discussion with Liang-Liang Wan. This work is supported by the National Key R&D Program of China (Grants No. 2022YFA1404103 and No. 2018YFA0307200), NSFC (Grant No. 12274196), and a fund from Guangdong province (Grant No. 2019ZT08X324).

Appendix A Effective Hamiltonian near Weyl points

A general anisotropic Weyl Hamiltonian can be written as $\mathcal{H}_{W}(\mathbf{k})=A_{0}+\sum_{i}A_{i}k_{i}\sigma_{0}+\sum_{i,j}k_{i}v_{ij}\sigma_{j}$, where $\sigma_{j}$ are Pauli matrices and $i,j=x,y,z$. To confirm the Weyl points and their charges in the excitation spectra, we try to obtain the low-energy $2\times 2$ Hamiltonians near the Weyl points from the original $6\times 6$ bosonic BdG Hamiltonian $\mathcal{H}_{\mathrm{BdG}}(\mathbf{k})$, the final results are expected to be in the form of $\mathcal{H}_{W}(\mathbf{k})$. Due to the constraint of bosonic commutation relation, the eigenmodes of a $2n\times 2n$ Hermitian matrix $H_{\mathrm{BdG}}(\mathbf{k})$ are solved by a congruence transformation $T_{\mathbf{k}}^{\dagger}H_{\mathrm{BdG}}(\mathbf{k})T_{\mathbf{k}}$. The transformation matrix $T_{\mathbf{k}}$ is paraunitary and satisfies

[TABLE]

with $\tau_{z}=\sigma_{z}\otimes 1_{n\times n}$. The congruence transformation can be rewritten as a similarity transformation of a non-Hermitian matrix $\tau_{z}H_{\mathrm{BdG}}(\mathbf{k})$

[TABLE]

The matrix $\tau_{z}H_{\mathrm{BdG}}(\mathbf{k})$ and the transformation $T_{\mathbf{k}}$ are Krein-Hermitian and Krein-unitary with respect to $\tau_{z}$, respectively. In the following contents of this section, we define $H\equiv\tau_{z}H_{\mathrm{BdG}}(\mathbf{k})$.

We then apply the perturbation theory to approximately describe Weyl points. We decompose $H$ into a zeroth order diagonal matrix $H^{0}$ and a first order perturbation $H^{\prime}$

[TABLE]

where $H^{1}$ is block-diagonal with two subsets and $H^{2}$ is block off-diagonal. $H$ can be block-diagonalized by the Schrieffer-Wolff transformation

[TABLE]

where $e^{W}$ is a paraunitary matrix and $W$ is block off-diagonal like $H^{2}$. Consider the first order of the block off-diagonal part of $\tilde{H}$ is zero, we find

[TABLE]

where $W$ only up to first order can be obtained by

[TABLE]

where $E_{m}^{0}$ are the diagonal elements of $H^{0}$. Up to second order, $\tilde{H}$ is given by

[TABLE]

The matrix elements of $\tilde{H}$ can be expressed as

[TABLE]

This result is obtained by taking up to second order perturbation. In fact, this perturbation method can give rise to the effective model in an arbitrary order Massarelli et al. (2022).

By defining $\hat{p}_{\pm}^{\dagger}=\hat{p}_{x}^{\dagger}\pm i\hat{p}_{y}^{\dagger}$, we can rewrite the bosonic BdG Hamiltonian via

[TABLE]

where $\hat{\bm{\phi}}_{\mathbf{k}}^{\prime}=(\hat{p}_{+,\mathbf{k}},\hat{p}_{-,\mathbf{k}},\hat{p}_{z,\mathbf{k}},\hat{p}_{+,2\mathbf{M}-\mathbf{k}}^{\dagger}\hat{p}_{-,2\mathbf{M}-\mathbf{k}}^{\dagger},\hat{p}_{z,2\mathbf{M}-\mathbf{k}}^{\dagger})$ and

[TABLE]

Since 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 spectrum. 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 second order perturbation method is applied for $H=\tau_{z}(U^{\dagger}\mathcal{H}_{\mathrm{BdG}}(\mathbf{k})U)$. The energy spectra of the effective two-band models are shown in Fig 5. We further obtain all effective $\mathbf{k}\cdot\mathbf{p}$ Hamiltonians near Weyl points up to first order of $\mathbf{k}$ and find that the topological charges from $c=\text{sgn}[\det(v_{ij})]$ are consistent with the numerically calculated Chern numbers. For a typical Weyl Hamiltonian near the crossing point at $(k_{x},k_{y},k_{z})=(\pi/a_{0},0,1.639/a_{0})$ due to the hybridization among $p_{+}$ and $p_{z}$ orbitals, we find $A_{0}=4.0620$, $A_{z}=-0.8439a_{0}$, $v_{xx}=1.4110a_{0}$, $v_{yy}=-1.4110a_{0}$, $v_{zz}=3.1469a_{0}$, and other elements of $A_{i}$ and $v_{ij}$ are zero. For another typical Weyl Hamiltonian near the crossing point of $(k_{x},k_{y},k_{z})=(\pi/a_{0},0,1.1512/a_{0})$ due to the hybridization among $p_{-}$ and $p_{z}$ orbitals, we find $A_{0}=4.5556$, $A_{z}=-0.9985a_{0}$, $v_{xx}=1.4120a_{0}$, $v_{yy}=1.4120a_{0}$, $v_{zz}=2.9954a_{0}$, and other elements of $A_{i}$ and $v_{ij}$ are zero.

Appendix B Edge state densities of the bosonic excitation modes

We consider a semi-infinite lattice structure as shown in the left panel of Fig. 6. The periodic boundaries are applied along $(\hat{x}+\hat{y})/\sqrt{2}$ and $\hat{z}$ directions, their corresponding quasimomenta are denoted by $k_{\parallel}$ and $k_{z}$. Open boundary condition is applied along the $(\hat{x}-\hat{y})/\sqrt{2}$ direction. One unit cell of the semi-infinite lattice contains $N_{L}$ unit cells of the original lattice as shown in the right panel of Fig. 6. Consider that the quasimomentum of the condensate is $\mathbf{M}$, the bosonic BdG Hamiltonian is diagonalized as

[TABLE]

Here, $\hat{\bm{\beta}}_{\mathbf{k}}=(\hat{\beta}_{1,\mathbf{k}},\hat{\beta}_{2,\mathbf{k}},\cdots)^{\mathrm{T}}$ are the bosonic annihilation operators. $T_{\mathbf{k}}^{\dagger}\mathcal{H}_{\mathrm{BdG}}(\mathbf{k})T_{\mathbf{k}}$ is diagonal and $\hat{\bm{\gamma}}_{\mathbf{k}}=(\hat{\gamma}_{1,\mathbf{k}},\hat{\gamma}_{2,\mathbf{k}},\cdots)^{\mathrm{T}}$ are the annihilation operators of the Bogoliubov excitation modes. The paraunitary matrix $T_{\mathbf{k}}$ has the form of

[TABLE]

where the $3N_{L}\times 3N_{L}$ matrices $\mathbf{u}_{\mathbf{k}}$ and $\mathbf{v}_{\mathbf{k}}$ satisfy the normalization condition $\mathbf{u}_{\mathbf{k}}^{\dagger}\mathbf{u}_{\mathbf{k}}-\mathbf{v}_{\mathbf{k}}^{\dagger}\mathbf{v}_{\mathbf{k}}=1_{3N_{L}\times 3N_{L}}$.

The ground state $|\Omega\rangle$ should be the vacuum of all quasiparticles with $\hat{\gamma}_{j,\mathbf{k}}|\Omega\rangle=0$. Then we consider an excitation $\hat{\gamma}_{j,\mathbf{k}}^{\dagger}|\Omega\rangle$ which is corresponding to the $j$-th column vector of the matrix $T_{\mathbf{k}}$. After one quasiparticle being excited, the particle number changes by an amount

[TABLE]

Here, the particle number operator is

[TABLE]

with $N_{0}$ being the particle number of condensate.

We rewrite the $j$-th column vector of the matrix $T_{\mathbf{k}}$ as $|T_{\mathbf{k}}^{j}\rangle=[\cdots,u_{i\alpha}^{j},\cdots,v_{i\alpha}^{j},\cdots]^{\mathrm{T}}$, where $i=1,\cdots,N_{L}$ denote the original unit cells and $\alpha=x,y,z$. In Fig. 3(c) of the main text, we show the population of each original unit cell as $\langle\hat{n}_{i,\mathbf{k}}\rangle=\sum_{\alpha}(|u_{i\alpha}^{j}|^{2}+|v_{i\alpha}^{j}|^{2})$ with $\hat{n}_{i,\mathbf{k}}=\sum_{\alpha}\hat{p}_{i\alpha,\mathbf{k}}^{\dagger}\hat{p}_{i\alpha,\mathbf{k}}$. The populations of the topological edge states mainly concentrate on the neighborhood of the open boundaries, which illustrates the bulk-boundary correspondence.

Appendix C Band structure of the 3D cubic optical lattice

The 3D cubic optical lattice potential is given by $V(\mathbf{r})=V_{\mathrm{2D}}(x,y)+V^{\prime}(\mathbf{r})$, where

[TABLE]

and

[TABLE]

The recoil energy $E_{\mathrm{R}}=\hbar^{2}k_{L}^{2}/(2m_{a})$ is used as the unit of the energy, where $m_{a}$ is the mass of the bosonic atom. We choose lattice parameters as $(V_{0},\theta,V^{\prime})=(11E_{\mathrm{R}},0.484\pi,4.1E_{\mathrm{R}})$ and numerically calculate the exact band structure by plane wave expansion. Fig. 3(a) show the the fourth, fifth, and sixth bands related to $p_{x,y}$ orbitals in deep potential wells and $p_{z}$ orbitals in shallow potential wells. We find that the fourth and fifth bands are only degenerate at $\Gamma$ point ($k_{x}=k_{y}=k_{z}=0$) and $M$ point ($k_{x}=k_{L}$ and $k_{y}=k_{z}=0$) on the plane with $k_{z}=0$. $M$ point is also the position of the energy minima for the fourth single-particle energy band. Thus, we assume that atoms prefer to condensate at $M$ point in the week interaction limit.

To futher verify the theoretical model in the main text, we consider more hoppings that makes the tight-binding model more realistic. The energy spectra of the tight-binding model with more hoppings are shown in Fig. 7(a), which is very close to the exact band structure. The energy minimum is only located at $M$ point. We further use the tight-binding model with more hoppings to study the Bogoliubov excitations. The obtained Bogoliubov excitation spectra and the charges of Weyl points are shown in Fig. 7(b, c), which are consistent with the results obtained from the simplified tight-binding model of Eq. (1) in the main text.