A unified perspective of complex band structure: interpretations, formulations, and applications

At its core, conventional band structure is rooted in the translational symmetry of a crystalline material, where periodicity of the potential allows invocation of Bloch's theorem. We thus consider wavefunctions of the form

$$\begin{eqnarray}&&{{\psi}_{n,\vec{k}}}\left(\vec{r}\right)={{\text{e}}^{\text{i}\vec{k}\cdot \vec{r}}}{{u}_{n,\vec{k}}}\left(\vec{r}\right),\end{eqnarray} \tag{ 1 }$$

where $\vec{k}$ is a wavevector (usually in the first Brillouin zone), n is the band index, and u is a function with periodicity of the lattice. When used in the Schrödinger equation, $\vec{k}$ can be treated as a parameter (i.e. each $\vec{k}$ produces an independent system), and the eigenstates for a given $\vec{k}$ are calculated using standard techniques. The ensuing dispersion relation for the nth band, ${{E}_{n}}(\vec{k})$, details the material's band structure.

Despite the numerous applications of this conventional band structure, many systems of interest lack perfect translational symmetry. Surfaces and interfaces, for instance, break translational symmetry in (at least) one direction. These systems still possess a 'repeat unit' similar to the periodic unit cell, and equation (1) may remain applicable if solutions with complex $\vec{k}$ are also admitted [3–6].

Imaginary components of $\vec{k}$, when present, indicate evanescence; that is, the state grows or decays in magnitude from one repeat unit to the next. These states have sometimes been called 'generalized Bloch functions'. Conversely, bulk states that propagate throughout the material have real $\vec{k}$. For physicality, we restrict our attention to complex wavevectors that describe states with real energies [4–6], and the set of such $\vec{k}$ constitutes the material's complex band structure (CBS). In other words, CBS extends the dispersion relation to complex $\vec{k}$, such that it describes both the propagating bulk states from conventional band structure and also any evanescent states that are forbidden by translational symmetry [5].

Given CBS's broad applicability, it has found use in myriad contexts: Surfaces and interfaces [2, 5, 8–18, 19–29, 30–40, 41–51], the construction of Wannier functions [3, 6, 20, 52–56], impurities [54, 55, 57–67], high-energy electron diffraction [68], superlattices [40, 41, 66, 69–77], heterostructures [35, 46, 78–88], quantum wells [40, 41, 85, 89–91], magnetic systems [92–100], electron transport [48, 64, 66, 85, 92, 101–130], nanomaterials [64, 66, 77, 104, 117, 131–137], topological materials [67, 77, 138–144], solar cells [145], quantum size effects [15, 27, 39, 133–136, 146–149], and many others. However, these numerous applications often independently redevelop the main concepts and results of CBS, frequently with different notations or nomenclature. As a result, there are several formulations of CBS that appear unrelated at a quick inspection. Concrete examples of these formulations include the transfer or companion matrices [28–30, 33, 35, 42, 57, 67, 73, 74, 77, 83, 87, 90, 102, 105, 122, 141, 150–156], the propagation matrix [11, 16, 23], bivariational methods [157], wavefunction matching techniques [2, 7, 9, 38, 39, 44, 82, 89], and Green function-based approaches [24, 26, 27, 34, 41, 43, 135, 158, 159]. Note that these terms are very broad; e.g. several different 'transfer matrices' have been reported. To exacerbate these inconsistencies, many of the CBS derivations are unnecessarily tied to the choice of basis set and/or are laden with application-specific details. Consequently, the handful of studies that compare different CBS formulations [33, 38, 40, 41, 44, 105, 115, 119, 153, 155, 157, 160, 161] laid crucial groundwork but have yet to illuminate the full generality of CBS.

In this work we present (section 3) a unified perspective of CBS that delineates the key underpinnings of CBS that are often unstated or unclear. From these we develop an interpretation for CBS that relates to the amount of information in the material's Hamiltonian. We then use this realization to review previous studies (sections 4–9) and demonstrate the commonalities among the various CBS formulations.

As seen in our examples from section 2, the first step is to identify each 'layer' within the system. The layer doesn't necessarily have to be an atomic layer, but is simply a way to partition the system and define the repeat units. Mathematically, this is accomplished using a complete set of orthogonal (Hermitian) projection operators [167]. For example, let ${{\mathbf{P}}_{j}}$ be the projector for layer j such that $\mathbf{I}=\sum\nolimits_{j}{{\mathbf{P}}_{j}}$. Then, ${{\mathbf{H}}_{j,j}}\equiv {{\mathbf{P}}_{j}}\mathbf{H}{{\mathbf{P}}_{j}}$ denotes the 'block' of the Hamiltonian corresponding to layer j; likewise,

$$\begin{eqnarray}&&{{\mathbf{H}}_{j,k}}\equiv {{\mathbf{P}}_{j}}\mathbf{H}{{\mathbf{P}}_{k}}\end{eqnarray} \tag{ 3 }$$

is the coupling between layers j and k. Because it is built with operators, this projector-based notion of a 'matrix block' sidesteps any concerns about basis set locality or nonorthogonality and permits a discussion of CBS that is unencumbered by the specific choice of basis set [152]. Briefly returning to our examples in section 2, figure 1 specified the projectors

$$\begin{eqnarray}&&{{\mathbf{P}}_{n}}=\mathbf{a}_{n,+}^{\dagger}{{\mathbf{a}}_{n,+}}+\mathbf{a}_{n,-}^{\dagger}{{\mathbf{a}}_{n,-}},\end{eqnarray}$$

whereas the larger layers in figure 2 used

$$\begin{eqnarray}&&\begin{array}{*{20}{l}}{{\mathbf{P}}_{n}}=\mathbf{a}_{2n+1,+}^{\dagger}{{\mathbf{a}}_{2n+1,+}}+\mathbf{a}_{2n,-}^{\dagger}{{\mathbf{a}}_{2n,-}}\\ \,\,\,\,\,\,\,\,\,\,+\mathbf{a}_{2n,+}^{\dagger}{{\mathbf{a}}_{2n,+}}+\mathbf{a}_{2n-1,-}^{\dagger}{{\mathbf{a}}_{2n-1,-}}.\end{array}\end{eqnarray}$$

We also require (for now) the projectors to be chosen such that each layer only couples to its nearest neighbors. The size of each layer can always be made larger to satisfy this condition, and we will discuss relaxing it in section 9. Mathematically, the Hamiltonian then assumes a block tridiagonal form,

$$\begin{eqnarray}\mathbf{H}=\left[\begin{array}{c} {{\mathbf{H}}_{1,1}} & {{\mathbf{H}}_{1,2}} & \mathbf{0} & \mathbf{0} & \cdots \\ {{\mathbf{H}}_{2,1}} & {{\mathbf{H}}_{2,2}} & {{\mathbf{H}}_{2,3}} & \mathbf{0} & \cdots \\ \mathbf{0} & {{\mathbf{H}}_{3,2}} & {{\mathbf{H}}_{3,3}} & {{\mathbf{H}}_{3,4}} & \cdots \\ \mathbf{0} & \mathbf{0} & {{\mathbf{H}}_{4,3}} & {{\mathbf{H}}_{4,4}} & \cdots \\ \vdots & \vdots & \vdots & \vdots & \ddots \end{array}\right],\end{eqnarray} \tag{ 4 }$$

with ${{\mathbf{H}}_{j,k}}=\mathbf{H}_{k,j}^{\dagger}$ (by Hermiticity). As enumerated in equation (4), the system has a 'surface' to the left of layer 1 as there is nothing to the left of layer 1. Accordingly, layer 1 is the surface layer, layer 2 the first subsurface layer, and so forth to the 'bulk limit' at $\infty$ [135, 136, 149]. In the absence of a surface, the indices would extend to $-\infty$.

If we further neglect other forms of disorder (such as surface reconstructions or impurities), each layer becomes identical to the others. This is the repeat unit. The Hamiltonian is thus also block Toeplitz; that is, blocks are the same along a particular diagonal. We denote ${{\mathbf{H}}_{j,j}}={{\mathbf{H}}_{\text{D}}}$ and ${{\mathbf{H}}_{j,j-1}}={{\mathbf{H}}_{\text{S}}}$ (${{\mathbf{H}}_{j-1,j}}=\mathbf{H}_{\text{S}}^{\dagger}$) such that

$$\begin{eqnarray}\mathbf{H}=\left[\begin{array}{c} {{\mathbf{H}}_{\text{D}}} & \mathbf{H}_{\text{S}}^{\dagger} & \mathbf{0} & \mathbf{0} & \cdots \\ {{\mathbf{H}}_{\text{S}}} & {{\mathbf{H}}_{\text{D}}} & \mathbf{H}_{\text{S}}^{\dagger} & \mathbf{0} & \cdots \\ \mathbf{0} & {{\mathbf{H}}_{\text{S}}} & {{\mathbf{H}}_{\text{D}}} & \mathbf{H}_{\text{S}}^{\dagger} & \cdots \\ \mathbf{0} & \mathbf{0} & {{\mathbf{H}}_{\text{S}}} & {{\mathbf{H}}_{\text{D}}} & \cdots \\ \vdots & \vdots & \vdots & \vdots & \ddots \end{array}\right].\end{eqnarray} \tag{ 5 }$$

(The '$\text{D}$' and '$\text{S}$' subscripts stand for 'diagonal' and 'sub-/super-diagonal', respectively.) For the time being, we assume that ${{\mathbf{H}}_{\text{S}}}$ is invertible. Section 8 will discuss relaxing this condition.

For computational purposes (e.g. the examples in section 2), we assume that each layer can be sufficiently described by M basis functions, such that the operators ${{\mathbf{H}}_{\text{D}}}$ and ${{\mathbf{H}}_{\text{S}}}$ can be represented by $M\times M$ matrices. Matrix representations for the small layers in figure 1 are $2\times 2$,

$$\begin{eqnarray}{{\mathbf{H}}_{\text{D}}}=\left[\begin{array}{c} \alpha & \beta _{1}^{\ast} \\ {{\beta}_{1}} & -\alpha \end{array}\right]~\text{and}~{{\mathbf{H}}_{\text{S}}}=\left[\begin{array}{c} 0 & {{\beta}_{2}} \\ {{\beta}_{2}} & 0 \end{array}\right],\end{eqnarray}$$

whereas the larger layers in figure 2 require $4\times 4$ matrices,

$$\begin{eqnarray}{{\mathbf{H}}_{\text{D}}}=\left[\begin{array}{c} \alpha & 0 & \beta _{2}^{\ast} & 0 \\ 0 & \alpha & \beta _{1}^{\ast} & \beta _{2}^{\ast} \\ {{\beta}_{2}} & {{\beta}_{1}} & -\alpha & 0 \\ 0 & {{\beta}_{2}} & 0 & -\alpha \end{array}\right]~\text{and}~{{\mathbf{H}}_{\text{S}}}=\left[\begin{array}{c} 0 & 0 & {{\beta}_{2}} & {{\beta}_{1}} \\ 0 & 0 & 0 & {{\beta}_{2}} \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{array}\right].\end{eqnarray}$$

Note that ${{\mathbf{H}}_{\text{S}}}$ for the larger layers is not invertible, and provides an example for the discussion in section 8. These basis set and matrix details are only necessary for computation and (unless otherwise noted) will not be considered when formulating or discussing the theory of CBS.

This block tridiagonal and block Toeplitz structure in the Hamiltonian is foundational to CBS, as it mathematically leads to CBS. We now detail this idea.

Most applications aim to either diagonalize or (essentially) invert the Hamiltonian. The former appears in the context of the time-independent Schrödinger equation,

$$\begin{eqnarray}&&\mathbf{H}|\psi \rangle =E|\psi \rangle,\end{eqnarray} \tag{ 6 }$$

and the latter in finding the (retarded/advanced) Green function (GF) [168],

$$\begin{eqnarray}&&\mathbf{G}(E)=\underset{\eta \to {{0}^{\pm}}}{{\lim}}{{\left[\left(E+\text{i}\eta \right)\mathbf{I}-\mathbf{H}\right]}^{-1}}.\end{eqnarray} \tag{ 7 }$$

As foreshadowed, the block tridiagonal and block Toeplitz structure of $\mathbf{H}$ greatly aids in these processes.

To begin, consider the amount of information in $\mathbf{H}$. All that is needed to completely describe our system is ${{\mathbf{H}}_{\text{D}}}$, ${{\mathbf{H}}_{\text{S}}}$, and the number of layers (i.e. the number of block rows/ columns). With respect to the size of our system (the number of layers) this is $\mathcal{O}(1)$ information; that is, the amount of information does not change with system size. Because $\mathbf{H}$ uniquely determines the eigenvalues (eigenvectors) and the GF, there must only be $\mathcal{O}(1)$ information in these quantities. CBS is this $\mathcal{O}(1)$ information. All reported formulations of CBS access this information to solve specific problems, e.g. diagonalizing $\mathbf{H}$ for the time-independent Schrödinger equation. In this manner, CBS describes, and can be used to understand, all of a material's static and dynamic electronic properties.

Such an analysis of a matrix's information content is not new to the mathematics community, and is one facet of studying (block) quasi-separable (sometimes called rank-structured) matrices [169, 170]. As it suffices for our purposes (a more rigorous definition is presented in the footnote¹), a block quasi-separable matrix has certain sub-parts with specific rank. For the present discussion, the sub-parts that contain a sub-diagonal (or super-diagonal) block and the bottom-left (top-right) corner, but do not contain a diagonal block, will have rank equal to that of ${{\mathbf{H}}_{\text{S}}}$. Graphical examples of these sub-parts are

Clearly, block tridiagonal matrices (with or without the block Toeplitz structure) belong to the class of block quasi-separable matrices.

Ramifications of this block quasi-separable structure have been discussed in numerous contexts [169, 170]. In many cases, they include analytical characterizations and/or fast numerical algorithms for inverting or diagonalizing a matrix. References [112, 158, 171–180] overview specific results for block tridiagonal matrices. In the end, these results stem from finding and exploiting the 'generators' of the matrix [181]; that is, the minimal information needed to generate the entire matrix, its spectral decomposition, or its inverse. The generators are not unique, and some generators may be preferential to others for specific operations (such as inverting or diagonalizing the matrix).

Returning to our block tridiagonal and block Toeplitz Hamiltonian, we see that the set $\left\{{{\mathbf{H}}_{\text{D}}},{{\mathbf{H}}_{\text{S}}},N\right\}$ is a generator for $\mathbf{H}$ (where N is the number of layers). We will show in the following sections that the various formulations of CBS produce different generators, as needed for specific applications, but that they all represent the same information (i.e. CBS). For example, the transfer matrix [30] and the companion matrix [33] help in diagonalizing $\mathbf{H}$, whereas generators for (in effect) inverting $\mathbf{H}$ to obtain the GF have also been developed [26, 159]. In the following sections we derive some of these generators (that is, formulate CBS) and show how they aid in diagonalizing or inverting $\mathbf{H}$. Along the way, we prove that these generators embody the same information (CBS) and then discuss some uses and applications of CBS.

The MMT [159, 185] is a generalization of the Möbius (bilinear) transformation from complex variables [162] to matrices. Let ${{\mathbf{m}}_{11}}$, ${{\mathbf{m}}_{12}}$, ${{\mathbf{m}}_{21}}$, ${{\mathbf{m}}_{22}}$, and $\mathbf{z}$ be complex $n\times n$ matrices; a MMT is a mapping with the form

$$\begin{eqnarray}&&\mathbf{M}\bullet \mathbf{z}\equiv \left({{\mathbf{m}}_{11}}\mathbf{z}+{{\mathbf{m}}_{12}}\right){{\left({{\mathbf{m}}_{21}}\mathbf{z}+{{\mathbf{m}}_{22}}\right)}^{-1}}.\end{eqnarray} \tag{ 14 }$$

Similar to the complex Möbius transformation [162], we can associate with $\mathbf{M}$ the $2n\times 2n$ matrix [159]

$$\begin{eqnarray}\mathbf{M}=\left[\begin{array}{c} {{\mathbf{m}}_{11}} & {{\mathbf{m}}_{12}} \\ {{\mathbf{m}}_{21}} & {{\mathbf{m}}_{22}} \end{array}\right].\end{eqnarray}$$

A MMT and its matrix representation will be used interchangeably throughout this work. Note that the '•' in equation (14) emphasizes that $\mathbf{M}\bullet \mathbf{z}$ is not a canonical matrix-matrix product, but the action of a MMT. Lastly, it is easily verified that the application of MMTs is associative; that is, for any two MMTs ${{\mathbf{M}}_{1}}$ and ${{\mathbf{M}}_{2}}$,

$$\begin{eqnarray}&&{{\mathbf{M}}_{1}}\bullet \left({{\mathbf{M}}_{2}}\bullet \mathbf{z}\right)=\left({{\mathbf{M}}_{1}}{{\mathbf{M}}_{2}}\right)\bullet \mathbf{z},\end{eqnarray} \tag{ 15 }$$

where ${{\mathbf{M}}_{1}}{{\mathbf{M}}_{2}}$ is the usual matrix-matrix product.

As we will see in section 5.3, the GF for a block tridiagonal Hamiltonian can be formulated in terms of surface GFs. This section is, therefore, devoted to surface GFs, which have been thoroughly discussed elsewhere [123, 133, 135, 158, 159, 186, 187]. Herein we develop a layer-by-layer approach for surface GFs that parallels the transfer matrix for wavefunctions.

Before proceeding, we need to define the term 'surface GF' and specify our notation. Using the layer enumeration in equation (4), the surface GF is the (1, 1) block of the total system GF; that is,

$$\begin{eqnarray}&&\widetilde{\mathbf{G}}_{\infty}^{\text{R}}(E)={{\mathbf{P}}_{1}}\mathbf{G}(E){{\mathbf{P}}_{1}}\equiv {{\mathbf{G}}_{1,1}}(E).\end{eqnarray} \tag{ 16 }$$

Three comments are needed to clarify this notation for surface GFs. First, the tilde signifies a surface GF, as opposed to the full system GF in equation (7). Second, the surface GF depends on the orientation of the material relative to the surface. For example, the system described in equation (4) places the material to the right of the surface, where the surface is effectively between layers 0 and 1. We use the superscript '$\text{R}$' to label this case; '$\text{L}$' is used when the material is to the left of the surface. For reference, a left surface GF is the bottom-right block of the total system GF. Third, surface GFs for systems containing a finite number of layers will also be important; $\widetilde{\mathbf{G}}_{N}^{\text{L}/\text{R}}(E)$ denotes the left/right surface GF for a system with N layers. As demonstrated in equation (16), $\widetilde{\mathbf{G}}_{\infty}^{\text{L}/\text{R}}(E)$ is the left/right surface GF for a semi-infinite system.

We now show that MMTs can build surface GFs in a layer-by-layer fashion. To overview the procedure, we use Löwdin partitioning [188] (which is essentially an embedding technique [51, 189]) to construct an effective Hamiltonian for the surface layer that (i) only exists in the surface layer and (ii) is rigorously equivalent to the total system Hamiltonian within the surface layer. This happens at the expense of having an energy-dependent effective Hamiltonian. Because the surface GF is the block of $\mathbf{G}$ for the surface layer, the surface GF is identical to the GF of this effective Hamiltonian. Equation (21) formally states this result, and also relates the surface GF for a system with N layers to the surface GF for a system with N  −  1 layers. MMTs then exploit equation (21) to obtain the surface GF for an arbitrary number of layers; this realization leads to CBS in the next section. In what follows we detail this idea for right surface GFs, noting that a similar procedure is used for left surface GFs.

For concreteness, figure 3 displays various surface DOSs,

$$\begin{eqnarray}&&\widetilde{\rho}_{n}^{\text{L}/\text{R}}(E)=\frac{-1}{\pi}\text{Im}(\text{Tr}[\widetilde{\mathbf{G}}_{n}^{\text{L}/\text{R}}(E)]),\end{eqnarray} \tag{ 17 }$$

for our ongoing example (using the layers from both figures 1 and 2). When there are only a few layers in the system, the surface DOS is essentially a collection of δ functions (one for each eigenvalue of $\mathbf{H}$). As the system gets larger, the δ functions coalesce into bands, ultimately giving rise to the continuous surface DOSs seen for the semi-infinite systems. When there is a mirror plane between the layers (as in figure 1), the valence and conduction bands are anti-symmetrically shaped. There is no such similarity between the bands when the mirror plane is absent (e.g. the layers of figure 2).

Zoom In Zoom Out Reset image size

As in the derivation of the transfer matrix, we begin with the time-independent Schrödinger equation. To set up Löwdin partitioning, suppose that ${{\mathbf{P}}_{1}}$ is the projector for our surface layer and ${{\mathbf{Q}}_{1}}\equiv \mathbf{I}-{{\mathbf{P}}_{1}}={{\mathbf{P}}_{2}}+{{\mathbf{P}}_{3}}+\ldots$ is the projector for the rest of the system. If we substitute $\mathbf{I}={{\mathbf{P}}_{1}}+{{\mathbf{Q}}_{1}}$ before and after $\mathbf{H}$ in equation (6), distribute, and rearrange terms, we get the matrix equation

$$\begin{eqnarray}\left[\begin{array}{c} {{\mathbf{P}}_{1}}\mathbf{H}{{\mathbf{P}}_{1}} & {{\mathbf{P}}_{1}}\mathbf{H}{{\mathbf{Q}}_{1}} \\ {{\mathbf{Q}}_{1}}\mathbf{H}{{\mathbf{P}}_{1}} & {{\mathbf{Q}}_{1}}\mathbf{H}{{\mathbf{Q}}_{1}} \end{array}\right]\left[\begin{array}{c} {{\mathbf{P}}_{1}}|\psi \rangle \\ {{\mathbf{Q}}_{1}}|\psi \rangle \end{array}\right]=E\left[\begin{array}{c} {{\mathbf{P}}_{1}}|\psi \rangle \\ {{\mathbf{Q}}_{1}}|\psi \rangle \end{array}\right].\end{eqnarray} \tag{ 18 }$$

The bottom row can be rearranged to produce

$$\begin{eqnarray}&&\left(E{{\mathbf{Q}}_{1}}-{{\mathbf{Q}}_{1}}\mathbf{H}{{\mathbf{Q}}_{1}}\right){{\mathbf{Q}}_{1}}|\psi \rangle ={{\mathbf{Q}}_{1}}\mathbf{H}{{\mathbf{P}}_{1}}|\psi \rangle.\end{eqnarray}$$

We observe that the term in parentheses on the left-hand side is the inverse of the GF for the isolated system of all non-surface layers. If we invert this expression² to obtain the GF for that subsystem (denoted by ${{\mathbf{G}}_{{{\mathbf{Q}}_{1}}}}$), we get

$$\begin{eqnarray}&&{{\mathbf{Q}}_{1}}|\psi \rangle ={{\mathbf{G}}_{{{\mathbf{Q}}_{1}}}}(E){{\mathbf{Q}}_{1}}\mathbf{H}{{\mathbf{P}}_{1}}|\psi \rangle.\end{eqnarray} \tag{ 19 }$$

This gives us the wavefunction outside the surface layer (${{\mathbf{Q}}_{1}}|\psi \rangle$) in terms of the wavefunction in the surface layer (${{\mathbf{P}}_{1}}|\psi \rangle \equiv |{{\psi}_{1}}\rangle$). Substituting equation (19) into the top row of equation (18) and rearranging terms yields

$$\begin{eqnarray}&&{{\mathbf{P}}_{1}}\mathbf{H}{{\mathbf{P}}_{1}}|{{\psi}_{1}}\rangle +{{\mathbf{P}}_{1}}\mathbf{H}{{\mathbf{Q}}_{1}}{{\mathbf{G}}_{{{\mathbf{Q}}_{1}}}}(E){{\mathbf{Q}}_{1}}\mathbf{H}{{\mathbf{P}}_{1}}|{{\psi}_{1}}\rangle =E|{{\psi}_{1}}\rangle.\end{eqnarray}$$

We now have a time-independent Schrödinger equation exclusively in the surface layer, where the effective Hamiltonian is

$$\begin{eqnarray}&&{{\mathbf{P}}_{1}}\mathbf{H}{{\mathbf{P}}_{1}}+{{\mathbf{P}}_{1}}\mathbf{H}{{\mathbf{Q}}_{1}}{{\mathbf{G}}_{{{\mathbf{Q}}_{1}}}}(E){{\mathbf{Q}}_{1}}\mathbf{H}{{\mathbf{P}}_{1}}.\end{eqnarray} \tag{ 20 }$$

Mathematically, this effective Hamiltonian is related to Schur complements [190]. The second, energy-dependent term is sometimes called a self-energy, and can be interpreted as an embedding potential or as an open-system boundary condition.

The surface GF for our total system is then obtained by inverting this effective Hamiltonian within the surface layer's subspace (where ${{\mathbf{P}}_{1}}$ is the identity). If our system has N layers, we would obtain the surface GF for a system with N layers,

$$\begin{eqnarray}&&\widetilde{\mathbf{G}}_{N}^{\text{R}}(E)={{\left[E{{\mathbf{P}}_{1}}-{{\mathbf{P}}_{1}}\mathbf{H}{{\mathbf{P}}_{1}}-{{\mathbf{P}}_{1}}\mathbf{H}{{\mathbf{Q}}_{1}}{{\mathbf{G}}_{{{\mathbf{Q}}_{1}}}}(E){{\mathbf{Q}}_{1}}\mathbf{H}{{\mathbf{P}}_{1}}\right]}^{-1}}.\end{eqnarray}$$

The block tridiagonal structure of $\mathbf{H}$ now allows significant simplifications to this expression. First, ${{\mathbf{P}}_{1}}\mathbf{H}{{\mathbf{Q}}_{1}}={{\mathbf{H}}_{1,2}}{{\mathbf{P}}_{2}}$ because the surface layer (e.g. layer 1) only couples to layer 2. A similar statement holds for ${{\mathbf{Q}}_{1}}\mathbf{H}{{\mathbf{P}}_{1}}$. Then,

$$\begin{eqnarray}&&\widetilde{\mathbf{G}}_{N}^{\text{R}}(E)={{\left[E{{\mathbf{P}}_{1}}-{{\mathbf{H}}_{1,1}}-{{\mathbf{H}}_{1,2}}{{\mathbf{P}}_{2}}{{\mathbf{G}}_{{{\mathbf{Q}}_{1}}}}(E){{\mathbf{P}}_{2}}{{\mathbf{H}}_{2,1}}\right]}^{-1}}.\end{eqnarray}$$

Second, layer 2 of the total system is a surface layer of the isolated ${{\mathbf{Q}}_{1}}$ subsystem, meaning that ${{\mathbf{P}}_{2}}{{\mathbf{G}}_{{{\mathbf{Q}}_{1}}}}(E){{\mathbf{P}}_{2}}$ is a surface GF for that subsystem, which has one fewer layer. Thus,

$$\begin{eqnarray}&&\widetilde{{\mathbf{G}}}_{N}^{\text{R}}(E)={{[E{{\mathbf{P}}_{1}}-{{\mathbf{H}}_{1,1}}-{{\mathbf{H}}_{1,2}}\widetilde{{\mathbf{G}}}_{N-1}^{\text{R}}(E){{\mathbf{H}}_{2,1}}]}^{-1}}^{{}}.\end{eqnarray}$$

Finally, exploiting the block Toeplitz structure and, for simplicity, replacing ${{\mathbf{P}}_{1}}$ by the identity, we get

$$\begin{eqnarray}&&\widetilde{{\mathbf{G}}}_{N}^{\text{R}}(E)={{[E\mathbf{I}-{{\mathbf{H}}_{\text{D}}}-\mathbf{H}_{\text{S}}^{\dagger}\widetilde{{\mathbf{G}}}_{N-1}^{\text{R}}(E){{\mathbf{H}}_{\text{S}}}]}^{-1}}^{{}}.\end{eqnarray} \tag{ 21 }$$

Equation (21) relates the surface GF for a system with N layers to the surface GF for a system with N  −  1 layers. After some minor rearrangement, we can also show that a MMT describes this relationship between $\widetilde{\mathbf{G}}_{N}^{\text{R}}(E)$ and $\widetilde{\mathbf{G}}_{N-1}^{\text{R}}(E)$. Starting from equation (21),

$$\begin{eqnarray}\begin{array}{*{20}{l}}\widetilde{{\mathbf{G}}}_{N}^{\text{R}}(E)&=\mathbf{H}_{\text{S}}^{-1}{{[\left(E\mathbf{I}-{{\mathbf{H}}_{\text{D}}}\right)\mathbf{H}_{\text{S}}^{-1}-\mathbf{H}_{\text{S}}^{\dagger}\widetilde{{\mathbf{G}}}_{N-1}^{\text{R}}(E)]}^{-1}}^{{}}\\ &=\left[\begin{array}{c} \mathbf{0} & \mathbf{H}_{\text{S}}^{-1} \\ -\mathbf{H}_{\text{S}}^{\dagger} & \left(E\mathbf{I}-{{\mathbf{H}}_{\text{D}}}\right)\mathbf{H}_{\text{S}}^{-1} \end{array}\right]\bullet \widetilde{{\mathbf{G}}}_{N-1}^{\text{R}}(E),\end{array}\end{eqnarray} \tag{ 22a }$$

$$\begin{eqnarray}\begin{array}{*{20}{l}}={{\left[\begin{array}{c} \mathbf{0} & \mathbf{H}_{\text{S}}^{-1} \\ -\mathbf{H}_{\text{S}}^{\dagger} & \left(E\mathbf{I}-{{\mathbf{H}}_{\text{D}}}\right)\mathbf{H}_{\text{S}}^{-1} \end{array}\right]}^{N}}\bullet \widetilde{\mathbf{G}}_{0}^{\text{R}}(E),\end{array}\end{eqnarray} \tag{ 22b }$$

where this last step utilizes associativity of MMTs and we define $\widetilde{\mathbf{G}}_{0}^{\text{R}}(E)\equiv \mathbf{0}$. For notational convenience moving forward, we give this MMT the symbol

$$\begin{eqnarray}{{\mathbf{M}}_{\text{R}}}=\left[\begin{array}{c} \mathbf{0} & \mathbf{H}_{\text{S}}^{-1} \\ -\mathbf{H}_{\text{S}}^{\dagger} & \left(E\mathbf{I}-{{\mathbf{H}}_{\text{D}}}\right)\mathbf{H}_{\text{S}}^{-1} \end{array}\right].\end{eqnarray}$$

Similar to the transfer matrix, ${{\mathbf{M}}_{\text{R}}}$ depends on E.

This process can be repeated for left surface GFs. The end result is similar,

$$\begin{eqnarray}\begin{array}{*{35}{l}} \widetilde{\mathbf{G}}_{N}^{\text{L}}(E) & ={{\mathbf{M}}_{\text{L}}}\bullet \widetilde{\mathbf{G}}_{N-1}^{\text{L}}(E), \end{array}\end{eqnarray} \tag{ 23a }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & =\mathbf{M}_{\text{L}}^{N}\bullet \widetilde{\mathbf{G}}_{0}^{\text{L}}(E), \end{array}\end{eqnarray} \tag{ 23b }$$

where $\widetilde{\mathbf{G}}_{0}^{\text{L}}(E)\equiv \mathbf{0}$ and

$$\begin{eqnarray}{{\mathbf{M}}_{\text{L}}}=\left[\begin{array}{c} \mathbf{0} & \mathbf{H}_{\text{S}}^{-\dagger} \\ -{{\mathbf{H}}_{\text{S}}} & \left(E\mathbf{I}-{{\mathbf{H}}_{\text{D}}}\right)\mathbf{H}_{\text{S}}^{-\dagger} \end{array}\right].\end{eqnarray}$$

From equation (7), calculating the total system GF is tantamount to inverting the Hamiltonian. Mathematically, we thus want to consider the inverse of a block tridiagonal and block Toeplitz matrix. Theories for inverting a block tridiagonal matrix have been discussed numerous times in various contexts [169, 172–177], and very recently [179, 180] the block Toeplitz structure with them. In the end, the block quasi-separable structure of the block tridiagonal and block Toeplitz Hamiltonian leads to a structured GF, which we now formulate. As we will soon see, surface GFs play a key role, and, from the previous discussion (equations (22) and (23)), MMTs can be used to obtain all of the surface GFs. In this sense, the MMTs ${{\mathbf{M}}_{\text{L}}}$ and ${{\mathbf{M}}_{\text{R}}}$ are generators for $\mathbf{G}(E)$; that is, they embody the material's CBS.

To see this, let us first consider the diagonal blocks of the GF, ${{\mathbf{G}}_{n,n}}(E)$. A surface GF is one such block, but there many others unless the system is trivially small. As with the surface GFs above, we use Löwdin partitioning [188] to obtain an expression for ${{\mathbf{G}}_{n,n}}(E)$. Suppose ${{\mathbf{P}}_{n}}$ is the projector for our layer of interest (which is not a surface layer). Similar to before, we define ${{\mathbf{Q}}_{n}}=\mathbf{I}-{{\mathbf{P}}_{n}}$ to be the projector for all other layers of the system. A comparable statement to equation (18) can then be made where all 1 subscripts are replaced by n. Following the same logic, we obtain an effective Hamiltonian for layer n and ultimately arrive at

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\mathbf{G}}_{n,n}}(E) & ={{\left[E{{\mathbf{P}}_{n}}-{{\mathbf{P}}_{n}}\mathbf{H}{{\mathbf{P}}_{n}}-{{\mathbf{P}}_{n}}\mathbf{H}{{\mathbf{Q}}_{n}}{{\mathbf{G}}_{{{\mathbf{Q}}_{n}}}}(E){{\mathbf{Q}}_{n}}\mathbf{H}{{\mathbf{P}}_{n}}\right]}^{-1}} \\ {} & ={{\left[E\mathbf{I}-{{\mathbf{H}}_{\text{D}}}-{{\mathbf{P}}_{n}}\mathbf{H}{{\mathbf{Q}}_{n}}{{\mathbf{G}}_{{{\mathbf{Q}}_{n}}}}(E){{\mathbf{Q}}_{n}}\mathbf{H}{{\mathbf{P}}_{n}}\right]}^{-1}}, \end{array}\end{eqnarray}$$

where we have again replaced ${{\mathbf{P}}_{n}}$ with the identity in the energy term.

The key difference here is that there are two disconnected parts of the ${{\mathbf{Q}}_{n}}$ subsystem. Layers $1,2,\ldots,n-1$ are one set, and layers $n+1,n+2,\ldots$ are the other. Because layer n only couples to layers n  −  1 and n  +  1 in the total system,

$$\begin{eqnarray}&&{{\mathbf{P}}_{n}}\mathbf{H}{{\mathbf{Q}}_{n}}={{\mathbf{H}}_{n,n-1}}{{\mathbf{P}}_{n-1}}+{{\mathbf{H}}_{n,n+1}}{{\mathbf{P}}_{n+1}},\end{eqnarray}$$

and similarly for ${{\mathbf{Q}}_{n}}\mathbf{H}{{\mathbf{P}}_{n}}$. Thus,

$$\begin{eqnarray}&&\begin{array}{*{20}{l}}{{\mathbf{G}}_{n,n}}(E)=\left[E\mathbf{I}-{{\mathbf{H}}_{\text{D}}}\,-\,\left({{\mathbf{H}}_{n,n-1}}{{\mathbf{P}}_{n-1}}+{{\mathbf{H}}_{n,n+1}}{{\mathbf{P}}_{n+1}}\right)\right. \\ \left. {{\mathbf{G}}_{{{\mathbf{Q}}_{n}}}}(E)\left({{\mathbf{P}}_{n-1}}{{\mathbf{H}}_{n-1,n}}+{{\mathbf{P}}_{n+1}}{{\mathbf{H}}_{n+1,n}}\right)\right]^{-1}.\end{array}\end{eqnarray}$$

After distributing, each of these '$\mathbf{PGP}$' terms simplifies. First (and second), ${{\mathbf{P}}_{n\pm 1}}{{\mathbf{G}}_{{{\mathbf{Q}}_{n}}}}(E){{\mathbf{P}}_{n\mp 1}}=\mathbf{0}$ because layers n  −  1 and n  +  1 do not belong to the same part of the ${{\mathbf{Q}}_{n}}$ subsystem (i.e. they are decoupled in the subsystem). Third, layer n  −  1 is a surface layer for the subsystem containing layers $1,2,\ldots,n-1$; therefore (as in the surface GF discussion), ${{\mathbf{P}}_{n-1}}{{\mathbf{G}}_{{{\mathbf{Q}}_{n}}}}(E){{\mathbf{P}}_{n-1}}=\widetilde{\mathbf{G}}_{n-1}^{\text{L}}(E)$. Fourth, and similarly, layer n  +  1 is a surface layer for the subsystem containing layers $n+1,n+2,\ldots$ such that ${{\mathbf{P}}_{n+1}}{{\mathbf{G}}_{{{\mathbf{Q}}_{n}}}}(E){{\mathbf{P}}_{n+1}}=\widetilde{\mathbf{G}}_{N-n}^{\text{R}}(E)$, where we assume there are N layers in the total system. Putting these results together,

$$\begin{eqnarray}&&{{\mathbf{G}}_{n,n}}(E)={{[E\mathbf{I}-{{\mathbf{H}}_{\text{D}}}-{{\mathbf{H}}_{\text{S}}}\widetilde{{\mathbf{G}}}_{n-1}^{\text{L}}(E)\mathbf{H}_{\text{S}}^{\dagger}-\mathbf{H}_{\text{S}}^{\dagger}\widetilde{{\mathbf{G}}}_{N-n}^{\text{R}}(E){{\mathbf{H}}_{\text{S}}}]}^{-1}}^{{}}.\end{eqnarray} \tag{ 24 }$$

Using our convention that $\widetilde{\mathbf{G}}_{0}^{\text{L}/\text{R}}(E)=\mathbf{0}$, this equation also handles the cases when n labels a surface layer; compare to equation (21).

It is clear from equation (24) that calculating a diagonal block of the GF essentially reduces to calculating surface GFs. Note also that the appropriate surface GF(s) should be replaced by $\widetilde{\mathbf{G}}_{\infty}^{\text{L}/\text{R}}(E)$ if the material has an infinite or semi-infinite number of layers.

Tying this result back to our running example, figure 4 shows projected DOSs for various layers of our model system, using both choices of layers from figures 1 and 2. Each layer's projected DOS is an instance of equation (13),

$$\begin{eqnarray}&&{{\rho}_{n}}(E)=\frac{-1}{\pi}\text{Im}(\text{Tr}[{{\mathbf{G}}_{n,n}}(E)]),\end{eqnarray} \tag{ 25 }$$

and has been used [15, 27, 39, 134–136, 142, 149] to study quantum size effects. The projected DOS for layer 1 (i.e. the surface layer) is, as expected, the surface DOS for the semi-infinite system. As we move away from the surface, oscillations appear in the projected DOSs, which eventually converge to the bulk DOS (essentially layer $\infty$) at a rate determined by CBS (point (VII.11) below). The van Hove singularities [191] at the band edges of the bulk DOS signify of a one-dimensional material. Although the choice of layer is critically important for surface-related properties, it does not matter for bulk properties.

Zoom In Zoom Out Reset image size

Let us now consider the off-diagonal blocks of the GF. Expressions for these blocks can be straightforwardly obtained from Dyson's equation [168], which, similar to Löwdin partitioning, has us break up the system. Unlike Löwdin partitioning, one of our 'partitions' will only contain the coupling between two layers; the other 'partition' contains everything else. The idea is best presented through example. Suppose we want to calculate ${{\mathbf{G}}_{n+1,n}}(E)$. We regard the coupling between layers n and n  +  1,

$$\begin{eqnarray}&&\mathbf{V}={{\mathbf{P}}_{n}}\mathbf{H}{{\mathbf{P}}_{n+1}}+{{\mathbf{P}}_{n+1}}\mathbf{H}{{\mathbf{P}}_{n}},\end{eqnarray}$$

as a perturbation on the Hamiltonian ${{\mathbf{H}}_{0}}$, which describes layers $1,2,\ldots,n$ decoupled from layers $n+1,n+2,\ldots$. In other words, ${{\mathbf{H}}_{0}}=\mathbf{H}-\mathbf{V}$. Dyson's equation relates the GF for the total system ($\mathbf{G}$) to the GF for the unperturbed system (${{\mathbf{G}}_{0}}$),

$$\begin{eqnarray}&&\mathbf{G}(E)={{\mathbf{G}}_{0}}(E)+\mathbf{G}(E)\mathbf{V}{{\mathbf{G}}_{0}}(E).\end{eqnarray} \tag{ 26 }$$

Applying this to the (n  +  1, n) block,

$$\begin{eqnarray}&&\begin{array}{*{20}{l}}{{\mathbf{G}}_{n+1,n}}(E)={{\mathbf{P}}_{n+1}}{{\mathbf{G}}_{0}}(E){{\mathbf{P}}_{n}}+{{\mathbf{P}}_{n+1}}\mathbf{G}(E)\left({{\mathbf{P}}_{n}}\mathbf{H}{{\mathbf{P}}_{n+1}}\right. \\ \left. \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,+{{\mathbf{P}}_{n+1}}\mathbf{H}{{\mathbf{P}}_{n}}\right){{\mathbf{G}}_{0}}(E){{\mathbf{P}}_{n}}.\end{array}\end{eqnarray}$$

Because layers n and n  +  1 are decoupled in ${{\mathbf{H}}_{0}}$, ${{\mathbf{P}}_{n+1}}{{\mathbf{G}}_{0}}(E){{\mathbf{P}}_{n}}=\mathbf{0}$. Layer n is also a surface layer in ${{\mathbf{H}}_{0}}$ such that ${{\mathbf{P}}_{n}}{{\mathbf{G}}_{0}}(E){{\mathbf{P}}_{n}}=\widetilde{\mathbf{G}}_{n}^{\text{L}}(E)$. Thus,

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\mathbf{G}}_{n+1,n}}(E) & ={{\mathbf{P}}_{n+1}}\mathbf{G}(E){{\mathbf{P}}_{n+1}}\mathbf{H}{{\mathbf{P}}_{n}}{{\mathbf{G}}_{0}}(E){{\mathbf{P}}_{n}} \\ {} & ={{\mathbf{G}}_{n+1,n+1}}(E){{\mathbf{H}}_{\text{S}}}\widetilde{\mathbf{G}}_{n}^{\text{L}}(E). \end{array}\end{eqnarray}$$

This procedure of isolating the coupling between two layers and applying Dyson's equation can be iterated, leading to the general result

$$\begin{eqnarray}&&{{\mathbf{G}}_{m,n}}(E)={{\mathbf{G}}_{m,n+1}}(E){{\mathbf{H}}_{\text{S}}}\widetilde{\mathbf{G}}_{n}^{\text{L}}(E)~\text{when}~m>n.\end{eqnarray} \tag{ 27a }$$

Similar logic can be used when m  <  n, and results in

$$\begin{eqnarray}&&{{\mathbf{G}}_{m,n}}(E)={{\mathbf{G}}_{m,n-1}}(E)\mathbf{H}_{\text{S}}^{\dagger}\widetilde{\mathbf{G}}_{N-n+1}^{\text{R}}(E)~\text{when}~m<n.\end{eqnarray} \tag{ 27b }$$

Finally, Dyson's equation can equivalently be written as $\mathbf{G}(E)={{\mathbf{G}}_{0}}(E)+{{\mathbf{G}}_{0}}(E)\mathbf{VG}(E)$, which produces the equivalent expressions

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\mathbf{G}}_{m,n}}(E) & =\widetilde{\mathbf{G}}_{N-m+1}^{\text{R}}(E){{\mathbf{H}}_{\text{S}}}{{\mathbf{G}}_{m-1,n}}(E)~\text{when}~m>n, \end{array}\end{eqnarray} \tag{ 27c }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\mathbf{G}}_{m,n}}(E) & =\widetilde{\mathbf{G}}_{m}^{\text{L}}(E)\mathbf{H}_{\text{S}}^{\dagger}{{\mathbf{G}}_{m+1,n}}(E)~\text{when}~m<n. \end{array}\end{eqnarray} \tag{ 27d  }$$

These equations can also be extended to infinite or semi-infinite systems by replacing the appropriate surface GFs with $\widetilde{\mathbf{G}}_{\infty}^{\text{L}/\text{R}}(E)$.

Equations (24) and (27) demonstrate that every block of the GF can be straightforwardly obtained from surface GFs. Then, because the surface GFs are described by the MMTs ${{\mathbf{M}}_{\text{L}}}$ and ${{\mathbf{M}}_{\text{R}}}$, these MMTs transitively generate the entire GF [180].

This last statement suggests that ${{\mathbf{M}}_{\text{L}}}$ and ${{\mathbf{M}}_{\text{R}}}$ access the material's CBS in a way that facilitates matrix inversion; that is, ${{\mathbf{M}}_{\text{L}}}$ and ${{\mathbf{M}}_{\text{R}}}$ provide a GF-based route to CBS. We demonstrate this by showing some common structure between these MMTs and the transfer matrix $\mathbf{T}$ in equation (11). Mathematically, it is straightforward to prove that ${{\mathbf{M}}_{\text{L}}}$ is similar to $\mathbf{T}$,

$$\begin{eqnarray}{{\mathbf{M}}_{\text{L}}}=\left[\begin{array}{c} \mathbf{0} & \mathbf{I} \\ \mathbf{H}_{\text{S}}^{\dagger} & \mathbf{0} \end{array}\right]\mathbf{T}\left[\begin{array}{c} \mathbf{0} & \mathbf{H}_{\text{S}}^{-\dagger} \\ \mathbf{I} & \mathbf{0} \end{array}\right];\end{eqnarray}$$

meaning that ${{\mathbf{M}}_{\text{L}}}$ has the same eigenvalues as $\mathbf{T}$. From section 4, these eigenvalues are essentially the material's CBS. Likewise, ${{\mathbf{M}}_{\text{R}}}$ is similar to ${{\mathbf{T}}^{-1}}$,

$$\begin{eqnarray}{{\mathbf{M}}_{\text{R}}}=\left[\begin{array}{c} \mathbf{I} & \mathbf{0} \\ \mathbf{0} & {{\mathbf{H}}_{\text{S}}} \end{array}\right]{{\mathbf{T}}^{-1}}\left[\begin{array}{c} \mathbf{I} & \mathbf{0} \\ \mathbf{0} & \mathbf{H}_{\text{S}}^{-1} \end{array}\right].\end{eqnarray}$$

In other words, $\lambda ={{\text{e}}^{\text{i}ka}}$ is an eigenvalue of $\mathbf{T}$ (or ${{\mathbf{M}}_{\text{L}}}$) if and only if $1/\lambda ={{\text{e}}^{-\text{i}ka}}$ is an eigenvalue of ${{\mathbf{M}}_{\text{R}}}$. For reference,

$$\begin{eqnarray}{{\mathbf{T}}^{-1}}=\left[\begin{array}{c} \mathbf{0} & \mathbf{I} \\ -\mathbf{H}_{\text{S}}^{-1}\mathbf{H}_{\text{S}}^{\dagger} & \mathbf{H}_{\text{S}}^{-1}\left(E\mathbf{I}-{{\mathbf{H}}_{\text{D}}}\right) \end{array}\right].\end{eqnarray}$$

This relationship between $\mathbf{T}$, ${{\mathbf{M}}_{\text{L}}}$, and ${{\mathbf{M}}_{\text{R}}}$ is the multi-band generalization of a likewise statement for one-band materials, where $\mathbf{T}$, ${{\mathbf{M}}_{\text{L}}}$, and ${{\mathbf{M}}_{\text{R}}}$ are $2\times 2$ matrices and complex Möbius transformations were used [156].

There are a handful of properties that apply regardless of the material's dimensionality. More detailed discussions of these results can be found in [4, 6, 7, 163].

• (VI.1)  
  The dispersion relation $E(\vec{k})$ is analytic [162] except at $\vec{k}$ where multiple bands intersect (i.e. branch points). Proofs are presented in [6, 7, 52, 163].
• (VI.2)  
  If $\vec{g}$ is a reciprocal lattice vector (and thus real), then $E(\vec{k})=E(\vec{k}+\vec{g})$. This result is easily verified by using Bloch's theorem in the time-independent Schrödinger equation.
• (VI.3)  
  $E(\vec{k})=E{{({{\vec{k}}^{\ast}})}^{\ast}}$, where ^* denotes complex conjugation. Similar to the previous point, this result can also be verified using Bloch's theorem in the time-independent Schrödinger equation. Furthermore, if E is real, then both $\vec{k}$ and ${{\vec{k}}^{\ast}}$ are in the material's CBS. This symmetry between $\vec{k}$ and ${{\vec{k}}^{\ast}}$ can be observed in figures 1 and 2 from our example.
• (VI.4)  
  Branch points occur at non-real $\vec{k}$ when the potential (and thus Hamiltonian) are Hermitian [7]. However, branch points can be shifted to real $\vec{k}$ if the potential is non-Hermitian, as might happen in (for example) high-energy electron diffraction [68].

Heine [4] derived several results about the 'lines of real E' when $\vec{k}$ is only complex in one dimension that were subsequently classified by Chang [32]. These results are all visible in our example from section 2 (figures 1 and 2), which may be a useful companion.

• (VI.5)  
  Away from branch points and band edges, lines of real energy do not merge, split, begin, or terminate. Heine [4] demonstrated these results by considering the winding number [162] of the function E(k)  −  E(k₀) for some wavevector k₀ that is not a branch point or band edge. There is only one zero of this function inside some infinitesimally-small, closed contour around k₀, meaning that the function accumulates a phase of $2\pi$ going around the contour. Consequently, there are two points with real energies on the contour: One for the line to enter and one for it to leave. If a line merged, split, began, or terminated, there would have to be an odd number of real energies on the contour.
• (VI.6)  
  At a band edge, two lines of real energy meet and take ${{90}^{\circ}}$ turns into the complex plane. A proof is straightforward. Suppose k₀ is a band edge. E(k) has an extremum along the real axis at k₀, meaning ${{E}^{\prime}}\left({{k}_{0}}\right)=0$. Then, because E(k) is analytic at k₀, we know that

  $$\begin{eqnarray}\begin{array}{*{35}{l}} E(k) & =E({{k}_{0}})+{{E}^{\prime}}({{k}_{0}})(k-{{k}_{0}})+\frac{{{E}^{\prime \prime}}({{k}_{0}})}{2}{{(k-{{k}_{0}})}^{2}}+\mathcal{O}[{{(k-{{k}_{0}})}^{3}}] \\ {} & \approx E({{k}_{0}})+\frac{{{E}^{\prime \prime}}({{k}_{0}})}{2}{{(k-{{k}_{0}})}^{2}} \end{array}\end{eqnarray}$$

  for k near k₀. Because k₀ is in a real band, E(k₀) and ${{E}^{\prime \prime}}({{k}_{0}})$ are both real, implying that ${{(k-{{k}_{0}})}^{2}}$ must be real for E(k) to be real. Thus, k  −  k₀ is either real (meaning k is also in the conventional band) or purely imaginary (corresponding to a ${{90}^{\circ}}$ turn). From this reasoning, we also see that a band edge k₀ is a saddle point of E(k).
• (VI.7)  
  Lines of real energy do not generally go through branch points. In other words, if k_bp is a branch point, $E({{k}_{\text{bp}}})$ is not usually real. See the example in figure 2. The one major exception [4] to this statement occurs when there is a mirror plane between between the layers, as demonstrated in figure 1. $E({{k}_{\text{bp}}})$ will be real in this case. Either way [4],

  $$\begin{eqnarray}&&E(k)\approx E({{k}_{\text{bp}}})+C{{(k-{{k}_{\text{bp}}})}^{1/2}},\end{eqnarray} \tag{ 28 }$$

  for k near k_bp and some (generally complex) constant C. Note that equation (28) is only valid in the likely case that two bands meet at the branch point (see point (VII.2)). If m bands were to intersect at k_bp, the exponent 1/2 would be replaced by 1/m [6].

Because one-dimensional CBS is related to the eigenvalues of the transfer matrix $\mathbf{T}$ (alternatively and equivalently, the MMTs ${{\mathbf{M}}_{\text{L}}}$ and ${{\mathbf{M}}_{\text{R}}}$), we focus on the spectral decomposition of $\mathbf{T}$.

• (VI.8)  
  The transfer matrix is not Hermitian. As desired, this means that $\mathbf{T}$ can have complex eigenvalues.
• (VI.9)  
  The eigenvalues of $\mathbf{T}$ come in pairs: If λ is an eigenvalue, so is $1/{{\lambda}^{\ast}}$. To see this, we verify that $\mathbf{T}$ is similar to ${{\mathbf{T}}^{-\dagger}}$ (mathematically, $\mathbf{T}$ is complex symplectic),

  $$\begin{eqnarray}\mathbf{T}=\left[\begin{array}{c} \mathbf{0} & -\mathbf{H}_{\text{S}}^{-\dagger} \\ \mathbf{H}_{\text{S}}^{-1} & \mathbf{0} \end{array}\right]{{\mathbf{T}}^{-\dagger}}\left[\begin{array}{c} \mathbf{0} & {{\mathbf{H}}_{\text{S}}} \\ -\mathbf{H}_{\text{S}}^{\dagger} & \mathbf{0} \end{array}\right].\end{eqnarray}$$
• (VI.10)  
  $\mathbf{T}$ has an even (possibly zero) number of eigenvalues with unit magnitude, and an equal number of eigenvalues with magnitudes greater than and less than 1. As follows from fact 2.14.13 of [192],

  $$\begin{eqnarray}&&\text{det}(\mathbf{T})=\text{det}(\mathbf{H}_{\text{S}}^{-\dagger}{{\mathbf{H}}_{\text{S}}})=\frac{\text{det}({{\mathbf{H}}_{\text{S}}})}{\text{det}{{({{\mathbf{H}}_{\text{S}}})}^{\ast}}}.\end{eqnarray}$$

  Thus, $|\text{det}(\mathbf{T})|=1$, meaning the product of all eigenvalues is 1 (in magnitude). Appealing to the eigenvalue pairing in point (VI.9), it is straightforward to see the conclusion. Suppose that λ is an eigenvalue of $\mathbf{T}$ with unit magnitude. Then, $1/{{\lambda}^{\ast}}$ has unit magnitude and is also an eigenvalue. Similar reasoning reveals that there are an equal number of eigenvalues with $|\lambda |<1$ and $|\lambda |>1$.
• (VI.11)  
  Eigenvalues with unit magnitude (see points (VII.4) and (VII.5) below) are forbidden when E is not real. This situation arises, for example, when numerically evaluating the limit in equation (7) and a small imaginary component is present in the energy. In essence, the imaginary energy breaks symmetry so that the paired eigenvalues are no longer degenerate. Half of the states will have $|\lambda |<1$ and the other half will have $|\lambda |>1$ [187]. Note that, as expected, $|\lambda |\to 1$ as $\text{Im}(E)\to {{0}^{\pm}}$ for these eigenvalues.A proof of this result is straightforward and proceeds by contradiction. Suppose that λ is an eigenvalue of $\mathbf{T}$ ($|\lambda |=1$) when $\text{Im}(E)\ne 0$. Then, from equation (11) and fact 2.14.13 in [192],

  $$\begin{eqnarray}\begin{array}{*{20}{l}}0&=&\text{det}\left[\lambda \mathbf{I}-\mathbf{T}\right]\\ &=&\text{det}\left[\begin{array}{c} \lambda \mathbf{I}-\mathbf{H}_{\text{S}}^{-\dagger}\left(E\mathbf{I}-{{\mathbf{H}}_{\text{D}}}\right) & \mathbf{H}_{\text{S}}^{-\dagger}{{\mathbf{H}}_{\text{S}}} \\ -\mathbf{I} & \lambda \mathbf{I} \end{array}\right]\\ &=&\text{det}[-\lambda \mathbf{H}_{\text{S}}^{-\dagger}]\text{det}\left[E\mathbf{I}-\lambda \mathbf{H}_{\text{S}}^{\dagger}-{{\mathbf{H}}_{\text{D}}}-\frac{1}{\lambda}{{\mathbf{H}}_{\text{S}}}\right].\end{array}\end{eqnarray}$$

  The first determinant is nonzero because $|\lambda |\,=1$ and $\mathbf{H}_{\text{S}}^{-\dagger}$ is invertible. Considering the second determinant, we can write $\lambda ={{\text{e}}^{\text{i}\theta}}$ (θ is real) such that

  $$\begin{eqnarray}&&0=\text{det}[E\mathbf{I}-{{\text{e}}^{\text{i}\theta}}\mathbf{H}_{\text{S}}^{\dagger}-{{\mathbf{H}}_{\text{D}}}-{{\text{e}}^{-\text{i}\theta}}{{\mathbf{H}}_{\text{S}}}].\end{eqnarray} \tag{ 29 }$$

  Recalling that ${{\mathbf{H}}_{\text{D}}}$ is Hermitian, it is easily verified that ${{\text{e}}^{\text{i}\theta}}\mathbf{H}_{\text{S}}^{\dagger}+{{\mathbf{H}}_{\text{D}}}+{{\text{e}}^{-\text{i}\theta}}{{\mathbf{H}}_{\text{S}}}$ is also Hermitian. Then, because equation (29) can only hold when E is an eigenvalue of ${{\text{e}}^{\text{i}\theta}}\mathbf{H}_{\text{S}}^{\dagger}+{{\mathbf{H}}_{\text{D}}}+{{\text{e}}^{-\text{i}\theta}}{{\mathbf{H}}_{\text{S}}}$, we reach a contradiction. Hermitian operators must have real eigenvalues and E is not real. Thus, eigenvalues of $\mathbf{T}$ cannot have unit magnitude when E is not real.
• (VI.12)  
  $\mathbf{T}$ may not be diagonalizable. In general, the transfer matrix is not normal, meaning $\mathbf{T}{{\mathbf{T}}^{\dagger}}\ne {{\mathbf{T}}^{\dagger}}\mathbf{T}$. Consequently, $\mathbf{T}$ is not guaranteed to have a complete set of eigenvectors and any eigenvectors it does possess are not required to be mutually orthogonal [193]. Rather than diagonalizing $\mathbf{T}$, a Jordan or Schur decomposition [193] (which both generalize the notion of 'diagonalization') may be more appropriate for analytical or numerical considerations, respectively. Additional structure of the eigenvectors can be found in [154].

• (VII.1)  
  Each branch of the multi-valued dispersion relation corresponds to a band of the material. In this context, a 'band' includes the real band from conventional band structure, as well as some complex wavevectors. Prodan [7] further discusses, with examples of the various Riemann surfaces, branches and branch points of E(k) when k is strictly one-dimensional.
• (VII.2)  
  Branch points of $E(\vec{k})$ occur at the boundaries between two bands. Consequently, they will only appear in band gaps. It may be possible that branch points will appear inside a band for higher-dimensional (not one-dimensional) materials [45, 47]; however, in the absence of general tools for higher-dimensional CBS, heuristics were used to calculate those branch point energies³.
• (VII.3)  
  If ${{\vec{k}}_{\text{bp}}}$ is a branch point (and therefore in a band gap from point (VII.2)), states near ${{\vec{k}}_{\text{bp}}}$ with $E<E({{\vec{k}}_{\text{bp}}})$ will have more 'valence band character', and similarly for 'conduction band character' when $E>E({{\vec{k}}_{\text{bp}}})$ [19, 20, 55]. Using chemistry nomenclature, states near the branch point with $E=E({{\vec{k}}_{\text{bp}}})$ are 'non-bonding'.

• (VII.4)  
  Eigenvalues of $\mathbf{T}$ with unit magnitude correspond to the bulk-propagating states from conventional band structure. Briefly, if $\lambda ={{\text{e}}^{\text{i}ka}}$ and $|\lambda |\,=1$, then $\text{Im}(k)=0$. Because k is real, λ and $1/{{\lambda}^{\ast}}$ are degenerate eigenvalues. If E is outside of a band, there will not be any of these states or eigenvalues.
• (VII.5)  
  Eigenvalues of $\mathbf{T}$ with non-unit magnitudes correspond to the evanescent states that are excluded from conventional band structure. Depending on orientation, states corresponding to $|\lambda |<1$ grow or decay from one layer to the next, whereas states with $|\lambda |>1$ exhibit the opposite behavior. From point (VI.10), there are an equal number of states with $|\lambda |>1$ and $|\lambda |<1$.
• (VII.6)  
  The states corresponding to defective eigenvalues of $\mathbf{T}$ (i.e. the degenerate eigenvalues that lack a complete set of eigenvectors, thereby prohibiting $\mathbf{T}$ from being diagonalized) are not well understood. Some results [18, 68, 136, 147, 149, 150, 154, 163, 194–196] suggest that $\mathbf{T}$ is defective at a band edge or a branch point, or when a surface state exists (point (VII.10)). More work needs to be performed to understand the significance of defective $\mathbf{T}$.
• (VII.7)  
  CBS determines the dynamic properties of electrons in a material [52], which are described by the off-diagonal blocks of the GF. To see this, we apply general properties of MMTs [180]. Equation (27c) takes a lower-triangular block of the GF and moves one layer down, and can be approximated as

  $$\begin{eqnarray}&&{{\mathbf{G}}_{j,k}}(E)\approx \widetilde{\mathbf{G}}_{\infty}^{\text{R}}(E){{\mathbf{H}}_{\text{S}}}{{\mathbf{G}}_{j-1,k}}(E),\end{eqnarray}$$

  when j is sufficiently large (see point (VII.9) below). Then, using equation (30),

  $$\begin{eqnarray}\begin{array}{*{35}{l}} {{\mathbf{G}}_{j,k}}(E) & = & \left({{\mathbf{M}}_{\text{R}}}\bullet \widetilde{{\mathbf{G}}}_{\infty}^{\text{R}}(E)\right){{\mathbf{H}}_{\text{S}}}{{\mathbf{G}}_{j-1,k}}(E) \\ {} & = & \mathbf{H}_{\text{S}}^{-1}{{\left[-\mathbf{H}_{\text{S}}^{\dagger}\widetilde{{\mathbf{G}}}_{\infty}^{\text{R}}(E)+\left(E\mathbf{I}-{{\mathbf{H}}_{\text{D}}}\right)\mathbf{H}_{\text{S}}^{-1}\right]}^{-1}}{{\mathbf{H}}_{\text{S}}}{{\mathbf{G}}_{j-1,k}}(E). \end{array}\end{eqnarray}$$

  We now invoke a general property of MMTs (theorem 2 of [180]): The eigenvalues of $-\mathbf{H}_{\text{S}}^{\dagger}\widetilde{\mathbf{G}}_{\infty}^{\text{R}}(E)+\left(E\mathbf{I}-{{\mathbf{H}}_{\text{D}}}\right)\mathbf{H}_{\text{S}}^{-1}$ are exactly the eigenvalues of ${{\mathbf{M}}_{\text{R}}}$ with $|\lambda |>1$ (assuming $\text{Im}(E)\ne 0$). In this sense, the material's CBS determines the behavior of the GF as we move off the diagonal; these are the dynamical properties of electrons. If E is outside of a band, the electron transition probabilities fall off rapidly with each additional layer because there are no propagating states. Conversely, if E is in a band, there is at least one eigenvalue satisfying $|\lambda |\approx 1$, and the off-diagonal blocks of the GF are slow to decay. As a final comment, a similar analysis can be performed for each part of equation (27), and leads to the same result.

Surfaces were one of the first applications of CBS. Notably, and as discussed in section 5, CBS is closely linked to the surface GF for a semi-infinite system [12, 26, 29, 30, 50, 123, 135, 159, 187]. Herein we describe some additional results related to surfaces.

• (VII.8)  
  The surface GF (for a semi-infinite system) can be directly obtained from the eigenvectors of the transfer matrix [29] or the MMTs ${{\mathbf{M}}_{\text{L}}}$ and ${{\mathbf{M}}_{\text{R}}}$ [135, 159]; see equation (31). (Generalized eigenvectors may also be needed if the transfer matrix is not diagonalizable.) This explains why the surface GFs in figure 3 are different for the two choices of layers. Even though the eigenvalues of ${{\mathbf{M}}_{\text{R}}}$, ${{\mathbf{M}}_{\text{L}}}$, or $\mathbf{T}$ (i.e. the dispersion relation) are mostly independent of the choice of layer (see figures 1 and 2), the eigenvectors are much more sensitive.To elaborate, consider a semi-infinite material with (right) surface GF $\widetilde{\mathbf{G}}_{\infty}^{\text{R}}(E)$. When applied to a surface GF, the MMT ${{\mathbf{M}}_{\text{R}}}$ produces the surface GF for a system with an additional layer. This new system, however, is identical to the original semi-infinite system; thus,

  $$\begin{eqnarray}&&\widetilde{\mathbf{G}}_{\infty}^{\text{R}}(E)={{\mathbf{M}}_{\text{R}}}\bullet \widetilde{\mathbf{G}}_{\infty}^{\text{R}}(E).\end{eqnarray} \tag{ 30 }$$

  Mathematically, the surface GF for a semi-infinite system is a fixed point of ${{\mathbf{M}}_{\text{R}}}$.We previously showed [180] that the fixed points of a MMT are related to the MMT's invariant subspaces [197]. In other words, a fixed point corresponds to a collection of eigenvectors of the MMT (and possibly generalized eigenvectors if the MMT is not diagonalizable). Physically, we are interested in the half of the eigenvectors that correspond to right- (or left-)propagating states; that is, those with eigenvalues greater than (or less than) 1 in magnitude. Note that, if there are eigenvalues with unit magnitude (i.e. E is in a band) such that half of the eigenvalues are not greater (less) than 1 in magnitude (point (VI.9)), it is well known that the surface GF is not well-defined. Adding a small imaginary part to the energy (point (VI.11)) eliminates this problem.A proof of the upcoming result is lengthy and requires other properties of MMTs that we will not discuss here. For brevity, we simply state the result and direct the interested reader to [180] for more details. We also have to impose a basis set and assume that ${{\mathbf{H}}_{\text{D}}}$ and ${{\mathbf{H}}_{\text{S}}}$ are accurately represented by $M\times M$ matrices⁴. Let $\mathbf{U}$ be the $2M\times M$ matrix whose columns are the M (generalized) eigenvectors of ${{\mathbf{M}}_{\text{R}}}$ with eigenvalues greater than 1 in magnitude (recall that ${{\mathbf{M}}_{\text{R}}}$ will be a $2M\times 2M$ matrix). Then,

  $$\begin{eqnarray}&&\widetilde{\mathbf{G}}_{\infty}^{\text{R}}(E)={{\mathbf{U}}_{1}}\mathbf{U}_{2}^{-1},\end{eqnarray} \tag{ 31 }$$

  where ${{\mathbf{U}}_{1}}$ is the first M rows of $\mathbf{U}$ and ${{\mathbf{U}}_{2}}$ is the bottom M rows of $\mathbf{U}$. Likewise, if $\mathbf{V}$ is the $2M\times M$ matrix of eigenvectors of ${{\mathbf{M}}_{\text{L}}}$ with eigenvalues greater than 1 in magnitude, then

  $$\begin{eqnarray}&&\widetilde{\mathbf{G}}_{\infty}^{\text{L}}(E)={{\mathbf{V}}_{1}}\mathbf{V}_{2}^{-1},\end{eqnarray}$$

  where ${{\mathbf{V}}_{1}}$ and ${{\mathbf{V}}_{2}}$ are similarly defined.
• (VII.9)  
  CBS determines the rate at which $\widetilde{\mathbf{G}}_{N}^{\text{L}/\text{R}}(E)$ converges to $\widetilde{\mathbf{G}}_{\infty}^{\text{L}/\text{R}}(E)$ with increasing N (if it converges), which is part of an examination of quantum size effects. This result is also relevant to numerical methods that approximate $\widetilde{\mathbf{G}}_{\infty}^{\text{L}/\text{R}}(E)$ by $\widetilde{\mathbf{G}}_{N}^{\text{L}/\text{R}}(E)$ for N sufficiently large [186, 187]. Specifically, $\widetilde{\mathbf{G}}_{N}^{\text{L}/\text{R}}(E)$ converges to $\widetilde{\mathbf{G}}_{\infty}^{\text{L}/\text{R}}(E)$ as $\mathcal{O}(|{{\lambda}_{M}}{{|}^{2N}})$, where ${{\lambda}_{M}}$ is the largest (in magnitude) eigenvalue of ${{\mathbf{M}}_{\text{L}/\text{R}}}$ satisfying $|{{\lambda}_{M}}|\leqslant 1$. As implied, $\widetilde{\mathbf{G}}_{N}^{\text{L}/\text{R}}(E)$ does not converge to the semi-infinite limit when $|{{\lambda}_{M}}|=1$. In this case, E is in a band and $\widetilde{\mathbf{G}}_{\infty}^{\text{L}/\text{R}}(E)$ is not well-defined.Starting from equation (22b) and using the Jordan form of ${{\mathbf{M}}_{\text{R}}}={{\mathbf{U}}_{\text{R}}}\mathbf{JU}_{\text{R}}^{-1}$,

  $$\begin{eqnarray}&&\widetilde{\mathbf{G}}_{N}^{\text{R}}(E)=({{\mathbf{U}}_{\text{R}}}{{\mathbf{J}}^{N}}\mathbf{U}_{\text{R}}^{-1})\bullet \widetilde{\mathbf{G}}_{0}^{\text{R}}(E).\end{eqnarray}$$

  Similar to the discussion in point (VII.8), we assume that $\text{Im}(E)>0$ such that half of the eigenvalues in $\mathbf{J}$ are greater than 1 (in magnitude) and the other half are less than 1 (in magnitude); see points (VI.10) and (VI.11). We can then order the eigenvalues in $\mathbf{J}$ so that each element of $\mathbf{U}_{\text{R}}^{-1}\bullet \widetilde{\mathbf{G}}_{0}^{\text{R}}(E)$ is multiplied by an eigenvalue that is smaller than 1 (in magnitude) and simultaneously divided by an eigenvalue that is larger than 1 (in magnitude) when ${{\mathbf{J}}^{N}}$ is applied. Consequently, each element of $({{\mathbf{J}}^{N}}\mathbf{U}_{\text{R}}^{-1})\bullet \widetilde{\mathbf{G}}_{0}^{\text{R}}(E)$ exponentially converges to 0 with a rate determined by CBS. The slowest element to approach zero converges as $\mathcal{O}(|{{\lambda}_{M}}{{|}^{2N}})$.
• (VII.10)  
  CBS and the (semi-infinite) surface GF detail the existence of surface states and their asymptotic decay rates. Surface states appear at energies where there are no bulk-propagating states (that is, no λ have unit magnitude) and $\widetilde{\mathbf{G}}_{\infty}^{\text{L}/\text{R}}(E)$ has a pole [26, 28]. The latter condition is satisfied when, mathematically (with a momentary abuse of notation),

  $$\begin{eqnarray}&&0=\text{det}[{{(\widetilde{\mathbf{G}}_{\infty}^{\text{L}/\text{R}}(E))}^{-1}}].\end{eqnarray}$$

  From equation (31), this condition rigorously becomes

  $$\begin{eqnarray}&&0=\text{det}({{\mathbf{U}}_{2}})~\text{or}~0=\text{det}({{\mathbf{V}}_{2}}),\end{eqnarray} \tag{ 32 }$$

  with ${{\mathbf{U}}_{2}}$ and ${{\mathbf{V}}_{2}}$ defined as in point (VII.8). Last, if there exists a surface state at energy E, the state asymptotically decays (in magnitude) at a rate of $-\ln |{{\lambda}_{M}}|$ per layer [29], where ${{\lambda}_{M}}$ is the largest eigenvalue (in magnitude) satisfying $|{{\lambda}_{M}}|<1$. Physically, ${{\lambda}_{M}}$ corresponds to the most slowly decaying evanescent state from the material's CBS.
• (VII.11)  
  CBS describes the general decay of surface effects [135, 136], which helps explore quantum size effects. As depicted in figure 4, the diagonal blocks of the GF (which provide access to the projected DOSs) are one way to investigate this behavior. From equation (24), diagonal blocks of the GF are easily calculated from surface GFs for systems with the appropriate numbers of layers. More specifically, ${{\mathbf{G}}_{j,j}}(E)$ requires the surface GFs for the system with all layers to the left of layer j ($\widetilde{\mathbf{G}}_{j-1}^{\text{L}}(E)$) and for the system with all layers to the right of layer j ($\widetilde{\mathbf{G}}_{N-j}^{\text{R}}(E)$). We thus see that surface effects persist as long as these two surface GFs do not resemble the surface GFs for semi-infinite systems; the 'bulk' material would see semi-infinite systems on both sides. Because point (VII.9) describes how CBS determines the convergence of these surface GFs to their semi-infinite limits, the subsurface properties of a material are also tied to CBS.

Many of the CBS results for surfaces readily generalize to interfaces [27, 35, 158], where two materials come into contact with each other.

• (VII.12)  
  CBS describes the band 'lineup' between the materials on both sides of the interface when neither material is metallic [46, 81, 198, 199]. In almost every case, the two materials will have different workfunctions, meaning that electrons will spontaneously flow from one material to the other when the interface forms. The Fermi energies of the 'donor' and 'acceptor' materials will locally fall and rise near the interface, respectively, until there is no more net charge transfer across the interface. This results in an interface dipole. Tersoff [81] posited that charge is transferred from 'conduction-like' states of the donor to 'valence-like' states of the acceptor, implying that equilibration occurs when 'non-bonding' states from both materials are the next to be used. Using point (VII.3), the band lineup is thus determined by the branch point energies of the materials. This idea seems to work in many, but not all, cases [199].
• (VII.13)  
  Metal-induced gap states at metal-semiconductor interfaces arise from semiconductor states with complex k [36]. In the band gap, these states penetrate several layers into the metal, but decay according to $\text{Im}(k)$. Using similar reasoning to point (VII.12), the metal's Fermi energy is pinned at or near the branch point energy of the semiconductor. This interpretation has also been used to describe the level alignment between an oligomeric molecule and metal surfaces for molecular electronics applications [106, 107, 112, 116], although it may not always be appropriate [121].

If large layers cause ${{\mathbf{H}}_{\text{S}}}$ to be (near-)singular (see our example in figure 2), an obvious solution is to choose smaller layers. The problem is that smaller layers could lead to deviations from the block tridiagonal and block Toeplitz structure of the Hamiltonian, which were essential to deriving CBS in sections 4 and 5. We defer our discussion of going beyond the block tridiagonal structure to section 9 and focus here on relaxing the block Toeplitz structure.

In such a case, the Hamiltonian would instead have a periodic block Toeplitz structure [33, 35, 171, 201], where the blocks along each diagonal cycle through some finite number of distinct blocks. Such a system is essentially a superlattice. For example, there are two different small layers that alternate in figure 2; the Hamiltonian is structured as

$$\begin{eqnarray}\mathbf{H}=\left[\begin{array}{cc|cc|c} {{\mathbf{H}}_{\text{D},1}} & \mathbf{H}_{\text{S},1}^{\dagger} & \mathbf{0} & \mathbf{0} & \cdots \\ {{\mathbf{H}}_{\text{S},1}} & {{\mathbf{H}}_{\text{D},2}} & \mathbf{H}_{\text{S},2}^{\dagger} & \mathbf{0} & \cdots \\ \hline \mathbf{0} & {{\mathbf{H}}_{\text{S},2}} & {{\mathbf{H}}_{\text{D},1}} & \mathbf{H}_{\text{S},1}^{\dagger} & \cdots \\ \mathbf{0} & \mathbf{0} & {{\mathbf{H}}_{\text{S},1}} & {{\mathbf{H}}_{\text{D},2}} & \cdots \\ \hline \vdots & \vdots & \vdots & \vdots & \ddots \end{array}\right],\end{eqnarray}$$

where the lines show that the block Toeplitz structure is restored when two small layers are grouped together. We can address this structure by defining two transfer matrices (similarly for MMTs): One moves from a layer of type 1 to a layer of type 2 and vice versa. Using the associativity of matrix multiplication (MMT action), we can produce an 'aggregate' transfer matrix (MMT) for the periodic structure [35, 74, 136].

We can approximate the very evanescent states, which lead to the (near-)singularity, by states that are less evanescent [72, 89, 123, 135]. As long as the new states still decay (grow) much more rapidly than the physically-important states, the error can be negligible (or at least controllable). In essence, we compute the singular value decomposition [202] of ${{\mathbf{H}}_{\text{S}}}$ and artificially increase the small singular values. This results in an approximate coupling operator ${{\overline{\mathbf{H}}}_{\text{S}}}$ that has a condition number determined by the amount of increase. ${{\overline{\mathbf{H}}}_{\text{S}}}$ is thus invertible. The derivations in sections 4 and 5 are then readily applicable using ${{\overline{\mathbf{H}}}_{\text{S}}}$ instead of ${{\mathbf{H}}_{\text{S}}}$.

The invertibility of ${{\mathbf{H}}_{\text{S}}}$ was not needed until the final steps of deriving the transfer matrix (between equations (8) and (10)). Instead, we can rearrange equation (8) to obtain

$$\begin{eqnarray}&&\mathbf{H}_{\text{S}}^{\dagger}|{{\psi}_{n+1}}\rangle =\left(E\mathbf{I}-{{\mathbf{H}}_{\text{D}}}\right)|{{\psi}_{n}}\rangle -{{\mathbf{H}}_{\text{S}}}|{{\psi}_{n-1}}\rangle,\end{eqnarray}$$

where we stop just before inverting $\mathbf{H}_{\text{S}}^{\dagger}$. As in section 4, we combine this equation with the tautology $|{{\psi}_{n}}\rangle =\,|{{\psi}_{n}}\rangle$ to produce [85, 90, 113, 114, 120, 123, 153, 187]

$$\begin{eqnarray}&&\left[\begin{array}{c} \mathbf{H}_{\text{S}}^{\dagger}|{{\psi}_{n+1}}\rangle \\ |{{\psi}_{n}}\rangle \end{array}\right]=\left[\begin{array}{c} E\mathbf{I}-{{\mathbf{H}}_{\text{D}}} & -{{\mathbf{H}}_{\text{S}}} \\ \mathbf{I} & \mathbf{0} \end{array}\right]\left[\begin{array}{c} |{{\psi}_{n}}\rangle \\ |{{\psi}_{n-1}}\rangle \end{array}\right].\end{eqnarray}$$

Converting to the 'supercell' wavefunction,

$$\begin{eqnarray}\begin{array}{*{20}{l}}\left[\begin{array}{c} \mathbf{H}_{\text{S}}^{\dagger} & \mathbf{0} \\ \mathbf{0} & \mathbf{I} \end{array}\right]|{{\Psi}_{n+1}}\rangle =\left[\begin{array}{c} E\mathbf{I}-{{\mathbf{H}}_{\text{D}}} & -{{\mathbf{H}}_{\text{S}}} \\ \mathbf{I} & \mathbf{0} \end{array}\right]|{{\Psi}_{n}}\rangle,\end{array}\end{eqnarray}$$

and invoking Bloch's theorem (equation (1)), we get

$$\begin{eqnarray}{{\text{e}}^{\text{i}ka}}\left[\begin{array}{c} \mathbf{H}_{\text{S}}^{\dagger} & \mathbf{0} \\ \mathbf{0} & \mathbf{I} \end{array}\right]|{{\Psi}_{n}}\rangle =\left[\begin{array}{c} E\mathbf{I}-{{\mathbf{H}}_{\text{D}}} & -{{\mathbf{H}}_{\text{S}}} \\ \mathbf{I} & \mathbf{0} \end{array}\right]|{{\Psi}_{n}}\rangle.\end{eqnarray} \tag{ 33 }$$

We thus avoid inverting $\mathbf{H}_{\text{S}}^{\dagger}$ at the expense of having to solve a generalized eigenvalue problem [193] for CBS.

Generalized eigenvalue problems are more sophisticated than standard eigenvalue problems in that they admit both finite and infinite eigenvalues. For the present discussion, the generalized eigenvalue problem has a zero or infinite eigenvalue when

$$\begin{eqnarray}\left[\begin{array}{c} E\mathbf{I}-{{\mathbf{H}}_{\text{D}}} & -{{\mathbf{H}}_{\text{S}}} \\ \mathbf{I} & \mathbf{0} \end{array}\right]~\text{or}~\left[\begin{array}{c} \mathbf{H}_{\text{S}}^{\dagger} & \mathbf{0} \\ \mathbf{0} & \mathbf{I} \end{array}\right]\end{eqnarray}$$

is singular, respectively [193]. Both of these cases arise when ${{\mathbf{H}}_{\text{S}}}$ is singular, giving our justification that $\mathbf{T}$ tends to singularity as ${{\mathbf{H}}_{\text{S}}}$ becomes singular. Note that, when ${{\mathbf{H}}_{\text{S}}}$ is near-singular (such that the transfer matrix is analytically well-defined), $\mathbf{T}$ will have very large (in magnitude) eigenvalues that correspond to the nearly infinite eigenvalues of equation (33). In practice, these eigenvalues may lead to numerical instabilities when applying CBS [72, 83, 87], which partly motivates the next strategy.

Finally, we can build on the previous idea and remove the 0 and infinite eigenvalues. In essence, we deflate [200] the generalized eigenvalue problem so that only the states of interest remain [67, 89, 114, 123, 124, 126]. In other words, we project the highly evanescent states out of the problem. Unlike in section 8.2, these states are not approximated; rather, they are completely removed from consideration. After deflation, the resulting (and approximate) ${{\mathbf{H}}_{\text{S}}}$ is invertible such that the tools in sections 4 and 5 are again applicable. Compared to the generalized eigenvalue problem in section 8.3, this procedure requires a bit of extra effort to perform the deflation, but ultimately produces smaller matrices for the CBS computations.