Extended Bloch theorem for topological lattice models with open boundaries

TL;DR

This paper derives exact analytical solutions for boundary states in various topological lattice models with open boundaries, revealing the role of spectral mirror symmetry in obtaining full spectra including bulk and boundary states.

Contribution

It introduces a method based on spectral mirror symmetry to analytically solve for the entire spectrum of topological lattice models with open boundaries, including models previously lacking exact solutions.

Findings

Exact solutions for boundary states in nodal-line semimetals, spin-orbit coupled graphene, and topological insulators.

Identification of spectral mirror symmetry as key for full spectrum solutions.

Analytical spectrum derivation for the Lieb lattice with open boundaries in both directions.

Abstract

While the Bloch spectrum of translationally invariant noninteracting lattice models is trivially obtained by a Fourier transformation, diagonalizing the same problem in the presence of open boundary conditions is typically only possible numerically or in idealized limits. Here we present exact analytic solutions for the boundary states in a number of lattice models of current interest, including nodal-line semimetals on a hyperhoneycomb lattice, spin-orbit coupled graphene, and three-dimensional topological insulators on a diamond lattice, for which no previous exact finite-size solutions are available in the literature. Furthermore, we identify spectral mirror symmetry as the key criterium for analytically obtaining the entire (bulk and boundary) spectrum as well as the concomitant eigenstates, and exemplify this for Chern and $\mathcal Z_2$ insulators with open boundaries of co-dimension one. In the case of the two-dimensional Lieb lattice, we extend this further and show how to analytically obtain the entire spectrum in the presence of open boundaries in both directions, where it has a clear interpretation in terms of bulk, edge, and corner states.

Figures

Equations

$$\hat{H}$$

$$+\sum _{m,i}{c}_{A_i,m}^{†}{H}_{A\leftarrow B}^i{c}_{B,m}+H.c.$$

$$+\sum _{m,i}{c}_{A_i,m+1}^{†}{H}_{B\to A}^i{c}_{B,m}+H.c.,$$

$$|{\psi }_i⟩={\mathcal{N}}_i\sum _{m=1}^M{r}_i^m{c}_{A_i,m}^{†}|0⟩,$$

$$r_i=-\frac{{\left(H_{A\leftarrow B}^{†}\right)}^{1,i}}{{\left(H_{B\to A}^{†}\right)}^{1,i}}.$$

$$|{\psi }_i⟩={\mathcal{N}}_i\sum _{m=1}^M{r}_i^m\left(\sum _{\bar{n}=1}^n{c}_{A_{\bar{n}},m}^{†}{\varphi }_{i,\bar{n}}\right)|0⟩$$

$$r_i=-\frac{H_{A\leftarrow B}^{†}{\varphi }_i}{H_{B\to A}^{†}{\varphi }_i},$$

$$c_{A_i,k_{\perp }}^{†}$$

$$c_{B,k_{\perp }}^{†}$$

$$H(k_{\perp })$$

$${\Psi }_l(k_{\perp },m):=e^{i{k}_{\perp }m}{\Psi }_l(k_{\perp }).$$

$${\Psi }_{\mathrm{Bulk},l,i}(k_{\perp ,j},m)=\mathcal{A}{\Psi }_{l,i}(k_{\perp ,j},m)+\mathcal{B}{\Psi }_{l,i}(-k_{\perp ,j},m)$$

$$=\mathcal{A}{e}^{i{k}_{\perp ,j}m}{\Psi }_{l,i}(k_{\perp ,j})+\mathcal{B}{e}^{-i{k}_{\perp ,j}m}{\Psi }_{l,i}(-k_{\perp ,j}),$$

$${\Psi }_{\mathrm{Bulk},l,n+1}(k_{\perp ,j},m=0)={\Psi }_{\mathrm{Bulk},l,n+1}(k_{\perp ,j},m=M)$$

$$=\mathcal{A}{\Psi }_{l,n+1}(k_{\perp ,j})+\mathcal{B}{\Psi }_{l,n+1}(-k_{\perp ,j})=0$$

$$\Rightarrow \mathcal{B}/\mathcal{A}=-{\Psi }_{l,n+1}(k_{\perp ,j})/{\Psi }_{l,n+1}(-k_{\perp ,j}),$$

$${\Psi }_{\mathrm{Bulk},l,i}(k_{\perp ,j},m)=$$

$$e^{i{k}_{\perp ,j}m}{\Psi }_{l,i}(k_{\perp ,j})-e^{-i{k}_{\perp ,j}m}\frac{{\Psi }_{l;n+1}(k_{\perp ,j})}{{\Psi }_{l;n+1}(-k_{\perp ,j})}{\Psi }_{l,i}(-k_{\perp ,j}).$$

$$H_{B\to A}{\Psi }_B(m-1)+H_A{\Psi }_A(m)+H_{A\leftarrow B}$$

$$=E{\Psi }_A(m),$$

$$H_{B\to A}{r}^{m-1}{\Psi }_B+H_A{r}^m{\Psi }_A+H_{A\leftarrow B}{r}^m{\Psi }_B=E{r}^m{\Psi }_A.$$

$$H_A{\Psi }_A(1)+H_{A\leftarrow B}{\Psi }_B(1)=E{\Psi }_A(1)$$

$$\Rightarrow H_A{\Psi }_A+H_{A\leftarrow B}{\Psi }_B=E{\Psi }_A,$$

$$H_A(k_{||})$$

$$H_{A\leftarrow B}$$

$$H_{B\to A}$$

$$r_i(k_{||})=-\frac{{\left(H_{A\leftarrow B}^{†}\right)}^{i,i}}{{\left(H_{B\to A}^{†}\right)}^{i,i}}.$$

$$H_A=\left(\begin{array}{ccc}0& t& 0\\ t& 0& t(1+{\text{e}}^{i{k}_1})\\ 0& t(1+{\text{e}}^{-i{k}_1})& 0\end{array}\right),H_B=0,$$

$$H_{A\leftarrow B}=\left(\begin{array}{c}0\\ 0\\ t\end{array}\right),H_{B\to A}=\left(\begin{array}{c}2t{\text{e}}^{-\frac{i}{2}{k}_1}\text{cos}(k_2/2)\\ 0\\ 0\end{array}\right),$$

$$k_{\perp }=3k_3-k_2+k_1,$$

$$H=\left(\begin{array}{ccccc}{H}_A& H_{A\leftarrow B}& 0& 0& 0\\ H_{A\leftarrow B}^{†}& H_B& H_{B\to A}^{†}& 0& 0\\ 0& H_{B\to A}& H_A& \cdots & 0\\ 0& 0& ⋮& \ddots & H_{B\to A}^{†}\\ 0& 0& 0& H_{B\to A}& H_A\end{array}\right).$$

$$r_i({\vec{k}}_{||})=-\frac{{\varphi }_{i,3}}{2{\text{e}}^{\frac{i}{2}{k}_1}\text{cos}(k_2/2){\varphi }_{i,1}},$$

$$H_A=-H_B=t_2\sin (k_2-k_1),H_{B\to A}=t_1{\text{e}}^{\frac{i}{2}(k_2-k_1)},$$

$$H_{A\leftarrow B}=t_1(1+{\text{e}}^{i(k_2-k_1)}),k_{\perp }=-(k_1+k_2)/2,$$

$$H_A=-H_B=t_2\text{sin}(k_2-k_1){\tau }_z,H_{B\to A}=t_1{\text{e}}^{\frac{i}{2}(k_2-k_1)}{\tau }_0,$$

$$H_{A\leftarrow B}=t_1(1+e^{i(k_2-k_1)}){\tau }_0,k_{\perp }=-(k_1+k_2)/2,$$

$$H_A=-H_B=\lambda \left\{-\sqrt{2}\right[\text{sin}\left(k_1\right)+\text{sin}\left(k_2\right)]{\tau }_x$$

Full text

I Introduction

One of the main subjects of investigation in the condensed matter community are topological phases of matter, which already in the single-particle limit offer a plethora of fascinating properties such as the appearance of dissipationless, robust boundary states [1, 2, 3]. The boundary modes are not only a signature feature of first-order topological phases—e.g., Chern insulators, $\mathcal{Z}_{2}$ insulators, and strong, three-dimensional topological insulators feature chiral edge modes, helical edge modes, and two-dimensional Dirac states, respectively, on their boundaries [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]—higher-order topological phases are also marked by the appearance of boundary modes, which now localize to terminations with a codimension larger than one [17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32].

To access these boundary modes, one needs to consider a system with open boundary conditions, which greatly enlarges the unit cell and consequently renders the problem of finding analytical solutions more challenging and generally intractable in the limit of large system sizes. Despite this increased difficulty, we showed in previous works that exact solutions for the boundary-mode wave functions of any codimension $d$ can be straightforwardly obtained on a large family of $D$-dimensional lattices by making use of the naturally present destructive interference and without the need for fine-tuned tight-binding parameters [33, 28, 29]. The lattices we investigate consist of repeated structures that are subdivided in $A$ and $B$-type motifs, and have $A$ motifs at their boundaries, e.g., Fig. 1(c). On these lattices, we find $n$ exact solutions, where $n$ is the number of degrees of freedom of the $A$ motifs, which exponentially localize to the boundaries with a completely disappearing amplitude on the $B$-type motifs and whose eigenenergies correspond to the eigenenergies of the local Hamiltonian on $A$. Other methods also exist [34, 35, 36, 37, 38, 39, 40, 41], which are restricted to finding solutions for models with open boundary conditions in one direction only, while a generalization to higher-order phases is also made in an approximate fashion in Ref. 42. Despite this recent progress, there are several important questions regarding the exact treatment of the eigensystem of lattice models with open boundary conditions that remain unanswered, and we here systematically address them building on our framework introduced in Refs. 33, 28, 29.

Firstly, we investigate the premises for when it is possible to obtain, in addition to recovering boundary states, all remaining bulk states, and we find that this is the case when there is a mirror symmetry in the spectrum of the periodic Bloch Hamiltonian, $E(\vec{k}_{||},k_{\perp})=E(\vec{k}_{||},-k_{\perp})$, where the mirror relates the momenta $k_{\perp}$ in the direction of the open boundary conditions. In the presence of this symmetry, the problem of finding the complete eigensystem reduces to the problem of solving for the standing waves in an open string. Using this method, we are also able to show that the boundary wave functions that we have previously recovered are the only boundary modes in the system. We also provide a case study of the Lieb lattice where we obtain the full spectrum with open boundary conditions in both directions, where the solutions have a particularly illuminating interpretation in terms of bulk, edge, and corner states, respectively.

This leads to the second question, which addresses the issue of completeness of our boundary eigensystem, namely, in the absence of a spectral mirror symmetry, are the boundary-mode solutions we have found unique, or are there other boundary modes in the system that we have not found? This conundrum is resolved upon using a result derived using transfer matrices regarding the structure of edge states, which explicitly allows us to show that the boundary-mode solutions we find are the only solutions that exist.

Finally we note that certain mechanisms such as spin-orbit coupling (SOC) have have not been accounted for in previous studies. Here we investigate when it is possible to still solve for the boundary eigensystem, or the entire eigensystem in case of spectral mirror symmetry, in the presence of SOC. We find that as long as the spin on the $A$ motifs is well defined, and the coupling between different $A$ and $B$ motifs is invariant under spin rotation, we can still use our formalism to find analytical solutions. This is a particularly interesting problem to address in the case of three-dimensional strong topological insulators, where the spin sectors need to be explicitly coupled via SOC to obtain Dirac surface states.

We illustrate our findings by introducing explicit examples. To showcase our previously developed method, we introduce a “vanilla” example, which does not have a spectral mirror symmetry nor includes SOC, by considering a simple nearest-neighbor hopping model on the hyperhoneycomb lattice [43, 44, 45], which is the natural structure of $\beta$-Li2IrO3 [46]. This model realizes a nodal-line semimetal, which are systems that feature a degeneracy in the spectrum that forms a closed curve that is explicitly protected by the presence of spatial symmetries or time-reversal symmetry [47]. We find that our solution coincides with the wave function of the drumhead surface state, which corresponds to an exactly flat band in the spectrum. We thus solve for the drumhead state in a nodal-line semimetal in an exact fashion.

We continue by considering a variation to the Haldane model on the honeycomb lattice [5], which realizes a Chern-insulating phase and has the desired spectral mirror symmetry. We show that we are indeed not only able to find solutions to the chiral edge states but also solve for the complete eigensystem. Next, we introduce spin degree of freedom in this model by considering two time-reversed copies of the pseudo-Haldane model by which we realize a $\mathcal{Z}_{2}$ insulator in the form of a pseudo-Kane-Mele model [10, 11]. This expansion of our Hilbert space does not spoil the spectral mirror symmetry such that we can still solve for the full eigensystem with our boundary-mode solutions now corresponding to the helical edge states. For this model, the SOC does not connect the two spin sectors, and it is thus straight forward to solve this model.

However, we note that while the strategy of gluing two Chern insulators with opposite Chern number together indeed leads to the retrieval of a $\mathcal{Z}_{2}$ insulator, this approach does not suffice to obtain a strong topological insulator in three dimensions and we need to explicitly include a SOC, which connects the different spin sectors, to accomplish this phase. We implement a variation to the model introduced by Fu et al. [15], which realizes a strong topological insulator on the diamond lattice, and we show that we can find the boundary states both in the topologically trivial phase, the topologically weak phase, and the topologically strong phase. Moreover, we connect this phase to a nodal-line semimetal phase by which a direct connection between Dirac surface states and drumhead states is established.

In a final case study, we dissect the band spectrum of the two-dimensional Lieb lattice with open boundary conditions in two directions. By making use of the Bloch Hamiltonian as well as the two Hamiltonians relating to the two different types of edges that appear, we are able to solve for the entire spectrum analytically.

This paper is ordered as follows. We introduce our general method in Sec. II, where we summarize our formalism in Refs. [33, 28, 29] in Sec. II.1 for the sake of completion. In Sec. II.2, we show that the full eigensystem can be found in the presence of a spectral mirror symmetry, in Sec. II.3 we address the completeness of our boundary mode solutions, and we consider the case of SOC in Sec. II.4. The examples are introduced in Sec. III, and the discussion can be found in Sec. IV.

II Exact solutions

In this section, we study lattices with open boundary conditions in one direction of the form shown in Fig. 1(c), and show that in addition to finding exact solutions for the boundary modes [33, 28, 29], we can also fully diagonalize the entire system if the spectral symmetry $E(\vec{k}_{||},k_{\perp})=E(\vec{k}_{||},-k_{\perp})$ is present in the spectrum of the periodic Bloch Hamiltonian. Moreover, we show that this method still works in the presence of SOC by satisfying weak constraints.

We start by considering a $D$-dimensional lattice with $M$ unit cells, which consists of two distinguishable $(D-1)$-dimensional motifs, $A$ [red in Fig. 1(c)] and $B$ (blue). The $A$ motifs have $n$ degrees of freedom while the $B$ motifs have a single degree of freedom, and electrons are forbidden to hop between $A$ motifs directly. For simplicity, we set $D=1$ in the following discussion but note that the formalism can easily be generalized to higher dimensions by dimensional extension. We assume that hopping between $A$ and $B$ motifs takes place in a nearest-neighbor fashion, such that the real-space Hamiltonian reads

[TABLE]

where $c^{\dagger}_{A_{i},m}$ ($c^{\dagger}_{B,m}$) creates an electron in the $m$th unit cell within motif $A$ ($B$), and $i$, $i^{\prime}$ corresponds to the internal degrees of freedom on the $A$ motif, which may be a sublattice or orbital degree of freedom. The Hamiltonian is depicted in Fig. 1(c).

II.1 Boundary modes

In Refs. 33, 28, 29 we showed that if these lattices terminate with $A$ motifs on either side [cf. Fig. 1(c)], the following $n$ solutions exist:

[TABLE]

whose eigenvalues correspond to the eigenvalues $E_{i}$ of the Hamiltonian $H_{A}$ on motif $A$, and where $i$ labels the degrees of freedom on each $A$ motif, $\mathcal{N}_{i}$ is the normalization constants, and $c^{\dagger}_{A_{i},m}$ creates an electron on motif $A$ in unit cell $m$ with energy $E_{i}$. This wave function has zero weight on each of the $B$ motifs due to exact destructive interference, and we find

[TABLE]

Note that these simplified forms of Eqs. (2) and (3) assume that $H_{A}$ is written in its diagonal basis. When this is not the case, the exact wave functions take the following more general form:

[TABLE]

with

[TABLE]

where $\phi_{i,\bar{n}}$ is the $\bar{n}$th component of the $i$th eigenfunctions $\phi_{i}$ of $H_{A}$. Note that $r_{i}$ is a scalar since $H_{A\leftarrow B}^{\dagger}$ and $H_{B\rightarrow A}^{\dagger}$ are $1\times n$ matrices and $\phi_{i}$ is an $n$-component vector. The state $\ket{\psi_{i}}$ corresponds to a bulk mode when $|r_{i}|=1$, while localizing to the $A$ motif $m=1$ ($M$) when $|r_{i}|<1$ ($>1$). In the latter case, these solutions thus correspond to boundary states.

If the lattice instead terminates with an $A$ motif on one end and a $B$ motif on the other, we find that our solution in Eq. (2) is still relevant. Indeed, in this case, boundary modes appear on both ends of the lattice and when the lattice is sufficiently long, the solution in Eq. (2) may be mapped onto either end of the lattice as long as $|r_{i}|<1$. In this case, the solution on the boundary terminating with an $A$ ($B$) motif decays exponentially into the bulk and has an approximate zero weight on the $B$ ($A$) motifs. When $|r_{i}|>1$, on the other hand, the boundary states interfere with each other, such that the wave function in Eq. (2) is no longer a good approximation. Nevertheless, $|r_{i}|=1$ predicts when the boundary states enter the bulk and gap out up to finite size corrections and thus also provides key information for open systems without broken unit cells.

We emphasize that this method can also be used to find the edge and surface states of two- and three-dimensional crystals, respectively, as they can be seen as a collection of one-dimensional Hamiltonians parametrized by $\vec{k}_{||}$. In these cases, the function $r_{i}(\vec{k}_{||})$ is momentum dependent, which increases the likelihood of the states to localize to the edges or surfaces. Indeed, by setting the dimension $D=2$ and $3$ and implementing the relevant lattice Hamiltonian, we showed in Ref. 33 that these solutions correspond to the chiral, edge modes of a Chern insulator and the Fermi arcs in Weyl semimetals, respectively. To showcase this formalism explicitly, we show an explicit example in Sec. III.1 of a three-dimensional nodal-line semimetal, where the drumhead state is captured by our exact solutions.

We point out that in Refs. 28, 29, we generalized our formalism to higher-order boundaries by considering lattices that consist of a superposition of the chains like the one in Fig. 1(c) into an arrangement on which solutions of the form of $\ket{\psi_{i}}$ exist with additional prefactors $r^{\prime}_{i},r^{\prime\prime}_{i},\ldots$, which were found to capture corner modes as well as chiral hinge states. Indeed, in Sec. III.5 we show an explicit example of this generalization.

II.2 Complete diagonalization of mirror-symmetric lattice models

So far, we have only looked at retrieving boundary-mode solutions. In this section, we show we can find all energy eigenvalues and wave functions in the presence of the spectral symmetry $E(k_{\perp})=E(-k_{\perp})$ of the Bloch Hamiltonian. We start by considering a periodic chain with $2M$ unit cells, where each unit cell consists of $A$ and $B$ motifs as shown in Fig. 1(a). We find the Bloch Hamiltonian for this model by Fourier transforming the creation and annihilation operators in the real-space Hamiltonian in Eq. (1)

[TABLE]

such that $\hat{H}(k_{\perp})=(c^{\dagger}_{A_{i},k_{\perp}},c^{\dagger}_{B,k_{\perp}})H(k_{\perp})(c_{A_{i},k_{\perp}},c_{B,k_{\perp}})^{T}$ with

[TABLE]

We denote the eigenvectors of $H(k_{\perp})$ with $\Psi_{l}(k_{\perp})$, where $l=1,\ldots,n+1$ is the band index. The corresponding real-space wave function in the $m$th unit cell is then given by

[TABLE]

We proceed by cutting the chain in half by removing the $B$ sites in the [math]th (or $2M$th) and $M$th unit cell, such that two open chains remain with each $M$ A motifs and $(M-1)$ $B$ sites as shown in Fig. 1(b). Analogously to the standing waves in an open string, almost all wave functions for the geometry of such a chain can be obtained. If we now imagine that the $B$ sites at the end of the chain are reattached, and we suppose that the wave functions $\Psi_{l}(k_{\perp,j},m)$ and $\Psi_{l}(-k_{\perp,j},m)$ with $k_{\perp,j}=\pi j/M,\,j=1,\ldots,M-1$ have the same energy, i.e., $E_{l}(k_{\perp})=E_{l}(-k_{\perp})$, we can form a new eigenstate by considering the following linear combination

[TABLE]

where $i=1,\ldots,n+1$ refers to the component of the $l$th Bloch wave $\Psi_{l}(k_{\perp})$. Imposing that this linear combination, the so-called “generalized bulk state” $\Psi_{{\rm Bulk},l,i}(k_{\perp,j},m)$, vanishes on the “reattached” zeroth and $M$th $B$ motifs, we find, both for $m=0$ and $m=M$, that

[TABLE]

where $\Psi_{l,n+1}(k_{\perp,j})$ refers to the $(n+1)$th component of the $l$th Bloch wave $\Psi_{l}(k_{\perp})$, which coincides with the amplitude of $\Psi_{l}(k_{\perp})$ on the $B$ site. Setting $\mathcal{A}=1$ and subsequently $\mathcal{B}=-\Psi_{l,n+1}(k_{\perp,j})/\Psi_{l,n+1}(-k_{\perp,j})$, we thus find the following solution for the generalized bulk state:

[TABLE]

It follows that Eq. (9) must also be solutions for the open chain with $M$ $A$ motifs and $(M-1)$ $B$ sites. Moreover, for $1\leq j,j^{\prime}\leq M-1$, we find that these waves are orthogonal as long as $E_{l}(k_{\perp,j})\neq E_{l}(k_{\perp,j^{\prime}})$. We have thus found $\left(M-1\right)(n+1)$ standing wave solutions. The open chain, however, has $Mn+M-1=M(n+1)-1$ degrees of freedom, such that there are $n$ states unaccounted for. These missing states exactly correspond to the solutions in Eq. (2).

In the presence of the spectral symmetry, $E(k_{\perp})=E(-k_{\perp})$, the solutions in Eq. (2) are thus unique in the sense that no other boundary modes exist. If the spectral symmetry is broken, i.e., $E(k_{\perp})\neq E(-k_{\perp})$, the solutions in Eqs. (2) and (3) remain valid and the boundary modes persist, while the solutions in Eq. (9) no longer exist.

The construction presented here is applicable to a broad range of systems, and can be straightforwardly applied in multi-dimensional cases by requiring $E(\vec{k}_{||},k_{\perp})=E(\vec{k}_{||},-k_{\perp})$. In one-dimensional systems, this spectral symmetry is generated by mirror, inversion or time-reversal symmetry, whereas in higher-dimensional models only mirror symmetry leads to the desired spectral behavior. One might naively expect that other symmetries in two and three dimensions, e.g., inversion or time-reversal symmetry, would lead to a similar approach. This is prevented, however, by the fact that now the internal unit-cell labeling in the $A$ motifs starts to play a role and the unit cells in the $A$ motifs that are related to each other via these symmetries are not coupled to each other via the same site in $B$. The approach to find all bulk bands is, however, not destroyed by the presence of symmetries such as inversion or time reversal, as long as the desired spectral mirror symmetry is also present. Note, moreover, that the only requirement we impose is that the spectrum has to exhibit the mirror symmetry, while it is not necessary that the Hamiltonian itself is also mirror symmetric. Indeed, in Sec. III.2, we present an example closely related to the Haldane model, and find that the spectrum has the desired symmetry while the Hamiltonian does not.

II.3 Uniqueness

In the previous sections, we found $n$ boundary-mode solutions [cf. Eq. (2)], and showed they are unique in the case of a spectral mirror symmetry. However, even in the absence of such a spectral symmetry, it can be shown that no additional boundary states exist on the family of lattices that we study [cf. Fig. 1(c)]. Indeed, the results by Dwivedi and Chua derived for transfer matrices [36] imply that in our models all boundary states are exponentially decaying waves with only a single exponent, i.e., $\Psi(m)=r^{m}\Psi$, where $\Psi(m)$ denotes the wave function in the $m$th unit cell. For $1<m<M$, we thus find

[TABLE]

such that

[TABLE]

Similarly, for $m=1$, the following holds:

[TABLE]

where we note an equivalent expression is retrieved for $m=M$. Plugging this into Eq. (10), we immediately find $\Psi_{B}=0$. We thus show that all exponentially decaying wave functions that may exist on our lattices have a disappearing amplitude on all $B$ sites. The solutions in Eq. (2) are thus the only boundary mode solutions that exist.

II.4 Spin degree of freedom

We now end our general section on exact solutions by showing that the constraint of having a single degree of freedom on each $B$ motif can be relaxed. As an explicit example, we consider spin degree of freedom, and realize that in the most simple case we may restrict ourselves to models where the spin-up and down sectors are decoupled. These systems can be straightforwardly solved with our formalism, and indeed we consider such an example in the form of a two-dimensional $\mathcal{Z}_{2}$ insulator in Sec. III.3. While these $\mathcal{Z}_{2}$ insulators can be seen as being equal to the sum of two time-reversed Chern-insulator copies, their three-dimensional cousins cannot be interpreted as consisting of two separable systems, and spin-mixing terms generated by SOC are explicitly necessary to generate a three-dimensional topological insulator. These terms, however, may render the solutions to $r_{i}$ ambiguous.

There are, nevertheless, only two requirements that need to be satisfied for our formalism to be applicable also to this setting. Firstly, the eigenstates on the $A$ motifs need to have a well-defined spin, which can be established with a Hamiltonian of the form

[TABLE]

where $\vec{\tau}$ is the vector of Pauli matrices acting in spin space, and the spin-mixing terms are thus captured by the first term. Secondly, the coupling between the $A$ and $B$ motifs has to be invariant under spin rotations, such that generally

[TABLE]

This means that for a given $k_{||}$, we have two decoupled spins polarized in the $\pm\vec{d}(k_{||})$ directions, and we can solve the boundary modes of this system using our method, where Eq. (3) for finding $r_{i}(k_{||})$ is modified to

[TABLE]

The dependence of the spin-polarization on $k_{||}$ is particularly important in the case of the three-dimensional topological insulator, because it enables us to obtain exact solutions for the surface Dirac cones as we explicitly show in an example in Sec. III.4.

III Explicit examples

In this section, we present five explicit examples to showcase the formalism developed in the previous section. We start in Sec. III.1 by considering a nearest-neighbor model on the three-dimensional hyperhoneycomb lattice, which realizes a nodal-line semimetal and show that we can find the exact solution for the drumhead state as well as two Fermi-arc-like states. We continue by considering a variation to the Haldane model [5] in Sec. III.2, which we call the pseudo-Haldane model, and show that we can diagonalize the whole system by making use of the spectral symmetry $E(k_{x},k_{y})=E(k_{x},-k_{y})$ of the Bloch Hamiltonian to find the full spectrum. Next, we introduce spin-orbit coupling in Sec. III.3 by taking two time-reversed copies of the psuedo-Haldane model by which we realize a pseudo-Kane-Mele model [10, 11], which is a two-dimensional $\mathcal{Z}_{2}$ topological insulator. Again, we can solve the whole system exactly, and we find exact solutions to the helical edge states. In Sec. III.4, we allow for spin mixing terms inside the $A$ motif and show that we can solve for the surface Dirac modes by studying a three-dimensional topological insulator on the diamond lattice. Lastly, we treat the two-dimensional Lieb lattice with open boundary conditions in two directions in Sec. III.5, and show we can solve the full spectrum.

III.1 Nodal-line semimetal

We start by considering the three-dimensional hyperhoneycomb lattice shown in Fig. 2(a). This lattice has four sites in the unit cell, and we group them such that each $A$ [red in Fig. 2(a)] and $B$ (blue) motif contains three and one sublattice site, respectively, and we do not consider additional degrees of freedom such that $n=3$. The nearest-neighbor vectors $\vec{\delta}_{i}$ are explicitly indicated in Fig. 2(a), and read $\vec{\delta}_{1}=(\sqrt{3},\,0,\,1)/2$, $\vec{\delta}_{2}=(0,\,0,\,-1)$, $\vec{\delta}_{3}=(-\sqrt{3},\,0,\,1)/2$, $\vec{\delta}_{4}=(0,\,\sqrt{3},\,1)/2$, $\vec{\delta}_{5}=(0,\,-\sqrt{3},\,1)/2$, and the lattice vectors are defined as $\vec{a}_{1}=\vec{\delta}_{1}-\vec{\delta}_{3}$, $\vec{a}_{2}=\vec{\delta}_{4}-\vec{\delta}_{5}$, and $\vec{a}_{3}=\vec{\delta}_{3}-2\vec{\delta}_{2}+\vec{\delta}_{4}$. The Bloch Hamiltonian for this model is given in Eq. (7) with

[TABLE]

where $k_{i}\equiv\vec{k}\cdot\vec{a}_{i}$ and $t$ is the nearest-neighbor hopping parameter. Assuming open boundary conditions in $z$, and terminating with $A$ motifs on either end such that there is a bearded edge at $m=1$ and a zigzag edge at $m=M$, we find the Hamiltonian $\hat{H}=\Psi^{\dagger}H\Psi$ with $\Psi^{\dagger}=(\psi_{A,1}^{\dagger},\,\psi_{B,1}^{\dagger},\,\psi_{A,2}^{\dagger},\,\ldots,\,\psi_{B,M-1}^{\dagger},\,\psi_{A,M}^{\dagger})$, where $\alpha$ and $m$ in $\psi_{\alpha,m}^{\dagger}$ label the motif and unit cell, respectively, and

[TABLE]

The band spectrum for this model is shown in Fig. 2(b), where we have rescaled $k_{x}\rightarrow k_{x}/\sqrt{3}$. As $H_{A}$ is not written in a diagonal form, the exact solutions for the wave functions are given in Eq. (4) with $r_{i}(\vec{k}_{||})$ derived from Eq. (5)

[TABLE]

where $\phi_{i,\bar{n}}$ are the eigenfunctions of $H_{A}$. The eigenvalues $E_{i}$ associated with the wave functions in Eq. (4) are equivalent to the eigenvalues of $H_{A}$, i.e., $H_{A}\phi_{i}=E_{i}\phi_{i}$, with $E_{i}=0,\,\pm t\textrm{e}^{-ik_{1}}\sqrt{\textrm{e}^{ik_{1}}+3\textrm{e}^{2ik_{1}}+\textrm{e}^{3ik_{1}}}$, and are shown in red in Fig. 2(b). Plotting the inverse localization length $\xi^{-1}_{i}(\vec{k}_{||})=\textrm{ln}|r_{i}(\vec{k}_{||})|$, such that $\xi^{-1}_{i}(\vec{k}_{||})<0$ ($>0$) when $|r_{i}(\vec{k}_{||})|<1$ ($>1$), for the flat band as shown in the middle panel of Fig. 2(c), we clearly see that it corresponds to a drumhead state, and that the bulk-gap closing at which the band switches localization forms a closed loop (in black) revealing this is a nodal-line semimetal. The other two red bands also switch localization as shown in the upper and lower panel of Fig. 2(c). We have plotted energy contours (black, dashed lines) on top of the inverse localization length such that it is clearly visible that, following these contours, both these states switch surfaces. These bands thus have a Fermi-arc-like character in the sense that they are broken and live on opposite surfaces.

III.2 Chern insulator

We continue by introducing a variation to the Haldane model [5], which we call the pseudo-Haldane model, which has sublattice degree of freedom and realizes a Chern-insulating phase with chiral, one-dimensional edge modes on the two-dimensional honeycomb lattice shown in Fig. 3(a). In this case, the $A$ and $B$ motifs simply correspond to periodic chains consisting of the $A$ and $B$ sublattices of the honeycomb lattice only. Here we find the exact solution corresponding to the chiral edge modes, as was previously demonstrated for a lattice with chiral hinge states in Ref. 28, and also all solutions for the bulk states, which is a new result. The nearest-neighbor vectors $\vec{\delta}_{i}$ and lattice vectors $\vec{a}_{i}$ are explicitly indicated in Fig. 3(a), and read $\vec{\delta}_{1}=(-\sqrt{3},\,1,\,0)/2$, $\vec{\delta}_{2}=(0,\,-1,\,0)$, $\vec{\delta}_{3}=(\sqrt{3},\,1,\,0)/2$, $\vec{a}_{1}=\vec{\delta}_{3}-\vec{\delta}_{2}$ and $\vec{a}_{2}=\vec{\delta}_{1}-\vec{\delta}_{2}$. The Bloch Hamiltonian is given by Eq. (7) with

[TABLE]

where $k_{i}\equiv\vec{k}\cdot\vec{a}_{i}$, and $t_{1}$ ($t_{2}$) are the (next-)nearest-neighbor hopping parameters. We point out that this pseudo-Haldane model differs from the original Haldane model in that we have turned off all the next-nearest-neighbor hopping terms with a $k_{y}$-dependence, such that there is no long-range hopping connecting different $A$ motifs. The eigenvalues of the Bloch Hamiltonian have the spectral symmetry $E(k_{x},k_{y})=E(k_{x},-k_{y})$ such that we can use the method in section II.2 to fully diagonalize this model. Note that this model has additional spectral symmetries, $E(k_{x},k_{y})=E(-k_{x},k_{y})=E(-k_{x},-k_{y})$, which are in principle not required for the solvability of the model.

Assuming open boundary conditions in $y$ while keeping periodicity in $x$, and terminating with $A$ motifs on both sides, such that there is a zigzag and bearded edge at $m=1$ and $m=M$, respectively, we find the Hamiltonian given in Eq. (13) with the entries given in Eq. (14). The bulk bands computed using $E=E_{l}(k_{x},2\pi j/M)$ with $j=1,\ldots,M-1$ are shown in Fig. 4(a), while the chiral band with energy $H_{A}=t_{2}\,{\rm sin}(k)$ with $k=k_{1}-k_{2}=\sqrt{3}k_{x}$ is shown in Fig. 4(b). Combining these two spectra, we find Fig. 4(c), which corresponds exactly to results obtained via exact diagonalization (ED). Indeed, in Fig. 4(d) we compare the energies obtained via ED (in orange) and to those obtained analytically (in cyan) for a specific cut in $k$, and find that they overlap up to machine precision. The wave function corresponding to this chiral mode is given by Eq. (2) with $r(k)$ defined in Eq. (3) such that we find $r(k)=-2\,\textrm{cos}(k/2)$, and the mode necessarily attaches to the bulk when $k=\pm 2\pi/3\,\textrm{mod}\,2\pi$, i.e., when $|r(k)|=1$ [black dashed lines in Fig. 4(c)]. Moreover, from the specific form of $r(k)$ we deduce that the right (left) mover is localized to the $A$ lattice at the edge $m=M$ ($m=1$).

III.3 Two-dimensional topological insulator

We now increase the number of degrees of freedom on the $A$ and $B$ motifs by introducing spin, and study a modified version of the Kane-Mele model [10, 11], which we call the pseudo-Kane-Mele model. This system explicitly includes SOC terms, and realizes the quantum spin Hall effect on the graphene lattice, which features helical states on its edges. The Bloch Hamiltonian is given by Eq. (7) with

[TABLE]

with $\tau_{i}$ the Pauli matrices acting in spin space. Note that the terms proportional to $t_{2}$ now play the role of mirror-symmetric SOC [10, 11]. This model is closely related to the pseudo-Haldane model in Eq. (14) in the sense that the pseudo-Kane-Mele model realizes two time-reversed copies of the pseudo-Haldane model. There is only SOC inside each $A$ and $B$ motif, which leaves the spin up and down sectors decoupled. This means that, as in the previous section, we again have a spectral symmetry $E(k_{x},k_{y})=E(k_{x},-k_{y})$, and we can diagonalize the full model.

Taking open boundary conditions in $y$, the Hamiltonian is given in Eq. (13) with the entries in Eq. (15), and the accompanying band spectrum is shown in Fig. 5. The exact wave-function solutions are given in Eq. (2) with the eigenenergies $E_{1}(k)=-E_{2}(k)=t_{2}\,\textrm{sin}(k)$ (red bands in Fig. 5) the eigenvalues of $H_{A}$, where $k=k_{1}-k_{2}=\sqrt{3}k_{x}$. We use Eq. (11) to find $r_{i}(k)$ such that $r(k)=r_{1}(k)=r_{2}(k)=-2\,\textrm{cos}(k/2)$ as before. We have indicated the value of $|r(k)|$ explicitly in Fig. 5, and we see that two helical states localize to both edges.

III.4 Three-dimensional topological insulator

Here, we study a three-dimensional extension of the previous model in the form of a strong topological insulator in the diamond lattice [15]. We investigate the lattice with $(111)$ surfaces, such that the model consists of stacked, buckled honeycomb lattices as shown in Fig. 3(b)-3(d). Each of the triangular $A$ [red in Fig. 3(b)] and $B$ (in blue) sublattices of a buckled honeycomb sheet makes up the $A$ and $B$ motifs of the full model [cf. Figs. 3(c)-3(d)]. The nearest-neighbor vectors now read $\vec{\delta}_{1}=(-\sqrt{2}/\sqrt{3},\,\sqrt{2}/3,\,1/3)$, $\vec{\delta}_{2}=(0,\,-2\sqrt{2}/3,\,1/3)$, $\vec{\delta}_{3}=(\sqrt{2}/\sqrt{3},\,\sqrt{2}/3,\,1/3)$ and $\vec{\delta}_{4}=(0,\,0,\,-1)$, while the next-nearest-neighbor vectors $\vec{a}_{i}$ are defined as in Sec. III.2 with the additional vector $\vec{a}_{3}=\vec{\delta}_{2}-\vec{\delta}_{4}$. We adopt the Hamiltonian introduced by Fu et al. in Ref. 15, such that the Bloch Hamiltonian is given by Eq. (7) with

[TABLE]

with $k_{i}\equiv\vec{k}\cdot\vec{a}_{i}$, $t_{i}=t+\delta t_{i}$ the nearest-neighbor hopping parameter, $\lambda$ the next-nearest-neighbor SOC term, and $\tau_{i}$ acting in spin space as before. As opposed to the model in Ref. 15, we only include SOC inside each $A$ ($B$) layer and not between $A$ ($B$) layers. Taking open boundary conditions in $z$, the Hamiltonian is given in Eq. (13) with the entries in Eq. (16) and the energy spectrum is shown in Fig. 6. The exact wave-function solutions are given in Eq. (4), and the eigenvectors $\phi_{i}(\vec{k})$ and eigenvalues $E_{i}(\vec{k})$ (red bands in Fig. 6) are found by diagonalizing $H_{A}$. We thus find that there are always four (surface) Dirac cones in the system when $\lambda\neq 0$. As $H_{A\leftarrow B}$ and $H_{B\rightarrow A}$ are diagonal matrices, we find that, even though $H_{A}$ is not written in its diagonal form, we can use Eq. (11) to find $r_{i}(\vec{k})$

[TABLE]

Depending on the choice of $t_{i}$, $r(\vec{k})$ switch between $|r(\vec{k})|<1$ and $|r(\vec{k})|>1$ as a function of $\vec{k}$. This is indeed shown in Fig. 6, where we find that the localization of the Dirac cones changes depending on the choice for $t_{i}$. Indeed, when we turn on $\delta t_{1}$ only, we find that the system is a trivial insulator as the four cones are localized to the same surface [cf. Fig. 6(a)]. However, when we turn on $\delta t_{2}$ instead, we find we can realize both a weak [Fig. 6(b)] and strong topological insulator [Fig. 6(c)]. In the first case, this phase is found by tuning $\delta t_{2}$ such that the hopping $t_{2}=0$, which effectively reduces the diamond lattice into a stacking of decoupled honeycomb lattices with edges in the open surfaces. Each of these lattices realizes a quantum-spin Hall insulator with helical edge states, such that we find two Dirac cones on either of the surfaces. In Fig. 6(c), we tune $\delta t_{2}$ such that $t_{2}>t_{i},\,i=1,3,4$, and we find that one Dirac cone localizes to the bottom surface $m=1$ whereas the three other cones localize to the top surface $m=M$. We thus recover a strong topological insulator.

We note that even though in the case of the trivial insulator [cf. Fig. 6(a)] the surface state enters the bulk bands, it does not hybridize with the bands at the same energy. Indeed, the degeneracy in energy in this case is accidental and we could add local operators to the Hamiltonian to separate the solvable surface state from the continuous bulk spectrum. This can always be done whenever $|r(k)|\neq 1$, i.e., whenever Eq. (2) describes a surface state.

Lastly, we point out that if $\lambda=0$, this model [cf. Eq. (16)] reduces to a nodal-line semimetal with a doubly degenerate drumhead state. This is indeed already visible in Fig. 6(b), where the SOC is only turned on weakly. We have thus established a direct link between Dirac surface states and drumhead states. Moreover, we note that they are captured by the exact same wave function, and only their eigenenergies are affected by the change in $\lambda$.

III.5 Generalization to the Lieb lattice

Lastly, we turn to the two-dimensional Lieb lattice shown schematically in Fig. 7(a). This lattice is a two-dimensional generalization of the chain in Fig. 1(c), and consists of $A$ (in red), $B$ (in blue), and $B^{\prime}$ (in green) motifs with each one degree of freedom. The real-space Hamiltonian is depicted in Fig. 7(a), and reads

[TABLE]

In Refs. 28, 29, we showed that

[TABLE]

with $r=-t_{1}/t_{2}$ and $r^{\prime}=-t_{3}/t_{4}$, is an exact solution on this lattice with energy $E=m_{A}$, which corresponds to a corner state when $|r|\lessgtr 1$ and $|r^{\prime}|\lessgtr 1$, while it represents an edge state localized to the bottom (top) edge formed by $A$ and $B$ motifs when $|r|=1$ and $|r^{\prime}|<1$ ($>1$), and localized to the left (right) edge formed by $A$ and $B^{\prime}$ motifs when $|r^{\prime}|=1$ and $|r|<1$ ($>1$).

The Bloch Hamiltonian for this model reads $\hat{H}(\vec{k})=(c^{\dagger}_{A,\vec{k}},c^{\dagger}_{B,\vec{k}},c^{\dagger}_{B^{\prime},\vec{k}})H(\vec{k})(c_{A,\vec{k}},c_{B,\vec{k}},c_{B^{\prime},\vec{k}})^{T}$ with

[TABLE]

We find the eigenstates of this Bloch Hamiltonian have the spectral mirror symmetry $E(k_{x},k_{y})=E(-k_{x},k_{y})=E(k_{x},-k_{y})=E(-k_{x},-k_{y})$ regardless of our choice for parameters. We can also write the Hamiltonians corresponding to the edges formed by $A$ and $B$, and $A$ and $B^{\prime}$ motifs as

[TABLE]

respectively. In Fig. 7(b), we show the spectrum for this model with open boundary conditions in $x$ and $y$. In addition to the corner state solution (in red), whose eigenstate is given in Eq. (17) and whose energy is corresponds to $m_{A}$, we find we can also accurately compute the remainder of the band spectrum using $E(\pi j/M,\pi j^{\prime}/M^{\prime})$ with $E(k_{x},k_{y})$ the eigenstates of the Bloch Hamiltonian, to find the bulk bands (in blue), and $E_{\textrm{edge},AB}(\pi j/M)$ and $E_{\textrm{edge},AB^{\prime}}(\pi j^{\prime}/M^{\prime})$ to retrieve the edge bands (in orange and green, respectively), where $j=1,\ldots,M-1$ and $j^{\prime}=1,\ldots,M^{\prime}-1$. Indeed, in Fig. 7(c), we explicitly show the overlap of the results computed using ED (in orange) and those retrieved with the above description (in cyan) for a specific choice of parameters, and find an overlap up to machine precision.

IV Discussion

In this work, we address several open questions concerning the method we developed in Refs. 33, 28, 29 to find exact boundary state solutions. We show that in the presence of a spectral mirror symmetry, which relates the open boundaries, we can solve for the complete eigensystem, and that our boundary states are unique, also in the absence of this symmetry. Moreover, we consider the possibility of including SOC in our models, and show that only weak conditions need to be satisfied for us to be able to still solve the model. We emphasize that to our knowledge no exact solutions have been found that capture drumhead surface state, nor has the close relation between Dirac surface states and doubly-degenerate drumhead states—they have the exact same wave function—been established before. Moreover, to the best of our knowledge, we are the first to provide exact analytical solutions for the entire band spectrum of the two-dimensional Lieb lattice with open boundary conditions in two directions.

We note that other methods exist for finding exact wave-function solutions, such as the aforementioned transfer-matrix formalism [34, 35, 36], a rather involved method based on the generalization of Bloch’s theorem [37, 38], a construction specifically catered to time-reversal-symmetry preserving systems [39], and a formalism making use of boundary Green’s functions leading to effective boundary theories [40, 41, 42]. In case of the transfer-matrix approach, we point out that while transfer matrices specifically allow for the computation of bulk and boundary states, eigenvalues are treated on an equal footing with the other parameters in the model, and only solutions that are allowed to exist on the lattice come out as valid. This is in sharp contrast to the approach we have presented here, where the energy eigenvalues for the boundary states are computed directly from the Hamiltonian on the $A$ motif, $H_{A}$, while the bulk energies come about by solving the eigenvalue equation for the Bloch Hamiltonian $H(\vec{k}_{||},k_{\perp})$. We, moreover, note that the approach for finding bulk modes presented in Ref. 39 for systems with open boundary conditions in one direction carries some resemblance to our method in Sec. II.2. Whereas the focus in Ref. 39 is on time-reversal-symmetric models, we observe that the minimal generic requirement for using this method is in fact the presence of a spectral mirror symmetry relating the open boundaries.

Our work can presumably also be extended to other physical contexts where exact solutions are likely to be illuminating. One such interesting new area of application for our results in the field of non-Hermitian physics, relevant for a range of dissipative systems, in the context of topology, which has attracted much recent attention due to several exotic phenomena ranging from the breakdown of conventional bulk-boundary correspondence [49, 48] to the appearance of exceptional points, which correspond to a degeneracy both in the space of energies and the space of eigenvectors [50]. The failure of conventional bulk-boundary correspondence is strongly related to a phenomenon dubbed the “skin effect,” where bulk states pile up at the boundary [49, 50, 51, 52, 53, 54]. In this case, the continuum states thus behave as boundary modes, and a study of such bands following our method in Sec. II.2 is thus particularly interesting. Indeed, initial insight was recently enabled by insights from exact solutions in the simplest one and two-dimensional settings [49]. Our new way of obtaining the full spectrum analytically appears as particularly promising for understanding this rich class of systems better.

Acknowledgements.

We would like to thank M. Trescher for related collaborations, and V. Dwivedi for useful discussions. F.K.K. and E.J.B. are supported by the Swedish research council (VR) and the Wallenberg Academy Fellows program of the Knut and Alice Wallenberg Foundation. G.v.M. is supported by the research program of the Foundation for Fundamental Research on Matter (FOM), which is part of the Netherlands Organization for Scientific Research (NWO).