Topological Photonics: A Mathematical Perspective

1. Introduction

Topological photonics is a framework that follows both condensed matter physics and topology. It refers to designing the guiding properties of the propagating medium (e.g., a photonic crystal or a waveguide lattice) in such a way that the transport of electromagnetic energy is realized in unique, robust, and sometimes unexpected ways. Consider a simple thought experiment: imagine first the two-dimensional wave equation on a square domain, and assume homogeneous Dirichlet boundary conditions. We know that the accessible modes extend in periodic form throughout the whole domain and, in time, waves can propagate in all directions. This behavior is in response to the inherent symmetries of the medium. Imagine instead that we engineer the medium in such a way that all the energy concentrates in the boundary of the medium and propagates in only one direction. (In the language of optics, this would be seen as inhibiting back reflection and making the bulk medium act like an insulator).

Typically, in describing a photonic system, we refer to physical quantities such as frequency, wave vector, polarization, and dispersion. Instead, in the relatively new field of topological photonics, the term “topology” refers to a property of a photonic material that characterizes global behavior of the wavefunctions on their entire dispersion map; most importantly, this property takes quantized values. Think of this as characterizing the “genus” of an object, like a doughnut, with the “object” being described in wave vector space rather than in physical space. There are analogues in photonics to the topological fact that continuous deformations will not change the genus of an object. As an example, photonic topological insulators that are designed using artificial materials can support topologically nontrivial unidirectional states of light. These states are characterized by a particular “genus-like” number. Since this number is quantized, this unidirectional property will be robust to perturbations in the underlying photonic structure.

Photonics research often parallels or aims at explaining phenomena in other physical contexts. Bose-Einstein condensation in condensed matter physics is governed by the Gross-Pitaevskii equation, which is identical to the nonlinear Schrödinger equation that governs intense laser beam propagation in a dielectric medium such as air. In the quantum realm, nontrivial states of two-dimensional matter (e.g., a periodic lattice of atoms) with broken time-reversal symmetry can have the property that the bulk is an insulator while states (modes) exist that carry current along the sample edges without dissipation. The characteristic “genus-like” integer is called the Chern number (see section 3 for an example), which arises out of topological properties of the material’s band structure (see the discussion in section 2 and section 3 below).

In photonic crystals, a periodic variation of the dielectric properties of the medium affects photons in the same manner as solids modulate electrons (with the caveat that photons are bosons, while electrons are fermions). The question is whether the topological features are replicated in the analogous photonic system. In two foundational papers by Haldane and Raghu HR08RH08, the authors highlight the photonics analogue to quantum properties. They demonstrate the ingredients necessary to create a “one-way waveguide” which exhibits properties similar to the Quantum Hall Effect. While the model in HR08 has not been experimentally realized, it motivated further work by Wang, Chong, Joannopoulos, and Soljačić, in which they first predicted the existence of edge states in a magneto-optical crystal in the microwave regime WCJS08 and then demonstrated these experimentally WCJS09. Experiments by Rechtsman et al. RZP$^{+}$13 found topological edge states without the need for an external magnetic field by using a photonic crystal comprising helical waveguides.

Since then, the field of topological photonics has matured and continues to be very active, both in theory and experiments, as well as in the linear and nonlinear regimes. While we have briefly discussed its origins, it is not our purpose to give a detailed history of the field (for this purpose, we point the interested reader to the review articles JS14LCG$^{+}$22). Instead, we will focus our discussion on three prototypical examples: the one-dimensional Su-Schrieffer-Heeger model, the two-dimensional Haldane model, and the model of a photonic Floquet topological insulator from RZP$^{+}$13. We hope that this article highlights why this is a fertile area for mathematicians to explore and contribute to with their expertise.

2. SSH Model

The Su-Schrieffer-Heeger (SSH) model SSH79, is the simplest lattice model that exhibits topological features. It was devised to describe the electrical conductivity in a doped polyacetylene polymer chain. The lattice comprises repeating, two-node unit cells, where the couplings within and between unit cells are given by $t_1$ and $t_2$, respectively (Figure 1, top). The optical analogue is a linear lattice of waveguides in which the nearest-neighbor couplings are staggered (this can be implemented, e.g., by altering the physical spacings between the fibers). Mathematically, the SSH model can be described by the discrete nonlinear Schrödinger equation

$$\begin{equation} \begin{aligned} &i \dot{a}_n + t_1 b_n + t_2 b_{n-1} + \gamma |a_n|^2 a_n = 0 \\ &i \dot{b}_n + t_1 a_n + t_2 a_{n+1} + \gamma |b_n|^2 b_n = 0, \end{aligned} {\tag{1}} \end{equation}$$

where $(a_n, b_n)$ is the $n$th unit cell, and $\gamma >0$ is the strength of the cubic nonlinearity. (A rigorous mathematical derivation can be found in AC22).

Figure 1.

Top: schematic of SSH model, unit cell $(a_n, b_n)$ in dotted box. Coupling constant is $t_1$ within unit cell, and $t_2$ between unit cells. Bottom: reciprocal lattice, first BZ $[-\pi , \pi ]$ shown in red.

Our analysis follows that of AOP16, Chapter 1 and AC22, Section 8. The topological features of the optical SSH model can be understood by studying the linear model ($\gamma = 0$). As a first step, we look for plane wave solutions of the form

$$\begin{equation} {\mathbf{w}}(k) = \begin{pmatrix} a_n \\ b_n \end{pmatrix} = \begin{pmatrix} a \\ b \end{pmatrix}e^{i(k n - \lambda t)}, {\tag{2}} \end{equation}$$

where $\lambda$ is the frequency and $k$ is the wavenumber. Equation (2) is periodic in $k$ with period $2\pi$, since ${\mathbf{w}}(k+2m\pi ) = {\mathbf{w}}(k)$ for any integer $m$. The points $\{ 2 m \pi : m \in {\mathbb{Z}}\}$ define another linear lattice, which is called the reciprocal lattice (Figure 1, bottom). The first Brillouin zone (BZ) is the set of points closer to the origin than any other point of the reciprocal lattice, which in this case is the interval $[-\pi , \pi ]$. Due to the $2\pi$-periodicity, the BZ is topologically equivalent to the unit circle $S^1$.

Substituting the ansatz (2) into (1) and simplifying, we obtain the $k$-dependent eigenvalue problem $H(k) {\mathbf{v}}= -\lambda (k) {\mathbf{v}}$, where ${\mathbf{v}}= (a, b)^\top$, and $H(k)$ is the Hermitian matrix

$$\begin{equation*} H(k) = \begin{pmatrix} 0 & h(k) \\ h^*(k) & 0 \end{pmatrix}, \qquad h(k) = t_1 + t_2 e^{-ik}. \end{equation*}$$

Since $H(k)$ is $2\pi$-periodic, we only need to consider $k \in [-\pi , \pi ]$, i.e., in the first BZ. We note that since we are posing the problem on the full integer lattice, $k$ can take any value in $[-\pi , \pi ]$. The eigenvalues of $H(k)$ are

$$\begin{equation} \lambda (k) = \pm |h(k)| = \sqrt { t_1^2 + t_2^2 + 2 t_1 t_2 \cos k}, {\tag{3}} \end{equation}$$

which is the dispersion relation $\lambda (k)^2 = |h(k)|^2$ relating the frequency $\lambda$ and the wavenumber $k$. Each eigenvalue $\lambda (k)$ is a continuous function of the wavenumber on the first BZ, and is called a band. Since $H(k)$ is a $2 \times 2$ matrix, the SSH model has two bands. All of the bands of the system form its band structure, which is illustrated in the left column of Figure 2. (These terms are borrowed from solid state physics, where the band structure describes the energy levels that electrons can occupy in a solid).

Figure 2.

Left: band structure of the SSH model for $0 < t_2 < t_1$ (top) and $0 < t_1 < t_2$ (bottom). Right: circle in the complex plane traced counterclockwise by $h(k)$ for $k \in [-\pi , \pi ]$.

When $t_1 \neq t_2$, there is a space between the upper and lower bands, which is known as a band gap. This band gap has size $2 \Delta$, where $\Delta = |t_1 - t_2|$. The band gap closes when $t_1 = t_2$. Since the eigenvalues (3) are unchanged if $t_1$ and $t_2$ are exchanged, it appears at first glance that the cases $t_1 > t_2$ and $t_1 < t_2$ are identical, i.e., that the problem is symmetric about $t_1 = t_2$. Interestingly, this is not the case. For a complete picture, we need to look at the eigenvectors of $H(k)$ as well.

The normalized eigenvectors corresponding to the eigenvalues $\pm \lambda (k)$ are given by

$$\begin{equation} {\mathbf{v}}_\pm (k) = \frac{1}{\sqrt {2}}\begin{pmatrix} \pm \hat{h}(k) \\ 1 \end{pmatrix}, \quad \hat{h}(k) = \frac{h(k)}{|h(k)|}. {\tag{4}} \end{equation}$$

Since $\hat{h}(k)$ is a complex number of unit modulus whose argument is the same as that of $h(k)$, we can write ${\mathbf{v}}_\pm (k)$ as

$$\begin{equation} \begin{aligned} {\mathbf{v}}_\pm (k) &= \frac{1}{\sqrt {2}}\begin{pmatrix} \pm e^{i\phi (k)} \\ 1 \end{pmatrix} \\ \tan \phi &= \frac{\operatorname {\text{Im}}h(k)}{\operatorname {\text{Re}}h(k)} = -\frac{t_2 \sin k}{t_1 + t_2 \cos k}. \end{aligned} {\tag{5}} \end{equation}$$

As the wavenumber $k$ varies from $-\pi$ to $\pi$ over the BZ, the complex number $h(k)$ traces a clockwise circle in the complex plane with center $(t_1,0)$ and radius $t_2$. This circle encloses the origin when $t_1 < t_2$, but does not when $t_1 > t_2$ (Figure 2, right column). The topological invariant is the winding number of $h(k)$, which is the number of times $h(k)$ travels counterclockwise around the origin. We can see from Figure 2 that

$$\begin{equation} \operatorname {Ind}_h(0) = \begin{cases} -1 & t_1 < t_2 \\ 0 & t_1 > t_2 \end{cases}, {\tag{6}} \end{equation}$$

where a winding number of $-1$ represents a single clockwise trip around the origin. (The winding number is undefined if $t_1 = t_2$).

The same topological information can be obtained in a different way by computing a quantity known as the Berry phase Ber84 (also known as the Zak phase Zak89 in 1D). Intuitively, the Berry phase is the phase angle accumulated by a complex vector, e.g., one of the eigenvectors ${\mathbf{v}}_\pm (k)$, around a closed contour in $k$-space. We will take a brief digression to discuss these concepts, following AOP16, Chapter 2 and Van18, Chapter 3, before computing them for the SSH model.

Let ${\mathbf{v}}(k)$ be a normalized eigenvector of $H(k)$. For any wavenumbers $k_1$ and $k_2$, we define the relative phase between ${\mathbf{v}}(k_1)$ and ${\mathbf{v}}(k_2)$ by

$$\begin{equation} \gamma _{12} = \arg \langle {\mathbf{v}}(k_2), {\mathbf{v}}(k_1) \rangle , {\tag{7}} \end{equation}$$

where $\arg z$ is the phase, or argument, of the complex number $z$. (We are using the Hermitian inner product $\langle {\mathbf{u}}, {\mathbf{v}}\rangle = {\mathbf{u}}^\dagger {\mathbf{v}}= \sum _j u_j^* v_j$, where the dagger symbol denotes the conjugate transpose; the complex conjugation is placed on the first component to be consistent with the Dirac notation of quantum mechanics). The relative phase $\gamma _{12}$ satisfies the equation

$$\begin{equation} e^{i \gamma _{12}} = \frac{\langle {\mathbf{v}}(k_2), {\mathbf{v}}(k_1) \rangle }{|\langle {\mathbf{v}}(k_2), {\mathbf{v}}(k_1) \rangle |}. {\tag{8}} \end{equation}$$

It is important to note that the eigenvector ${\mathbf{v}}(k)$ is not unique. In particular, since it is a unit vector, it is specified only up to multiplication by a constant unit complex number $e^{i \theta }$. The transformation ${\mathbf{v}}(k) \mapsto e^{i\theta } {\mathbf{v}}(k)$ is called a gauge transformation. The relative phase $\gamma _{12}$ is not invariant under a gauge transformation, since if we take ${\mathbf{v}}(k_j) \mapsto e^{i\theta _j} {\mathbf{v}}(k_j)$, $\langle {\mathbf{v}}(k_2), {\mathbf{v}}(k_1) \rangle$ transforms to

$$\begin{equation*} \langle e^{i \theta _2} {\mathbf{v}}(k_2), e^{i \theta _1} {\mathbf{v}}(k_1) \rangle = e^{i(\theta _1-\theta _2)} \langle {\mathbf{v}}(k_2), {\mathbf{v}}(k_1) \rangle , \end{equation*}$$

thus $\gamma _{12} \mapsto \gamma _{12} + (\theta _1 - \theta _2)$.

We wish to define the change of phase of ${\mathbf{v}}(k)$ in such a way as to be gauge invariant. To do this, we take a sequence $(k_1, k_2, \dots , k_N)$ of $N$ points in $k$-space ordered in a loop. We then define the discrete Berry phase by

$$\begin{equation*} \begin{aligned} \gamma &= \arg e^{i(\gamma _{12} + \gamma _{23} + \dots + \gamma _{N,1})} \\ & = \arg \left( \langle {\mathbf{v}}(k_2), {\mathbf{v}}(k_1) \rangle \langle {\mathbf{v}}(k_3), {\mathbf{v}}(k_2) \rangle \cdots \langle {\mathbf{v}}(k_1), {\mathbf{v}}(k_N) \rangle \right), \end{aligned} \end{equation*}$$

which is the phase accumulated by ${\mathbf{v}}(k)$ around the loop. Unlike the relative phases $\gamma _{jk}$, the Berry phase is gauge invariant; if we take the gauge transformations ${\mathbf{v}}(k_j) \mapsto e^{\theta _j} {\mathbf{v}}(k_j)$, the Berry phase transforms to $\gamma + (\theta _1 - \theta _2) + (\theta _2 - \theta _3) + \dots + (\theta _N - \theta _1)$, which is equal to $\gamma$, since all of the $\theta _j$ cancel. We note that the Berry phase is only unique up to an integer multiple of $2 \pi$ unless we take the principal value of the argument, i.e., restrict $\arg z$ to $(-\pi , \pi ]$.

We now move from discrete to continuous. In particular, we wish to compute the phase accumulated by ${\mathbf{v}}(k)$ along a continuous, closed path. For small $\Delta k$, let $\Delta \gamma$ be the relative phase accumulated between ${\mathbf{v}}(k)$ and ${\mathbf{v}}(k+\Delta k)$. Following (8), $\Delta \gamma$ satisfies the equation

$$\begin{equation} e^{i \Delta \gamma } = \frac{\langle {\mathbf{v}}(k+\Delta k) , {\mathbf{v}}(k) \rangle }{|\langle {\mathbf{v}}(k+\Delta k), {\mathbf{v}}(k) \rangle |}. {\tag{9}} \end{equation}$$

Since $\Delta k$ is small and ${\mathbf{v}}$ is a unit vector, the denominator in (9) is approximately 1, thus

$$\begin{equation*} e^{i \Delta \gamma } \approx \langle {\mathbf{v}}(k+\Delta k),{\mathbf{v}}(k) \rangle . \end{equation*}$$

Expanding both sides in a Taylor series to first order in $\Delta \gamma$ and $\Delta k$ and simplifying, we find that $\Delta \gamma$ is approximately given by

$$\begin{equation} \Delta \gamma \approx i \left\langle {\mathbf{v}}(k), \frac{d}{dk} {\mathbf{v}}(k) \right\rangle \Delta k. {\tag{10}} \end{equation}$$

We define the Berry connection by

$$\begin{equation} A(k) = i \left\langle {\mathbf{v}}(k), \frac{d}{dk} {\mathbf{v}}(k) \right\rangle = i {\mathbf{v}}(k)^\dagger \frac{d}{dk} {\mathbf{v}}(k), {\tag{11}} \end{equation}$$

which is the coefficient of $\Delta k$ on the RHS of (10). We then define the Berry phase to be the integral of Berry connection around a closed contour $\mathcal{C}$:

$$\begin{equation} \gamma = \oint _\mathcal{C} A(k) dk. {\tag{12}} \end{equation}$$

As in the discrete case, the Berry connection is not gauge invariant, while the Berry phase is invariant under gauge transformations (modulo integer multiples of $2 \pi$).

Returning to the SSH model, we first compute the Berry connection using the eigenvector ${\mathbf{v}}_+(k)$:

$$\begin{align*} A(k) &= i \left\langle \frac{1}{\sqrt {2}}\begin{pmatrix} e^{i \phi (k)} \\ 1 \end{pmatrix}, \frac{1}{\sqrt {2}}\begin{pmatrix} i \phi '(k) e^{i \phi (k)} \\ 0 \end{pmatrix} \right\rangle \\ &= -\frac{1}{2} \phi '(k). \end{align*}$$

We then compute the Berry phase by integrating the Berry connection from $-\pi$ to $\pi$. This is a closed contour in the BZ since the endpoints of $[-\pi , \pi ]$ correspond to the same point on the unit circle $S^1$. The Berry phase is

$$\begin{align*} \gamma &= -\frac{1}{2} \int _{-\pi }^\pi \phi '(k) dk = \begin{cases} \pi & t_1 < t_2 \\ 0 & t_1 > t_2, \end{cases} \end{align*}$$

which is the change in phase of the eigenvector ${\mathbf{v}}_+(k)$ over the BZ. Since it is a constant multiple of (6), it conveys the same information as the winding number.

Figure 3.

Eigenvalues of $B$ for $0 < t_2 < t_1$ (top left), $t_1 = t_2$ (top right), and $0 < t_1 < t_2$ (bottom left). Edge mode eigenvectors ${\mathbf{v}}= ({\mathbf{a}}, {\mathbf{b}})$ for $0 < t_1 < t_2$ (bottom right); ${\mathbf{a}}$ (solid blue line) and ${\mathbf{b}}$ (dotted orange line) are either in-phase (top) or out-of-phase (bottom). $N=25$ unit cells with Dirichlet boundary conditions.

The fundamental topological difference between the two asymmetric lattice configurations ($t_1<t_2$ and $t_1>t_2$) becomes evident when we consider a finite lattice. Specifically, we take a lattice comprising $2N$ waveguides ($N$ unit cells) with Dirichlet boundary conditions at the two ends ($b_0 = 0$ and $a_{N+1} = 0$). For the linear system, solutions are standing waves of the form $({\mathbf{a}}, {\mathbf{b}})^\top e^{-i\lambda t}$, where ${\mathbf{a}}= (a_1, \dots , a_N)$ and ${\mathbf{b}}= (b_1, \dots , b_N)$ represent the $a$ and $b$ sublattices of the system. Substituting this ansatz into (1) yields the eigenvalue problem $B{\mathbf{v}}= -\lambda {\mathbf{v}}$, where ${\mathbf{v}}= ({\mathbf{a}}, {\mathbf{b}})^\top \in \mathbb{C}^{2N}$ and $B$ is the off-diagonal block matrix

$$\begin{equation*} B = \begin{pmatrix} 0 & B_0 \\ B_0^\top & 0 \end{pmatrix}, \quad B_0 = \begin{pmatrix} t_1 \\ t_2 & t_1 \\ & \ddots & \ddots \end{pmatrix}. \end{equation*}$$

An intuitive understanding of the difference between the two cases can be gained by considering the two extreme configurations, where one of the coupling constants is set to 0. If $t_1 > 0$ and $t_2 = 0$, the lattice comprises $N$ independent dimers with internal coupling constant $t_1$. The eigenvalues of $B$ are $\pm t_1$, each with multiplicity $N$. On the other hand, if $t_1 = 0$ and $t_2 > 0$, the lattice instead comprises $N-1$ independent dimers (staggered from the ones in the previous case) with internal coupling constant $t_2$, as well as two unconnected nodes at the ends of the lattice. In addition to eigenvalues at $\pm t_2$, each with multiplicity $N-1$, the matrix $B$ has two eigenvalues at $0$. The eigenvectors corresponding to these zero eigenvalues are $(1,0,\dots ,0)^\top$ and $(0, \dots ,0,1)^\top$. These are called edge modes, since they are localized at the ends of the lattice.

Since $B$ is a $2N\times 2N$ matrix, its spectrum is a discrete set of $2N$ eigenvalues, as opposed to the two continuous bands of eigenvalues found from the dispersion relation in the infinite lattice case. The eigenvalues of $B$ are shown in Figure 3 for $0 < t_2 < t_1$, $t_1 = t_2$, and $0 < t_1 < t_2$. The two asymmetric configurations contain a “gap,” which closes when $t_1=t_2$. This eigenvalue gap is analogous to the band gap in the infinite lattice case. When $0 < t_2 < t_1$, there are no eigenvalues in this gap, and all of the eigenvectors are nonlocalized. When $0 < t_1 < t_2$, however, there are two eigenvalues close to (but not exactly at) 0 which lie within this gap (these eigenvalues approach 0 in the limit $N \rightarrow \infty$). As in the case where $t_1=0$, these eigenvalues correspond to edge modes, since ${\mathbf{a}}$ and ${\mathbf{b}}$ are localized to the left and right edges of the lattice, respectively (Figure 3, bottom right). All of the remaining modes are nonlocalized.

Finally, we briefly comment on what occurs when a cubic nonlinearity ($\gamma > 0)$ is present (see MS21 for a more thorough treatment). Standing wave solutions of the form $({\mathbf{a}}, {\mathbf{b}})^\top e^{-i\lambda t}$ solve the equation $B {\mathbf{v}}+ \gamma {\mathbf{v}}^3 + \lambda {\mathbf{v}}= 0$, where ${\mathbf{v}}= (a_1, \dots , a_N, b_1, \dots , b_N)$. Numerical continuation experiments show that the edge modes from the linear model persist for small $\gamma$.

3. Haldane Model

We now turn to a two-dimensional model. We start with a honeycomb lattice (Figure 4), which is constructed from a two-site unit cell, with sites labeled $a$ and $b$. These unit cells tile the plane periodically along the two primitive lattice vectors ${\mathbf{v}}_1 = \left( \frac{3}{2}, \frac{\sqrt {3}}{2}\right)$ and ${\mathbf{v}}_2 = \left( \frac{3}{2}, -\frac{\sqrt {3}}{2}\right)$ to obtain a hexagonal lattice. We note that the $a$-sites and $b$-sites form two offset, triangular sublattices. This structure is similar to that of the material graphene, which is a hexagonal lattice constructed entirely from carbon atoms. The spatial location of a unit cell is specified by the vector ${\mathbf{r}}_{\mathbf{n}}= m {\mathbf{v}}_1 + n {\mathbf{v}}_2$, where ${\mathbf{n}}= (m,n) \in {\mathbb{Z}}^2$. It is therefore natural to index the unit cells by the vector ${\mathbf{n}}$; the locations of the lattice sites $a_{\mathbf{n}}$ and $b_{\mathbf{n}}$ in unit cell ${\mathbf{n}}$ are $r_{\mathbf{n}}$ and $r_{\mathbf{n}}+ (1,0)$, respectively. Each node in the honeycomb lattice is connected to its three nearest neighbors with coupling strength $t_1>0$, which is a coupling between sublattices. The directions of the nearest-neighbor (NN) couplings are given by the vectors $\boldsymbol{\delta }_1 = \left( \frac{1}{2}, -\frac{\sqrt {3}}{2}\right)$, $\boldsymbol{\delta }_2 = \left( \frac{1}{2}, \frac{\sqrt {3}}{2}\right)$, and $\boldsymbol{\delta }_3 = -(\boldsymbol{\delta }_1+\boldsymbol{\delta }_2) = \left( -1, 0 \right)$, which are depicted in Figure 4. The resulting linear model can be written as

$$\begin{equation} \begin{aligned} & i \dot{a}_{\mathbf{n}}+ t_1 \sum _{\langle {\mathbf{m}}\rangle } b_{{\mathbf{m}}} = 0, \quad & i \dot{b}_{\mathbf{n}}+ t_1 \sum _{\langle {\mathbf{m}}\rangle } a_{{\mathbf{m}}} = 0, \end{aligned} {\tag{13}} \end{equation}$$

where the angle brackets indicate that the sum is taken over nearest neighbors. The system (13) obeys time-reversal symmetry, i.e., is invariant under the transformation $t \mapsto -t$, $(a_{\mathbf{n}}, b_{\mathbf{n}}) \mapsto (a_{\mathbf{n}}^*, b_{\mathbf{n}}^*)$.

The Haldane model Hal88 adds two more terms to the honeycomb model. We will see that this results in a band gap in the spectrum, similar to what occurs in the SSH model with unequal couplings. First, the Haldane model has an on-site energy term of magnitude $t_0$, which takes opposite signs on the $a$ and $b$ sublattices. In addition, there is an imaginary, next-nearest neighbor (NNN) coupling term with strength $i t_2$. (In the original Haldane model, this coupling term has the complex strength $t_2 e^{i\phi }$; we take $\phi = \pi /2$ here for simplicity). Each node is coupled to its six next-nearest neighbors, which is a coupling within sublattices. These couplings are staggered so that there is no net flux into or out of a single lattice site. The directions of the next-nearest neighbor couplings are given by the vectors $\boldsymbol{\delta }'_1 = {\mathbf{v}}_1$, $\boldsymbol{\delta }'_2 = -{\mathbf{v}}_2$, and $\boldsymbol{\delta }'_3 = {\mathbf{v}}_2 - {\mathbf{v}}_1$ (Figure 4). Using this notation, the linear model for the Haldane lattice is

$$\begin{equation} \begin{aligned} & i \dot{a}_{\mathbf{n}}+ t_0 a_{\mathbf{n}}+ t_1 \sum _{\langle {\mathbf{m}}\rangle } b_{{\mathbf{m}}} \pm i t_2 \sum _{\lAngle {\mathbf{m}}\rAngle } a_{{\mathbf{m}}} = 0 \\ & i \dot{b}_{\mathbf{n}}- t_0 b_{\mathbf{n}}+ t_1 \sum _{\langle {\mathbf{m}}\rangle } a_{{\mathbf{m}}} \pm i t_2 \sum _{\lAngle {\mathbf{m}}\rAngle } b_{{\mathbf{m}}} = 0, \end{aligned} {\tag{14}} \end{equation}$$

where the double angle brackets indicate that the sum is taken over next-nearest neighbors. The signs of the NNN couplings are indicated by the arrows in Figure 4, where outward and inward pointing arrows denote couplings of $+i t_2$ and $-i t_2$, respectively. The arrangement of the arrows in two staggered, counterclockwise triangles ensures no net flux results from the NNN term. Time-reversal symmetry is broken when $t_2 \neq 0$, but is unaffected by the on-site term $t_0$.

Figure 4.

Schematic of the Haldane lattice. Rhombus is unit cell with sites $a$ and $b$. Primitive lattice vectors ${\mathbf{v}}_1$ and ${\mathbf{v}}_2$. Nearest neighbor coupling vectors $\boldsymbol{\delta }_1$, $\boldsymbol{\delta }_2$, and $\boldsymbol{\delta }_3$. Next-nearest neighbor coupling vectors $\boldsymbol{\delta }'_1$, $\boldsymbol{\delta }'_2$, and $\boldsymbol{\delta }'_3$.

As in the SSH model, the first step is to compute the band structure, which is found by looking for plane wave solutions to (14) of the form

$$\begin{equation} {\mathbf{w}}({\mathbf{k}}) = \begin{pmatrix} a_{\mathbf{n}}\\ b_{\mathbf{n}}\end{pmatrix} = \begin{pmatrix} a \\ b \end{pmatrix}e^{i({\mathbf{k}}\cdot {\mathbf{r}}_{\mathbf{n}}- \lambda t)}, {\tag{15}} \end{equation}$$

where