A generalization of Bloch's theorem for arbitrary boundary conditions: Theory

TL;DR

This paper generalizes Bloch's theorem for finite-range lattice systems with arbitrary boundary conditions, providing exact solutions, boundary indicators, and insights into topological edge modes and localized excitations.

Contribution

It introduces an exact, analytic framework for solving lattice Hamiltonians with arbitrary boundaries, extending Bloch's theorem and linking bulk-boundary properties to topological phenomena.

Findings

Exact solutions for energy eigenvalues and eigenstates with arbitrary boundaries.

Identification of conditions for robust zero-energy edge modes.

Demonstration of fractional Josephson effect without fermionic parity switch.

Abstract

We present a generalization of Bloch's theorem to finite-range lattice systems of independent fermions, in which translation symmetry is broken only by arbitrary boundary conditions, by providing exact, analytic expressions for all energy eigenvalues and eigenstates. By transforming the single-particle Hamiltonian into a corner-modified banded block-Toeplitz matrix, a key step is a bipartition of the lattice, which splits the eigenvalue problem into a system of bulk and boundary equations. The eigensystem inherits most of its solutions from an auxiliary, infinite translation-invariant Hamiltonian that allows for non-unitary representations of translation symmetry. A reformulation of the boundary equation in terms of a boundary matrix ensures compatibility with the boundary conditions, and determines the allowed energy eigenstates. We show how the boundary matrix captures the interplay between bulk and boundary properties, leading to efficient indicators of bulk-boundary correspondence. Remarkable consequences of our generalized Bloch theorem are the engineering of Hamiltonians that host perfectly localized, robust zero-energy edge modes, and the predicted emergence, e.g. in Kitaev's chain, of localized excitations whose amplitudes decay exponentially with a power-law prefactor. We further show how the theorem yields diagonalization algorithms for the class of Hamiltonians under consideration, and use the proposed bulk-boundary indicator to characterize the topological response of a multi-band time-reversal invariant s-wave superconductor under twisted boundary conditions, showing how a fractional Josephson effect can occur without a fermionic parity switch. Finally, we establish connections to the transfer matrix method and demonstrate, using the paradigmatic Kitaev's chain example, that a non-diagonalizable transfer matrix signals the presence of solutions with a power-law prefactor.

Figures

Equations

$$\displaystyle\widehat{\bm{H}}=\sum_{{\mathbf{r}}}\sum_{{\mathbf{j}}}\big{[}\hat{\Phi}^{\dagger}_{{\mathbf{j}}}K_{{\mathbf{r}}}\hat{\Phi}^{\;}_{{\mathbf{j}}+{\mathbf{r}}}+\frac{1}{2}(\hat{\Phi}^{\dagger}_{{\mathbf{j}}}\Delta_{{\mathbf{r}}}\hat{\Phi}_{{\mathbf{j}}+{\mathbf{r}}}^{\dagger}+\text{h.c.})\big{]},\quad$$

$${\hat{\Phi }}_{\mathbf{k}}^{†}\equiv \sum _{\mathbf{j}}\frac{e^{i\mathbf{k}\cdot \mathbf{j}}}{\sqrt{M}}{\hat{\Phi }}_{\mathbf{j}}^{†}$$

$$\begin{gathered} \hat{\mathbf{H}}=\frac{1}{2}\sum _{\mathbf{k}}[{\hat{\Phi }}_{\mathbf{k}}^{†}{K}_{\mathbf{k}}{\hat{\Phi }}_{\mathbf{k}}^{}+{\hat{\Phi }}_{-\mathbf{k}}^{†}{K}_{-\mathbf{k}}^{\ast }{\hat{\Phi }}_{-\mathbf{k}}^{} \\ +{\hat{\Phi }}_{\mathbf{k}}^{†}{\Delta }_{\mathbf{k}}{\hat{\Phi }}_{-\mathbf{k}}^{†}+{\hat{\Phi }}_{\mathbf{k}}{\Delta }_{-\mathbf{k}}^{\ast }{\hat{\Phi }}_{-\mathbf{k}}], \end{gathered}$$

$$K_{\mathbf{k}}\equiv \sum _{\mathbf{r}}{e}^{i\mathbf{k}\cdot \mathbf{r}}{K}_{\mathbf{r}},{\Delta }_{\mathbf{k}}\equiv \sum _{\mathbf{r}}{e}^{i\mathbf{k}\cdot \mathbf{r}}{\Delta }_{\mathbf{r}}.$$

$$\mathbf{j}=j{\mathbf{a}}_1+{\mathbf{j}}_{\perp },j=1,\dots ,N=N_1,{\mathbf{j}}_{\perp }=\sum _{\mu =2}^D{j}_{\mu }{\mathbf{a}}_{\mu }.$$

$${\hat{\Phi }}_{{\mathbf{k}}_{\perp }}^{†}=\sqrt{N}\sum _{{\mathbf{j}}_{\perp }}\frac{e^{i{\mathbf{k}}_{\perp }\cdot {\mathbf{j}}_{\perp }}}{\sqrt{M}}{\hat{\Phi }}_{{\mathbf{j}}_{\perp }}^{†},{\mathbf{k}}_{\perp }=\sum _{\mu =2}^D\frac{k_{\mu }}{N_{\mu }}{\mathbf{b}}_{\mu },$$

$${\hat{\Phi }}_{{\mathbf{j}}_{\perp }}^{†}\equiv \left[\begin{array}{ccc}{\hat{\Phi }}_{{\mathbf{a}}_1+{\mathbf{j}}_{\perp }}^{†}& {\hat{\Phi }}_{2{\mathbf{a}}_1+{\mathbf{j}}_{\perp }}^{†}& {\hat{\Phi }}_{N{\mathbf{a}}_1+{\mathbf{j}}_{\perp }}^{†}\end{array}\right].$$

$$\begin{gathered} {\hat{H}}_{N,{\mathbf{k}}_{\perp }}=\frac{1}{2}({\hat{\Phi }}_{{\mathbf{k}}_{\perp }}^{†}{K}_{{\mathbf{k}}_{\perp }}{\hat{\Phi }}_{{\mathbf{k}}_{\perp }}^{}+{\hat{\Phi }}_{-{\mathbf{k}}_{\perp }}^{†}{K}_{-{\mathbf{k}}_{\perp }}^{\ast }{\hat{\Phi }}_{-{\mathbf{k}}_{\perp }}^{} \\ +{\hat{\Phi }}_{{\mathbf{k}}_{\perp }}^{†}{\Delta }_{{\mathbf{k}}_{\perp }}{\hat{\Phi }}_{-{\mathbf{k}}_{\perp }}^{†}+{\hat{\Phi }}_{{\mathbf{k}}_{\perp }}{\Delta }_{-{\mathbf{k}}_{\perp }}^{\ast }{\hat{\Phi }}_{-{\mathbf{k}}_{\perp }}), \end{gathered}$$

$$[K_{{\mathbf{k}}_{\perp }}{]}_{j{j}^{\prime }}=K_{j^{\prime }-j,{\mathbf{k}}_{\perp }}\equiv \sum _{{\mathbf{r}}_{\perp }}{e}^{i\mathbf{k}\cdot \mathbf{r}}{K}_{\mathbf{r}},\mathbf{r}=(j^{\prime }-j){\mathbf{a}}_1+{\mathbf{r}}_{\perp },$$

$$K_{r,{\mathbf{k}}_{\perp }}={\Delta }_{r,{\mathbf{k}}_{\perp }}=0,\forall {\mathbf{k}}_{\perp }\text{if}|r|>R.$$

$$\mathbf{j}=b{\mathbf{a}}_1+{\mathbf{j}}_{\perp },b=1,\dots ,R,N-R+1,\dots ,N;$$

$$\begin{gathered} [{\hat{W}}_{{\mathbf{k}}_{\perp }}{]}_{b{b}^{\prime }}=\frac{1}{2}({\hat{\Phi }}_{b,{\mathbf{k}}_{\perp }}^{†}{W}_{{\mathbf{k}}_{\perp }}^{(K)}{\hat{\Phi }}_{b^{\prime },{\mathbf{k}}_{\perp }}^{}+{\hat{\Phi }}_{b,-{\mathbf{k}}_{\perp }}^{†}(W_{-{\mathbf{k}}_{\perp }}^{(K)}{)}^{\ast }{\hat{\Phi }}_{b^{\prime },-{\mathbf{k}}_{\perp }}^{} \\ +{\hat{\Phi }}_{b,{\mathbf{k}}_{\perp }}^{†}{W}_{{\mathbf{k}}_{\perp }}^{(\Delta )}{\hat{\Phi }}_{b^{\prime },-{\mathbf{k}}_{\perp }}^{†}+{\hat{\Phi }}_{b,{\mathbf{k}}_{\perp }}(W_{-{\mathbf{k}}_{\perp }}^{(\Delta )}{)}^{\ast }{\hat{\Phi }}_{b^{\prime },-{\mathbf{k}}_{\perp }}), \end{gathered}$$

$$\hat{H}={\hat{H}}_N+\hat{W}=\sum _{{\mathbf{k}}_{\perp }}{\hat{H}}_{{\mathbf{k}}_{\perp }},{\hat{H}}_{{\mathbf{k}}_{\perp }}={\hat{H}}_{N,{\mathbf{k}}_{\perp }}+{\hat{W}}_{{\mathbf{k}}_{\perp }}.$$

$${\hat{H}}_{{\mathbf{k}}_{\perp }}=\frac{1}{2}{\hat{\Psi }}_{{\mathbf{k}}_{\perp }}^{†}{\tilde{H}}_{{\mathbf{k}}_{\perp }}{\hat{\Psi }}_{{\mathbf{k}}_{\perp }}+\frac{1}{2}\text{tr}(K_{{\mathbf{k}}_{\perp }}+W_{{\mathbf{k}}_{\perp }}^{(K)})$$

$${\tilde{H}}_{N,{\mathbf{k}}_{\perp }}=\left[\begin{array}{cc}{K}_{{\mathbf{k}}_{\perp }}& {\Delta }_{{\mathbf{k}}_{\perp }}\\ -{\Delta }_{-{\mathbf{k}}_{\perp }}^{\ast }& -K_{-{\mathbf{k}}_{\perp }}^{\ast }\end{array}\right],$$

$${\tilde{W}}_{{\mathbf{k}}_{\perp }}=\left[\begin{array}{cc}{W}_{{\mathbf{k}}_{\perp }}^{(K)}& W_{{\mathbf{k}}_{\perp }}^{(\Delta )}\\ {-W_{-{\mathbf{k}}_{\perp }}^{(\Delta )}}^{\ast }& {-W_{-{\mathbf{k}}_{\perp }}^{(K)}}^{\ast }\end{array}\right].$$

$${\hat{\Psi }}_{{\mathbf{k}}_{\perp }}^{†}\equiv \left[\begin{array}{ccc}{\hat{\Psi }}_{1,{\mathbf{k}}_{\perp }}^{†}& \cdots & {\hat{\Psi }}_{N,{\mathbf{k}}_{\perp }}^{†}\end{array}\right],{\hat{\Psi }}_{j,{\mathbf{k}}_{\perp }}^{†}\equiv \left[\begin{array}{cc}{\hat{\Phi }}_{j,{\mathbf{k}}_{\perp }}^{†}& {\hat{\Phi }}_{j,-{\mathbf{k}}_{\perp }}^{}\end{array}\right],$$

$${\tilde{H}}_{{\mathbf{k}}_{\perp }}\mapsto H_{{\mathbf{k}}_{\perp }}\equiv H_{N,{\mathbf{k}}_{\perp }}+W_{{\mathbf{k}}_{\perp }},$$

$$h_r=\left[\begin{array}{cc}{K}_{r,{\mathbf{k}}_{\perp }}& {\Delta }_{r,{\mathbf{k}}_{\perp }}\\ -{\Delta }_{r,-{\mathbf{k}}_{\perp }}^{\ast }& -K_{r,-{\mathbf{k}}_{\perp }}^{\ast }\end{array}\right],\text{vspace}-0.5cm$$

$$[W{]}_{b{b}^{\prime }}=\left[\begin{array}{cc}{W}_{b{b}^{\prime },{\mathbf{k}}_{\perp }}^{(K)}& W_{b{b}^{\prime },{\mathbf{k}}_{\perp }}^{(\Delta )}\\ -(W_{b{b}^{\prime },-{\mathbf{k}}_{\perp }}^{(\Delta )}{)}^{\ast }& -(W_{b{b}^{\prime },-{\mathbf{k}}_{\perp }}^{(K)}{)}^{\ast }\end{array}\right].$$

$$H_N=\left[\begin{array}{cccccccc}{h}_0& \dots & h_R& & & 0& \cdots & 0\\ ⋮& \ddots & & \ddots & & & \ddots & ⋮\\ h_R^{†}& & \ddots & & \ddots & & & 0\\ & \ddots & & & & & & \\ & & & & & & \ddots & \\ 0& & & \ddots & & \ddots & & h_R\\ ⋮& \ddots & & & \ddots & & \ddots & ⋮\\ 0& \cdots & 0& & & h_R^{†}& \cdots & h_0\end{array}\right],$$

$$W=\left[\begin{array}{ccccccccc}{w}_{11}^{(l)}& \dots & w_{1R}^{(l)}& & 0& & w_{11}& \dots & w_{1R}\\ ⋮& \ddots & ⋮& & ⋮& & ⋮& \ddots & ⋮\\ w_{R1}^{(l)}& \dots & w_{RR}^{(l)}& & ⋮& & w_{R1}& \dots & w_{RR}\\ & & & & & & & & \\ & & & & & & & & \\ 0& \cdots & \cdots & & 0& & \cdots & \cdots & 0\\ & & & & & & & & \\ & & & & & & & & \\ w_{11}^{†}& \dots & w_{1R}^{†}& & ⋮& & w_{11}^{(r)}& \dots & w_{1R}^{(r)}\\ ⋮& \ddots & ⋮& & ⋮& & ⋮& \ddots & ⋮\\ w_{R1}^{†}& \dots & w_{RR}^{†}& & 0& & w_{R1}^{(r)}& \dots & w_{RR}^{(r)}\end{array}\right],$$

$$w_{b{b}^{\prime }}^{(l)}\equiv W_{b{b}^{\prime }},w_{b{b}^{\prime }}^{(r)}\equiv W_{N-b+1,N-b^{\prime }+1},$$

$$w_{b{b}^{\prime }}\equiv W_{b,N-b^{\prime }+1}.$$

$$\widehat{H}=\frac{1}{2}\sum_{j=1}^{N}\hat{\Psi}^{\dagger}_{j}h_{0}\hat{\Psi}^{\;}_{j}+\frac{1}{2}\sum_{r=1}^{R}\Big{(}\sum_{j=1}^{N-r}\hat{\Psi}^{\dagger}_{j}h_{r}\hat{\Psi}^{\;}_{j+r}+\text{h.c.}\Big{)}\\ +\frac{1}{2}\sum_{b,b^{\prime}}\hat{\Psi}^{\dagger}_{b}W_{bb^{\prime}}\hat{\Psi}^{\;}_{b^{\prime}}+\frac{1}{2}\text{tr}(K+W^{(K)}),$$

$$\widehat{W}=\frac{1}{2}\sum_{r=1}^{R}\sum_{b=N-R+1}^{N}\Big{(}\hat{\Psi}^{\dagger}_{b}\,g_{r}\,\hat{\Psi}^{\;}_{b+r-N}+\text{h.c.}\Big{)},$$

$$\displaystyle d\equiv\left\{\begin{array}[]{lcl}d_{\rm int}&\mbox{if}&\Delta=0\;\,\text{(number-conserving)}\\ 2d_{\rm int}&\mbox{if}&\Delta\neq 0\;\,\text{(number-non-conserving)}\end{array}\right..$$

$${\hat{H}}_N=-t\sum _{j=1}^{N-1}(c_j^{†}{c}_{j+1}+c_{j+1}^{†}{c}_j),$$

$$\displaystyle H_{N}=-t\sum_{j=1}^{N-1}\big{(}|j\rangle\langle j+1|+|j+1\rangle\langle j|\big{)},$$

$$\displaystyle H_{N}|k\rangle=-2t\cos k|k\rangle+\frac{t}{\sqrt{N}}\Big{(}|1\rangle+e^{ik(N+1)}|N\rangle\Big{)},$$

$$H_{N}|-k\rangle=-2t\cos k|-k\rangle+\frac{t}{\sqrt{N}}\Big{(}|1\rangle+e^{-ik(N+1)}|N\rangle\Big{)}.$$

Full text

I Introduction

Modern electronic transport theory in crystalline solids relies on two fundamental tenets. On the one hand, because of the Pauli exclusion principle, electrons satisfy Fermi-Dirac statistics; on the other, Bloch’s theorem allows labeling of the one-electron wave-functions in terms of their crystal momenta. The set of allowed momenta, defining the so-called Brillouin zone, is determined by symmetry and the fact that Born-von-Karman (periodic) boundary conditions (BCs) are enforced on the systemashcroftmermin . It is the organization of electrons within the Brillouin zone that is key to defining its conduction properties. While the assumption of a perfect crystal with a unit cell that is periodically repeated emphasizes the (discrete) symmetry of translation, the torus topological constraint imposed by the Born-von-Karman condition further eliminates the potential emergence of edge or boundary electronic states in a real, finite crystal. Although much of the transport properties are determined by bulk electrons, technologically relevant processes on the surface of solids are known to lead to intriguing phenomena, such as surface superconductivity gor'kov or Kondo screening of magnetic impurities resulting in exotic surface spin textures kondoleonid . Early theoretical investigations by Tamm and Shockley Tamm ; Shockley initiated the systematic study of surface state physics, that witnessed a landmark achievement with the discovery of the quantum Hall effect hall , and that today finds its most striking applications in topological insulating and superconducting materials bernevig .

The organization of bulk electrons comes with a twist. The quantum electronic states labeled by crystal momenta organize in ways subject to classification according to integer values of topological invariants defined over the entire Brillouin zone zak ; chiu16 . The first Chern number, determined in terms of the Berry connection, is one of those topological invariants, defining a topologically non-trivial electronic phase whenever its value differs from zero bernevig . For instance, the transverse conductivity of a quantum Hall fluid is proportional to such a Chern number. Perhaps surprisingly, there appears to be a connection between a non-vanishing value of the topological invariant, a bulk property, and the emergence of “robust” boundary states, an attribute of the surface. This principle is known as the bulk-boundary correspondence bernevig ; Graf13 ; chiu16 . At first, this relation seems odd, since surface properties are totally independent from those of the bulk; for example, one can deposit impurities, generate strain and reconstruction, or add externally applied electric fields only on the surface. Nonetheless, it seems reasonable to assume that as long as the symmetry protecting the surface states is not broken by external means, a bulk-boundary correspondence will still hold, although the quantum surface state will, in general, get transformed isaev11 . In other words, although the mere existence of a boundary mode may be robust, only classical information may be protected in general isaev11 ; ips .

It is apparent that Bloch’s theorem and its consequences pertain to the realm of bulk physics. A crystal without boundaries is required to establish it. But, can one generalize Bloch’s theorem for independent electrons to arbitrary BCs, so that bulk and surface states can be handled on an equal footing, and physical insight about the interplay between bulk and boundary may be gained? In light of our previous discussion, it is clear that to accomplish such a task one needs to give up on some concepts, such as the notion of a Brillouin zone. If possible, such a generalization would allow us to formulate a bulk-boundary correspondence principle that makes use of both bulk and boundary information. It is tempting to argue that the relative importance of BCs diminishes as the size of the crystal grows. Notwithstanding, for example, recent work shows that BCs impact the quasi-conserved local charges of one-dimensional systems, with important consequences for bulk quench dynamics fagotti ; gluza16 . More generally, the statistical mechanics of topologically nontrivial systems begs some answers directly relevant to the questions above quelle15 ; quelle16 .

In this paper, we generalize Bloch’s theorem to systems of independent electrons subject to arbitrary BCs. Intuitively speaking, one may expect such a result on the basis that translation symmetry is only mildly broken by BCs – namely, clean (disorder-free) systems are translationally-invariant away from the boundary. Our generalized Bloch theorem makes this idea precise, by providing an exact (often in fully closed-form) description of the eigenstates of the system’s Hamiltonian in terms of generalized eigenstates of non-unitary representations of translation symmetry in infinite space, that is, with boundaries at infinity and no torus topology commentAM . As a result, both exponentially decaying edge modes and more exotic modes with power-law prefactors can emerge, provided the BCs allow them. Our generalized Bloch theorem leverages the bulk-boundary separation of the Schrödinger equation we introduced in Ref. [abc, ] and the full solution of the bulk equation rigorously established in Ref. [JPA, ]. It extends the diagonalization procedure described in Ref. [abc, ], and recently used in Ref. [KatsuraTwisted, ], to a more general class of Hamiltonians and BCs, which in particular allows for different modifications to be imposed on different boundaries. A unifying theme behind these results is an effective analytic continuation to the complex plane of the standard Bloch’s Hamiltonian off the Brillouin zone. This analytic continuation is remarkably useful because the original problem reduces to a matrix polynomial function JPA . Interestingly, a recent study made use of similar polynomial structures for the purpose of topological classification read16 .

The outline of this paper is as follows. In Sec. II we discuss a re-arrangement of the fermionic basis that allows us to reduce the diagonalization of the original many-electron finite-range quadratic Hamiltonian in second quantization, subject to specified BCs, to the one of a single-particle Bogoliubov-de Gennes Hamiltonian that has the structure of a corner-modified block-Toeplitz matrix, as introduced in Ref. [JPA, ]. Section III develops a structural characterization of the energy eigenstates for the many-electron systems under consideration, culminating into our generalization of Bloch’s theorem. Like the usual Bloch’s theorem, such a generalization is first and foremost a practical tool for calculations, granting direct access to exact energy eigenvalues and eigenstates. In Sec. IV, we provide two new procedures – one numerical and another algebraic – for carrying out the exact diagonalization of the single-particle Hamiltonian, based on the generalized Bloch theorem. The algebraic procedure, which may provide closed-form solutions to the problem, is explicitly illustrated through a number of examples in Sec. V. While, in order to illustrate our methodology, we focus largely on one-dimensional systems here, we anticipate that additional applications to higher-dimensional problems will be addressed in a companion paper PRB2 . Remarkably, while mid-gap modes with power-law prefactors have been predicted for systems with long-range couplings, we show analytically that they can also prominently manifest in short-range tight-binding models of topological insulators and superconductors pientka13 ; degottardi13 ; vodola14 ; ortiz14 ; yazdani14 .

Crucially, our generalized Bloch theorem also allows derivation of a boundary indicator for the bulk-boundary correspondence, which contains information from both the bulk and the BCs and, as remarked in Ref. [abc, ], is computationally more efficient than other indicators also applicable in the absence of translational symmetry kitaev01 . This is the subject of Sec. VI. In the same section, we expand on the analysis of the two-band time-reversal invariant $s$-wave topological superconducting wire we introduced previously swavePRL ; swavePRB , by employing our newly defined indicator of bulk-boundary correspondence – constructed by using the generalized Bloch theorem, as opposed to the simplified Ansatz we presented in Ref. [abc, ]. Specifically, this indicator is employed in the analysis of the Josephson response of the $s$-wave superconductor in a bridge configuration, sharply diagnosing the occurrence of a fractional $4\pi$-periodic Josephson effect. Remarkably, we find that this is possible without a conventional fermionic parity switch, which we explain based on a suitable transformation into two decoupled systems, each undergoing a parity switch. Section  VII establishes some important connections between our generalized Bloch theorem and the widely employed transfer matrix approach Lee . Interestingly, from the standpoint of computing energy levels, our bulk-boundary separation is in many ways complementary to the transfer matrix method. While the latter can handle bulk disorder (at a computational cost), it does not, a priori, lend itself to investigating the space of arbitrary BCs in a transparent way. On the contrary, our generalized Bloch theorem can handle arbitrary BCs efficiently, as long as the bulk respects translational invariance – with arbitrary (finite-range) disorder on the boundary being permitted. Looking afresh at the transfer matrix approach from the generalized Bloch theorem’s perspective yields a remarkable result: the generalized eigenvectors of the transfer matrix, whose role has been appreciated only recently DwivediBB , describe energy eigenstates with power-law corrections to an otherwise exponential behavior. Our generalized Bloch theorem further suggests a way to extend the transfer matrix approach to a disordered bulk and arbitrary BCs. A discussion of the main implications of our work, along with outstanding research questions, concludes in Sec. VIII, whereas additional technical material is included in separate appendixes.

II From independent fermions to Toeplitz matrices

We begin by describing the class of model Hamiltonians investigated in this and the companion paper PRB2 . The upshot of this section will be a non-conventional re-ordering of the physical subsystems’ labels that allows recasting the single-particle (Bogoliubov-de Gennes (BdG)) Hamiltonians in Toeplitz form, essential for the exact diagonalization procedure we will describe.

Consider a $D$-dimensional, translation-invariant infinite system of independent fermions. Such a system is described in full generality by a quadratic, not necessarily particle-number-conserving, Hamiltonian in Fock space. In a lattice approximation, the vector position of a given fermion in the regular crystal lattice can be written as the sum of a Bravais lattice vector and a basis vector ashcroftmermin . We will include these basis vectors as part of the internal labels, and denote Bravais lattice vectors as ${\mathbf{j}}\equiv\sum_{\mu=1}^{D}j_{\mu}{\mathbf{a}}_{\mu}$, with ${\mathbf{a}}_{1},\dots,{\mathbf{a}}_{D}$ primitive vectors and each $j_{\mu}\in\mathbb{Z}$. An orthonormal basis of the Hilbert space of single-particle states is thus labeled by Bravais lattice vectors ${\mathbf{j}}$, and a finite number of internal labels $m=1,\dots,d_{\rm int}$. We denote by $c^{\;}_{{\mathbf{j}}m}$ ($c^{\dagger}_{{\mathbf{j}}m}$) the fermionic annihilation (creation) operator corresponding to lattice vector ${\mathbf{j}}$ and internal state $m$. The Hamiltonian of a translation-invariant system can then be written as

[TABLE]

with $\hat{\Phi}^{\dagger}_{{\mathbf{j}}}\equiv\begin{bmatrix}c^{\dagger}_{{\mathbf{j}}1}&\cdots&c^{\dagger}_{{\mathbf{j}}d_{\rm int}}\end{bmatrix}$, ${\mathbf{r}}$ a Bravais lattice vector, and the $d_{\rm int}\times d_{\rm int}$ hopping and pairing matrices $K_{{\mathbf{r}}}$, $\Delta_{{\mathbf{r}}}$ satisfying $K_{-{\mathbf{r}}}=K_{{\mathbf{r}}}^{\dagger}$, $\Delta_{-{\mathbf{r}}}=-\Delta_{{\mathbf{r}}}^{\rm T},$ where the superscript ${\rm T}$ denotes the transpose operation. For arrays, such as $\hat{\Phi}^{\dagger}_{{\mathbf{j}}}$ and $\hat{\Phi}_{{\mathbf{j}}}^{\;}$, we stick to the convention that those appearing on the left (right) of a matrix are row (column) arrays.

Since the infinite system is translation-invariant in all $D$ directions, it is customary to introduce the volume containing the electrons by imposing Born-von Karman (periodic) BCs over a macroscopic volume commensurate with the primitive cell of the underlying Bravais lattice. If the allowed ${\mathbf{j}}$’s in the macroscopic volume correspond to $j_{\mu}=1,\dots,N_{\mu}$, then,

[TABLE]

defines the Fourier-transformed array of creation operators of real Bloch wavevector (or crystal momentum), ${\mathbf{k}}\equiv\sum_{\mu=1}^{D}\frac{k_{\mu}}{N_{\mu}}{\mathbf{b}}_{\mu}$, with $k_{\mu}$ integers such that ${\mathbf{k}}$ lies inside the Brillouin zone. The total number of primitive cells is given by $M=N_{1}N_{2}\dots N_{D}$, and ${\mathbf{b}}_{\mu}$ defines the reciprocal lattice vectors satisfying ${\mathbf{a}}_{\mu}\cdot{\mathbf{b}}_{\nu}=2\pi\delta_{\mu\nu}$, with $\delta_{\mu\nu}$ representing Kronecker’s delta ashcroftmermin . Finally, by letting $*$ denote complex conjugation, one can express the Hamiltonian of Eq. (1) in momentum space as

[TABLE]

which has a block structure in terms of the matrices

[TABLE]

Now let us turn our attention to systems that are periodic along $(D-1)$ directions and terminated by two parallel hyperplanes perpendicular to the direction ${\mathbf{a}}_{1}$. We then write the allowed values of ${\mathbf{j}}$ as

[TABLE]

In this scenario, each Bloch wavevector ${\mathbf{k}}$ is no longer a good quantum number. However, we can still block-diagonalize the Hamiltonian in the partial basis

[TABLE]

where $\hat{\Phi}_{{\mathbf{j}}_{\perp}}^{\dagger}$ is defined to be the array

[TABLE]

A system with sudden termination at hyperplanes correspondig to $j=1$ and $j=N$ is modeled by open (or hardwall) BCs, in which case the Hamiltonian can be expressed as $\widehat{H}_{N}\equiv\sum_{{\mathbf{k}}_{\perp}}\widehat{H}_{N,{\mathbf{k}}_{\perp}}$,

[TABLE]

in terms of $Nd_{\rm int}\times Nd_{\rm int}$ matrices

[TABLE]

and analogously defined matrices $\Delta_{{\mathbf{k}}_{\perp}}$. We will henceforth assume that the range $R$ of hopping and pairing along the ${\mathbf{a}}_{1}$ direction is finite. This means that

[TABLE]

In this paper, we are interested in BCs more general than open BCs. They are modeled by a Hermitian many-body operator $\widehat{W}$ on Fock space which satisfies the following restrictions (see also Appendix A):

• •

$\widehat{W}$ has no effect beyond the “boundary slab”, containing basis vectors

[TABLE]

• •

$\widehat{W}$ is periodic along the $D-1$ directions ${\mathbf{a}}_{2},\dots,{\mathbf{a}}_{D}$, and has a decomposition analogous to that of $\widehat{H}_{N}$.

Because of the latter restriction, $\widehat{W}\equiv\sum_{{\mathbf{k}}_{\perp}}\widehat{W}_{{\mathbf{k}}_{\perp}}$ with

[TABLE]

where $b,b^{\prime}\in\{1,\dots,R,\ N-R+1,\dots,N\}$, $W_{{\mathbf{k}}_{\perp}}^{(K)}$ is Hermitian and $W_{{\mathbf{k}}_{\perp}}^{(\Delta)}$ is antisymmetric for each ${\mathbf{k}}_{\perp}$. Then, the model Hamiltonian, with arbitrary BCs, becomes

[TABLE]

From now on, we will focus on diagonalizing one such block $\widehat{H}_{{\mathbf{k}}_{\perp}}$, for a fixed value of ${\mathbf{k}}_{\perp}$. We will investigate the interplay between ${\mathbf{k}}_{\perp}$ and our diagonalization algorithm, (and, more generally, disordered BCs), in Ref. [PRB2, ].

The next step consists of deriving the BdG Hamiltonian for this block. The conventional way blaizot is to use the (Nambu) basis $\hat{\Psi}_{{\mathbf{k}}_{\perp}}^{\dagger}\equiv\begin{bmatrix}\hat{\Phi}_{{\mathbf{k}}_{\perp}}^{\dagger}&\hat{\Phi}_{-{\mathbf{k}}_{\perp}}\end{bmatrix}$, with $\hat{\Phi}^{\dagger}_{{\mathbf{k}}_{\perp}}$ defined in Eq. (2), so that $\widehat{H}_{{\mathbf{k}}_{\perp}}$ can be expressed in the form,

[TABLE]

in terms of a Hermitian matrix $\widetilde{H}_{{\mathbf{k}}_{\perp}}$ (note that the matrix $W_{{\mathbf{k}}_{\perp}}^{(K)}$ has entries $[W_{{\mathbf{k}}_{\perp}}^{(K)}]_{jj^{\prime}}=0$ if any of $j,j^{\prime}$ take values from the set $\{R+1,\dots,N-R\}$). This relation leads us to a BdG Hamiltonian $\widetilde{H}_{{\mathbf{k}}_{\perp}}\equiv\widetilde{H}_{N,{{\mathbf{k}}_{\perp}}}+\widetilde{W}_{{\mathbf{k}}_{\perp}}$ with

[TABLE]

The diagonalization of the BdG Hamiltonian $\widetilde{H}_{{\mathbf{k}}_{\perp}}$ implies that of $\widehat{H}_{{\mathbf{k}}_{\perp}}$, as detailed for example in Ref. [blaizot, ].

The $2\times 2$ block-structure of $\widetilde{H}_{{\mathbf{k}}_{\perp}}$ emphasizes the intrinsic charge-conjugation symmetry under the anti-unitary operator $\mathcal{C}\equiv(\mathds{1}_{Nd_{\sf int}}\tau_{x})\,{\mathcal{C}}_{\text{cc}}$, i.e., $\mathcal{C}\widetilde{H}_{{\mathbf{k}}_{\perp}}\mathcal{C}^{-1}=-{\widetilde{H}_{-{\mathbf{k}}_{\perp}}},$ where $\tau_{x}$ is the Pauli $\sigma_{x}$-matrix in the Nambu basis, and ${\mathcal{C}}_{\text{cc}}$ denotes complex conjugation. Such a block-structure, however, does not explicitly highlight the role of translation invariance. For this reason, we reorder the (Nambu) basis according to abc

[TABLE]

so that the BdG Hamiltonian transforms to

[TABLE]

in terms of a banded block-Toeplitz matrix $H_{N,{{\mathbf{k}}_{\perp}}}=H_{N}$, with entries $[H_{N}]_{jj^{\prime}}=h_{j^{\prime}-j}$ along the diagonals, and a block matrix $W_{{\mathbf{k}}_{\perp}}=W$, where

[TABLE]

[TABLE]

Explicitly, in array form, we have:

[TABLE]

[TABLE]

where we have used the notation

[TABLE]

Here, the superscript $(l)$ [or $(r)$] indicates the entries that allow hoppings only near the left [or right] boundary, whereas the ones without superscript allow hoppings from the left to the right boundary slabs. The matrix $H=H_{N}+W$ is a corner-modified banded block-Toeplitz matrix as defined in Ref. [JPA, ], and is amenable to the exact solution approach described therein RemarkSymm .

This transformed BdG Hamiltonian allows us to write the second-quantized Hamiltonian $\widehat{H}_{{\mathbf{k}}_{\perp}}$ in the form

[TABLE]

where we have dropped the label ${\mathbf{k}}_{\perp}$ everywhere. In particular, for one-dimensional systems ($D$=1), we recover (up to a constant) the class of Hamiltonians considered in Ref. [abc, ], provided that $W$ is expressible as

[TABLE]

for some $2d_{\rm int}\times 2d_{\rm int}$ matrices $g_{r}$.

Notice that for particle number-conserving systems ($\Delta=0=W^{(\Delta)}$), the single-particle Hamiltonian is just $H=K+W^{(K)}$, which is already a corner-modified, banded block-Toeplitz matrix. In such cases, the re-ordering of the basis is not required, and one may directly apply the diagonalization procedure described in the following sections to $H$, with internal blocks of dimension $d_{\rm int}$. In order to have a uniform notation, we shall use

[TABLE]

III Algebraic characterization of energy eigenstates

A main goal of this work is to diagonalize the single-particle Hamiltonian $H=H_{N}+W$, which is a corner-modified, banded block-Toeplitz matrix. In this section, we investigate the structure of its energy eigenstates, which will culminate in a generalization of Bloch’s theorem to systems described by such model Hamiltonians. Our analysis will illustrate, in particular, that for non-generic parameter values, Hamiltonians may display a finite number of exceptional (singular) energies corresponding to dispersionless, flat bands. The latter represent a macroscopic number of energy eigenstates that are localized in the bulk and, thus, are completely insensitive to BCs. It is remarkable that the analytic continuation of the Bloch Hamiltonian can still encompass this situation. We will show how to use it to construct the localized flat band energy eigenstates directly in real space.

III.1 An impurity problem as a motivating example

Consider the simple tight-binding Hamiltonian

[TABLE]

defined on an open chain of $N$ (even) lattice sites with nearest-neighbor hopping strength $t$, and lattice constant $a=1$. The corresponding single-particle Hamiltonian is

[TABLE]

and breaks translation-invariance due to the presence of the boundary, so that the crystal momentum is not a good quantum number. In fact, for any $k\in(0,2\pi]$, the state $|k\rangle=\frac{1}{\sqrt{N}}\sum_{j=1}^{N}e^{ikj}|j\rangle$ (labeled by $k$) obeys

[TABLE]

with a similar relation holding for $-k$

[TABLE]

The first term on the right-hand side of Eqs. (7)-(8) indicates that $|k\rangle$ and $|-k\rangle$ “almost” (for large $N$) satisfy the eigenvalue relation with energy $-2t\cos k$, while the two terms in the brackets show that the eigenvalue relation is violated near the two edges of the chain. Under periodic BCs, $-2t\cos k$ is the actual energy eigenvalue of the eigenstate $|k\rangle$ (and $|-k\rangle$), and $k$ is the crystal momentum, given by $k=2\pi q/N,\ q=1,\dots,N\in(0,2\pi]$ashcroftmermin .

Because of the identical first term $-2t\cos k$ in Eqs. (7) and (8), the states $|k\rangle$ and $|-k\rangle$ can be linearly combined in order to cancel off the similar-looking boundary contributions. For $\alpha,\beta\in\mathds{C}$, the eigenvalue relation

[TABLE]

is recovered provided that the constraint

[TABLE]

is satisfied. For this to hold, the coefficients of both $|1\rangle$ and $|N\rangle$ must vanish, which leads to the kernel equation

[TABLE]

The determinant of the above “boundary matrix” $B$ must vanish, which happens if the condition $e^{i2k(N+1)}=1$ is satisfied, that is, when $k={\pi q}/{(N+1)},$ $q=1,\dots,N$. For each of these values of $k$, $\alpha=-\beta=1/\sqrt{2}$ provides the required kernel vector of the boundary matrix, with the resulting $N$ eigenvectors

[TABLE]

of energy $\epsilon_{k}=-2t\cos k$. Notice that the allowed values of $k$ differ from the case of periodic BCs Mikeska .

Encouraged by these results, let us change the Hamiltonian by adding an on-site potential at the edges,

[TABLE]

so that the total single-particle Hamiltonian becomes $H=H_{N}+W$. The boundary matrix $B$ changes to

[TABLE]

While it is harder to predict analytically the values of $k$ for which it has a non-trivial kernel, it is interesting to examine the limit $w\gg t$. Then, we can approximate the relevant kernel condition as

[TABLE]

showing nontrivial solutions if $e^{i2k(N-1)}=1$. There are now $(N-2)$ $k$-values yielding stationary eigenstates as before. The two missing eigenstates are localized at the edges, and can be taken to be $|1\rangle$ and $|N\rangle$, to leading order in $t/w\ll 1$. These localized states are reminiscent of Tamm-Shockley modes Tamm ; Shockley .

In hindsight, it is natural to ask whether this approach to diagonalization may be improved and extended to more general Hamiltonians. The answer is Yes, and this paper provides the appropriate tools.

III.2 The bulk-boundary system of equations

The above motivating example suggests that it may be possible to isolate the extent to which boundary effects prevent bulk eigenstates from becoming eigenstates of the actual Hamiltonian. Consider Eqs. (7) and (8) in particular. We may condense them into a single relative eigenvalue equation, $P_{B}H_{N}|\pm k\rangle=(-2t\cos k)P_{B}|\pm k\rangle,$ in terms of the projector $P_{B}\equiv\sum_{j=2}^{N-1}|j\rangle\langle j|.$ The extension of this observation to the general class of Hamiltonians $H=H_{N}+W$ requires only knowledge of the range $R$ in Eq. (4). The block-structure of $H_{N}$ defines a subsystem decomposition of the single-particle state space abc ,

[TABLE]

where $\mathcal{H}_{L}$ and $\mathcal{H}_{I}$ are lattice and internal state spaces of dimensions $N$ and $d$, respectively. Let $\{|j\rangle,\ j=1,\dots,N\}$ and $\{|m\rangle,\ m=1,\dots,d\}$ be their respective orthonormal bases. Define bulk and boundary projectors,

[TABLE]

with $\mathds{1}\equiv\mathds{1}_{N}\otimes\mathds{1}_{d}$ the identity matrix on $\mathcal{H}$, and $\mathds{1}_{N}$, $\mathds{1}_{d}$ the identity matrices on $\mathcal{H}_{L}$ and $\mathcal{H}_{I},$ respectively (see Fig. 1). The defining property of the bulk projector is that it annihilates any boundary contribution $W$, that is, $P_{B}W=0$. Because $P_{B}+P_{\partial}=\mathds{1}$, the bulk-boundary system of equations,

[TABLE]

may be seen to be completely equivalent to the standard eigenvalue equation, $H|\epsilon\rangle=\epsilon|\epsilon\rangle$ JPA .

This bulk-boundary separation of the eigensystem problem is advantageous because the bulk equation is, in a well-defined sense, translation-invariant. Let us define a left-shift operator $T\equiv\sum_{j=1}^{N-1}|j\rangle\langle j+1|$ on the lattice space ${\cal H}_{L}$ (see Appendix B). Then, one may verify that

[TABLE]

By extending $T$ infinitely on both directions, we obtain a translation-invariant auxiliary Hamiltonian,

[TABLE]

where $\bm{T}\equiv\sum_{j\in\mathbb{Z}}|j\rangle\langle j+1|$ now denotes the generator of discrete translations on the (infinite-dimensional) vector space spanned by $\{|j\rangle\}_{j\in\mathds{Z}}$, and $\mathbf{1}$ the corresponding identity operator. The subtle difference between Hamiltonians $H_{N}$ and $\bm{H}$ is that while $T$ is not invertible, $\bm{T}$ is, and in fact $\bm{T}^{-1}=\bm{T}^{\dagger}$. This difference is decisive in solving the corresponding eigenvalue problems. On the one hand, the eigenvalue equation $\bm{H}|\Psi_{\epsilon}\rangle=\epsilon|\Psi_{\epsilon}\rangle$ is equivalent to the infinite system of linear equations

[TABLE]

where $|\Psi_{\epsilon}\rangle\equiv\sum_{j\in\mathbb{Z}}|j\rangle\otimes|\psi_{j}\rangle.$ On the other, the bulk equation $P_{B}H_{N}|\epsilon\rangle=\epsilon P_{B}|\epsilon\rangle$, with $|\epsilon\rangle\equiv\sum_{j=1}^{N}|j\rangle\otimes|\psi_{j}\rangle$ is equivalent to Eq. (15) but restricted to the finite domain $R<j\leq N-R$. Hence, the bulk equation is underdetermined (there are $2R$ more vector variables than constraints). In particular, if $|\Psi_{\epsilon}\rangle$ is an eigenstate of the infinite Hamiltonian as above, then

[TABLE]

is a solution of the bulk equation. It is in this sense of shared solutions with $\bm{H}$ that the bulk equation is, as anticipated, translation-invariant.

III.3 Exact solution of the bulk equation

Let us revisit the energy eigenvalue equation, Eq. (15). If the goal were to diagonalize the infinite-system Hamiltonian $\bm{H}$, then one should focus on finding energy eigenvectors associated to normalized states in Hilbert space. However, our model systems are of finite extent, and we are only interested in using $\bm{H}$ as an auxiliary operator for finding the translation-invariant solutions of the bulk equation. Hence, we will allow $\bm{H}$ to act on arbitrary vector sequences of the form $\Psi=\sum_{j\in\mathds{Z}}|j\rangle|\psi_{j}\rangle,$ possibly “well outside” the Hilbert state space, and so we will drop Dirac’s ket notation. From the standpoint of solving the bulk equation, every sequence that satisfies $\bm{H}\Psi=\epsilon\Psi$ is acceptable, so one must find them all. In the space of all sequences, the translation symmetry $\bm{T}$ remains invertible but is no longer unitary, because the notion of adjoint operator is not defined. This is important, because it means that translations need not have their eigenvalues on the unit circle, or be diagonalizable. Nonetheless, $[\bm{T},\bm{H}]=0$, and so both features have interesting physical consequences for finite systems.

We will refer to the space of solutions of the bulk equation as the bulk solution space and denote it by

[TABLE]

for any fixed energy $\epsilon$. Let $\mathcal{M}_{-\infty,\infty}(\epsilon)\equiv{\rm Ker}\ (\bm{H}-\epsilon{\bf 1})$ denote the space of eigenvectors of $\bm{H}$ of energy $\epsilon$ within the space of all sequences. In terms of these spaces, our arguments in Sec. III.2 establish the relation

[TABLE]

where we dropped the argument $\epsilon$. Translation invariance is equivalent to the properties $\bm{T}\mathcal{M}_{-\infty,\infty}\subseteq\mathcal{M}_{-\infty,\infty}$ and $\bm{T}^{-1}\mathcal{M}_{-\infty,\infty}\subseteq\mathcal{M}_{-\infty,\infty}$ footnote2 . If the matrix $h_{R}$ is invertible, Eq. (16) becomes $\bm{P}_{1,N}\mathcal{M}_{-\infty,\infty}=\mathcal{M}_{1,N}$ JPA .

Since $\bm{T}$ commutes with $\bm{T}^{-1}$, the generator of translations to the right, these two symmetries share eigenvectors of the form $\Phi_{z,1}|u\rangle\equiv\sum_{j\in\mathds{Z}}z^{j}|j\rangle|u\rangle,$ with $z$ an arbitrary non-zero complex number and $|u\rangle$ any internal state: there are $d$ linearly independent eigenvectors of translations for each $z\neq 0$. As a simple but important consequence of the identities

[TABLE]

one finds that

[TABLE]

where the linear operator

[TABLE]

acts on the internal space ${\cal H}_{I}$ only. This $H(z)$ is precisely the reduced bulk Hamiltonian $h_{B}(z)$ of Ref. [abc, ], obtained here by way of a slightly different argument. Since $H_{k}=H(z=e^{ik})$ is the usual Bloch Hamiltonian of a one-dimensional system with Born-von-Karman BCs, $H(z)$ is the analytic continuation of $H_{k}$ off the Brillouin zone.

One can similarly continue the energy dispersion relation off the Brillouin zone, by relating $\epsilon$ to $z$ via

[TABLE]

In practice, it is advantageous to use the polynomial

[TABLE]

We will say that $\epsilon$ is regular if $P(\epsilon,z)$ is not the zero polynomial, and singular otherwise. That is, $P(\epsilon,z)=0$ identically for all $z$ if $\epsilon$ is singular. Such a (slight) abuse of languageJPA is permitted since we are interested in varying $\epsilon$ for a fixed Hamiltonian. For any given Hamiltonian of finite range $R$, there are at most a finite number of singular energies. Physically, singular energies correspond to flat bands, as one can see by restriction to the Brillouin zone. We can now state a first useful result, whose formal proof follows from the general arguments in Ref. [JPA, ]:

Theorem 1. If $\epsilon$ is regular, the number of independent solutions of the bulk equation is $\dim\mathcal{M}_{1,N}(\epsilon)=2Rd$, for any system size $N>2R$.

This result ties well with the physical meaning of the number $2Rd=\text{dim(Range}\,P_{\partial})$ as counting the total number of degrees of freedom on the boundary, which is equal to the dimension of the boundary subspace. The condition $N>2R$ implies that the system is big enough to contain at least one site in the bulk.

III.3.1 Extended-support bulk solutions at regular energies

The solutions of the bulk equation that are inherited from $\bm{H}$ have non-vanishing support on the full lattice space ${\cal H}_{L}$, and are labeled by the eigenvalues of ${\bm{T}}$, possibly together with a second “quantum number” that appears because ${\bm{T}}$ is not unitary on the space of all sequences. For any $z\neq 0$, if $|u\rangle$ satisfies the eigenvalue equation $H(z)|u\rangle=\epsilon|u\rangle,$ then Eq. (17) implies that $\Phi_{z,1}|u\rangle$ is an eigenvector of $\bm{H}$ with eigenvalue $\epsilon$. In order to be more systematic, let $\{z_{\ell}\}_{\ell=1}^{n}$ denote the $n$ distinct non-zero roots of Eq. (19), and $\{s_{\ell}\}_{\ell=1}^{n}$ their respective multiplicities. For generic values of $\epsilon$, $H(z_{\ell})$ has exactly $s_{\ell}$ eigenvectors $\{|u_{\ell s}\rangle\}_{s=1}^{s_{\ell}}$ in ${\cal H}_{I}$, satisfying

[TABLE]

Since $\bm{H}\Phi_{z_{\ell},1}|u_{\ell s}\rangle=\epsilon\Phi_{z_{\ell},1}|u_{\ell s}\rangle,$ the states

[TABLE]

are solutions of the bulk equation. Intuitively, these states are “eigenstates of the Hamiltonian up to BCs.”

For a few isolated values of $\epsilon$, $H(z_{\ell})$ can have less than $s_{\ell}$ eigenvectors. However, the number of eigenvectors of $\bm{H}$ is still $s_{\ell}$ JPA , as we illustrate here by example. Suppose for concreteness that

[TABLE]

Since $R=1$ and $d=1$, we expect two eigenvectors for each value of $\epsilon$. One concludes that the eigenspace of energy $\epsilon$ is spanned by the sequences $\Phi_{z_{\ell},1},\ \ell=1,2,$ if $z_{1}\neq z_{2}$. But, if $\epsilon=\pm t$, then $z_{1}=z_{2}=\mp 1$, and

[TABLE]

How can one get two independent solutions in this case? The answer is that, in addition to $\Phi_{z_{1},1}$, the factor $(\bm{T}-z_{1})^{2}$ contributes another sequence to the kernel of $\bm{H}-\epsilon\bm{1}$, namely, $\Phi_{z_{1},2}=\sum_{j\in\mathds{Z}}jz_{1}^{j-1}|j\rangle.$ There are two eigenvectors in total, even though there is only one root.

Returning to the general case, the sequences JPA ; normalization

[TABLE]

span the kernel of $(\bm{T}-z)^{s}$ for $v=1,\dots,s$. In other words, $\Phi_{z,v}$ is a generalized eigenvector of the translational symmetry $\bm{T}$ of rank $v$ with eigenvalue $z$. We refer to eigenvectors with $v>1$ as the power-law solutions of the bulk equation (solutions with a power-law prefactor). They exist because translations are not diagonalizable in the full space of sequences (as opposed to the Hilbert space of square-summable sequences), leading to the new quantum number $v$.

The power-law solutions of the bulk equation may be found from the action of $\bm{H}$ on the generalized eigenvectors of $\bm{T}$. For arbitrary internal state $|u_{x}\rangle$, we have:

[TABLE]

Then one can show from Eqs. (21) and (22) that the action of $\bm{H}$ on the vector sequence $\Psi=\sum_{x=1}^{v}\Phi_{z,x}|u_{x}\rangle$, where $\{|u_{x}\rangle\}$ are arbitrary internal states, is given by

[TABLE]

Here, $H_{v}(z)$ is an upper triangular block-Toeplitz matrix with non-trivial blocks

[TABLE]

In matrix form, by letting $H^{(x)}\equiv\partial_{z}^{x}H(z)$, we have

[TABLE]

We refer to $H_{v}(z)$ as the generalized reduced bulk Hamiltonian of order $v$. Notice that $H_{1}(z)=H(z)$. In the partial basis

[TABLE]

organized as a row vector, the entries of $|u\rangle=\begin{bmatrix}|u_{1}\rangle&\dots&|u_{v}\rangle\end{bmatrix}^{\rm T}$ are the vector-valued coordinates of $\Psi$, $\Psi=\Phi_{z}|u\rangle=\sum_{x=1}^{v}\Phi_{z,x}|u_{x}\rangle.$ Then, Eq. (23) can be rewritten as

[TABLE]

Now it becomes clear that for $\Psi$ to be an eigenvector of $\bm{H}$, the required condition is $H_{v}(z)|u\rangle=\epsilon|u\rangle$, which is analogous to the condition derived for the generic case $v=1$. If a root $z_{\ell}$ of Eq. (18) has multiplicity $s_{\ell}$, then $\bm{H}$ has precisely $s_{\ell}$ linearly independent eigenvectors corresponding to $z_{\ell}$. This provides a characterization of the eigenstates of $\bm{H}$, which may be regarded as extending Bloch’s theorem to $\bm{H}$ viewed as a linear transformation on the space of all vector-valued sequences, and whose rigorous justification follows from Ref. [JPA, ]:

Theorem 2. For fixed, regular $\epsilon$, let $\{z_{\ell}\}_{\ell=1}^{n}$ denote the distinct non-zero roots of Eq. (18), with respective multiplicities $\{s_{\ell}\}_{\ell=1}^{n}$. Then, the eigenspace of $\bm{H}$ of energy $\epsilon$ is a direct sum of $n$ vector spaces spanned by generalized eigenstates of $\bm{T}$ of the form

[TABLE]

where the linearly independent vectors $\{|u_{\ell s}\rangle\}_{s=1}^{s_{\ell}}$ are chosen in such a way that $H_{s_{\ell}}(z_{\ell})|u_{\ell s}\rangle=\epsilon|u_{\ell s}\rangle$, and $|u_{\ell s}\rangle=\begin{bmatrix}|u_{\ell s1}\rangle&\dots&|u_{\ell ss_{\ell}}\rangle\end{bmatrix}^{\rm T}$.

Once the eigenvectors of $\bm{H}$ are calculated, the bulk solutions of extended support are readily obtained by projection. Let, for $v\geq 1$,

[TABLE]

be the projections of generalized eigenvectors of $\bm{T}$. Then

[TABLE]

describes a basis of the translation-invariant solutions of the bulk equation, where

[TABLE]

Remark.— The bulk equation bears power-law solutions only at a few isolated values of $\epsilon$ Trench85 . However, linear combinations of $v=1$ solutions show power-law-like behavior, as soon as two or more of the roots of Eq. (18) are sufficiently close to each other. Suppose, for instance, that for some value of energy $\epsilon$, two of the roots of Eq. (18) coincide at $z_{*}$. For energy differing from $\epsilon$ by a small amount $\delta\epsilon$, the double root $z_{*}$ bifurcates into two roots slightly away from each other, with values $z_{*}\pm\delta z$. The relevant bulk solution space is spanned by

[TABLE]

showing that the second vector has indeed a close resemblance to the power-law solution $|z_{*},2\rangle$. Similar considerations apply if $d>1$, as it is typically the case in physical applications. Assuming that the relevant bulk solutions at energy $\epsilon+\delta\epsilon$ are described by analytic vector functions $|\psi(z_{*}+\delta z)\rangle$ and $|\psi(z_{*}-\delta z)\rangle$, then, from the above analysis, it is clear that for energy $\epsilon$, the power-law bulk solution will be proportional to

[TABLE]

We will make use of this observation for the calculation of power-law solutions in Sec. V.2.

III.3.2 Emergent solutions at regular energies

While the extended solutions of the bulk equation correspond to the nonzero roots of Eq. (18), the polynomial $P(\epsilon,z)$ defined in Eq. (19) may also include $z_{0}=0$ as a root of multiplicity $s_{0}$, that is, we may generally write

[TABLE]

However, $|z=0\rangle|u\rangle=0$ does not describe any state of the system. This observation suggests that the extended solutions of the bulk equation may fail to account for all $2Rd$ solutions we expect for regular $\epsilon$. That this is indeed the case follows from a known result in the theory of matrix polynomials teran15 , implying that $2Rd=2s_{0}+\sum_{\ell=1}^{n}s_{\ell}$ for matrix polynomials associated to Hermitian Toeplitz matrices JPA . Hence, the number of solutions of the bulk equation of the form given in Eq. (26) is

[TABLE]

We call the missing $2s_{0}$ solutions of the bulk equation emergent, because they are no longer controlled by $\bm{H}$ and (nonunitary) translation symmetry, but rather they appear only because of the truncation of the infinite lattice down to a finite one, and only if $\det h_{R}=0$ JPA . Emergent solutions are a direct, albeit non-generic, manifestation of translation-symmetry-breaking; nonetheless, remarkably, they can also be determined by the analytic continuation of the Bloch Hamiltonian, in a precise sense.

While full technical detail is provided in Appendix C, the key to computing the emergent solutions is to relate the problem of solving the bulk equation to a half-infinite Hamiltonian, rather than the doubly-infinite $\bm{H}$ we have exploited thus far. Let us define the unilateral shifts

[TABLE]

The Hamiltonian

[TABLE]

is then the half-infinite counterpart of $\bm{H}$. The corresponding half-infinite bulk projector is

[TABLE]

Suppose there is a state $\Upsilon^{-}$, that solves the equation $\bm{P}^{-}_{B}(\bm{H}_{-}-\epsilon\bm{1}_{-})\Upsilon^{-}=0.$ Then one can check that $|\psi\rangle=\bm{P}_{1,N}\Upsilon^{-}$ is a solution of the bulk equation, Eq. (12). Clearly, some of the bulk solutions we arrive at in this way using $\bm{H}_{-}$ will coincide with those obtained from $\bm{H}$. These are precisely the extended solutions we already computed in Sec. III.3.1. In contrast, the emergent solutions are obtained only from $\bm{H}_{-}$.

Since $\bm{T}_{-}\bm{T}_{-}^{\star}=\bm{1}_{-}$, we may write $\bm{P}^{-}_{B}(\bm{H}_{-}-\epsilon\bm{1}_{-})=\bm{T}_{-}^{\star R}K^{-}(\epsilon,\bm{T}_{-}),$ in terms of the matrix polynomial

[TABLE]

Half of the emergent solutions, namely, the ones localized on the left edge, are determined by the kernel of $K_{s_{0}}^{-}(\epsilon,z_{0}=0)\equiv K^{-}(\epsilon)$, with $[K^{-}_{v}(\epsilon,z)]_{xx^{\prime}}$ constructed as in Eq. (24). Explicitly, such a matrix, which was obtained by different means in Ref. [JPA, ], takes the form

[TABLE]

for systems with fairly large $s_{0}>2R+1$. Let $\{|u^{-}_{s}\rangle\}_{s=1}^{s_{0}}$ denote a basis of the kernel of $K^{-}(\epsilon)$, with

[TABLE]

Then,

[TABLE]

are the emergent solutions with support on the first $s_{0}$ lattice sites, with $s_{0}$ obeying Eq. (28).

We are still missing $s_{0}$ emergent solutions for the right edge. They may be constructed from the kernel of the lower-triangular block matrix $K^{+}(\epsilon)\equiv[K^{-}(\epsilon)]^{\dagger}=[K^{-}_{s_{0}}(\epsilon,z_{0}=0)]^{\dagger}.$ Let $\{|u^{+}_{s}\rangle\}_{s=1}^{s_{0}}$ denote a basis of the kernel of $K^{+}(\epsilon)$, with

[TABLE]

Then,

[TABLE]

are the emergent bulk solutions associated to the right edge, supported on the lattice sites $N-s_{0}+1,\ldots,N$. Again, for mathematical justifications, see Appendix  C.

In what follows, we shall denote the spaces spanned by left- and right- localized emergent bulk solutions by $\mathcal{F}_{1}^{-}$ and $\mathcal{F}_{N}^{+}$, and their bases by $\mathcal{B}^{-}\equiv\{|\psi_{s}^{-}\rangle\}_{s=1}^{s_{0}}$ and $\mathcal{B}^{+}\equiv\{|\psi_{s}^{+}\rangle\}_{s=1}^{s_{0}}$, respectively.

III.3.3 Bulk-localized states at singular energies

If $h_{R}$ is not invertible, there can be at most a finite number of singular energy values (usually referred to as flat bands), leading to bulk-localized solutions: these solutions are finitely-supported and appear everywhere in the bulk. Hence, a singular energy cannot be excluded from the physical spectrum of a finite system by way of BCs. In contrast, emergent solutions are also finitely-supported but necessarily “anchored” to the edges (and only appearing for regular values of $\epsilon$).

Recall that if $\epsilon$ is singular, then $\det(H(z)-\epsilon\mathds{1})=0$ for any $z$. Thus, there exists an analytic vector function,

[TABLE]

satisfying $H(z)|v(z)\rangle=\epsilon|v(z)\rangle$ for all $z$. To obtain $|v(z)\rangle$, one can construct the adjugate matrix of $(H(z)-\epsilon\mathds{1}_{d})$. (Recall that the adjugate matrix $\text{adj}(M)$ associated to a square matrix $M$ is constructed out of the signed minors of $M$ and satisfies $\text{adj}(M)M=\det(M)\mathds{1}$.) Hence,

[TABLE]

and so one can use any of the non-zero columns of $\text{adj}(H(z)-\epsilon\mathds{1}_{d})$, suitably pre-multiplied by a power of $z$, for the vector polynomial $|v(z)\rangle$. By matching powers of $z$, this equation becomes

[TABLE]

The idea now is to use the linearly independent solutions of Eq. (35) to construct finite-support solutions of the bulk equation. Let us denote such solutions by $|v_{\mu}\rangle\equiv\begin{bmatrix}|v_{\mu 0}\rangle&|v_{\mu 1}\rangle&\dots&|v_{\mu\delta_{0}}\rangle\end{bmatrix}^{\rm T}$, for $\mu=1,\dots,\mu_{0}$. One can check directly that the finitely-supported sequences

[TABLE]

all satisfy $(\bm{H}-\epsilon\bm{1})\Psi_{js}=0$ because $|v_{\mu}\rangle$ obeys Eq. (35). Hence, the states $\bm{P}_{1,N}\Psi_{j\mu}$ provide finitely-supported solutions of the bulk equation. In addition, as long as $2R<j<N-2R-\delta_{0}$, the boundary equation is also satisfied trivially, and so all such states become eigenvectors of $H_{N}+W$ with the singular energy $\epsilon$. This is why singular energies, if present for the infinite system, are necessarily also part of the spectrum of the finite system and display macroscopic degeneracy of order $\mathcal{O}(N)$.

Let us further remark that the sequences $\Psi_{j\mu}$ and associated solutions of the bulk equation need not be linearly independent. To obtain a complete (rather than overcomplete), set of solutions for flat bands, one would require a technical tool, the Smith normal form Gohberg , which is beyond the scope of this paper. See Ref. [JPA, ] for details.

III.4 The boundary matrix

For regular energies, the bulk solutions determine a subspace of the full Hilbert space [Theorem 1], whose dimension $2Rd\ll dN$ for typical applications. While not all bulk solutions are eigenstates of the Hamiltonian $H=H_{N}+W$, the actual eigenstates must necessarily appear as bulk solutions. Hence, the bulk-boundary separation in Eqs. (12), and, in particular, the bulk equation, identifies by way of a translational symmetry analysis a small search subspace. In order to find the energy eigenstates efficiently, one must solve the boundary equation on this search subspace. Since the boundary equation is linear, its restriction to the space of bulk solutions can be represented by a matrix, the boundary matrix abc . The latter is a square matrix that combines our basis of bulk solutions with the relevant BCs.

Let ${\cal B}\equiv{\cal B}_{\text{ext}}\cup\mathcal{B}^{-}\cup\mathcal{B}^{+}$ be a basis for $\mathcal{M}_{1,N}$. Then, building on the previous section, the Ansatz state

[TABLE]

represents the solutions of the bulk equation parametrized by the $2Rd$ amplitudes $\bm{\alpha}$, where

[TABLE]

Moreover, let as before $b=1,\dots,R,N-R+1,\dots,N$ label the boundary sites. Then,

[TABLE]

In particular, the boundary equation is equivalent to the requirement that $\langle b|(H_{N}+W-\epsilon\mathds{1})|\Psi_{\mathcal{B}}\rangle\bm{\alpha}=0$ for all boundary sites. Since $\langle b|(H_{N}+W-\epsilon\mathds{1})|\Psi_{\mathcal{B}}\rangle\equiv\langle b|H_{\epsilon}|\Psi_{\mathcal{B}}\rangle$ denotes a row array of internal states, it is possible to organize these arrays into the boundary matrix

[TABLE]

By construction, the boundary matrix $B$ is a block matrix of block-size $d\times 1$. In terms of this matrix, Eq. (40) provides the useful identity

[TABLE]

One may write an analogous equation in Fock space by defining an array

[TABLE]

Then Eq. (42) translates into

[TABLE]

It is interesting to notice that this (many-body) relation remains true even if $\epsilon$ is allowed to be a complex number.

III.5 The generalized Bloch theorem

The bulk-boundary separation of the energy eigenvalue equation shows that actual energy eigenstates are necessarily linear combinations of solutions of the bulk equation. This observation leads to a generalization of Bloch’s theorem for independent fermions under arbitrary BCs:

Theorem 3 (Generalized Bloch theorem). Let $H=H_{N}+W$ denote the single-particle Hamiltonian of a clean system subject to BCs described by $W=P_{\partial}W$. If $\epsilon$ is a regular energy eigenvalue of $H$ of degeneracy $\cal{K}$, the associated eigenstates can be taken to be of the form

[TABLE]

where $\{\bm{\alpha}_{\kappa},\ \kappa=1,\dots,\cal{K}\}$ is a basis of the kernel of the boundary matrix $B(\epsilon)$ at energy $\epsilon$.

In short, $(H_{N}+W)|\epsilon,\bm{\alpha}\rangle=\epsilon|\epsilon,\bm{\alpha}\rangle$ if and only if $B\bm{\alpha}=0$, in which case it also follows from Eq. (43) that $\eta_{\epsilon,\bm{\alpha}}^{\dagger}$ is a normal fermionic mode of the many-body Hamiltonian $\widehat{H}$. From now on, we will refer to energy eigenstates of the form $|\Psi_{\mathcal{B}}\rangle\bm{\alpha}_{\kappa}$ as generalized Bloch states. Recall that $H$ acts on ${\cal H}={\mathbb{C}}^{N}\otimes{\mathbb{C}}^{d}$, with couplings of finite range $R$. A lower bound on $N$ should be obeyed, in order for the above theorem to apply. If $\det h_{R}\neq 0$, since there are no emergent solutions nor flat bands, generalized Bloch states describe the allowed energy eigenstates as soon as $N>2R$, independently of $d$. If $h_{R}$ fails to be invertible, we should require that $N>2\,{\rm max}(s_{0},R)$ to ensure that emergent solutions on opposite edges do not overlap, and are thus independent. Since $s_{0}\leq Rd$, this condition is satisfied for any $N>2Rd$. In general, $N>2R(d+1)$ always suffices for generalized Bloch states to describe generic energy eigenstates JPA .

We further note that if $\epsilon$ is not an energy eigenvalue, the kernel of $B(\epsilon)$ is trivial. Thus, the degeneracy of a single-particle energy level coincides with the dimension of the kernel of $B(\epsilon)$. Let $\rho(\omega)$ denote the single-particle density of states. Combining its definition with the generalized Bloch theorem, we then see that

[TABLE]

an alternative formula to the usual

[TABLE]

from the theory of Green’s functions ballentine . Another interesting and closely related formula is

[TABLE]

for the partition function of the single-particle Hamiltonian, with the dependence on BCs highlighted quelle15 .

We conclude this section by showing how, for periodic BCs, one consistently recovers the conventional Bloch’s theorem. In this case, the appropriate matrix $W$ reads

[TABLE]

since then one can check that

[TABLE]

in terms of the fundamental circulant matrix

[TABLE]

Physically, $V$ is the generator of translations (to the left) for a system displaying ring (1-torus) topology.

The Bloch states are the states that diagonalize $H_{p}$ and $V$ simultaneously. Theorem 3 guarantees that we can choose the eigenstates of $H_{p}$ to be linear combinations of translation-invariant and emergent solutions. Thus, we only need to check if these linear combinations include eigenstates of $V$. There is no hope of retaining the emergent solutions, because they are localized and too few in number (at most $2Rd$) to be rearranged into eigenstates of $V$. The same holds for translation-invariant solutions with a power-law prefactor. Hence, the search subspace that is compatible with the translational symmetry $V$ is described by the simplified Ansatz abc

[TABLE]

Now, $V|\psi_{\ell 1}\rangle=z_{\ell}|\psi_{\ell 1}\rangle-z_{\ell}(1-z_{\ell}^{N})|N\rangle|u_{\ell s_{\ell}1}\rangle,$ and so the generalized Bloch states can only be eigenstates of $V$ if $e^{ik_{\ell}N}=1$ with $z_{\ell}=e^{ik_{\ell}}$, and all but one entry in $\bm{\alpha}$ vanish. That is, $|\epsilon,\bm{\alpha}\rangle\equiv|\epsilon,k_{\ell}\rangle=|z_{\ell},1\rangle|u_{\ell 1,1}\rangle.$ As one may verify, $H_{p}|\epsilon,k_{\ell}\rangle=|z_{\ell},1\rangle H(z_{\ell})|u_{\ell 1,1}\rangle=\epsilon|\epsilon,k_{\ell}\rangle,$ showing that $|\epsilon,k_{\ell}\rangle$ is indeed compatible with the boundary matrix. Manifestly, $|\epsilon,k_{\ell}\rangle$ is an eingenstate of $H_{p}$ in the standard Bloch form – thereby recovering the conventional Bloch’s theorem for periodic BCs, as desired.

IV The bulk-boundary algorithms

The results of Sec. III can be used to develop diagonalization algorithms for the relevant class of single-particle Hamiltonians. We will describe two such algorithms. The first treats $\epsilon$ as a parameter for numerical search. The second is inspired by the algebraic Bethe Ansatz, as suggested by comparing our Eq. (42) to Eq. (28) of Ref. [ortiz05, ].

IV.1 Numerical “scan-in-energy” diagonalization

The procedure described in this section is a special instance of the Eigensystem Algorithm described in Ref. [JPA, ], specialized to Hermitian matrices. It employs a search for energy eigenvalues along the real line, and takes advantage of the results of Sec. III to determine whether a given number is an eigenvalue. The overall procedure is schematically depicted in Fig. 2.

The first part of the algorithm finds all eigenvectors of $H$ that correspond to the flat (dispersionless) energy band, if any exists. Two steps are entailed:

Find all real values of $\epsilon$ for which $\det(H(z)-\epsilon\mathds{1}_{d})$ vanishes for any $z$. Output these as singular eigenvalues of $H$. 2. 2.

For each of the eigenvalues found in step (1), find and output a basis of the corresponding eigenspace of $H$ using any conventional algorithm.

In implementing step (2) above, one can leverage the analysis of Sec. III.3.3. The following part of the algorithm, which repeats until all eigenvectors of $H$ are found, proceeds according to the following steps:.

Choose a seed value of $\epsilon$, different from those eigenvalues found already. 2. 4.

Find all $n$ distinct non-zero roots of the equation $\det(H(z)-\epsilon\mathds{1}_{d})=0$. Let these roots be $\{z_{\ell},\ \ell=1,\dots,n\}$, and their respective multiplicities $\{s_{\ell},\ \ell=1,\dots,n\}$. 3. 5.

For each such roots, construct the generalized reduced bulk Hamiltonian $H_{s_{\ell}}(z_{\ell})$ [Eq. (24)]. 4. 6.

Find a basis of the eigenspace of $H_{s_{\ell}}(z_{\ell})$ with eigenvalue $\epsilon$. Let the basis vectors be $\{|u_{\ell s}\rangle,\ s=1,\dots,s_{\ell}\}$. The bulk solution corresponding to $(\ell,s)$ is $|\psi_{\ell s}\rangle=|z_{\ell},1\rangle|u_{\ell s}\rangle$, with $\Phi_{z_{\ell}}$ defined in Eq. (25). 5. 7.

If $h_{R}$ is non-invertible, find $s_{0}=Rd-\sum_{\ell=1}^{n}s_{\ell}/2$. Construct matrices $K^{-}(\epsilon)$ as described in Eq. (III.3.2), and $K^{+}(\epsilon)=[{K^{-}}(\epsilon)]^{\dagger}$. 6. 8.

Find bases of the kernels of $K^{-}(\epsilon)$ and $K^{+}(\epsilon)$. Let the basis vectors be $\{|u_{s}^{-}\rangle,\ s=1,\dots,s_{0}\}$ and $\{|u_{s}^{+}\rangle,\ s=1,\dots,s_{0}\}$, respectively. The emergent bulk solutions corresponding to each $s$ are follow from Eqs. (32) and (33). 7. 9.

Construct the boundary matrix $B(\epsilon)$ [Eq. (41)]. 8. 10.

If $\det B(\epsilon)=0$, output $\epsilon$ as an eigenvalue. Find a basis $\{\bm{\alpha}_{\kappa},\ \kappa=1,\dots,{\cal K}\}$ of the kernel of $B(\epsilon)$. Then a basis of the eigenspace of $H$ corresponding to energy $\epsilon$ is $\{|\epsilon_{\kappa}\rangle=|\Psi_{\mathcal{B}}\rangle\bm{\alpha}_{\kappa},\ \kappa=1,\dots,{\cal K}\}$, with $|\Psi_{\mathcal{B}}\rangle$ being defined in Eqs. (39). If all $2dN$ eigenvectors are not yet found, then go back to step (3). 9. 11.

If $\det B(\epsilon)\neq 0$, choose a new value of $\epsilon$ as dictated by the relevant root-finding algorithm NoteRoot . Go back to step (4).

Some considerations are in order, in regard to the fact that the determinant of $B(\epsilon)$ plotted as a function of energy $\epsilon$ may display finite-precision inaccuracies, that appear as fictitious roots. Such issues arise at those $\epsilon$ where two (or more) of the roots of Eq. (18) cross as a function of $\epsilon$, due to the non-orthogonality of the basis $\mathcal{B}$ that results from the procedure described in Sec. III.3. Let $\epsilon_{*}$ be a value of energy for which this happens, so that the bulk equation bears a power-law solution. For $\epsilon\approx\epsilon_{*}$ (except $\epsilon_{*}$ itself), Eq. (18) has two roots that are very close in value, so that the corresponding bulk solutions overlap almost completely. This results in a boundary matrix having two nearly identical columns, with determinant vanishing in the limit $\epsilon\rightarrow\epsilon_{*}$, irrespective of $\epsilon_{*}$ being an eigenvalue of $H$ (hence, a physical solution). However, if we calculate $B(\epsilon)$ exactly at $\epsilon_{*}$, then the basis $\mathcal{B}$ contains power-law solutions, and accurately indicates whether $\epsilon_{*}$ is an eigenvalue. This also means that the function $\det B(\epsilon)$ has a discontinuity at $\epsilon=\epsilon_{*}$.

A simple way to identify those fictitious roots is as follows. Rewrite the polynomial in Eq. (19) as

[TABLE]

which is treated as a polynomial in $z$ with coefficients depending on $\epsilon$ (if $s_{0}$ changes with $\epsilon$, we use the smallest possible value of $s_{0}$ in Eq. (44)). $P(\epsilon,z)$ has double roots at $\epsilon_{*}$ if and only if the discriminant $D(P(\epsilon_{*},z))=0$ Gelfand08 . The latter gives a polynomial expression in $\epsilon$, of degree $\mathcal{O}(dR)$. By finding the roots of this equation, one can obtain all the values of $\epsilon$ for which fictitious roots of $\det B(\epsilon)$ may appear. To check whether these roots are true eigenvalues, one then needs to construct $B(\epsilon)$ by including the power-law solutions in the Ansatz.

We further note that, while the Ansatz is not continuous at such values of $\epsilon$, the fact that the bulk solution space is the kernel of the linear operator $P_{B}(H_{N}+W-\epsilon)$ implies that it must change smoothly with $\epsilon$. A way to improve numerical accuracy would be to construct an orthonormal basis (e.g., via Gram-Schmidt orthogonalization) of ${\cal M}_{1,N}(\epsilon)$ at each $\epsilon$, and use this basis to construct a modified boundary matrix $\tilde{B}(\epsilon)$. In practice, one may directly compute the new determinant by using

[TABLE]

where $\mathcal{G}\equiv\langle\Psi_{\mathcal{B}}|\Psi_{\mathcal{B}}\rangle$ is the Gramian matrix grammatrix of the basis of bulk solutions obtained in steps (4) to (8) of the algorithm, with entries $\mathcal{G}_{ss^{\prime}}\equiv\langle\psi_{s}|\psi_{s^{\prime}}\rangle,$ $s,s^{\prime}=1,\dots,2Rd.$ In fact, it can be checked that the bulk solutions

[TABLE]

form an orthonormal basis of the bulk solution space ${\mathcal{M}}_{1,N}$. The calculation of the entries of the Gramian is straightforward thanks to the analytic result

[TABLE]

In regard to the time and space complexity of the algorithm, the required resources depend entirely on those needed to compute the boundary matrix. For generic $\epsilon$, regardless of the invertibility of $h_{R}$, the size of $B(\epsilon)$ is $2Rd\times 2Rd$, independently of $N$. Calculation of each of its entries is also simple from the point of view of complexity, thanks to the fact that $H=H_{N}+W$ is symmetrical RemarkSymm ; JPA . Accordingly, both the number of steps and the memory space used by this algorithm do not scale with the system size $N$, making this approach computationally more efficient than conventional methods of diagonalization of generic Hermitian matrices Demmel07 .

IV.2 Algebraic diagonalization

The scan-in-energy algorithm can be further developed into an algorithm that yields an analytic solution (often closed-form), in the same sense as the Bethe Ansatz method does for a different class of (interacting) quantum integrable systems. The idea is to obtain, for generic values of $\epsilon$, an analytic expression for $B(\epsilon)$, since its determinant will then provide a condition for $\epsilon$ to be an eigenvalue, and the corresponding eigenvectors can be obtained from its kernel. As mentioned, for generic $\epsilon$, the extended bulk solutions do not include any power-law solutions. This property can be exploited to derive an analytic expression for $B(\epsilon)$ in such a generic setting. The values of $\epsilon$ for which power-law solutions appear, or the analytic expression fails for other reasons, can be dealt with on a case-by-case basis.

By the Abel-Ruffini theorem, a completely closed-form solution by radicals in terms of $\epsilon$ can be achieved if the degree in $z$ of the characteristic polynomial of the reduced bulk Hamiltonian is at most four. If this is not the case, the roots $\{z_{\ell}\}$ do not possess an algebraic expression in terms of $\epsilon$ and entries of $H$. The workaround is then to consider $\{z_{\ell}\}$ as free variables, with the constraint that each of them satisfy the characteristic equation of $H(z)$. With these tools in hand, the following procedure can be used to find an analytical solution for generic values of $\epsilon$:

Construct the polynomial $P(\epsilon,z)$ in Eq. (44), which is a bivariate polynomial in $\epsilon$ and $z$. Determine $s_{0}$ using $s_{0}=2Rd-\text{deg}(P(\epsilon,z)),$ where $\text{deg}(.)$ denotes the degree of the polynomial in $z$. 2. 2.

Assuming that $\epsilon$ and $z$ satisfy $P(\epsilon,z)=0$, find an expression for the eigenvector $|u(\epsilon,z)\rangle$ of $H(z)$ with eigenvalue $\epsilon$. 3. 3.

Consider variables $\{z_{\ell},\ \ell=1,\dots,2Rd-2s_{0}\}$, each satisfying $P_{\epsilon}(z_{\ell})=0$. Each of these corresponds to a bulk solution $|z_{\ell},1\rangle|u(\epsilon,z_{\ell})\rangle$. 4. 4.

If $h_{R}$ is not invertible, construct matrices $K^{-}(\epsilon)$ and $K^{+}(\epsilon)=[K^{-}(\epsilon)^{\dagger}]$ [Eq. (III.3.2)]. 5. 5.

Find bases for their kernels, each of which contains $s_{0}$ vectors. Let these be $\{|u_{s}^{-}(\epsilon)\rangle,\ s=1,\dots,s_{0}\}$ and $\{|u_{s}^{+}(\epsilon)\rangle,\ s=1,\dots,s_{0}\}$. These correspond to finite-support solutions of the bulk equation. 6. 6.

Construct the boundary matrix $B(\epsilon)\equiv B(\epsilon,\{z_{\ell}\})$ [Eq. (41)]. 7. 7.

The condition for $\epsilon$ being an eigenvalue of $H$ is $\det B(\epsilon,\{z_{\ell}\})=0$. Therefore, a complete characterization of eigenvalues is

[TABLE] 8. 8.

If $\text{deg}(P(\epsilon,z))\leq 4$, substitute for each $z_{\ell}$ the closed-form expression of the corresponding root $z_{\ell}(\epsilon)$. The eigenvalue condition in step (7) simplifies to a single equation, $\det B(\epsilon,\{z_{\ell}(\epsilon)\})=0.$ 9. 9.

For every eigenvalue $\epsilon$, the kernel vector $\bm{\alpha}(\epsilon,\{z_{\ell}\})$ of $B(\epsilon,\{z_{\ell}\})$ provides the corresponding eigenvector of $H$.

In steps (2), (5) and (9), we need to obtain an analytic expression for the basis of the kernel of a square symbolic matrix of fixed kernel dimension in terms of its entries. This can be done in many different ways, and often is possible by inspection. One possible way was described in Sec. III.3.3 in connection to evaluating Ker($H(z)-\epsilon\mathds{1}_{d})$ for singular values of $\epsilon$. The above analysis does not hold when $\epsilon$ satisfies any of the following conditions:

1. (i)

$\det(H(z)-\epsilon\mathds{1})=0$ has one or more double roots. This is equivalent to $D(P(\epsilon,z))=0$, as discussed in Sec. IV.1. This is a polynomial equation in terms of $\epsilon$, the roots of which yield all required values of $\epsilon$. 2. (ii)

The coefficient $p_{s_{0}}(\epsilon)$ of $z^{s_{0}}$ in $P(\epsilon,z)$ vanishes, or equivalently, $\epsilon$ is a root of $p_{s_{0}}(\epsilon)=0$. 3. (iii)

Each entry of $|u(\epsilon,z)\rangle$ vanishes. Such points are identified by solving simultaneously the equations $\langle m|u(\epsilon,z)\rangle=0,\ m=1,\dots,d$ and $P(\epsilon,z)=0$, Since a necessary and sufficient condition for these polynomials (in $z$) to have a common root is that their resultant vanishes Gelfand08 , we find the relevant values of $\epsilon$ by equating the pairwise resultants to zero. 4. (iv)

$\{|u_{s}^{-}(\epsilon)\rangle,\ s=1,\dots,s_{0}\}$ or $\{|u_{s}^{+}(\epsilon)\rangle,\ s=1,\dots,s_{0}\}$ are linearly dependent. To find such values of $\epsilon$, one may form the corresponding Gramian matrix and equate its determinant to zero.

For all the values of $\epsilon$ thus identified, $B(\epsilon)$ is calculated by following steps (4)-(10) in the scan-in-energy algorithm. To summarize, this algebraic procedure achieves diagonalization in analytic form: the upshot is a system of polynomial equations, whose simultaneous roots are the eigenvalues, and an analytic expression for the eigenvectors, with parametric dependence on the eigenvalue.

V Illustrative examples

This section contains three paradigmatic examples illustrating the use of our generalized Bloch theorem, along with the resulting algebraic procedure of diagonalization.

V.1 The impurity model revisited

Let us first reconsider the impurity model of Sec. III.1. The single-particle Hamiltonian is the corner-modified, banded block-Toeplitz matrix $H=H_{N}+W$, with

[TABLE]

The boundary consists of two sites, so that $P_{\partial}=|1\rangle\langle 1|+|N\rangle\langle N|$, for any $N>2$. Likewise, $R=1=d$. The first step in diagonalizing $H$ is solving the bulk equation. Since the reduced bulk Hamiltonian $H(z)=-t(z+z^{-1})$,

[TABLE]

Thus, every value of $\epsilon$ is regular and yields two ($=$ the number of boundary degrees of freedom) solutions of the bulk equation. If $\epsilon\neq\pm 2t$, the solutions are $|z_{\ell},1\rangle$, with

[TABLE]

with $z_{1}z_{2}=1$ and $\epsilon=-t(z_{1}+z_{2})$. The special values $\epsilon=\pm 2t$ for which $H_{N}$ yields only one of the two bulk solution have an interpretation as the edges of the energy band. If $\epsilon=2t$, then $H(z)$ yields only $|z_{1}=-1,1\rangle$, whereas if $\epsilon=-2t$, it yields only $|z_{1}=1,1\rangle$. In order to obtain the missing bulk solution in each case, one must consider the effective Hamiltonian [Eq. (24)]

[TABLE]

One may check that $H_{2}(z_{1})-\epsilon\mathds{1}\equiv 0$ if $\epsilon=\pm 2t,\ z_{1}=\mp 1$. Thus, the two linearly independent solutions of the bulk equation at these energies are $|z_{1}=1,v\rangle,$ $v=1,2$, if $\epsilon=-2t$, and $|z_{1}=-1,v\rangle,$ $v=1,2$, if $\epsilon=2t$.

For the purpose of solving the boundary equation, and hence the full diagonalization problem, it is convenient to organize the solutions of the bulk equation as

[TABLE]

For comparison with Sec. III.1, one should think of $z_{1}=e^{ik}$ and $z_{2}=e^{-ik}$. Because the Ansatz is naturally broken into three pieces, so is the boundary matrix. For instance, when $\epsilon\neq\pm 2t$, direct calculation yields

[TABLE]

However, from Eq. (45) it follows that

[TABLE]

This allows a simpler form to be obtained, by effectively changing the argument of the boundary matrix from $\epsilon$ to $z_{\ell}$ (or $k$). The complete final expression reads:

[TABLE]

Notice that if $\epsilon$ approaches $\pm 2t$, the two distinct roots collide at $z_{1}=z_{2}=\mp 1$, and the boundary matrix becomes, trivially, a rank-one matrix, signaling the discontinuous behavior anticipated in Sec. IV.1. Furthermore, it follows from Eq. (27) that the power-law solution at $\epsilon=\pm 2t$ may be written as $\partial_{z}(|z_{1},1\rangle)=|z_{1},2\rangle$. The entries of the second column of the corresponding boundary matrices satisfy $\langle b|H_{\epsilon}|z_{1},2\rangle=\partial_{z_{2}}\langle b|H_{\epsilon}|z_{2},1\rangle|_{z_{2}=z_{1}}$, where $z_{1}$ is the double root. Thus, the entries in the second column of the boundary matrix for $\epsilon=\pm 2t$ can be obtained by differentiating with respect to $z_{2}$ the second column of the boundary matrix for other (generic) values of $\epsilon$, an observation we will use in other examples as well (see e.g. Sec. V.2.2). We now analyze separately different regimes (see also Fig. 3 for illustration).

V.1.1 Vanishing impurity potential

If $w=0$, then $B(\epsilon=2t)$ and $B(\epsilon=-2t)$ have a trivial kernel; the exotic states $|\epsilon=\pm 2t\rangle$ cannot possibly arise as physical eigenvectors. For other energies, we find that the kernel of the boundary matrix

[TABLE]

is nontrivial only if $z_{1}^{N+1}=z_{2}^{N+1}$, in which case we can take $\alpha_{1}=1$ and $\alpha_{2}=-1$. From Eq. (45), it also follows that $z_{1}z_{2}=1$. Hence, there are $2N+2$ solutions,

[TABLE]

Of the associated $2N+2$ (un-normalized) Ansatz vectors

[TABLE]

two vanish identically ($q=-N-1$ and $q=0$). For $q=\pm 1,\dots,\pm N$, it is immediate to check that $|\epsilon_{-q}\rangle=-|\epsilon_{q}\rangle$. This means that the Ansatz yields exactly $N$ linearly independent energy eigenvectors, of energy

[TABLE]

This is precisely the result of Sec. III.1, where the solutions were labelled in terms of allowed quantum numbers $k=\pi q/(N+1),\ q=1,\dots,N$.

According to our general theory, the eigenspaces of $H$ are in one-to-one correspondence with the zeroes of $\det B(\epsilon)$. For this system then, there should be at most $N$ zeroes. The reason we find $2N+2$ zeroes is due to the above-mentioned (quadratic) change of argument in the boundary matrix from $\epsilon$ to $k$. Such a change of variables is advantageous for analytic work, and the associated redundancy is always rectified at the level of the Ansatz.

V.1.2 Power-law solutions

What would it take for $|\epsilon=\pm 2t\rangle$ to become eigenvectors? The kernel of $B(\epsilon=2t)$ is nontrivial only if

[TABLE]

These two values coincide up to corrections of order $1/N$, but remember that our analysis is exact for any $N>2$. Similarly, the kernel of $B(\epsilon=-2t)$ is nontrivial only if

[TABLE]

Only one of these conditions can be met: for fixed $w$, either $|\epsilon=2t\rangle$ is an energy eigenstate or $|\epsilon=-2t\rangle$ is, but not both. Let us look more closely at the state at the bottom of the energy band. As we just noticed, this state will be a valid eigenstate for either of the two values of $w$. Let us pick $w\equiv w_{N}=-t({N+1})/({N-1})$, since it yields the most interesting ground state. Then,

[TABLE]

so that one can set $\alpha_{1}=1/(w_{N}+t)$, $\alpha_{2}=-1/w_{N}$, and

[TABLE]

Notice that $\langle j|\epsilon=-2t\rangle=-\langle N-j+1|\epsilon=-2t\rangle$; that is, the power-law eigenvector of the impurity problem is an eigenstate of inversion symmetry.

V.1.3 Strong impurity potential

Lastly, consider the regime where $t\ll|w|$, for large $N$. Then, the values $\epsilon=\pm 2t$ are excluded from the physical spectrum, and the eigenstates of the system can be determined from $\det B(\epsilon)=0$. We expect bound states of energy $w$ to leading order and well-localized at the edges, so that $0<|z_{1}|<1<|z_{2}|$, say, with $z_{1}$ ($z_{2}$) associated to the left (right) edge. It is convenient to take advantage of this feature and modify the original Ansatz to

[TABLE]

so that $|z_{1},1\rangle$ ($z_{2}^{-N}|z_{2},1\rangle$) peaks at the left (right) edge, respectively. The boundary matrix becomes

[TABLE]

since $|z_{1}|^{N}\approx 0\approx|z_{2}|^{-N}$. Keeping in mind that $z_{1}z_{2}=1$, we see that the kernel of $\tilde{B}(\epsilon)$ is two-dimensional for

[TABLE]

and otherwise trivial. The corresponding energy eigenstates can be chosen to be

[TABLE]

Notice that $|\epsilon_{b},2\rangle$ is the mirror image of $|\epsilon_{b},1\rangle$, up to normalization. The large-$N$ approach to boundary modes exemplified by the preceding calculation can be made systematic, as we will further explain in Sec. VI.1.

The remaining $(N-2)$ eigenstates consist of standing waves. They can be computed from the original boundary matrix, approximated for $t\ll|w|$ as

[TABLE]

This boundary matrix has a nontrivial kernel only if

[TABLE]

in which case one may choose $\alpha_{1}=z_{2},\ \alpha_{2}=-z_{1}$. Then,

[TABLE]

Moreover, $|\epsilon_{s}\rangle=-|\epsilon_{N-1+s}\rangle,$ $s=1,\dots,N-2.$ Hence, as needed, we have obtained $(N-2)$ linearly independent eigenvectors of energy $\epsilon_{s}=-2t\cos[\pi s/(N-1)]$.

The above discussion is further illustrated in Fig. 3, where the determinant of the exact boundary matrix is displayed as a function of energy.

V.2 Engineering perfectly localized

zero-energy modes: A periodic Anderson model

Having illustrated the algebraic diagonalization method on a simple impurity model, we illustrate next its usefulness toward Hamiltonian engineering. In this section, we will design from basic principles a “comb” model, see Fig. 4, with the peculiar property of exhibiting a perfectly localized mode at zero energy while all other modes are dispersive. The zero mode is distributed over two sites on the same end of the comb, with weights determined by a ratio of hopping amplitudes.

The starting point is the single-particle Hamiltonian

[TABLE]

In order to have perfectly localized eigenvectors at zero energy, the bulk equation must bear emergent solutions. Therefore, we assume that $h_{1}$ is non-invertible. Let $|u^{-}\rangle$ be in the kernel of $h_{1}^{\dagger}$. Since $T$ annihilates $|j=1\rangle$,

[TABLE]

Similarly, if $|u^{+}\rangle$ is in the kernel of $h_{1}$, then $|j=N\rangle|u^{+}\rangle$ is also in the kernel of $H$. Therefore, $|j=1\rangle|u^{-}\rangle$ and $|j=N\rangle|u^{+}\rangle$ are perfectly localized zero energy modes.

A concrete example may be obtained by choosing

[TABLE]

whose kernel is spanned by

[TABLE]

respectively. This example corresponds to a many-body Hamiltonian of two coupled fermionic chains, as illustrated in Fig. 4:

[TABLE]

where $c_{j}$ and $f_{j}$ denote the $j$th fermions in the upper and lower chain, $t_{0}$ denotes intra-ladder hopping in one of the chains, and $t_{1}$ is the diagonal hopping strength between the two chains of the ladder, respectively. Physically, this “topological comb model” is closely related to the one-dimensional periodic Anderson model in its non-interacting (spinless) limit, see Ref. [Dagotto, ].

V.2.1 Zero-energy modes

The perfectly localized zero-energy modes in this case are $|j=1\rangle|u^{-}\rangle$ and $|j=N\rangle|u^{+}\rangle$, that translate, after normalization, into the fermionic operators

[TABLE]

The operator $\eta_{2}^{\dagger}$ trivially describes a zero-energy mode, since it corresponds to the last fermion on the lower chain, that is decoupled from the rest. However, $\eta_{1}^{\dagger}$ corresponds to a non-trivial zero energy mode, localized over the first sites of the two chains. For large values of $|t_{0}/t_{1}|$, $\eta_{1}^{\dagger}$ is localized mostly on the $f$-chain, whereas for small values it is localized mostly on the $c$-chain.

Remarkably, such a non-trivial zero-energy mode is robust against arbitrary fluctuations in hopping strengths, despite the absence of a protecting chiral symmetry. Imagine that in Eq. (54) the hopping strengths $t_{0,j}$ and $t_{1,j}$ are position-dependent. Then, $\widehat{H}$ may be written as

[TABLE]

where $\widehat{G}$ does not contain terms involving $c_{1}$ and $f_{1}$, so that $[\widehat{G},c_{1}]=0=[\widehat{G},f_{1}]$. Then it is easy to verify that the expression for the zero-energy mode is obtained from $\eta_{1}^{\dagger}$ in Eq. (55) after substituting $t_{0}\mapsto t_{0,1}$ and $t_{1}\mapsto t_{1,1}$. We conclude that the zero-energy edge mode is protected by an “emergent symmetry”, that has a non-trivial action only on the sites corresponding to $j=1$. Likewise, assume for concreteness that $t_{0}=\pm t_{1}$, and consider the inter-chain perturbation described by

[TABLE]

In this case, the corresponding single-particle Hamiltonian becomes $H=\mathds{1}_{N}\otimes h_{0}+T\otimes h_{1}+T^{\dagger}\otimes h_{1}^{\dagger}$ with

[TABLE]

Nevertheless, the zero-energy mode corresponding to $|1\rangle|u^{-}\rangle$ is still an emergent solution for $\epsilon=0$, and can be verified to satisfy the boundary equation as well. The topological nature of this zero-energy mode is confirmed by its non-trivial Berry phase bernevig at half-filling. Under periodic BCs, the Hamiltonian in momentum space is

[TABLE]

leading to the following eigenvectors for the two bands:

[TABLE]

Direct calculation shows that the Berry phase has the non-trivial value $\pi$ $(\text{mod}\,2\pi)$, as long as $t_{1}\neq 0$.

V.2.2 Complete closed-form solution

We now obtain a complete closed-form solution of the eigenvalue problem corresponding to Eq. (54) (open BCs). The reduced bulk Hamiltonian is

[TABLE]

with the associated polynomial ($R=1,d=2$)

[TABLE]

The model has two energy bands with a gap containing $\epsilon=0$, and no chiral symmetry. Because $H$ is real, this enforces the symmetry $z\leftrightarrow z^{-1}$ of the non-zero roots of $P(\epsilon,z)$ that satisfy $z_{1}z_{2}=1$. For generic $\epsilon\neq 0$, there are two distinct non-zero roots and, therefore, two extended bulk solutions. The eigenvector of $H(z)$ may be generically expressed as

[TABLE]

Using Eq. (28), the number of emergent bulk solutions is $2Rd-2=2=2s_{0}$, one localized on each edge. As $K^{-}(\epsilon)=h_{1}^{\dagger}$ and $K^{+}(\epsilon)=h_{1}$, such solutions are found from their kernels, spanned by $|u^{-}\rangle$ and $|u^{+}\rangle$, independently of $\epsilon$. The boundary matrix

[TABLE]

whose kernel is nontrivial only if

[TABLE]

In this case, since $z_{1}z_{2}=1$, we may reduce this system to one variable by substituting $z_{2}=z_{1}^{-1}$, which then yields the polynomial equation

[TABLE]

The algebraic system of equations (56) and (57) determine the “dispersing” extended-support bulk modes of the system. When these equations are both satisfied, the kernel of the boundary matrix is spanned by

[TABLE]

and the corresponding eigenvectors of $H$ are given by

[TABLE]

which, upon substituting $z_{1}=e^{ik}$, can be recast as remarkAlternative

[TABLE]

To check whether $|\epsilon\rangle$ in Eq. (58) indeed satisfies the eigenvalue equation, notice that

[TABLE]

Using the expression for $|\epsilon\rangle$, $\langle N|H-\epsilon\mathds{1}|\epsilon\rangle$ vanishes trivially, while, for $j=1$,

[TABLE]

which is seen to vanish from the relation

[TABLE]

The first term on the right hand-side is equal to $P(\epsilon,e^{ik})=0$, whereas the second term vanishes due to Eq. (57). Finally, for $2\leq j\leq N-1$, we get

[TABLE]

which equals zero, completing the argument.

Next, we find the values of $\epsilon$ for which Eq. (56) has a double root. The discriminant of $P(\epsilon,z)$ is $D(P(\epsilon,z))=(\epsilon^{2}-t_{1}^{2})^{2}-4\epsilon^{2}t_{0}^{2}$, and vanishes for $\epsilon=-t_{0}\pm\sqrt{t_{0}^{2}+t_{1}^{2}}$ and $\epsilon=t_{0}\pm\sqrt{t_{0}^{2}+t_{1}^{2}}$, for which the corresponding double roots are $z_{1}=+1$ and $z_{1}=-1$, respectively. In these cases, the bulk equation may have power-law solutions. While one could construct the reduced bulk Hamiltonian $H_{2}(z)$ to identify these solutions, another quick way to proceed is suggested by Eq. (27), as already remarked in Sec. V.1. A power-law solution may now be written as

[TABLE]

where $z_{1}$ is the double root corresponding to $\epsilon$. The first column of the new boundary matrix remains the same as the original one, while its second column is determined from the derivative of the second column of the original boundary matrix with respect to $z_{2}$, computed at $z_{2}=z_{1}$. For $\epsilon=-t_{0}\pm\sqrt{t_{0}^{2}+t_{1}^{2}}$, we have $z_{1}=1$ and

[TABLE]

Some algebra reveals that $\det B(\epsilon)\neq 0$, so that these values of $\epsilon$ do not appear in the spectrum of $H$ for any values of parameters $t_{0},t_{1}$. Similar analysis for $\epsilon=t_{0}\pm\sqrt{t_{0}^{2}+t_{1}^{2}}$ yields the same conclusion. Therefore, there are no power-law solutions compatible with open BCs.

We now derive the perfectly localized zero energy modes described in Sec. V.2.1. Notice that for $\epsilon=0$, the only possible roots of $P(\epsilon,z)$ are $z_{0}=0$, and from its degree it follows that there are $s_{0}=2$ emergent solutions on each edge. In this case,

[TABLE]

with its kernel spanned by

[TABLE]

Similarly, the kernel of $K^{+}(0)$ is spanned by

[TABLE]

Thus, the Ansatz for $\epsilon=0$ consists of all four perfectly localized solutions (see Eqs. (32) and (33)). The boundary matrix in this case is

[TABLE]

which has a two-dimensional kernel, spanned by

[TABLE]

The corresponding two zero-energy edge modes are then

[TABLE]

consistent with the results of Sec. V.2.1. The eigenvector $|\epsilon=0,\bm{\alpha}_{1}\rangle$ has support only on the first site of the two band chain. Since $|N\rangle|u^{+}\rangle=|N\rangle[0\quad 1]^{\rm T}$, the eigenvector $|\epsilon=0,\bm{\alpha}_{2}\rangle$ represents the decoupled degree of freedom at the right end of the chain, as shown in Fig. 4 (a) and (b).

V.3 The Majorana Chain

Kitaev’s Majorana chainkitaev01 is a prototypical model of $p$-wave topological superconductivity alicea12 ; beenakker13 . In terms of spinless fermions, the relevant many-body Hamiltonian in the absence of disorder and under open BCs reads

[TABLE]

where $\mu,t,\Delta\in{\mathbb{R}}$ denote the chemical potential, hopping amplitude, and pairing strengths, respectively. This Hamiltonian, expressed in spin language via a Jordan-Wigner transformation, describes the well-known anisotropic XY spin chain, which has a long history in quantum magnetism, including analysis of boundary effects for both open BCs and periodic liebschultzmattis ; Mikeska ; pfeuty ; XYremark .

Expressed in the form of Eq. (13), the corresponding single-particle Hamiltonian is

[TABLE]

Thus, $R=1$, $d=2d_{\text{int}}=2$, and $h_{R}=h_{1}$ (hence the model) is invertible in the generic parameter regime $|t|\neq|\Delta|$, for arbitrary $\mu$. We have already characterized in detail both the invertible regime abc and the non-invertible regime JPA for generic, regular energy values. While, given the importance of the model, we will summarize some of these results in what follows, our emphasis here will be on (i) addressing singular energy values, in particular, by directly computing compactly-supported eigenstates of flat-band eigenvectors directly in real space; (ii) uncovering the existence of zero-energy Majorana modes with a power-law prefactor, emerging in an invertible but non-generic parameter regime recently discussed in the context of transfer-matrix analysis Hegde16 .

V.3.1 The parameter regime $|t|=|\Delta|$, $\mu\neq 0$

We briefly recall some key steps and results presented in Sec. 5.2 of Ref. [JPA, ]. For concreteness, we assume $t=\Delta$, but a similar analysis may be repeated for the case $t=-\Delta$. The reduced bulk Hamiltonian in this case is

[TABLE]

with associated polynomial

[TABLE]

As in the topological comb example, for generic values of $\epsilon$ the above has two distinct non-zero roots $z_{1}$ and $z_{2}$, which implies a two-dimensional space of extended bulk solutions and one emergent solution on each edge. Let the two extended solutions be labeled by $z_{1}$ and $z_{2}=z_{1}^{-1}$, with $|z_{1}|\leq 1$. Then, we get

[TABLE]

The two emergent solutions are obtained from the one-dimensional kernels of the matrices $K^{-}(\epsilon)=h_{1}^{\dagger}$ and $K^{+}(\epsilon)=h_{1}$, which are spanned by

[TABLE]

respectively. Following Eq. (41), the boundary matrix is

[TABLE]

Our analysis in Ref. [JPA, ] shows that open BCs do not allow any contributions from the emergent solutions in the energy eigenstates, which are linear combinations of the two extended solutions. The condition for $\epsilon$ to be an energy eigenvalue is $\det B(\epsilon)=0$, which simplifies to

[TABLE]

Explicitly, as long as $\epsilon\notin\mathcal{S}\equiv\{\mu\pm 2t,-\mu\pm 2t\}$, the corresponding eigenstate is

[TABLE]

The above equation is particularly interesting for zero energy, since it dictates the necessary and sufficient conditions for the existence of Majorana modes. For $\epsilon=0$, the root $z_{1}$ takes values

[TABLE]

In the large-$N$ limit, the factor $z_{1}^{N+1}$ in the right hand-side of Eq. (61) vanishes thanks to our choice of $|z_{1}|<1$. However, the left hand-side vanishes only in the topologically non-trivial regime characterized by $|\mu|<2|t|$, giving rise to a localized Majorana excitation. The unnormalized Majorana wavefunction in this limit is characterized by an exact exponential decay (see also Fig. 5), namely,

[TABLE]

For the analysis of the non-generic energy values in $\mathcal{S}$, we return to the finite system size $N$. For such $\epsilon$, $P(\epsilon,z)$ has double roots at $z_{1}=1$ and $z_{1}=-1$, so that the bulk equation has one power-law solution in each case JPA . These solutions are compatible with the BCs for certain points in the parameter space, determined by the condition $2tN+\mu(N+1)=0$. Explicitly, the eigenstates corresponding to eigenvalues $\epsilon=\pm(\mu+2t)$ are then

[TABLE]

V.3.2 The parameter regime $|t|=|\Delta|,\ \mu=0$

This regime, sometimes affectionately called the “sweet spot,” is remarkable. Since the analytic continuation of the Bloch Hamiltonian is

[TABLE]

one finds that $\det(H(z)-\epsilon\mathds{1}_{2})=\epsilon^{2}-4t^{2}.$ Thus, the energies $\epsilon=\pm 2t$ realize a flat band and its charge conjugate. From the point of view of the generalized Bloch theorem, these two energies are singular. According to Sec. III.3.3, they necessarily belong to the physical spectrum of the Kitaev chain regardless of BCs, each yielding $\mathcal{O}(N)$ corresponding bulk-localized eigenvectors.

In order to construct such eigenvectors, note that for $\epsilon=\pm 2t$, the adjugate of $H(z)-\epsilon\mathds{1}_{d}$ is the matrix

[TABLE]

which immediately provides two kernel vectors

[TABLE]

In this case, we see that the kernel vectors contain polynomials in $z^{-1}$ of degree $2<\delta_{0}=(d-1)2Rd=4$ (recall Eq. (34)). For a suitable range of lattice coordinates $j$s, the compactly-supported sequences

[TABLE]

yield non-zero solutions $|\Psi_{j\mu,\pm}\rangle=\bm{P}_{1,N}\Psi_{j\mu,\pm},$ $\mu=1,2$, of the bulk equation. However, it is not a priori clear how many of these are linearly independent. For example, it is immediate to check that

[TABLE]

In this case, a basis of compactly-supported solutions can be chosen from the states

[TABLE]

Out of these $N+1$ states, the ones corresponding to $j=1,\dots,N-1$ can be immediately checked to be eigenstates of energy $\epsilon\pm 2t$ fn . In contrast, $|\tilde{\Psi}_{0}\rangle$ and $|\tilde{\Psi}_{N}\rangle$ are not eigenstates: they do not satisfy the boundary equation trivially like other states localized in the bulk. We have thus found $2N-2$ eigenstates of the Hamiltonian, $N-1$ for each band $\epsilon=\pm 2t$.

The two missing eigenstates appear at $\epsilon=0$, which is a regular value of energy and so it is controlled by the generalized Bloch theorem. For $\epsilon=0$, there are four emergent solutions (two on each edge), out of which only

[TABLE]

are compatible with the BCs. Since these solutions are perfectly localized on the two edges, they exist for any $N>2$ (see also Fig. 5). Interestingly, the above states also appeared as solutions of the bulk equation at the singular energies $\epsilon=\pm 2t$, and failed to satisfy the BCs at those values of energy. We do not know whether this fact is just a coincidence or has some deeper significance.

V.3.3 Majorana wavefunction oscillations in the regime $t\neq\Delta$

Recently, it was shown Hegde16 that, inside the so-called “circle of oscillations”, namely, the parameter regime

[TABLE]

the Majorana wavefunction oscillates while decaying in space. Such oscillations in Majorana wavefunction are not observed outside this circle. This observation has consequences on the fermionic parity of the ground state degottardi13 . Because of duality, spin excitations in the XY chain show a similar behavior in the corresponding parameter regime XYremark $B_{z}^{2}=t^{2}-\Delta^{2}=J_{x}J_{y}.$ We now analyze this phenomenon by leveraging the analysis of Sec. III. For simplicity, we address directly the large-$N$ limit.

Clearly, whether a wavefunction oscillates in space depends on the nature of the extended bulk solutions that contribute to the wavefunction. In particular, let $|\psi\rangle=|z,1\rangle|u\rangle$ be one such bulk solution. For a wavefunction to be decaying asymptotically, we must have $|z|<1$. Further, if $z\in{\mathbb{R}}$, then $|\psi_{j}\rangle=z|\psi_{j-1}\rangle$ implies that the part of the wavefunction associated to this bulk solution simply decays exponentially without any oscillations. On the other hand, if $z\equiv|z|e^{i\phi}$ with non-zero phase, then a linear combination of vectors

[TABLE]

can show oscillatory behavior while decaying. This is precisely the phenomenon observed in this case. When $t\neq\Delta$, the reduced bulk Hamiltonian is

[TABLE]

with associated characteristic equation

[TABLE]

For $\epsilon=0$, the above admits four distinct roots in general, out of which two lie inside the unit circle and contribute to the Majorana mode on the left edge. Whether any of these two roots is complex decides if the Majorana wavefunction oscillates for those parameter values. Notice that the characteristic equation is quadratic in the variable $\omega=z+z^{-1}$. We get the two values of $\omega$ to be

[TABLE]

Likewise, notice that for $\mu^{2}<4(t^{2}-\Delta^{2})$, we get both $\omega_{+}$ and $\omega_{-}$ to be complex, which necessarily means that both $z_{1},z_{2}$ inside the unit circle are also necessarily complex. Further, the symmetry of Eq. (63) forces that $z_{2}=z_{1}^{*}$. This leads to the oscillatory behavior of the Majorana wavefunction in the regime $\mu^{2}<4(t^{2}-\Delta^{2})$, that is, inside the circle defined by Eq. (62). Thus, the spatial behavior of Majorana excitations in this regime is formally similar to the solution of an underdamped classical harmonic oscillator (see Fig. 5). Outside the circle, the roots $\omega_{\pm}$ are real. With some algebra, it can be shown that $|\omega_{\pm}|>2$ in this regime, which also means that both $z_{1},z_{2}$ are real roots. This is why oscillations are not observed in this parameter regime, in agreement with the results of Ref. [Hegde16, ]. The Majorana wavefunction in this case resembles qualitatively the solution of a overdamped harmonic oscillator.

The situation when the parameters lie precisely on the circle is particularly interesting. In this case, we find that $\omega_{+}=\omega_{-}\equiv\omega_{0}=-4t/\mu$. Let us assume $t/\Delta>0$ for simplicity. It then follows that $z_{1}=z_{2}=-2(t-\Delta)/\mu$, which rightly indicates appearance of a power-law solution. Let us specifically analyze the case of open BCs on one end (for $N\gg 1$ as stated). One of the two decaying bulk solutions is $|\psi_{1,1}\rangle=|z_{1},1\rangle|u(z_{1})\rangle$, where

[TABLE]

The other bulk solution is obtained from

[TABLE]

The relevant boundary matrix,

[TABLE]

may be computed by relating its second column to the partial derivative of the first column at $z=z_{1}$ as also done previously. Explicitly:

[TABLE]

where we also used Eq. (63) for simplification. Some algebra reveals that $B(0)$ has a one-dimensional kernel, spanned by the vector

[TABLE]

This leads to the power-law Majorana wavefunction

[TABLE]

which decays exponentially with a linear prefactor (see Fig. 5). In principle, the existence of such exotic Majorana modes could be probed in proposed Kitaev-chain realizations based on linear quantum dot arrays Fulga13 , which are expected to afford tunable control on all parameters.

VI An indicator of the bulk-boundary correspondence

As stated in the Introduction, a main motivation behind the development of the generalized Bloch theorem is to elucidate the bulk-boundary correspondence. In this section, we start presenting an indicator of bulk-boundary correspondence based on the results from Sec. III, generalizing the original definition in Ref. [abc, ]. The indicator is built out of the boundary matrix and, therefore, encodes information from the bulk and the BCs. We will then consider an application of the indicator to study the Josephson response of an $s$-wave two-band topological superconductor swavePRL ; swavePRB . Interestingly, and to the best of our knowledge, this system provides the first example of an unconventional (fractional) Josephson effect not accompanied by a fermionic parity switch. We explain the physical reasons behind such a result.

VI.1 Derivation of the indicator

For a system of size $N$, the existence of localized modes at energy $\epsilon$ reflects into a non-trivial kernel of the corresponding boundary matrix, which we now denote by $B_{N}(\epsilon)$ in order to emphasize the dependence on $N$ and $\epsilon$. As we increase $N$ without changing the BCs, the energy $\epsilon$ of the bound modes (that is, modes that remain asymptotically normalizable) attains a limiting value. For instance, in topologically non-trivial, particle-hole or chiral- symmetric systems under hard-wall BCs, the mid-gap bound modes attain zero energy in the large-$N$ limit. This convergence of bound modes and their energies is nicely captured by a modified version of the boundary matrix in the limit $N\gg 1$, which we now construct.

Consider a system of $N$ sites in a ring topology, as shown in Fig. 6(a), so as to allow non-zero contribution from the matrix $w_{bb^{\prime}}$ in the BCs described by $W$ (see Eq. (5)). Let us assume that the system hosts one or more bound modes near the junction formed by the two ends, which converge in the large-$N$ limit to energy $\epsilon$. The resulting modes are the bound modes of a bridge configuration that extends to infinity on both sides, and where the boundary region is shown in Fig. 6(b). For each $N$, we may express the bound eigenstate as in Eq. (36). Such bound states have contributions only from those bulk solutions that are normalizable for $N\gg 1$. The extended-support solutions corresponding to $|z_{\ell}|=1$ are not normalizable, and therefore must drop out from the Ansatz. Further, while the amplitude of those corresponding to $|z_{\ell}|>1$ blows up near $j=N$, they remain normalizable in the limit. This becomes apparent once we rescale such solutions by $z_{\ell}^{-N}$. These rescaled solutions almost vanish at $j=1$ for large $N$. Based on these considerations, we propose a modified Ansatz for finite $N$,

[TABLE]

expressed in terms of up most $2Rd$ amplitudes.

The above Ansatz may be used to compute a corresponding boundary matrix ${B}_{N}(\epsilon)$ in the same way as described in Sec. III.4. Note that ${B}_{N}(\epsilon)$ may not capture the bound modes appearing at finite $N$ since, by construction, it does not incorporate contributions from extended support solutions corresponding to $|z_{\ell}|=1$. However, ${B}_{\infty}(\epsilon)\equiv\lim_{N\rightarrow\infty}{B}_{N}(\epsilon)$ is now well-defined, and describes accurately the presence and exact form of bound modes in the limit. The condition for a non-trivial kernel becomes $\det[{B}_{N}^{\dagger}(\epsilon){B}_{N}(\epsilon)]=0$. Based on this condition, we define the quantity

[TABLE]

as an indicator of bulk-boundary correspondence. This captures precisely the interplay between the bulk properties and the BCs that may lead to the emergence of bound modes, in the sense that, as we parametrically change either or both of the reduced bulk Hamiltonian and the BCs, $\mathcal{D}_{\epsilon}$ shows a singularity at (and only at) the parameter value for which the system hosts bound modes at energy $\epsilon$. Unlike most other topological indicators that are derived from bulk properties (i.e., in a torus topology), our indicator is constructed from a boundary matrix, that incorporates the relevant properties of the bulk. In cases where the bound modes are protected by a symmetry, this allows for the indicator to be computed for arbitrary BCs that respect the symmetry, paving the way to characterizing the robustness of the bound modes against classes of boundary perturbations.

An interesting situation is that of $w_{bb^{\prime}}=0$, in which case the large-$N$ limit consists of two disjoint semi-infinite chains. Then ${B}_{\infty}(\epsilon)$ is block diagonal,

[TABLE]

where $B_{\infty}^{-}$ ($B_{\infty}^{+}$) may be interpreted as the boundary matrix of a semi-infinite chain, describing the edge modes at the left (right) edge, respectively.

While the indicator ${\cal D}_{\epsilon}$ of Eq. (66) signals the presence of bound states, it does not convey information about the degeneracy of that energy level, which is nevertheless contained in the boundary matrix. Therefore, it is often useful to also study the behavior of the degeneracy indicator as a function of $\epsilon$:

[TABLE]

In practice, the dimension of the kernel is obtained by counting the number of zero singular values of ${B}_{\infty}(\epsilon)$.

Remark.— With reference to the discussion in Sec. IV.1, recall that in numerical computations, ${B}_{\infty}(\epsilon)$ signals fictitious roots whenever the bulk equation has a power-law solution. In such cases, we once again remedy the issue by resorting to the Gramian. Then the corrected value of the indicator is given by

[TABLE]

Thus, the correct degeneracy of the energy is obtained by counting zero (within numerical accuracy) singular values of the matrix $\tilde{B}_{\infty}(\epsilon)=B_{\infty}(\epsilon)\mathcal{G}(\epsilon)^{-1/2}$.

VI.2 Application: An $s$-wave topological superconducting wire

The usefulness of the proposed indicator of bulk-boundary correspondence was demonstrated in the context of characterizing the Josephson response of a two-band time-reversal invariant $s$-wave topological superconducting wire in Ref. [abc, ]. While the calculations reported there employed a simplified Ansatz, including only extended-support solutions of the bulk equation, we now validate the analysis by using the complete Ansatz given in Eqs. (36) and (65), and further analyze and interpret our results in terms of fermionic parity switches.

The relevant $s$-wave, spin-singlet, two-band superconductor model swavePRL ; swavePRB derives its topological nature from the interplay between a Dimmock-type intra-band spin-orbit coupling and inter-band hybridization terms. Due to the spin degree of freedom in each of the two relevant orbitals, say, $c$ and $d$, the Nambu basis corresponding to an atom at position $j$ consists of 8 fermionic operators, that we write as the vector

[TABLE]

In this basis, the single-particle Hamiltonian under open BCs is given by

[TABLE]

where the real parameters $\mu,u_{cd},t,\lambda,\Delta$ denote the chemical potential, the interband hybridization, hopping, spin-orbit coupling and pairing potential strengths, respectively, and $\tau_{\alpha},\nu_{\alpha},\sigma_{\alpha}$, $\alpha=\{x,y,z\}$, are Pauli matrices in Nambu, orbital and spin spaces.

The topological properties of the above Hamiltonian were analyzed in Ref. [swavePRB, ]. The BdG Hamiltonian is time-reversal invariant, which places it in the symmetry class DIII. The topological phases may thus be distinguished by a ${\mathbb{Z}}_{2}$-invariant, given by the parity of the sum of the Berry phases for the two occupied negative bands in one of the Kramers’ sectors only swavePRB . For open BCs and for non-vanishing pairing, the system in its trivial phases was found to host zero or two pairs of Majoranas on each edge, in contrast to the topologically non-trivial phase supporting one pair of Majoranas per edge. Similar to the two-dimensional version of the model, one may see that the existence of such Majorana modes is protected by a non-trivial chiral symmetry, of the form $\tau_{y}\sigma_{z}$. The single-particle Hamiltonian $H_{N}$ for open BCs can be exactly diagonalized as described in Sec. IV. In the large $N$-limit, the boundary matrix $B_{\infty}(\epsilon=0)$ calculated by using the Ansatz in Eq. (65) yields degeneracy ${\cal K}_{0}=0,4,8$ in the no-pair, one-pair, and two-pair phases, respectively, verifying the bulk-boundary correspondence previously established through numerical diagonalization.

VI.2.1 Josephson response

In the Josephson ring configuration considered in Ref. [abc, ], the first and last sites of the open chain are coupled by the same hopping and spin-orbit terms as in the rest of the chain, only weaker by a factor of $1/w$. A flux $\phi$ is introduced between the two ends via this weak link. In the large-$N$ limit, this link acts as a junction, with the corresponding tunneling term in the many-body Hamiltonian being given by

[TABLE]

The total Hamiltonian is then $\widehat{H}(\phi)=\widehat{H}_{N}+\widehat{H}_{T}(\phi)$. It was demonstrated abc that the Hamiltonian displays fractional Josephson effect in the topologically non-trivial phase, as inferred from its $4\pi$-periodic many-body ground state energy [Fig. 7(a)], with the phenomenon being observed only if the open-chain Hamiltonian correspondingly hosts an odd number of Majorana pairs per edge. The physics behind the $4\pi$-periodicity was explained in terms of the crossing of a positive and a negative single-particle energy level happening at precisely zero energy as a function of flux $\phi$.

The singular behavior resulting at flux values $\phi=\pi,3\pi$ from the indicator ${\cal D}_{\epsilon=0}(\phi)$ computed using both the simplified Ansatz as in Ref. [abc, ] and the complete Ansatz of Eq. (65) is shown in Fig. 7(c). The qualitative features are clearly unchanged, indicating that in the large-$N$ limit the bound modes formed near the junction are linear combinations only of extended-support solutions, with no contributions from emergent ones. As seen in Fig. 7(d), at both $\phi=\pi$ and $\phi=3\pi$ the junction hosts a total of four Majoranas.

VI.2.2 Parity switch and decoupling transformation

Despite the $4\pi$-periodic Josephson response witnessed in the topologically non-trivial phase, it turns out that the ground state fermionic parity remains unchanged for all flux values. In the non-trivial regime of interest, we may focus on the three low-lying energy levels. Specifically, for values of $\phi<\pi$, let $|\Phi(\phi)\rangle$ denote the many-body ground state, with energy $E_{0}(\phi)$, as in Fig. 7(a). As we will show, there are two degenerate quasi-particle excitations, say, $\eta_{1}(\phi),\eta_{2}(\phi)$, with small positive energy $\epsilon_{0}(\phi)$. This results in a two-fold degenerate first excited many-body state, with energy $E_{1}(\phi)=E_{0}(\phi)+\epsilon_{0}(\phi)$, and a corresponding eigenspace is spanned by $\{\eta_{1}^{\dagger}(\phi)|\Phi(\phi)\rangle,\eta_{2}^{\dagger}(\phi)|\Phi(\phi)\rangle\}$. The second excited state, $\eta_{1}^{\dagger}(\phi)\eta_{2}^{\dagger}(\phi)|\Phi(\phi)\rangle$, is not degenerate and has energy $E_{2}(\phi)=E_{0}(\phi)+2\epsilon_{0}(\phi)$. Note that this state has the same (even) fermionic parity as the ground state. At $\phi=\pi$, the quasi-particle excitation has exactly zero energy, $\epsilon_{0}(\pi)=0$, causing all three energy levels to become degenerate. As $\phi$ crosses $\pi$, $\epsilon_{0}(\phi)$ becomes negative. Therefore, for $\pi<\phi<3\pi$, we find that $E_{2}(\phi)<E_{1}(\phi)<E_{0}(\phi)$. The continuation of the state $\eta_{1}^{\dagger}(\phi)\eta_{2}^{\dagger}(\phi)|\Phi(\phi)\rangle$ with energy $E_{2}(\phi)$ thus becomes the new ground state, whereas the continuation of the original ground state $|\Phi(\phi)\rangle$ now attains the maximum energy among these three levels. Since the new ground state has the same parity as the original one, the system shows no parity switch, with a similar analysis holding for the crossover at $\phi=3\pi$. We conclude that the absence of a fermionic parity switch originates from the twofold degeneracy of the single-particle energy levels.

While the system under open BCs is time-reversal invariant, away from $\phi=0,2\pi$ this symmetry is broken by the tunneling term $\widehat{H}_{T}(\phi)$. Therefore, Kramer’s theorem is not responsible in general for the degeneracy in the single-particle levels. Instead, we now explain the physical origin of this degeneracy in terms of a “decoupling transformation” in real space, thanks to which the system in the Josephson bridge configuration is mapped into two decoupled systems in the same configuration, each with half the number of internal degrees of freedom as the original one. Although each of these smaller systems does undergo a parity switch, the total parity being the sum of individual parities remains unchanged.

Observe that the Hamiltonian $\widehat{H}(\phi)$ is invariant under the unitary symmetries $\hat{S}_{1}$ and $\hat{S}_{2}$, defined by the action

[TABLE]

We can use the eigenbasis of $\hat{S}_{1}$ to decouple $\widehat{H}(\phi)$ into two independent Hamiltonians. Consider, for each site $j=1,\ldots,N$, the canonical transformation

[TABLE]

and let $\hat{U}_{1}$ be the unitary change of basis defined by $\hat{U}_{1}:\hat{\Psi}_{j}^{\dagger}\mapsto[\,\hat{\Psi}_{+,j}^{\dagger}\;\hat{\Psi}_{-,j}^{\dagger}]$, where

[TABLE]

By letting $\hat{\Psi}_{\pm}^{\dagger}\equiv[\hat{\Psi}_{\pm,1}^{\dagger}\;\dots\;\hat{\Psi}_{\pm,N}^{\dagger}]$, the action of $\hat{U}_{1}$ then decouples $\widehat{H}(\phi)$ according to

[TABLE]

where $\widehat{H}_{\pm}(\phi)$ describes two smaller systems, each in a Josephson ring configuration, with hopping and pairing amplitudes given by

[TABLE]

with $\tilde{\sigma}_{\alpha}$ denoting Pauli matrices in the modified spin basis. The decoupling transformation in Eq. (67) is close in spirit to the one already employed under periodic BCs swavePRL ; swavePRB . Indeed, it is worth remarking that $\hat{\Psi}_{+,j}$ and $\hat{\Psi}_{-,j}$ are still time-reversals of each other, in the sense that $\mathcal{T}\hat{\Psi}_{+,j}^{\dagger}\mathcal{T}^{-1}=\hat{\Psi}_{-,j}^{\dagger}$, with ${\cal T}$ being the anti-unitary time-reversal operator for the system. Because of the tunneling term, however, the two decoupled (commuting) Hamiltonians $\widehat{H}_{\pm}(\phi)$ are related by $\mathcal{T}\widehat{H}_{+}(\phi)\mathcal{T}^{-1}=\widehat{H}_{-}(4\pi-\phi)$.

It now remains to show that $\widehat{H}_{\pm}(\phi)$ have identical single-particle energy spectrum, and therefore lead to the desired degeneracy in the energy levels of $\widehat{H}(\phi)$. This follows by examining the symmetries of the single-particle BdG Hamiltonian $H(\phi)$. Corresponding to $\hat{S}_{1}$, $H(\phi)$ has a unitary symmetry $S_{1}=\mathds{1}_{N}\otimes\nu_{x}\sigma_{z}$, and thus gets block-diagonalized into two blocks, $H_{\pm}(\phi)$, upon the action of $U_{1}$. Similarly, corresponding to $\hat{S}_{2}$, $H(\phi)$ has another unitary symmetry $S_{2}=i\mathds{1}_{N}\otimes\tau_{z}\sigma_{x}$. Further, $S_{1}$ and $S_{2}$ satisfy the anti-commutation relation $\{S_{1},S_{2}\}=0$, which is responsible for the doubly degenerate eigenvalue spectrum fn6 . In fact, one can also verify directly that $\widehat{H}_{+}(\phi)$ and $\widehat{H}_{-}(\phi)$ satisfy $\hat{S}_{2}\widehat{H}_{+}(\phi)\hat{S}_{2}^{\dagger}=\widehat{H}_{-}(\phi)$. This explains the origin of the double degeneracy of each single-particle energy level, and hence of the absence of fermionic parity switch.

VII Transfer matrix in the light of the generalized Bloch theorem

Starting with the work in Refs. [HatsugaiPRL, ]-[HatsugaiPRB, ], the transfer matrix has remained the tool of choice for analytical investigations of the bulk-boundary correspondence DelplaceZakPhase ; MongEdge ; MaoAnalytic ; DwivediBB including, as mentioned, recent studies of Majorana wavefunctions in both clean and disordered Kitaev wires Hegde16 . In this section, we revisit the transfer matrix approach to band-structure determination in the light of our generalized Bloch theorem. In particular, we show how, in situations where the transfer matrix fails to be diagonalizable, our analysis makes it possible to give physical meaning to the generalized eigenvectors by relating them to the power-law solutions discussed in Sec. III.3.

VII.1 Basics of the standard transfer matrix method

While our conclusions apply more generally to arbitrary finite-range clean models, for concreteness we refer in our discussion to the simplest setting where both approaches are applicable, namely, a one-dimensional chain with nearest-neighbor hopping. We further focus on open (hard-wall) BCs, as most commonly employed in transfer-matrix studies. The relevant single particle-Hamiltonian $H_{N}$ is then a tridiagonal block-Toeplitz matrix, with entries $h_{1}^{\dagger},h_{0}$ and $h_{1}$ along the three diagonals. Generically, $h_{1}$ is assumed to be invertible. The starting point of the method entails obtaining the recurrence relation between eigenvector components. Specifically, if $|\epsilon\rangle=\sum_{j=1}^{N}|j\rangle|\psi_{j}\rangle$ is an eigenvector of $H$ with energy eigenvalue $\epsilon$ relative to the usual Hilbert-space factorization ${\cal H}={\cal H}_{L}\otimes{\cal H}_{I}$, the components $|\psi_{j}\rangle$ satisfy the recurrence relation

[TABLE]

In terms of the $2d\times 2d$ transfer matrix

[TABLE]

the above recurrence relation may be reformulated as

[TABLE]

where we have written $\bm{P}_{j,j+1}|\epsilon\rangle\equiv\begin{bmatrix}|\psi_{j}\rangle&|\psi_{j+1}\rangle\end{bmatrix}^{\rm T}$. Thus,

[TABLE]

which can be leveraged for obtaining the complete set of eigenvectors of $H_{N}$. We can define $|\psi_{0}\rangle,|\psi_{N+1}\rangle$ by using the relations

[TABLE]

so that $\bm{P}_{N,N+1}|\epsilon\rangle=T(\epsilon)\bm{P}_{0,1}|\epsilon\rangle$ in terms of the matrix $T(\epsilon)\equiv{t(\epsilon)}^{N}$. Hard-wall BCs enforce $|\psi_{0}\rangle=0=|\psi_{N+1}\rangle$. Substituting these boundary values leads to

[TABLE]

which has a non-trivial solution if and only if

[TABLE]

Therefore, all values of $\epsilon$ that obey the above condition are eigenvalues of $H_{N}$. For each eigenvalue, the corresponding $|\psi_{1}\rangle$ is obtained as the kernel of $T_{22}(\epsilon)$. In practice, $T(\epsilon)$ is calculated by first diagonalizing $t(\epsilon)$ by a similarity transformation, and then exponentiating the eigenvalues along its diagonal Lee .

As can be appreciated from this example, the standard version of the transfer matrix method relies on invertibility of certain matrices, although “inversion-free” Boykin ; Biczo or partially inversion-free DwivediBB modifications have also been suggested. In the standard case, the only prerequisite for constructing $t(\epsilon)$ at each step is the banded structure of the single-particle Hamiltonian and, most importantly, the resulting matrix $T(\epsilon)$ is assumed to be diagonalizable.

VII.2 Connections to the generalized Bloch theorem

In order to relate the above analysis to the generalized Bloch formalism, the key observation is to note that the set of equations in Eq. (68) constitute the complete bulk equation, as described in Sec. III.2. Consequently, Eq. (70) is satisfied by any bulk solution $|\psi\rangle\in\mathcal{M}_{1,N}$, where $\mathcal{M}_{1,N}$ denotes the bulk solution space as usual. It is insightful to recast Eq. (71) in the form

[TABLE]

suggesting that the action of the transfer matrix in the bulk solution space is closely related to the one of the left shift $T$. When restricted to $\mathcal{M}_{1,N}$, the above yields the following operator identity:

[TABLE]

with $z\in\mathds{C}$. This relation may be used to establish a direct connection between the basis of the bulk solution space described in the generalized Bloch theorem, and the Jordan structure of the transfer matrix. In the absence of power-law solutions, each bulk solution $|\psi_{\ell s}\rangle$ is annihilated by $\bm{P}_{1,2}(T-z_{\ell}\mathds{1}_{N})=\bm{P}_{1,2}[P_{B}(T-z_{\ell}\mathds{1}_{N})]$. In such cases, Eq. (73) reads

[TABLE]

implying that $\bm{P}_{1,2}|\psi_{\ell s}\rangle$ is an eigenvector of $t(\epsilon)$ with eigenvalue $z_{\ell}$. Naturally, a Bloch wave-like bulk solution corresponds to an eigenvalue on the unit circle, whereas an exponential solution corresponds to one inside or outside the unit circle, in agreement with the literature Lee .

While, as remarked, the transfer matrix is typically assumed to be diagonalizable, we now show that generalized eigenvectors of $t(\epsilon)$ are physically meaningful, and in fact related to the power-law solutions of the bulk equation. Let $\epsilon$ be a value of energy for which power-law solutions are present. We can then generalize our earlier calculation for the eigenvectors of the transfer matrix by noting that each $|\psi_{\ell s}\rangle$ is annihilated by $\bm{P}_{1,2}(T-z_{\ell}\mathds{1}_{N})^{s_{\ell}}$, where $s_{\ell}$ is the multiplicity of the root $z_{\ell}$ as usual. Then, a similar calculation reveals that $\bm{P}_{1,2}|\psi_{\ell s}\rangle$ is a generalized eigenvector of $t(\epsilon)$, satisfying

[TABLE]

Thus, generalized eigenvectors of the transfer matrix are projections of solutions with a power-law prefactor. In some non-generic scenarios, they indeed contribute to the energy eigenstates, as we discussed TMremark .

This analysis is vividly exemplified by the parameter regime corresponding to the circle of oscillations in the Majorana chain, Eq. (62), which we found to be associated to a zero-energy power-law Majorana wavefunction. Accordingly, we expect the corresponding transfer matrix to possess generalized eigenvectors of rank two, failing to be diagonalizable. Let us verify this explicitly. Except for the points $\mu=0,\Delta/t=\pm 1$ in this regime, the matrix $h_{1}$ in Eq. (V.3) is invertible. The transfer matrix is then

[TABLE]

where $\mu,t$ and $\Delta$ satisfy Eq. (62). It can be checked that $t(\epsilon=0)$ has only two eigenvalues, namely, $z_{\ell}={-2(t+(-1)^{\ell}\Delta)}/{\mu},$ $\ell=1,2,$ each of algebraic multiplicity two, and that both of these eigenvalues have only one eigenvector, given by

[TABLE]

hence geometric multiplicity equal to one. Both $z_{1},z_{2}$ are then defective, making $t(\epsilon=0)$ not diagonalizable. In fact, $t(\epsilon=0)$ has one generalized eigenvector of rank two corresponding to each eigenvalue, given by

[TABLE]

Returning to the general case, a number of additional remarks are worth making, in regard to points of contact and differences between the transfer matrix approach and our generalized Bloch theorem. First, the eigenstate Ansatz obtained from the analytic continuation of the Bloch Hamiltonian provides a global characterization of energy eigenvectors (and generalized eigenvectors), as opposed to the local characterization afforded within the transfer-matrix approach, whereby each eigenvector is reconstructed “iteratively” for any given eigenvalue. Further to that, the generalized Bloch theorem unveils the role of non-unitary representations of translational symmetry for finite systems. Perhaps most importantly, the two methods differ in the way BCs are handled. Clearly, in both approaches it is necessary to match BCs in order to obtain the physical energy spectrum. While open BCs are most commonly used in transfer-matrix calculations, the method has also been applied to relaxed surfaces Lee and generalized periodic BCs Molinari , all of which belong to the class of BCs considered in this paper. In this sense, it is tempting to compare Eq. (72) with the condition on the determinant of the boundary matrix, $\det B(\epsilon)=0$. However, the class of BCs to which the transfer matrix approach can be successfully applied is not a priori clear, thus whether such a condition can be established for as general a class of BCs as our theorem covers has not been investigated to the best of our knowledge.

From a numerical standpoint, the computational complexity of the standard transfer matrix method for clean systems (when applicable) is independent of the system size $N$, as is the case of our scan-in-energy algorithm in Sec. IV.1. In those cases where inversion of certain matrices is a difficulty and inversion-free approaches are used Boykin ; Biczo , the latter also have a comparable computational complexity to our method. Interestingly, all approaches so far that are truly inversion-free rely at some point or another on the solution of a non-linear eigenvalue problem TMremark . Thanks to the fact that, as noted, the construction of $t(\epsilon)$ in the generic case relies only on the banded structure of $H_{N}$, bulk disorder can be handled efficiently within transfer-matrix approaches, albeit for a limited class of BCs. For general BCs as we consider, it is thus natural to combine the transfer matrix approach with the bulk-boundary separation we have introduced, in order to still find solutions efficiently: the transfer matrix can be employed to find all possible solutions of the bulk equation in the presence of bulk disorder, and the latter can then be used as input for the boundary matrix, that provides a condition for energy eigenstates.

VIII Discussion and outlook

We have formulated a generalization of Bloch’s theorem applicable to clean systems of independent fermions on a lattice, subject to BCs that are arbitrary – other than respecting the finite-range nature of the overall Hamiltonian. This generalization, which leverages a reformulation of the problem in terms of corner-modified block-Toeplitz matrices, affords exact, analytical expressions for all the energy eigenvalues and eigenstates of the system – which consistently recovers the ones derived from the standard Bloch’s theorem for periodic BCs. As a key component to this theorem, one obtains an exact structural Ansatz, close in spirit to the Bethe Ansatz, for all (regular) energy eigenstates in dispersive bands. This Ansatz is easy to construct since it depends only on the energy eigenvalue and the bulk properties of the Hamiltonian. The individual components of this Ansatz reflect translation invariance in a way we have made precise and are, as such, determined by the analytic continuation of the Bloch Hamiltonian, as shown.

Based on the generalized Bloch theorem, we have provided both a numerical and an algebraic diagonalization algorithm for the class of quadratic Hamiltonians under consideration. For generic energy values, the former is computationally more efficient than existing ones in that its complexity is independent upon the system size; the latter is especially well-suited for symbolic computation or pen-and-paper solutions, as we explicitly demonstrated by solving in closed form a number of tight-binding Hamiltonians of interest, under various BCs. With an eye toward applications in synthetic quantum matter, we have also used the generalized Bloch theorem to engineer a quasi one-dimensional Hamiltonian that support a perfectly localized, robust zero-energy mode, notwithstanding the lack of chiral and charge-conjugation protecting symmetries.

Remarkably, our generalized Bloch theorem predicts the existence, under specific (non-generic) conditions, of edge states that decay exponentially in space with a power-law prefactor. Such exotic states were previously believed to arise only in systems with long-range couplings. In our framework, their origin may be traced back to the description of the system’s eigenstates in terms of non-unitary representations of translation symmetry “outside Hilbert space” – again capturing the fact that such a symmetry is only mildly broken by the BCs, in a precise sense. Notably, we have shown how the emergence of zero-energy Majorana modes with a linear prefactor is possible in the paradigmatic Kitaev chain by proper Hamiltonian tuning on the so-called “circle of oscillations”. Their “critical” spatial behavior separates the theoretically observed Majorana wavefunction oscillations inside such a circle from the simple exponential decay outside.

Our generalized Bloch theorem makes no prediction about the (singular) energy values which correspond to dispersionless, or flat, bands of eigenstates. We have nonetheless provided a prescription for identifying such energy values without diagonalizing the full Hamiltonian, and showed how such energy values necessarily enter the physical energy spectrum irrespective of the BCs. In such singular cases, we have further provided a procedure to effectively obtain a (possibly overcomplete) basis of perfectly localized states using an analytic continuation of the Bloch Hamiltonian, and explicitly illustrated such a procedure in the Kitaev’s Majorana chain Hamiltonian at its sweet spot.

Building on our proposal in Ref. [abc, ], we have rigorously derived and further explored a proposed boundary indicator for the bulk-boundary correspondence. This indicator leverages the other key component to our generalized Bloch theorem, the boundary matrix, and is unique in the sense that, unlike most other indicators in the literature, it combines information from both the bulk and the boundary. The utility of this indicator is seen from our analysis of the $4\pi$-periodic Josephson effect in a model of a $s$-wave topological superconductor. In the process, we show how, remarkably, the $4\pi$-periodicity that distinguishes a topologically nontrivial response is not accompanied by a fermionic parity switch in this system. We have provided a physical explanation of this behavior by exhibiting a decoupling transformation, which maps the relevant Hamiltonian to two uncoupled “virtual” wires – each undergoing a parity switch.

Finally, for systems where no bulk disorder is present, and subject to BCs for which the well-known transfer matrix approach is also applicable, we have shown how the generalized Bloch theorem may be used to obtain a physical interpretation of the transfer matrix’s generalized eigenvectors, in terms of bulk solutions with a power-law prefactor. An explicit example is seen, again, in the semi-infinite Kitaev’s chain with open BCs, precisely in the same circle-of-oscillations parameter regime that hosts power-law zero-energy Majorana modes. While, in this way, our method may be seen to provide yet another inversion-free alternative to the standard transfer-matrix approach, the connections we have identified in this work naturally point to further possibilities for fruitfully combining the two approaches. In particular, since the bulk-boundary separation we proposed remains useful in the presence of bulk disorder, one may envision a hybrid approach for solving disordered systems subject to arbitrary BCs, by employing transfer-matrix techniques to handle the resulting bulk equation.

The tools we have developed here may serve as the starting point for a number of additional studies and applications. As mentioned, in the companion paper PRB2 , we will provide a formulation of the bulk-solution Ansatz and the generalized Bloch theorem further accounting for the role played by the transverse momentum (${\mathbf{k}}_{\perp}$) in higher-dimensional systems with non-trivial boundaries – as opposed to the single ${\mathbf{k}}_{\perp}$-analysis presented here. We will show that topological power-law modes discussed in this paper are not just a feature of one-dimensional systems, and indeed are present in higher dimensions too. Beside exploring the interplay between ${\mathbf{k}}_{\perp}$, the boundary matrix, and the edge states in a number of paradigmatic model Hamiltonians, we will also demonstrate how the treatment of one-dimensional homogeneous systems can be effectively extended to those of interfaces. From a computational standpoint, we expect that the diagonalization algorithms emerging from our approach will be useful for large-scale electronic calculations in both one- and higher- dimensions, possibly in conjunction with perturbative approaches for incorporating interactions.

Towards a deeper understanding of bulk-boundary correspondence in topological insulators and superconductors, our approach can be instrumental in studying robustness against boundary perturbations. It is natural to start by asking how certain symmetries of the system influence the nature of the proposed indicator, or the boundary matrix from which the indicator itself is derived. This can possibly lead to identifying a symmetry principle which dictates the bulk-boundary correspondence, as well as an interpretation at the basic dynamical-system level in terms of stability theory. Likewise, the framework we have developed may also serve as a concrete starting point for rigorously deriving an effective boundary theory for lattice systems.

Lastly, while we have focused on fermions in this paper, the general foundation of our method laid out in Ref. [JPA, ] is equally valid for bosons and immediately applicable to non-Hermitian effective Hamiltonians with non-trivial boundaries, as often arising in semi-classical models of open quantum systems in various contexts nonhermitianrev ; MandalEP ; tayebi16 ; Chong . We plan to explore the corresponding generalized Bloch theorems in forthcoming publications, and to ultimately provide extensions to Markovian open quantum systems described by quadratic Lindblad master equations.

Acknowledgements

We gratefully acknowledge useful discussions with Smitha Vishveshwara. Work at Dartmouth was supported in part by the US NSF through Grant No. PHY-1620541 and the Constance and Walter Burke Special Projects Fund in Quantum Information Science.

Appendix A Further discussion on arbitrary BCs

Section II imposes two restrictions on the allowed form of BCs, described by $\widehat{W}$. The first restricts the non-trivial action of $\widehat{W}$ to the boundary hyperplanes. Since the corresponding single-particle operator $W$ satisfies the relation $P_{B}W=0$, with $P_{B}$ being the bulk projector associated to $H_{N}$, $W$ can be thought of as a corner-modification of the banded block-Toeplitz matrix $H_{N}$. The operators $H_{N}+W$ represent boundary value problems in such a way that a change of BCs is encoded in a change of $W$. The intuition behind these ideas comes from finite-difference methods for solving differential equations. We briefly illuminate this connection here.

Consider for concreteness the Schrödinger boundary value problem

[TABLE]

describing a particle in an infinite one-dimensional potential well. The discretization $x\mapsto x_{j}=j\Delta x$, with $j=0,1,\dots,N+1=L/\Delta x$, reduces this problem to the lattice boundary value problem

[TABLE]

in terms of the centered second difference approximation to the Laplacian. This set of linear equations is equivalent to the eigenvalue equation $(H_{N}-\epsilon\mathds{1}_{N})|\psi\rangle=0$, with

[TABLE]

By comparison, the more general BCs

[TABLE]

lead to the lattice boundary value problem

[TABLE]

together with Eq. (75). The system of linear equations in Eqs. (75)–(77) is equivalent to the eigenvalue problem $(H_{N}+W-\epsilon\mathds{1}_{N})|\psi\rangle=0$, with

[TABLE]

a corner modification of the lattice Laplacian $H_{N}$. For the special case $\alpha_{1}=\alpha_{2},\ \beta_{1}=-\beta_{2}$, we have discussed the exact diagonalization of $H_{N}+W$ in Sec. V.1.

Appendix B Algebras of shift operators

Consider the topologically inequivalent manifolds corresponding to the finite line segment, the circle (of finite or infinite radius), the semi-infinite line, and the infinite line, as illustrated in Fig. 8. Given a physical system whose state space has support on those manifolds, one can define distinct shift (or translation by a distance $a$) operators acting on the physical states. Certainly, those shift operators encode topological information that depending on the circumstances may have physical consequences. In the following we will study the algebra of those shift operators. The subtle difference between the various shift (or translation) operators is reflected in the fundamental discussions that led to the modern theory of macroscopic electric polarization in many-body systems in terms of Berry phases macroP1 ; macroP2 , and the concomitant definition of the position operator in extended systems positionop .

The finite line segment.— This section is based on Ref. [fockpfs, ], where the matrices we are about to consider appeared with a different physical meaning. Consider a line of finite length $L=Na$, written in terms of a characteristic length $a$, typically defined by a periodic potential or lattice. The left shift operator is given by $T=\sum_{j=1}^{N-1}|j\rangle\langle j+1|,$ in terms of the orthonormal lattice states $|j\rangle$. The lattice state $|1\rangle$ is annihilated by $T$, $T|1\rangle=0$, and $|N\rangle$ is annihilated by $T^{\dagger}$, mirroring the fact that the boundary of a line segment consists of two points. For states other than $|1\rangle,|N\rangle$, $T$ and $T^{\dagger}$ act as ordinary translations, to the left or right respectively, i.e., $T|j\rangle=|j-1\rangle$ and $T^{\dagger}|j\rangle=|j+1\rangle$.

While $T$ can be regarded as the generator of bulk translations, it is not a unitary transformation. Instead,

[TABLE]

and notice also that $T^{N}=0$. The commutator $[T,T^{\dagger}]=|1\rangle\langle 1|-|N\rangle\langle N|$ captures the extent of translation-symmetry breaking introduced by the BCs. The lattice-regularized position operator $X=\sum_{j=1}^{N}j\ |j\rangle\langle j|$ satisfies the commutation relation

[TABLE]

While this is formally analogous to $[x,e^{ip/\hbar}]=-e^{ip/\hbar}$, care must be exercised with such analogy, precisely because of issues of definition of the domains of functions where operators act upon.

The circle.— The other compact one-dimensional manifold is the circle. The standard (periodic) left shift operator in this case is given by

[TABLE]

No lattice state $|j\rangle$ is annihilated by either $V$ or $V^{\dagger}$, because the circle is a manifold with no boundary. One can further check that $VV^{\dagger}=\mathds{1}=V^{N}.$ The relation between periodic shifts and the position operator $X$ is better described in terms of $U\equiv e^{i\frac{2\pi}{N}X},$ since then we have the Heisenberg-Weyl relation

[TABLE]

This Heisenberg-Weyl algebra is well-known in statistical mechanics in connection to clock models pclock , but its relevance to tight-binding models appears to have gone unnoticed. The two generators are related by the discrete Fourier transform $F$ as $FUF^{\dagger}=V^{\dagger}$ and $FVF^{\dagger}=U$, see for example Ref. [pclock, ] for more details and references.

By comparing Eq. (78) to Eq. (79), one sees that the $U(1)$ symmetry of the shift algebra associated to the line segment is broken to a $\mathds{Z}_{N}$ symmetry for the circle. In practice, the full $U(1)$ symmetry is recovered by introducing twisted generalizations of the Heisenberg-Weyl algebra, $V_{\phi}U_{\phi}=e^{i\frac{2\pi}{N}}U_{\phi}V_{\phi},$ $V_{\phi}^{N}=e^{i\phi}\mathds{1},$ with $U_{\phi}$ and $V_{\phi}$ unitary. Their meaning is clear in terms of tight-binding models. Twisted Heisenberg-Weyl algebras describe physical problems subject to generalized Born-von-Karman BCs, needed for example for defining topological invariants such as the Chern number. A representation of these algebra is given by

[TABLE]

In statistical mechanics, our twisted Heisenberg-Weyl algebras are connected to chiral Potts models, but this connection seems to be unknown in the literature.

The semi-infinite line.— The left and right unilateral shifts $\bm{T}_{-},\ \bm{T}_{-}^{\star}$ were introduced in Sec. III.3.2. The commutator $[\bm{T}_{-},\bm{T}_{-}^{\dagger}]=|1\rangle\langle 1|$ captures in some sense the extent of translation symmetry breaking. The lattice position operator $X_{-}=\sum_{j=1}^{\infty}j\,|j\rangle\langle j|$ satisfies the commutation relations $[X_{-},\bm{T}_{-}]=-\bm{T}_{-},\quad[X_{-},\bm{T}_{-}^{\star}]=\bm{T}_{-}^{\star}.$ The relation $\bm{T}_{-}^{\star}=\bm{T}_{-}^{\dagger}$ holds if the domain of these linear transformations is restricted to the Hilbert space of square summable half-infinite sequences.

The real line.— The shift operator is $\bm{T}\equiv\sum_{j\in\mathds{Z}}|j\rangle\langle j+1|,$ and it is unitary when restricted to the Hilbert space of square-summable sequences, that is, $\bm{T}^{-1}=\bm{T}^{\dagger}$. We carefully refrained from restricting $\bm{T}$ so in Sec. III.3.1. With $X\equiv\sum_{j\in\mathds{Z}}j\,|j\rangle\langle j|$ (an unbounded Hermitian operator in Hilbert space), one can show that $[X,\bm{T}]=-\bm{T},\quad[X,\bm{T}^{-1}]=\bm{T}^{-1}$ both in and out of Hilbert space.

In summary, the shift operators associated to the finite and the semi-infinite line segment do not commute with their adjoints, reflecting the presence of boundary points for these topologies. In contrast, the shift operators defined on the circle and the line $V$ and $\bm{T}$ do commute with their adjoints (or inverses) and are unitary (or just invertible) – which is why they can represent translation symmetry. As a consequence, $V,V^{\dagger}$ can be diagonalized simultaneously, and the same goes for $\bm{T},\bm{T}^{\dagger}$ footapp2 . Their eigenvalues lay on the unit circle due to unitarity. The key difference between these two types of translation symmetry stems from their interplay with lattice position operators. For all the shift algebras but the one associated to the circle, the position operators generate $U(1)$ rotations of the shift operators. For the Heisenberg-Weyl algebra, this $U(1)$ symmetry appears instead as a family of inequivalent unitary irreducible representations of the defining relation Eq. (79).

Appendix C Emergent solutions at regular energies

This appendix provides further mathematical detail on the procedure for computing emergent bulk solutions outlined in Sec. III.3.2. Specifically, we pick up the discussion where we left it therein, right after the definition of the matrix polynomial $K^{-}(\epsilon,\bm{T}_{-})$ in Eq. (30).

Left-localized emergent bulk solutions.— In analogy to the sequences $\Phi_{z,v}$ associated to ${\bm{T}}$ in Eq. (21), let us define states

[TABLE]

in such a way that $\Upsilon^{-}_{0,v}=|j=v\rangle$ and, also,

[TABLE]

It is then immediate to verify that

[TABLE]

Moreover, using Eq. (80), one also obtains the more general relation

[TABLE]

in terms of the upper-triangular $v\times v$ block matrix

[TABLE]

It will be crucial for later use to notice that $K^{-}_{v}(\epsilon,z)$ is a block-Toeplitz matrix.

Both $K^{-}_{v}(\epsilon,z)$ and $H_{v}(z)$ are defined by the same formula, recall Eq. (24). The key difference between the two is that $K^{-}_{v}(\epsilon,z)$ is well-defined also at $z=0$. So suppose that $z_{0}=0$ is a root of $P(\epsilon,z)$ of multiplicity $s_{0}>0$. Then, one can show using tools from Ref. [JPA, ], that there are precisely $s_{0}$ independent solutions of the equation

[TABLE]

The corresponding emergent bulk solutions are

[TABLE]

They are localized on the left edge over the first $s_{0}$ sites. For Hermitian Hamiltonians, $s_{0}\leq dR$ necessarily.

Right-localized emergent bulk solutions.— Left-localized emergent bulk solutions cannot appear alone; they can only appear in conjunction with a set of right-localized emergent bulk solutions. The reason is as follows. Consider the unitary, Hermitian operator

[TABLE]

which implements a mirror transformation of the lattice, by acting trivially on internal states. The transformed Hamiltonian is the Hermitian block-Toeplitz matrix

[TABLE]

in which the hopping matrices have been exchanged as $h_{r}\leftrightarrow h_{r}^{\dagger}$. Therefore, the left-localized emergent bulk solutions for $\widetilde{H}_{N}$ are dictated by the matrix $\widetilde{K}^{-}(\epsilon)$ with entries $[\widetilde{K}^{-}(\epsilon)]_{ij}=[K^{-}(\epsilon)]_{ij}^{\dagger}$. If $|\widetilde{\psi}^{-}\rangle$ denotes a left-localized emergent solution for $\widetilde{H}_{N}$, then

[TABLE]

implying that the state $U|\widetilde{\psi}_{s}^{-}\rangle=\sum_{j=1}^{s_{0}}|N-j+1\rangle|\tilde{u}^{-}_{sj}\rangle$ is an emergent bulk solution for $H_{N}$, localized on the right edge. Similarly, the left-localized emergent bulk solutions of $H_{N}$ are in one-to-one correspondence with the right-localized emergent solutions of $\widetilde{H}_{N}$. This conclusion relies havily on the commutation relation $P_{B}U=UP_{B}$, which is always necessarily true for closed systems (Hermitian Hamiltonians), as we considered here.

But how can we compute the right-localized emergent bulk solutions directly in terms of $H_{N}$? In Sec. III.3.2, we answered this question with the help of the matrix $K^{+}(\epsilon)\equiv K^{-}(\epsilon)^{\dagger}$. We will justify this answer here. Let $|\widetilde{\psi}_{s}^{-}\rangle=\sum_{j=1}^{s_{0}}|j\rangle|\tilde{u}^{-}_{sj}\rangle$, $s=1,\dots,s_{0}$, denote the left-localized emergent solutions associated to $\widetilde{H}_{N}$, and let

[TABLE]

denote the corresponding right-localized emergent solutions of $H_{N}$, so that $|u^{+}_{sj}\rangle\equiv|\tilde{u}^{-}_{s,s_{0}-j+1}\rangle$. Our goal is to show that the arrays

[TABLE]

are annihilated by $K^{+}(\epsilon)$. Because $|u^{+}_{s}\rangle=\tilde{U}|\tilde{u}^{-}_{s}\rangle$, with $\tilde{U}=\tilde{U}^{\dagger}=\sum_{j=1}^{s_{0}}|j\rangle\langle s_{0}-j+1|$, we conclude that $K^{+}(\epsilon)$ is related to $\widetilde{K}^{-}(\epsilon)$ via $K^{+}(\epsilon)=\tilde{U}\widetilde{K}^{-}(\epsilon)\tilde{U}.$ This leads to the entries

[TABLE]

thanks to the fact that $K^{-}(\epsilon)$ is a block-Toeplitz matrix. Hence, $K^{+}(\epsilon)=[K^{-}(\epsilon)]^{\dagger}$, as desired.