Of Bulk and Boundaries: Generalized Transfer Matrices for TightBinding Models
Abstract
We construct a generalized transfer matrix corresponding to noninteracting tightbinding 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 noninvertible. 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 noncontractible 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.
pacs:
71.15.Ap, 73.20.At, 02.10.UdI 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 Blochlike formalism? As it turns out, the answer is yes, if we allow the quasimomentum to be complex.
Let be a wavefunction for lattice sites indexed by . Then, an imaginary part of the quasimomentum, 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 quasimomentum 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 quasi1D systemMong and Shivamoggi (2011) parametrized by the transverse quasimomentum , 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 timereversal 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 kspace 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 2dimensional 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 bulkboundary correspondence. In particular, using the methods of noncommutative 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 timereversal invariant topological insulatorsProdan (2011). Complementary to this are approaches based on Green’s functionsVolovik (2009); Essin and Gurarie (2011), which have also demonstrated the bulkboundary correspondence from a field theoretic perspective. In addition, aspects of this correspondence have also been discussed from the viewpoint of quantum transport using SmatricesFulga 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 tightbinding 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 quasi1D 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 2torus, 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 2torus, 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 1particle state, as we demonstrate in the following.
ii.1 Outline for tightbinding 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 finitedimensional, momentumdependent 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 quasimomentum (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 tightbinding Hamiltonian of such an arbitrary 1D model parametrized by the transverse quasimomentum , with finite range hopping along the finite direction. In the most general form, we can write such a Hamiltonian as
(1) 
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.
By considering the action of this Hamiltonian on a 1particle state
(2) 
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
(3) 
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
(4) 
Here, is the hopping matrix connecting nearest neighbor supercells and is the onsite matrix, which encodes the hopping between degrees of freedom inside the supercell as well as the onsite energies. The wavefunction for a supercell, , is defined as
(5) 
We shall denote the standard basis of with , where .
For a nonsingular , the transfer matrix construction works by noticing thatAvila et al. (2013)
(6) 
can be rewritten as
(7) 
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 sites
We shall seek to compute the transfer matrix in a basisindependent 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
(8) 
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 degreesoffreedom 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 bands
We perform a reduced singular value decompositionStrang (2009) (SVD) of ,
(9) 
where
(10) 
with
(11) 
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
(12) 
with
(13) 
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
(14) 
with , or, equivalently,
(15) 
with . We have defined analogous to and , so that
(16) 
Also,
(17) 
with and defined in a similar fashion.
We can rewrite the recursion relation in eq. (8), in terms of the Green’s function , as
(18) 
But
(19) 
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 , , etc
(20) 
where the is a matrix. After some matrix gymnastics (see Appendix A.1), the first two equations can be reorganized as
(21) 
with
(22) 
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
(23) 
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
(24) 
Using these and eq. (156), an explicit computation shows that
(25) 
which we can write as
(26) 
However, we can gauge this phase away by the gauge transform
(27) 
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.

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
(28) if
(29) 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 complexsymplectic (). The spectral properties of unitary operators have been studied in great detail in the mathematics literatureSchulzBaldes (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 diagonal
^{4} . 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
(30) 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
(31) 
and the bulk gap as
(32) 
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
(33) 
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
(34) 
which can be solved to get
(35) 
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
(36) 
so that the eigenvalues of the transfer matrix are
(37) 
For instance, for , the explicit expression for the Floquet discriminants is
(38) 
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 bandgap as
(39) 
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
(40) 
Now, span, so that any state can be written as a linear combination of the eigenvectors of the transfer matrix, and
(41) 
We deduce that decays as only if the contribution of the growing eigenvalue is 0, i.e,
(42) 
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)
(43) 
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
(44) 
and similarly, and for and , respectively. Clearly,
(45) 
Then a sufficient condition for to be a left edge state is
(46) 
and similarly for a right edge,
(47) 
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
(48) 
will satisfy this Dirichlet boundary condition on the left edge. Similarly, on the right edge, as , we have
(49) 
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
(50) 
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
(51) 
Then a sufficient condition for to be a left edge state is
(52) 
while for the right edge, we have
(53) 
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 
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
(54) 
or, alternatively,
(55) 
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
(56) 
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
(57) 
where span the cokernel 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
(58) 
Alternatively, in terms of the left eigenvectors, we can also write .
In general, the analogue of this expression is the nonorthogonal representationMeyer (2000) of
(59) 
Hence, the decay condition (eq. (56)) implies , which, using eq. (57), can be written explicitly as
(60) 
which constitutes linear equations for variables . Note that is unique up to a nonzero complex scalar since the righthand 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
(61) 
which is essentially a Cramer’s condition. The analogous right edge conditions reads as
(62) 
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
(63) 
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
(64) 
Substituting the integral representation of from eq. (63), we get
(65) 
where denotes the submatrix 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.
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 relation
(67) 
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
(68) 
Also, the conditions on the Green’s function in eq. 24 reduce to
(69) 
As is real and , . Hence, all transfer matrices for are symplectic, by construction
We can write out the Floquet discriminant, the trace of the transfer matrix, as
(70) 
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 1dimensional, the decay condition requires that be an eigenvector:
(71) 
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
(72) 
The eigenvalue condition () can then be equivalently expressed as
(73) 
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
(74) 
where we have exercised our right to scale by an arbitrary complex number. This then suggests that we choose
(75) 
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.
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
(76) 
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
(77) 
and
(78) 
Going through the machinery above, we obtain the Floquet discriminant as
(79) 
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 flux
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.
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
(80) 
with , where we have defined , then we can write the recursion relation as
(81) 
where and are periodic with period . Following Hatsugai and othersLee and Joannopoulos (1981); Chang and Schulman (1982), we can compute the transfer matrix as
(82) 
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,
(83) 
As , in this basis, , and all other matrix elements of are zero, which is what we demand in eq. (80). Hence, we set and