Boundary Green functions of topological insulators and superconductors

TL;DR

This paper introduces a recursive boundary Green function method to analyze topological insulators and superconductors, revealing their phase diagrams, edge mode localization, and topological indices across various models.

Contribution

It develops a novel recursive approach to boundary Green functions that captures topological phases and boundary physics in a unified framework.

Findings

Provides a recursive method for boundary Green functions

Identifies fixed points corresponding to topological phases

Applies method to models like SSH, Kitaev chain, and Chern insulator

Abstract

Topological insulators and superconductors are characterized by their gapless boundary modes. In this paper, we develop a recursive approach to the boundary Green function which encodes this nontrivial boundary physics. Our approach describes the various topologically trivial and nontrivial phases as fixed points of a recursion and provides direct access to the phase diagram, the localization properties of the edge modes, as well as topological indices. We illustrate our approach in the context of various familiar models such as the Su-Schrieffer-Heeger model, the Kitaev chain, and a model for a Chern insulator. We also show that the method provides an intuitive approach to understand recently introduced topological phases which exhibit gapless corner states.

Tables

Table 1. Table 1: Periodic table of topological insulators and superconductors in d = 0 , … , 7 𝑑 0 … 7 d=0,\dots,7 dimensions. The first column denotes the ten symmetry classes of fermionic Hamiltonians, charactered by the absence ( 0 0 ) or presence ( ± plus-or-minus \pm or 1 1 1 ) of time-reversal ( 𝒯 ^ ^ 𝒯 \hat{\mathcal{T}} ), particle-hole ( 𝒞 ^ ^ 𝒞 \hat{\mathcal{C}} ), and chiral symmetries ( 𝒮 ^ ^ 𝒮 \hat{\mathcal{S}} ). Time reversal and particle hole symmetry exist in two types denoted by ± plus-or-minus \pm . Classes which support only trivial phases are denoted by ”-”, while classes with nontrivial topological classifications are indicated by the type of topological invariant ( ℤ , 2 ​ ℤ , ℤ 2 ℤ 2 ℤ subscript ℤ 2 \mathbb{Z},2\mathbb{Z},\mathbb{Z}_{2} ).

Class	$\widehat{𝒯}$	$\widehat{𝒞}$	$\widehat{𝒮}$	0	1	2	3	4	5	6	7
A	0	0	0	$\mathbb{Z}$	-	$\mathbb{Z}$	-	$\mathbb{Z}$	-	$\mathbb{Z}$	-
AIII	0	0	1	-	$\mathbb{Z}$	-	$\mathbb{Z}$	-	$\mathbb{Z}$	-	$\mathbb{Z}$
AI	+	0	0	$\mathbb{Z}$	-	-	-	$2\mathbb{Z}$	-	${\mathbb{Z}}_2$	${\mathbb{Z}}_2$
BDI	+	+	1	${\mathbb{Z}}_2$	$\mathbb{Z}$	-	-	-	$2\mathbb{Z}$	-	${\mathbb{Z}}_2$
D	0	+	0	${\mathbb{Z}}_2$	${\mathbb{Z}}_2$	$\mathbb{Z}$	-	-	-	$2\mathbb{Z}$	-
DIII	-	+	1	-	${\mathbb{Z}}_2$	${\mathbb{Z}}_2$	$\mathbb{Z}$	-	-	-	$2\mathbb{Z}$
AII	-	0	0	$2\mathbb{Z}$	-	${\mathbb{Z}}_2$	${\mathbb{Z}}_2$	$\mathbb{Z}$	-	-	-
CII	-	-	1	-	$2\mathbb{Z}$	-	${\mathbb{Z}}_2$	${\mathbb{Z}}_2$	$\mathbb{Z}$	-	-
C	0	-	0	-	-	$2\mathbb{Z}$	-	${\mathbb{Z}}_2$	${\mathbb{Z}}_2$	$\mathbb{Z}$	-
CI	+	-	1	-	-	-	$2\mathbb{Z}$	-	${\mathbb{Z}}_2$	${\mathbb{Z}}_2$	$\mathbb{Z}$

Equations

$$H({\mathbf{k}}_{\perp })=\sum _{n=1}^N\left({\psi }_n^{†}{h}_n{\psi }_n+{\psi }_{n+1}^{†}{V}_n{\psi }_n+{\psi }_n^{†}{V}_n^{†}{\psi }_{n+1}\right).$$

$$\mathbf{H}=\left(\begin{array}[]{ccccc}h_{1}&V_{1}^{\dagger}\\ V_{1}&h_{2}&V_{2}^{\dagger}\\ &\ddots&\ddots&\ddots\\ &&V_{N-2}&h_{N-1}&V_{N-1}^{\dagger}\\ &&&V_{N-1}&h_{N}\end{array}\right).$$

$$(\omega \mathbb{I}-\mathbf{H})\mathbf{G}(\omega )=\mathbb{I}.$$

$$\left(g_N^{-1}-V_{N-1}{G}_{N-1}{V}_{N-1}^{†}\right)G_N=\mathbb{I}.$$

$$\left(g^{-1}-VG{V}^{†}\right)G=\mathbb{I}.$$

$$H=\sum _n\left(-t_1{c}_{2n}^{†}{c}_{2n-1}-t_2{c}_{2n+1}^{†}{c}_{2n}+h.c.\right),$$

$$V_N=\{\begin{array}{ll}{t}_1& \mathrm{odd}N\\ t_2& \mathrm{even}N\end{array}$$

$$G_{2N}=[\omega -t_1^2(\omega -t_2^2{G}_{2N-2}{)}^{-1}{]}^{-1}.$$

$$G_{2N}-G_{2N-2}=\beta (G_{2N-2})$$

$$\beta (x)=[\omega -t_1^2(\omega -t_2^{2}x{)}^{-1}{]}^{-1}-x.$$

$$G_{\mathrm{triv}}$$

$$G_{\mathrm{top}}$$

$${\beta }^{\prime }(G_{\mathrm{top}})$$

$${\beta }^{\prime }(G_{\mathrm{triv}})$$

$$\xi =\frac{4[\max (t_1,t_2){]}^2}{|t_1^2-t_2^2|},$$

$$r_N=\mathbb{I}-2iW\frac{1}{G_N^{-1}+i{W}^{†}W}{W}^{†}.$$

$$r_N=\frac{\mathbb{I}-iW{G}_N{W}^{†}}{\mathbb{I}+iW{G}_N{W}^{†}}.$$

$$H_{\mathrm{bound}}=-{\left[G(\omega =0)\right]}^{-1}.$$

$$g_n^{-1}(\omega )\psi (n)-V_{n-1}\psi (n-1)-V_n^{†}\psi (n+1)=0.$$

$$M_{n}(\omega)=\left(\begin{array}[]{cc}(V_{n}^{\dagger})^{-1}g_{n}^{-1}(\omega)&-(V_{n}^{\dagger})^{-1}V_{n-1}\\ \mathbb{I}&0\end{array}\right).$$

$${\mathcal{M}}_N=M_N{M}_{N-1}\dots M_1,$$

$$\left(\begin{array}[]{c}0\\ \psi(N)\end{array}\right)=\mathcal{M}_{N}(\omega_{0})\left(\begin{array}[]{c}\psi(1)\\ 0\end{array}\right).$$

$$\det {\mathcal{M}}_{N,11}({\omega }_0)=0,$$

$$\ker {\mathcal{M}}_{N,11}({\omega }_0)=\mathrm{span}\{{\psi }_1(1),\dots ,{\psi }_k(1)\}.$$

$$(V_N^{†}{)}^{-1}{g}_N^{-1}{X}_N-(V_N^{†}{)}^{-1}{V}_{N-1}{X}_{N-1}{X}_N=\mathbb{I},$$

$$\left(\begin{array}[]{c}\mathbb{I}\\ X_{N}\end{array}\right)=M_{N}\left(\begin{array}[]{c}\mathbb{I}\\ X_{N-1}\end{array}\right)X_{N},$$

$$\displaystyle\left(\begin{array}[]{c}\mathbb{I}\\ X_{N}\end{array}\right)=\mathcal{M}_{N}\left(\begin{array}[]{c}\mathbb{I}\\ 0\end{array}\right)X_{1}\dots X_{N}.$$

$$X_1\dots X_N={\mathcal{M}}_{N,11}^{-1}.$$

$$X_N={\mathcal{M}}_{N,21}{X}_1\dots X_N={\mathcal{M}}_{N,21}{\mathcal{M}}_{N,11}^{-1}.$$

$$G_N={\mathcal{M}}_{N,21}{\mathcal{M}}_{N,11}^{-1}(V_N^{†}{)}^{-1}$$

$$M_n^{†}{\Omega }_n{M}_n={\Omega }_{n-1}$$

Full text

I Introduction

There is currently considerable interest in topological phases of matter and topological quantum phase transitions. The best studied and simplest variety of topological phases are the independent-fermion phases Hasan and Kane (2010); Qi and Zhang (2011); Bernevig and Hughes (2013) collected and classified in the periodic table and its extensions Schnyder et al. (2008); Kitaev (2009); Ryu et al. (2010); Teo and Kane (2010); Essin and Gurarie (2011); Fulga et al. (2012); Chiu et al. (2016). These topological phases can be characterized in terms of bulk topological indices and, by the bulk-boundary correspondence, exhibit gapless modes localized at the sample boundaries.

Here, we present a simple and intuitive recursive approach which provides access to these topological phases and their transitions. Our approach focuses on the boundary Green function Sancho et al. (1985); Umerski (1997) which encodes the presence or absence of gapless surface modes in an immediate manner. Focusing on tight-binding models, we imagine that the bulk system is constructed by iteratively adding boundary sites (in one dimension) or boundary layers (in higher dimensions). We can then exploit that the boundary Green function can be computed recursively and, by virtue of the bulk gap, becomes invariant under this recursion in the thermodynamic limit. Thus, the boundary Green function can be obtained as fixed points of the recursion.

For a typical model which exhibits topologically trivial and nontrivial phases, the recursion will generate several fixed-point boundary Green functions. The physical boundary Green function and hence the phase diagram can be identified from a stability analysis. The various phases correspond to the stable fixed points of the recursion and the approach to the fixed point contains information about the spatial localization of the boundary mode away from the boundary as well as the bulk gap. This has some similarities with a real space renormalization group approach.

The fixed-point boundary Green function provides access to various observables of interest such as the tunneling density of states. Most importantly, it immediately contains information on the topological invariants by virtue of the fact that it encodes the full boundary spectrum. We can relate our approach more systematically to the general classifications of free-fermion topological phases by establishing a connection to their classification in terms of scattering matrices Fulga et al. (2012). This can be done by expressing the reflection matrix describing reflection from a gapped phase in terms of its boundary Green function. In this way, we immediately obtain explicit expressions for the topological invariants in terms of the boundary Green function.

Our approach can be applied to a wide variety of models. Here, we illustrate the method on a number of popular models. For simple cases such as the Su-Schrieffer-Heeger model Su et al. (1979), the recursion is conveniently implemented in a direct manner. To make the approach more systematic and to simplify its implementation, we find it helpful to introduce an approach based on the transfer matrix Pichard and Sarma (1981); Lee and Joannopoulos (1981); Hatsugai (1993a, b). This is illustrated in one dimension for the Kitaev chain Kitaev (2001) and in two dimensions for a model of a Chern insulator Bernevig and Hughes (2013).

Very recently, it was pointed out Benalcazar et al. (2016) that some 2D and 3D models exhibit topologically protected excitations at corners. These models have boundary Green functions which are regular, reflecting the fact that the boundary is gapped. Nevertheless, it is natural to interpret the inverse of the boundary Green function at the Fermi energy as an effective boundary Hamiltonian which can by itself be topologically trivial or nontrivial. We show for a model discussed in Ref. Benalcazar et al. (2016) that this effective boundary Hamiltonian is indeed topologically nontrivial when the model exhibits topological corner states. Our boundary Green function method may thus open the path to a more systematic analysis of such higher-order topological phases.

The paper is organized as follows. In Sec. II, we introduce the principal idea of this work and illustrate it with the Su-Schrieffer-Heeger model. We also show how to relate the boundary Green function to the reflection matrix which connects our approach with the systematic classification of free-fermion systems. In Sec. III, we show how the transfer matrix can be employed for computing the fixed-point boundary Green function. This is used in Sec. IV to compute the boundary Green function for the Kitaev chain and a model of a Chern insulator. For both models, we show how the topological invariants can be extracted from the fixed-point boundary Green function. A 2D model with topological corner states is investigated in Sec. V. We conclude in Sec. VI. Some aspects are relegated to appendices.

II Basic considerations

II.1 Boundary Green function

Consider a $d$-dimensional tight-binding Hamiltonian, which is translationally invariant along $d-1$ directions. Assuming periodic boundary conditions in these directions, we can define a $(d-1)$-dimensional crystal momentum $\mathbf{k}_{\perp}$ to label the Hamiltonian. In the remaining dimension, we apply open boundary conditions and retain a real-space representation. By choosing a sufficiently large unit cell, a Hamiltonian with a finite hopping range can always be brought into a form which couples only neighboring unit cells. Thus, without loss of generality, we can consider a family of quasi-1D Hamiltonians which are labelled by the transverse momentum $\mathbf{k}_{\perp}$ and take the form

[TABLE]

Here, $n=1,2\ldots N$ labels the unit cells along the direction with open boundary conditions, and $\psi_{n}$ ($\psi^{\dagger}_{n}$) is a column (row) vector of electron annihilation (creation) operators in the $n$th unit cell. These vectors have $M$ entries, reflecting the dimension of the local Hilbert space associated with each unit cell and a given $\mathbf{k}_{\perp}$. The intra- and inter-unit-cell couplings are described by $M\times M$ matrices $h_{n}$ and $V_{n}$, respectively. Note that we suppressed the $\mathbf{k}_{\perp}$ dependence of $\psi_{n}$, $h_{n}$, and $V_{n}$ for notational simplicity. The structure of this quasi-1D tight-binding Hamiltonan is illustrated in Fig. 1.

More explicitly, the corresponding first-quantized Hamiltonian can be written as an $NM\times NM$ block tridiagonal Hamiltonian matrix $\mathbf{H}(\mathbf{k}_{\perp})$ which takes the form

[TABLE]

Both the spectrum and the eigenfunctions of the system are encoded in the Green function (or resolvent operator) $\mathbf{G}(\omega)$ of this Hamiltonian, which is defined through

[TABLE]

Here, $\mathbb{I}$ is the identity matrix of dimension $NM$.

The low-energy excitations of topological insulators and superconductors are localized near their boundary. This motivates us to consider the boundary Green functions $G_{1}=\mathbf{G}_{11}$ and $G_{N}=\mathbf{G}_{NN}$, which correspond to the $1,1$ and $N,N$ blocks of the full Green function matrix $\mathbf{G}$. Thus, these boundary Green functions are $M\times M$ matrices.

By virtue of the tridiagonal structure of the first-quantized Hamiltonian (2), we can compute the boundary Green function recursively, extending the system unit cell by unit cell in the direction with open boundary conditions. In this way, we can relate the boundary Green function $G_{N}$ of a chain with $N$ unit cells to the boundary Green function $G_{N-1}$ of a chain with $N-1$ unit cells. Indeed, the two boundary Green functions are simply related by the Dyson equation

[TABLE]

Here, $g_{N}^{-1}(\omega)=\mathbb{I}\omega-h_{N}$ is the bare Green function of the $N$th unit cell. Similar recursions were used in the context of topological phases in Refs. Dai et al. (2008); Kim et al. (2015).

When considering insulators and superconductors with a gapped bulk, be they topological or not, we expect that the boundary Green function becomes independent of the number of layers in the limit of large $N$. Hence, the boundary Green function should approach a fixed point $G$ which satisfies the fixed-point equation

[TABLE]

Here, we have assumed that the system is uniform along the chain, so that we can drop the indices of $g$ and $V$. This closed equation for the boundary Green function (a matrix quadratic equation) is the central starting point of our considerations.

II.2 The Su-Schrieffer-Heeger (SSH) model: An illustrative example

It should be evident that the boundary Green function encodes much information of interest in the context of topological insulators and superconductors. Most importantly, it knows about gapless end, edge, or surface states, and can thus be used to derive their properties (including their extension into the bulk) as well as the topological phase diagram of the model. It is useful to illustrate this statement in the context of the Su-Schrieffer-Heeger model. This model is well known to exhibit topologically distinct phases and is sufficiently simple so that we can implement the recursive approach in a rather straight-forward manner, without relying on the more powerful methods developed in the next sections.

The Su-Schrieffer-Heeger (SSH) model describes a 1D chain of spinless fermions with one orbital per site and alternating hopping strengths between neighbouring sites. For a chain with $2N$ sites and open boundary conditions, the Hamiltonian takes the form

[TABLE]

where $c_{j}$ ($c_{j}^{\dagger}$) annihilates (creates) a fermion at site $j$, and $t_{1}$ and $t_{2}$ are the hopping amplitudes (taken to be real for simplicity). While the model is gapless for $\left|t_{1}\right|=\left|t_{2}\right|$, the dimerization opens a gap around $E=0$ when the hopping amplitudes differ in magnitude. The gapless point separates a topological phase for $\left|t_{1}/t_{2}\right|<1$ from a trivial phase for $\left|t_{1}/t_{2}\right|>1$. In the topological phase, there is a midgap state localized at each end of the chain in the thermodynamic limit.

For the SSH model, we have $g_{N}^{-1}=\omega$ and

[TABLE]

from Eq. (6). In view of the dimerization of the hopping amplitudes, it is convenient to iterate the recursion Eq. (4) and to obtain a recursion for chains with an even number of sites,

[TABLE]

To analyze this recursion, we write it as

[TABLE]

with the $\beta$-function

[TABLE]

The fixed-point boundary Green function $G$ follows from the zeroes of the $\beta$-function. Solving the corresponding quadratic equation, we find the two solutions

[TABLE]

where we have restricted attention to the limit of small $\omega$.

These two fixed-point boundary Green functions correspond to the topological and nontopological phases of the model. In the topological phase, the midgap state makes the boundary Green function singular for $\omega\to 0$. In contrast, there is no such singularity in the trivial phase where the end remains gapped and the boundary Green function is regular. These statements are an immediate consequence of the Lehmann representation for the Green function.

The phase diagram of the model emerges when studying the stability of the fixed-point Green functions under the recursion (8). A fixed point is stable when $\beta^{\prime}(x)<0$. A simple calculation yields

[TABLE]

for $\omega=0$. This implies that the singular boundary Green function $G_{\rm top}$ corresponding to the topological phase is stable for $\left|t_{1}/t_{2}\right|<1$, while the regular boundary Green function $G_{\rm triv}$ corresponding to the trivial phase is stable for $\left|t_{1}/t_{2}\right|>1$.

The fixed-point analysis also provides information about the gap, the correlation length of the model, and their behavior near the topological critical point $t_{1}=t_{2}$. The correlation length of the phase with fixed-point Green function $G$ is given by $\xi=4/|\beta^{\prime}(G)|$. (Notice that we define $\xi$ as the decay length of a wavefunction; the diagonal elements of the Green function decay on half of this scale.) Thus, we find the result

[TABLE]

which is written in a way that applies to both phases. The correlation length diverges at the topological critical point $|t_{1}|=|t_{2}|$. It is also interesting to note that this result is consistent with the quasiparticle weight $Z$ of the fixed-point Green function. In the topological phase, one expects $G_{\rm top}\simeq Z/\omega$ for small $\omega$ and $Z=4/\xi$ reflecting the fact that the midgap state is rapidly oscillating and decays exponentially into the bulk on the scale of the correlation length.

To extract the spectral gap $\Delta$, we note that it is related to the correlation length through $\Delta=\hbar v_{F}/\xi$ (assuming that the gap is small compared to the bandwidth, i.e., $|t_{1}-t_{2}|\ll t_{1},t_{2}$). Here, the Fermi wavelength should be taken at the center of the band with $t=t_{1}=t_{2}$, i.e., $v_{F}=2t/\hbar$. This yields $\Delta=|t_{1}-t_{2}|$ in agreement with a direct evaluation of the spectrum.

II.3 Topologial invariants

The boundary Green function contains direct information on the topological nature of the phases through the gapless end, edge, or surface states described by the boundary Green function. As for the SSH model above, this can be read off directly from the boundary Green function at energies $\omega$ below the bulk gap. In 1D, the topological phase is signaled by midgap end states and hence by a contribution to the Green function which is singular for $\omega\to 0$. In higher dimensions, there are edge or surface modes with a linear dispersion at low energies. We will discuss this in more detail below.

A systematic connection to the periodic table of topological phases can be established by relating the boundary Green function to the reflection matrix. Ref. Fulga et al. (2012) derived the periodic table based on the reflection matrix $r_{N}$ and provided explicit expressions for the topological invariants in terms of the $r_{N}$. With the relation between the reflection matrix and the boundary Green function, the topological invariants can alternatively be written in terms of the latter.

Start with a tight-binding model of the kind discussed above and attach a lead to one end with $M$ propagating channels for each ${\mathbf{k}}_{\perp}$. In the spirit of the tight-binding model, the channels are only coupled to the last site of the chain, with the matrix $W$ denoting the coupling between the lead channels and the adjacent site of the system. By the gapped nature of the system, the incoming wave is fully reflected. Based on the Mahaux-Weidenmüller formula, the corresponding scattering – or reflection – matrix $r_{N}$ can then be written as Mahaux and Weidenmüller (1969)

[TABLE]

Note that the Green function $G_{N}$ appearing in this expression is indeed the boundary Green function since the coupling is to the last site of the chain only. Some rearranging then leads to the relation Aleiner et al. (2002)

[TABLE]

We include an alternative and more detailed derivation of this relation in App. A. Moreover, we couch the scattering-matrix approach to the periodic table and the topological invariants in terms of the boundary Green function in App. B.

Before concluding this section, we briefly consider systems with a trivial bulk topological invariant and thus without a gapless boundary mode. In this case, the boundary Green function evaluated at $\omega=0$ is generally invertible. This can be used to define an effective boundary Hamiltonian

[TABLE]

Even if the original model does not have a nontrivial bulk topological invariant, this boundary Hamiltonian can still be topologically nontrivial, leading to a second order topological phase. We use this observation to extend our approach to describe higher-order topological phases in Sec. V.2.

III Boundary Green function and transfer matrix

In this section, we introduce a more efficient method to compute the fixed-point boundary Green function which comes in handy when discussing models which are more complicated than the SSH model. The method is based on the transfer matrix Pichard and Sarma (1981); Lee and Joannopoulos (1981) and effectively provides an explicit construction of the fixed-point boundary Green function.

We start by reviewing the transfer matrix based on the Schrödinger equation for the Hamiltonian in Eq. (1),

[TABLE]

Here, ${\Psi}(n)$ is an $M$-component wave function. The site index $n$ takes on values $n=1,\dots,N$, and we impose the open boundary conditions $\psi_{j}(0)=0=\psi_{j}(N+1)$.

As usual for tight-binding Hamiltonians, we can relate the two-component quantities $\Psi(n)=[\psi(n+1)^{T},\psi(n)^{T}]^{T}$ via $\Psi(n)=M_{n}(\omega)\Psi(n-1)$, where the matrix $M_{n}(\omega)$ takes the form

[TABLE]

Here and in the following, we assume $V_{n}$ to be invertible. Notice that $V_{0}$ and $V_{N}$ do not appear in the Hamiltonian and can be chosen as any invertible matrix in order to define $M_{0}$ and $M_{N}$. The matrix $M_{n}(\omega)$ is sometimes called a transfer matrix. Here, we reserve this terminology for the matrix

[TABLE]

connecting the two ends of the chain.

If $\omega_{0}$ is an eigenenergy of the Schrödinger equation (19), the transfer matrix $\mathcal{M}_{N}(\omega_{0})$ connects the eigenfunctions $\psi(n)$ at sites $n=1$ and $n=N$,

[TABLE]

If the state corresponding to the eigenenergy is $k$-fold degenerate, ${\Psi}_{p}(n)$ with $p=1,\dots,k$, then

[TABLE]

and

[TABLE]

Here, the subscript $11$ denotes the blocks under the partition of $\mathcal{M}_{N}$ into a $2\times 2$ block matrix. One immediate consequence of this is that $k\leq M$, namely the degeneracy is bounded by the number of internal degrees of freedom per unit cell.

We now show how the boundary Green function $G_{N}$ can be expressed in terms of the transfer matrix $\mathcal{M}_{N}$, starting directly from the recursion in Eq. (4). Defining $X_{n}=G_{n}V^{\dagger}_{n}$, we can rewrite the recursion as

[TABLE]

or equivalently

[TABLE]

where $M_{N}$ is defined in Eq. (20).

This relation can be applied iteratively to obtain

[TABLE]

The first row of this equation implies that

[TABLE]

Then, the second row yields

[TABLE]

Hence, we obtain the desired expression

[TABLE]

for the boundary Green function in terms of the transfer matrix. This formula was also derived in the broader context of block Toeplitz matrices in Ref. Boffi et al. (2014).

It is useful to exploit the symplectic structure of the transfer matrix. Notice that the matrix $M_{n}$ has the property

[TABLE]

in terms of the anti-hermitian matrix

[TABLE]

As a consequence, the transfer matrix $\mathcal{M}_{N}$ obeys the relation

[TABLE]

Remember that the choice of matrices $V_{N}$ and $V_{0}$ at the boundary is arbitrary as long as they are invertible. If we choose $V_{0}=V_{N}$, then the transfer matrix obtains a symplectic structure as it satisfies the relation

[TABLE]

As usual for symplectic matrices, this implies that the eigenvalues of $\mathcal{M}_{N}$ appear in pairs of $\lambda$ and $1/\lambda^{*}$. In particular, if $v$ is a right eigenvector of $\mathcal{M}_{N}$ with eigenvalue $\lambda$,

[TABLE]

then $(\Omega_{0}v)^{\dagger}$ is a left eigenvector of $\mathcal{M}_{N}$ with eigenvalue $1/\lambda^{*}$,

[TABLE]

This separates the $2M$ eigenvalues of the transfer matrix into $M$ eigenvalues $\lambda_{j}$ ($j=1,\ldots,M$) with $|\lambda_{j}|>1$, and another $M$ eigenvalues $1/\lambda_{j}^{*}$ with modulus smaller than unity. The Lyapunov exponents $1/\xi_{j}$ are then defined through

[TABLE]

and are strictly positive at the Fermi energy ($\omega=0$) for a gapped bulk Beenakker (1997). At the topological phase transition, the bulk gap closes and therefore at least one Lyapunov exponent vanishes for some $j$, corresponding to $\left|\lambda_{j}\right|=1$ or a diverging $\xi_{j}$. This is explicitly shown in App. C.

So far, we did not assume that the system is uniform along the direction with open boundary conditions. Let us now consider models with translational invariance of period $a$. (In the SSH model discussed above, this period was equal to 2 due to the dimerization of the hopping.) Considering a chain of length $Na$, the transfer matrix becomes

[TABLE]

where we defined the matrix

[TABLE]

associated with one period of length $a$.

In the translationally invariant case, the evaluation of the boundary Green function can be simplified by diagonalizing the matrix $T$. Again, $T$ satisfies a symplectic structure and its eigenvalues can be separated into $\Lambda=\operatorname{diag}(\lambda_{1},\dots,\lambda_{M})$ with $\left|\lambda_{j}\right|>1$, and a second set $(\Lambda^{-1})^{*}$. Thus, we can write

[TABLE]

We can use this to express $T$ and hence $\mathcal{M}_{aN}$, and insert this expression into the Eq. (34) for the boundary Green function. Upon taking the thermodynamic limit $N\to\infty$ and using that the system is gapped (see App. D for details), this yields the desired explicit expression

[TABLE]

for the fixed-point boundary Green function (see also Ref. Boffi et al. (2014)). While we have suppressed the argument $\omega$ in this expression, it is valid for any $\omega$ which falls inside the bulk gap.

For a homogeneous system in which the on-site potential $h_{n}$ and the hopping amplitude $V_{n}$ are independent of the site index $n$ along the chain, i.e., a system with $a=1$, the evaluation of the boundary Green functions can be simplified further. In this case, the columns of the matrix $U$ in Eq. (45) consist of the right eigenvectors of

[TABLE]

The simple structure of $T$ implies that the right eigenvector with eigenvalue $\lambda_{i}$ can be written as $(\lambda_{i}z_{i}^{T},z_{i}^{T})^{T}$, where the vector $z_{i}$ obeys

[TABLE]

Let us also denote the eigenvectors corresponding to eigenvalues $|\lambda_{i}|>1$ as $x_{i}$ and the corresponding eigenvectors with eigenvalues $1/\lambda_{i}^{*}$ as $y_{i}$. We can now identify $U_{21}=(x_{1},\ldots,x_{M})$ and $U_{11}=(\lambda_{1}x_{1},\ldots,\lambda_{M}x_{M})$. Then, Eq. (45) implies

[TABLE]

for the fixed-point boundary Green function. The essential steps in evaluating $G$ thus reduce to solving the nonlinear eigenvalue problem Eq. (47) and inserting the result into Eq. (48).

For later use, it is useful to collect expressions for all right and left eigenvectors of $T$ for $a=1$. We write the right eigenvectors of $T$ as $v_{i}=(x_{i}^{T},x_{i}^{T}/\lambda_{i})^{T}$ for $|\lambda_{i}|>1$ and as $v_{M+i}=(y_{i}^{T},y_{i}^{T}\lambda_{i}^{*})^{T}$ for the corresponding eigenvalues $1/\lambda_{i}^{*}$ with modulus less than unity. The symplectic nature of $T$ implies then that the corresponding left eigenvectors take the form $w_{i}^{T}=(-y_{i}^{\dagger}V^{\dagger}\lambda_{i},y_{i}^{\dagger}V)$ with eigenvalue $\lambda_{i}$ and $w^{T}_{M+i}=(-x_{i}^{\dagger}V^{\dagger}/\lambda_{i}^{*},x_{i}^{\dagger}V)$ with eigenvalue $1/\lambda_{i}^{*}$.

IV Application to 1D and 2D models

In this section, we apply the formalism developed in the previous section to two homogeneous tight-binding models with $a=1$. We consider the Kitaev chain as an example of a 1D model and a Chern insulator as an example of a 2D model.

IV.1 Perturbation theory

We are mainly interested in the fixed-point boundary Green function for small $\omega$ as this limit captures the topological properties and the effective low-energy excitations. This limit can be most directly treated by means of perturbation theory. Denoting $T_{0}=T(\omega=0)$, we have

[TABLE]

For small $\omega$, we can treat $\omega A$ in first-order perturbation theory (properly generalized to nonhermitian matrices, for which left and right eigenvectors are no longer identical).

If we denote by $w_{i}^{T}$ and $v_{i}$ the left and right eigenvectors of $M_{0}$ with eigenvalues $\lambda_{i}$ (not restricted in magnitude), we obtain the first-order corrections

[TABLE]

These expressions are evident generalizations of the corresponding results in quantum mechanics.

We can now use the explicit structure of the $v_{i}$ discussed at the end of the previous section as well as the simple form of $A$. This yields

[TABLE]

for the shift in the eigenvalues with $\left|\lambda_{i}\right|>1$ and

[TABLE]

with

[TABLE]

for the corresponding first-order corrections to the right eigenvectors. In the last expression, the sum is only over those indices $j$ with $\left|\lambda_{j}\right|>1$.

At small $\omega$, the fixed-point boundary Green function in Eq. (48) then becomes

[TABLE]

Notice that we retain the correction due to a nonzero $\omega$ only in the matrix inverse because this is where a singular $\omega$ dependence can originate from. Notice also that $\delta s_{i}$ does not enter this expression so that we refrain from giving an explicit expression.

IV.2 Kitaev chain

IV.2.1 Model and boundary Green function

As a first example, consider the Kitaev chain which is a lattice model for a spinless $p$-wave superconductor. The Hamiltonian for a chain of $N$ sites takes the form

[TABLE]

where $c_{j}$ ($c_{j}^{\dagger}$) is the fermion annihilation (creation) operator at site $j$. We assume below that both the hopping $t$ and the pairing $\Delta$ are real.

This model has trivial and topological phases, with the latter hosting Majorana zero modes at both ends of the chain. The trivial phase occurs for $\left|\mu\right|>2\left|t\right|$ when the chemical potential $\mu$ falls outside the normal-state band. The topological phase occurs for $\left|\mu\right|<2\left|t\right|$. The model supports a chiral symmetry $\mathcal{U}_{S}=\tau_{x}$, so that it falls into class BDI. Indeed, in addition to the topological phase transition lines at $\left|\mu\right|=2\left|t\right|$, the model also becomes gapless in the absence of pairing correlations, $\Delta=0$, and the topological phases with $\Delta>0$ and $\Delta<0$ have topological indices $\pm 1$.

To apply the approach developed in Sec. III, we introduce a Nambu spinor $\psi_{j}^{\dagger}=(c_{j}^{\dagger},c_{j})$. Then, we can identify

[TABLE]

where $\tau_{x,y,z}$ are Pauli matrices in Nambu space.

Applying the perturbative approach developed in Sec. IV.1, we first determine the eigenvalues and eigenvectors of $T$ for $\omega=0$. The eigenvalues follow from the condition

[TABLE]

which yields

[TABLE]

Denote the pair of roots for the upper sign as $\lambda_{1,2}$, the roots for the lower sign as $\lambda_{3,4}$, and fix the labeling by imposing $\lambda_{1}\lambda_{4}^{*}=1=\lambda_{2}\lambda_{3}^{*}$. The corresponding eigenvectors fulfilling Eq. (47) are

[TABLE]

To determine the two eigenvectors corresponding to eigenvalues with $|\lambda|>1$, we consider the quantity

[TABLE]

Thus, $\left|\lambda_{1}\right|$ and $\left|\lambda_{2}\right|$ are both larger than or both smaller than unity when $\left|2t\right|>\left|\mu\right|$.

Let us first consider the situation that $\left|\lambda_{1}\right|$ and $\left|\lambda_{2}\right|$ are both larger than unity. This occurs for $\left|2t\right|>\left|\mu\right|$ and $\Delta t<0$. Then, the matrix $(x_{1},x_{2})=1+\tau_{x}$ is noninvertible and the fixed-point boundary Green function $G$ becomes singular in the limit $\omega\to 0$. This is in agreement with expectations as the model exhibits a topological phase for these parameters.

Including the first-order corrections in $\omega$, we find

[TABLE]

where

[TABLE]

Thus, the fixed-point boundary Green function takes a form which implies the existence of a zero-energy end state. The rank of the numerator $1+\tau_{x}$ confirms that there is one zero-energy end state at each end, and the matrix structure imposes that the end-state wavefunction has particle-hole symmetry as it should for a Majorana bound state. The boundary Green function of the Kitaev chain was also derived in Ref. Zazunov et al. (2016) by a different method.

As for the SSH model, the quasiparticle weight of the fixed-point boundary Green function (62) encodes information on the localization length of the Majorana end state. In particular, one expects that the Majorana localization length diverges at the topological phase transition. For the parameters of Eq. (62), we have either $\left|\lambda_{1}\right|\to 1$ or $\left|\lambda_{2}\right|\to 1$, so that $a_{12}$ is dominated by one of the two terms and does indeed diverge.

Similarly, we find $\left|\lambda_{3,4}\right|>1$ for $\left|2t\right|>\left|\mu\right|$ and $\Delta t>0$. Then, the matrix $(x_{1},x_{2})=1-\tau_{x}$ (exploiting that the eigenvectors are only defined up to an overall prefactor) and we find

[TABLE]

where $a_{34}$ follows from the expression for $a_{12}$ by the replacements $(1,2)\leftrightarrow(3,4)$. As expected, this expression is again consistent with the particle-hole symmetry of Majorana end states.

Now, consider the nontopological phase in the complementary parameter range $\left|2t\right|<\left|\mu\right|$. Assuming for definiteness that $\left|\lambda_{1,3}\right|>1$, we have the relation

[TABLE]

and find that all $\lambda_{j}$ are real. We then observe that the matrix $(x_{1},x_{2})$ is invertible and we can evaluate the fixed-point boundary Green function for $\omega=0$. This yields

[TABLE]

which is nonsingular indicating that there are no zero-energy end states.

IV.2.2 Topological index

The Kitaev chain with the Hamiltonian given in Eq. (55) belongs to symmetry class BDI, as it has a chiral symmetry operator $\mathcal{U}_{S}=\tau_{x}$. For 1D systems, class BDI is characterized by a topological $\mathbb{Z}$-invariant which is given by

[TABLE]

Here, $M$ is even and $\nu(F)$ counts the number of negative eigenvalues of a matrix $F$. This expression follows from Ref. Fulga et al. (2012) in conjunction with the relation (17) between reflection matrix and boundary Green function (see App. B for a summary).

The existence of zero-energy boundary modes implies that for $\omega\to 0$, the boundary Green function has the form

[TABLE]

with a singular Hermitian matrix $K$ whose rank corresponds to the number of zero modes. This general form for symmetry class BDI is consistent with the expressions (62) and (64) which we found above for the Kitaev chain. The chiral symmetry implies that $\mathcal{U}_{S}$ anticommutes with the Hamiltonian, so that $[U_{S},VKV^{\dagger}]=0$ (see App. B for more details). Thus, we can choose a basis in which the Hermitian matrices $\mathcal{U}_{S}$ and $VKV^{\dagger}$ are simultaneously diagonalized and order the eigenvalues such that

[TABLE]

with $\mathbb{I}_{n}$ the identity matrix of dimension $n$, and

[TABLE]

with $D_{\pm,p}$ a diagonal matrix containing the nonzero eigenvalues of $VKV^{\dagger}$. Here, $p$ ($q$) denotes the number of zero-energy modes with positive (negative) chirality. We also use the notation $0_{n}$ for an $n\times n$ zero matrix.

Using the asymptotic form of the boundary Green function in Eq. (68), we thus find

[TABLE]

and therefore

[TABLE]

Thus, the topological index of symmetry class BDI can be computed as the difference between the numbers of zero-energy modes with positive and negative chiralities.

In the absence of zero modes, $\lim_{\omega\to 0}G(\omega)=G(0)$ is generically an invertible, full-rank matrix. Then, chiral symmetry implies $\{\mathcal{U}_{S},VG(0)V^{\dagger}\}=0$. (Notice that unlike for the singular boundary Green function in the topological phase, the leading contribution to $G$ for $\omega\to 0$ is independent of, and thus even in, $\omega$.) Hence, we can choose a basis in which

[TABLE]

and

[TABLE]

where the diagonal $M/2\times M/2$ matrix $D_{M/2}$ collects eigenvalues of $VG(0)V^{\dagger}$. We obtain the matrix

[TABLE]

As $\{\sigma_{z}\otimes\mathbb{I}_{M/2},F\}=0$, the spectrum of $F$ must be symmetric with respect to zero. Hence, we conclude $\mathcal{Q}=0$ as expected.

For the Kitaev chain in the topological phase, the asymptotic form of the boundary Green function is $G\sim(1\pm\tau_{x})/\omega$. This leads to the topological invariant $\mathcal{Q}=\pm 1$, corresponding to the existence of one Majorana zero mode with chirality $\pm 1$. When this model is in its trivial phase, we have $G\sim\tau_{z}$ and the topological invariant becomes $\mathcal{Q}=0$.

In realistic 1D topological superconductors, the chiral symmetry will typically be absent. This lowers the symmetry class from BDI to D, which has a $\mathbb{Z}_{2}$ topological invariant. In this case, there can be at most one zero-energy boundary mode when the system is in the topological phase. Then, the asymptotic Green function $G\sim K/\omega$ involves a matrix $K$ of rank one. The topological invariant for this case is (see App. B for details)

[TABLE]

which takes the value $+1$ in the trivial phase and $-1$ in the topological phase.

In the presence of one boundary zero mode, diagonalizing the matrix $VKV^{\dagger}$ results in only one non-zero eigenvalue. In this diagonal basis, one finds

[TABLE]

which leads to $\mathcal{Q}=-1$.

In the trivial phase, $\lim_{\omega\to 0}G(\omega)=G(0)$ is generically a full rank matrix. Because of particle-hole symmetry (see App. B for further details), the spectrum of $VG(0)V^{\dagger}$ is symmetric about zero, which leads to $\mathcal{Q}=1$.

As an aside, this implies that the topological invariant for class D signals the presence or absence of Majoranas in the Kitaev chain despite its being in class BDI.

IV.3 Chern insulator

IV.3.1 Model and boundary Green function

In this section, we consider a Chern insulator as an example of a 2D topological phase. A simple model Hamiltonian for a Chern insulator on a square lattice is given by Bernevig and Hughes (2013)

[TABLE]

The corresponding spectrum is gapped except for $M=0,2,4$, and the model exhibits topological phases for $0<M<2$ – with edge states around $k=0$ – and for $2<M<4$ – with edge states around $k=\pi$.

To compute the boundary Green function, we rewrite the Hamiltonian in real space along the $x$-direction (and retain the momentum representation in the $y$-direction with periodic boundary conditions). This yields

[TABLE]

with $\Lambda(k_{y})=2-M-\cos k_{y}$. We can then readily read off

[TABLE]

It is useful to compute the fixed-point boundary Green function in a perturbation theory in $(\omega-\sin k_{y}\sigma_{y})$, starting with

[TABLE]

where $T_{0}=T(\sin k_{y}=0,\omega=0)$. For $\omega-\sin k_{y}\sigma_{y}=0$, the eigenvalues satisfy

[TABLE]

We denote the roots for the upper sign as $\lambda_{1,2}$ and for the lower sign as $\lambda_{3,4}$, and fix the labeling by imposing $\lambda_{1}\lambda_{4}^{*}=1=\lambda_{2}\lambda_{3}^{*}$. The corresponding eigenvectors are

[TABLE]

Similar to the Kitaev chain, we expect a topological phase when $|\lambda_{1}|$ and $|\lambda_{2}|$ are either both larger or both smaller than unity. The relevant parameter ranges can again be deduced from the quantity

[TABLE]

Thus, the topological phase requires $\Lambda^{2}<1$, which implies $\left|\lambda_{1,2}\right|>1$ for $B<0$ and $\left|\lambda_{3,4}\right|>1$ for $B>0$. Using that $\cos k_{y}=\pm 1$ when $\sin k_{y}=0$, we recover that the model is topological for $0<M<4$.

We can now evaluate the fixed-point boundary Green function in perturbation theory by following essentially the same steps as for the Kitaev chain. When $\left|\lambda_{1,2}\right|>1$, a straight-forward calculation yields

[TABLE]

with

[TABLE]

Similarly, when $\left|\lambda_{3,4}\right|>1$, we have

[TABLE]

where

[TABLE]

Thus, we find that when $0<M<2$, the system has a gapless chiral boundary mode near $k_{y}=0$ with dispersion $\omega=\operatorname{sgn}(B)\sin k_{y}$. When $2<M<4$, the gapless chiral boundary mode occurs around $k_{y}=\pi$ with dispersion $\omega=\operatorname{sgn}(B)\sin k_{y}$.

In the nontopological phase for $M<0$ or $M>4$, we have $\Lambda^{2}>1$ and can thus choose $\left|\lambda_{1}\right|,\left|\lambda_{3}\right|>1$. As for the Kitaev chain, the $\lambda_{j}$ are all real and obey the relation

[TABLE]

The boundary Green function is nonsingular and can hence be evaluated approximately for $\omega-\sin k_{y}=0$. This yields

[TABLE]

IV.3.2 Topological index

The Chern insulator with Hamiltonian (78) belongs to class A. The topological invariant for a two-dimensional system in this symmetry class takes the form Fulga et al. (2012) (see also App. B)

[TABLE]

where the zero-energy reflection matrix is expressed in terms of the boundary Green function as

[TABLE]

In line with Laughlin’s argument Laughlin (1981) and the familiar expression for quantum pumping of charge Brouwer (1998), this can be viewed as the pumped charge in units of the electron charge when threading the effectively cylinder-shaped sample (periodic boundary conditions in the $y$-direction) by one flux quantum.

When $0<M<2$, the boundary Green function diverges at $k_{y}=0$ and has rank one,

[TABLE]

Away from $k_{y}=0$ the boundary Green function is invertible and finite. This yields a winding number $\mathcal{Q}=1$.

When $2<M<4$, the boundary Green function diverges at $k_{y}=\pi$ and has rank one,

[TABLE]

Away from $k_{y}=\pi$ the boundary Green function is well behaved. We thus obtain the opposite winding number $\mathcal{Q}=-1$.

When $M<0$ or $M>4$, the boundary Green function is well behaved at all $k_{y}$ and the winding number is zero.

V Higher-order topological insulators

V.1 Model

Recent work introduced higher-order topological insulators Benalcazar et al. (2016). In 2D, these phases have gapped edges but gapless corner states. In this section, we show that our approach readily accommodates such systems and provides a rather transparent picture.

The model investigated in Ref. Benalcazar et al. (2016) takes the form

[TABLE]

This square-lattice Hamiltonian has dimerized hopping amplitudes in the $x$ and $y$-directions and an added $\pi$-flux per plaquette, as illustrated in Fig. 2. The hopping amplitudes $t_{\alpha,j}$ ($\alpha=x,y$ and $j=1,2$) are assumed to be real. As a consequence of the $\pi$ flux, the model has an insulating bulk which makes it amenable to our boundary Green function approach.

While the edges of this model are gapped, there are zero-energy corner states Benalcazar et al. (2016). In the presence of $C_{4}$ symmetry, these corner states can be characterized by a quadrupole moment which is a $\mathbb{Z}_{2}$ topological index. In the absence of $C_{4}$ symmetry, there is a $\mathbb{Z}_{2}\times\mathbb{Z}_{2}$ topological index which can be associated with the edge polarizations.

The main observation underlying the current section is that even if the fixed-point boundary Green function is gapped, it defines a low-energy boundary Hamiltonian in a natural manner. This boundary Hamiltonian can be either trivial or topological. Gapless corner states (or corresponding generalizations to higher dimensions) appear if this boundary Hamiltonian is topological.

V.2 Effective boundary Hamiltonian

We now derive the fixed-point boundary Green function of this model. The dimerization of the hopping amplitudes makes this a model with period $a=2$, with $T=M_{2}M_{1}$ being the product of two matrices. As described in Sec. III, we can deduce the fixed-point boundary Green function from the eigenvectors and eigenvalues of $T$.

Consider a system with $2N\times 2N$ sites and focus on the Green function for a boundary parallel to the $x$-direction (in which we thus assume periodic boundary conditions). Introducing a two-component spinor $\psi_{n,j}^{T}=(c_{2n-1,j}^{T},c_{2n,j}^{T})$ and its Fourier transform

[TABLE]

the Hamiltonian can be rewritten as

[TABLE]

Here, we define

[TABLE]

and

[TABLE]

The matrices $\tau_{x,y,z}$ are Pauli matrices in the sublattice space of even and odd sites along the $x$-direction, with $\tau_{\pm}=(\tau_{x}\pm i\tau_{y})/2$.

With these ingredients, we can construct the transfer matrix $T=M_{2}M_{1}$ using Eq. (20). As the edge is gapped, we can set $\omega=0$ from the outset. This yields

[TABLE]

The matrix $T$ decomposes into two independent $2\times 2$ blocks, which we can treat separately. Both blocks have the same eigenvalues $\lambda$ (with $|\lambda|>1$) and $1/\lambda^{*}$. We find that $\lambda$ fulfills the equation

[TABLE]

Evaluating the corresponding eigenvectors, we can then compute the fixed-point boundary Green function based on Eq. (45). This yields

[TABLE]

which is invertible and purely real, reflecting the fact that this describes the gapped edge of a gapped bulk.

We can now define an effective boundary Hamiltonian

[TABLE]

from the fixed-point boundary Green function. Here, the prefactor

[TABLE]

has a fixed sign as a function of $k$. It is our central observation that up to the prefactor $f_{k}$, this is just the Hamiltonian of the SSH model (6). Thus, the boundary Hamiltonian is topologically equivalent to the bulk Hamiltonian of the SSH model, with $t_{1}$ and $t_{2}$ replaced by $t_{x,1}$ and $t_{x,2}$.

A similar analysis can be performed for the edge parallel to the $y$-direction. In view of the fact that the $\pi$-flux can be equally included via the hopping amplitudes in the $x$-direction, the calculation is entirely equivalent and just results in exchanging the subscripts $x$ and $y$. Thus, the edge along the $y$-direction can also be characterized by a $\mathbb{Z}_{2}$ invariant. Altogether, model is thus characterized by a $\mathbb{Z}_{2}\times\mathbb{Z}_{2}$ invariant, which encodes the polarizations of the edges in the two directions Benalcazar et al. (2016).

The various cases can be illustrated by plotting the wave function (modulus squared) of the excitation which is closest to zero energy, see Fig. 3. When the boundary Hamiltonians for both the $x$- and the $y$-edges are in the topological phase, the model exhibits corner states [panel (a)]. The excitation corresponds to a pure bulk state when both edge Hamiltonians are in the nontopological phase [panel (b)]. When only one of the two boundary Hamiltonians is in the topological phase, the excitations are localized at the boundary in the topological direction, but delocalized in the other one [panels (c) and (d)].

VI Conclusion

The boundary Green function is a natural quantity to characterize topological insulators and superconductors as it directly encodes the boundary modes and hence the topological phase diagram as well as the topological indices. In this paper, we have developed a systematic approach to compute boundary Green functions of topological insulators and superconductors by relying on a recursive approach. We show for several familiar models that our approach can be applied in a rather straight-forward manner.

Though simpler, our recursive approach has some similarities with a real-space renormalization group approach. For the noninteracting fermion lattice problems discussed in this paper, we find that the recursion relation can be solved explicitly in terms of the transfer matrix which greatly facilitates analytical calculations. When attaching leads to the sample, we can relate the boundary Green function to the reflection matrix. In view of earlier work on deriving the periodic table of topological phases of free fermions from the reflection matrix Fulga et al. (2012), this provides us with explicit expressions for the topological indices in terms of the boundary Green functions.

As an interesting application, we show that for insulators with trivial bulk topology, the boundary Green function defines a boundary Hamiltonian in a natural manner. Remarkably, even if the bulk is nontopological, this boundary Hamiltonian can still be topologically nontrivial which, say in 2D, is reflected in gapless corner states. We show this explicitly for a model which was recently proposed. This not only provides an intuitive approach to these recently introduced topological systems, but also has some potential to serve as a starting point for a more systematic investigation of these higher-order topological insulators.

While we restricted ourselves to clean and noninteracting models, the recursive approach should be extendable to include disorder and possibly interactions. For instance, in the presence of disorder, the recursion for the boundary Green function remains valid. While for a fixed disorder configuration, the boundary Green function would no longer approach a fixed-point Green function in the thermodynamic limit, the distribution of boundary Green functions over the disorder ensemble should still have such a fixed point.

Acknowledgements.

We thank Piet Brouwer and Andrei Bernevig for a critical reading of the manuscript, and acknowledge financial support by the Deutsche Forschungsgemeinschaft CRC 183. Yimu Bao thanks the Dahlem Center for hospitality during a research internship.

Appendix A Boundary Green function and reflection matrix

In this appendix, we present an alternative derivation of the relation (17) between reflection matrix and boundary Green function.

We start by noting that if we multiply $\psi_{n}^{\dagger}$ from the left in Eq. (19), we find that $\mathrm{Im}\psi_{n}^{\dagger}V_{n}^{\dagger}\psi_{n+1}$ is independ of $n$. Thus, we can define the current operator

[TABLE]

so that the current $\mathcal{I}=\Psi(n)^{\dagger}\hat{\mathcal{I}}_{n}\Psi(n)$ becomes independent of $n$. Here, we use the notation $\Psi(n)=[\psi(n+1)^{T},\psi(n)^{T}]^{T}$.

It is more convenient to have a current operator which is independent of cell index $n$. To this end, we introduce the transformation

[TABLE]

with $O_{n}=\operatorname{diag}(1,V_{n})$. Here, the Pauli matrices $\sigma_{x,y,z}$ act in the two-component space of $\Psi(n)$. In the transformed basis, the matrix $M_{n}$ defined in Eq. (20) becomes

[TABLE]

so that the relation $\Psi(n+1)=M_{n}\Psi(n)$ remains unchanged.

We can further diagonalize the current operator by introducing a unitary transformation

[TABLE]

such that

[TABLE]

Then, the current operator becomes

[TABLE]

In this new basis, the upper and lower components (of length $M$ each) of $\Psi(n)$ describe the right- and left-moving states at unit cell $n$.

If we attach normal metal leads from left and right to the first and the last unit cell of the quasi-1D system of $N$ cells, we can define the scattering matrix

[TABLE]

where the $2\times 2$ block structure reflects the two leads, $r_{N},r^{\prime}_{N}$ are reflection matrices and $t_{N},t^{\prime}_{N}$ are transmission matrices. The scattering matrix relates the outgoing states on left and right to the incoming states. After the basis transformations, we have

[TABLE]

Here the upper and lower components are outgoing (incoming) and incoming (outgoing) states on the right (left).

In the new basis, the transfer matrix can be written as Nazarov and Blanter (2009)

[TABLE]

If we add one unit cell to the system, we have $\mathcal{M}_{N+1}=M_{N+1}\mathcal{M}_{N}$, where $M_{N}$ in the new basis becomes

[TABLE]

by straightforward algebra. We then find the recurrence relations

[TABLE]

for the reflection and transmission matrices. By adding and subtracting these two equations, we find

[TABLE]

From the last equation, we have

[TABLE]

Inserting this expression back into the previous one, we find

[TABLE]

Multiplying $t_{N+1}$ from the right, we obtain

[TABLE]

Defining

[TABLE]

we have

[TABLE]

Comparing with Eq. (4), we find that the boundary Green function can be written as

[TABLE]

Solving for $r_{N}$, we recover the relation (17) between the reflection matrix $r_{N}$ and the boundary Green function $G_{N}$.

Appendix B Topological invariants

Single particle Hamiltonians can be classified based on nonspatial and nonunitary symmetries, i.e., the time-reversal, particle-hole, and chiral symmetries. This defines the ten Altland-Zirnbauer symmetry classes Altland and Zirnbauer (1997). Ref. Fulga et al. (2012) discussed the topological classification of the ten symmetry classes in terms of the reflection matrix at the Fermi energy. In this section, we adapt this discussion to the language of the boundary Green functions and provide explicit expressions for the topological indices.

B.1 Bott periodicity

B.1.1 Symmetries

Time-reversal symmetry requires the existence of an antiunitary operator $\hat{\mathcal{T}}=\mathcal{U}_{T}\hat{\mathcal{K}}$ with a unitary matrix $\mathcal{U}_{T}$ and complex conjugation operator $\hat{\mathcal{K}}$, such that

[TABLE]

Correspondingly, the boundary Green function obeys the relation

[TABLE]

which can be proven iteratively using the Dyson equation [see Eq. (4)]. If we apply the time-reversal operator twice, we obtain

[TABLE]

By Schur’s lemma, we conclude that $\mathcal{U}_{T}^{*}\mathcal{U}_{T}=\exp(i\alpha)\mathbb{I}$ is a multiple of the identity matrix. Moreover, because of the unitarity of $\mathcal{U}_{T}$, we find $\exp(2i\alpha)=\pm 1$. Thus, there are two types of time-reversal symmetries,

[TABLE]

or $\mathcal{T}^{2}=\pm 1$.

Particle-hole symmetry requires the existence of an antiunitary operator $\hat{\mathcal{C}}=\mathcal{U}_{C}\hat{\mathcal{K}}$ such that

[TABLE]

Here, $\mathcal{U}_{C}$ is a unitary matrix. The constraint of particle-hole symmetry on the boundary Green function is

[TABLE]

Similar to the case of time-reversal symmetry, we find that there are two types of particle-hole symmetries given by

[TABLE]

or $\hat{\mathcal{C}}^{2}=\pm 1$.

Chiral symmetry requires the existence of a unitary operator $\hat{\mathcal{S}}=\mathcal{U}_{S}$ represented by a unitary matrix $\mathcal{U}_{S}$, such that

[TABLE]

Because of chiral symmetry, the boundary Green function fulfills

[TABLE]

Applying chiral symmetry twice, we find that $\mathcal{U}_{S}^{2}=\exp(i\alpha)\mathbb{I}$. Redefining $\mathcal{U}_{S}\to\mathcal{U}_{S}\exp(-i\alpha/2)$, we have

[TABLE]

which implies $\mathcal{U}_{S}=\mathcal{U}_{S}^{\dagger}$. When both time-reversal and particle-hole symmetry exist, the system also has chiral symmetry with $\mathcal{U}_{S}=\mathcal{U}_{T}\mathcal{U}_{C}^{*}$. If either time-reversal or particle-hole symmetry is absent, chiral symmetry is also absent. However, chiral symmetry can exist when both time-reversal and particle-hole symmetry are absent.

The absence or presence of time-reversal ($\hat{\mathcal{T}}$), particle-hole ($\hat{\mathcal{C}}$), and chiral symmetries ($\hat{\mathcal{S}}$) of the various kinds according to Eqs. (133), (136), and (139) leads to the ten Altland-Zirnbauer symmetry classes, as summarized in the first two columns of Table 1.

B.1.2 Dimensional reduction

In all symmetry classes with a chiral symmetry, the reflection matrix (the subscript $N$ is omitted for simplicity) satisfies the constraint

[TABLE]

which follows from Eqs. (17) and (138). This implies that in the limit $\omega\to 0$, the matrix $\mathcal{U}_{S}r$ is hermitian where we use the shorthand $r=\lim_{\omega\to 0}r(\omega)$. This enables one to define an effective Hamiltonian in $d-1$ dimensions as Fulga et al. (2012)

[TABLE]

This Hamiltonian is again gapped as it has eigenvalues $\pm 1$ and has no chiral symmetry.

In the absence of chiral symmetry, the eigenvalues of the reflection matrix $r$ at the Fermi energy are in general complex. One can then construct an effective Hamiltonian in $d-1$ dimensions by doubling the degrees of freedom, Fulga et al. (2012)

[TABLE]

Because of the identity

[TABLE]

the effective Hamiltonian acquires a chiral symmetry with $\mathcal{U}_{S}=\sigma_{z}\otimes\mathbb{I}$.

*Classes A and AIII.—*Hamiltonians in classes A and AIII are thus transformed into one another under this process of dimensional reduction. This is referred to as complex Bott periodicity (mod $2$).

Now we take time-reversal and particle-hole symmetries into account, in order to reproduce the real Bott periodicity for the remaining eight classes of the ten-fold table.

*Classes AI and AII.—*We first focus on classes AI and AII with only time-reversal symmetry. Let us denote $r(\mathbf{k}_{\perp})=\lim_{\omega\to 0}r(\mathbf{k}_{\perp},\omega)$. By Eqs. (130) and (131), time-reversal symmetry restricts the reflection matrix at the Fermi energy to obey

[TABLE]

Taking the hermitian conjugate and multiplying by $\mathcal{U}_{T}$ from both sides, we obtain

[TABLE]

According to Eq. (147), the effective Hamiltonian $H_{d-1}$ has chiral symmetry. In addition, we can confirm using Eqs. (149) and (150) that is also acquires time-reversal symmetry with $\hat{\mathcal{T}}=(\sigma_{x}\otimes\mathcal{U}_{T})\hat{\mathcal{K}}$ since

[TABLE]

and particle-hole symmetry with $\hat{\mathcal{C}}=(i\sigma_{y}\otimes\mathcal{U}_{T})\hat{\mathcal{K}}$ since

[TABLE]

Notice that time-reversal and particle-hole symmetries combine into the chiral symmetry as expected. The type of time-reversal and particle-hole symmetry of the effective Hamiltonian $H_{d-1}$ can be determined from

[TABLE]

Hence, dimensional reduction transforms the classes AI and AIII into CI and DIII, respectively.

*Classes C and D.—*Next, consider classes C and D with only particle-hole symmetry. Similar to the previous case, particle-hole symmetry implies that the reflection matrix at the Fermi energy obeys

[TABLE]

The effective Hamiltonian $H_{d-1}$ according to Eq. (147) has chiral symmetry. In addition, we can establish using Eqs. (155) and (156) that $H_{d-1}$ has time-reversal symmetry with $\hat{\mathcal{T}}=(\mathbb{I}_{2}\otimes\mathcal{U}_{T})\hat{\mathcal{K}}$ since

[TABLE]

and particle-hole symmetry with $\hat{\mathcal{C}}=(\sigma_{z}\otimes\mathcal{U}_{T})\hat{\mathcal{K}}$ since

[TABLE]

Time-reversal and particle-hole symmetries combine as

[TABLE]

and hence produce the chiral symmetry as expected.

The types of the new time-reversal and particle-hole symmetry are determined from

[TABLE]

Hence, classes C and D transform into CII and BDI respectively under the procedure of dimensional reduction.

*Classes CI, CII, DIII, and BDI.—*Finally, we focus on the symmetry classes where time-reversal, particle-hole, and chiral symmetries are all present. In this situation, we have $\mathcal{U}_{S}=\mathcal{U}_{T}\mathcal{U}_{C}^{*}$. Then, $\mathcal{U}_{S}^{2}=1$ implies $\mathcal{U}_{T}^{*}\mathcal{U}_{C}\mathcal{U}_{S}^{*}=1$. This can be used to show that

[TABLE]

Notice that $\mathcal{U}_{C}^{*}\mathcal{U}_{C}=\pm 1$ and $\mathcal{U}_{T}^{*}\mathcal{U}_{T}=\pm 1$ are just numbers.

The effective Hamiltonian $H_{d-1}$, in this case defined by Eq. (141), has no chiral symmetry and can thus have either time-reversal or particle-hole symmetry, but not both. Using Eq. (162), we find that $H_{d-1}$ has the property

[TABLE]

Thus, $\mathcal{U}_{C}\hat{\mathcal{K}}$ defines a time-reversal symmetry or particle-hole symmetry for $H_{d-1}$, depending on whether the quantity $(\mathcal{U}_{C}^{*}\mathcal{U}_{C})(\mathcal{U}_{T}^{*}\mathcal{U}_{T})$ takes on the value $+1$ or $-1$, respectively. Whether this symmetry squares to $+1$ or $-1$ is determined by $\mathcal{U}_{C}^{*}\mathcal{U}_{C}$. This implies that under the procedure of dimensional reduction, classes CI, CII, DIII, and BDI transform into C, AII, D, and AI, respectively.

Using this procedure of dimensional reduction, a $d$-dimensional Hamiltonian in one symmetry class is related to a $(d-1)$-dimensional Hamiltonian $H_{d-1}$ in another symmetry class. Yet, both Hamiltonians have the same topological invariants Fulga et al. (2012). This reproduces the Bott periodicity of the topological classification of symmetry classes Schnyder et al. (2008); Kitaev (2009); Ryu et al. (2010). In particular, classes A and AIII are transformed into one another, producing the complex Bott periodicity (mod $2$). The remaining eight classes with antiunitary symmetries are shifted by one in Table 1, giving rise to the real Bott periodicity (mod $8$). This is summarized in Table 1.

B.2 Topological invariants

From 1D Hamiltonians, dimensional reduction produces 0D Hamiltonians. For classes without symmetry between positive and negative eigenenergies (classes A, AI, and AII), the topological invariant in 0D is given by the imbalance between positive and negative eigenenergies (measured relative to the Fermi energy). This corresponds to a $\mathbb{Z}$ index, except in class AII, where this number is always even because of Kramers degeneracy ($2\mathbb{Z}$ index). Of the remaining classes with symmetry between positive and negative-energy spectra, D and BDI are special as there is no level repulsion between the levels of a pair with positive and negative energies Altland and Zirnbauer (1997); Beenakker (2015). These classes have a $\mathbb{Z}_{2}$ topological invariant, which is associated with fermion parity and can be expressed as the sign of a Pfaffian

[TABLE]

where the Hamiltonian $H_{0}$ is the zero-dimensional Hamiltonian, which is antisymmetric in an appropriate basis. The remaining five symmetry classes are always topologically trivial.

B.2.1 1D systems

The dimensional reduction implies that in one dimension, there are topological phases in symmetry classes AIII, BDI, CII, DIII, and D. The topological invariants in 0D immediately provide explicit expressions for boundary topological invariants in 1D.

Dimensional reduction connects the classes AIII, BDI, and CII in 1D to classes with $\mathbb{Z}$ or $2\mathbb{Z}$ indices in 0D. Using Eq. (141), the 1D classes are therefore characterized by the boundary topological invariant

[TABLE]

where $\nu(A)$ denotes the number of negative eigenvalues of the hermitian matrix $A$. The constant $M/2$ ensures that the trivial phase has topological invariant $\mathcal{Q}=0$. This invariant indeed counts the imbalance between positive and negative eigenenergies of $H_{d-1}$.

By the bulk-boundary correspondence, these topological invariants can also be computed from the bulk 1D Hamiltonian $H(k)$. Due to the presence of chiral symmetry, the bulk Hamiltonian can be written as

[TABLE]

and the corresponding bulk topological invariant is given by Zak (1989); Ryu and Hatsugai (2002)

[TABLE]

which is a winding number.

In class DIII, $H_{d-1}$ has particle-hole symmetry realized via a unitary matrix $\mathcal{U}_{C}$ with $\mathcal{U}_{C}^{*}\mathcal{U}_{C}=1$. This implies that $\mathcal{U}_{C}=\mathcal{U}_{C}^{T}$, so that it is possible to write $\mathcal{U}_{C}=\mathcal{V}_{C}\mathcal{V}_{C}^{T}$ with a unitary matrix $\mathcal{V}_{C}$. By the unitary transformation

[TABLE]

$H_{d-1}$ becomes antisymmetric and purely imaginary. According to Eq. (164), the boundary topological invariant for class DIII can thus be expressed as

[TABLE]

Again, this invariant can also be computed from the bulk Hamiltonian $H(k)$ as Qi et al. (2010)

[TABLE]

where $h(k)$ is defined as in Eq. (166).

Class D also allows us to transform $H_{d-1}$ as

[TABLE]

where $\mathcal{V}_{C}^{\dagger}r\mathcal{V}_{C}$ is purely real. Thus, the boundary topological invariant takes the form

[TABLE]

according to Eq. (164). The corresponding bulk topological invariant is Kitaev (2001)

[TABLE]

B.2.2 2D systems

By the procedure of dimensional reduction, the boundary invariants of 2D topological phases correspond to bulk invariants of 1D Hamiltonians. Using the 1D bulk invariant for classes AIII, BDI, and CII in Eq. (167), we obtain the boundary topological invariants for 2D systems in classes A, C, and D expressed in terms of the boundary Green function,

[TABLE]

These boundary topological invariants can again alternatively be computed from a bulk topological index expressed in terms of the bulk Hamiltonian. This bulk topological index takes the form of a Chern number Chiu et al. (2016),

[TABLE]

where the Berry curvature is defined as

[TABLE]

Here, $\ket{u_{\alpha}(\mathbf{k})}$ denotes the Bloch wave function of the $\alpha$ band. Note that in computing the Berry curvature, one sums over occupied bands only.

For 2D system in class DIII, the boundary topological invariant follows from the 1D bulk invariant for class D in Eq. (178). This yields

[TABLE]

The corresponding bulk topological invariant of the bulk Hamiltonian $H(k_{1},k_{2})$ can be written as a product of topological invariants of 1D Hamiltonians Qi et al. (2010)

[TABLE]

where $\mathcal{Q}(H(k_{1},k_{2}))$ with $k_{2}=0,\pi$ is given in Eq. (169).

Finally, the topological invariant of 2D system in class AII is given by

[TABLE]

according to the 1D bulk topological invariant for class DIII in Eq. (169). This topological invariant can be computed from the bulk via Fu and Kane (2006)

[TABLE]

where the matrix element of $w$ is defined as

[TABLE]

with occupied Bloch bands $\ket{u_{n}(\mathbf{k})}$.

B.2.3 3D systems

The 3D topological invariant in class AII, which transforms into DIII under dimensional reduction, takes the form

[TABLE]

with

[TABLE]

The topological invariants in classes AIII, CI, DIII, and CII can be written in terms of 2D bulk topological invariants. For simplicity, we do not write them down explictly here.

Appendix C Condition for a vanishing bulk gap

In this appendix, we show that when the bulk band gap vanishes, the matrix $T$ defined in Eq. (46) has at least one eigenvalue with unit modulus and thus one vanishing Lyapunov exponent. This implies that the topological phase transition with its associated gap closing is signalled by at least one vanishing Lyapunov exponent.

Consider a homogeneous, quasi-1D system with periodic boundary conditions. According to the Bloch theorem, we can write the wave function as as $\psi(n)=e^{ikn}\mathbf{u}$, where $\mathbf{u}$ is an $M$-component column vector. Using Eq. (19), we find a zero-energy eigenstate when

[TABLE]

This implies

[TABLE]

Thus, $\lambda=e^{-ik}$ is an eigenvalue of $T(\omega=0)$ defined in Eq. (46), which implies that the corresponding Lyapunov exponent vanishes.

Appendix D Derivation of Eq. (45)

From Eq. (42) and (44), we have

[TABLE]

which gives

[TABLE]

Hence

[TABLE]

Taking the thermodynamic limit $N\to\infty$, we have $\Lambda^{-2N}\to 0$ and thus arrive at Eq. (45).