Excitation band topology and edge matter waves in Bose–Einstein condensates in optical lattices

Shunsuke Furukawa; Masahito Ueda

doi:10.1088/1367-2630/17/11/115014

1. Introduction

Topological insulators and superconductors have attracted a great deal of attention in recent years for their rich variety of quantized responses and robust gapless edge states originating from nontrivial topology of bulk Bloch bands [1–3]. A prototype of topological insulators is an integer quantum Hall system [4], which exhibits a quantized Hall conductivity ${\sigma }_{{xy}}=-({e}^{2}/h)C,$ where C is the sum of the Chern numbers of occupied bands [5]. The gapless edge states are characterized by $| C|$ sets of chiral modes propagating clockwise (counterclockwise) along the system's edge for $C\gt 0$ ( $C\lt 0$ ) [6, 7]. Here a nonzero value of C is caused by the breaking of time-reversal symmetry due to a magnetic field. The discovery of ${{\mathbb{Z}}}_{2}$ topological insulators [8–14] has opened up a new avenue for realizing a topologically nontrivial structure in Bloch bands through spin–orbit coupling, without breaking time-reversal symmetry. Topological superconductors have been shown to exhibit exotic edge states consisting of Majorana fermions, which are protected by particle–hole symmetry [15, 16]. A unified understanding of these topological phases has been achieved with a topological periodic table, where such phases are systematically classified for quadratic fermionic Hamiltonians in different dimensions and symmetry classes [17, 18].

Ultracold atomic systems have recently emerged as a new platform for exploring the physics of topological phases, especially owing to ongoing experimental developments for engineering synthetic gauge fields [19, 20] which can be used to produce such states. Different schemes have been proposed and implemented to create nearly uniform magnetic fields in continuum [21], optical lattices [22–25], and synthetic dimensions [26–28]. Furthermore, the Haldane model [29], in which non-uniform fluxes pierce through the system, has been realized using fermionic atoms in a periodically modulated honeycomb optical lattice [30] (see [31–35] for early theoretical proposals for realizing the same or related models). The Haldane model is a prototypical example of Hamiltonians that exhibit topologically distinct regimes characterized by the Chern numbers ${C}_{\pm }$ associated with upper (+) and lower (−) Bloch bands. The phase diagram of this model has been vindicated experimentally using momentum-resolved interband transitions [30]. Another recent remarkable achievement has been an interferometric measurement of the π Berry flux in the momentum space of a honeycomb lattice [36].

A notable feature of atomic systems is that one can study the effect of quantum statistics. For example, by implementing the technique of [30] for bosonic atoms, one can realize a bosonic counterpart of the Haldane model. In the noninteracting case, topological properties of Bloch bands do not depend on quantum statistics. For weakly interacting bosons in optical lattices, Bogoliubov excitation bands give the elementary excitations of Bose–Einstein condensates (BECs). It is then interesting to ask how the topological properties of Bloch bands are carried over to those of Bogoliubov excitation bands in the interacting case.

Band topology of bosonic or classical vibrational modes has been studied in photonic [37–41], phononic [42–44], magnonic [45–47], and polaritonic [48] excitations. Nontrivial Chern numbers of bulk excitation bands give rise to in-gap chiral edge modes, as observed experimentally in photonic systems [39–41]. Ultracold bosonic atoms in optical lattices are expected to offer a unique platform for the studies of band topology because of the high controllability of such systems combined with the macroscopic quantum nature of BECs.

In this paper, we study the topological properties of Bogoliubov excitation bands in BECs in optical lattices with broken time-reversal symmetry, using a Bose–Hubbard extension of the Haldane model (Haldane–Bose–Hubbard model). We show that the topological properties of the Bloch bands in the noninteracting case [29] are smoothly carried over to those of the Bogoliubov excitation bands in the interacting case. Furthermore, the parameter ranges that exhibit nontrivial band topology enlarge with increasing the Hubbard interaction or the particle density (see figure 3). In the presence of sharp boundaries, chiral edge modes appear in the gap between topologically nontrivial excitation bands. We demonstrate that by coherently transferring a portion of the condensate into an edge mode, a density wave is formed along the edge owing to an interference with the background condensate. This property can be used as a macroscopically enhanced signature of an edge mode.

We note that Vasic et al [49] have recently studied the ground-state phase diagram of the Haldane–Bose–Hubbard model, predicting the emergence of uniform and chiral BEC phases and plaquette Mott insulators with loop currents. While the Bogoliubov excitation bands in the uniform and chiral BEC phases have also been studied, their topological properties such as Chern numbers and associated edge modes have not been analyzed in detail. Our work addresses such topological properties of excitations in the uniform BEC phase, clarifies the parameter ranges showing topologically nontrivial bands in the presence of interactions, and demonstrates a unique method of detecting an edge mode though a macroscopic quantum interference. We also note that strong correlation effects on the band topology have been studied in the spin- $\frac{1}{2}$ fermionic Haldane–Hubbard model in [50–52].

The rest of the paper is organized as follows. In section 2, we describe our model, and review the band structure in the noninteracting case. In section 3, we present a Bogoliubov theory for homogeneous condensates with weak repulsive interactions, and analyze the topology of Bogoliubov excitation bands. In sections 4 and 5, we analyze the ground state and excitations of inhomogeneous condensates by using a Bogoliubov–de Gennes (BdG) theory. After describing the basic formalism in section 4, we present numerical results for box and harmonic traps in section 5. In section 6, we present a summary of this paper and discuss an outlook for future studies.

2. Model and band structure

In this section, we first describe our model—a Bose–Hubbard version of the Haldane model in a honeycomb lattice. We then review the band structure in the noninteracting case, and discuss the single-particle ground state into which bosons condense at zero temperature.

2.1. Bose–Hubbard Hamiltonian

We consider bosonic atoms in an optical lattice, which are well described in the tight-binding limit by a Bose–Hubbard model. The Hamiltonian of our system is given by

$\begin{eqnarray}&&H=-\displaystyle \sum _{{\boldsymbol{r}},{{\boldsymbol{r}}}^{\prime }}J({\boldsymbol{r}},{{\boldsymbol{r}}}^{\prime }){a}^{\dagger }({\boldsymbol{r}})a({{\boldsymbol{r}}}^{\prime })+\displaystyle \frac{U}{2}\displaystyle \sum _{{\boldsymbol{r}}}{a}^{\dagger }{({\boldsymbol{r}})}^{2}a{({\boldsymbol{r}})}^{2},\end{eqnarray} \tag{ 1 }$

where ${\boldsymbol{r}}$ and ${{\boldsymbol{r}}}^{\prime }$ run over all the site positions of the lattice, $a({\boldsymbol{r}})$ is the bosonic annihilation operator at the site ${\boldsymbol{r}},$ and U describes the on-site Hubbard interaction. The diagonal element $J({\boldsymbol{r}},{\boldsymbol{r}})$ gives a potential energy at the site ${\boldsymbol{r}},$ and the off-diagonal element $J({\boldsymbol{r}},{{\boldsymbol{r}}}^{\prime })$ with ${\boldsymbol{r}}\ne {{\boldsymbol{r}}}^{\prime }$ describes the (generally complex) hopping amplitude between the two sites and satisfies $J({{\boldsymbol{r}}}^{\prime },{\boldsymbol{r}})={J}^{*}({\boldsymbol{r}},{{\boldsymbol{r}}}^{\prime })$ in order for H to be hermitian. We set the total number of particles to N:

$\begin{eqnarray}&&\displaystyle \sum _{{\boldsymbol{r}}}{a}^{\dagger }({\boldsymbol{r}})a({\boldsymbol{r}})=N.\end{eqnarray} \tag{ 2 }$

We focus on the case in which the hopping terms in equation (1), which will be denoted by H_kin, are given by the Haldane model in a honeycomb lattice [29]; see figure 1. A honeycomb lattice consists of A and B sublattices, and it is spanned by the primitive vectors ${{\boldsymbol{a}}}_{\mathrm{1,2}}$ ( ${{\boldsymbol{a}}}_{3}$ is also introduced for convenience). We also introduce the vectors ${{\boldsymbol{\delta }}}_{\mathrm{1,2,3}},$ which are directed from a B site to the three neighboring A sites. These vectors are given by

$\begin{eqnarray}{{\boldsymbol{\delta }}}_{1} & = & d\left(\displaystyle \frac{1}{2},-\displaystyle \frac{\sqrt{3}}{2}\right),{{\boldsymbol{\delta }}}_{2}=d\left(\displaystyle \frac{1}{2},+\displaystyle \frac{\sqrt{3}}{2}\right),{{\boldsymbol{\delta }}}_{3}=d(-1,0),\\ {{\boldsymbol{a}}}_{1} & = & {{\boldsymbol{\delta }}}_{2}-{{\boldsymbol{\delta }}}_{3},{{\boldsymbol{a}}}_{2}={{\boldsymbol{\delta }}}_{3}-{{\boldsymbol{\delta }}}_{1},{{\boldsymbol{a}}}_{3}={{\boldsymbol{\delta }}}_{1}-{{\boldsymbol{\delta }}}_{2},\end{eqnarray} \tag{ 3 }$

where d is the length between the neighboring sites. We also introduce the reciprocal lattice vectors

$\begin{eqnarray}&&{{\boldsymbol{a}}}_{1}^{*}=\displaystyle \frac{2\pi }{d}\left(\displaystyle \frac{1}{3},\displaystyle \frac{1}{\sqrt{3}}\right),{{\boldsymbol{a}}}_{2}^{*}=\displaystyle \frac{2\pi }{d}\left(-\displaystyle \frac{1}{3},\displaystyle \frac{1}{\sqrt{3}}\right),\end{eqnarray} \tag{ 4 }$

which satisfy ${{\boldsymbol{a}}}_{i}\cdot {{\boldsymbol{a}}}_{j}^{*}=2\pi {\delta }_{{ij}}(i,j=1,2).$

**Figure 1.** (a) Haldane model on a honeycomb lattice. The kinetic part of the model, H_kin, consists of the nearest-neighbor hopping J₁, the next-nearest-neighbor hopping J₂, and the potential difference 2Δ between the two sublattices A and B. Furthermore, atoms acquire an Aharanov–Bohm phase Φ when hopping along every dashed line in the arrowed direction (see equation (5)). (b) The first Brillouin zone of the honeycomb lattice. Here ${{\boldsymbol{a}}}_{1}^{*}$ and ${{\boldsymbol{a}}}_{2}^{*}$ represent the reciprocal lattice vectors (see equation (4)).
Download figure:
Standard image High-resolution image

**Figure 1.** (a) Haldane model on a honeycomb lattice. The kinetic part of the model, H_kin, consists of the nearest-neighbor hopping J₁, the next-nearest-neighbor hopping J₂, and the potential difference 2Δ between the two sublattices A and B. Furthermore, atoms acquire an Aharanov–Bohm phase Φ when hopping along every dashed line in the arrowed direction (see equation (5)). (b) The first Brillouin zone of the honeycomb lattice. Here ${{\boldsymbol{a}}}_{1}^{*}$ and ${{\boldsymbol{a}}}_{2}^{*}$ represent the reciprocal lattice vectors (see equation (4)).
Download figure:
Standard image High-resolution image

The Haldane model consists of the real nearest-neighbor hopping J₁, the complex next-nearest-neighbor hopping ${J}_{2}{{\rm{e}}}^{\pm {\rm{i}}{\rm{\Phi }}},$ and the potential difference 2Δ between the two sublattices. Nonzero values of $J({\boldsymbol{r}},{{\boldsymbol{r}}}^{\prime })$ are thus given by

$\begin{eqnarray}-J({\boldsymbol{r}}+{{\boldsymbol{\delta }}}_{j},{\boldsymbol{r}}) & = & -J({\boldsymbol{r}},{\boldsymbol{r}}+{{\boldsymbol{\delta }}}_{j})=-{J}_{1}{\rm{for}}{\boldsymbol{r}}\in B;\\ -J({\boldsymbol{r}}+{{\boldsymbol{a}}}_{j},{\boldsymbol{r}}) & = & -{J}^{*}({\boldsymbol{r}},{\boldsymbol{r}}+{{\boldsymbol{a}}}_{j})=-{J}_{2}{{\rm{e}}}^{-{\rm{i}}{\epsilon }_{X}{\rm{\Phi }}}{\rm{for}}{\boldsymbol{r}}\in X=A,B;\\ -J({\boldsymbol{r}},{\boldsymbol{r}}) & = & {\epsilon }_{X}{\rm{\Delta }}+V({\boldsymbol{r}}){\rm{for}}{\boldsymbol{r}}\in X=A,B,\end{eqnarray} \tag{ 5 }$

where $j=1,2,3,$ and ${\epsilon }_{A,B}=\pm 1.$ Here, ${\boldsymbol{r}}\in X$ indicates that the site ${\boldsymbol{r}}$ belongs to the X sublattice, and $V({\boldsymbol{r}})$ describes an external potential which depends on the setting of our system. We assume ${J}_{1},{J}_{2}\gt 0$ in the following.

2.2. Band structure in the noninteracting case

Here we set U = 0, and review the band structure of the Haldane model in the noninteracting case [29]. Assuming the periodic boundary conditions in the two directions of the honeycomb lattice, we perform the Fourier expansion

$\begin{eqnarray}&&a({\boldsymbol{r}})=\displaystyle \frac{1}{\sqrt{{N}_{\mathrm{uc}}}}\displaystyle \sum _{{\boldsymbol{k}}}{a}_{X}({\boldsymbol{k}}){{\rm{e}}}^{{\rm{i}}{\boldsymbol{k}}\cdot {\boldsymbol{r}}}({\boldsymbol{r}}\in X=A,B),\end{eqnarray} \tag{ 6 }$

where the sum is taken over the discrete momenta ${\boldsymbol{k}}$ in the first Brillouin zone, and ${N}_{\mathrm{uc}}$ is the total number of unit cells in the system (i.e., half of the total number of sites). The kinetic part of the Hamiltonian (1) is then rewritten as

$\begin{eqnarray}&&{H}_{\mathrm{kin}}=\displaystyle \sum _{{\boldsymbol{k}}}\left({a}_{A}^{\dagger }({\boldsymbol{k}}),{a}_{B}^{\dagger }({\boldsymbol{k}})\right){\mathcal{H}}({\boldsymbol{k}})\left(\begin{array}{c}{a}_{A}({\boldsymbol{k}})\\ {a}_{B}({\boldsymbol{k}})\end{array}\right).\end{eqnarray} \tag{ 7 }$

Here the 2 × 2 hermitian matrix ${\mathcal{H}}({\boldsymbol{k}})$ can be written in the form

$\begin{eqnarray}&&{\mathcal{H}}({\boldsymbol{k}})={h}_{0}({\boldsymbol{k}})I+{\boldsymbol{h}}({\boldsymbol{k}})\cdot {\boldsymbol{\sigma }},\end{eqnarray} \tag{ 8 }$

where I is the identity matrix, ${\boldsymbol{\sigma }}=({\sigma }_{1},{\sigma }_{2},{\sigma }_{3})$ are the Pauli matrices. The coefficients ${h}_{0}({\boldsymbol{k}})$ and ${\boldsymbol{h}}({\boldsymbol{k}})=({h}_{1}({\boldsymbol{k}}),{h}_{2}({\boldsymbol{k}}),{h}_{3}({\boldsymbol{k}}))$ are calculated as

$\begin{eqnarray}{h}_{0}({\boldsymbol{k}}) & = & -2{J}_{2}\mathrm{cos}({\rm{\Phi }})\displaystyle \sum _{j}\mathrm{cos}({\boldsymbol{k}}\cdot {{\boldsymbol{a}}}_{j}),{h}_{1}({\boldsymbol{k}})=-{J}_{1}\displaystyle \sum _{j}\mathrm{cos}({\boldsymbol{k}}\cdot {\delta }_{j}),\\ {h}_{2}({\boldsymbol{k}}) & = & -{J}_{1}\displaystyle \sum _{j}\mathrm{sin}({\boldsymbol{k}}\cdot {\delta }_{j}),{h}_{3}({\boldsymbol{k}})={\rm{\Delta }}+2{J}_{2}\mathrm{sin}({\rm{\Phi }})\displaystyle \sum _{j}\mathrm{sin}({\boldsymbol{k}}\cdot {{\boldsymbol{a}}}_{j}).\end{eqnarray} \tag{ 9 }$

The two energy bands are obtained through the diagonalization of equation (8) as

$\begin{eqnarray}&&{e}_{\pm }({\boldsymbol{k}})={h}_{0}({\boldsymbol{k}})\pm h({\boldsymbol{k}}),h({\boldsymbol{k}}):= | {\boldsymbol{h}}({\boldsymbol{k}})| .\end{eqnarray} \tag{ 10 }$

An example of energy bands is presented in figure 2(a). For noninteracting fermions, complete filling of the lower band leads to a band insulator.

**Figure 2.** (a) An example of a band structure in the noninteracting case, calculated along the path ${M}_{1}\to {K}_{-}\to {K}_{+}\to {M}_{1}$ in the first Brillouin zone shown in figure 1(b). The parameters are chosen to be in the trivial phase with ${C}_{\pm }=0.$ (b) and (c) Examples of Bogoliubov excitation bands in the presence of interaction. A transition in the band topology occurs at ${Un}/{J}_{1}={({Un}/{J}_{1})}_{{\rm{c}}}\simeq 0.321$ (a solution to equation (36)), at which the band gap closes at the ${K}_{-}$ point (see (b)). For ${Un}/{J}_{1}\gt {({Un}/{J}_{1})}_{{\rm{c}}},$ the higher band acquires a nontrivial Chern number ${C}_{+}=+1.$
Download figure:
Standard image High-resolution image

For noninteracting bosons, Bose–Einstein condensation into the lowest-energy single-particle state occurs at zero temperature³ . When ${J}_{2}={\rm{\Delta }}=0,$ the bottom of the lower band ${e}_{-}({\boldsymbol{k}})$ is located at ${\boldsymbol{k}}={\bf{0}}.$ To find whether the position of the bottom can change owing to finite J₂ or Δ, we expand equation (10) around ${\boldsymbol{k}}={\bf{0}}$ as

$\begin{eqnarray}&&{e}_{\pm }({\boldsymbol{k}})={h}_{0}({\bf{0}})+\displaystyle \frac{9}{2}{J}_{2}{k}^{2}{d}^{2}\mathrm{cos}{\rm{\Phi }}\pm \left(h({\bf{0}})-\displaystyle \frac{9}{4}\displaystyle \frac{{J}_{1}^{2}{k}^{2}{d}^{2}}{\sqrt{9{J}_{1}^{2}+{{\rm{\Delta }}}^{2}}}\right)+{\mathcal{O}}({k}^{3}).\end{eqnarray} \tag{ 11 }$

Therefore, ${e}_{-}({\boldsymbol{k}})$ is minimized at ${\boldsymbol{k}}={\bf{0}}$ if the following condition is met:

$\begin{eqnarray}&&\displaystyle \frac{{J}_{1}^{2}}{\sqrt{9{J}_{1}^{2}+{{\rm{\Delta }}}^{2}}}+2{J}_{2}\mathrm{cos}{\rm{\Phi }}\gt 0.\end{eqnarray} \tag{ 12 }$

This condition is satisfied in the regime ${J}_{2},| {\rm{\Delta }}| \ll {J}_{1},$ which is relevant to the Haldane model realized in the scheme of [30] ⁴ .

To determine the single-particle ground state, we parametrize ${\boldsymbol{h}}({\bf{0}})$ using the polar coordinates as

$\begin{eqnarray}&&{\boldsymbol{h}}({\bf{0}})=(-3{J}_{1},0,{\rm{\Delta }})=h({\bf{0}})(\mathrm{sin}{\theta }_{0}\mathrm{cos}\pi ,\mathrm{sin}{\theta }_{0}\mathrm{sin}\pi ,\mathrm{cos}{\theta }_{0}).\end{eqnarray} \tag{ 13 }$

The ${\boldsymbol{k}}={\bf{0}}$ part of H_kin is then diagonalized by means of the transformation

$\begin{eqnarray}&&\left(\begin{array}{c}{a}_{A}({\bf{0}})\\ {a}_{B}({\bf{0}})\end{array}\right)=U({\theta }_{0},\pi )\left(\begin{array}{c}{a}_{+}({\bf{0}})\\ {a}_{-}({\bf{0}})\end{array}\right),\end{eqnarray} \tag{ 14 }$

where we have defined the 2 × 2 unitary matrix

$\begin{eqnarray}U(\theta ,\varphi ):= {{\rm{e}}}^{-{\rm{i}}\varphi (I+{\sigma }_{3})/2}{{\rm{e}}}^{-{\rm{i}}\theta {\sigma }_{2}/2}=\left(\begin{array}{cc}{{\rm{e}}}^{-{\rm{i}}\varphi }\mathrm{cos}(\theta /2) & -{{\rm{e}}}^{-{\rm{i}}\varphi }\mathrm{sin}(\theta /2)\\ \mathrm{sin}(\theta /2) & \mathrm{cos}(\theta /2)\end{array}\right).\end{eqnarray} \tag{ 15 }$

For noninteracting bosons, Bose–Einstein condensation occurs in the mode created by ${a}_{-}^{\dagger }({\bf{0}}).$ For interacting bosons, the condensate wave function is gradually modified with increasing the interaction, as discussed in the next section.

3. Bogoliubov theory and excitation band topology for homogeneous condensates

In this section, we present the Bogoliubov theory [55, 56] for homogeneous condensates with weak repulsive interactions $U\gt 0,$ and determine the band structure of Bogoliubov excitations. Here by 'homogeneous', we refer to the situation in which the system has the periodicity of the honeycomb lattice (we do not require the equivalence of the two sublattices). We then analyze the topology of the Bogoliubov excitation bands, and determine the parameter ranges that exhibit nontrivial topology.

3.1. Bogoliubov theory

To formulate the Bogoliubov theory for the present system, we first need to determine a condensate wave function by using the Gross–Pitaevskii (GP) theory. In the GP theory, we introduce the GP energy functional E by replacing $(a({\boldsymbol{r}}),{a}^{\dagger }({\boldsymbol{r}}))$ by $(\psi ({\boldsymbol{r}}),{\psi }^{*}({\boldsymbol{r}}))$ in the Hamiltonian (1), and minimize it with respect to $(\psi ({\boldsymbol{r}}),{\psi }^{*}({\boldsymbol{r}}))$ under the constraint ${\displaystyle \sum }_{{\boldsymbol{r}}}{| \psi ({\boldsymbol{r}})| }^{2}=N.$ Since the single-particle ground state is formed at ${\boldsymbol{k}}={\bf{0}}$ (as discussed in section 2.2), we introduce the following homogeneous ansatz for the interacting case:

$\begin{eqnarray}&&\psi ({\boldsymbol{r}})=\displaystyle \frac{{\psi }_{X}}{\sqrt{{N}_{\mathrm{uc}}}}({\boldsymbol{r}}\in X).\end{eqnarray} \tag{ 16 }$

We also introduce the chemical potential μ as a Lagrange multiplier to satisfy the particle-number constraint. The functional to be minimized is then given by

$\begin{eqnarray}&&E-\mu N=({\psi }_{A}^{*},{\psi }_{B}^{*})[{\mathcal{H}}({\bf{0}})-\mu I]\left(\begin{array}{c}{\psi }_{A}\\ {\psi }_{B}\end{array}\right)+\displaystyle \frac{U}{2{N}_{\mathrm{uc}}}({| {\psi }_{A}| }^{4}+{| {\psi }_{B}| }^{4}).\end{eqnarray} \tag{ 17 }$

Minimizing this with respect to ${\psi }_{X}^{*}(X=A,B)$ gives a homogeneous version of the GP equations:

$\begin{eqnarray}&&[{\mathcal{H}}({\bf{0}})-\mu I]\left(\begin{array}{c}{\psi }_{A}\\ {\psi }_{B}\end{array}\right)+\displaystyle \frac{U}{{N}_{\mathrm{uc}}}\left(\begin{array}{c}{\psi }_{A}^{*}{\psi }_{A}^{2}\\ {\psi }_{B}^{*}{\psi }_{B}^{2}\end{array}\right)=0.\end{eqnarray} \tag{ 18 }$

Since the single-particle ground state is created by ${a}_{-}^{\dagger }({\bf{0}})$ in equation (14), it is convenient to parametrize ${({\psi }_{A},{\psi }_{B})}^{T}$ as

$\begin{eqnarray}&&\left(\begin{array}{c}{\psi }_{A}\\ {\psi }_{B}\end{array}\right)=\sqrt{N}\left(\begin{array}{c}{f}_{A}\\ {f}_{B}\end{array}\right)=\sqrt{N}\left(\begin{array}{c}\mathrm{sin}(\theta /2)\\ \mathrm{cos}(\theta /2)\end{array}\right),\end{eqnarray} \tag{ 19 }$

where $\theta ={\theta }_{0}$ when U = 0. Multiplying equation (18) by $({f}_{A},{f}_{B})$ or $(-{f}_{B},{f}_{A})$ from the left, we obtain

$\begin{eqnarray}&&{h}_{0}({\bf{0}})-h({\bf{0}})\mathrm{cos}(\theta -{\theta }_{0})+2{Un}({f}_{A}^{4}+{f}_{B}^{4})=\mu ,\end{eqnarray} \tag{ 20a }$

$\begin{eqnarray}&&h({\bf{0}})\mathrm{sin}(\theta -{\theta }_{0})+2{{Unf}}_{A}{f}_{B}({f}_{A}^{2}-{f}_{B}^{2})=0,\end{eqnarray} \tag{ 20b }$

where $n:= N/(2{N}_{\mathrm{uc}})$ is the average number of particles per site. The parameter θ is determined by solving equation (20b); the chemical potential μ is determined by substituting the obtained θ into the left-hand side of (20a). For the obtained θ, we consider the unitary transformation

$\begin{eqnarray}&&\left(\begin{array}{c}{a}_{A}({\bf{0}})\\ {a}_{B}({\bf{0}})\end{array}\right)=U(\theta ,\pi )\left(\begin{array}{c}{a}_{+}\\ {a}_{-}\end{array}\right)=\left(\begin{array}{cc}-{f}_{B} & {f}_{A}\\ {f}_{A} & {f}_{B}\end{array}\right)\left(\begin{array}{c}{a}_{+}\\ {a}_{-}\end{array}\right).\end{eqnarray} \tag{ 21 }$

Bose–Einstein condensation occurs in the lower-band mode created by ${a}_{-}^{\dagger }.$ When ${\rm{\Delta }}=0,$ the two sublattices are equivalent, and thus equation (20b) gives $\theta =\pi /2$ [49]. For $| {\rm{\Delta }}| \ll {J}_{1},$ we can expand equation (20) in terms of $\theta -\pi /2,$ obtaining

$\begin{eqnarray}&&\theta =\displaystyle \frac{\pi }{2}-\displaystyle \frac{{\rm{\Delta }}}{3{J}_{1}+{Un}}+{\mathcal{O}}\left({{\rm{\Delta }}}^{2}/{J}^{2}\right),\end{eqnarray} \tag{ 22a }$

$\begin{eqnarray}&&\mu =-3{J}_{1}+{Un}-6{J}_{2}\mathrm{cos}{\rm{\Phi }}+{\mathcal{O}}\left({{\rm{\Delta }}}^{2}/J\right).\end{eqnarray} \tag{ 22b }$

As seen in equation (22a), the potential difference 2Δ induces a density imbalance between the two sublattices; this imbalance is reduced by a repulsive interaction $U\gt 0.$

We now discuss excitations from the condensate ground state by using the Bogoliubov theory. To this end, using equations (6) and (21), we decompose $a({\boldsymbol{r}})$ into the condensate and noncondensate parts as

$\begin{eqnarray}&&a({\boldsymbol{r}})=\displaystyle \frac{{f}_{X}}{\sqrt{{N}_{\mathrm{uc}}}}{a}_{-}+\tilde{a}({\boldsymbol{r}})({\boldsymbol{r}}\in X).\end{eqnarray} \tag{ 23 }$

Here, the noncondensate part is given by

$\begin{eqnarray}&&\tilde{a}({\boldsymbol{r}})=\displaystyle \frac{1}{\sqrt{{N}_{\mathrm{uc}}}}\left[-{\epsilon }_{X}{f}_{\bar{X}}{a}_{+}+\displaystyle \sum _{{\boldsymbol{k}}\ne {\bf{0}}}{a}_{X}({\boldsymbol{k}}){{\rm{e}}}^{{\rm{i}}{\boldsymbol{k}}\cdot {\boldsymbol{r}}}\right]({\boldsymbol{r}}\in X)\end{eqnarray} \tag{ 24 }$

with $\bar{A}=B$ and $\bar{B}=A.$ Following the Bogoliubov approximation, we replace both ${a}_{-}$ and ${a}_{-}^{\dagger }$ by $\sqrt{N},$ substitute equation (23) into $H-\mu N,$ and expand $H-\mu N$ up to quadratic order in $\tilde{a}({\boldsymbol{r}}).$ The terms linear in $\tilde{a}({\boldsymbol{r}})$ or ${\tilde{a}}^{\dagger }({\boldsymbol{r}})$ disappear because of the stability condition of the condensate in equation (20b), and we arrive at the Bogoliubov Hamiltonian

$\begin{eqnarray}&&H-\mu N=\displaystyle \frac{1}{2}\displaystyle \sum _{{\boldsymbol{k}}\ne {\bf{0}}}{{\boldsymbol{\alpha }}}^{\dagger }({\boldsymbol{k}}){\mathcal{M}}({\boldsymbol{k}}){\boldsymbol{\alpha }}({\boldsymbol{k}})+\displaystyle \frac{1}{2}\left({a}_{+}^{\dagger },{a}_{+}\right){{\mathcal{M}}}_{+}\left(\begin{array}{c}{a}_{+}\\ {a}_{+}^{\dagger }\end{array}\right)+\mathrm{const}.\end{eqnarray} \tag{ 25 }$

with

$\begin{eqnarray}&&{{\boldsymbol{\alpha }}}^{\dagger }({\boldsymbol{k}}):= \left({a}_{A}^{\dagger }({\boldsymbol{k}}),{a}_{B}^{\dagger }({\boldsymbol{k}}),{a}_{A}(-{\boldsymbol{k}}),{a}_{B}(-{\boldsymbol{k}})\right).\end{eqnarray} \tag{ 26 }$

Here, we have introduced the 4 × 4 matrix ${\mathcal{M}}({\boldsymbol{k}})$ and the 2 × 2 matrix ${{\mathcal{M}}}_{+}$ as

$\begin{eqnarray}{\mathcal{M}}({\boldsymbol{k}})=\left(\begin{array}{cc}{\mathcal{H}}({\boldsymbol{k}})-\mu I+4{{UnF}}^{2} & 2{{UnF}}^{2}\\ 2{{UnF}}^{2} & {{\mathcal{H}}}^{T}(-{\boldsymbol{k}})-\mu I+4{{UnF}}^{2}\end{array}\right),\end{eqnarray} \tag{ 27 }$

$\begin{eqnarray}&&{{\mathcal{M}}}_{+}=\left[{h}_{0}({\bf{0}})+h({\bf{0}})\mathrm{cos}(\theta -{\theta }_{0})-\mu +8{{Unf}}_{A}^{2}{f}_{B}^{2}\right]I+4{{Unf}}_{A}^{2}{f}_{B}^{2}{\sigma }_{1}\end{eqnarray} \tag{ 28 }$

with $F:= \mathrm{diag}({f}_{A},{f}_{B}).$

To diagonalize the Bogoliubov Hamiltonian (25), we perform generalized Bogoliubov transformations

$\begin{eqnarray}&&{\boldsymbol{\alpha }}({\boldsymbol{k}})=W({\boldsymbol{k}}){\boldsymbol{\beta }}({\boldsymbol{k}}),\left(\begin{array}{c}{a}_{+}\\ {a}_{+}^{\dagger }\end{array}\right)={W}_{+}\left(\begin{array}{c}{b}_{+}({\bf{0}})\\ {b}_{+}^{\dagger }({\bf{0}})\end{array}\right)\end{eqnarray} \tag{ 29 }$

with

$\begin{eqnarray}&&{{\boldsymbol{\beta }}}^{\dagger }({\boldsymbol{k}}):= \left({b}_{+}^{\dagger }({\boldsymbol{k}}),{b}_{-}^{\dagger }({\boldsymbol{k}}),{b}_{+}(-{\boldsymbol{k}}),{b}_{-}(-{\boldsymbol{k}})\right).\end{eqnarray} \tag{ 30 }$

Here, $W({\boldsymbol{k}})$ and ${W}_{+}$ are paraunitary matrices which satisfy

$\begin{eqnarray}&&{W}^{\dagger }({\boldsymbol{k}}){\tau }_{3}W({\boldsymbol{k}})=W({\boldsymbol{k}}){\tau }_{3}{W}^{\dagger }({\boldsymbol{k}})={\tau }_{3},{W}_{+}^{\dagger }{\sigma }_{3}{W}_{+}={W}_{+}{\sigma }_{3}{W}_{+}^{\dagger }={\sigma }_{3}\end{eqnarray} \tag{ 31 }$

with ${\tau }_{3}:= \mathrm{diag}(1,1,-1,-1).$ These equations ensure the invariance of the bosonic commutation relations. If the matrices $W({\boldsymbol{k}})$ and ${W}_{+}$ are chosen to satisfy

$\begin{eqnarray}{W}^{\dagger }({\boldsymbol{k}}){\mathcal{M}}({\boldsymbol{k}})W({\boldsymbol{k}}) & = & \mathrm{diag}({E}_{+}({\boldsymbol{k}}),{E}_{-}({\boldsymbol{k}}),{E}_{+}(-{\boldsymbol{k}}),{E}_{-}(-{\boldsymbol{k}})),\\ {W}_{+}^{\dagger }{{\mathcal{M}}}_{+}{W}_{+} & = & {E}_{+}({\bf{0}})I,\end{eqnarray} \tag{ 32 }$

the Bogoliubov Hamiltonian (25) is diagonalized as

$\begin{eqnarray}&&H-\mu N=\displaystyle \sum _{{\boldsymbol{k}}}{E}_{+}({\boldsymbol{k}}){b}_{+}^{\dagger }({\boldsymbol{k}}){b}_{+}({\boldsymbol{k}})+\displaystyle \sum _{{\boldsymbol{k}}\ne {\bf{0}}}{E}_{-}({\boldsymbol{k}}){b}_{-}^{\dagger }({\boldsymbol{k}}){b}_{-}({\boldsymbol{k}})+\mathrm{const}..\end{eqnarray} \tag{ 33 }$

The paraunitary matrix $W({\boldsymbol{k}})$ satisfying equation (32) can be constructed numerically by using the method described in [45, 57]. Examples of the calculated Bogoliubov excitation bands ${E}_{\pm }({\boldsymbol{k}})$ are presented in figures 2(b) and (c).

3.2. Band gap

Before discussing the topology of the Bogoliubov excitation bands, let us analyze the gap between the two bands ${E}_{\pm }({\boldsymbol{k}}).$ The band topology cannot change as far as the band gap remains open. On the other hand, the closing of the gap may signal a change in topology.

As seen in figure 2 and known in the noninteracting case [29], the smallest gap is found at one of the two points ${\boldsymbol{k}}=\pm {\boldsymbol{K}}:= \pm \left({{\boldsymbol{a}}}_{1}^{*}+{{\boldsymbol{a}}}_{2}^{*}\right)/3$ in the Brillouin zone; see the ${K}_{\pm }$ points in figure 1. At these points, because of ${h}_{1}(\pm {\boldsymbol{K}})={h}_{2}(\pm {\boldsymbol{K}})=0,$ the matrix ${\mathcal{M}}(\pm {\boldsymbol{K}})$ in equation (27) is decoupled into A and B sublattice blocks, each of which can be diagonalized by the standard 2 × 2 Bogoliubov transformation. The excitation energies at ${\boldsymbol{k}}=\pm {\boldsymbol{K}}$ are then calculated to be

$\begin{eqnarray}&&{E}_{X}(\pm {\boldsymbol{K}})={\lambda }_{X}\pm 3\sqrt{3}{\epsilon }_{X}{J}_{2}\mathrm{sin}{\rm{\Phi }}(X=A,B),\end{eqnarray} \tag{ 34 }$

$\begin{eqnarray}&&{\lambda }_{X}:= {\left[{(3{J}_{2}\mathrm{cos}{\rm{\Phi }}+{\epsilon }_{X}{\rm{\Delta }}-\mu +4{{Unf}}_{X}^{2})}^{2}-{(2{{Unf}}_{X}^{2})}^{2}\right]}^{1/2}.\end{eqnarray} \tag{ 35 }$

The higher (lower) of these energies gives ${E}_{+}(\pm {\boldsymbol{K}})$ ( ${E}_{-}(\pm {\boldsymbol{K}})$ ) in equation (33). The gap-closing conditions at ${\boldsymbol{k}}=\pm {\boldsymbol{K}}$ are thus obtained as

$\begin{eqnarray}&&0={E}_{A}(\pm {\boldsymbol{K}})-{E}_{B}(\pm {\boldsymbol{K}})={\lambda }_{A}-{\lambda }_{B}\pm 6\sqrt{3}{J}_{2}\mathrm{sin}{\rm{\Phi }}.\end{eqnarray} \tag{ 36 }$

For ${J}_{2},| {\rm{\Delta }}| \ll {J}_{1},$ by using (22), we can expand ${\lambda }_{A}-{\lambda }_{B}$ as

$\begin{eqnarray}&&{\lambda }_{A}-{\lambda }_{B}=\displaystyle \frac{2{\rm{\Delta }}}{G({Un}/{J}_{1})}+{\mathcal{O}}\left(\displaystyle \frac{{J}_{2}^{2}}{{J}_{1}},\displaystyle \frac{{{\rm{\Delta }}}^{2}}{{J}_{1}}\right),G(s):= {\left(1+\displaystyle \frac{2s}{3}\right)}^{\displaystyle \frac{1}{2}}\left(1+\displaystyle \frac{s}{3}\right).\end{eqnarray} \tag{ 37 }$

Equation (36) is then rewritten into a simple form:

$\begin{eqnarray}&&{\rm{\Delta }}\pm 3\sqrt{3}{J}_{2}G({Un}/{J}_{1})\mathrm{sin}{\rm{\Phi }}+{\mathcal{O}}\left({J}_{2}^{2}/{J}_{1},{{\rm{\Delta }}}^{2}/{J}_{1}\right)=0.\end{eqnarray} \tag{ 38 }$

For ${Un}=0,$ this reduces to the exact phase boundaries ${\rm{\Delta }}=\mp 3\sqrt{3}{J}_{2}\mathrm{sin}{\rm{\Phi }}$ in the noninteracting case [29]. Equation (38) indicates that with increasing the strength of interaction ${Un}/{J}_{1},$ these boundaries are shifted to larger $| {\rm{\Delta }}|$ by a factor of $G({Un}/{J}_{1});$ see figure 3. Here, $G(s)$ is an increasing function of s, and expanded for $| s| \ll 1$ as

$\begin{eqnarray}&&G(s)=1+2\;s/3+{s}^{2}/18+{\mathcal{O}}({s}^{4}).\end{eqnarray} \tag{ 39 }$

Across the gap-closing lines, the topology of the Bogoliubov excitation bands changes, as we discuss next.

**Figure 3.** Chern number ${C}_{+}$ of the higher excitation band in the Haldane–Bose–Hubbard model with ${J}_{2},| {\rm{\Delta }}| \ll {J}_{1}.$ The boundaries between regions of different ${C}_{+}$ values are given by equation (38), and indicated for ${Un}/{J}_{1}=0,0.5,1.$ The solid curves correspond to the noninteracting case [29]. The regions with nontrivial topology ( ${C}_{+}\ne 0$ ) enlarge as ${Un}/{J}_{1}$ increases.
Download figure:
Standard image High-resolution image

**Figure 3.** Chern number ${C}_{+}$ of the higher excitation band in the Haldane–Bose–Hubbard model with ${J}_{2},| {\rm{\Delta }}| \ll {J}_{1}.$ The boundaries between regions of different ${C}_{+}$ values are given by equation (38), and indicated for ${Un}/{J}_{1}=0,0.5,1.$ The solid curves correspond to the noninteracting case [29]. The regions with nontrivial topology ( ${C}_{+}\ne 0$ ) enlarge as ${Un}/{J}_{1}$ increases.
Download figure:
Standard image High-resolution image

3.3. Chern number

We now analyze the topology of the Bogoliubov excitation bands. We first note that technically, ${\mathcal{M}}({\boldsymbol{k}})$ in equation (27) (more specifically, ${\mathcal{H}}({\boldsymbol{k}})$ in it) needs a modification for such a purpose because it is not periodic in the first Brillouin zone in the present representation. This is because the Fourier expansion (6) is based on the real-space positions ${\boldsymbol{r}}$ —while this treatment is useful in making ${\mathcal{H}}({\boldsymbol{k}})$ possess the C₃ symmetry of the original lattice, the spacing between the two sublattices introduces an additional phase factor, which breaks the periodicity in the first Brillouin zone. To recover the periodicity, we define $\tilde{{\mathcal{M}}}({\boldsymbol{k}})$ by replacing ${\mathcal{H}}({\boldsymbol{k}})$ by $\tilde{{\mathcal{H}}}({\boldsymbol{k}})={{\rm{e}}}^{-{\rm{i}}{\boldsymbol{k}}\cdot {{\boldsymbol{\delta }}}_{3}(I-{\sigma }_{3})/2}{\mathcal{H}}({\boldsymbol{k}}){{\rm{e}}}^{{\rm{i}}{\boldsymbol{k}}\cdot {{\boldsymbol{\delta }}}_{3}(I-{\sigma }_{3})/2}$ in equation (27). We then introduce the paraunitary matrix $\tilde{W}({\boldsymbol{k}})$ as the matrix that 'diagonalizes' $\tilde{{\mathcal{M}}}({\boldsymbol{k}})$ in the sense of equation (32).

The 4 × 4 paraunitary matrix $\tilde{W}({\boldsymbol{k}})$ can be parametrized as

$\begin{eqnarray}\tilde{W}({\boldsymbol{k}})=\left(\begin{array}{cc}{\mathcal{U}}({\boldsymbol{k}}) & {{\mathcal{V}}}^{*}(-{\boldsymbol{k}})\\ {\mathcal{V}}({\boldsymbol{k}}) & {{\mathcal{U}}}^{*}(-{\boldsymbol{k}})\end{array}\right)\end{eqnarray} \tag{ 40 }$

with

$\begin{eqnarray}{\mathcal{U}}({\boldsymbol{k}})=\left(\begin{array}{cc}{u}_{A+}({\boldsymbol{k}}) & {u}_{A-}({\boldsymbol{k}})\\ {u}_{B+}({\boldsymbol{k}}) & {u}_{B-}({\boldsymbol{k}})\end{array}\right),{\mathcal{V}}({\boldsymbol{k}})=\left(\begin{array}{cc}{v}_{A+}({\boldsymbol{k}}) & {v}_{A-}({\boldsymbol{k}})\\ {v}_{B+}({\boldsymbol{k}}) & {v}_{B-}({\boldsymbol{k}})\end{array}\right).\end{eqnarray} \tag{ 41 }$

To discuss the topology of each excitation band, we introduce the vectors

$\begin{eqnarray}&&| {w}_{\gamma }({\boldsymbol{k}})\rangle ={({u}_{A\gamma }({\boldsymbol{k}}),{u}_{B\gamma }({\boldsymbol{k}}),{v}_{A\gamma }({\boldsymbol{k}}),{v}_{B\gamma }({\boldsymbol{k}}))}^{T}(\gamma =\pm ),\end{eqnarray} \tag{ 42 }$

where $\gamma =\;+\;$ and − correspond respectively to the first and second columns of $W({\boldsymbol{k}}).$ It follows from equations (31) and (32) that these vectors satisfy the eigen equation

$\begin{eqnarray}&&\tilde{{\mathcal{M}}}({\boldsymbol{k}})| {w}_{\gamma }({\boldsymbol{k}})\rangle ={E}_{\gamma }({\boldsymbol{k}}){\tau }_{3}| {w}_{\gamma }({\boldsymbol{k}})\rangle \end{eqnarray} \tag{ 43 }$

and the orthonormality condition

$\begin{eqnarray}&&\langle {w}_{\gamma }({\boldsymbol{k}})| {\tau }_{3}| {w}_{{\gamma }^{\prime }}({\boldsymbol{k}})\rangle ={\delta }_{\gamma {\gamma }^{\prime }}.\end{eqnarray} \tag{ 44 }$

For each band, we introduce the Berry curvature [37, 45]

$\begin{eqnarray}&&{{\mathcal{B}}}_{\gamma }({\boldsymbol{k}})={\rm{i}}{\epsilon }_{{ij}}\langle {\partial }_{i}{w}_{\gamma }({\boldsymbol{k}})| {\tau }_{3}| {\partial }_{j}{w}_{\gamma }({\boldsymbol{k}})\rangle \end{eqnarray} \tag{ 45 }$

with ${\partial }_{i}:= \frac{\partial }{\partial {k}_{i}}.$ The Chern number can then be introduced as [5, 37, 45]

$\begin{eqnarray}&&{C}_{\gamma }={\displaystyle \int }_{\mathrm{BZ}}\displaystyle \frac{{{\rm{d}}}^{2}{\boldsymbol{k}}}{2\pi }{{\mathcal{B}}}_{\gamma }({\boldsymbol{k}}).\end{eqnarray} \tag{ 46 }$

In the noninteracting case ${Un}=0,$ both ${C}_{\pm }$ are quantized to integers, and satisfy the zero sum rule ${\displaystyle \sum }_{\gamma }{C}_{\gamma }=0.$ In the interacting case ${Un}\gt 0,$ however, $| {w}_{-}({\boldsymbol{k}})\rangle$ is not defined at⁵ ${\boldsymbol{k}}={\bf{0}}$ , and thus there is an ambiguity in the definition of ${C}_{-}.$ Nevertheless, ${C}_{+}$ is still well-defined, and quantized to an integer⁶ . We numerically calculate ${C}_{+}$ using the manifestly gauge-invariant method proposed in [58].

The 'phase diagram' based on the Chern number ${C}_{+}$ is presented in figure 3 (we note that this diagram is not based on ground-state transitions). We have numerically confirmed that the boundaries between regions of different ${C}_{+}$ values are given precisely by the gap-closing condition (36) (or equation (38) for ${J}_{2},| {\rm{\Delta }}| \ll {J}_{1}$ ) obtained in section 3.2. The obtained results indicate that the topology of the Bloch bands known for the noninteracting case ${Un}=0$ [29] are smoothly carried over to that of the Bogoliubov excitation bands for the interacting case ${Un}\gt 0,$ and that the regions displaying nontrivial topology ${C}_{+}\ne 0$ enlarge with increasing ${Un}/{J}_{1}.$ When ${C}_{+}\ne 0,$ the bulk-edge correspondence [7] dictates that in-gap chiral edge modes intervening between the upper and lower bulk bands appear when the system has a boundary. We numerically demonstrate the emergence of such modes in section 5.

In closing this section, two remarks are in order.

The first remark is on the reason why the regions displaying topological bands expand with increasing the interaction U. In the noninteracting case U = 0, the potential difference 2Δ between the two sublattices drives a transition from topological to trivial bands as it favors sublattice-separated Bloch wave functions, which have trivial topology. Algebraically, this potential difference induces a finite difference between ${\lambda }_{A}$ and ${\lambda }_{B}$ defined in equation (35), and the transition occurs when $| {\lambda }_{A}-{\lambda }_{B}| =6\sqrt{3}{J}_{2}| \mathrm{sin}{\rm{\Phi }}| .$ The interaction $U\gt 0$ has the effect of obscuring this difference (through ${{Unf}}_{X}^{2}$ in equation (35)), and thus a larger $| {\rm{\Delta }}|$ is required to drive the topological-to-trivial transition. It will be interesting to investigate whether a similar stabilization of topological bands occurs in a wider variety of interacting systems.

The second remark is on the case of attractive interactions $U\lt 0.$ While an attractive Bose gas is unstable against collapse in the thermodynamic limit, it can form a metastable condensate in a finite system if N is below a certain critical value [59]. As far as a quasi-homogeneous condensate is realized, we can perform the same Bogoliubov analysis as in this section, and obtain the phase boundaries of topologically nontrivial regions as in equation (38); these regions gradually shrink with increasing $| U| n$ for $U\lt 0.$ In these regions, in-gap chiral edge modes discussed in section 5 are also expected to be formed.

4. Ground state and excitations in trapped condensates: formalism

In this section, we present the formalism for calculating the ground state and excitations of trapped condensates. We first describe the BdG theory [55, 56] for inhomogneous condensates on lattices. We then apply this theory to a strip geometry, which is convenient for discussing edge modes. We also describe an extended Thomas–Fermi approximation which can give a simple analytic expression for the density profile in a given trap potential. Numerical results obtained using the formalism are presented in section 5.

4.1. BdG theory for inhomogeneous condensates

The BdG theory for inhomogeneous condensates can be derived by linearizing a time-dependent GP equation. We start from the Heisenberg equation of motion for the time-dependent operator $a({\boldsymbol{r}},t):$

$\begin{eqnarray}&&{\rm{i}}{\rm{\hslash }}{\partial }_{t}a({\boldsymbol{r}},t)=[a({\boldsymbol{r}},t),H]=-\displaystyle \sum _{{{\boldsymbol{r}}}^{\prime }}J({\boldsymbol{r}},{{\boldsymbol{r}}}^{\prime })a({{\boldsymbol{r}}}^{\prime },t)+{{Ua}}^{\dagger }({\boldsymbol{r}},t)a{({\boldsymbol{r}},t)}^{2}.\end{eqnarray} \tag{ 47 }$

Replacing $(a({\boldsymbol{r}},t),{a}^{\dagger }({\boldsymbol{r}},t))$ by classical fields $(\psi ({\boldsymbol{r}},t),{\psi }^{*}({\boldsymbol{r}},t)),$ we obtain the GP equation

$\begin{eqnarray}&&i{\rm{\hslash }}{\partial }_{t}\psi ({\boldsymbol{r}},t)=-\displaystyle \sum _{{{\boldsymbol{r}}}^{\prime }}J({\boldsymbol{r}},{{\boldsymbol{r}}}^{\prime })\psi ({{\boldsymbol{r}}}^{\prime },t)+U{\psi }^{*}({\boldsymbol{r}},t)\psi {({\boldsymbol{r}},t)}^{2}.\end{eqnarray} \tag{ 48 }$

Equation (2) imposes the normalization condition ${\displaystyle \sum }_{{\boldsymbol{r}}}{| \psi ({\boldsymbol{r}},t)| }^{2}=N.$ Inserting a stationary-state ansatz $\psi ({\boldsymbol{r}},t)=\psi ({\boldsymbol{r}}){{\rm{e}}}^{-{\rm{i}}\mu t/{\rm{\hslash }}}$ into equation (48), we obtain the time-independent GP equation

$\begin{eqnarray}&&\mu \psi ({\boldsymbol{r}})=-\displaystyle \sum _{{{\boldsymbol{r}}}^{\prime }}J({\boldsymbol{r}},{{\boldsymbol{r}}}^{\prime })\psi ({{\boldsymbol{r}}}^{\prime })+U{\psi }^{*}({\boldsymbol{r}})\psi {({\boldsymbol{r}})}^{2}.\end{eqnarray} \tag{ 49 }$

The solution $\psi ({\boldsymbol{r}})$ with the lowest frequency $\mu /{\rm{\hslash }}$ gives the GP ground state.

We now discuss small fluctuations around the GP ground state: $\psi ({\boldsymbol{r}},t)=\psi ({\boldsymbol{r}}){{\rm{e}}}^{-{\rm{i}}\mu t/{\rm{\hslash }}}+\phi ({\boldsymbol{r}},t).$ Expanding the GP equation (48) to first order in $\phi ({\boldsymbol{r}}),$ we obtain a linear differential equation

$\begin{eqnarray}&&{\rm{i}}{\rm{\hslash }}{\partial }_{t}\phi ({\boldsymbol{r}},t)=-\displaystyle \sum _{{{\boldsymbol{r}}}^{\prime }}J({\boldsymbol{r}},{{\boldsymbol{r}}}^{\prime })\phi ({{\boldsymbol{r}}}^{\prime },t)+2U{| \psi ({\boldsymbol{r}})| }^{2}\phi ({\boldsymbol{r}},t)+U{{\rm{e}}}^{-2{\rm{i}}\mu t/{\rm{\hslash }}}\psi {({\boldsymbol{r}})}^{2}{\phi }^{*}({\boldsymbol{r}},t).\end{eqnarray} \tag{ 50 }$

Assuming a solution of the form

$\begin{eqnarray}&&\phi ({\boldsymbol{r}},t)={{\rm{e}}}^{-{\rm{i}}\mu t/{\rm{\hslash }}}\left[u({\boldsymbol{r}}){{\rm{e}}}^{-{\rm{i}}{Et}/{\rm{\hslash }}}+{v}^{*}({\boldsymbol{r}}){{\rm{e}}}^{{\rm{i}}{Et}/{\rm{\hslash }}}\right],\end{eqnarray} \tag{ 51 }$

we obtain BdG equations

$\begin{eqnarray}&&+{Eu}({\boldsymbol{r}})=-\displaystyle \sum _{{{\boldsymbol{r}}}^{\prime }}J({\boldsymbol{r}},{{\boldsymbol{r}}}^{\prime })u({{\boldsymbol{r}}}^{\prime })+(2U{| \psi ({\boldsymbol{r}})| }^{2}-\mu )u({\boldsymbol{r}})+U\psi {({\boldsymbol{r}})}^{2}v({\boldsymbol{r}}),\end{eqnarray} \tag{ 52a }$

$\begin{eqnarray}&&-{Ev}({\boldsymbol{r}})=-\displaystyle \sum _{{{\boldsymbol{r}}}^{\prime }}{J}^{*}({\boldsymbol{r}},{{\boldsymbol{r}}}^{\prime })v({{\boldsymbol{r}}}^{\prime })+(2U{| \psi ({\boldsymbol{r}})| }^{2}-\mu )v({\boldsymbol{r}})+U{\psi }^{*}{({\boldsymbol{r}})}^{2}u({\boldsymbol{r}}).\end{eqnarray} \tag{ 52b }$

If the condensate is stable, these equations admit ${N}_{s}-1$ (linearly independent) sets of solutions $({u}_{j}({\boldsymbol{r}}),{v}_{j}({\boldsymbol{r}}))$ with positive frequencies ${E}_{j}/{\rm{\hslash }},$ where N_s is the total number of lattice sites. Such solutions can be chosen to satisfy the orthonormality condition

$\begin{eqnarray}&&\displaystyle \sum _{{\boldsymbol{r}}}[{u}_{j}^{*}({\boldsymbol{r}}){u}_{l}({\boldsymbol{r}})-{v}_{j}^{*}({\boldsymbol{r}}){v}_{l}({\boldsymbol{r}})]={\delta }_{{jl}}.\end{eqnarray} \tag{ 53 }$

The present formulation of the BdG theory is based on the linearization of the GP equation (48), which is an equation of motion for the classical field. However, we can check the consistency of this classical-field formulation with the operator formulation for the homogeneous case in section 3. Indeed, substituting the ansatz (16) into the GP equation (49) reproduces equation (18). Furthermore, substituting

$\begin{eqnarray}&&u({\boldsymbol{r}})=\displaystyle \frac{{u}_{X\gamma }({\boldsymbol{k}})}{\sqrt{{N}_{\mathrm{uc}}}}{{\rm{e}}}^{{\rm{i}}{\boldsymbol{k}}\cdot {\boldsymbol{r}}},v({\boldsymbol{r}})=\displaystyle \frac{{v}_{X\gamma }({\boldsymbol{k}})}{\sqrt{{N}_{\mathrm{uc}}}}{{\rm{e}}}^{{\rm{i}}{\boldsymbol{k}}\cdot ({\boldsymbol{r}}-\displaystyle \frac{1-{\epsilon }_{X}}{2}{{\boldsymbol{\delta }}}_{3})}({\boldsymbol{r}}\in X)\end{eqnarray} \tag{ 54 }$

into the BdG equation (52) reproduces equation (43). It is known that the operator formulation for the inhomogeneous case also leads to the same set of BdG equations as in equation (52) [55, 56].

4.2. Strip geometry

We apply the BdG theory to analyze excitations in trapped condensates. We describe site positions by two-dimensional coordinates ${\boldsymbol{r}}=(x,y).$ For simplicity, we consider a strip geometry, which is periodic only along the y direction as in figure 4 and analogous to a graphene nanoribbon. This geometry can be used to describe the central part of an elongated condensate around which the system is approximately uniform in the elongated direction. Along the x direction, we introduce a box trap with sharp zigzag edges (figure 4) or a harmonic trap [60] in section 5; however, we do not assume a specific trap potential in this subsection.

**Figure 4.** Strip (or 'nanoribon') geometry with zigzag edges. A periodic boundary condition is imposed only along the y direction. The red rectangle indicates a unit cell which contains ${N}_{x}=16$ lattice sites. We denote the number of unit cells by N_y. The origin of the x coordinate is placed at the center of the system.
Download figure:
Standard image High-resolution image

**Figure 4.** Strip (or 'nanoribon') geometry with zigzag edges. A periodic boundary condition is imposed only along the y direction. The red rectangle indicates a unit cell which contains ${N}_{x}=16$ lattice sites. We denote the number of unit cells by N_y. The origin of the x coordinate is placed at the center of the system.
Download figure:
Standard image High-resolution image

Exploiting the translation invariance along the y direction, we make the following ansatz for the GP ground state:

$\begin{eqnarray}&&\psi (x,y)=\displaystyle \frac{{\psi }_{x}}{\sqrt{{N}_{y}}}=\sqrt{{n}_{y}}{f}_{x},{n}_{y}:= \displaystyle \frac{N}{{N}_{y}},\end{eqnarray} \tag{ 55 }$

where N_y is the number of unit cells, and ${\psi }_{x}=\sqrt{N}{f}_{x}$ satisfies the normalization condition ${\displaystyle \sum }_{x}{| {\psi }_{x}| }^{2}=N{\displaystyle \sum }_{x}{| {f}_{x}| }^{2}=N.$ Substituting this into the GP equation (49), we obtain

$\begin{eqnarray}&&\mu {f}_{x}=\displaystyle \sum _{{x}^{\prime }}{{\mathcal{H}}}_{{{xx}}^{\prime }}(0){f}_{{x}^{\prime }}+{{Un}}_{y}{f}_{x}^{*}{f}_{x}^{2}.\end{eqnarray} \tag{ 56 }$

Here we have introduced

$\begin{eqnarray}&&{{\mathcal{H}}}_{{{xx}}^{\prime }}({k}_{y})=-\displaystyle \sum _{y}J(x,y,{x}^{\prime },{y}^{\prime }){{\rm{e}}}^{-{\rm{i}}{k}_{y}(y-{y}^{\prime })},\end{eqnarray} \tag{ 57 }$

where $({x}^{\prime },{y}^{\prime })$ is a particular site position, and the sum over y is restricted to the site positions for fixed x; the translation invariance along the y direction ensures that the right-hand side (rhs) of equation (57) does not depend on y'. An accurate solution to equation (56) with the lowest frequency $\mu /{\rm{\hslash }}$ can be obtained by numerically performing the imaginary time evolution with the replacement $\mu \to -{\partial }_{\tau }.$

After obtaining the GP ground state (55), we solve the BdG equations (52) to calculate excitations. We introduce the following ansatz with momentum k_y in the y direction:

$\begin{eqnarray}&&u(x,y)=\displaystyle \frac{{u}_{x}}{\sqrt{{N}_{y}}}{{\rm{e}}}^{{\rm{i}}{k}_{y}y},v(x,y)=\displaystyle \frac{{v}_{x}}{\sqrt{{N}_{y}}}{{\rm{e}}}^{{\rm{i}}{k}_{y}y}.\end{eqnarray} \tag{ 58 }$

Substituting these into the BdG equation (52), we obtain an eigen equation

$\begin{eqnarray}&&{\mathcal{M}}({k}_{y}){\boldsymbol{w}}=E{\tau }_{3}{\boldsymbol{w}}.\end{eqnarray} \tag{ 59 }$

Here we have introduced a $2{N}_{x}$ -component vector ${\boldsymbol{w}}={(\{{u}_{x}\},\{{v}_{x}\})}^{T},$ and a $2{N}_{x}\times 2{N}_{x}$ matrix

$\begin{eqnarray}{\mathcal{M}}({k}_{y})=\left(\begin{array}{cc}{\mathcal{H}}({k}_{y})-\mu I+2{{Un}}_{y}{F}^{*}F & {{Un}}_{y}{F}^{2}\\ {{Un}}_{y}{F}^{*2} & {{\mathcal{H}}}^{T}(-{k}_{y})-\mu I+2{{Un}}_{y}{F}^{*}F\end{array}\right),\end{eqnarray} \tag{ 60 }$

with $F:= \mathrm{diag}(\{{f}_{x}\}).$ The eigen equation (59) can be solved numerically by using the method of [57]. We denote positive-frequency solutions by ${{\boldsymbol{w}}}_{j}({k}_{y})={(\{{u}_{{xj}}({k}_{y})\},\{{v}_{{xj}}({k}_{y})\})}^{T}$ and ${E}_{j}({k}_{y})$ with $j=1,2,...,{N}_{x}$ in descending order in energy. In the operator formulation of the BdG theory, this leads to the fact that the Hamiltonian is diagonalized as

$\begin{eqnarray}&&H-\mu N=\displaystyle \sum _{({k}_{y},j)\ne (0,{N}_{x})}{E}_{j}({k}_{y}){b}_{j}^{\dagger }({k}_{y}){b}_{j}({k}_{y})+\mathrm{const}.,\end{eqnarray} \tag{ 61 }$

where the new bosonic operators $\{{b}_{j}({k}_{y})\}$ are related to the original ones as

$\begin{eqnarray}&&a(x,y)=\psi (x,y)+\displaystyle \frac{1}{\sqrt{{N}_{y}}}\displaystyle \sum _{({k}_{y},j)\ne (0,{N}_{x})}{{\rm{e}}}^{{\rm{i}}{k}_{y}y}\left[{u}_{{xj}}({k}_{y}){b}_{j}({k}_{y})+{v}_{{xj}}^{*}(-{k}_{y}){b}_{j}^{\dagger }(-{k}_{y})\right].\end{eqnarray} \tag{ 62 }$

To discuss excitations in a strip geometry, it is useful to introduce the spectral weight at zero temperature:

$\begin{eqnarray}&&\rho (x,{x}^{\prime },{k}_{y},\omega )={\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}t\displaystyle \sum _{y}{{\rm{e}}}^{-{\rm{i}}{k}_{y}(y-{y}^{\prime })+{\rm{i}}\omega t}\langle [a(x,y,t),{a}^{\dagger }({x}^{\prime },{y}^{\prime },0)]\rangle ,\end{eqnarray} \tag{ 63 }$

where $(x,y)$ and $({x}^{\prime },{y}^{\prime })$ are taken similarly as in equation (57), and the average $\langle \cdot \rangle$ is taken over the vacuum of the Bogoliubov excitations in equation (61). Using equations (61) and (62), equation (63) can be rewritten as

$\begin{eqnarray}\rho (x,{x}^{\prime },{k}_{y},\omega ) & = & 2\pi \displaystyle \sum _{j}\left[{u}_{{xj}}({k}_{y}){u}_{{x}^{\prime }j}^{*}({k}_{y})\delta (\omega -{E}_{j}({k}_{y})/{\rm{\hslash }})\right.\\ & & +\left.{v}_{{xj}}^{*}(-{k}_{y}){v}_{{x}^{\prime }j}(-{k}_{y})\delta (\omega +{E}_{j}(-{k}_{y})/{\rm{\hslash }})\right].\end{eqnarray} \tag{ 64 }$

In section 5, we calculate $\rho (x,x,{k}_{y},\omega )$ to discuss excitations that can be probed at each position x. We note that similar calculations are performed for a fermionic Hofstadter model in various traps by Buchhould et al [60].

4.3. Extended Thomas–Fermi approximation for density profiles

Here we provide an approximate solution to the scaled GP equation (56). This helps us tune the potential and the interaction to obtain the desired density profiles later. We extend the Thomas–Fermi approximation [55, 56], where the condensate wave function is assumed to vary slowly in space, so that the oscillations of f_x between the two sublattices due to ${\rm{\Delta }}\ne 0$ are also taken into account. We assume

$\begin{eqnarray}&&{f}_{x}=\bar{f}(x)+{\epsilon }_{x}\delta f(x),\end{eqnarray} \tag{ 65 }$

where $\bar{f}(x)$ and $\delta f(x)$ are slowly varying real functions, and ${\epsilon }_{x}=+1$ ( $-1$ ) if x belongs to the A (B) sublattice. As can be seen in figure 4, all the sites having the same x belong to the same sublattice. Assuming $| {\rm{\Delta }}| \ll {J}_{1},$ we can expect $| \delta f(x)| \ll \bar{f}(x).$ Then the rhs of equation (56) can be approximated as

$\begin{eqnarray}&&\left[\displaystyle \sum _{{x}^{\prime }=x-3/2}^{x+3/2}{{\mathcal{H}}}_{{{xx}}^{\prime }}(0)\right]\;\bar{f}(x)+\left[\displaystyle \sum _{{x}^{\prime }=x-3/2}^{x+3/2}{\epsilon }_{{x}^{\prime }}{{\mathcal{H}}}_{{{xx}}^{\prime }}(0)\right]\;\delta f(x)+{{Un}}_{y}{[\bar{f}(x)+\delta f(x)]}^{3}\\ &&\quad \approx \left[V(x)-3{J}_{1}-6{J}_{2}\mathrm{cos}{\rm{\Phi }}+{{Un}}_{y}\bar{f}{(x)}^{2}\right]\;\bar{f}(x)\\ &&\quad +{\epsilon }_{x}\left[{\rm{\Delta }}\bar{f}(x)+(V(x)+3{J}_{1}-6{J}_{2}\mathrm{cos}{\rm{\Phi }}+3{{Un}}_{y}\bar{f}{(x)}^{2})\delta f(x)\right],\end{eqnarray} \tag{ 66 }$

where we have ignored higher-order terms in Δ and $\delta f(x).$ Here, the first and second lines of the last expression can be viewed as uniform and staggered components, respectively, because ${\epsilon }_{x}$ oscillates rapidly and other functions of x vary slowly. Requiring that equation (56) hold separately for different components (because they cannot cancel each other), we obtain

$\begin{eqnarray}&&{{Un}}_{y}\bar{f}{(x)}^{2}=\mathrm{max}({{Un}}_{\mathrm{TF}}^{\mathrm{max}}-V(x),0),{{Un}}_{\mathrm{TF}}^{\mathrm{max}}:= \mu +3{J}_{1}+6{J}_{2}\mathrm{cos}{\rm{\Phi }},\end{eqnarray} \tag{ 67a }$

$\begin{eqnarray}&&\delta f(x)=-\displaystyle \frac{{\rm{\Delta }}\bar{f}(x)}{V(x)-\mu +3{J}_{1}-6{J}_{2}\mathrm{cos}{\rm{\Phi }}+3{{Un}}_{y}\bar{f}{(x)}^{2}}.\end{eqnarray} \tag{ 67b }$

Therefore, the density profile scaled by the interaction U is obtained as

$\begin{eqnarray}&&{{Un}}_{y}{| {f}_{x}| }^{2}\approx \mathrm{max}({{Un}}_{\mathrm{TF}}^{\mathrm{max}}-V(x),0)\left(1-\displaystyle \frac{{\epsilon }_{x}{\rm{\Delta }}}{{{Un}}_{\mathrm{TF}}^{\mathrm{max}}-V(x)+3{J}_{1}}\right).\end{eqnarray} \tag{ 68 }$

This expression agrees well with the density profiles of the numerically calculated GP ground states presented in section 5. As seen in this expression, ${n}_{\mathrm{TF}}^{\mathrm{max}}$ can be interpreted as the maximal uniform-component density achieved at the potential minimum with $V(x)=0.$

5. Ground state and excitations in trapped condensates: numerical results

In this section, we present numerical results on the ground state and excitations of trapped condensates based on the formalism described in section 4. We exploit the strip geometry described in section 4.2, and introduce a box trap or a harmonic trap along the x direction. While harmonic traps are used more commonly in experiments of ultracold atoms, a box trap with sharp boundaries has also been realized recently [61]. Sharp boundaries can also be realized in synthetic dimensions [26–28]. We demonstrate that chiral edge states reflecting the nontrivial bulk band topology do appear in a box trap. For a harmonic trap, by contrast, our results show that such edge states are substantially obscured and difficult to observe, as opposed to what is expected from a semiclassical picture.

5.1. Box trap

We first consider the case of a box trap

$\begin{eqnarray}V(x)=\left\{\begin{array}{cc}0 & \left(| x| \lt \displaystyle \frac{3}{8}{N}_{x}d\right);\\ \infty & (\mathrm{otherwise}),\end{array}\right.\end{eqnarray} \tag{ 69 }$

which has sharp boundaries of zigzag type as in figure 4. The extended TF result (68) leads to a uniform density in each sublattice inside the box. Summing equation (68) over the lattice sites $| x| \lt \frac{3}{8}{N}_{x}d,$ we find ${{Un}}_{y}\approx {{Un}}_{\mathrm{TF}}^{\mathrm{max}}{N}_{x}.$ We note that in the BdG calculation for a strip geometry, the product ${{Un}}_{y}$ is an input parameter, i.e., if ${{Un}}_{y}$ is fixed, the result does not depend on individual values of U and n_y. We can thus tune ${{Un}}_{y}$ to obtain a desired scaled average density ${{Un}}_{\mathrm{TF}}^{\mathrm{max}}/{J}_{1}\approx {{Un}}_{y}/({N}_{x}{J}_{1}).$

We set $({J}_{2}{{\rm{e}}}^{{\rm{i}}{\rm{\Phi }}}/{J}_{1},{\rm{\Delta }}/{J}_{2},{{Un}}_{y}/({N}_{x}{J}_{1}))=(0.1i,1.2\times {3}^{3/2},1)$ so that the case of figure 2(c) is realized in a quasi-homogeneous region in the bulk. The scaled density profile ${{Un}}_{y}{| {f}_{x}| }^{2}/{J}_{1}$ of the GP ground state in figure 5 is indeed almost uniform in each sublattice, and agrees well with the extended TF result (68), except over a few sites near each boundary.

**Figure 5.** Density profile (scaled by $U/{J}_{1}$ ) of the ground state along the x direction in a box trap (69) with ${N}_{x}=80$ (see figure 5). The data for the GP ground state are compared with the extended TF result (68). The parameters of the model are chosen such that the case of figure 2(c) with interaction-induced nontrivial topology ${C}_{+}=+1$ is realized in a quasi-homogeneous region in the bulk.
Download figure:
Standard image High-resolution image

**Figure 5.** Density profile (scaled by $U/{J}_{1}$ ) of the ground state along the x direction in a box trap (69) with ${N}_{x}=80$ (see figure 5). The data for the GP ground state are compared with the extended TF result (68). The parameters of the model are chosen such that the case of figure 2(c) with interaction-induced nontrivial topology ${C}_{+}=+1$ is realized in a quasi-homogeneous region in the bulk.
Download figure:
Standard image High-resolution image

We now discuss excitations calculated by the BdG theory described in section 4.2. For an interval I of the x coordinate, we consider the integrated spectral weight ${\displaystyle \sum }_{x\in I}\rho (x,x,{k}_{y},\omega ),$ where $\rho (x,x,{k}_{y},\omega )$ is defined in equation (63). The results for (a) $I=(-30d,-15d)$ and (b) $I=(+15\;d,+30d)$ are presented in figure 6. In both of these results, the continuum of excitations corresponding to the two bulk particle (hole) excitation bands is found with high (low) spectral density. Furthermore, inside the band gaps, chiral modes with negative and positive velocities are clearly formed in (a) and (b), respectively⁷ . These results are consistent with the formation of chiral edge modes propagating in the −y (+y) direction at the left (right) boundary as expected from the bulk topological number ${C}_{+}=+1.$

**Figure 6.** Integrated spectral weight ${\displaystyle \sum }_{x\in I}\rho (x,x,{k}_{y},\omega )$ for (a) x in $I=(-30d,-15d)$ and (b) x in $I=(+15d,+30d),$ which contain the left and right edges, respectively. The BdG calculation is performed using the GP ground state shown in figure 5, and the spectral weight $\rho (x,x,{k}_{y},\omega )$ is calculated with equation (63). Positive (negative) energies correspond to particle (hole) excitations.
Download figure:
Standard image High-resolution image

**Figure 6.** Integrated spectral weight ${\displaystyle \sum }_{x\in I}\rho (x,x,{k}_{y},\omega )$ for (a) x in $I=(-30d,-15d)$ and (b) x in $I=(+15d,+30d),$ which contain the left and right edges, respectively. The BdG calculation is performed using the GP ground state shown in figure 5, and the spectral weight $\rho (x,x,{k}_{y},\omega )$ is calculated with equation (63). Positive (negative) energies correspond to particle (hole) excitations.
Download figure:
Standard image High-resolution image

To detect the chiral edge modes between the two excitation bands, a high-frequency probe is required. This contrasts with the case of fermionic topological insulators, where edge modes cross the Fermi level and can be excited with infinitesimal energies. In ultracold-atom experiments, stimulated Raman transitions can be used to create excitations with desired momentum and frequency. Since edge modes are isolated from bulk modes in momentum and frequency in figure 6, Raman transitions can transfer a portion of a condensate selectively into a particular edge mode with momentum k_y and frequency $\omega ={E}_{j}({k}_{y})/{\rm{\hslash }},$ realizing an 'edge matter wave'⁸ . The resulting condensate wave function is a coherent superposition of the background condensate and the edge matter wave, which is given by

$\begin{eqnarray}&&\psi (x,y,t)\approx \sqrt{{n}_{y}}{{\rm{e}}}^{-{\rm{i}}\mu t/{\rm{\hslash }}}\left[{f}_{x}+\alpha {u}_{{xj}}({k}_{y}){{\rm{e}}}^{{\rm{i}}({k}_{y}y-\omega t)}+{\alpha }^{*}{v}_{{xj}}^{*}({k}_{y}){{\rm{e}}}^{-{\rm{i}}({k}_{y}y-\omega t)}\right],\end{eqnarray} \tag{ 70 }$

where α is the complex amplitude of the edge mode. We assume ${| \alpha | }^{2}\ll 1$ to ensure that the linearization done in equation (50) to derive the BdG theory works. We relate the amplitude α to the microscopic process later. Since the edge mode has its weight mainly around an edge, we can expect that a density wave is formed along the edge as a result of the interference with the background condensate. Using equation (70), the scaled density profile relative to the ground state is calculated as

$\begin{eqnarray}&&\displaystyle \frac{U}{{J}_{1}}\left[{| \psi (x,y,t)| }^{2}-{n}_{y}{| {f}_{x}| }^{2}\right]\approx \displaystyle \frac{{{Un}}_{y}}{{J}_{1}}\left[\alpha {z}_{{xj}}({k}_{y}){{\rm{e}}}^{{\rm{i}}({k}_{y}y-\omega t)}+{\rm{c}}.{\rm{c}}.\right]\end{eqnarray} \tag{ 71 }$

with ${z}_{{xj}}({k}_{y}):= {f}_{x}^{*}{u}_{{xj}}({k}_{y})+{f}_{x}{v}_{{xj}}({k}_{y}).$ We shift the origin of t such that α becomes real. We plot equation (71) for different times in figure 7. A density wave is indeed formed along the right edge, and it propagates in the negative y direction with the phase velocity ${E}_{40}({k}_{y})/({k}_{y}{J}_{1})=-1.83\times \sqrt{3}d$ of the edge mode. We expect that such a propagating density wave can be used as a macroscopically enhanced experimental signature of an edge mode.

**Figure 7.** Interference patterns of an edge matter wave with the background condensate at different times ${J}_{1}t/{\rm{\hslash }}=0,0.5,1,1.5,2.$ The color shows the scaled density relative to that of the ground state, $\frac{U}{{J}_{1}\alpha }\left[{| \psi (x,y,t)| }^{2}-{n}_{y}{| {f}_{x}| }^{2}\right],$ calculated from equation (71). The center of each triangular pixel corresponds to a site of the honeycomb lattice. The edge matter wave is created by transferring a portion of the condensate into the edge mode with the momentum ${k}_{y}\left(\sqrt{3}d\right)/(2\pi )=-0.35$ (equivalent to 0.65 due to the periodicity of the Brillouin zone) and the energy ${E}_{40}({k}_{y})/{J}_{1}=4.03$ as can be seen from figure 6(b). During the time evolution (a)–(e), the pattern propagates in the (negative) y direction along the right edge with the phase velocity ${E}_{40}({k}_{y})/({k}_{y}{J}_{1})=-1.83\times \sqrt{3}d;$ see, e.g., the propagation of an antinode with a positive variation (red) as indicated by the arrows.
Download figure:
Standard image High-resolution image

In the above argument, we keep the amplitude α of the edge matter wave undetermined. Here we determine α based on the microscopic process. This would help design an experimental setup for creating and observing an edge matter wave. A pair of Raman lasers with wave vectors ${{\boldsymbol{k}}}_{\mathrm{1,2}}$ and frequencies ${\omega }_{\mathrm{1,2}}$ are prepared in such a manner that ${\boldsymbol{k}}={{\boldsymbol{k}}}_{1}-{{\boldsymbol{k}}}_{2}$ points in the y direction [i.e., ${\boldsymbol{k}}=(0,{k}_{y},0)$ ], and that ${\rm{\hslash }}\omega ={\rm{\hslash }}({\omega }_{1}-{\omega }_{2})$ is resonant with the excitation energy ${E}_{j}({k}_{y})$ of the edge mode. These lasers bring about the following time-dependent perturbation:

$\begin{eqnarray}&&{H}^{\prime }(t)=\displaystyle \sum _{(x,y)}{\rm{\hslash }}{\rm{\Omega }}\mathrm{cos}({k}_{y}y-\omega t){a}^{\dagger }(x,y)a(x,y),\end{eqnarray} \tag{ 72 }$

where the sum runs over positions $(x,y)$ of lattice sites, and Ω describes the strength of the atom–light coupling. This perturbation adds the following term to the rhs of the linearized GP equation (50):

$\begin{eqnarray}&&{\rm{\hslash }}{\rm{\Omega }}\mathrm{cos}({k}_{y}y-\omega t)\psi (x,y,t)\approx {\rm{\hslash }}{\rm{\Omega }}\mathrm{cos}({k}_{y}y-\omega t)\sqrt{{n}_{y}}{f}_{x}{{\rm{e}}}^{-{\rm{i}}\mu t/{\rm{\hslash }}}.\end{eqnarray} \tag{ 73 }$

This describes a transfer of the condensate particles to the target edge mode; higher-order processes in which particles are further transferred to higher-energy states are neglected. In the presence of this term, we assume a solution of the form

$\begin{eqnarray}&&\phi (x,y,t)=\sqrt{{n}_{y}}{{\rm{e}}}^{-{\rm{i}}\mu t/{\rm{\hslash }}}\left[\alpha (t){u}_{{xj}}({k}_{y}){{\rm{e}}}^{{\rm{i}}({k}_{y}y-\omega t)}+{\alpha }^{*}(t){v}_{{xj}}^{*}({k}_{y}){{\rm{e}}}^{-{\rm{i}}({k}_{y}y-\omega t)}\right].\end{eqnarray} \tag{ 74 }$

The amplitude $\alpha (t)$ of the edge matter wave is found to obey

$\begin{eqnarray}&&\left[{\rm{i}}{\rm{\hslash }}{\partial }_{t}\alpha (t)\right]{\tau }_{3}{{\boldsymbol{w}}}_{j}({k}_{y})=\displaystyle \frac{{\rm{\hslash }}{\rm{\Omega }}}{2}{\boldsymbol{f}},\end{eqnarray} \tag{ 75 }$

where ${{\boldsymbol{w}}}_{j}({k}_{y})={(\{{u}_{{xj}}({k}_{y})\},\{{v}_{{xj}}({k}_{y})\})}^{T}$ as defined in section 4.2 and ${\boldsymbol{f}}:= {(\{{f}_{x}\},\{{f}_{x}^{*}\})}^{T}.$ If Raman lasers are illuminated in the time interval $[0,\delta t],$ the integrated amplitude of the edge matter wave is obtained as

$\begin{eqnarray}&&\alpha (\delta t)=-{\rm{i}}\displaystyle \frac{{\rm{\Omega }}\delta t}{2}{{\boldsymbol{w}}}_{j}^{\dagger }({k}_{y}){\boldsymbol{f}},\end{eqnarray} \tag{ 76 }$

where we have used ${{\boldsymbol{w}}}_{i}^{\dagger }({k}_{y}){\tau }_{3}{{\boldsymbol{w}}}_{j}({k}_{y})={\delta }_{{ij}}$ in solving equation (75). Since the GP ground state ${\boldsymbol{f}}$ is extended over the strip and the edge-mode wave function ${{\boldsymbol{w}}}_{j}({k}_{y})$ is localized only over a few sites around the edge, their overlap scales as ${{\boldsymbol{w}}}_{j}^{\dagger }({k}_{y}){\boldsymbol{f}}\sim 1/\sqrt{{N}_{x}}.$ Using this result, we can tune Ω or $\delta t$ to achieve a desired amplitude α.

5.2. Harmonic trap

We next consider the case of a harmonic trap $V(x)=\frac{\kappa }{2}{x}^{2}.$ The extended TF result (68) suggests that the condensate extends over the range $| x| \lt {R}_{\mathrm{TF}}:= \sqrt{2{{Un}}_{\mathrm{TF}}^{\mathrm{max}}/\kappa }.$ Integrating equation (68) over this range, we find

$\begin{eqnarray}&&{{Un}}_{y}\approx {\displaystyle \int }_{-{R}_{\mathrm{TF}}}^{{R}_{\mathrm{TF}}}\displaystyle \frac{4}{3d}{\rm{d}}x({R}_{\mathrm{TF}}^{2}-{x}^{2})=\displaystyle \frac{8\kappa {R}_{\mathrm{TF}}^{3}}{9d}=\displaystyle \frac{16{R}_{\mathrm{TF}}}{9d}{{Un}}_{\mathrm{TF}}^{\mathrm{max}}.\end{eqnarray} \tag{ 77 }$

To achieve the desired values of ${{Un}}_{\mathrm{TF}}^{\mathrm{max}}$ and ${R}_{\mathrm{TF}},$ we can thus set the input parameters as

$\begin{eqnarray}&&\displaystyle \frac{\kappa }{2}=\displaystyle \frac{{{Un}}_{\mathrm{TF}}^{\mathrm{max}}}{{R}_{\mathrm{TF}}^{2}},{{Un}}_{y}=\displaystyle \frac{16{R}_{\mathrm{TF}}}{9d}{{Un}}_{\mathrm{TF}}^{\mathrm{max}}.\end{eqnarray} \tag{ 78 }$

We set $({J}_{2}{{\rm{e}}}^{{\rm{i}}{\rm{\Phi }}}/{J}_{1},{\rm{\Delta }}/{J}_{2},{{Un}}_{\mathrm{TF}}^{\mathrm{max}}/{J}_{1})=(0.1i,1.2\times {3}^{3/2},1)$ so that the case of figure 2(c) (with nontrivial topology ${C}_{+}=+1$ ) is realized around the center of the trap. The scaled density profile ${{Un}}_{y}{| {f}_{x}| }^{2}/{J}_{1}$ of the GP ground state in figure 8(a) indeed shows the maximum of near unity at the center; the oscillating pattern of the profile agrees well with the extended TF result (68). As we move away from the center, the density decreases towards zero, and the case of figure 2(a) (with trivial topology ${C}_{+}=0$ ) is expected to be realized for $| x| \gt {R}_{\mathrm{TF}}.$ This indicates that topological boundaries appear inside the condensate, assuming local homogeneity as in the semiclassical approach. To locate such boundaries, we plot in figure 8(b) the local band edges ${E}_{X}(-{\boldsymbol{K}})$ with $X=A,B,$ which are calculated by substituting the density profile of figure 8(a) into equation (34). The two energies ${E}_{A/B}(-{\boldsymbol{K}})$ indeed crosses at $x/d\approx -26$ and 24.5, where topological boundaries are expected to be formed.

**Figure 8.** (a) Density profile (scaled by $U/{J}_{1}$ ) of the ground state in a harmonic trap. The data for the GP ground state are compared with the extended TF result (68). The parameters of the model are chosen so that the case of figure 2(c) is realized around the center of the trap. Calculations are performed in a sufficiently wide strip (with ${N}_{x}=160$ ) so that the effect of the strip edges is negligible. (b) Local band edges ${E}_{A/B}(-{\boldsymbol{K}})$ calculated from the density profile in (a). For each $x\in X,$ ${E}_{X}(-{\boldsymbol{K}})$ is calculated by substituting the local density ${{Un}}_{y}{| {f}_{x}| }^{2}$ and the local potential $V(x)-\mu$ into $2{{Unf}}_{X}^{2}$ and −μ, respectively, in equation (34); ${E}_{\bar{X}}(-{\boldsymbol{K}})$ is calculated with a similar procedure using the average density at the two neighboring sites.
Download figure:
Standard image High-resolution image

**Figure 8.** (a) Density profile (scaled by $U/{J}_{1}$ ) of the ground state in a harmonic trap. The data for the GP ground state are compared with the extended TF result (68). The parameters of the model are chosen so that the case of figure 2(c) is realized around the center of the trap. Calculations are performed in a sufficiently wide strip (with ${N}_{x}=160$ ) so that the effect of the strip edges is negligible. (b) Local band edges ${E}_{A/B}(-{\boldsymbol{K}})$ calculated from the density profile in (a). For each $x\in X,$ ${E}_{X}(-{\boldsymbol{K}})$ is calculated by substituting the local density ${{Un}}_{y}{| {f}_{x}| }^{2}$ and the local potential $V(x)-\mu$ into $2{{Unf}}_{X}^{2}$ and −μ, respectively, in equation (34); ${E}_{\bar{X}}(-{\boldsymbol{K}})$ is calculated with a similar procedure using the average density at the two neighboring sites.
Download figure:
Standard image High-resolution image

To see whether edge modes of the topological origin appear at such boundaries, we plot in figure 9(a) the local spectral weight $\rho (x,x,\omega )={\displaystyle \sum }_{{k}_{y}}\rho (x,x,{k}_{y},\omega ).$ The formation of an excitation gap with the vanishing spectral weight is seen around the center x = 0. As we move away from the center, the gap gradually closes as expected from the semiclassical approach. However, the reopening of the gap as in figure 8(b) is not seen in this figure. This is more clearly seen in the integrated spectral weight for small intervals shown in figure 9(b). For I = (0, 3d), the formation of a gap with the vanishing spectral weight can be seen. For $I=(24d,27d)$ and $(30d,33d),$ by contrast, the spectral weight has a tiny but non-vanishing value at the valley, indicating a gapless nature. This indicates the breakdown of the semiclassical approach. A possible reason for it is as follows. The local band edges calculated with the semiclassical approach in figure 8(b) indicate that not only the size of the gap but also its location (in energy) changes as a function of the position x. Therefore, beyond the semiclassical picture, the energy gap for a particular position x is easily penetrated by the states in the same energy range in the surrounding region. The observation of edge modes in a harmonic trap thus remains a challenging issue.

**Figure 9.** (a) Local spectral weight $\rho (x,x,\omega )={\displaystyle \sum }_{{k}_{y}}\rho (x,x,{k}_{y},\omega ).$ A BdG calculation is performed using the GP ground state displayed in figure 8, and the spectral weight $\rho (x,x,{k}_{y},\omega )$ is calculated with equation (63). (b) Integrated local spectral weight ${\displaystyle \sum }_{x\in I}\rho (x,x,\omega )$ for small intervals I = (0, 3d), (24d, 27d), and (30d, 33d).
Download figure:
Standard image High-resolution image

6. Summary and outlook

We have studied the topological properties of Bogoliubov excitation bands in BECs in optical lattices on the basis of a Bose–Hubbard extension of the Haldane model. We have shown that the topological properties of the Bloch bands in the noninteracting case are smoothly carried over to those of Bogoliubov excitation bands in the interacting case, and that the parameter regions showing nontrivial topology enlarge with increasing the Hubbard interaction or the particle density. In the presence of sharp boundaries, chiral edge modes appear in the gap between topologically nontrivial excitation bands. We propose that Raman transitions can be used to coherently transfer a portion of the condensate into an edge mode, and that a density wave is formed along the edge due to an interference with the background condensate. This can be used as a macroscopically enhanced experimental signature of the edge mode. By contrast, our results for a harmonic trap show that edge states are substantially obscured and difficult to observe, as opposed to what is expected from a semiclassical picture.

Our study demonstrates that BECs in optical lattices offer a unique playground in the studies of band topology. The macroscopic nature of BECs can enhance signatures of a topological edge mode. The high controllability of BECs allows various methods of exciting particles to such edge modes. While we have considered Raman transitions in this paper, a trap quench can provide another useful method. While both the bulk and edge modes are excited by such a quench, the edge excitations may exhibit a distinct time evolution because of their chiral nature (see [63] for related numerical demonstrations for fermions). It will also be interesting to exploit the high controllability of optical lattices to design bosonic systems with different symmetries or dimensionality, where different topological classes can be explored as expected from the studies of fermions [17, 18].

Note added.—Recently, we became aware of two independent works [64, 65], where nontrivial topology of Bogoliubov excitation bands and associated edge states were discussed in different bosonic systems. Engelhardt and Brandes [64] have considered a one-dimensional system with inversion symmetry, while Bardyn et al [65] have considered a kagome vortex lattice with potential realization in nonlinear optical systems or exciton-polariton condensates.

Acknowledgments

The authors thank Yusuke Horinouchi and Ryuichi Shindou for useful discussions. SF acknowledges the enlightening lecture by Shuichi Murakami on Berry Phase Physics and Topological Insulators (Department of Physics, University of Tokyo, 2013), from which this work was partly motivated. This work was supported by KAKENHI Grant Nos. 25800225 and 26287088 from the Japan Society for the Promotion of Science, a Grant-in-Aid for Scientific Research on Innovative Areas 'Topological Materials Science' (KAKENHI Grant No. 15H05855), the Photon Frontier Network Program from MEXT of Japan, and the Mitsubishi Foundation.