A Mathematical preliminaries

Of Bulk and Boundaries: Generalized Transfer Matrices for Tight-Binding Models


We construct a generalized transfer matrix corresponding to noninteracting tight-binding lattice models, which can subsequently be used to compute the bulk bands as well as the edge states. Crucially, our formalism works even in cases where the hopping matrix is non-invertible. Following Hatsugai [PRL 71, 3697 (1993)], we explicitly construct the energy Riemann surfaces associated with the band structure for a specific class of systems which includes systems like Chern insulator, Dirac semimetal and graphene. The edge states can then be interpreted as non-contractible loops, with the winding number equal to the bulk Chern number. For these systems, the transfer matrix is symplectic, and hence we also describe the windings associated with the edge states on and interpret the corresponding winding number as a Maslov index.

71.15.Ap, 73.20.At, 02.10.Ud

I Introduction

The topological phases of matter have been a subject of considerable interestBernevig (2013); Shen (2013); Hasan and Kane (2010). In their simplest form, for noninteracting periodic systems, they are characterized by a topological invariant on the Brillouin zone. For instance, the two dimensional topological states may have a nontrivial Chern number of the bundle of the phase of the Bloch eigenstates over the Brillouin zoneThouless et al. (1982).

One of the remarkable features of the topological phases is the existence of boundary states – modes localized on the surface/edge at energies that reside in the bulk band gap, which cannot be gapped out by local perturbations. For topological phases that are insulating in the bulk, the existence of these modes is related to the fact that the topological phase is topologically different from vacuum, a trivial insulator, and hence they cannot be connected without closing the gap.

The topological phases are commonly described by noninteracting lattice Hamiltonians, which can be studied in the framework of Bloch theory (also known as Floquet theory in the differential equation literatureEastham (1973); Brown et al. (2012)). But one of the central tenets of Bloch theory is that we only consider eigenstates to the Hamiltonian that are translation invariant(up to a phase), which is desirable if the system in question is truly periodic, so that the (quasi-)momentum is well defined. However, the presence of an edge naturally breaks that symmetry. Can we still account for the edge modes in a Bloch-like formalism? As it turns out, the answer is yes, if we allow the quasi-momentum to be complex.

Let be a wavefunction for lattice sites indexed by . Then, an imaginary part of the quasi-momentum, Im leads to , a growing or decaying exponential along the direction , which is not normalizable in the case of an infinite system. However, an edge state can be naturally thought of as a decaying exponential (indeed, we can calculate the penetration depth of the edge states), which corresponds to a quasi-momentum with a nonzero imaginary part.

A system with boundary can be naturally described in the language of transfer matricesHatsugai (1993a, b); Lee and Joannopoulos (1981); Chang and Schulman (1982); Tauber and Delplace (2015). By assuming the system to be periodic in directions parallel to the edge/surface, we can reduce it to a family of quasi-1D systemMong and Shivamoggi (2011) parametrized by the transverse quasi-momentum , for which a transfer matrix gives the wavefunction of a block in terms of previous block(s). The eigenvalues of the transfer matrix decide whether the state is periodic or decaying: an eigenvalue with magnitude unity corresponds to Bloch state, while the magnitude being less(greater) than unity corresponds to a decaying(growing) exponential as .

Transfer matrices have been studied in diverse contexts, for instance, electronic band structureHatsugai (1993a); Chang and Schulman (1982); Lee and Joannopoulos (1981); Avila et al. (2013); Agazzi et al. (2014), disorderMacKinnon and Kramer (1983) conductivityAndriotis et al. (2002), Majorana fermionsHegde et al. (2015) and wave motion in electromechanical systemsStephen (2006). They have also been studied by mathematicians as monodromy matrices under the banner of Floquet theoryEastham (1973); Brown et al. (2012); Yakubovich and Starzhinskii (1975); Bronski and Rapti (2011). In condensed matter, they have been used to compute the invariants for time-reversal invariant systemsAvila et al. (2013). However, all of these constructions have been limited to very specific models and/or invertible hopping matrices.

As we will see, the perspective of transfer matrices and a complex k-space is not a mere curiosity or a tool to simplify calculations; instead, it offers interesting insights into the geometry associated with the edge states of the system. For one, as the eigenvalue condition for the transfer matrix, expressed as the characteristic polynomial, is algebraic in energy and , one can associate algebraic curves with the characteristic polynomials. A natural thing to do is to complexify , as the algebraic equation always has roots in the complex plane, which follows from the fundamental theorem of algebra. The characteristic polynomial defines an algebraic variety of codimension 1 on for a -dimensional system, often termed as a Bloch varietyGieseker et al. (1992). We shall not delve much into this picture, however, algebraic curves in complex spaces can also be naturally thought of as Riemann surfaces, a perspective that turns out to be particularly useful, which we shall consider in this article.

One of the first steps in that direction was taken by Y HatsugaiHatsugai (1993b, a), who showed that for the Hofstadter model with flux per plaquette, the edge states correspond to nontrivial windings around the holes of the -Riemann surface, which is a 2-dimensional surface of genus , where is the periodicity of the lattice in presence of the magnetic field, or, equivalently, the number of bulk bands (hence being the number of band gaps). Furthermore, he proves that the winding number so associated with the edge states is equal to the bulk Chern number. The Riemann surface picture was first introduced in condensed matter physics by W KohnKohn (1959), however, it has been investigated in substantial detail in mathematics literatureMcKean and van Moerbeke (1975).

Substantial progress has been made since, with regards to the bulk-boundary correspondence. In particular, using the methods of non-commutative geometry, the equality of the Chern number and the Hall conductivity at finite disorder have been rigorously addressed within the mathematical physics literatureBellissard et al. (1994); Kellendonk et al. (2002), and with generalizations to time-reversal invariant topological insulatorsProdan (2011). Complementary to this are approaches based on Green’s functionsVolovik (2009); Essin and Gurarie (2011), which have also demonstrated the bulk-boundary correspondence from a field theoretic perspective. In addition, aspects of this correspondence have also been discussed from the viewpoint of quantum transport using S-matricesFulga et al. (2012); Bräunlich et al. (2010).

In this article we seek to study the bulk and edge spectra of and the complex geometry associated with generic tight-binding Hamiltonians using transfer matrices, as a continuation of Hatsugai’s analysis beyond the Hofstadter model. We work out a general construction of transfer matrices for quasi-1D systems, including the cases when the hopping matrices are singular. The size of the transfer matrices turns out to be twice the rank of the hopping matrix. We work out analytic computations in some detail for various systems where the hopping matrix is of rank 1, e.g, Hofstadter model, Chern insulator, Dirac semimetal and graphene, as well as for a model with the a rank 2 hopping matrix, viz, the topological crystalline insulator model proposed by FuFu (2011). In all these case, we derive explicit analytic expressions for the bulk band edges as well as the edge state spectra.

For the rank 1 systems, we also describe the topological winding associated with the edge states on the energy Riemann surface, following Hatsugai. The case of Chern insulator is particularly nice, as the energy Riemann surface turns out to be a 2-torus, and we work out Hatsugai’s construction explicitly using elliptic functions. Furthermore, as the transfer matrix for these models is real symplectic matrix, and the associated Lie group is homeomorphic to a solid 2-torus, we plot the transfer matrix as a function of the transverse momentum for a given edge spectrum using an explicit parametrization of . Using the fact that , we identify the corresponding winding number as a Maslov index, which, in all the cases discussed, turns out to be equal to the bulk Chern number.

Finally, we use our transfer matrix formalism to study Chern insulator on a rectangular geometry in presence of diagonal disorder, where we marry our general transfer matrix formalism with the conventional numerical methodsKramer and MacKinnon (1993) to compute localization lengths and the their scaling. We then demonstrate the existence of edge states for a disordered Chern insulator.

The rest of this article is organized as follows: in §II, we describe our general construction of the transfer matrix and discuss their properties and applications. In §III, we describe the computations for rank 1 systems, taking Chern insulator as the prototypical model. In §IV, we construct the -Riemann sheet and describe the windings associated with edge states. In §V, we analytically compute the transfer matrix for a rank 2 system, the TCI. In §VI, we apply our formalism to Chern insulator with diagonal disorder. Finally, we conclude and make general comments in §VII. Details of some of the lengthy calculations and mathematical facts relevant to the work are relegated to the appendices.

Ii Transfer matrices

Transfer matrices arise naturally in discrete calculus as a representation of finite order linear difference equations, where it is an operator that implements a first order shift on a block. In general, a transfer matrix will depend on the independent variable; however, if the system is periodic, the transfer matrix that translates by a single period acquires a special significance, as its repeated application can propagate a solution as far as one wishes. This matrix is often known as the monodromy matrix in dynamical systems literatureEastham (1973); Brown et al. (2012).

As noninteracting lattice models are essentially composed of hopping (i.e, shift) operators, which act on the wavefunctions, the Schrödinger equation can alternatively be written as difference equations (or recursion relations) by action on a 1-particle state, as we demonstrate in the following.

ii.1 Outline for tight-binding lattice Hamiltonians

The tight binding lattice Hamiltonians, in presence of (discrete) translation symmetry, are diagonal in the momentum basis, thereby reducing the computation of the bulk spectrum to the diagonalization of a finite-dimensional, momentum-dependent Bloch Hamiltonian. However, as we are interested in the edge states, the translation symmetry is naturally broken in the direction normal to the edge, as the system is finite in that direction. We shall assume that for a system in space dimensions, we have periodic boundary conditions(PBC) along directions parallel to the edge, so that the corresponding quasi-momentum (the dimensional torus) is a good quantum numberMong and Shivamoggi (2011); Ramamurthy and Hughes (2014), which reduces a given -dimensional system into a family of 1D system parametrized by .

Consider, then, the tight-binding Hamiltonian of such an arbitrary 1D model parametrized by the transverse quasi-momentum , with finite range hopping along the finite direction. In the most general form, we can write such a Hamiltonian as


where is the range of the hopping and we have internal degrees of freedom (spin/ orbital/ sublattice) per site of the lattice. In the second line, we have bundled up the creation/annihilation matrices corresponding to orbitals in the -vectors , while is the corresponding hopping matrix with nearest neighbors. We shall suppress the explicit dependence on in the following equations to avoid notational clutter; however, all parameters should be assumed to depend on , unless stated otherwise.

Figure 1: (a) A schematic depiction of the recursion relation, with internal degrees of freedom, range of interaction and Dirichlet boundary conditions at the left edge. We can form blocks (supercells) of such sites with 2 sites each, so the there is only nearest neighbor hopping between them. (b) A simplified depiction of the reduced recursion relation, with , , corresponding to the coefficients of , and subspaces(introduced in §II.2), respectively. (c) We club together with to obtain , which is translated by one step using the transfer matrix.

By considering the action of this Hamiltonian on a 1-particle state


where is the fermionic vacuum state and is the wavefunction for each physical site, the eigenvalue problem can be written as a recursion relation in , as


We now construct blocks consisting of these sites, so that the system is periodic in these blocks and the hopping between such blocks is restricted to nearest neighborLee and Joannopoulos (1981), as shown in Fig. 1). These blocks form the sites of a superlattice. We shall hereafter refer to those blocks as supercells. Note that this is always possible as the hopping has a finite range, hence we can always choose a supercell consisting of sites. In terms of these supercells, each containing degrees of freedom, the recursion relation becomes


Here, is the hopping matrix connecting nearest neighbor supercells and is the on-site matrix, which encodes the hopping between degrees of freedom inside the supercell as well as the on-site energies. The wavefunction for a supercell, , is defined as


We shall denote the standard basis of with , where .

For a nonsingular , the transfer matrix construction works by noticing thatAvila et al. (2013)


can be rewritten as


Hence, we have a transfer matrix . However, this does not work for a singular , which is often the case.

What exactly does mean? Think of the degrees of freedom inside each supercell as sites1. Then, is the number of linearly independent rows in , and hence the number of degrees of freedom that enter in the recursion relation, when expressed in a suitable basis. In more physical terms, denotes the number of bonds between adjacent supercells. Hence, physically, the singularity of means that there are sites in the supercell from which one cannot hop directly to a site in another supercell. In other words, if the degrees of freedom in a supercell are thought of as nodes of a graph, where and encode the connectivity of the graph, then there are nodes in a supercell that are not connected to any nodes in other supercells, and is the number of links between the supercells. We seek to compute a transfer matrix for singular , where we, in some sense, mod out the redundant degrees of freedom, thereby inverting on a reduced subspace to get a reduced transfer matrix.

We shall seek to compute the transfer matrix in a basis-independent fashion, i.e., without referring to the explicit forms of the and matrices. However, we present an explicit calculation for the case of Chern insulator in Appendix D which motivated us for the following general construction.

ii.2 Constructing the transfer matrix

We begin with the recursion relation


Let . We will see that the corresponding transfer matrix will be . Indeed, if had full rank (), we could have inverted it to get a transfer matrix, as computed in §II.1. In the following, we shall also assume a big enough supercell that is nilpotent of degree 2, i.e, , so that . Physically, for , this simply means that in a given supercell, the nodes in a supercell that are connected to the right neighboring supercell and the left neighboring supercell are not directly connected to each other.

Now consider . We note that is the resolvent (or the Green’s function) of a single supercell. Clearly, is singular when is an eigenvalue of . What does that mean? Consider a system with the uncoupled -degrees-of-freedom supercells, obtaining by tuning to in the recursion relation of eq. (8). The corresponding spectrum consists of degenerate levels. As we turn on the hopping , these degenerate levels broaden into bands. Hence, the eigenvalues of can be interpreted as the centers of the bands. Since we are primarily concerned with the band gaps and the edge states therein, we can take to be nonsingular as far as we do not venture deep inside the bands2.

We perform a reduced singular value decompositionStrang (2009) (SVD) of ,






The first two expressions follow from the definition of SVD, while implies the third. Hence, is required to ensure that the and subspaces are orthogonal and the corresponding coefficients can be extracted by taking inner products.

The SVD can equivalently be written as




We shall hereby refer to these vector pairs as channels. As we can still change the phases of and without violating the orthonormality, we choose their phases such that all the singular values are positive, i.e, . Clearly, .

Now, morally speaking, we claim that the only directions in relevant for the problem are ’s and ’s, i.e, span and span. Take a basis of as , where , and expand as


with , or, equivalently,


with . We have defined analogous to and , so that




with and defined in a similar fashion.

We can rewrite the recursion relation in eq. (8), in terms of the Green’s function , as




which follows from the SVD, eq. (11) and eq. (16). We can now premultiply eq. (18) by , and to extract the coefficients , and , respectively. In order to simplify notation, we denote the restriction of to and subspaces by , , etc3. Then,


where the is a matrix. After some matrix gymnastics (see Appendix A.1), the first two equations can be reorganized as




Hence, we have managed to construct a closed form expression for a transfer matrix explicitly for the given recursion relation. This is one of our central results.

Defining , we can also express this result as


This expression is somewhat cleaner, but it obscures the different physical significance associated with and , as well as properties of , which we now state. As the Green’s function is Hermitian, i.e, , we have


Using these and eq. (156), an explicit computation shows that


which we can write as


However, we can gauge this phase away by the gauge transform


In the following, whenever we refer to the transfer matrix, we shall assume that we have gauged away the phase of the determinant of so that .

ii.3 Properties

Before we go on to compute physically relevant quantities from the transfer matrix, we discuss a few features of our construction:

  • The transfer matrix propagates and degrees of freedom. Given one of the ’s, say, , we can compute for all , and, using the expression for in eq. (20), compute the wavefunction . Furthermore, as is nonsingular by construction, we can also use to propagate backwards.

  • The transfer matrix is basis independent, as we have never referred to the explicit form of the and matrices. It reduces the computation of transfer matrix for a system to the identification of the and matrices, as everything else can be mechanized. We shall illustrate that with examples later on (see Table 2).

  • The size and spectral properties of the transfer matrix are independent of the size of the supercell chosen, once it is above a certain size. Hence, we can define a minimal supercell, which is a block consisting of the minimum number of sites so that the hopping between the supercells is nearest neighbor and the corresponding hopping matrix is nilpotent. In Appendix B, we show that if we take a supercell that is times the minimal supercell, the transfer matrix is simply exponentiated by , i.e, , but its size, which is twice the rank of the hopping matrix, stays invariant under this operation. But as in computing the band structures, we are concerned only with the behavior of for large (See §II.4), the band structure, as expected, stays invariant under such a transformation. Hence, we can always make the supercell bigger than the minimal supercell, while leaving the bands and edge states invariant. We shall use this property in certain proofs.

  • As are simply restrictions of the Green’s functions, they are propagators connecting the and degrees of freedom for each supercell, while encodes the tunneling probabilities, or the relative strength of each channel. In fact, the recursion equation in terms of and (eq. (20)) has a simple diagrammatic interpretation as superpositions of possible nearest neighbor hopping processes. This is illustrated in Fig. 2, where the Green’s functions express the propagation within a block and the tunneling between blocks. The transfer matrix equation (21) can then be seen as an equation of constraint that respects these hopping processes.

    Figure 2: Diagrammatic representation of the recursion relations (eq. (20)) for (a) and (b) .
  • In Floquet theory for a continuous independent variable, the monodromy matrix is symplectic if the system is HamiltonianBronski and Rapti (2005). Is that also true for the discrete case? An explicit computation using eq. (24) shows that




    Physically, this condition implies that the various channels that connect the nearest neighbor supercells need to be independent, so that the order of tunneling () and propagation () is irrelevant. We term this condition as being -unitary or complex-symplectic (). The spectral properties of -unitary operators have been studied in great detail in the mathematics literatureSchulz-Baldes (2014). Furthermore, if is real, we say that is symplectic (. In the discussion on bulk bands, we show that if the transfer matrix is symplectic, it can effectively be decomposed into a set of chains, one corresponding to each channel. The conditions on obtained above are physical manifestations of that fact. As is, by definition, a diagonal matrix, in order for it to commute with another matrix , , in general, must also be diagonal4. Hence, for to be symplectic, and must also be diagonal.

  • Recall that a complex square matrix A is termed normal if it commutes with its adjoint, i.e, if


    The matrix is diagonalizable by a unitary matrix if and only if it is normal. In other words, the eigenvectors of a matrix form an orthonormal basis if and only if it is normal, a condition often ignored in physics literature. However, as it turns out, the transfer matrices are almost never normal. Hence, in the subsequent arguments, we shall not assume that is normal in general, which, naturally, makes them somewhat more involved.

ii.4 Using the transfer matrix

We now discuss the computation and interpretation of spectra from the transfer matrices.

Bulk bands

The transfer matrix can be used to propagate a state spatially into the system. The eigenstates with an eigenvalue propagates indefinitely if , i.e, the lies on the unit circle in the complex plane, while it grows/decays as is lies outside/inside the unit circle. Hence, a given lies in the bulk band if all eigenvalues of lie on the unit circle, while it lies in the gap if all of them lie off the unit circle.

Formally, let be the spectrum of for a point . We define the bulk band, , as


and the bulk gap as


For , the possibility exists that there can be points for which some eigenvalues are on and some off the unit circle. We shall term such points partial gaps, , defined as


By construction, each falls in one of these sets.

To compute the bulk bands, one needs to compute the eigenvalues of the transfer matrix. This can always be done numerically in a given case, however, if the transfer matrix is symplectic, its characteristic polynomial has further structure(see Appendix C) which allows us to compute the eigenvalues analytically.

Let us start off with , where implies that is symplectic. It also implies that the product of eigenvalues is unity, so that the eigenvalues are reciprocals of each other, so that


which can be solved to get


Hence, either both the eigenvalues are on the unit circle or both on the real line. In turn, a given either belongs to or , so that .

For , if the transfer matrix is symplectic, the eigenvalues always come in reciprocal pairs, i.e, if is an eigenvalue, so is . Hence, given a transfer matrix, we construct the Floquet discriminants using the traces of powers of (see Appendix C for details), where


so that the eigenvalues of the transfer matrix are


For instance, for , the explicit expression for the Floquet discriminants is


Hence, if the transfer matrix is symplectic, we can essentially decompose it into a set of systems, which are independent of each other!

From the expression for the eigenvalues, we deduce that we have an oscillating state for and a growing/decaying state for . Hence, we can alternatively define the bulk band and the band-gap as


Furthermore, the band edges are simply given by the conditions . Hence, we can simply solve this conditions for , numerically if needed, to compute the band edges, without having to diagonalize for all possible , which one would need to do in general.

Decay conditions

The edge states are typically the states that reside outside the bulk bands, , which implies that they are growing/decaying as . In order to be normalizable, they are taken to be decaying into the bulk away from the edges. Typically, one is interested in the existence of these states, and, should they exist, in the edge spectrum, i.e, the energy of the edge state, , as a function of the transverse momentum .

Given for an arbitrary site , we can use the transfer matrix to compute . Hence, given a , we are concerned with the asymptotics of for , where is the vector norm over . A can be a legitimate left edge state if as . Similarly, can be a legitimate right edge state if as . In this subsection, we seek the conditions imposed on by the transfer matrix (i.e, the bulk) for it to be a legitimate decaying edge state, while we defer the implications of the boundary condition to the next subsection. In the following, we consider the left edge, the arguments for the right edge being their exact analogues.

Consider first the case when is normal and satisfies the conventional eigenvalue equation


Now, span, so that any state can be written as a linear combination of the eigenvectors of the transfer matrix, and


We deduce that decays as only if the contribution of the growing eigenvalue is 0, i.e,


This restricts to a subspace of corresponding to .

However in general, is not a normal matrix and hence it cannot have an eigenvalue equation in the usual sense and cannot be diagonalized by a unitary matrix, a fact that is often overlooked in the physics literature. Nevertheless, it can always be brought to a Jordan canonical form Kato (1995). The behavior of is still dictated by the generalized eigenvalues of , so that to ensure exponential spatial decay of as , we now require that contain no generalized eigenvectors with eigenvalues . Similarly, for the right edge, we want to decay as . Hence, the corresponding condition demands that contain no generalized eigenvectors with eigenvalues .

In order to obtain a precise mathematical condition for the subspaces, we express in its Jordan canonical form Kato (1995)


where the sum extends over all generalized eigenvalues , project in the eigenspace of and are nilpotent matrices. However, it remains true that the determinant of is equal to the product of its generalized eigenvalues , as can be seen by applying a similarity transform to (43). We define the projector to the subspace as


and similarly, and for and , respectively. Clearly,


Then a sufficient condition for to be a left edge state is


and similarly for a right edge,


A rigorous proof of this statement is provided in Appendix A.3. We shall term these decay conditions.

Boundary conditions

We now discuss the boundary conditions required to compute the physical edge states of the system, as observed in an exact diagonalization of the lattice models on finite size lattices. Most of the following is a restatement of the results by Lee and JoannopoulosLee and Joannopoulos (1981) in our formalism, which, we believe, is more general. In the following, we shall only consider the system that is terminated abruptly at layer 0 and , hereafter termed a hard boundary condition.

We mention in passing that as the edge state spectrum is strongly dependent on the boundary conditions, it can get modified quite drastically by local terms at the boundary, an effect commonly known as edge reconstruction. In order to consider the most general case, we should take a Hamiltonian , where is an operator localized at the edge, which can account for the edge reconstruction, for instance, due to an impuritySong et al. (2010) or lattice deformationPark et al. (2013). Such a boundary condition imposes additional conditionsLee and Joannopoulos (1981); Dang et al. (2014); Mao et al. (2010) on the eigenvectors of the transfer matrix. As our purpose in this article is to expound the geometry and topology associated with the band structure which is independent of such local deformations, we shall not discuss such cases in detail.

For concreteness, we consider just a left edge since the right edge is analogous. The hard boundary conditionHatsugai (1993b) is the simplest Dirichlet boundary condition at an edge, whereby we simply demand that , which leads to . Hence, any initial state vector


will satisfy this Dirichlet boundary condition on the left edge. Similarly, on the right edge, as , we have


Note that are still undetermined for on their respective edges. We shall use the decay conditions to fix these in the next section.

Formally, we can also define projectors to write the boundary condition in a way similar to the decay conditions. We begin by defining the matrices


as the injectors into the and subspaces, respectively. In terms of these operators, the Dirichlet condition on the left edge is equivalent to the statement that , while the right edge is equivalent to . Finally, define the projectors


Then a sufficient condition for to be a left edge state is


while for the right edge, we have


These are our boundary conditions.

Physical edge states

We have obtained two sets of conditions, viz, the decay conditions and the boundary conditions, that we need to solve simultaneously in order to obtain the physical edge states. However, before we attempt to do so, we can ask a somewhat perverse question, which turns out to have important consequences: What if we chose the wrong decay condition for a given boundary? We tabulate the situation as follows:

Left edge Unphysical
Unphysical Right edge
Table 1: Boundary(rows) vs decay(column) conditions.

The wrong choice of decay condition implies that the corresponding state grows (instead of decaying) exponentially in the bulk, and is hence not normalizable and unphysical. However, we shall see that in order to account for all the windings corresponding to the edge state, we shall need to take the unphysical states into account. Furthermore, these should not be thought of as a complete fantasy, as they can be revealed by changing the boundary condition, as we shall demonstrate explicitly in §III.4

At this point, we can compute the physical edge states by solving the decay conditions and the boundary conditions simultaneously. For instance, for the left edge state, we seek to simultaneously solve


or, alternatively,


Note that as , this is a homogeneous linear system of up to equations for the variables, viz, the coefficients of . But for a nontrivial state, we demand that , from which we can obtain a Cramer’s condition, which can be numerically solved to obtain the physical edge spectrum.

In the following, we also analytically construct a closed form expression combining the decay and the boundary conditions for the case when there are an equal number of eigenvalues are inside and outside the unit circle in the complex plane, which corresponds to an , i.e, in the bulk gap. This implies that , so that and


We seek to represent in terms of the (generalized) eigenvectors of . Let be the generalized eigenvalues of with corresponding left and right generalized eigenvectors being ’s and ’s. Furthermore, let us assume that that lies outside the unit circle for while it lies inside the the unit circle for . Then, we define the left and right subspaces corresponding to as


where span the co-kernel and range of , respectively.

If were normal, i.e, diagonalizable by a unitary transform, then the right eigenvectors ’s form an orthonormal basis of . As projects along a subset of these eigenvectors, it is an orthogonal projection, which can be written as


Alternatively, in terms of the left eigenvectors, we can also write .

In general, the analogue of this expression is the non-orthogonal representationMeyer (2000) of


Hence, the decay condition (eq. (56)) implies , which, using eq. (57), can be written explicitly as


which constitutes linear equations for variables . Note that is unique up to a non-zero complex scalar since the right-hand sides are all zero. Thus the space of unique solutions really is the complex projective valued. The equations (60) have a nontrivial solution if and only if


which is essentially a Cramer’s condition. The analogous right edge conditions reads as


These conditions incorporate both the boundary and decay conditions and can be solved numerically to obtain as a function of to obtain the edge spectrum, .

Equation (61) is very convenient for numerical computations, but we also present an alternative characterization which is more explicit in terms of ’s projection. The general spectral decomposition of the resolvent of Kato (1995) yields


Essentially, we use the fact that the integrand has poles whenever equals an eigenvalue of so that .

Now, in the simpler case of a normal , we have , so that


Substituting the integral representation of from eq. (63), we get


where denotes the sub-matrix of the argument and we have expressed the dependence of explicitly. Such an equation, though impractical for numerical computations, make explicit the analytic properties of an edge dispersion in open neighborhoods where it exists as a solution.

In the most general case where is oblique(non-orthogonal), the analogue of eq. (64) is


where is still given by the integral equation (63).

Iii The case of

Now that we have a hammer, we look for a nail. The simplest nontrivial case for our formalism corresponds to . In the following, we shall see that this case offers further simplifications, as well as additional structure that is not present, or at least not immediately obvious, in the higher rank cases.

iii.1 Transfer matrix

Let us start with an explicit calculation of the transfer matrix using eq. (22). For , is a matrix, i.e, a number, which we can set to by a suitable normalization of the recursion relation5. We write the transfer matrix as


where we have defined the restricted determinant as

The prefactor becomes after we gauge away the phase of by the gauge transform from eq. (27), as


Also, the conditions on the Green’s function in eq. 24 reduce to


As is real and , . Hence, all transfer matrices for are symplectic, by construction6.

We can write out the Floquet discriminant, the trace of the transfer matrix, as


The band edges are given by , which can be used to solve for , at least locally. Note that enters the calculation only as , which is a rational function of , so that solving for the band edges is equivalent to finding the zeros of a polynomial in .

The case of the edge states is also particularly simple for . As the subspaces corresponding to are either or 1-dimensional, the decay condition requires that be an eigenvector:


where is the eigenvector of corresponding to the eigenvalue , i.e, in the bulk gap. We get the analogous condition for the right edge.

We start off with a somewhat less restrictive condition which can be represented in a neat geometric way. We demand simply that , be eigenvectors of (hereby referred to as the eigenvalue condition). We note that for any antisymmetric and , we have


The eigenvalue condition () can then be equivalently expressed as


In dynamical systems literature, the function is often referred to as the Evans functionBronski and Rapti (2011). This is equivalent to the statement that satisfies either the left or the right decay condition, i.e, it lies entirely in the growing or the decaying subspace. We can later check whether these states are physical by computing the corresponding eigenvalues.

The hard boundary condition (eq. (48) and (49)) implies that


where we have exercised our right to scale by an arbitrary complex number. This then suggests that we choose


such that is the symplectic form with respect to the and subspaces. Such a choice of automatically incorporates the hard boundary data into the Evans function.

We term a solution of the equations (73),(74),(75) for either as the edge spectrum. Note that this includes the physical as well as the unphysical edge states, as defined in §II.4. These conditions describe a curve in the space, which has an associated winding number.

iii.2 Hofstadter model

We start off by repeating Hatsugai’sHatsugai (1993a) calculation in our formalism. The Hofstadter Hamiltonian, after a partial Fourier transform along , the direction with PBC, is given by


where . The system is periodic with period , and we get a gapped system with edge states for odd . We club together physical sites to make a supercell, so that has all entries equal to zero except , while has as its diagonal entries while it has 1’s on the first diagonal. For instance, for the simplest nontrivial case of , these matrices are




Going through the machinery above, we obtain the Floquet discriminant as


In general, for arbitrary , is a polynomial in of order . The edge state calculation is identical to Hatsugai’s, so we shall not discuss it in any detail. However, we emphasize his remark that if the total number of sites is commensurate with the flux , i.e, a multiple of , then for a given , we either get edge states on both left and right edges, or no edge states at all. In order to have an edge state for all , which will have an associated winding, we need to consider a system with the number of sites incommensurate with the flux7 (see fig. 3).

In our picture, for the latter case, the number of supercells is not an integer. This makes physical sense for a Hofstadter model as the degrees of freedoms per supercell are physical sites for the Hofstadter model, so that we can remove those sites. In general, the degrees of freedom inside a supercell are not physical sites. However, we shall see that the number of supercells being fractional still formally makes sense, and hence we can contrive a (potentially unphysical) boundary conditions for those cases which will exhibit the winding of the edge states. We shall hereafter use the word incommensurate (with the superlattice) to refer to the cases where the number of supercells is not an integer.


Figure 3: (color online) The spectrum of Hofstadter model, with the band edges (dark blue) computed using the transfer matrix formalism and the left and right edge state dispersion (dashed and dashed-dot) from the Evans equation (73), overlaid on the spectrum computed using exact diagonalization for a (top) commensurate and (bottom) incommensurate system. Note that in the latter case, the edge states seen in exact diagonalization exactly follow the winding right edge state obtained from the transfer matrix.

iii.3 Natural basis and “unfolding”

Hatsugai’s calculation of the transfer matrix worked because of the fact that the system had nearest neighbor hopping. Before we proceed to further examples, we stop to consider the implications of nearest neighbor hopping inside a single supercell, which implies, in our notation, that is tridiagonal, and , as in the Hofstadter model. Explicitly, if


with , where we have defined , then we can write the recursion relation as


where and are periodic with period . Following Hatsugai and othersLee and Joannopoulos (1981); Chang and Schulman (1982), we can compute the transfer matrix as


where is the transfer matrix from site to site with periodicity . This construction always results in the transfer matrix being polynomial in , as it simply involves a product of matrices linear in . Subsequently, the Floquet discriminant, is a polynomial in , a fact which we shall use when discussing the winding in §IV.

Now, an interesting aspect of our computation of the transfer matrix is that it is basis independent, as we have not used the explicit form of or anywhere in this calculation. Can we choose a basis for our system where and are of the form in eq. (80)? Let us assume so, and let such a basis of be .

We start off by noting that a natural orthonormal basis for , viz,


As , in this basis, , and all other matrix elements of are zero, which is what we demand in eq. (80). Hence, we set and