Self-Organized Bosonic Domain Walls

TL;DR

This paper explores how hardcore bosons on honeycomb lattice ribbons form self-organized domain walls that separate topologically distinct charge-density-wave regions, leading to superfluid transport along the domain walls, with implications for bosonic topological phenomena.

Contribution

It demonstrates the formation of self-organized bosonic domain walls with topological properties in honeycomb lattice ribbons, a novel phenomenon in bosonic systems.

Findings

Charge domain walls are energetically favorable.

Superfluid transport occurs along the domain walls.

The system exhibits topologically distinct charge-density-wave regions.

Abstract

Hardcore bosons on honeycomb lattice ribbons with zigzag edges are studied using exact numerical simulations. We map out the phase diagrams of ribbons with different widths, which contain superfluid and insulator phases at various fillings. We show that charge domain walls are energetically favorable, in sharp contrast to the more typical occupation of a set of sites on a single sublattice of the bipartite geometry at $\rho=\frac{1}{2}$ filling. This `self-organized domain wall' separates two charge-density-wave (CDW) regions with opposite Berry curvatures. Associated with the change of topological properties, superfluid transport occurs down the domain wall. Our results provide a concrete context to observe bosonic topological phenomena and can be simulated experimentally using bosonic cold atoms trapped in designed optical lattices.

Equations

$$H=-t\sum _{⟨i,j⟩}(b_i^{†}{b}_j^{}+\text{H.c.})+V\sum _{⟨i,j⟩}{n}_i{n}_j-\mu \sum _i{n}_i.$$

$$H=$$

$$\displaystyle{\cal H}({\bf k})=\left[\begin{array}[]{cc}\mu&f({\bf k})\\ f^{*}({\bf k})&3V-\mu\\ \end{array}\right],$$

$${\Omega }_{\lambda }(\mathbf{k})=\frac{\partial A_y(\mathbf{k})}{\partial k_x}-\frac{\partial A_x(\mathbf{k})}{\partial k_y},$$

$$\displaystyle H_{2}(k_{x})=\left[\begin{array}[]{cccc}0&-t\gamma_{k}&0&0\\ -t\gamma^{*}_{k}&0&-t&0\\ 0&-t&0&-t\gamma_{k}\\ 0&0&-t\gamma^{*}_{k}&0\\ \end{array}\right],$$

$$H=$$

$$\displaystyle H_{2}^{h}(k_{x})=\left[\begin{array}[]{cccc}-2V&-t\gamma_{k}&0&0\\ -t\gamma^{*}_{k}&-3V&-t&0\\ 0&-t&-3V&-t\gamma_{k}\\ 0&0&-t\gamma^{*}_{k}&-2V\\ \end{array}\right].$$

$$\displaystyle H_{3}(k_{x})=\left[\begin{array}[]{ccc}h_{1}(k_{x})&h_{2}&0\\ h^{\dagger}_{2}&h_{1}(k_{x})&h_{2}\\ 0&h^{\dagger}_{2}&h_{1}(k_{x})\\ \end{array}\right],$$

$$\displaystyle h_{1}(k_{x})=\left(\begin{array}[]{cc}0&-t\gamma_{k}\\ -t\gamma^{*}_{k}&0\\ \end{array}\right),h_{2}=\left(\begin{array}[]{cc}0&0\\ -t&0\\ \end{array}\right).$$

$$\displaystyle H_{4}(k_{x})=\left[\begin{array}[]{cccc}h_{1}(k_{x})&h_{2}&0&0\\ h^{\dagger}_{2}&h_{1}(k_{x})&h_{2}&0\\ 0&h^{\dagger}_{2}&h_{1}(k_{x})&h_{2}\\ 0&0&h^{\dagger}_{2}&h_{1}(k_{x})\\ \end{array}\right].$$

$$H=$$

$$S_{A,i}^{+}=(\sqrt{2S-a_{i,A}^{†}{a}_{i,A}^{}})a_{i,A}^{},$$

$$S_{A,i}^{-}=a_{i,A}^{†}(\sqrt{2S-a_{i,A}^{†}{a}_{i,A}^{}})$$

$$S_{A,i}^z=S-a_{i,A}^{†}{a}_{i,A}^{}.$$

$$S_{B,i}^{+}=a_{i,B}^{†}(\sqrt{2S-a_{i,B}^{†}{a}_{i,B}^{}}),$$

$$S_{B,i}^{-}=(\sqrt{2S-a_{i,B}^{†}{a}_{i,B}^{}})a_{i,B}^{}$$

$$S_{B,i}^z=a_{i,B}^{†}{a}_{i,B}^{}-S.$$

$$H=$$

$$\displaystyle{\cal H}({\bf k})=\left[\begin{array}[]{cc}\mu&f({\bf k})\\ f^{*}({\bf k})&3V-\mu\\ \end{array}\right]$$

$$\displaystyle\left[\begin{array}[]{cc}\mu&f({\bf k})\\ -f^{*}({\bf k})&-3V+\mu\\ \end{array}\right].$$

$$\displaystyle{\cal U}_{\bf k}=\left[\begin{array}[]{cc}\cosh\theta_{\bf k}e^{i\phi_{\bf k}}&-\sinh\theta_{\bf k}\\ -\sinh\theta_{\bf k}&\cosh\theta_{\bf k}e^{-i\phi_{\bf k}}\\ \end{array}\right],$$

$${\Omega }_{\lambda }(\mathbf{k})=\frac{\partial A_y(\mathbf{k})}{\partial k_x}-\frac{\partial A_x(\mathbf{k})}{\partial k_y},$$

Full text

††preprint: APS/123-QED

Introduction.- One of the most interesting properties of condensed matter systems is their condensation into ordered low temperature phases, breaking an underlying symmetry of the Hamiltonian. Such phases typically minimize the free energy $F$; coexistence of the distinct ordered patterns involves a domain wall, increasing $F$. Nevertheless, domain walls often exist in practice in experiments (or in simulations) as a consequence of long annealing times. This is especially the case in the presence of disorder which can pin their motion.

In addition to being manifest as meta-stable states, domain walls can also arise in other ways. An important example is provided by doping away from the commensurate antiferromagnetic (AF) filling of the cuprate superconductorsTranquada (2013), or the Hubbard and $t$-$J$ models that describe themZaanen and Gunnarsson (1989); White and Scalapino (2003); Hellberg and Manousakis (1999). Dopants do not spread uniformly, but instead form “charge stripes”. Across these stripes there is a ‘$\pi$-phase shift’ of the AF orderFradkin et al. (2015). The up-spin occupied sublattice interchanges across the stripe, realizing a domain wall.

In model Hamiltonian studies on ‘ladder’ geometries using the density matrix renormalization group, the charge patterns are found to be ‘vertical stripes’, i.e. the doped holes lie parallel to the short direction of the clusterWhite and Scalapino (2003). These charge patterns are fundamentally connected not only to magnetism, e.g. the $\pi$ phase shift, but also to charge density wave and $d$-wave pairing order. Studies of stripe physics and the associated domain walls remain of great interestZheng et al. (2017); Huang et al. (2017), with the possible coexistence of Luther-Emery liquid states in which the spin excitations are gapped, and quasi-long range superconducting correlations being a key issue Jiang and Devereaux (2018).

In this Letter we study bosonic particles on honeycomb ribbons. We discuss four novel features of this geometry. First, we argue that charge domain walls are energetically favorable compared to occupation of a set of sites on a single sublattice of the bipartite geometry, even at half-filling. This is a rather unique feature compared to situations in which domain walls are excitations rather than the ground state. Second, the low density sites of the domain wall are arranged ‘horizontally’ (parallel to the long axis), rather than vertically. These ‘self-organized domain walls’ open the possibility of superfluid transport down the chain. Third, associated with this physics is a non-trivial Berry curvature, which changes sign across the domain wall. Finally, the system realizes an exotic one-dimensional supersolid.

The model and method.- We consider spinless hardcore bosons on zigzag ribbons of a honeycomb lattice, described by the extended Bose-Hubbard model Bloch et al. (2008); Cazalilla et al. (2011)

[TABLE]

$b^{\phantom{\dagger}}_{i}(b_{i}^{\dagger}$) are hardcore boson annihilation(creation) operators, $n_{i}=b_{i}^{\dagger}b^{\phantom{\dagger}}_{i}$ is the number operator. Hardcore bosons obey commutation relations $[b^{\phantom{\dagger}}_{i},b_{j}^{\dagger}]=0$ for sites $i\neq j$ and on-site anticommutation relations $\{b^{\phantom{\dagger}}_{i},b_{i}^{\dagger}\}=1$. The first term in Eq. (1) describes nearest-neighbor (NN) hopping, with amplitude $t$ taken as the unit of energy $(t=1)$. The second term in Eq. (1) is the NN interaction $V$. Finally, $\mu$ denotes the chemical potential, which controls the number of bosons in the system. The model in Eq. (1) has a $U(1)$ symmetry, and is invariant under the transformation $b_{j}\rightarrow e^{{\rm i}\theta}b_{j}$ where $\theta$ is a real-valued phase. This symmetry is spontaneously broken in a superfluid phase. The Hamiltonian is invariant under the inversion transformation center of honeycomb lattice ribbons with zigzag edges.

The band structure of the honeycomb lattice consists of two inequivalent Dirac points, which are characterized by $\pm\pi$ Berry phasesCastro Neto et al. (2009). As a result of the nontrivial topological property, localized flat bands connecting the two Dirac points appear on the zigzag edgesFujita et al. (1996); Nakada et al. (1996); Son et al. (2006); Yao et al. (2009). Figure 1(b) shows the low-energy bands. As the widths of the ribbons are reduced, the lengths of the flat bands are shortened. The band bottom corresponds to the chemical potential at which the hardcore bosons begin to fill into the system, which determines the phase boundary for $\rho=0$.

We employ the stochastic series expansion (SSE) quantum Monte Carlo (QMC) method Bauer et al. (2011) with directed loop updates to study Eq.(1). SSE expands the partition function in a power series and the trace is written as a sum of diagonal matrix elements. The directed loop updates and the fact that the discrete configuration space can be sampled without floating point operations make the approach very efficient Syljuåsen and Sandvik (2002); Alet et al. (2005); Pollet et al. (2004). Our simulations are on finite lattices with the total number of sites $N=2\times W\times L$ with $W$ the width and $L$ the length of a ribbon[see Fig.1(a)]. The temperature is set to be low enough to obtain the ground-state properties. We also use the exact diagonalization (ED) method, which is numerically exact, but has strong size limitations.

The phase diagram.- Phase diagrams of ribbons with $W=2$ are shown in Fig. 2(a) (see others in Ref.[28]). In the atomic limit ($t/V=0$), the density abruptly jumps from empty ($\rho=0$) to half-filled ($\rho=\frac{1}{2}$) at $\mu/V=0$, with bosons placed so they never occupy adjacent sites. At $\rho=\frac{1}{2}$, no further bosons can be added without being neighbors, costing energy $\propto V$: there is a jump in $\mu$. If the half-filled bosons are placed so that they occupy only a single sublattice, the empty sites of one of the boundaries are special: they interact with only two neighboring occupied sites. Thus the $\rho=\frac{1}{2}$ CDW insulator terminates at $\mu/V=2$cdw . Once these special sites are completely occupied, the increase in density pauses again until $\mu/V=3$, at which point bosons are added to the remaining empty sites with three occupied neighbors, completely filling the lattice. This atomic limit picture explains the positions of the insulating lobes bases in Fig. 2(a).

The $\rho=\frac{1}{2}$ insulator for $0<\mu/V<2$ has a $(W+1)$ fold degeneracy, including two single sublattice CDWs and $(W-1)$ configurations with a domain wall arranged along each row of vertical bonds. The key observation for the presence of self-organized domain walls of Fig. 2(a) is that a domain wall has a pair of ‘edges’ where empty sites have only two occupied neighbors. This greater multiplicity of special sites compared to a single sublattice leads to a lower energy at finite hopping $t/V\neq 0$: the second order energy decrease when a boson hops onto a special site is $-t^{2}/V$ compared to the $-t^{2}/2V$ for hopping within the CDW. Fig. 1(a) illustrates these different hopping processes, and the larger overall energy decrease, $-5t^{2}/2V$, of a site adjacent to a domain wall compared to $-3t^{2}/2V$ gained by a boson inside the CDW. Filling half of the special sites in the domain wall phase also explains the densities $\rho=\frac{1}{2}+\frac{1}{2W}$ of the $2<\mu/V<3$ insulating lobe of Fig. 2(a).

The atomic insulator phases initially persist at small $t/V$, but the range in chemical potential over which they are stable decreases. They completely disappear beyond a critical value of $t/V$. We note that there is an additional valence-bond insulator at $\rho=\frac{1}{4}$ for $W=2$ which has no atomic counterpartHou et al. (2007); Bergman (2013). At non-zero $t/V$, all these insulators are separated by incommensurate superfluid regions.

Quantum phases suggested by these strong coupling arguments can be precisely determined using the SSE by measuring the average density $\rho=\frac{1}{N}\sum_{i}\langle n_{i}\rangle$ and the superfluid density $\rho_{s}=\langle W^{2}\rangle/4\beta t$, where $W$ is the winding number counting the net number of times the paths of the particles have wound around the periodic cellPollock and Ceperley (1987); Rousseau (2014), and $\beta$ is the inverse temperature. Insulating behavior is characterized by $\rho_{s}=0$ and a plateau of $\rho$ representing the persistence of the atomic limit steps in the chemical potential to finite $t/V$. Conversely, the superfluid has nonzero $\rho_{s}$ and finite compressibility $\kappa=\partial\rho/\partial\mu$. These features are clearly seen in the SSE results of Fig. 2(b) for $W=2$. A collection of plots like Fig. 2(b) for different $t/V$ generates the phase diagram in Fig. 2(a).

We analyze the lifting of the atomic limit $(W+1)$-fold degeneracy of the $\rho=\frac{1}{2}$ insulator by the hopping non-perturbatively on small lattices using ED. As shown for $W=2$ in Fig. 3(a), as $t/V$ increases, $W+1=3$ distinct eigenenergy curves emerge from the degenerate $t/V=0$ limit. The site densities and the density-density correlations at Figs. 3(c) and 3(e) confirm the lowest is the unique domain wall phase. When $t$ is non-zero, the single-sublattice CDW states form linear combinations and in the resulting state superpositions, all sites have average filling around 0.5. However the domain wall state, because of inversion symmetry, has two inequivalent sets of sites, half of which have densities which are low, and half which are high. In the domain wall phase, both edge sites of the unit cell have high densities, while one of them is occupied and the other is not in the CDW phase. The density-density correlations between such sites have distinct values, which are large for the domain wall phase, and small for the CDW [see Fig. 3(e)]. Thus the presence of two well-separated $\rho_{i}(t/V)$ trajectories and large density-density correlation between opposite-edge sites are ‘smoking guns’ that the ground state has a domain wall, as already suggested by the strong coupling argument.

Figures 3(b,d,f) demonstrate the existence of a domain wall in the $W=3$ ground state. In Fig. 3(b), the $W+1=4$-fold atomic limit degeneracy is lifted by hopping. In the larger lattice results, the two lowest states are nearly degenerate, as the two upper ones SM . If one averages the $\rho_{i}$ of the two domain wall states, one finds two ‘occupied’ sites which have $\rho_{i}$ large in both states, two ‘empty’ sites which have $\rho_{i}$ small in both states, and two sites which exchange $\rho_{i}$ small and large. This leads to three $\rho_{i}$ trajectories: large, small, and approximately half-filled. Meanwhile, each of the CDW states has six inequivalent sites, three high and three low density, which all exchange between the two degenerate cases. When averaged, all sites would therefore have $\rho_{i}\sim 0.5$. The density-density correlation between opposite-edge sites in the domain wall phase is much larger than that in the CDW phase [See Fig. 3(d)]. Thus the densities and density-density correlations observed in Fig. 3(f) and 3(d) for $W=3$ offer compelling evidences that the ground state manifests a charge domain wall.

The Berry curvature and the domain-wall superfluid.- Having established the phase diagram and the existence of domain walls in the $\rho=1/2$ insulator, we now focus on behavior at the interface. The topological property of bulk CDW phase of hardcore bosons is first investigated. The Bose-Hubbard Hamiltonian in Eq.(1) is equivalent to a spin$-1/2$ $XXZ$ model through a mapping $S^{+}_{i}=b^{{\dagger}}_{i}$ and $S^{z}_{i}=n_{i}-\frac{1}{2}$. The associated Holstein-Primakoff transformation yields,

[TABLE]

where $a^{\dagger}_{i,\alpha},a^{\phantom{\dagger}}_{i,\alpha}$ ($\alpha=A,B$ denoting the sublattice) are the bosonic creation and annihilation operatorsbos . In momentum space, $H=\sum_{\bf k}\psi^{\dagger}_{\bf k}{\cal H}({\bf k})\psi^{\phantom{\dagger}}_{\bf k}$, where $\psi_{\bf k}=\{a^{\phantom{\dagger}}_{A,{\bf k}},a^{\dagger}_{B,-{\bf k}}\}^{T}$, and

[TABLE]

with $f({\bf k})=-t(1+e^{-{\rm i}{\bf k}\cdot{\bf a}_{1}}+e^{-{\rm i}{\bf k}\cdot{\bf a}_{2}})$. The magnon band structure, i.e., the excitation of the mapped spin$-1/2$ $XXZ$ model, has two branches: $E^{\pm}_{\bf k}=\pm(\mu-\frac{3V}{2})+\epsilon({\bf k})$, where $\epsilon({\bf k})=\sqrt{(\frac{3V}{2})^{2}-|f({\bf k})|^{2}}$ [see Fig.4(a)]. The Berry curvature associated with each magnon band is given by

[TABLE]

where $A_{i}=-{\rm i}\langle u_{\lambda,{\bf k}}|\frac{\partial}{\partial k_{i}}|u_{\lambda,{\bf k}}\rangle$ ($i=x,y$) is the Berry potential. $\lambda=\pm$ denotes the two magnon bandsOwerre (2016); Guo et al. (2016); Fukui et al. (2005). The Berry curvature is peaked at the Brillouin zone (BZ) corners, Fig. 4(a), and is antisymmetric with respect to the inversion center ${\bf k}=(0,0)$. The sum of the Berry curvature of each band in the BZ (the Chern number) vanishes identically. The Berry curvatures for the two $\rho=\frac{1}{2}$ CDW insulators differ by an overall sign, thus the Berry curvature changes the sign across the domain wall. Due to bulk-boundary correspondence, there should appear in-gap domain-wall statesSemenoff et al. (2008); Cheng et al. (2018). As shown in Fig.4(b), there are two such branches associated with the domain wall. One of them is at the bottom of the magnon excitation spectrum, and it corresponds to the superfluid above the $\rho=\frac{1}{2}$ CDW insultor, which is localized near the domain wall.

For the density $\rho=\frac{1}{2}$, the domain wall is thus formed by a quasi-1D region with negligible occupancy, confined by robust CDW phases. As the chemical potential is further increased, bosons are preferentially added to the empty sites forming the domain wall, since they only interact with two occupied neighbors in contrast with three neighbors of an empty CDW site. A single extra boson can hop freely along the chain, lowering its kinetic energy without changing the interacting energy. Additional added particles behave effectively as interacting bosons in 1D, which condense to superfluid transport down the domain wall. The values of the superfluid density, Fig. 2(b), follow a dome shape, and are maximal at density $\sim\frac{1}{2}+\frac{1}{4W}$ when half of such empty sites are occupied. After the domain wall is full, the superfluid vanishes, and the system becomes a $\rho=\frac{1}{2}+\frac{1}{2W}$ insulator. Figure 5 shows the density for $W=2,4$. The presence of two sets of well-separated local traces is consistent with the $\rho=\frac{1}{2}$ insulator being a domain-wall phase. (See discussion of Fig. 3(e).)

While the local density of the sites on the domain wall with small occupation increases, that of the high-occupancy sites first decreases even as $\mu$ grows. This anomalous behavior is a signature of the flow of bosons onto the domain wall and the appearance of a superfluid localized near the domain wall. It is noteworthy that the domain wall superfluid coexists with diagonal (density) order: the zigzag honeycomb nanoribbon realizes an exotic one-dimensional supersolidBatrouni et al. (2006). The emergence of this superfluid is a manifestation of change in topological properties when crossing the domain wall.

To verify the localization of the superfluid near the domain wall, we show the single-particle correlator $\langle b^{\dagger}_{0}b^{\phantom{\dagger}}_{r}\rangle$ in Fig. 6. The correlator along the zigzag chain on the domain wall is slower than a power-law decay with distance, which is characteristic of a gapless quasi-1D superfluid. In contrast, the excitation is gapped for the $\rho=\frac{1}{2}$ domain-wall insulator, and the correlator decays exponentially. As one moves away from the domain wall, the correlator becomes increasingly short-ranged, and $\rho_{s}$ decreases. For wide ribbons, the superfluid density decays exponentially with the distance away from the domain wallSM .

Conclusions.- We studied hardcore bosons on zigzag edge honeycomb lattice ribbons using exact simulations. The phase diagram contains superfluid and insulating phases and, remarkably, at $\rho=\frac{1}{2}$ filling the ground state contains a charge domain wall rather than occupation of a single sublattice. This ‘self-organized domain wall’ separates CDW regions with opposite Berry curvature, and supports superfluid transport in coexistence with diagonal (density) order. Our results demonstrate that honeycomb ribbons provide a concrete geometry for the observation of bosonic topological phenomenon.

This physics can be explored experimentally. Cold atoms in optical lattices provide a well-established means to emulate the Bose-Hubbard modelJaksch et al. (1998), large values of $U$ which achieve the hard-core limit can be attained, and the honeycomb geometry has been generated Polini et al. (2013); Zhu et al. (2007); Soltan-Panahi et al. (2011). The use of a synthetic dimension may also provide a simple realization of honeycomb ribbonsSuszalski and Zakrzewski (2016); Ozawa and Price (2019). New tools based on quantum gas microscopes allow observation of the density profile at the level of individual atomsGross and Bloch (2017); Bloch et al. (2012); Bakr et al. (2009); Sherson et al. (2010), and hence direct comparison with our real-space measurements. The Berry curvature can also be obtained via interferometric techniquesDuca et al. (2015); Fläschner et al. (2016).

Acknowledgments.- The authors thank G. Batrouni, G. Chen and F. Nori for helpful discussions. H.G. acknowledges support from the NSFC grant No. 11774019. X.Z. and S.F. are supported by the National Key Research and Development Program of China under Grant No. 2016YFA0300304, and NSFC under Grant Nos. 11574032 and 11734002. R.M. acknowledges support from NSAF-U1530401, and from NSFC Grant No. 11674021 and No. 11851110757. R.T.S. was supported by DOE grant DE-SC0014671.

Appendix A 1. Band structure

The band structure of honeycomb lattice ribbons can be determined analytically. One fermion has exactly the same energy as one hardcore boson due to the absence of exchange statistics. For $W=2$, the Hamiltonian in the momentum space is,

[TABLE]

with $\gamma_{k}=1+e^{ik_{x}}$. The energy spectrum contains four branches, $E_{i}=(\pm 1\pm\sqrt{9+8\cos k_{x}})t/2$ ($i=1,2,3,4.$). The band bottom is located at $k_{x}=0$, and the corresponding eigenvalue is $-(1+\sqrt{17})t/2$ (See Fig. S7). Thus the lower boundary of the phase diagram, where the density first begins to become nonzero, is a straight line $\mu/V=-\frac{1+\sqrt{17}}{2}t/V$.

Under a particle-hole transformation $b^{\dagger}_{i}\rightarrow h_{i}$, the Hamiltonian, Eq. (1) of the main text, becomes,

[TABLE]

where $n^{h}_{i}=h_{i}^{\dagger}h_{i}$ is the hole number operator, and $E_{0}=\frac{3}{2}VN-LV-\mu N$ ($N=2L\times W$ the total number of sites) is a constant. In momentum space, the hole Hamiltonian is

[TABLE]

The energy spectrum $E_{1,2}=-\frac{t}{2}-\frac{5V}{2}\pm\frac{1}{2}P_{k,+};E_{3,4}=\frac{t}{2}-\frac{5V}{2}\pm\frac{1}{2}P_{k,-}$ with $P_{k,\pm}=\sqrt{9t^{2}+8t^{2}\cos k_{x}+V^{2}\pm 2tV}$. The band bottom determines the upper boundary of the phase diagram, which is described by the curve $\frac{\mu}{V}=\frac{1}{2}\frac{t}{V}+\frac{1}{2}\sqrt{17\frac{t^{2}}{V^{2}}+2\frac{t}{V}+1}+\frac{5}{2}$.

For the $W=3$ case, the Hamiltonian is,

[TABLE]

where

[TABLE]

The lower boundary is described by $\mu/V=-2.76\,t/V$.

For $W=4$, the Hamiltonian is

[TABLE]

From the eigenvalue of the band bottom, we have $\mu/V=-2.851\,t/V$ describing the phase boundary for $\rho=0$. The phase boundary for $\rho=1$ can be obtained numerically in the hole representation.

Appendix B 2. The phase diagrams for several other widths

The phase diagrams of the ribbons with several other widths are shown in Fig.S8, which contain superfluid regions and insulator phases at specific fillings. Compared to that of the periodic honeycomb lattice, the regions of the $\rho=\frac{1}{2}$ insulators are greatly reduced, and there appear new insulating regions at the filling $\rho=\frac{1}{2}+\frac{1}{2W}$. There is an additional valence-bond insulators at $\rho=\frac{3}{4}$ for $W=4$ which has no atomic counterpart.

Appendix C 3. Strong Coupling

The Exact Diagonalization and Density Matrix Renormalization methods can be used to obtain further details of the evolution of the low energy spectrum and site densities away from the atomi ($t=0$) limit. The $W+1$-fold degeneracy of the atomic limit $\rho=\frac{1}{2}$ insulator is lifted by the hopping. For the $W=2$ case, the ground state is the unique domain wall phase. The first- and second- excited states are linear combinations of the single-sublattice CDW states to preserve the symmetries of the Hamiltonian. As shown in Fig. S9, the $N=16$ site ED eigenenergies of the two mixed single-sublattice CDW states lie above the domain wall state and are slightly split by finite size effects. This finite size splitting is significantly reduced in the $N=32$ site ED and the $N=48$ site DMRG calculations. We plot the site densities of the mixed CDW states in Fig. S10, which are the first- and second-excited states. All sites have average filling around $0.5$. As the value of $t/V$ increases, the density difference between the two inequivalent sites becomes larger due to stronger quantum fluctuations.

For the $W=3$ case, the lowest four eigenenergies are grouped into two sets, each of which contains two mixed domain-wall or CDW states. Due to the finite-size effect, the eigenenergies are split. As the lattice size is increased (see Fig. S9), the splitting of the lower (upper) two eigenenergies is significantly reduced, suggesting they are degenerate in the thermodynamic limit. The site densities of the lower (upper) two states are plotted in Fig. S11 (Fig. S12), which are consistent with those of the domain-wall (CDW) phases. Thus the domain-wall phase is the ground state of the $W=3$ ribbon.

Appendix D 4. The Berry curvature

The model in Eq. (1) is equivalent to a spin$-1/2$ $XXZ$ model through a mapping $S^{+}_{i}=b^{{\dagger}}_{i}$ and $S^{z}_{i}=n_{i}-\frac{1}{2}$,

[TABLE]

Using the Holstein-Primakoff transformation, the spin operators are expressed in terms of bosonic creation and annihilation operators. The honeycomb lattice is bipartite. The transformation on sublattice $A$ is defined as

[TABLE]

Conversely, on sublattice $B$, the spin is in the opposite direction for antiferromagnetic order. Thus the spin operators are defined as

[TABLE]

Expanding the square root in Eq.(34) in powers of $1/S$, the zeroth order terms are kept in the linear spin-wave theory. Then the bosonic tight binding Hamiltonian becomes

[TABLE]

We ignore the four-operator terms and a constant. Under a Fourier transformation, $H=\sum_{\bf k}\psi^{\dagger}_{\bf k}{\cal H}({\bf k})\psi^{\phantom{\dagger}}_{\bf k}$, where $\psi_{\bf k}=\{a^{\phantom{\dagger}}_{A,{\bf k}},a^{\dagger}_{B,-{\bf k}}\}^{T}$ is the basis, and

[TABLE]

with $f({\bf k})=-t(1+e^{-{\rm i}{\bf k}\cdot{\bf a}_{1}}+e^{-{\rm i}{\bf k}\cdot{\bf a}_{2}})$ [${\bf a}_{1}=(\sqrt{3},0),{\bf a}_{2}=(\sqrt{3}/2,3/2)$ are the primitive vectors]. To diagonalize the above Hamiltonian, we consider the following non-Hermitian matrix,

[TABLE]

The eigenvalues are given by $E^{\pm}_{\bf k}=\mu-\frac{3V}{2}\pm\epsilon({\bf k})$ with $\epsilon({\bf k})=\sqrt{(\frac{3V}{2})^{2}-|f({\bf k})|^{2}}$. The eigenvector matrix is

[TABLE]

where $\sinh 2\theta_{\bf k}=\frac{|f({\bf k})|}{\epsilon({\bf k})}$, $\tan\phi_{\bf k}=\frac{Imf({\bf k})}{Ref({\bf k})}$. The first (second) column is the eigenvector $u_{+,{\bf k}}$ ($u_{-,{\bf k}}$) corresponding to $E^{+}_{\bf k}$ ($E^{-}_{\bf k}$). The Hamiltonian is thus diagonalized by the transformation: ${\cal U}^{\dagger}_{\bf k}{\cal H}({\bf k}){\cal U}_{\bf k}=\textrm{diag}(E^{+}_{\bf k},-E^{-}_{\bf k})$. In Fig.S13, we plot the magnon band structure at $\mu=\frac{3V}{2}=6t$, which consist of two gapped branches.

The above magnon bands of the mapped spin$-1/2$ XXZ model correspond to the excitation spectrum above the CDW insulator of the Bose-Hubbard model, which is the superfluid. When the spectrum becomes gapless, i.e., $E^{\pm}_{\bf k}=0$, superfluid begins to appear in the system. Thus $E^{\pm}_{\bf k}=0$ determine the phase boundary between the CDW and superfluid phases, which is $\mu=\frac{3V}{2}-\sqrt{(\frac{3V}{2})^{2}-(3t)^{2}}$ (the maximum value of $|f({\bf k}|)$ is $3t$). We plot the phase boundary from the spin-wave approximation in Fig.S14. It is qualitatively consistent with the exact phase diagram from the QMC method except that the range in the chemical potential is slightly increased. We also determine the phase boundary of the $\rho=\frac{1}{2}$ domain-wall phase on a $W=12$ zigzag ribbon, which greatly shrinks from that of a periodic system.

The Berry curvature associated with each magnon band is given by

[TABLE]

where $A_{i}=-{\rm i}\langle u_{\lambda,{\bf k}}|\frac{\partial}{\partial k_{i}}|u_{\lambda,{\bf k}}\rangle$ ($i=x,y$) is the Berry potential, and $\lambda=\pm$ denotes the two magnon bands.

Appendix E 5. The domain-wall superfluid

To demonstrate the domain-wall superfluid clearly, we perform QMC simulations on a $W=12$ and $L=24$ ribbon. Since there are more approximately degenerate $\rho=\frac{1}{2}$ phases on wider ribbons, a small pinning field is used to select a specific configurationAssaad and Herbut (2013). The pinning field is a staggered potential with the value $-|\Delta|$ ($|\Delta|$) on each occupied (unoccupied) site of the targeted domain-wall phase. In the following simulations, the strength of the pinning field is set to $|\Delta|=0.1$, and we focus on the phase with the domain wall in the middle.

Figure S15 shows the average density and superfluid density as a function of chemical potential at $V/t=4$. As the chemical potential increases, bosons are continuously added to the $\rho=\frac{1}{2}$ domain-wall phase until the domain wall is full, and the ribbon becomes a $\rho=\frac{1}{2}+\frac{1}{2W}$ insulator. Between the two insulators, the superfluid density is nonzero, implying the system is a superfluid.

Next we demonstrate the superfluid is localized near the domain wall, forming a one-dimensional superfluid channel. Figure S16 shows the local densities as functions of chemical potential. For the $\rho=\frac{1}{2}$ insulator, the profile of the local density indicates it is a domain-wall phase with the domain wall in the middle of the ribbon. From about $\mu=5$, additional bosons are added to the system. The value on the low-occupation sites of the domain wall increases significantly. In the $\rho=\frac{1}{2}+\frac{1}{2W}$ insulator, the site densities approach $\rho_{i}\sim 0.5$, which corresponds to the case the domain wall is full. The local densities are only slightly affected for sites more than $2a$ ($a$ the lattice constant) away from the domain wall. Thus the added bosons mainly reside near the domain wall. It is also noted that the local densities on the highly-occupied sites near the domain wall first decrease. This implies the bosons begin to flow between high- and low-occupation sites, which is consistent with the appearance of superfluid. Since the decrease is most significant on the domain wall, the gapless superfluid is mainly around it.

We also calculate the single-particle correlator. As the distance $r_{\bot}$ away from the domain wall increases, the decay of the correlator becomes rapid, suggesting the superfluid density decreases. For the zigzag chain on the domain wall, the decay of the correlator is slower than power law. However the decay becomes exponential from the fourth zigzag chain, where the superfluid begins to vanish. Moreover the curves remain almost unchanged for the far zigzag chains, manifesting the uniform bulk CDW order there. Inset shows $\langle b^{\dagger}_{0}b^{\phantom{\dagger}}_{r}\rangle$ as a function of $r_{\bot}$ (the distance away from the domain wall) at fixed $r=6$. The data are well fit using an exponential decay. Thus the superfluid is localized near the domain wall, and decays exponentially into the bulk.

Appendix F 6. The domain wall with turns

It is interesting to see whether it is possible to have domain walls with turns. The domain wall organized by the NN interactions should be straight and along the ribbon. The open zigzag edges of the ribbons are created by removing a line of vertical bonds. Then both sites connected by such a bond can be occupied. To be at half filling, we have a line of vertical bonds empty in the bulk, which form the ”self-organized” domain wall. Once the domain wall has turns, non-vertical bonds crossed should be empty. To maintain the number of bosons at half filling, there must appear adjacent occupied sites connected by periodic boundary condition, which induce interacting energy. As shown in Fig.S18(a), a ”Z” shaped domain wall crosses one non-vertical bond, and there appear a pair of adjacent bosons generating an interaction $V$. A domain wall along a direction other than the ribbon is shown in Fig.S18(b). It crosses three non-vertical bonds, and there appear three pairs of adjacent bosons with an interacting energy $3V$. Generally when a domain wall crosses the non-vertical bonds for $n$ times, there appear $n$ pairs of adjacent bosons, and the interacting energy increases by $nV$. So for large interactions, a domain wall with turns is not energetically favored.

Although the domain walls with turns can not be generated by the NN interactions, they can be designed using staggered potentials. We create a ”Z” shaped domain-wall using staggered potentials with the strength $|\Delta|=4t$ on a $W=12$ and $L=24$ ribbon, which is periodic along the longer axis. Figure S19(a) plots the average density and the superfluid density as a function of $\mu$. The $\rho=\frac{1}{2}$ plateau with vanishing superfluid corresponds to the domain-wall insulator. As the chemical potential further increases, bosons are added to the $\rho=\frac{1}{2}$ insulator. The superfluid density has finite values, and the system becomes superfluid. We calculate the distribution of the added bosons, $\delta\rho_{i}=\rho_{i}(\mu_{2})-\rho_{i}(\mu_{2})$, with $\mu_{2},\mu_{1}$ marked in Fig.S19(a). As shown in Fig.S19(b), the added bosons mainly distribute near the $Z$ shaped domain wall, implying the superfluid above the $\rho=\frac{1}{2}$ transporting down such a domain wall. With this method, complex shaped domain wall can be designed, realizing circuits of superfluid.