Boundary edge networks induced by bulk topology

TL;DR

This paper develops an effective edge network theory to analyze boundary topology in topological insulators, accurately predicting bound state properties and clarifying complex edge phenomena in higher-order and fullerene systems.

Contribution

It introduces a novel edge network approach that captures boundary topology and bound state characteristics in coupled edge states of topological insulators.

Findings

Edge networks predict energy and fractional charge of bound states.

The theory applies to second-order topological insulators and topological fullerenes.

It simplifies understanding of complex boundary phenomena.

Abstract

We introduce an effective edge network theory to characterize the boundary topology of coupled edge states generated from various types of topological insulators. Two examples studied are a two-dimensional second-order topological insulator and three-dimensional topological fullerenes, which involve multi-leg junctions. As a consequence of bulk-edge correspondence, these edge networks can faithfully predict properties such as the energy and fractional charge related to the bound states (edge solitons) in the aforementioned systems, including several aspects that were previously complicated or obscure.

Figures

Equations

$$H_{\text{edge}}=\sum _i\int d{x}_i{\Psi }^{†}(x_i)(-i{v}_i{\partial }_{x_i}{\sigma }_z+{\mathcal{M}}_i({\theta }_i))\Psi (x_i).$$

$${\tilde{H}}_{\text{edge}}=\sum _i(H_0^i+H_{\text{int}}^i).$$

$$H_{0,1}^i$$

$$H_{0,2}^i$$

$$H_{\text{int,2}}^i$$

$$H_{\text{int,4}}^i$$

$$H_{B,0}^i$$

$$H_{B,1}^i$$

$$\hat{Q}\equiv \int dx{j}_{\mu }(x)\propto \frac{{ϵ}^{01}}{2\pi }[\theta (x_2=+\infty )-\theta (x_1=-\infty )].$$

$$\phi =|{\theta }_{i+1}-{\theta }_i|/2,E=\text{sgn}({\theta }_{i+1}-{\theta }_i)\cos \phi ,N_s=-\frac{\phi }{\pi }.$$

$${\mathcal{M}}_1(0)=+{\sigma }_x,{\mathcal{M}}_2(\frac{\pi }{2})=+{\sigma }_y,$$

$${\mathcal{M}}_3(\pi )=-{\sigma }_x,{\mathcal{M}}_4(\frac{3\pi }{2})=-{\sigma }_y.$$

$$H$$

$$\mathcal{H}(\vec{k})=$$

$$+(m_z-2t_x\cos k_x-2t_y\cos k_y){\tau }_3{\sigma }_3$$

$$+\lambda {\tau }_1{\sigma }_1+\lambda {\tau }_2{\sigma }_2,$$

$$\left\{\begin{aligned} &\Psi_{1}(0)=\big{(}\frac{e^{-\varphi-i\theta_{1}}}{\sqrt{1+e^{-2\varphi}}}+r_{1}\frac{e^{\varphi-i\theta_{1}}}{\sqrt{1+e^{2\varphi}}}\big{)}\psi_{\alpha}({x_{1}=0})+\big{(}\frac{1}{\sqrt{1+e^{-2\varphi}}}+\frac{r_{1}}{\sqrt{1+e^{2\varphi}}}\big{)}\psi_{\beta}({x_{1}=0}),\\ &\Psi_{2}(0)=t_{12}\frac{e^{\varphi-i\theta_{2}}}{\sqrt{1+e^{2\varphi}}}\psi_{\alpha}({x_{2}=0})+t_{12}\frac{1}{\sqrt{1+e^{2\varphi}}}\psi_{\beta}({x_{2}=0}),\\ &\Psi_{3}(0)=t_{13}\frac{e^{\varphi-i\theta_{3}}}{\sqrt{1+e^{2\varphi}}}\psi_{\alpha}({x_{3}=0})+t_{13}\frac{1}{\sqrt{1+e^{2\varphi}}}\psi_{\beta}({x_{3}=0}).\end{aligned}\right.$$

$$\left\{\begin{aligned} &\frac{e^{-\varphi-i\theta_{1}}}{\sqrt{1+e^{-2\varphi}}}+r_{1}\frac{e^{\varphi-i\theta_{1}}}{\sqrt{1+e^{2\varphi}}}=t_{13}\frac{1}{\sqrt{1+e^{2\varphi}}}e^{i\alpha_{1}},\\ &t_{12}\frac{e^{\varphi-i\theta_{2}}}{\sqrt{1+e^{2\varphi}}}=\big{(}\frac{1}{\sqrt{1+e^{-2\varphi}}}+r_{1}\frac{1}{\sqrt{1+e^{2\varphi}}}\big{)}e^{i\alpha_{2}},\\ &t_{13}\frac{e^{\varphi-i\theta_{3}}}{\sqrt{1+e^{2\varphi}}}=t_{12}\frac{1}{\sqrt{1+e^{2\varphi}}}e^{i\alpha_{3}}.\end{aligned}\right.$$

$$S=\frac{1}{e^{3\phi }-e^{i\Lambda }}\left(\begin{array}{ccc}{\tilde{r}}_1& {\tilde{t}}_{12}& {\tilde{t}}_{13}\\ {\tilde{t}}_{21}& {\tilde{r}}_2& {\tilde{t}}_{23}\\ {\tilde{t}}_{31}& {\tilde{t}}_{32}& {\tilde{r}}_3\end{array}\right),\Lambda =\sum _i({\theta }_i+{\alpha }_i).$$

$$\frac{E}{m}=\cosh\varphi=\cos\bigg{(}\frac{\Lambda+2n\pi}{3}\bigg{)},~{}n\in\mathbb{Z},~{}N_{s}=-\frac{\sum_{i}\theta_{i}}{2\pi}.$$

$$\frac{E}{m}=\cosh\varphi=\cos\bigg{(}\frac{\Lambda}{l}\bigg{)},~{}N_{s}=-\frac{\sum_{i}\theta_{i}}{2\pi},$$

$$\mathcal{H}=-t_0\sum _{⟨i,j⟩}(c_i^{†}{c}_j+\text{H.c.})-t_1\sum _{⟨⟨i,j⟩⟩}(e^{-i{\varphi }_{ij}}{c}_i^{†}{c}_j+\text{H.c.}).$$

$$\phi =\frac{2{\sum }_i{\theta }_i-\pi /2}{3},E=\cos \phi ,N_s=-\frac{{\sum }_i{\theta }_i}{2\pi }.$$

$$\left\{\begin{aligned} \Psi_{1}(0)&=\bigg{(}\frac{e^{-\varphi-i\theta_{1}}}{\sqrt{1+e^{-2\varphi}}}+r\frac{e^{\varphi-i\theta_{1}}}{\sqrt{1+e^{2\varphi}}}\bigg{)}\psi_{\alpha}(x_{1}=0)+\bigg{(}\frac{1}{\sqrt{1+e^{-2\varphi}}}+\frac{r}{\sqrt{1+e^{2\varphi}}}\bigg{)}\psi_{\beta}(x_{1}=0),\\ \Psi_{2}(0)&=t\frac{e^{\varphi-i\theta_{2}}}{\sqrt{1+e^{2\varphi}}}\psi_{\alpha}(x_{2}=0)+t\frac{1}{\sqrt{1+e^{2\varphi}}}\psi_{\beta}(x_{2}=0).\end{aligned}\right.$$

$$\{\begin{array}{rl}\frac{e^{-\phi -i{\theta }_1}}{\sqrt{1+e^{-2\phi }}}+r\frac{e^{\phi -i{\theta }_1}}{\sqrt{1+e^{2\phi }}}& =t\frac{1}{\sqrt{1+e^{2\phi }}},\\ \frac{1}{\sqrt{1+e^{-2\phi }}}+r\frac{1}{\sqrt{1+e^{2\phi }}}& =t\frac{e^{\phi -i{\theta }_2}}{\sqrt{1+e^{2\phi }}}.\end{array}$$

$$\{\begin{array}{rl}r& =r^{\prime }=e^{\phi }\frac{e^{i({\theta }_1+{\theta }_2)}-1}{e^{2\phi }-e^{i({\theta }_1+{\theta }_2)}},\\ t& =e^{i{\theta }_2}\frac{e^{2\phi }-1}{e^{2\phi }-e^{i({\theta }_1+{\theta }_2)}},\\ t^{\prime }& =e^{i{\theta }_1}\frac{e^{2\phi }-1}{e^{2\phi }-e^{i({\theta }_1+{\theta }_2)}}.\end{array}$$

$$|r{|}^2+|t{|}^2=\frac{e^{2\phi }(2-2\cos ({\theta }_1+{\theta }_2))+e^{4\phi }+1-2e^{2\phi }}{e^{4\phi }+1-2e^{2\phi }\cos ({\theta }_1+{\theta }_2)}=1.$$

$$S=\left(\begin{array}{cc}t& r\\ r^{\prime }& t^{\prime }\end{array}\right)=\frac{1}{e^{2\phi }-e^{i({\theta }_1+{\theta }_2)}}\left(\begin{array}{cc}\tilde{t}& \tilde{r}\\ {\tilde{r}}^{\prime }& {\tilde{t}}^{\prime }\end{array}\right).$$

$$e^{i{\theta }_1+i{\theta }_2}-e^{2\phi }=0,2\phi =i({\theta }_1+{\theta }_2+2n\pi ),n\in \mathbb{Z},$$

$$E=\cosh\varphi=\cos\bigg{(}\frac{\varphi}{2}\bigg{)},~{}N_{s}=-\frac{|\theta_{2}+\theta_{1}|}{2\pi}.$$

$$\eta =\frac{r}{t}=\frac{(e^{i{\theta }_1}-e^{-i{\theta }_2})}{(e^{\phi }-e^{-\phi })}.$$

$$\arg(\eta)=\arctan\bigg{(}\frac{\sin\theta_{1}+\sin\theta_{2}}{\cos\theta_{1}-\cos\theta_{2}}\bigg{)}=\frac{\pi}{2}-\frac{\theta_{1}-\theta_{2}}{2},$$

$$|\eta {|}^2=\frac{2-2\cos {\theta }_1\cos {\theta }_2+2\sin {\theta }_1\sin {\theta }_2}{e^{2\phi }+e^{-2\phi }-2}=\frac{{\sin }^2(\frac{{\theta }_1+{\theta }_2}{2})}{{\sinh }^2\phi }.$$

Full text

I Introduction

A central feature of topological insulators (TI) is the bulk-edge correspondence: a $d$-dimensional TI with given symmetries has a bulk energy gap but symmetry protected gapless $d-1$ dimensional boundary excitations Kane and Mele (2005); Fu et al. (2007); Moore and Balents (2007); Fu and Kane (2007); Qi and Zhang (2011); Hasan and Kane (2010); Moore (2010). Recent studies on higher-order TIs generalized this bulk-edge correspondence. An $n$-th order TI has protected gapless modes of co-dimension $n$ Benalcazar et al. (2017a, b); Langbehn et al. (2017); Ezawa (2018); Schindler et al. (2018a); Song et al. (2017); Khalaf et al. (2017); Khalaf (2018); Fang and Fu (2017). A two-dimensional (2d) second order topological insulator (2d SOTI), for instance, is an insulator with gapped edge but gapless corners Benalcazar et al. (2017a, b); Langbehn et al. (2017); Ezawa (2018), i.e., there are localized in-gap states at corners under open boundary conditions. The higher order TIs can be derived from gapping out boundary Hamiltonian Langbehn et al. (2017); Khalaf et al. (2017); Khalaf (2018); Fang and Fu (2017). More specifically, to obtain a 2d SOTI, one can gap out a single helical edge state Benalcazar et al. (2017a, b); Langbehn et al. (2017), or alternately a pair of coupled counter-propagating helical edge states Langbehn et al. (2017); Zhu (2018); Yan et al. (2018). The point of this paper is to develop an effective theory to describe coupled edge states more generally and their dependence on the topology of the system boundary, which allows a description of the domain-wall states that remain at the intersection of edges for various types of edge junctions.

Meanwhile, one can think of the connected problem of higher order TIs. If we put an ordinary 2d TI on a closed surface of some 3d manifold, is it possible to have gapped 2d faces and 1d edges, but gapless 0d corner modes? Topological fullerenes Rüegg et al. (2013) are an example of this kind of system. They are polyhedral surfaces wrapped by the Haldane honeycomb lattice model Haldane (1988), leaving wedge disclination defects at the vertices Rüegg and Lin (2013); Rüegg et al. (2013). While these fullerenes do not currently exist in nature, very recent experiments indicate that twisted bilayer graphene at small twist angle supports a network of domain walls with threefold junctions (“Y-junctions”) Huang et al. (2018); Wu et al. (2018). These domain walls Ju et al. (2015) are not strictly topologically protected but conductance is expected to be high at the length scales of this network. If the planar system has non-vanishing Chern number, these topological fullerenes have gapped bulk and hinge states (here a “hinge state” is localized at the intersection of two 2d surfaces), but characteristic corner-localized in-gap states. These corner states can be related to the existence of nontrivial defect states bound to isolated wedge disclinations Teo and Hughes (2013); Gopalakrishnan et al. (2013); Benalcazar et al. (2014). The connection between the fullerene problem and the 2d SOTI can be viewed as follows: the classification of 2d SOTI is derived from that of TIs in 1d, which is identical to the classification of co-dimension 2 topological defects Teo and Kane (2010); Chiu et al. (2016); Langbehn et al. (2017). This implies that the topological fullerenes and certain classes of 2d SOTI should be describable in the same framework. The emergence of states bound to defects (such as disclination or dislocation) has previously been explained in several cases by edge soliton theory, i.e., the effective theory for a pair of coupled counter-propagating helical edge states Lee et al. (2007); Qi et al. (2008); Ran et al. (2009); Rüegg and Lin (2013); Klinovaja and Loss (2015). Although this theory is able to predict the fractional charge bound to the (edge) soliton Su et al. (1980); Goldstone and Wilczek (1981); Jackiw and Semenoff (1983) in those examples, one needs to extend the approach in order to incorporate crystalline symmetries in more complicated systems and obtain faithful bound state energies. (Note that in a system of noninteracting fermions, fractional charge should be thought of as an offset or displacement of the charge density, rather than as a property of elementary excitations.)

In this article, we propose a generic edge network theory to capture the boundary topology of coupled edge states. As a consequence of the bulk-edge correspondence, the edge states carry the necessary information of their topological insulator parents. By assigning proper boundary conditions on edge states at their vertices, the edge networks correctly predict the existence of bound states (edge solitons) and other information. We further considered edge states living on the hinges of varies 3d manifolds, where the edge states are generated from topological insulators attached on corresponding surfaces. Such edge networks can faithfully predict the energy and fractional charge of bound states located at the vertices, going beyond previous edge soliton theories. These edge networks are shown here to capture the key properties of topological fullerenes as well as some 2d SOTIs, and it is hoped that they will be useful for other problems as well.

The rest of the paper is organized as follows: In Sec. II, we briefly review basic facts and notation for an edge network made from multiple pairs of coupled helical edge states. In Sec. III, we discuss the minimal edge network constrained to lie on a closed 1d loop, and show the existence of bound state with fractional charge in the presence of certain symmetries. Based on this we further propose an AI class 2d SOTI that can be easily realized with atoms in an optical lattice. In Sec. IV, we consider edge networks with a multi-leg vertex. We first derive the bound state energy and charge for a Y-junction via a scattering matrix approach in Sec. IV.1. Then, in Sec. IV.2, we apply the results in Sec. IV.1 to the tetrahedral topological fullerene as an example. Starting from edge networks, we map the tetrahedral topological fullerene to the 2d SOTI we proposed. We summarize the main results in Sec. V with an eye toward future developments and applications of this picture.

II Description of edge network

We start with several pairs of coupled helical edge states, e.g., living on the hinges of the 3d manifold shown in Fig.[1.(a)]. The network is described by the effective Hamiltonian:

[TABLE]

Here, $i$ labels the hinges, and $x_{i}$ is the coordinate set along a specific hinge. The two component wave-function $\Psi(x_{i})=(\psi_{\alpha}(x_{i}),\psi_{\beta}(x_{i}))^{T}$ denotes a pair of coupled counter-propagating helical edge states living on $i$-th hinge, and varies smoothly on the scale of the lattice constant. The magnitude of edge velocity $v$ is set identical for all edge states, and their directions should be compatible with the positive direction of $x_{i}$. The mass term ${\mathcal{M}_{i}(\theta_{i})}=m\cos\theta_{i}\sigma_{x}+m\sin\theta_{i}\sigma_{y}$ describes the coupling on hinge $x_{i}$, where $\sigma_{x,y,z}$ are Pauli matrixes. Without loss of generality, we assume that $m\geq 0$ and $0\leq\theta_{i}\leq 2\pi$. If $m=0$, the helical edge states are decoupled and their energy spectrum is gapless. A non-zero mass term can locally gap out a pair of edge states, which is the situation that we are interested in.

To the Hamiltonian we need to add proper boundary conditions for these edge states at vertices where two or more edges come together. The boundary condition describe the scattering process at the junction. By doing so we can solve Eq.[1] and predict the existence of localized edge solitons that lie in the (bulk and edge) gaps, as well as their properties.

Before discussing edge network on specific configuration, we point out that the Hamiltonian Eq.[1] may be generalized into the case of Helical Luttinger liquid  Wu et al. (2006); Xu and Moore (2006); Hou et al. (2009); Giamarchi :

[TABLE]

The noninteracting Hamiltonian $H_{0}^{i}$ on each hinge can be divided into two parts: the linearized free Dirac field $H_{0,1}^{i}$ and their coupling ($H_{0,2}^{i}$) with two real-valued classical scalar field $\lambda_{1,2}(x_{i})$ Chamon et al. :

[TABLE]

Here $\psi_{\alpha(\beta),i}$ ($\lambda_{1(2),i}$) is short for $\psi_{\alpha(\beta)}(x_{i})$ ($\lambda_{1,2}(x_{i})$). Compared with ${\mathcal{M}}_{i}(\theta_{i})$ in Hamiltonian Eq.[1], we find that $\lambda_{1,i}=m\cos\theta_{i}$ and $\lambda_{2,i}=m\sin\theta_{i}$. For helical Luttinger liquid, we only need to consider the forward scattering $H_{\text{int},2}^{i}$ and chiral interaction $H_{\text{int},4}^{i}$ Wu et al. (2006); Xu and Moore (2006); Hou et al. (2009); Giamarchi :

[TABLE]

where $g_{2,i}$ and $g_{4,i}$ are interacting constants. One can conduct the standard bosonization procedure for Hamiltonian Eq.[2] by defining bosonic field $\partial_{x_{i}}\tilde{\varphi}_{i}=-\pi[\rho_{\alpha}(x_{i})+\rho_{\beta}(x_{i})]$ and $\partial_{x_{i}}\tilde{\theta}_{i}=\pi[\rho_{\alpha}(x_{i})-\rho_{\beta}(x_{i})]$, where $\rho_{\alpha(\beta)}(x_{i})$ stands for the density for counter-propagating edge states, i.e. $\rho_{\alpha(\beta)}(x_{i})=\psi^{\dagger}_{\alpha(\beta),i}\psi_{\alpha(\beta),i}$. The $\tilde{\theta}_{i}$ here should be distinguished from $\theta_{i}$ in effective mass ${\mathcal{M}}$. The Bosonized Hamiltonian for each hinge, $H_{B}^{i}=H_{B,0}^{i}+H_{B,1}^{i}$ reads:

[TABLE]

Here, $u\equiv v\sqrt{(1+g_{4}/2)^{2}-(g_{2}/2)^{2}}$ is the velocity, $K=\sqrt{(1+g_{4}/2-g_{2}/2)/(1+g_{4}/2+g_{2}/2)}$ is the Luttinger parameter, and $a$ is the lattice constant whose inverse stands for the momentum cut off of vacuum Giamarchi ; Chamon et al. ; Shankar (2017). The Hamiltonian $H_{B}^{i}$ is also interacting, and the interaction $H_{B,1}^{i}$ can be minimized by set $\tilde{\varphi}(x_{i})=\theta(x_{i})+\pi$. Referring to the bonsonized conserved current $j^{\mu}_{i}=\epsilon^{\mu\nu}\partial_{\nu}\tilde{\varphi}(x_{i})/2\pi\approx\epsilon^{\mu\nu}\partial_{\nu}\theta(x_{i})/2\pi$, for the simplest two terminal junction with two legs $x_{1,2}$ (see in Fig.[1.(b)]), the topological charge $\hat{Q}$ is given by Chamon et al. ; Fröhlich and Marchetti (1988):

[TABLE]

A mass kink of ${\mathcal{M}}_{i}(\theta_{i})$ implies nonzero topological charge $\hat{Q}$, see in Fig.[1.(b)]. This is in accordance with the soliton charge $N_{s}$ derived from non-interacting Fermionic theory Su et al. (1980); Goldstone and Wilczek (1981); Jackiw and Semenoff (1983), see also Eq.[7] in later on Sec. III. For simplicity, in the rest of our article we will focus on the non-interacting model Eq.[1]. It is reasonable to believe that the value of soliton charge remains unchanged when turning on interaction because it can be calculated from properties away from the junction. However, the response of bound state energy with respect to external flux may be modified by interaction, and may need a deeper description, e.g., by boundary conformal field theory Oshikawa et al. (2006); Hou et al. (2012).

III Edge states on closed 1d loop

We first consider the minimal example of an edge network, a pair of helical edge states living on the boundary of a closed 1d loop, as shown in Fig.[1.(c)]. The point is to determine how symmetries fix the free coefficients introduced in the previous discussion. The basis is chosen as $\Psi(x_{i})=(\psi_{\alpha}(x_{i}),\psi_{\beta}(x_{i}))^{T}$, where $\psi_{\alpha}(x_{i})$ ($\psi_{\beta}(x_{i})$) denotes the chiral edge states propagating in the clockwise (anti-clockwise) direction. We set four coordinates $x_{i=1,2,3,4}\geq 0$, and define $x_{5}=x_{1}$. The coupling for edge states on each leg is given by an effective mass ${\mathcal{M}}_{i}(\theta_{i})$, where we have set $v=m=1$ for simplicity. We use a set of trial wave functions $\Phi(x_{i})_{o(e)}=\frac{1}{\sqrt{N_{o(e)}}}\exp(-|(x_{i}-x_{i}^{o(e)})\sin\varphi|)\chi(x_{i})_{o(e)}$ to look for bound states localized the origin (o) and end (e) of $i$-th edge, with $\chi(x_{i})_{o(e)}=(a_{i}^{o(e)},b_{i}^{o(e)})^{T}$. Here $a_{o(e)},b_{o(e)}$, $\varphi$ and normalization constant $1/\sqrt{N_{\pm}}$ are coefficients to be determined.

Substituting the trial wave function $\Phi{(x_{i})}_{o(e)}$ into Eq.[1] for each individual edge, we find modes localized at two ends of $i$-th edge. For the states at the origin of $i$-th edge, we have the wave function $\chi(x_{i})_{o}=e^{i\delta_{i}^{o}}(e^{i(\varphi-\theta_{i})},1)^{T}$ with energy $\epsilon_{i}^{o}=\cos\varphi$. For the states at the end, we have $\chi(x_{i})_{e}=e^{i\delta^{e}_{i}}(e^{-i(\varphi+\theta_{i})},1)^{T}$ with energy $\epsilon_{i}^{e}=\cos\varphi$. Here, $\delta_{i}^{o,e}$ are overall phase factors. The wave-function we solved previously should satisfy the boundary condition at the corner, i.e., $\Phi({x_{i+1}\rightarrow x_{i+1}^{o}})_{o}=\Phi({x_{i}\rightarrow x_{i}^{e}})_{e}$ and $\epsilon_{i+1}^{o}=\epsilon_{i}^{e}$. If $\theta_{i}=\theta_{i+1}$, the only allowed solution is $\varphi=0$, which means that the localization length $\xi=1/|\sin\varphi|\rightarrow\infty$ and no bound state exists. If $\theta_{i}\neq\theta_{i+1}$, we have a mass kink at the intersection of $i$-th and $i+1$-th edge. The solution corresponds to an un-paired edge soliton Lee et al. (2007) localized at the intersection, with energy and fractional fermion number Su et al. (1980); Goldstone and Wilczek (1981); Jackiw and Semenoff (1983) given by:

[TABLE]

Since we measure the charge with respect to the vacuum, there is a minus sign for the soliton charge $N_{s}$. Eq.[7] predicts the existence of a domain wall state for any two adjoint edges. More specifically, the edge soliton derived from the aforementioned effective theory can be used to explain fractional charge in varies systems, such as the bound states induced by magnetic domain wall in the quantum spin hall effect Qi et al. (2008), or the localized state bound to 2d disclination (dislocation) defect in topological insulators Rüegg and Lin (2013); Ran et al. (2009).

The minimal edge network can explain the corner states in at least some kinds of 2d SOTI. The 2d SOTIs have gapped bulk and edges, but gapless corners. They can be derived from gapping out topological edge states. Heuristically, one potential way to get a 2d SOTI is by stacking 1d TIs, making the 0d boundaries of these 1d TIs form another set of 1d TIs in the perpendicular direction. This is one way to obtain the quadrupole insulator Benalcazar et al. (2017a, b). Alternatively, one can couple a pair of (or more) counter-propagating helical edge states living on the boundary of 2d TI and gap them out. Here we will use the latter picture extensively. Crystalline symmetries with unitary symmetry operator $U$, such as reflection Langbehn et al. (2017), inversionKhalaf (2018) and rotation symmetry Song et al. (2017); Fang and Fu (2017); Khalaf et al. (2017), can constrain the distribution of effective mass term ${\mathcal{M}}_{i}(\theta_{i})$ on the boundary. On the edges compatible with crystalline symmetry, $[{\mathcal{M}}_{i}(\theta_{i}),U]=0$. If two adjoint edges are related by crystalline symmetry with operator $U$, then $U^{\dagger}{\mathcal{M}}_{i}(\theta_{i})U={\mathcal{M}}_{i+1}(\theta_{i+1})$. If ${\mathcal{M}}_{i}(\theta_{i})\neq{\mathcal{M}}_{i+1}(\theta_{i+1})$, a domain wall state emerges at the intersection of two adjoining edges, as demonstrated before.

Distinct from corner-localized zero modes in a 2d second-order topological superconductor Langbehn et al. (2017); Zhu (2018), we find that, in the absence of particle-hole symmetry and chiral symmetry Langbehn et al. (2017); Shiozaki and Sato (2014); Lau et al. (2016); Trifunovic and Brouwer (2017), one can have corner modes with non-zero energy. The system we consider has two reflection-symmetric axes, as shown in Fig.[1.(c)]. The reflection operator for the red axis is $U_{b}=\sigma_{x}$, while the reflection operator for the blue axis is $U_{r}=\sigma_{y}$. Edge $AB(x_{1})$ and $CD(x_{3})$ are reflection symmetric edges for $U_{b}$, thus the only symmetry-allowed mass term is $\pm\sigma_{x}$. Similarly edge $AD(x_{4})$ and $BC(x_{2})$ are reflection symmetric edges for $U_{r}$, thus the only symmetry-allowed mass term is $\pm\sigma_{y}$. In summary, the effective mass terms on four edges $x_{1,2,3,4}$ are:

[TABLE]

Referring to Eq.[7], we find $\varphi={\pi}/{4},~{}E=\cos\varphi=1/\sqrt{2}$, and $N_{s}=-1/4$ for each corner, corresponding to four edge solitons on the loop.

SOTIs have been claimed to be appear in various systems Serra-Garcia et al. (2018); Peterson et al. (2018); Imhof et al. (2017), including bismuth Schindler et al. (2018b). Based on recent progress of two-dimensional spin-orbit coupling in cold atom system Liu et al. (2014); Wu et al. (2016), we provide a feasible experimental proposal of 2d SOTI with edge mass distribution as Eq.[8]. By stacking two Chern insulator layers with opposite Chern numbers (which can be easily realized in experiments by adding a magnetic field with gradient), the 2d tight-binding Hamiltonian for our 2d SOTI model is:

[TABLE]

Here, $s=\pm$ stands for layer index. The positive $t_{\alpha=x,y}$ and $t_{\text{so}}$ denotes, respectively, the inner-layer spin conserved and spin-flip hopping. The $m_{z}^{s}$ represents an effective Zeeman term, with $m_{z}^{+}=m_{z}$ and $m_{z}^{-}=-m_{z}$, which can be realized by a magnetic field with gradient. The spin-flip hopping $t_{\text{so}}$ and $\lambda$ comes from the spin-orbit coupling induced by effective inner-layer and inter-layer Raman coupling, respectively. Transforming $H$ into the momentum space yields $H=\sum_{k,\sigma,\sigma^{\prime}}\hat{c}_{\bf{k},\sigma}^{\dagger}\mathcal{H}_{\sigma,\sigma^{\prime}}(k)\hat{c}_{\bf{k},\sigma^{\prime}}$, with

[TABLE]

where $\bm{\tau}$ and $\bm{\sigma}$ are Pauli matrices in spin space and layer space, respectively. If $\lambda=0$, the Hamiltonian Eq.[10] has particle-hole symmetry ${\mathcal{P}}=\tau_{1}\sigma_{3}{\mathcal{K}}$, time-reversal symmetry ${\mathcal{T}}=\tau_{2}\sigma_{2}{\mathcal{K}}$, and chiral symmetry ${\mathcal{S}}=\tau_{3}\sigma_{1}$, where ${\mathcal{K}}$ stands for complex conjugate. With $|m_{z}|<2t_{x}+2t_{y}$, the system can be viewed as a robust index spin hall effect Zhou et al. (2008). It also has two reflection symmetric axes along $x-,y-$ directions, denoted by operator $\hat{U}_{x}=\tau_{2}\sigma_{2}$, $\hat{U}_{y}=\tau_{1}\sigma_{1}$. A small but non-zero $\lambda$ breaks the chiral and particle hole symmetry, and gaps out the helical edge states from the original index spin hall effect. By projecting the low energy Hamiltonian of Eq.[10] into the helical edge states derived from $\lambda=0$, one can get the effective edge Hamiltonian identical to Eq.[8], leading to the similar set of gapped edges and gapless corners.

We further confirm the analytic results by numerically diagonalizing the Hamiltonian Eq.[9] for two different boundary conditions, as shown in Fig.[2.(a)]. We find four corner modes with non-zero energy for both patterns. Fig.[2.(b)] shows the energy spectrum close to the Fermi surface. The inter-layer coupling $\lambda$ opens a gap $E_{\text{gap}}\approx 2\lambda$ at the boundary, and we can see clearly four corner-localized in gap states. At half-filling, one out of four in-gap states is filled, which compensates the $-1/4$ defect charge at each corner. In realistic cold atom experiments, the detection of the fractional charge at the corner can be conducted by conventional single site resolution. By turning on an $s$-wave onsite interaction for atomsLiu et al. (2014), this model becomes a 2d second order topological superfluid.

IV Edge states in multi-leg junction

We now turn to study edge networks with multiple pairs of edge states coming together at a vertex (or equivalently a junction). Fig.[3.(a)] shows a Y-junction with six edge states living on three legs. For each semi-infinite axis $x_{i}(i=1,2,3)$, we use the $\psi_{\alpha}(x_{i})$ and $\psi_{\beta}(x_{i})$ to denote the outgoing and incoming chiral edge states for the $i$-th leg, respectively. Instead of matching wave functions by hand as in the previous minimal 1D edge network, here, we introduce a more generic scattering matrix approach: injecting a mode along a specified leg will lead to reflection and transmission after scattering at the junction, and the poles of scattering matrix implies the existence of bound states. We make the following assumptions to capture the scattering process: (1) Away from the junction, each chiral edge state should be identical to that of an isolated Chern insulator layer, at most up to a global phase factor; (2) During the scattering process, the edge states from the same Chern insulator layer should maintain their amplitude, but could capture a phase shift. The value of the phase shift depends on the details of scattering, but will be constrained by symmetries in specific examples.

IV.1 Scattering theory for Y-junction

For an isolated junction, the incoming and outgoing scattering modes can be described by combining incoming and outgoing chiral edge states. Different from localized state, for scattering state with momentum $k$, we denote $\eta=v/m$ and set $k\eta=\sinh\varphi>0$. Then for $x_{i}>0$, under the basis \big{(}\psi_{\alpha}(x_{i}),\psi_{\beta}(x_{i})\big{)}^{T}, for each individual leg, from Eq.[1] we can derive normalized wave function of incoming and outgoing modes as: $\psi_{\text{in}}^{T}(x_{i})=(e^{-\varphi-i\theta_{i}},1)^{T}/\sqrt{1+e^{-2\varphi}},\quad\psi_{\text{out}}^{T}(x_{i})=(e^{\varphi-i\theta_{i}},1)^{T}/\sqrt{1+e^{+2\varphi}}$, with corresponding energy $E/m=+\cosh\varphi$. If we inject a mode along negative $x_{1}$ direction, the wave function on leg $x_{1}$ is given by $\Psi_{1}(x_{1})=e^{-ikx_{1}}\psi_{\text{in}}(x_{1})+r_{1}e^{ikx_{1}}\psi_{\text{out}}(x_{1})$. Meanwhile, the wave function on leg $x_{2}$ is given by $\Psi_{2}(x_{2})=t_{12}e^{ikx_{2}}\psi_{\text{out}}(x_{2})$, and $\Psi_{3}(x_{3})=t_{13}e^{ikx_{3}}\psi_{\text{out}}(x_{3})$ for wave function on leg $x_{3}$. We have used $r_{1}$ for reflection coefficient on leg $x_{1}$, and $t_{12}~{}(t_{13})$ for transmission coefficient for the scattering from $x_{1}$ to $x_{2}~{}(x_{3})$. With this we can expand the wave function around the intersection as:

[TABLE]

In Fig.[3.(a)], the edge states in same color are from the same Chern insulator layer. Due to the continuity of edge state wave function for each individual layer, we have $\psi_{\alpha}(x_{1}\rightarrow 0^{+})=\psi_{\beta}(x_{3}\rightarrow 0^{+})$, $\psi_{\beta}(x_{1}\rightarrow 0^{+})=\psi_{\alpha}(x_{2}\rightarrow 0^{+})$, and $\psi_{\beta}(x_{2}\rightarrow 0^{+})=\psi_{\alpha}(x_{3}\rightarrow 0^{+})$. During the scattering process they can capture an additional phase factor $e^{i\alpha_{i}}$, which depends on the details of the scattering process. This leads to:

[TABLE]

With this we can solve $r_{1}$, $t_{12}$ and $t_{13}$ in the term of $\varphi,\alpha_{i}$, and $\theta_{i}$. By injecting modes along the negative directions of rest two legs (see in Appendix), we can derive whole coefficients for the scattering matrix $S$:

[TABLE]

For arbitrary scattering process, $\Psi_{\text{out}}=S\Psi_{\text{in}}$, where \Psi_{\text{in(out)}}^{T}=\big{(}\psi(x_{1}),~{}\psi(x_{2}),~{}\psi(x_{3})\big{)}_{\text{out(in)}}. The pole of the scattering matrix, $e^{3\varphi}-e^{i\Lambda}=0$, implies the existence of bound state localized at the junction. Note that, in the presence of edge soliton, each of these semi-infinite legs contributes a fractional charge $-\theta_{i}/2\pi$. With these we find:

[TABLE]

As we mentioned before, $\Lambda=\sum_{i}(\theta_{i}+\alpha_{i})$, which depends on the details of scattering. The energy-phase relation Eq.[14] for 3-leg Y-junction can be easily generalized to $l$-leg junction:

[TABLE]

where we have let $2n\pi$ be absorbed into $\Lambda$ for latter convenience.

IV.2 Application to topological fullerenes

The multi-leg edge junction can be used to describe the bound state in an isolated wedge disclination Rüegg and Lin (2013); Rüegg et al. (2013); Teo and Hughes (2013); Gopalakrishnan et al. (2013); Benalcazar et al. (2014), which is the building block of topological fullerenes Rüegg and Lin (2013); Rüegg et al. (2013). More specifically, the Y-junction edge network mentioned above can be used to analyze one vertex of tetrahedral topological fullerenes (as shown in Fig.[1.(a)]), which is a wedge disclination defect with Frank index $f=3$ (or $180^{\circ}$ Frank angle). The Frank index $f$ here stands for the number of $60^{\circ}$ Chern insulator layers taken away from the complete Haldane lattice. In order to build the edge network for such a disclination, let us first consider three $60^{\circ}$ semi-infinite triangular layers (A,B,C) of Haldane honeycomb lattice coming together, as shown in Fig.[3.(d)]. Each layer is coupled with its two neighbors across the seams. The tight-binding Hamiltonian for such a disclination is given byHaldane (1988); Rüegg and Lin (2013); Rüegg et al. (2013):

[TABLE]

Here, $c^{\dagger}_{i}$ ($c_{i}$) is creation (annihilation) operator for spinless fermion on $i$-th site. The $t_{0}$ and $t_{1}$ denotes, respectively, the nearest-neighbor hopping and next-nearest-neighbor hopping amplitudes. The $e^{i\phi_{ij}}$ provides an additional phase factor for next-nearest-neighbor hopping. Within the topological region, each individual layer can provide chiral edge states surrounding the bulk. The local Chern vectorBianco and Resta (2011) for each layer points outside the plane of the paper, which ensures six edge states propagating according to the pattern in the figure.

These six edge states are not independent. The blue (red, gray) edge states $1,2$ ($3,4$; $5,6$) come from the same triangular layer, and they are connected by $\psi_{2j}(x_{j}\rightarrow 0^{+})=\psi_{2j-1}(x_{j+2}\rightarrow 0^{+})$. If an edge state is coupled with its time-reversal counterpart across the seam, we say this seam does not have phase mismatch. The total wave function on a lattice site across the seam is given by $\varphi_{\text{edge},\alpha}(x_{i})=e^{ik_{E}a}\psi_{\alpha}(x_{i})$ and $\varphi_{\text{edge},\beta}(x_{i})=e^{-ik_{E}a}\psi_{\beta}(x_{i})$, respectively, where $k_{E}$ denotes the edge momentum and $a$ stands for the lattice constant. The $\psi_{\alpha,\beta}(x_{i})$ here should be understood as the edge states on corresponding sub lattice. The effective coupling between two states is $\int d\tau\lambda\varphi^{*}_{\text{edge},\alpha}(x_{i})\varphi_{\text{edge},\beta}(x_{i})$, with $\lambda$ stands for the bond across the seam. The integral is done within a unit cell. For an isolated disclination, the total phase mismatch $\sum_{i}\theta_{i}$ for all legs (seams) is fixed in the absence of external flux. Due to the quantization of charge pumping, the function $\sum_{i}\alpha_{i}$ should be linear to $\sum_{i}\theta_{i}$, i.e. $\sum_{i}\alpha_{i}=A\sum_{i}\theta_{i}+B$. The coefficients $A,B$ are related to the parameters from the tight-binding model, such as the effective radius $\rho$ and the Haldane gap $m=3\sqrt{3}t_{1}$Rüegg and Lin (2013). By comparing with the results from exact diagonalizing the tight-binding Hamiltonian Eq.[16] (see in Appendix), we find that, for Haldane gap $m\approx t_{0}=1$,

[TABLE]

Similarly, for the vertex of an octahedral topological fullerene, the number of legs is $l=6-2=4$, and we further have $\varphi=2\sum_{i}\theta_{i}/4$, $E=\cos\varphi$, and $N_{s}=-\sum_{i}\theta_{i}/2\pi$ (with $1\leq i\leq 4$). For the vertex of an icosahedral topological fullerene, the number of legs is $l=6-1=5$, and we further have $\varphi=(2\sum_{i}\theta_{i}+\pi/2)/5$, $E=\cos\varphi$, $N_{s}=-\sum_{i}\theta_{i}/2\pi$ (with $1\leq i\leq 5$).

We now turn to determine the value of $\theta_{i}$ for each legRüegg and Lin (2013), especially for the cases with external flux. Let us first consider the process of combining two smaller $60^{\circ}$-layers in Fig.[3.(b)] to a larger $120^{\circ}$-layer in Fig.[3.(c)]. The two smaller layers are cut from the same Haldane honeycomb lattice model, and they are next to each other in the original lattice. With the open boundary condition, both of them can hold chiral edge states, which are denoted by red and blue arrows in Fig.[3.(b)]. We can set a simultaneous coordinate for both layers across the seam, thus the total wave function on a lattice site on the blue (red) edge is $\varphi_{\text{edge},\alpha}(\xi)=e^{-ik_{E}\xi}\psi_{\alpha}(\xi)$ ($\varphi_{\text{edge},\beta}=e^{ik_{E}\xi}\psi_{\beta}(\xi)$). In the presence of inversion symmetry, $k_{E}a=\pi$ for Haldane honeycomb lattice model Rüegg and Lin (2013). Thus the base functions $e^{ik_{E}\xi}$ oscillates with a period of two sites. In order to glue two layers back to a larger layer without phase mismatch across the seam, the amplitudes should be in the pattern in Fig.[3.(b)]. The edge states on the decoupled two branches can be written separately as $\varphi_{\text{edge},\alpha}=e^{-ik_{E}\xi}\psi_{\alpha}$ for the lower branch of blue edge, and as $\varphi_{\text{edge},\beta}=e^{ik_{E}\xi}\psi_{\beta}$ for the right branch of the red edge. However, as shown in Fig.[3.(c)], the edge states has an additional phase shift when bypassing the corner. This leads to $i\psi_{\alpha}=-\psi_{\beta}$ or $\psi_{\alpha}=i\psi_{\beta}$. Thus we have that the proper phase difference across the seam should be $\pm i$. To avoid any ambiguity induced by the gauge chosen for wave functions, we define the effective mass term ${\mathcal{M}}_{i}(\theta_{i})$ on each leg with respect to the scenario without phase mismatch. Thus if there is no phase mismatch on a certain leg, then ${\mathcal{M}}_{i}(\theta_{i}=0)=m\sigma_{x}$.

Note that, for a wedge disclination with Frank index $f=3$ in Fig.[3.(d)], if we glue AB and BC across the seam as shown in Fig.[3.(b)], the system can be viewed as a Haldane honeycomb lattice on the half plane. The gluing process means that we have chosen to measure the relevant phase factor of edge states on all layers from $\psi_{1}$ with a fixed reference point. Thus the coupling across the seams $AB$ and $BC$ should not have a phase mismatch, thus ${\mathcal{M}}_{1}(0)={\mathcal{M}}_{2}(0)=m\sigma_{x}$. However, the lower boundaries of A and C has phase mismatch and ${\mathcal{M}}_{3}(\pi/2)=m\sigma_{y}$ Rüegg and Lin (2013). Finally, referring to Eq.[17], we have $\varphi=\pi/6$, $E=\cos\pi/6$ and $N_{s}=-1/4$ for the vertex of Tetrahedral topological fullerene. Eq.[17] also stands in the presence of external flux. Adding an external flux $\Phi$ opposite to local Chern vector at the center of junction is equivalent to change the coupling pattern with additional phase factor $e^{i\Phi}$ for the bond across the Dirac string Rüegg and Lin (2013); Rüegg et al. (2013). For simplicity we can put the Dirac string along $x_{3}$, thus $\theta_{3}=\phi+\pi/2$ and Eq.[17] can be written as $\varphi=2\Phi/3+\pi/6$. More specifically, if $\Phi=\pi/2$, we have $\sum_{i}\theta_{i}=\pi$ and $\varphi=\pi/2$. Thus the external flux $\Phi=\pi/2$ moves the bound state energy to $E=\cos\varphi=0$, as well as the fractional charge to $N_{s}=-\sum_{i}\theta_{i}/2=-1/2$. This is consisted with the analysis from symmetry: an external flux with $\Phi=\pi/2$ can restore the particle hole symmetry of the system Rüegg et al. (2013). Thus the bound state energy should be [math] and the fractional charge should be $-1/2$. Similar results apply for vertices of octahedral and icosahedral topological fullerenes (see in Appendix), and are in accordance with numerical results Rüegg and Lin (2013); Rüegg et al. (2013).

The corner states in topological fullerenes can be further explained by the edge networks with a group of multi-leg junctions. In Fig.[4.(a)] we plot the edge network for the tetrahedral topological fullerene, with ${\mathcal{M}}_{1}(0)=m\sigma_{x}$, ${\mathcal{M}}_{2}(\pi/2)=m\sigma_{y}$, ${\mathcal{M}}_{3}(\pi)=-m\sigma_{y}$, ${\mathcal{M}}_{4}(3\pi/2)=-m\sigma_{x}$, and ${\mathcal{M}}_{5}(0)={\mathcal{M}}_{6}(0)=m\sigma_{x}$. However, Eq.[17] is derived for an isolated vertex with all coordinates point outward, which is slightly different from the settings in Fig.[4.(a)]. Note that, for a pair of helical edge states living on a $i$-th hinge with effective mass ${\mathcal{M}}_{i}(\theta_{i})$, changing the direction of coordinates is equivalent to changing the mass term to \tilde{\mathcal{M}}_{i}\big{(}(-\theta_{i})~{}\text{mod}~{}2\pi\big{)}. Thus for each individual vertex, we can first flip the coordinates to the pattern in Fig.[3.(a)], by then using Eq.[17] we find four corner localized states with $E=\cos\pi/6$ and $N_{s}=-1/4$.

V Conclusion

We have constructed a generic edge network theory and shown its ability to capture the boundary topology of coupled edge states with different geometric constraints. We first discussed the minimal edge network on a closed 1d loop, and demonstrated that crystalline symmetry can produce spatial-dependent mass term, leading to the domain wall states at the intersection of adjoint edges. After discussing a model 2d second-order TI, we constructed edge networks for multi-leg junctions, which can faithfully reflect the properties of bound states in disclination defects. The edge network can include polyhedral hinges, which allows determination of the corner states in topological fullerenes. These results can help to understand the origin of topologically generated localized states in a variety of situations.

We can view the similarities between the 2D second-order TI and the 3D topological fullerine as reflecting the fact that the classification of 2d SOTI is derived from that of TIs in 1d, which is the same as classification of co-dimension 2 topological defects Teo and Kane (2010); Chiu et al. (2016); Langbehn et al. (2017), including point defects in surfaces. In this sense, the 2d SOTI we proposed is in the same topological class as a corresponding system with wedge disclination defects. Based on effective edge theory, we can map the topological fullerenes to the 2d SOTI Eq.[10] derived from gapping out helical edge states in Sec. III. For any polyhedron, one can traverse all the vertices along hinges without repeats. The traversal forms a closed 1d loop (see the blue thick arrows in Fig.[4.(a)]). We can cut the polyhedron into two congruent Chern insulator layers along the traversal, as shown in Fig.[4.(b)]. The two Chern insulator layers can be viewed as a “twisted” index spin hall effect. The edges on the closed 1d loop are gapped out by the gluing process, and the effective mass changes after bypassing each corner due to crystalline symmetries, leading to an edge soliton with fractional charge located at the corner. This is identical to the generation of fractional charge in our 2d SOTI model. Similarly, we can also map the octahedral and icosahedral topological fullerenes to (less natural) 2d SOTIs.

More generally, the networks of edges discussed here could be generalized to incorporate proximity-induced superconductivity or Luttinger liquid corrections, or conceivably to include additional localized degrees of freedom such as boundary Majorana states or spins as in previous studies of the Kondo effect in Y-junctions Oshikawa et al. (2006). In the cases discussed here, there are enough symmetries or other physical constraints to determine the key properties of the localized states in an edge network quite directly, while in other situations the properties such as fractional offset charges might be actively tuned by symmetry-breaking perturbations. Planar networks of helical edges and three-leg junctions have recently been discovered in bilayer graphene at small twist angles, which suggests that the study of edge networks is likely to become increasingly relevant to experiment.

ACKNOWLEDGEMENTS

We thank Daniel E. Parker and Takahiro Morimoto for useful conversations. This work was supported as part of the Center for Novel Pathways to Quantum Coherence in Materials, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences. J.E.M. acknowledges additional support from a Simons Investigatorship.

APPENDIX

0.1 Edge network for two-terminal junction

In this section, we look into the scattering theory of the simplest two-terminal junction. In order to keep in accordance with the scattering theory in Sec. IV.1, we set the positive direction of the two legs being opposite to each other and pointing outside the junction. This switches the $\theta_{1}$ to $-\theta_{1}$ compared with the notation in Sec. III.

We construct the conventional scattering theory as following: suppose we have a wave injected along the negative $x_{1}$ direction. The wave function on leg $x_{1}$ is given by $\Psi_{1}=e^{-ikx_{1}}\psi_{1,\text{in}}+re^{ikx_{1}}\psi_{1,\text{out}}$. Meanwhile, the wave function on leg $x_{2}$ is given by $\Psi_{2}=te^{ikx_{2}}\psi_{2,\text{out}}$. Here $r$ and $t$ stand for the reflection and transmission coefficients, respectively. Different from a localized state, for a scattering state with momentum $k$, we denote $\eta=v/m$ and set $k\eta=\sinh\varphi>0$. Then for $x_{i}>0$, using the basis \big{(}\psi_{\alpha}(x_{i}),\psi_{\beta}(x_{i})\big{)}^{T}, for each individual leg, from Eq.[1] we can derive a normalized wave function of incoming and outgoing modes as: $\psi_{\text{in}}^{T}(x_{i})=(e^{-\varphi-i\theta_{i}},1)^{T}/\sqrt{1+e^{-2\varphi}},\quad\psi_{\text{out}}^{T}(x_{i})=(e^{\varphi-i\theta_{i}},1)^{T}/\sqrt{1+e^{+2\varphi}}$, with corresponding energy $E/m=+\cosh\varphi$. We can expand the wave function around the junction by the combination of incoming and outgoing edge states:

[TABLE]

As shown in Fig.[A1], $\psi_{\alpha}(x_{1})=\psi_{\beta}(x_{2})$ and $\psi_{\beta}(x_{1})=\psi_{\alpha}(x_{2})$ since they are from the same Chern insulator. For the SOTI Eq.[10], the wave function should be continuous at the junction:

[TABLE]

By solving this we derive:

[TABLE]

The reflection and transmission coefficients $r$ and $t$ satisfy conservation of probability current:

[TABLE]

Finally we have the scattering matrix for two terminal junction as:

[TABLE]

One can easily check that the scattering matrix is unitary $S^{\dagger}S={\mathbf{1}}$. The coefficients of scattering matrix, see in Eq.[A3] has simultaneous poles:

[TABLE]

which stands for bound states localized at the junction with energy and fractional charge as:

[TABLE]

Remember that $\theta_{1}$ here is equal to $-\theta_{1}$ in Sec. III due to the flipping of $x_{1}$-leg’s direction, the above results is in accordance with Eq.[7]. We further define $\eta$ as:

[TABLE]

The argument and the absolute value of $\eta$ are:

[TABLE]

[TABLE]

Thus the bound state energy can also be parametrized by reflection and transmission coefficents as:

[TABLE]

0.2 Edge networks for a Y-junction

In this section, we provide more details about the edge network description of a three-leg junction (“Y-junction”).

1.2.1 Bound states from matching wave function

Different from the scattering matrix approach in Sec. IV.1, here we get the same results by matching the trial wave function and validate that the poles of scattering states do correspond to localized states. We can derive the trial wave function by using the similar method in Sec. III. Substitute the trial wave function $\chi{(x_{i})}$ in to Eq.[1] for each leg independently, we find the modes localized at two ends of $i$-th edge, with $\chi{(x_{i})}=e^{i\delta_{i}}(e^{i(\varphi-\theta_{i})},1)^{T}$ for energy $\epsilon_{i}^{o}=\cos\varphi$. This gives the relation between $a_{i}$ and $b_{i}$ on the same leg. More specifically: for the leg 1, we have $\chi{(x_{1})}=e^{i\delta_{1}}(e^{i(\varphi-\theta_{1})},1)^{T}$, with the basis $\Psi(x_{1})=(\psi_{2}(x_{1}),\psi_{3}(x_{1}))^{T}$; for the leg 2, we have $\chi{(x_{2})}=e^{i\delta_{2}}(e^{i(\varphi-\theta_{2})},1)^{T}$, with the basis $\Psi(x_{2})=(\psi_{4}(x_{2}),\psi_{5}(x_{2}))^{T}$; for the leg 3, we have $\chi{(x_{3})}=e^{i\delta_{3}}(e^{i(\varphi-\theta_{3})},1)^{T}$, with the basis $\Psi(x_{3})=(\psi_{6}(x_{3}),\psi_{1}(x_{3}))^{T}$.

Due to the continuity of the bound state wave function, the boundary conditions are:

[TABLE]

where $\alpha_{i=1,2,3}$ are phase factors acquired across the junction as mentioned in main text. We also have $\psi_{1}(x_{3}\rightarrow 0^{+})=\psi_{2}(x_{1}\rightarrow 0^{+}),~{}\psi_{3}(x_{1}\rightarrow 0^{+})=\psi_{4}(x_{2}\rightarrow 0^{+}),~{}\psi_{5}(x_{2}\rightarrow 0^{+})=\psi_{6}(x_{3}\rightarrow 0^{+})$ since they are the edge states from the same Chern insulator layer. With these we have:

[TABLE]

which is equivalent to

[TABLE]

This is in accordance with Eq.[14] in main text. Thus the poles of the scattering matrix do correspond to the localized states at the junction. Similar results also apply for a vertex of octahedral or icosahedral topological fullerenes, as shown in Sec. IV.1.

2.2.2 Y-junction scattering matrix

In this section we provide more details about how to derive the full scattering matrix Eq.[13] in Sec. IV.1. Similarly to the two terminal junction, for the scattering states of Y-junction, we denote $\eta=v/m$. We set $k\eta=\sinh\varphi>0$. Then for $x_{i}>0$, under the basis \big{(}\psi_{\alpha}(x_{i}),\psi_{\beta}(x_{i})\big{)}^{T}, for each individual leg, from Eq.[1] we can derive normalized wave function of incoming and outgoing modes as: $\psi_{\text{in}}^{T}(x_{i})=(e^{-\varphi-i\theta_{i}},1)^{T}/\sqrt{1+e^{-2\varphi}},\quad\psi_{\text{out}}^{T}(x_{i})=(e^{\varphi-i\theta_{i}},1)^{T}/\sqrt{1+e^{+2\varphi}}$, with corresponding energy $E/m=+\cosh\varphi$. As mentioned in main text, we first inject the mode along negative $x_{1}$ direction. The wave function on leg $x_{1}$ is given by $\Psi_{1}=e^{-ikx_{1}}\psi_{1,\text{in}}+re^{ikx_{1}}\psi_{1,\text{out}}$. Meanwhile, the wave function on leg $x_{2}$ is given by $\Psi_{2}=t_{12}e^{ikx_{2}}\psi_{2,\text{out}}$, and the wave function on leg $x_{3}$ is given by $\Psi_{3}=t_{13}e^{ikx_{3}}\psi_{3,\text{out}}$. We can expand the wave function around the intersection as:

[TABLE]

Note that due to the continuity of edge state wave function for each individual layer, we have $\psi_{\alpha}(x_{1}\rightarrow 0^{+})=\psi_{\beta}(x_{3}\rightarrow 0^{+})$, $\psi_{\beta}(x_{1}\rightarrow 0^{+})=\psi_{\alpha}(x_{2}\rightarrow 0^{+})$, and $\psi_{\beta}(x_{2}\rightarrow 0^{+})=\psi_{\alpha}(x_{3}\rightarrow 0^{+})$. Following the assumption we made in Sec. IV, during the scattering process, the amplitude of the chiral edge states from same triangular Chern insulator is conserved, but they may acquire an additional phase factor $\alpha_{i}$ when by passing the junction. By matching the coeffients of $\psi_{\alpha(\beta),i}$ we have:

[TABLE]

From the above equation we derive that:

[TABLE]

One can check that the scattering is unitary:

[TABLE]

To derive the full scattering matrix, we can further inject the mode along negative $x_{2}$ ($x_{3}$) direction. By following the similar procedure for injecting along negative $x_{1}$ direction, we have:

[TABLE]

[TABLE]

Finally, we derive scattering matrix as:

[TABLE]

where $kv/m=\sinh\varphi,~{}\Lambda=\sum_{i}(\theta_{i}+\alpha_{i})$. The poles of $S$ denotes the existence of bound state with energy:

[TABLE]

which is the Eq.[14] in main text. It is easy to check that the Scattering matrix here is unitary, i.e., $S^{\dagger}S={\mathbf{1}}$. Eq.[14] can be generalized to $l$-leg junction: ${E}/{m}=\cosh\varphi=\cos[{(\Lambda+2n\pi)}/{l}],~{}n\in\mathbb{Z},~{}N_{s}=-{\sum_{i}\theta_{i}}/{2\pi}$. For latter convenience we let $2n\pi$ be absorbed into $\sum_{i}\alpha_{i}$.

3.2.3 Comparison with numerical results from exact diagonalization of tight-binding Hamiltonian

The bound-state energy Eq.[15] is depending on $\Lambda=\sum_{i}(\theta_{i}+\alpha_{i})$. As we showed in main text, $\sum_{i}\alpha_{i}=A\sum_{i}\theta_{i}+B$, substitute these into Eq.[15] we have:

[TABLE]

where $l=f-6$ is the number of legs for a disclination with Frank index $f$. In order to figure out the value of $A,B$ and derive the full response function as Eq.[17], in principle we need two data points (the bound state energy at two different flux value $\Phi$) from the exact diagonalizing tight-binding Hamiltonian Eq.[9]. In fact, we do take two data points directly for $m\approx t_{0}$ and get Eq.[17] in the main text. However, note that for Frank index $f=3$ ($f=1$), although the response of bound state energy with respect to external flux for different $m$ are different, adding an external flux $\Phi=-\pi/2$ ($\Phi=+\pi/2$) can restore the particle hole symmetry, and move the bound state energy to zero. Thus we can define $\Phi_{0}=|{B}/({1+A})|=(\sum_{i}\theta_{i})~{}\text{mod}~{}\pi$, which is fixed for given $f$. Note that $\sum_{i}\theta_{i}$ is the total phase mismatch at the junction. With these Eq.[A23] can be reduced to:

[TABLE]

The plus or minus sign here depends on whether the local Chern vector is align with or opposite to the direction of external flux. Now we only need one date point (for example, the energy of bound state in the absence of external flux) from exact diagonalization to get the value $A$ in Eq.[A24] and reproduce Eq.[17] directly. For $m=3\sqrt{3}t_{1}\approx t_{0}=1$, we derive $A$ first and Eq.[A24] is then simplified as:

[TABLE]

which is Eq.[17] in the presence of external flux $\Phi$.

We further compare the results from Eq.[17] with full numerical results derived from exactly diagonalizing the tight-binding Hamiltonian Eq.[16], as shown in Fig.[A2]. We plot the bound state energy with external flux ($\Phi$) under different Haldane mass $m=3\sqrt{3}t_{1}$ and different Frank index $f$. The direct fittings of numerical results do take the form of Eq.[A24], as shown in the left-bottom of each sub-figure.

4.2.4 Comparison to numerical results from continuous model

The bound state energy with respect to external flux from continuous model for conical singularities Rüegg and Lin (2013) is given by:

[TABLE]

where

[TABLE]

Here $m$ is the Haldane mass, $\rho$ stands for the radius of the hole in disclination, $E$ is the bound state energy, $j$ is half integer, $f$ is Frank index and stands for the number of $\pi/3$ wedges removed, $\Phi$ denotes the external flux, and $K(\kappa\rho)$ is modified Bessel functions of the second kind. We have set the positive direction of external flux opposite to local Chern vector. In practice, in order to derive full energy-flux relation for given $m$, one may need (at least) one data point (bound state energy at given $\Phi$) from exact diagonalizing Eq.[9] to get the value of effective radius $\rho$. After that we can derive the bound state energy with external flux from (numerically) solving Eq.[A26].

We have shown that our analytic results in Eq.[A24] fit quite well with the numerical results from diagonalizing tight-binding model in previous subsection. Our method also give the proper results from solving Eq.[A26] directly, as shown in Fig.[A3.(a)]. From here we know that, the phase shift $\sum_{i}\alpha_{i}$ should be a function of Haldane mass $m$ and effective radius $\rho$.

0.3 Boundary Hamiltonian for arbitrary edge

In this section we derive the effective edge Hamiltonian for an arbitrary edge. Note that the in gap state wave function distribution for our Tetrahedral type TI (Eq.[10]) is different from that of Quadrupole insulator, see in Fig.[A3.(b,c)]. We further show that our Tetrahedral type 2d SOTI can hold fractional charge at the corner of rectangular boundaries, regardless of the orientation of the rectangle.

The Bloch Hamiltonian for our Tetrahedral type TI is Eq.[10], as shown in main text. In the absence of inter-layer coupling, i.e. $\lambda=0$, the system can be viewed as index spin hall effect:

[TABLE]

Around $(k_{x}=0,k_{y}=0)$, the low energy version for Hamiltonian Eq.[A28] is given by:

[TABLE]

where $\tilde{m}_{z}=m_{z}-2t_{x}-2t_{y}$. For simplicity we assume $t_{x}=t_{y}=t_{0}$.

In order to figure out the edge states at the cut along $\vec{e}_{2}=\cos\theta\hat{e}_{x}+\sin\theta\hat{e}_{y}$ direction (see in Fig.[A4].(a)), we define a new set of basis in both spatial and momentum spaces:

[TABLE]

Substituting Eq.[A30] into Eq.[A29], the Low energy Hamiltonian can be written in the form of $k_{1,2}$:

[TABLE]

Consider the model Hamiltonian Eq.[A31] defined on the half-space $x_{1}>0$ in the $x_{1}-x_{2}$ plane. We replace $k_{1}\rightarrow-i\partial_{x_{1}},k_{2}\rightarrow 0$, and neglect the higher order terms in Eq.[A31]:

[TABLE]

By using the ansatz $\psi_{0}=e^{\eta x_{1}}\phi$, we can find a pair of counter-propagating edge states:

[TABLE]

where $N_{\uparrow(\downarrow)}$ is the normalization constant. This procedure Qi and Zhang (2011) leads to a $2\times 2$ effective Hamiltonian defined by $H^{\alpha,\beta}_{\text{edge}}(k_{2})=\bra{\Psi_{\alpha}}h(\vec{k})\ket{\Psi_{\beta}}$, to the leading order in $k_{2}$, we arrive at the effective Hamiltonian for helical edge states:

[TABLE]

Similarly, the inter-layer coupling $\lambda\tau_{1}\sigma_{1}$ (or $\lambda\tau_{2}\sigma_{2}$), under the basis $\Psi_{\alpha,\beta}$, gives birth to an additional term $-\lambda\sin\theta\sigma_{y}$ (or $-\lambda\cos\theta\sigma_{x}$). In summary, under the basis $\Psi_{\alpha,\beta}$, the total effective edge Hamiltonian is given by:

[TABLE]

We can define the effective mass term as:

[TABLE]

This related the effective mass term of $i$-th edge to its orientation $\hat{e}_{i}=\cos\theta_{i}\vec{e}_{x}+\sin\theta_{i}\vec{e}_{y}$. According to our previous results, the kink of effective mass term at the corner can give birth to corner localized charge. The value of the charge (edge soliton) $N_{s}$ is:

[TABLE]

For any rectangular boundary, $\delta\theta=\pi/2$ since two adjoint edges are perpendicular to each other. Thus the corner localized fractional charge should be $-1/4e$, regardless the orientation of rectangle. We have confirmed this by exact diagonalizing the tight-binding Hamiltonian, as shown in Fig.[A4.(b,c)]. This result can be generalized to the corner state with arbitrary fractional charge by tuning the angle $\theta$ between two adjoint edges.

0.4 Fractional charge for edge soliton

In the absence of particle hole symmetry, the domain wall state for a SSH chain can hold bound state with non-zero energy and fractional charge aside from $-1/2e$ Jackiw and Semenoff (1983); Goldstone and Wilczek (1981). In this section, we summarized and slightly modified their previous works Jackiw and Semenoff (1983) and derive the similar results for edge solitons. This is in accordance with the results Eq.[6] from bosonization in Sec. II.

Suppose we have a one-dimensional Dirac Hamiltonian in the external field $\varphi$:

[TABLE]

For simplicity we assume that $\epsilon>0$. Up to a global normalization constant and a unitary transformation this Hamiltonian can be connected to Hamiltonian Eq.[A35]. In the absence of $\epsilon\sigma_{x}$, the Hamiltonian respects the charge conjugation symmetry and can hold zero mode when $\varphi(x)$ has a kink. The presence of $\epsilon\sigma_{x}$ breaks the charge conjugation symmetry of the system. We will see later on that this Hamiltonian can hold bound state with nonzero energy and fractional charge.

In the vacuum where the system does not hold a soliton, $\varphi=\varphi_{0}=\text{const}$. We denote $\varphi_{0}=\mu$ for simplicity. In the presence of a soliton, $\varphi(x)=\varphi_{s}(x)$, and in principle the $\varphi_{s}(x)$ should have a kink. In order to compute the charge, we need to derive the eigenstates of this two situation. The Schrodinger equation for these two scenarios can be written as:

[TABLE]

where $\psi_{E}^{0}$ stands for the normal state without solitons, $\psi^{s}_{E}$ stand for the situation in the presence of soliton.

The charge density at level $E$ is $\rho_{E}(x)=\psi_{E}^{\dagger}(x)\times\psi_{E}(x)$, and the physical charge density is got by integrating $\rho_{E}$ over all negative $E$, since the negative energy levels are filled in the half-filling:

[TABLE]

Finally the soliton charge is obtained by integrating the charge density in the soliton field over all $x$, but to avoid an infinity, we must subtract a similar integral of the charge density when no soliton is present:

[TABLE]

We can calculate the exact value of $Q$ even if we do not know the exact form of $\varphi(x)$. All we need to know about $\varphi_{s}$ is that it interpolates between opposite “vacuum” values as $x$ passes from $-\infty$ to $+\infty$:

[TABLE]

We now study the eigenstates of Eq.[A39]. The vacuum problem is trivial: the wave functions are plane waves $\propto e^{ikx}$ and the spectrum is continuous $E^{0}=\pm(k^{2}+\mu^{2}+\epsilon^{2})^{1/2}$.

In the presence of soliton, we first assume that the wave-function of the eigenstate is $(u,v)^{T}$. Thus we have:

[TABLE]

which can be simplified as:

[TABLE]

In order to solve these two equations, we first add up two equations:

[TABLE]

and then subtract the second equation from the first one, such that we have:

[TABLE]

We define the new parameters:

[TABLE]

then we can rewrite the result as:

[TABLE]

From the second line of Eq.[A48] we know that

[TABLE]

Substitute this into the first line of Eq.[A48], we have:

[TABLE]

From Eq.[A48] and Eq.[A50] we can figure out a possible solution:

[TABLE]

corresponds to the energy $E=\epsilon$. Note that the $U$ is localized at the kink $x=0$ due to the form of $\varphi_{s}(x)$.

To calculate the particle density, we still need to know the eigenstate for all negative energy solutions. We assume that $U\propto e^{ikx}$ and $\varphi\approx\pm\mu$ at large $x$ limit, thus we have the normalized factor:

[TABLE]

from which we can figure out the normalized wave function for the negative energy:

[TABLE]

This gives the wave function in originally basis:

[TABLE]

where

[TABLE]

The wave function $\psi_{k}$ satisfies:

[TABLE]

The Charge-density at negative $E$ is given by:

[TABLE]

where the validity of second line comes from Eq.[A50].

The soliton charge is the integral over all $x$ and $k$ above evaluated with $\varphi=\varphi_{s}$, minus a similar integral in the vacuum; but in the vacuum, $|U_{k}|^{2}\equiv\varphi(x)=\mu$, such that the last two term in Eq.[A57] vanished. Thus we have the soliton charge:

[TABLE]

The double integral can be evaluated by completeness: The $U^{0}_{k}$ represent all the Schrodinger modes in the vacuum, while the $U_{k}^{s}$ are one short of being complete in the soliton sector, since the normalized bound state is not among them. Hence the first term contributes $-1$ to $Q$. To evaluate the second term in Eq.[A58], let us consider the wave function in the presence of a soliton when $x=\pm\infty$. These may be given in terms of transmission ($T$) and reflection coefficients ($R$):

[TABLE]

Thus, upon dropping oscillatory terms, we are left with the soliton charge:

[TABLE]

where the plus sign between the contributions at $x=+\infty$ and at $x=-\infty$ arises because of sign reversal in $\varphi_{s}(x)$. Unitarity, $|T|^{2}+|R|^{2}=1$, permits a final evaluation:

[TABLE]

Note that, if we denote $\epsilon=m\cos\theta$ and $\varphi_{s}(\pm\infty)=\pm\mu=m\sin(\mp\theta)$, Eq.[A61] is reduced to:

[TABLE]

which is in accordance with Eq.[6] derived from bosonization of the helical Luttinger liquid.