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

We show that Bose-Einstein condensates in optical lattices with broken time-reversal symmetry can support chiral edge modes originating from nontrivial bulk excitation band topology. To be specific, we analyze a Bose-Hubbard extension of the Haldane model, which can be realized with recently developed techniques of manipulating honeycomb optical lattices. The topological properties of Bloch bands known for the noninteracting case are smoothly carried over to Bogoliubov excitation bands for the interacting case. We show that the parameter ranges that display topological bands enlarge with increasing the Hubbard interaction or the particle density. In the presence of sharp boundaries, chiral edge modes appear in the gap between topological excitation bands. We demonstrate that by coherently transferring a portion of a condensate into an edge mode, a density wave is formed along the edge owing to an interference with the background condensate. This offers a unique method of detecting an edge mode through a macroscopic quantum phenomenon.


Introduction
Topological insulators and superconductors have attracted great 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,2,3]. A prototype of topological insulators is an integer quantum Hall system [4], which exhibits a quantized Hall conductivity σ 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 > 0 (C < 0) [6,7]. Here a nonzero value of C is caused by the breaking of timereversal symmetry due to a magnetic field. The discovery of Z 2 topological insulators [8,9,10,11,12,13,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,23,24,25], and synthetic dimensions [26,27,28]. Furthermore, the Haldane model [29], in which nonuniform fluxes pierce through the system, has been realized using fermionic atoms in a periodically modulated honeycomb optical lattice [30] (see Refs. [31,32,33,34,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 ± 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 Ref. [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 (BEC). 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 previously in photonic [37,38,39,40,41], phononic [42,43,44], magnonic [45,46,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,40,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 and a potential combination 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 Fig. 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 a 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-1 2 fermionic Haldane-Hubbard model in Refs. [50,51,52]. The rest of the paper is organized as follows. In Sec. 2, we describe our model, and review the band structure in the noninteracting case. In Sec. 3, we present a Bogoliubov theory for homogeneous condensates with weak repulsive interactions, and analyze the topology of Bogoliubov excitation bands. In Secs. 4 and 5, we analyze the ground state and excitations of inhomogeneous condensates by using a Bogoliubov-de Gennes theory. After describing the basic formalism in Sec. 4, we present numerical results for box and harmonic traps in Sec. 5. In Sec. 6, we present a summary of this paper and discuss an outlook for future studies.

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.

Bose-Hubbard Hamiltonian
We consider bosonic atoms in an optical lattice, which are well described in the tightbinding limit by a Bose-Hubbard model. The Hamiltonian of our system is given by where r and r ′ run over all the site positions of the lattice, a(r) is the bosonic annihilation operator at the site r, and U describes the on-site Hubbard interaction. The diagonal element J(r, r) gives a potential energy at the site r, and the off-diagonal element J(r, r ′ ) with r = r ′ describes the (generally complex) hopping amplitude between the two sites and satisfies J(r ′ , r) = J * (r, r ′ ) in order for H to be hermitian. We set the total number of particles to N: We focus on the case in which the hopping terms in Eq. (1), which will be denoted by H kin , are given by the Haldane model in a honeycomb lattice [29]; see Fig. 1. A honeycomb lattice consists of A and B sublattices, and is spanned by the primitive vectors a 1,2 (a 3 is also introduced for convenience). We also introduce the vectors δ 1,2,3 , which are directed from a B site to the three neighboring A sites. These vectors are given by where d is the length between the neighboring sites. We also introduce the reciprocal lattice vectors which satisfy a i · a * j = 2πδ ij (i, j = 1, 2). The Haldane model consists of the real nearest-neighbor hopping J 1 , the complex next-nearest-neighbor hopping J 2 e ±iΦ , and the potential difference 2∆ between the two sublattices. Nonzero values of J(r, r ′ ) are thus given by where j = 1, 2, 3, and ǫ A,B = ±1. Here, r ∈ X indicates that the site r belongs to the X sublattice, and V (r) describes an external potential which depends on the setting of our system. We assume J 1 , J 2 > 0 in the following.

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 where the sum is taken over the discrete momenta k in the first Brillouin zone, and N 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 Here the 2 × 2 hermitian matrix H(k) can be written in the form where I is the identity matrix, σ = (σ 1 , σ 2 , σ 3 ) are the Pauli matrices. The coefficients h 0 (k) and h(k) = (h 1 (k), h 2 (k), h 3 (k)) are calculated as The two energy bands are obtained through the diagonalization of Eq. (8) as An example of energy bands is presented in Fig. 2(a). For noninteracting fermions, complete filling of the lower band leads to a band insulator. For noninteracting bosons, Bose-Einstein condensation into the lowest-energy single-particle state occurs at zero temperature. ‡ When J 2 = ∆ = 0, the bottom of ‡ In the thermodynamic limit, Bose-Einstein condensation does not occur at finite temperatures in two dimensions. In finite-size systems, however, a large condensate fraction can still be achieved if coherence is formed over the system at sufficiently low temperatures. the lower band e − (k) is located at k = 0. To find whether the position of the bottom can change owing to finite J 2 or ∆, we expand Eq. (10) around k = 0 as Therefore, e − (k) is minimized at k = 0 if the following condition is met: This condition is satisfied in the regime J 2 , |∆| ≪ J 1 , which is relevant to the Haldane model realized in the scheme of Ref. [30]. § To determine the single-particle ground state, we parametrize h(0) using the polar coordinates as h(0) = (−3J 1 , 0, ∆) = h(0)(sin θ 0 cos π, sin θ 0 sin π, cos θ 0 ).
For noninteracting bosons, Bose-Einstein condensation occurs in the mode created by a † − (0). For interacting bosons, the condensate wave function is gradually modified with increasing the interaction, as discussed in the next section.

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 > 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. § When ∆ = 0 and Φ = π, the condition in Eq. (12) is written simply as J 1 > 6J 2 . By relating phases of bosons to angles of classical XY spins, we find that this condition is equivalent to the known stability condition of the ferromagnetic order in a honeycomb lattice magnet with competing ferromagnetic 2J 1 and antiferromagnetic −2J 2 couplings [53, 54].

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(r), a † (r)) by (ψ(r), ψ * (r)) in the Hamiltonian (1), and minimize it with respect to (ψ(r), ψ * (r)) under the constraint r |ψ(r)| 2 = N. Since the single-particle ground state is formed at k = 0 (as discussed in Sec. 2.2), we introduce the following homogeneous ansatz for the interacting case: We also introduce the chemical potential µ as a Lagrangian multiplier to satisfy the particle-number constraint. The functional to be minimized is then given by Minimizing this with respect to ψ * X (X = A, B) gives a homogeneous version of the GP equations: Since the single-particle ground state is created by a † − (0) in Eq. (14), it is convenient to parametrize (ψ A , ψ B ) T as where n := N/(2N uc ) is the average number of particles per site. The parameter θ is determined by solving Eq. (20b); the chemical potential µ is determined by substituting the obtained θ into the LHS of (20a). For the obtained θ, we consider the unitary transformation Bose-Einstein condensation occurs in the lower-band mode created by a † − . When ∆ = 0, the two sublattices are equivalent, and thus Eq. (20b) gives θ = π/2 [49]. For |∆| ≪ J 1 , we can expand Eq. (20) in terms of θ − π/2, obtaining As seen in Eq. (22a), the potential difference 2∆ induces a density imbalance between the two sublattices; this imbalance is reduced by a repulsive interaction U > 0. We now discuss excitations from the condensate ground state by using the Bogoliubov theory. To this end, using Eqs. (6) and (21), we decompose a(r) into the condensate and noncondensate parts as Here, the noncondensate part is given bỹ withĀ = B andB = A. Following the Bogoliubov approximation, we replace both a − and a † − by √ N, substitute Eq. (23) into H − µN, and expand H − µN up to quadratic order inã(r). The terms linear inã(r) orã † (r) disappear because of the stability condition of the condensate in Eq. (20b), and we arrive at the Bogoliubov Hamiltonian with Here, we have introduced the 4 × 4 matrix M(k) and the 2 × 2 matrix M + as To diagonalize the Bogoliubov Hamiltonian (25), we perform generalized Bogoliubov transformations with  38), and indicated for U n/J 1 = 0, 0.5, 1. The solid curves correspond to the noninteracting case [29]. The regions with nontrivial topology (C + = 0) enlarge as U n/J 1 increases.
Here, W (k) and W + are paraunitary matrices which satisfy with τ 3 := diag(1, 1, −1, −1). These equations ensure the invariance of the bosonic commutation relations. If the matrices W (k) and W + are chosen to satisfy the Bogoliubov Hamiltonian (25) is diagonalized as The paraunitary matrix W (k) satisfying Eq. (32) can be constructed numerically by using the method described in Refs. [45,57]. Examples of the calculated Bogoliubov excitation bands E ± (k) are presented in Fig. 2(b,c).

Band gap
Before discussing the topology of the Bogoliubov excitation bands, let us analyze the gap between the two bands E ± (k). The band topology cannot change as far as the band gap remains open. Thus the closing of the gap can signal a change in topology. As seen in Fig. 2 and known in the noninteracting case [29], the smallest gap is found at one of the two points k = ±K := ± (a * 1 + a * 2 ) /3 in the Brillouin zone; see the K ± points in Fig. 1. At these points, because of h 1 (±K) = h 2 (±K) = 0, the matrix M(±K) in Eq. (27) is decoupled into A and B sublattice blocks, each of which can be diagonalized by a standard 2 × 2 Bogoliubov transformation. The excitation energies at k = ±K are then calculated to be The higher (lower) of these energies gives E + (±K) (E − (±K)) in Eq. (33). The gapclosing conditions at k = ±K are thus obtained as For J 2 , |∆| ≪ J 1 , by using (22), we can expand λ A − λ B as Equation (36) is then rewritten into a simple form For Un = 0, this reduces to the exact phase boundaries ∆ = ∓3 √ 3J 2 sin Φ in the noninteracting case [29]. Equation (38) indicates that with increasing the strength of interaction Un/J 1 , these boundaries are shifted to larger |∆| by a factor of G(Un/J 1 ); see Fig. 3. Here, G(s) is an increasing function of s, and expanded for |s| ≪ 1 as Across the gap-closing lines, the topology of the Bogoliubov excitation bands changes, as we discuss next.

Chern number
We now analyze the topology of the Bogoliubov excitation bands. We first note that technically, M(k) in Eq. (27) (more specifically, H(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) was based on the real-space positions r -while this treatment was useful in making H(k) possess the C 3 symmetry of the original lattice, the spacing between the two sublattices introduced an additional phase factor, which broke the periodicity in the first Brillouin zone. To recover the periodicity, we defineM(k) by replacing H(k) byH(k) = e −ik·δ 3 (I−σ 3 )/2 H(k)e ik·δ 3 (I−σ 3 )/2 in Eq. (27). We then introduce the paraunitary matrix W (k) as the matrix which "diagonalizes"M(k) in the sense of Eq. (32).
The 4 × 4 paraunitary matrixW (k) can be parametrized as with To discuss the topology of each excitation band, we introduce the vectors where γ = + and − correspond respectively to the first and second columns of W (k). It follows from Eqs. (31) and (32) that these vectors satisfy the eigen equatioñ and the orthonormality condition For each band, we introduce the Berry curvature [37,45] with ∂ i := ∂ ∂k i . The Chern number can then be introduced as [5, 37, 45] In the noninteracting case Un = 0, both C ± are quantized to integers, and satisfy the zero sum rule γ C γ = 0. In the interacting case Un > 0, however, |w − (k) is not defined at k = 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 Ref. [58]. The "phase diagram" based on the Chern number C + is presented in Fig. 3 (we note that this diagram is not based on ground-state transitions). We numerically confirmed that the boundaries between regions of different C + values are given precisely by the gapclosing condition (36) (or Eq. (38) for J 2 , |∆| ≪ J 1 ) obtained in Sec. 3

.2. The obtained
This is because the Bogoliubov excitations consist only of modes orthogonal to the GP ground state [55,56]; see Eq. (24). ¶ This does not contradict the zero sum rule. In the interacting case U n > 0, the rule applies to the sum over all the particle and hole bands. Namely, C + + C − + C ′ + + C ′ − = 0, where C ′ ± are the Chern numbers associated with the hole bands [45]. One can easily show C + + C ′ + = 0 = C − + C ′ − , and thus the sum rule is trivially satisfied. Therefore, one cannot use the ill-defined nature of C − to change C + to an arbitrary value. 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 > 0, and that the regions displaying nontrivial topology C + = 0 enlarge with increasing Un/J 1 . When C + = 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 Sec. 5. In closing this section, two remarks are in order. The first remark is on the reason why the ranges 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 λ A and λ B defined in Eq. (35), and the transition occurs when |λ A − λ B | = 6 √ 3J 2 | sin Φ|. The interaction U > 0 has the effect of obscuring this difference [through Unf 2 X in Eq. (35)], and thus a larger |∆| 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 < 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 Eq. (38); these regions gradually shrink with increasing |U|n for U < 0. In these regions, in-gap chiral edge modes discussed in Sec. 5 are also expected to be formed.

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 Bogoliubov-de Gennes (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 Sec. 5.

Bogoliubov-de Gennes theory for inhomogeneous condensates
The BdG theory for inhomogeneous condensates can be derived by linearizing a timedependent GP equation. We start from the Heisenberg equation of motion for the time-dependent operator a(r, t): Replacing (a(r, t), a † (r, t)) by classical fields (ψ(r, t), ψ * (r, t)), we obtain the GP equation Equation (2) imposes the normalization condition r |ψ(r, t)| 2 = N. Inserting a stationary state ansatz ψ(r, t) = ψ(r)e −iµt/ into Eq. (48), we obtain the timeindependent GP equation The solution ψ(r) with the lowest frequency µ/ gives the GP ground state.
We now discuss small fluctuations around the GP ground state: ψ(r, t) = ψ(r)e −iµt/ + φ(r, t). Expanding the GP equation (48) to first order in φ(r), we obtain a linear differential equation Assuming a solution of the form φ(r, t) = e −iµt/ u(r)e −iEt/ + v * (r)e iEt/ , we obtain Bogoliubov-de Gennes (BdG) equations If the condensate is stable, these equations admit N s − 1 (linearly independent) sets of solutions (u j (r), v j (r)) with positive frequencies E j / , where N s is the total number of lattice sites. Such solutions can be chosen to satisfy the orthonormality condition 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 Sec. 3. Indeed, substituting the ansatz (16) into the GP equation (49) reproduces Eq. (18). Furthermore, substituting into the BdG equation (52) reproduces Eq. (43). It is known that the operator formulation for the inhomogeneous case also leads to the same set of BdG equations as in Eq. (52) [55,56].

Strip geometry
We apply the BdG theory to analyze excitations in trapped condensates. We describe site positions by two-dimensional coordinates r = (x, y). For simplicity, we consider a strip geometry, which is periodic only along the y direction as in Fig. 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 (Fig. 4) or a harmonic trap [60] in Sec. 5; however, we do not assume a specific trap potential in this subsection. Exploiting the translation invariance along the y direction, we make the following ansatz for the GP ground state: where N y is the number of unit cells, and ψ x = √ N f x satisfies the normalization condition x |ψ x | 2 = N x |f x | 2 = N. Substituting this into the GP equation (49), Here we have introduced where (x ′ , y ′ ) 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 RHS of Eq. (57) does not depend on y ′ . An accurate solution to Eq. (56) with the lowest frequency µ/ can be obtained by numerically performing the imaginary time evolution with the replacement µ → −∂ τ . 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: Substituting these into the BdG equation (52), we obtain an eigen equation Here with F := diag({f x }). The eigen equation (59) can be solved numerically by using the method of Ref. [57]. We denote positive-frequency solutions by 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 where the new bosonic operators {b j (k y )} are related to the original ones as a(x, y) = ψ(x, y) To discuss excitations in a strip geometry, it is useful to introduce the spectral weight at zero temperature where (x, y) and (x ′ , y ′ ) are taken similarly as in Eq. (57), and the average · is taken over the vacuum of the Bogoliubov excitations in Eq. (61). Using Eqs. (61) and (62), Eq. (63) can be rewritten as In Sec. 5, we calculate ρ(x, x, k y , ω) 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].

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 ∆ = 0 are also taken into account. We assume wheref (x) and δf (x) are slowly varying real functions, and ǫ x = +1 (−1) if x belongs to the A (B) sublattice. As can be seen in Fig. 4, all the sites having the same x belong to the same sublattice. Assuming |∆| ≪ J 1 , we can expect |δf (x)| ≪f (x). Then the RHS of Eq. (56) can be approximated as where we have ignored higher-order terms in ∆ and δf (x). Here, the first and second lines of the last expression can be viewed as uniform and staggered components because ǫ x oscillates rapidly and other functions of x vary slowly. Requiring that Eq. (56) holds separately for different components (because they cannot cancel each other), we obtain (67b) Therefore, the density profile scaled by the interaction U is obtained as This expression agrees well with the density profiles of the numerically calculated GP ground states presented in Sec. 5. As seen in this expression, n max TF can be interpreted as the maximal uniform-component density achieved at the potential minimum with V (x) = 0.

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 Sec. 4. We exploit the strip geometry described in Sec. 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,27,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 an expectation from a semiclassical picture.

Box trap
We first consider the case of a box trap which has sharp boundaries of zigzag type as in Fig. 4. The extended TF result (68) leads to a uniform density in each sublattice inside the box. Summing Eq. (68) over the lattice sites |x| < 3 8 N x d, we find Un y ≈ Un max TF 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 max TF /J 1 ≈ Un y /(N x J 1 ). We set (J 2 e iΦ /J 1 , ∆/J 2 , Un y /(N x J 1 )) = (0.1i, 1.2 × 3 3/2 , 1) so that the case of Fig. 2(c) is realized in a quasi-homogeneous region in the bulk. The scaled density  Figure 6. Integrated spectral weight x∈I ρ(x, x, k y , ω) 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 Fig. 5, and the spectral weight ρ(x, x, k y , ω) is calculated with Eq. (63). Positive (negative) energies correspond to particle (hole) excitations.
profile Un y |f x | 2 /J 1 of the GP ground state in Fig. 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.
We now discuss excitations calculated by the BdG theory described in Sec 4.2. For an interval I of the x coordinate, we consider the integrated spectral weight x∈I ρ(x, x, k y , ω), where ρ(x, x, k y , ω) is defined in Eq. (63). The results for (a) I = (−30d, −15d) and (b) I = (+15d, +30d) are presented in Fig. 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.
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 Fig. 6, Raman transitions can transfer a + In Fig. 6 [and Fig. 2(c)], the bulk band gap is relatively small because we examine the case in which nontrivial topology C + = +1 is induced by the interaction U . That is, the system is located in a narrow region between the lines for U n/J 1 = 0 and 1 in Fig. 3, where a small band gap closes and opens again as we change U n/J 1 between these values. If the system is located deep inside a region with nontrivial topology, a much larger band gap can be created, and edge modes can then be distinguished more clearly from bulk modes. Interference patterns of an edge matter wave with the background condensate at different times J 1 t/ = 0, 0.5, 1, 1.5, 2. The color shows the scaled density relative to that of the ground state, U J1α |ψ(x, y, t)| 2 − n y |f x | 2 , calculated from Eq. (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 ( √ 3d)/(2π) = −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 Fig. 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× √ 3d; see, e.g., the propagation of an antinode with a positive variation (red) as indicated by the arrows.
portion of a condensate selectively into a particular edge mode with momentum k y and frequency ω = E j (k y )/ , 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 where α is the complex amplitude of the edge mode. We assume |α| 2 ≪ 1 to ensure that the linearization done in Eq. (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 Eq. (70), the scaled density profile relative to the ground state is calculated as U J 1 |ψ(x, y, t)| 2 − n y |f x | 2 ≈ Un y J 1 αz xj (k y )e i(kyy−ωt) + c.c.
(71) * While the edge matter wave has an infinite lifetime within the BdG theory, it can acquire a finite lifetime due to collisions between quasiparticles and condensed particles (known as Beliaev damping [62]). The estimation of this lifetime however requires a detailed analysis of the collision processes, which is beyond the scope of the present paper.
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 Eq. (71) for different times in Fig. 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 × √ 3d 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.
In the above argument, we have kept 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 k 1,2 and frequencies ω 1,2 are prepared in such a manner that k = k 1 − k 2 points in the y direction [i.e., k = (0, k y , 0)], and that ω = (ω 1 − ω 2 ) is resonant with the excitation energy E j (k y ) of the edge mode. These lasers induce the following time-dependent perturbation to the Hamiltonian: 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): Ω cos(k y y − ωt)ψ(x, y, t) ≈ Ω cos(k y y − ωt) √ n y f x e −iµt/ .
This describes a transfer of the condensate particles to the target edge mode; higherorder 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 φ(x, y, t) = √ n y e −iµt/ α(t)u xj (k y )e i(kyy−ωt) + α * (t)v * xj (k y )e −i(kyy−ωt) .
The amplitude α(t) of the edge matter wave is found to obey where we have used w † i (k y )τ 3 w j (k y ) = δ ij in solving Eq. (75). Since the GP ground state f is extended over the strip and the edge-mode wave function w j (k y ) is localized only over a few sites around the edge, their overlap scales as w † j (k y )f ∼ 1/ √ N x . Using this result, we can tune Ω or δt to achieve a desired amplitude α. For each x ∈ X, E X (−K) is calculated by substituting the local density U n y |f x | 2 and the local potential V (x) − µ into 2U nf 2 X and −µ, respectively, in Eq. (34); EX (−K) is calculated with a similar procedure using the average density at the two neighboring sites.

Harmonic trap
We next consider the case of a harmonic trap V (x) = κ 2 x 2 . The extended TF result (68) suggests that the condensate extends over the range |x| < R TF := 2Un max TF /κ. Integrating Eq. (68) over this range, we find To achieve the desired values of Un max TF and R TF , we can thus set the input parameters as We set (J 2 e iΦ /J 1 , ∆/J 2 , Un max TF /J 1 ) = (0.1i, 1.2×3 3/2 , 1) so that the case of Fig. 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 Fig. 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 Fig. 2(a) (with trivial topology C + = 0) is expected to be realized for |x| > R 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 Fig. 8(b) the local band edges E X (−K) with X = A, B, which are calculated by substituting the density profile of Fig. 8(a) Figure 9.
To see whether edge modes of the topological origin appear at such boundaries, we plot in Fig. 9(a) the local spectral weight ρ(x, x, ω) = ky ρ(x, x, k y , ω). 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 Fig. 8(b) is not seen in this figure. This is more clearly seen in the integrated spectral weight for small intervals shown in Fig. 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 Fig. 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.

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 ranges showing nontrivial topology enlarges 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.
We expect 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 offers 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 Ref. [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.