Topological Edge-State Manifestation of Interacting 2D Condensed Boson-Lattice Systems in a Harmonic Trap

TL;DR

This paper demonstrates that interactions in a 2D bosonic lattice system under a harmonic trap can reveal topological edge states, with a clear band structure and increased edge modes, supported by analytical and numerical results.

Contribution

It shows that weak interactions enable the observation of topological edge states in bosonic systems within harmonic traps, which is typically challenging in noninteracting systems.

Findings

Edge states are clearly visible in the band structure with harmonic trapping.

The number of edge states crossing the gap increases due to interactions.

Analytical expressions for edge-state energies match numerical results.

Abstract

In this Letter, it is shown that interactions can facilitate the emergence of topological edge states of quantum-degenerate bosonic systems in the presence of a harmonic potential. This effect is demonstrated with the concrete model of a hexagonal lattice populated by spin-one bosons under a synthetic gauge field. In fermionic or noninteracting systems, the presence of a harmonic trap can obscure the observation of edge states. For our system with weakly interacting bosons in the Thomas--Fermi regime, we can clearly see a topological band structure with a band gap traversed by edge states. We also find that the number of edge states crossing the gap is increased in the presence of a harmonic trap, and the edge modes experience an energy shift while traversing the first Brillouin zone which is related to the topological properties of the system. We find an analytical expression for the edge-state energies and our comparison with numerical computation shows excellent agreement.

Figures

b of the main text. One can partly reproduce these oscillations. By introducing the particle density $\bar{\rho}_{j}=|\bar{\Psi}_{j}|^{2}$ one can write $\bar{\Psi}_{j}/\bar{\Psi}_{i}=\sqrt{\bar{\rho}_{j}/\bar{\rho}_{i}}\exp(i\Delta\theta_{i,j})$, where $\Delta\theta_{i,j}$ denotes the phase difference of the ground-state wave-function at sites $j$ and $i$. Since $V_{\rm eff}$ and $\bar{\mathcal{H}}_{\rm latt}$ are real the phase difference $\Delta\theta_{i,j}$ should not contribute to Eq. ([17](#A2.E17)). Indeed, due to time-reversal symmetry the ground-state wave-function $\bar{\Psi}_{j}$ ca

Equations

$${\hat{\mathcal{H}}}_{\mathrm{latt}}=-w\sum _{⟨j{j}^{\prime }⟩}{\hat{\mathbf{\Psi }}}_j^{†}{\hat{\mathbf{\Psi }}}_{j^{\prime }}+i\lambda \sum _{⟨⟨j{j}^{\prime }⟩⟩}{\nu }_{j{j}^{\prime }}{\hat{\mathbf{\Psi }}}_j^{†}{S}_z{\hat{\mathbf{\Psi }}}_{j^{\prime }},$$

$${\hat{V}}_{\omega }$$

$${\hat{\mathcal{H}}}_{\mathrm{int}}$$

$$\hat{\mathcal{H}}={\hat{\mathcal{H}}}_{\mathrm{latt}}+{\hat{\mathcal{H}}}_{\mathrm{int}}+{\hat{\mathcal{H}}}_Z+{\hat{V}}_{\omega }-\mu \sum _j{\hat{\rho }}_j,$$

$${\hat{\mathcal{H}}}_{\mathrm{B}}$$

$$V_{\mathrm{eff},j}=M{\omega }^2{x}_j^2/2+{\bar{\mathcal{H}}}_{\mathrm{latt}}+U{\bar{\rho }}_j-\mu .$$

$$\sum _j\left({\mathcal{H}}_{\mathrm{latt}}^{ij}-{\bar{\mathcal{H}}}_{\mathrm{latt}}{\delta }_{ij}+V_{\mathrm{eff},i}{\delta }_{ij}\right){\bar{\Psi }}_j=0.$$

$$\sum _{n^{\prime }}{H}_{k,n-n^{\prime }}{\mathbf{\Phi }}_{k,n^{\prime }}+V_n{\mathbf{\Phi }}_{k,n}=E_k{\mathbf{\Phi }}_{k,n},$$

$$E_k{\mathbf{\Phi }}_{k,\phi }=\left[H_{k,\phi }+V_{\bar{n}}+V_{\bar{n}}^{\prime }(i{\partial }_{\phi }-\bar{n})\right]{\mathbf{\Phi }}_{k,\phi },$$

$$i{\gamma }_{\bar{n}}{\partial }_{\phi }{\tilde{\mathbf{\Phi }}}_{k,\phi }=(E_k-{\tilde{H}}_{k,\phi }-{\tilde{V}}_{\bar{n}}+\text{mathds}1\bar{n}{\gamma }_{\bar{n}}){\tilde{\mathbf{\Phi }}}_{k,\phi }$$

$${\tilde{\mathbf{\Phi }}}_{k,\phi }=e^{-i\left(\bar{n}+\frac{E_k-\mathcal{V}(\bar{n})}{{\gamma }_{\bar{n}}}\right)\Delta \phi }\mathcal{P}{e}^{\frac{i}{{\gamma }_{\bar{n}}}{\int }_{{\phi }_0}^{\phi }d{\phi }^{\prime }{\tilde{H}}_{k,{\phi }^{\prime }}}{\tilde{\mathbf{\Phi }}}_{k,{\phi }_0}.$$

$${\mathbf{\Phi }}_{k,\phi }$$

$$\times \sum _{\nu }{e}^{i{\int }_{{\phi }_0}^{\phi }d{\phi }^{\prime }\left(\frac{{\epsilon }_{k,{\phi }^{\prime }}^{(\nu )}}{{\gamma }_{\bar{n}}}+A_{k,{\phi }^{\prime }}^{(\nu )}-\frac{s{Z}_{k,{\phi }^{\prime }}^{(\nu )}}{2}\right)}\left[{\mathbf{\varphi }}_{k,{\phi }_0}^{(\nu )†}{\mathbf{\Phi }}_{k,{\phi }_0}\right]{\mathbf{\varphi }}_{k,\phi }^{(\nu )}$$

$$E_k^{(\bar{n}\nu )}=$$

$$\frac{1}{2\pi }{\int }_0^{2\pi }d\phi \left({\epsilon }_{k,\phi }^{(\nu )}+{\gamma }_{\bar{n}}{A}_{k,\phi }^{(\nu )}-\frac{s{\gamma }_{\bar{n}}}{2}{Z}_{k,\phi }^{(\nu )}\right).$$

$${\mathcal{C}}_{\nu }=\frac{1}{2\pi }{\int }_0^{2\pi }d\phi {\int }_0^{2\pi }dk{\Omega }_{k,\phi }^{(\nu )},$$

$$\sum _j\left({\mathcal{H}}_{\mathrm{latt}}^{ij}-{\bar{\mathcal{H}}}_{\mathrm{latt}}{\delta }_{ij}+V_{\mathrm{eff},i}{\delta }_{ij}\right){\bar{\Psi }}_j=0,$$

$$V_{\mathrm{eff},j}=V_{\omega ,j}+{\bar{\mathcal{H}}}_{\mathrm{latt}}+U{\bar{\rho }}_j-\mu ,$$

$$V_{\mathrm{eff},i}={\bar{\mathcal{H}}}_{\mathrm{latt}}-\sum _j{\mathcal{H}}_{\mathrm{latt}}^{ij}\frac{{\bar{\Psi }}_j}{{\bar{\Psi }}_i}.$$

$$V_{\mathrm{eff},i}={\bar{\mathcal{H}}}_{\mathrm{latt}}-\sum _j{\mathcal{H}}_{\mathrm{latt}}^{ij}\sqrt{\frac{{\bar{\rho }}_j}{{\bar{\rho }}_i}}.$$

$$V_{\mathrm{eff},i}={\bar{\mathcal{H}}}_{\mathrm{latt}}-\sum _j{\mathcal{H}}_{\mathrm{latt}}^{ij}\sqrt{1-\frac{V_{{\omega }_j}-V_{{\omega }_i}-V_{\mathrm{eff},j}+V_{\mathrm{eff},i}}{\mu -{\bar{\mathcal{H}}}_{\mathrm{latt}}-V_{{\omega }_i}+V_{\mathrm{eff},i}}}.$$

$$H_{k,n-n^{\prime }}=R_k{\delta }_{n-n^{\prime },0}+M_k{\delta }_{n-1,n^{\prime }}+M_k^{†}{\delta }_{n+1,n^{\prime }},$$

$$\displaystyle R_{k}=\left(\begin{array}[]{cc}2\lambda\sin k&-w(1+e^{-ik})\\ -w(1+e^{ik})&-2\lambda\sin k\end{array}\right),$$

$$\displaystyle M_{k}=\left(\begin{array}[]{cc}i\lambda(1-e^{-ik})&0\\ -w&i\lambda(1-e^{-ik})\end{array}\right).$$

$$H_{k,\phi }=d_1{\sigma }_1+d_2{\sigma }_2+d_3{\sigma }_3,$$

$${\varphi }_{k,\phi }^{(\pm )}={\left(\frac{d\pm d_3}{\sqrt{2d(d\pm d_3)}},\pm \frac{d_1+i{d}_2}{\sqrt{2d(d\pm d_3)}}\right)}^{\mathrm{T}}.$$

$$A_{k,\phi }^{(\nu )}=\frac{-w^2\left(1+\cos \phi +\cos (\phi -k)\right)}{2{ϵ}_{k,\phi }^{(\nu )}\left[{ϵ}_{k,\phi }^{(\nu )}-2\lambda \left(\sin \phi -\sin (\phi -k)-\sin k\right)\right]}.$$

$$Z_{k,\phi }^{(\nu )}=d_3(k,\phi )/{ϵ}_{k,\phi }^{(\nu )}.$$

$$s=\frac{{\mathbf{a}}_1^2(\mathbf{\delta }\cdot {\mathbf{a}}_2)-({\mathbf{a}}_1\cdot {\mathbf{a}}_2)(\mathbf{\delta }\cdot {\mathbf{a}}_1)}{{\mathbf{a}}_1^2{\mathbf{a}}_2^2-({\mathbf{a}}_1\cdot {\mathbf{a}}_2{)}^2}.$$

Full text

Optical lattice experiments offer the possibility of simulating atomic crystal structures, creating a clean and well-controlled environment for probing many-body physics concepts and phenomena. In recent years, lattices reassembling the Hofstadter Williams et al. (2010); Aidelsburger et al. (2013); Miyake et al. (2013) and Haldane Jotzu et al. (2014) models were created in optical lattice setups. A central focus in such systems lies in the manifestation of surface states which is a result of the nontrivial bulk topological properties. In topological insulators Hasan and Kane (2010); Qi and Zhang (2011), for instance, electron excitations form a Fermi sea and provided the Fermi level is within the band gap the near-equilibrium dynamics is dominated by the states localized either at the sharp boundaries or at the interface with a system of a different topology. However, unlike solids, ultracold atoms are typically confined by a harmonic trapping potential. The influence of the harmonic trap on the band structure and its topology is therefore of significant importance.

Theoretical studies suggest the presence of a confining potential can modify both the bulk and edge energy spectrum significantly, leading not only to a change in group velocity of edge modes but also to an emergence of additional edge states Stanescu et al. (2010); Buchhold et al. (2012), their disappearance Yan et al. (2015), or localization and a shrinking of the bulk region Buchhold et al. (2012); Kolovsky et al. (2014). Possible ways of overcoming these difficulties include inducing topological interfaces Goldman et al. (2016); Leder et al. (2016) or creating so-called box traps Gaunt et al. (2013); Chomaz et al. (2015); Mukherjee et al. (2017). Creation of such trapping potentials represents a separate challenge and cost for experimental setup. The role of mean field interactions in the Haldane boson model was also considered in Ref. Furukawa and Ueda (2015) where it was shown that the bulk gap can close when the harmonic trap is taken into account. We, however, show that an interacting gas of spin-one bosons prepared in a polar ground state on a two-dimensional lattice can have a clear gap in the energy spectrum of the spin-$\pm 1$ excitations. Furthermore the gap is crossed by edge states that reflect the topological structure of the lattice.

Advances in the creation of synthetic gauge fields for atoms in optical lattices, as well as in photonic crystals, have allowed the investigation of topological bands populated by bosons Goldman et al. (2014); Lu et al. (2014). Of particular relevance is the realization of topological collective excitations in these systems Galilo et al. (2015); Bardyn et al. (2016); Peano et al. (2016). While fermions form a Fermi sea, bosons tend to occupy the lowest energy state available and population of the higher energy edge modes is a challenge. Ways to overcome this challenge were suggested. In particular, quantum quenches provide a tool to selectively induce dynamical instabilities in edge modes of topological bosonic systems Galilo et al. (2015) (quenches in similar fermionic systems were considered in Ref. Caio et al. (2015)). The use of periodic driving has also been proposed Engelhardt et al. (2016); Dehghani and Mitra (2016); Peano et al. (2016). In these studies, idealized (e.g. open) boundary conditions were adopted for simplicity.

In this Letter, we show that pristine edge modes, much like those occurring in systems with idealized open boundary conditions, can occur in certain harmonically trapped interacting systems. In particular, interactions treated at the mean-field level can screen out the confining potential near the center of the trap, leading to an effectively flat potential. Additionally in a spinor system one can have pairing terms which are either small or zero near the center of the trap. We demonstrate that, with these combined ingredients, robust topological edge states can exist in the collective excitation spectrum of such a system when the bulk gap is larger than the characteristic energy scale $\gamma\sim M\omega^{2}x_{\rm TF}a$ in the Thomas–Fermi regime where $\omega$ is the trapping frequency, $M$ is the mass of the constituent bosons, $x_{\rm TF}$ is the Thomas–Fermi radius of the system, and $a$ is the lattice constant. Additionally, for the case of $\gamma$ much smaller than the band gap, we will derive an analytical expression for the edge state dispersion, which is also a central result of this work. We illustrate our findings on a honeycomb Kane–Mele lattice model, although our results hold for other lattice geometries. The presence of the spin degree of freedom, on the other hand, is essential for an accurate reconstruction of the energy band structure of a general lattice model with idealized boundary conditions.

We focus on $S=1$ spinor condensates on a lattice. We introduce boson spinor operators $\hat{\bf\Psi}_{j}=(\hat{\Psi}_{j,1},\hat{\Psi}_{j,0},\hat{\Psi}_{j,-1})^{T}$ and consider a spin-one generalization of Kane–Mele model:

[TABLE]

where we denote the $3\times 3$ spin-one matrices as ${\bf S}=(S_{x},S_{y},S_{z})$ and the tunneling towards the left (right) next nearest-neighbor site is determined by $i\lambda\nu_{jj^{\prime}}=+i(-i)\lambda$ Kane and Mele (2005). Although we expect the results of the current work to exist for other models, we concentrate on Eq. (1) to be concrete, and because we considered this model in a previous work Galilo et al. (2015). The model exhibits two bands for each spin component. The next-nearest-neighbor term vanishes for spin 0 components and the Hamiltonian reduces to the graphene model. The other two spin components each have a gap of $4\lambda$ in energy spectrum, closed only by chiral topological edge modes. In the following, we will write this model more compactly as $\hat{\mathcal{H}}_{\rm latt}=\sum_{ij}\hat{\bf\Psi}^{\dagger}_{i}\mathcal{H}_{\rm{latt}}^{ij}\hat{\bf\Psi}_{j}$ where the $\mathcal{H}_{\rm{latt}}^{ij}$ matrix can be directly determined from E1. (1).

Since we are mainly interested in the manifestation of edge modes, we will consider the strip geometry. We keep the lattice periodicity along ${\bf a}_{1}$ primitive lattice vector, while considering harmonic trapping confinement in the perpendicular direction given by:

[TABLE]

where $\omega$ is the confining frequency, $\hat{\rho}_{j}=\hat{\bf\Psi}_{j}^{\dagger}\hat{\bf\Psi}_{j}$ is the local number operator, and $x_{j}$ is the location of the $j$th site with respect to the center of the trap taken along the $x$ direction as depicted in Fig. 1.

We now include on-site interaction terms that do not break spin rotational symmetry:

[TABLE]

where $\hat{\bf S}_{j}=\hat{\bf\Psi}_{j}^{\dagger}{\bf S}\hat{\bf\Psi}_{j}$ is the local spin operator, and $U$ and $U_{s}$ are the density and spin interaction strengths. We also consider a quadratic Zeeman term $\hat{\mathcal{H}}_{Z}=q_{Z}\sum_{j}\hat{\bf\Psi}_{j}^{\dagger}S_{z}^{2}\hat{\bf\Psi}_{j}$ which is experimentally adjustable using external magnetic fields or microwave fields Gerbier et al. (2006). The linear Zeeman term is omitted due to the $S_{z}$ symmetry of our model. The total Hamiltonian now reads

[TABLE]

where $\mu$ is the chemical potential.

At low enough temperatures atoms occupy the ground state $\bar{\bf\Psi}_{j}$ which we treat in mean field theory. We will focus on the case when the ground state is a polar state where $\bar{\bf S}_{j}=\bar{\Psi}_{j}^{\dagger}{\bf S}\bar{\Psi}_{j}=0$ and $\bar{\bf\Psi}_{j}=(0,\sqrt{\bar{\rho}_{j}},0)^{T}$. To achieve this, we require a positive $q_{Z}$ for positive $U_{s}$, and $q_{Z}>2|U_{s}|$ for negative $U_{s}$ Stenger et al. (1998). We consider a slowly varying confining potential, much slower than the lattice characteristic length scale, so the ground state population extends over a considerable number of lattice sites. Provided the density and spin fluctuations, $\delta\hat{\rho}=\hat{\rho}-\bar{\rho}$ and $\delta\hat{\bf S}=\hat{\bf S}-\bar{\bf S}$, are small, the fluctuations about the ground state $\hat{\bm{\psi}}=\hat{\Psi}-\bar{\Psi}$ up to quadratic order are described by the Bogoliubov–de Gennes (BdG) Hamiltonian:

[TABLE]

where we have introduced the mean kinetic energy per particle in the condensate ${\bar{\mathcal{H}}_{\rm latt}}=\sum_{i,j}\bar{\bf\Psi}_{i}^{\dagger}\ {\mathcal{H}_{\rm latt}^{ij}}\bar{\bf\Psi}_{j}/\sum_{j^{\prime}}\bar{\rho}_{j^{\prime}}$ and the effective potential:

[TABLE]

Prior to evaluating the fluctuation modes of the BdG Hamiltonian, we first discuss the effect of screening of confining potential. The ground state density profile $\bar{\Psi}$ is determined by the time-independent Gross–Pitaevskii equation:

[TABLE]

At large distances from the center of the trap the population density $\bar{n}_{j}$ vanishes and the confining potential takes the leading role in $V_{\rm eff}$. On the other hand, for distances closer to the center of the trap, the ground state population is larger. Under the Thomas–Fermi approximation for the mean-field density, applicable when $\sqrt{M\omega^{2}a^{2}w}\ll U{\max_{j}\bar{\rho}_{j}}$, the effective potential $V_{\rm eff}$ indeed will vanish identically inside the Thomas–Fermi radius given by $x_{\rm TF}=\sqrt{2(\mu-{\bar{\mathcal{H}}_{\rm latt}})/(M\omega^{2})}$. Shown in Fig.1 is the effective potential which is computed numerically and demonstrates the statements above. Such a screening occurs in standard scalar BECs. In this sense, interactions contribute to the screening of the harmonic trap inside the Thomas–Fermi region. Deviations from perfect Thomas–Fermi screening can be seen as oscillations on the scale of the lattice spacing in the effective potential $V_{\rm eff}$. These are due to variations of the harmonic potential within the unit cell that are not perfectly screened (see Supplemental Material). Although the magnitude of these variations is small in comparison to the band and gap widths controlled by $\omega$ and $\lambda$ we retain them in our numerical computation.

We now turn to the discussion of the fluctuation spectrum of the BdG Hamiltonian Eq. (5), which is the central focus of this work. Because of $S_{z}$ symmetry, the spin-$\pm 1$ components decouple from the spin-0 ones in $\hat{\mathcal{H}}_{B}$. The density fluctuation term in Eq. (5) carries only spin-0 components $\delta\hat{\rho}_{j}^{2}=\bar{\rho}_{j}(\hat{\psi}_{j,0}^{\dagger}\hat{\psi}_{j,0}+\hat{\psi}_{j,0}\hat{\psi}_{j,0}+{\rm H.c.})$, while the spin fluctuation term couples the spin-$\pm 1$ components: $\delta\hat{\bf S}_{j}^{2}=\bar{\rho}_{j}(\hat{\psi}_{j,+1}^{\dagger}\hat{\psi}_{j,+1}+\hat{\psi}_{j,-1}^{\dagger}\hat{\psi}_{j,-1}+\hat{\psi}_{j,+1}\hat{\psi}_{j,-1}+{\hat{\psi}_{j,+1}^{\dagger}\hat{\psi}_{j,-1}^{\dagger}}+{\rm H.c.})$. In typical spinor condensates, the spin-spin interaction is much smaller than the density-density interaction. For instance, for 87Rb, $U_{S}$ will be 2 orders of magnitude smaller than $U$. For this case, the Hamiltonian describing the $\pm 1$ components will be, to a good approximation, translationally invariant within the Thomas–Fermi radius. An interface is then encountered at $x=x_{\rm TF}$, where the effective potential rises (Fig. 1) and the system loses its translational invariance. Similar approximations have previously been employed to simplify the computation of the spatial structure of parametrically amplified spin modes Scherer et al. (2010).

To better understand the effect of trap screening on the manifestation of topological edge modes, we now turn to an effective description of the edge state emergence around the screening boundaries outside the Thomas–Fermi radius. As the edge states are composed entirely of $\pm 1$ components, we will restrict our attention to these modes. Outside of the Thomas–Fermi radius, we may neglect the last two terms in Eq. (5) as the mean-field particle density vanishes in this region, so that pairing terms are absent. We assume the boundary region to be sufficiently small in comparison to the lattice constant induced by optical lattice potential and will discuss the effect of these restrictions later.

In the following, we will describe a fairly general method to compute the dispersion of edge modes in a soft potential, but will give quantities relevant to our particular model in the Supplemental Material. Note that outside the Thomas–Fermi radius, the effective potential reduces to the trapping potential, and interactions are unimportant since the mean-field density vanishes in this region. As a result the spin components decouple and without loss of generality we will focus on one spin component. Let $k$ denote the wave number along the ${\bf a}_{1}$ periodic direction and $n$ label the unit cell along the ${\bf a}_{2}$ lattice vector. The two sublattices of our model are spatially separated and consequently the corresponding atoms at the same $n$th unit cell experience a different magnitude of the confining potential. While the excitations related to one sublattice experience an effective potential ${\cal V}(n)$, the excitations related to the other sublattice experience a potential ${\cal V}(n+s)$, where $s$ accounts for the relative lattice separation. For the specific mode Eq.l (1) of this Letter, $s=1/3$ (see Supplemental Material). The energy spectrum of the system outside of the Thomas–Fermi radius is given by the eigenvalues $E_{k}$ of a difference equation of the form

[TABLE]

where $H_{k,n-n^{\prime}}$ and $V_{n}={\rm diag}({\cal V}(n),{\cal V}(n+s))$ are $2\times 2$ matrices corresponding to the sublattice degrees of freedom (the treatment can also be generalized to larger matrices). In Eq. (8), $H_{k,n-n^{\prime}}$ follows from ${\mathcal{H}_{\rm latt}^{ij}}$ while $V_{n}$ follows from the effective potential Eq. (16).

We consider a lattice point ${\bar{n}}$ outside of the Thomas–Fermi radius and expand about it to first order: $V_{n}=V_{\bar{n}}+V_{\bar{n}}^{\prime}(n-{\bar{n}})$. We will search for solutions localized about $\bar{n}$. Such an approximation is valid since the effective potential is nearly linear in this region (see Fig. 1). We next introduce ${\bf\Phi}_{k,\varphi}=\sum_{n}e^{-in\varphi}{\bf\Phi}_{k,n}$. This is $2\pi$ periodic in $\varphi$ which follows directly from its definition. We now work with the wave function using this momentum representation for the motion in the $x$ direction. Then Eq. (8) in terms of ${\bf\Phi}_{k,\varphi}$ reduces to

[TABLE]

where $H_{k,\varphi}=\sum_{n}e^{-i(n-n^{\prime})\varphi}H_{k,n-n^{\prime}}$ is the Fourier transform of $H_{k,n-n^{\prime}}$. The following analysis is simplified by performing a unitary transformation so that matrices corresponding to the potential in the above are made to be proportional to the identity matrix. Such a transformation is achieved by ${\cal U}_{\varphi}={\rm diag}(1,e^{is\varphi})$. Under this transformation, Eq. (9) becomes

[TABLE]

where $\tilde{{\bf\Phi}}_{k,\varphi}={\cal U}^{\dagger}_{{\varphi}}{{\bf\Phi}}_{k,\varphi}$, $\tilde{H}_{k,{\varphi}}={\cal U}_{{\varphi}}^{\dagger}H_{k,{\varphi}}{\cal U}_{{\varphi}}$, $\tilde{V}_{\bar{n}}={\cal V}(\bar{n})\mathds{1}$, and $\gamma_{\bar{n}}={\cal V}^{\prime}({\bar{n}})$. The solution to this is readily found to be

[TABLE]

Here ${\cal P}$ is the path ordering operator for $\varphi$ and $\Delta\varphi=\varphi-\varphi_{0}$ where $\varphi_{0}$ is arbitrary. Note that the eigenvalues of $\tilde{H}_{k,{\varphi}}$ give the bulk band dispersions in the noninteracting limit.

For the case when $\gamma_{\bar{n}}$ is significantly smaller than the eigenvalue spacing of $\tilde{H}_{k,\varphi}$, i.e., much smaller than the bulk band gap, the adiabatic approximation can be used. To do so, we introduce the “instantaneous” eigenbasis: $\tilde{H}_{k,\varphi}\tilde{{\bm{\phi}}}_{k,\varphi}^{(\nu)}=\varepsilon_{k,\varphi}^{(\nu)}\tilde{{\bm{\phi}}}_{k,\varphi}^{(\nu)}$ where $\nu$ labels the eigenvectors or eigenenergies and let ${\bm{\phi}}_{k,\varphi}^{(\nu)}={\cal U}_{\varphi}\tilde{{\bm{\phi}}}_{k,\varphi}^{(\nu)}$. Here, $\varepsilon_{k,\varphi}^{(\nu)}$ are the bulk eigenstates of the noninteracting infinite system. Invoking the adiabatic approximation, and transforming back to the original variables, one finds

[TABLE]

where $A_{k,\varphi}^{(\nu)}=i{\bm{\phi}}_{k,\varphi}^{(\nu)\dagger}\partial_{\varphi}{\bm{\phi}}_{k,\varphi}^{(\nu)}$ is the Berry connection and $Z_{k,\varphi}^{(\nu)}={\bm{\phi}}_{k,\varphi}^{(\nu)\dagger}\sigma_{z}{\bm{\phi}}_{k,\varphi}^{(\nu)}$ ($\sigma_{z}$ is a Pauli matrix). Such derivations are standard for periodically driven quantum systems in the adiabatic limit (here $\varphi$ is analogous to time) Wilczek and Shapere (1989).

Requiring ${\bf\Phi}_{k,\varphi}$ to be $2\pi$ periodic in $\varphi$ enables one to determine the energies $E_{k}$ entering Eq. (12). Labeling these as $E_{k}^{(\bar{n}\nu)}$, one finds

[TABLE]

Different values of ${\bar{n}}$ correspond to edge states localized at different places along along the $x$ direction. This expression is remarkable in that it expresses the edge state dispersion purely in terms of the external potential and basic quantities of the bulk system. The (integer) Chern numbers of the bulk system are given by

[TABLE]

where $\Omega_{k,\varphi}^{(\nu)}=i\left(\partial_{k}{\bm{\phi}}_{k,\varphi}^{(\nu)\dagger}\partial_{\varphi}{\bm{\phi}}_{k,\varphi}^{(\nu)}-\partial_{\varphi}{\bm{\phi}}_{k,\varphi}^{(\nu)\dagger}\partial_{k}{\bm{\phi}}_{k,\varphi}^{(\nu)}\right)$ is the Berry curvature. Using this and Eq. (13), one sees that the edge state band energies change by $\gamma_{\bar{n}}{\cal C}_{\nu}$ as the one-dimensional Brillouin zone is traversed: $\Delta E_{k}^{(\bar{n}\nu)}=E_{k=2\pi}^{(\bar{n}\nu)}-E_{k=0}^{(\bar{n}\nu)}=\gamma_{\bar{n}}{\cal C}_{\nu}$.

A result similar to Eq. (13) was previously uncovered in Ref. Lee and Park (2015), though the methods used there are rather different from those used presently. The work Lee and Park (2015) considered a topological system under a strictly linear external field, the so-called Wannier–Stark ladder. A semiclassical theory of the system was developed using the methods of Ref. Sundaram and Niu (1999). The theory was quantized by the Bohr–Sommerfeld condition to deduce the eigenenergies. Our expression Eq. (13) reduces to the result of Ref. Lee and Park (2015) in the limit when ${\cal V}_{\bar{n}}$ is a strictly linear potential and the external field couples identically to sites with the same $n$ (i.e., $s=0$). We also point out that the above derivation can be extended to lattice systems having more than two sublattices. For this case, the term containing $Z_{k,\varphi}^{(\nu)}$ in Eq. (13) will be modified, reflecting the more complex lattice geometry.

We support our results by performing numerical computations. We evolve the Gross–Pitaevskii equation [Eq. (7)] in imaginary time to find the ground state in the Thomas–Fermi regime. We then compute the collective spin excitations. Comparison with the analytical result Eq. (13) shows excellent agreement (see Fig. 2). Deviations from Eq. (13) do occur for the edge states with energies close to the corresponding bulk band from which they emerge. For these energy levels the overlap between states Eq. (12) and bulk states cannot be ignored. The focus of this Letter, however, lies in describing the emergence of topological edge states with energies well within the gap and/or outside the bulk bands. Nevertheless, we point out that our method can be improved by requiring corresponding matching conditions at the screening radius and retaining higher order derivatives of the confining potential. We also find that the degeneracy of the edge states changes depending on its energy. As can be seen in Fig. 2, the chiral edge states can pass through a given energy level several times. However, the difference between the number of right and left movers, $N_{R}^{\nu}$ and $N_{L}^{\nu}$ correspondingly, at the same energy level is fixed by the topological structure of the bulk states: $N_{R}^{\nu}-N_{L}^{\nu}=C_{\nu}$. This is a direct indication of the topological nature of the bulk states and a consequence of the bulk-boundary correspondence Hasan and Kane (2010); Qi and Zhang (2011).

In summary, despite harmonic confinement being known for obscuring observation of topological edge states, we have shown that interactions can facilitate the emergence of these states in topological lattices populated by a spinor Bose condensate. In the Thomas–Fermi regime sharp boundaries emerge due to the screening of the harmonic trap inside the Thomas–Fermi radius. We have found localized states emerging outside the Thomas–Fermi radius. We have also shown that these states carry information about the lattice topology inside the screening radius. Our results are valid both for antiferromagnetic ($U_{s}>0$) and ferromagnetic ($U_{s}<0$) interactions. Though we have focused, for concreteness, on the spin-one Kane–Mele model, the same effect will be present in other spinor lattice systems, like a spinful Hofstadter model Aidelsburger et al. (2013). Our analysis gives new insight into the emergence of topological edge states at soft boundaries and opens new doors to the exploration of topological properties in optical lattice experiments with spin-1 bosons.

The authors would like to acknowledge Ari Turner for useful discussions. We are grateful for support from the Schrodinger Scholarship Scheme at Imperial College London (B. G.); the European Union’s Seventh Framework Programme for Research, Technological Development, and Demonstration under Grant No. PCIG-GA-2013-631002 and the Aspen Center for Physics under Grant No. PHYS-1066293 (R. B.)

Supplemental material

Appendix A Effects of harmonic potential in the absence of interactions

In order to emphasize the importance of the screening procedure described in the main text, we illustrate in Fig. 3 how the energy band structure is obscured by the harmonic potential when the interaction terms are switched off, i.e. $U=U_{s}=0$. The figure shows results of numerical direct diagonalization of the Kane–Mele Hamiltonian with the same lattice size and parameter choice as in Fig. 2 of the main text. The energy bands show a considerable deviation from the case of hard-wall boundary conditions which are widely accepted in solids. One important consequence of an external harmonic potential is that the region of extended bulk states is considerably reduced. Also, identification of a gap is rather subtle, since many states emerge and fill the gap region once the harmonic potential is switched on. The main text of this manuscript identifies how the bulk energy spectrum can be recovered by switching the interactions on.

Appendix B The origin of unscreened variations of effective potential

Here we explain the origin of the unscreened variations seen in the effective potential $V_{\rm eff}$ within the Thomas–Fermi radius as seen in Fig. 2 of the main text. Let us recall the time-independent Gross–Pitaevskii equation:

[TABLE]

and the effective potential which we defined as:

[TABLE]

where $V_{\omega,j}=M\omega^{2}x_{j}^{2}/2$ is the harmonic confining potential. In the above expression we have also introduced the mean kinetic energy per particle in the condensate ${\bar{\mathcal{H}}_{\rm latt}}=\sum_{i,j}\bar{\bf\Psi}_{i}^{\dagger}\ {\mathcal{H}_{\rm latt}^{ij}}\bar{\bf\Psi}_{j}/\sum_{j^{\prime}}\bar{\rho}_{j^{\prime}}$. Notice that inside the Thomas–Fermi radius the particle density does not vanish, i.e. $\bar{\Psi}_{j}\neq 0$. With this assumption we divide equation Eq. (15) by $\bar{\Psi}_{i}$ to obtain

[TABLE]

The above relation allows for small unscreened fluctuations of $V_{\rm eff}$ inside the Thomas–Fermi radius. This is due to the second term in the right hand-side of Eq. (17). This term is not constant and can fluctuate due to variation of the ground state function within the nearest neighbors. However, in the Thomas–Fermi regime the fluctuations due to the lattice Hamiltonian $\mathcal{H}_{\rm latt}^{ij}$ are small in comparison to the band and gap widths. This is reflected in the direct numerical computation shown in Fig. 1b of the main text. One can partly reproduce these oscillations. By introducing the particle density $\bar{\rho}_{j}=|\bar{\Psi}_{j}|^{2}$ one can write $\bar{\Psi}_{j}/\bar{\Psi}_{i}=\sqrt{\bar{\rho}_{j}/\bar{\rho}_{i}}\exp(i\Delta\theta_{i,j})$, where $\Delta\theta_{i,j}$ denotes the phase difference of the ground-state wave-function at sites $j$ and $i$. Since $V_{\rm eff}$ and $\bar{\mathcal{H}}_{\rm latt}$ are real the phase difference $\Delta\theta_{i,j}$ should not contribute to Eq. (17). Indeed, due to time-reversal symmetry the ground-state wave-function $\bar{\Psi}_{j}$ can be chosen to be real. Thus one has

[TABLE]

We next express $\bar{\rho}_{j}$ from the equation defining $V_{\rm eff}$ [Eq. (16)], and plug it into equation Eq. (18):

[TABLE]

Notice that so far we did not use any approximations, and the above relation is exact inside the Thomas–Fermi radius. Next, we regard Eq. (19) as a recursive relation for $V_{\rm eff}$. In what follows we compute the value of $V_{\rm eff}$ by plugging $V_{\rm eff}=0$ into the right-hand side of the relation Eq. (19). Then the obtained value we plug back into the right-hand side again, and repeat this procedure for a couple of times. The results after one iteration are shown in Fig.  4. There is excellent agreement with the exact value of $V_{\rm eff}$ in the middle of the trap. One can obtain a better approximation by including more terms in the expansion. However, this procedure gives poor convergence near the Thomas–Fermi radius. These results complement the screened value of the effective potential by taking into account the small variations within the nearest-neighbor hopping.

Appendix C Topological edge-states in presence of harmonic potential and interactions

Below we unveil the quantities used in (8) and thereafter in the main text in correspondence to the Kane–Mele model. Since in the region outside the Thomas–Fermi radius the system is effectively non-interacting the spin components decouple. Below we describe only the spin $m=1$ component. The other two components $m=-1$ and $m=0$ can be related to the $m=1$ case by setting $\lambda\rightarrow-\lambda$ and $\lambda\rightarrow 0$ with $q_{Z}\rightarrow 0$ respectively. The Hamiltonian matrix $H_{k,n-n^{\prime}}$ for the spin $m=1$ component is given by

[TABLE]

where $R_{k}$ and $M_{k}$ are $2\times 2$ matrices acting on the sublattice indices and are given by

[TABLE]

[TABLE]

Switching from $n-n^{\prime}$ to its conjugate $\varphi$ and introducing Pauli matrices $\sigma_{1}$, $\sigma_{2}$, and $\sigma_{3}$ reduces the Hamiltonian matrix to

[TABLE]

where $d_{1}=-w(1+\cos k+\cos\varphi)$, $d_{2}=-w(\sin k+\sin\varphi)$ and $d_{3}=-2\lambda\left(\sin\varphi-\sin(\varphi-k)-\sin k\right)$. The above Hamiltonian has two eigenvalues $\epsilon^{(\pm)}_{k,\varphi}=\pm\sqrt{d_{1}^{2}+d_{2}^{2}+d_{3}^{2}}=\pm d$ and the corresponding eigenvectors can be written as

[TABLE]

It is then straightforward to find an expression for the Berry connection $A_{k,\varphi}^{(\nu)}=i{\bm{\phi}}_{k,\varphi}^{(\nu)\dagger}\partial_{\varphi}{\bm{\phi}}_{k,\varphi}^{(\nu)}$ (in a particular gauge):

[TABLE]

This term causes the edge states to experience a shift in energy while wrapping them around the first Brillouin zone. There is also a term caused by the non-linearity of the potential:

[TABLE]

The trapping potential also lifts the energies of the edge mode excitations localized about lattice site $\bar{n}$. For reasons of computational convenience we have considered the trapping potential to be oriented in the direction perpendicular to the primitive lattice vector ${\bf a}_{1}$. Thus the excitations within the sublattice A experience a harmonic potential ${\cal V}(n)=M\omega^{2}n^{2}({\bf a}_{2}\cdot\hat{\bf x}-x_{c})^{2}/2+q_{Z}-\mu$ centered at position $x_{c}=({\bf a}_{2}\cdot\hat{\bf x})(N+1+s)/2$ and shifted by the chemical potential $\mu$ and quadratic Zeeman coupling $q_{Z}$, where $2N$ is the total number of sites in a unit cell strip highlighted in Fig. 1 of the main text. The sites of sublattice $B$, however, are displaced by a vector $\bm{\delta}$ and consequently experience a slightly different magnitude of the harmonic trap ${\cal V}(n+s)$. The value of the shift $s$ accounts for the displacement $\bm{\delta}$ and is related to the primitive lattice vectors ${\bf a}_{1}$ and ${\bf a}_{2}$. We find

[TABLE]

For the hexagonal lattice considered in this manuscript, $s=1/3$.