# Unified approach to Floquet lattices, topological insulators, and their
nonlinear dynamics

Mark Ablowitz Department of Applied Mathematics, University of Colorado,
Boulder, CO 80301, USA Justin T. Cole Department of Mathematics, University of
Colorado, Colorado Springs, CO 80918, USA S. D. Nixon22footnotemark: 2

###### Abstract

A unified method to analyze the dynamics and topological structure associated
with a class of Floquet topological insulators is presented. The method is
applied to a system that describes the propagation of electromagnetic waves
through the bulk of a two-dimensional lattice that is helically-driven in the
direction of propagation. Tight-binding approximations are employed to derive
reduced dynamical systems. Further asymptotic approximations, valid in the
high-frequency driving regime, yield a time-averaged system which governs the
leading order behavior of the wave. From this follows an analytic calculation
of the Berry connection, curvature and Chern number by analyzing the local
behavior of the eigenfunctions near the critical points of the spectrum.
Examples include honeycomb, Lieb and kagome lattices. In the nonlinear regime
novel equations governing slowly varying wave envelopes are derived. For the
honeycomb lattice, numerical simulations show that for relatively small
nonlinear effects a striking spiral patterns occurs; as nonlinearity increases
localized structures emerge and for somewhat higher nonlinearity the waves
collapse.

##  1 Introduction

The theoretical study of topological insulators began nearly 40 years ago with
the Quantum Hall Effect in electronics [1, 2]. Topological insulators are
materials that behave as an insulator in its interior, but conduct edge waves
along their boundaries. An effective model predicting the same phenomena for
periodic systems without any magnetic flux termed the Quantum Anomalous Hall
Effect), was introduced by Haldane in [3]. The underlying structure was later
reinterpreted to fit a photonic setting where light waves in an optical
lattice with suitable permittivity took the place of electrons in magnetic
field [4].

The quantization of the Hall conductance can be expressed in terms of
invariant integrals [5] and gives a unique classification for the geometric
structure of the problem [6]. At the same time, by interpreting the eigenspace
of a Hermetian operator as a line bundle over its parameter space [7], it
follows that these topological invariants are related to the Berry phase [8].
In an optical setting, this means the Floquet-Bloch spectral bands for a
periodic waveguide array, form a natural line bundle over the Brillouin zone.
The relevant characterization is the Chern class (and the related Chern
number) [9].

The broken time-reversal symmetry that characterizes Chern insulators and
protected edge states has since been realized experimentally using numerous
techniques. In [10], the first observation of an electromagnetic protected
edge state, capable of propagating around defects, was observed in a magneto-
optical setting. In this set up, the symmetry breaking was achieved by the
application of an external magnetic field to a square lattice of ferrite rods.

An experiment in the photonic regime that does not rely on magnetic fields was
reported in [11]. Here a honeycomb lattice is constructed with waveguides that
are helically varying in the direction of beam propagation, which acts as the
symmetry breaking mechanism and supports topologically protected edge states.
Under an appropriate change of variables, this has the same effect as a
magnetic field in the governing Schrödinger equation. A tight-binding system
in linear and nonlinear regimes was derived in [12]. The field of Topological
Photonics has since seen many important developments with novel experimental
designs that realize spectral bands with underlying nontrivial topological
invariants (see review articles [13, 14]).

Topological invariants can only change discretely; as a result, induced
behavior will be robust against perturbations/defects. The phase transition
points of the band structure, where the Chern number change, occur at
degeneracies of the linear spectrum. In a 2D optical lattice, the simplest
point degeneracies occur at conically shaped Dirac points in the dispersion
bands. Thus, the modulation of honeycomb, Lieb, or kagome lattices, all of
which exhibit Dirac/Dirac-like points in their band structure, can lead to
nontrivial topologies. The edges states for honeycomb [12, 15] and Lieb [16]
lattices with longitudinal driving have been investigated analytically using
tight-binding approximations in the lattice nonlinear Schrödinger equation
(NLS). Here, nonlinear edge solitons are also found to be unidirectional
solutions. A tight binding model associated with magnetized ferrite rods was
recently derived using Wannier functions in [17].

This paper presents a unified approach that starts with the lattice NLS
equations and derives averaged, constant coefficient, models suitable for the
analytical study of underlying Floquet-Bloch eigenfunctions, computation of
the Chern number and modeling nonlinear dynamics of the original equation. We
use a tight-binding approximation that forgoes the typical Peierls
transformation in favor of a more analytically tractable discrete system and
further reduce the model by averaging the effect of modulations in the
direction of propagation in the bulk (not gap) regime. We apply this method to
a wide range of examples including honeycomb, Lieb and kagome lattices in an
effort to both illustrate it’s ease of use and highlight the commonalities
among topological spectra. For honeycomb lattices, it is noteworthy that this
averaging technique connects with the well known Haldane model [3] in
Brillouin/Fourier space.

In the nonlinear regime, we derive envelope equations for the evolution of
wave packets near the degeneracy of the linear spectrum. There are two common
spectral phenomena we consider. The first is the typical Dirac point which
appears in the spectrum of the honeycomb and kagome lattices. The second is a
Dirac point with a flat band intersecting it as occurs in the Lieb and 1/5
depleted lattices and results in a novel third-order system. In both cases, we
investigate the relative effects of driving and nonlinear refraction.

The nonlinear evolution in the case of a honeycomb lattice is particularly
interesting. The equations are a Dirac system with an additional Dirac mass
term which was not present in prior conical diffraction studies [18, 19]. With
this extra term numerical computations show that with fixed mass at relatively
small nonlinearity there are striking spiral patterns that gradually disperse;
as nonlinearity increases a localized stationary soliton-like structure
appears and at yet higher nonlinearity collapse occurs. This similar to what
occurs with two dimensional NLS equations cf. [20].

###  1.1 Governing Wave Equation

The governing equation for intense paraxial light propagating in an optical
lattice is the nonlinear Schrödinger equation, which we write in normalized
form as

$\mathrm{i}\frac{\partial\Psi}{\partial
z}+\nabla^{2}\Psi-V(\mathrm{\mathbf{r}},z)\Psi+\sigma\left|\Psi\right|^{2}\Psi=0$
(1)

where $\Psi$ is the envelope of the electric field, $z$ is the propagation
distance, $\mathrm{\mathbf{r}}$ is the transverse plane, and the potential
$V(\mathrm{\mathbf{r}},z)$ is the index of refraction for a doubly periodic
waveguide array in the transverse direction subject to longitudinal driving in
the $z$ direction

$V(\mathrm{\mathbf{r}},z)={V}(\mathrm{\mathbf{r}}-\mathrm{\mathbf{h}}(z)).$
(2)

Optical lattices of this type have been constructed experimentally [11, 21]
and used in the construction of topological insulators. The broken time-
symmetry created by rotating the lattice induces an effective magnetic field
and yields nontrivial topology in the spectrum and the construction of
topologically protected edge modes.

The base lattice is doubly periodic in the transverse plane with periods
$\mathrm{\mathbf{v}}_{1}$ and $\mathrm{\mathbf{v}}_{2}$ and consists of $L$
sublattices. A unit cell of the lattice can be pictured as a parallelogram
with sides $\mathrm{\mathbf{v}}_{1}$ and $\mathrm{\mathbf{v}}_{2}$ and $L$
lattice sites in its interior, one from each sublattice; see Figures 1, 4(a),
6 and 7. The full lattice comes from tiling the plane with copies of the unit
cell, i.e.,

$V(\mathrm{\mathbf{r}})=V_{0}^{2}~{}-\sum_{(n,m)\in\mathbb{Z}^{2}}\sum_{\ell=1}^{L}\widetilde{V}(\mathrm{\mathbf{r}}-\mathrm{\mathbf{r}}_{nm}^{\ell})$
(3)

where
$\mathrm{\mathbf{r}}_{nm}^{\ell}=n~{}\mathrm{\mathbf{v}}_{1}+m~{}\mathrm{\mathbf{v}}_{2}+\mathrm{\mathbf{d}}_{\ell}$,
the vectors $\mathrm{\mathbf{d}}_{\ell}$ define the locations of the
sublattice sites within the unit cell and $(n,m)$ is a numbering of the unit
cell tiling. For simplicity, we assume that the individual lattice sites have
a uniform Gaussian shape,

$\widetilde{V}(\mathrm{\mathbf{r}})=V_{0}^{2}\mathrm{e}^{-\mathrm{\mathbf{r}}^{2}/w}$
(4)

and the lattice sites occur at the minimums of ${V}$. We take driving function
of the form

$\mathrm{\mathbf{h}}(z)=\frac{\kappa}{\omega}\Big{(}\cos(\omega
z),~{}\sin(\omega z)\Big{)}.$ (5)

It is now convenient to move to a frame of reference that follows the lattice,
$\widetilde{\mathrm{\mathbf{r}}}=\mathrm{\mathbf{r}}-\mathrm{\mathbf{h}}(z)$,
which transforms equation (1) into

$\mathrm{i}\frac{\partial\Psi}{\partial
z}-\mathrm{i}\mathrm{\mathbf{h}}^{\prime}(z)\cdot\widetilde{\nabla}\Psi+\widetilde{\nabla}^{2}\Psi-{V}\big{(}\widetilde{\mathrm{\mathbf{r}}}\big{)}\Psi+\sigma\left|\Psi\right|^{2}\Psi=0.$
(6)

In contrast to some previous discrete models for Floquet lattices, cf. [11,
15, 16] we do not employ the common Peierls transform. As a useful
consequence, the coefficients of the final averaged system can be expressed in
a closed form. For future work, this allows for the study of lattices where
the waveguides are not rotating in unison since the
$\mathrm{\mathbf{h}}^{\prime}$ term can be treated as the local rotation for
the purposes of tight-binding approximations. Results for this non-Peierls
approximation agree with those found in models that do take the Peierls
transform on lattices that have a moving frame of reference. In Section 3, we
consider a the fast oscillation, $\omega\gg 1$, regime. Alternate methods
instead consider an adiabatic Floquet problem [22].

##  2 Tight-binding approximation

The tight-binding approximation assumes a strong lattice regime where the
contrast between the waveguides and the background is large, i.e., $V_{0}\gg
1$ [19, 12, 15]. With the lattice sites well-separated, we replace equation
(1) by a suitable discrete model. Local behavior near each lattice site is
approximated by an exponentially localized orbital function that approximates
the discrete eigenmode for a single site lattice times a modulating function
that adjusts the local interaction magnitude. Due to the rapid decay of these
orbitals, only the interactions between nearest neighbor lattice sites are
considered in the final approximation.

In the neighborhood around $\mathrm{\mathbf{r}}_{nm}^{\ell}$, we define
$\phi_{nm}^{\ell}(\mathrm{\mathbf{r}})=\phi_{0}(\mathrm{\mathbf{r}}-\mathrm{\mathbf{r}}_{nm}^{\ell})$
where $\phi_{0}(\bm{r})$ is the orbital for a lattice site centered at the
origin satisfying the equation

$\left[\mu+\nabla^{2}+\widetilde{V}(\mathrm{\mathbf{r}})\right]\phi_{0}(\mathrm{\mathbf{r}})=0.$
(7)

where the orbital function has been normalized to have $L^{2}(\mathbb{R}^{2})$
norm one. The orbital expansion for the wave field, $\psi$, is now given by

$\psi(\mathrm{\mathbf{r}},z)\approx\sum_{(n,m)\in\mathbb{Z}^{2}}\left(\sum_{\ell=1}^{L}a_{nm}^{\ell}(z)\phi_{nm}^{\ell}(\mathrm{\mathbf{r}})\right)\mathrm{e}^{-\mathrm{i}\left(\mu+V_{0}^{2}\right)z}$
(8)

where each orbital function $\phi_{nm}^{\ell}(\mathrm{\mathbf{r}})$ is
modulated by a related coefficient $a_{nm}^{\ell}(z)$. A system of difference
equations is obtained by substituting anatz (8) into the modified Schrödinger
equation (1) and taking the inner product of the result with
$\phi_{nm}^{\ell}$ for each of lattice sites.

By combining the cluster of lattice sites within each unit cell,
$\mathrm{\mathbf{a}}_{nm}(z)=(a_{nm}^{1}(z),\ldots,a_{nm}^{\ell}(z))^{T}$, we
arrive at a nearest neighbor system of difference equations. Note: No
interaction occurs more than one unit cell away in any direction.

$\mathrm{i}\frac{\mathrm{d}\mathrm{\mathbf{a}}_{nm}}{\mathrm{d}z}+\sum_{i=-1}^{1}\sum_{j=-1}^{1}\mathcal{L}_{i}^{j}(z)\mathrm{\mathbf{a}}_{n+i,m+j}(z)+\sigma\rm{N}(\mathrm{\mathbf{a}}_{nm})\,\mathrm{\mathbf{a}}_{nm}=0$
(9)

where the $\mathcal{L}_{i}^{j}$ are a set of connectivity matrices,
$\rm{N}(\mathrm{\mathbf{a}})=\rm{diag}\left(|a^{1}|^{2},\ldots,|a^{L}|^{2}\right)$
is the onsite nonlinear response and $\widetilde{\sigma}$ is a constant. The
connectivity matrices are given by the definition

$\mathcal{L}_{i}^{j}(z)=\sum_{\ell=1}^{L}\sum_{\mathrm{\mathbf{w}}_{\ell
k}}\delta_{ij}(\mathrm{\mathbf{w}}_{\ell
k})\big{[}q(z,\mathrm{\mathbf{w}}_{\ell k})\big{]}_{\ell k}$ (10)

where the inner sum is taken over all the vectors $\mathrm{\mathbf{w}}$ that
point from a lattice site on the sublattice $\ell$ to a nearest neighbor on
any other sublattice $k$, the delta function
$\delta_{ij}(\mathrm{\mathbf{w}})$ is $1$ $\spadesuit$if the vector
$\mathrm{\mathbf{w}}$ with starting point $\mathrm{\mathbf{r}}_{00}^{\ell}$
has its end point in the lattice cell $(i,j)$ and $0$ otherwise$\spadesuit$,
and $\big{[}q(z,\mathrm{\mathbf{w}})\big{]}_{\ell k}$ is the matrix with the
function $q(z,\mathrm{\mathbf{w}})$ in row $\ell$ column $k$ and $0$
everywhere else. For example, the matrix $\mathcal{L}_{1}^{0}(z)$ will have
one element, of the form $q(z,\mathrm{\mathbf{w}})$, for each lattice site
that has a nearest neighbor in the lattice cell one over in the direction
$\mathrm{\mathbf{v}}_{1}$, whereas $\mathcal{L}_{0}^{0}(z)$ will have one
element for each nearest within the same lattice cell. The connectivity
matrices follow the symmetry properties,

$\displaystyle\mathcal{L}_{i}^{j}(z)$
$\displaystyle=\Big{[}\mathcal{L}_{-i}^{~{}j}(z)\Big{]}^{H}$ (11a)
$\displaystyle\mathcal{L}_{i}^{j}(z)$
$\displaystyle=\Big{[}\mathcal{L}_{~{}i}^{-j}(z)\Big{]}^{H}$ (11b)

where H is the conjugate transpose. See Appendix A for further details.

The topological properties of the lattice, such as the Chern Number, come from
the associated complex vector bundle constructed by taking the reciprocal unit
cell, or Brillouin zone, as a base space with attached vector space at each
point defined by the eigenspace for a particular band, see Section 4. To carry
this out, we take a discrete Fourier transform of the system (9), i.e, we
assume plane waves of the form

$\widetilde{\mathrm{\mathbf{a}}}_{nm}(z)=\bm{\alpha}(\mathrm{\mathbf{k}},z)\mathrm{e}^{\mathrm{i}\mathrm{\mathbf{k}}\cdot(n\mathrm{\mathbf{v}}_{1}+m\mathrm{\mathbf{v}}_{2})},$
(12)

where
$\mathrm{\mathbf{k}}=k_{x}\hat{\mathrm{\mathbf{i}}}+k_{y}\hat{\mathrm{\mathbf{j}}}$
and move the problem to the wavenumber domain. The end result of the tight-
binding approximations is a system of nonlinear differential equations,
parameterized by $\mathrm{\mathbf{k}}$, of the form

$\mathrm{i}\displaystyle\frac{\partial\bm{\alpha}}{\partial
z}+M\left(z,\mathrm{\mathbf{k}}\right)\bm{\alpha}+\widetilde{\sigma}\rm{N}(\bm{\alpha})\bm{\alpha}=0.$
(13)

where the matrix $M(z,\mathrm{\mathbf{k}})$ is $L\times L$, periodic in $z$,
doubly periodic in $\mathrm{\mathbf{k}}$, Hermitian, due to the symmetry of
the pairwise interactions, and has zeros along the diagonal, since any onsite
interactions can be removed by a phase transformation. The individual elements
of $M(z,\mathrm{\mathbf{k}})$ can be written as

$M_{ij}(z,\mathrm{\mathbf{k}})=e^{-\mathrm{i}\mathrm{\mathbf{k}}\cdot(\mathrm{\mathbf{d}}_{j}-\mathrm{\mathbf{d}}_{i})}\sum_{\mathrm{\mathbf{w}}}q(z,\mathrm{\mathbf{w}})e^{\mathrm{i}\mathrm{\mathbf{k}}\cdot\mathrm{\mathbf{w}}}$
(14)

where

$q(z,\mathrm{\mathbf{w}})=q_{0}-\mathrm{i}s_{0}\mathrm{\mathbf{h}}^{\prime}(z)\cdot\mathrm{\mathbf{w}},$
(15)

is a $z$-dependent coefficient function, with $q_{0}$ and $s_{0}$ constants,
derived in the Appendix A and the sum is taken over the set of vectors,
$\mathrm{\mathbf{w}}$, that point from a site on sublattice $i$ to a site on
sublattice $j$ (see 1, 4(a), 6 or 7).

##  3 Averaged System

Since equation (13) is a system of differential equations with periodic
coefficients, it usually does not have an explicit solution. So, we consider
the case where the frequency of the driving $\mathrm{\mathbf{A}}(z)$ occur on
a faster scale than the background. This leads to an averaged system,
independent of $z$. Employing a multiple scales approximation where the
lattice matrix is a function of the fast scale:
$\zeta=\displaystyle\frac{z}{\varepsilon}$, $0<\varepsilon\ll 1$,

$M\left(z\right)=M\left(\displaystyle\zeta\right),$ (16)

and the eigenvectors are expanded as
$\bm{\alpha}=\bm{\alpha}_{0}+\varepsilon\bm{\alpha}_{1}+\varepsilon^{2}\bm{\alpha}_{2}\ldots\,$.
In terms of the driven potential in Section 1.1, we define
$\omega=\varepsilon^{-1}$. Note: The dependence on $\mathrm{\mathbf{k}}$ is
not explicitly represented for this analysis, since the steps proceed
independently.

Matrix lattice equation (13) becomes

$\varepsilon^{-1}\mathrm{i}\frac{\partial\bm{\alpha}}{\partial\zeta}+\mathrm{i}\frac{\partial\bm{\alpha}}{\partial
z}+M\left(\zeta\right){\bm{\alpha}}=0.$ (17)

At $\mathrm{O}(\varepsilon^{-1})$,

$\mathrm{i}\frac{\partial\bm{\alpha}_{0}}{\partial\zeta}=0$ (18)

and as a result $\bm{\alpha}_{0}$ will be strictly a function of $z$, i.e.,
$\bm{\alpha}_{0}(\zeta,z)=\bm{\alpha}_{0}(z)$. At $\mathrm{O}(1)$,

$\mathrm{i}\frac{\partial\bm{\alpha}_{1}}{\partial\zeta}+\mathrm{i}\frac{\partial\bm{\alpha}_{0}}{\partial
z}+M(\zeta)\bm{\alpha}_{0}=0.$ (19)

To avoid secular growth in $\bm{\alpha}_{1}$, we require that the average
effect of the forcing from $\bm{\alpha}_{0}$ is at higher order. Define
$\left<f\right>=\displaystyle\frac{1}{T}\displaystyle\int_{0}^{T}f(\zeta)\,\mathrm{d}\zeta$
where $T$ is the period of $M(\zeta)$. We note that the terms:
$\mathrm{i}\partial_{z}\bm{\alpha}_{0}+\left<M\right>\bm{\alpha}_{0}$ grow
without bound in $\zeta$. Now, the secularity condition can be written as

$\mathrm{i}\frac{\partial\bm{\alpha}_{0}}{\partial
z}+\left<M\right>\bm{\alpha}_{0}=\varepsilon\mathrm{\mathbf{F}}_{0}.$ (20)

Substituting
$\mathrm{i}\partial_{z}\bm{\alpha}_{0}=-\left<M\right>\bm{\alpha}_{0}$ into
equation (19) allows us to solve for first correction as

$\bm{\alpha}_{1}=\mathrm{i}M_{1}(\zeta)\bm{\alpha}_{0}(z)$ (21)

where
$M_{1}(\zeta)=\displaystyle\int_{0}^{\zeta}\left(M\left(\widetilde{\zeta}\right)-\left<M\right>\right)\mathrm{d}\tilde{\zeta}$.

At $\mathrm{O}(\varepsilon)$,

$\mathrm{i}\frac{\partial\bm{\alpha}_{2}}{\partial\zeta}+\mathrm{i}\frac{\partial\bm{\alpha}_{1}}{\partial
z}+M(\zeta)\bm{\alpha}_{1}+\mathrm{\mathbf{F}}_{0}=0$ (22)

where the terms left over from the previous order must now be balanced as
well. We can use (21) to replace $\bm{\alpha}_{1}$ and once again substitute
$\mathrm{i}\partial_{z}\bm{\alpha}_{0}=-\left\langle
M\right\rangle\bm{\alpha}_{0}$ to simplify the correction terms we will be
deriving for the $\bm{\alpha}_{0}$ equation. These substitutions transform
equation (22) into

$\displaystyle\mathrm{i}\frac{\partial\bm{\alpha}_{2}}{\partial\zeta}$
$\displaystyle+\mathrm{i}\big{(}M(\zeta)M_{1}(\zeta)-M_{1}(\zeta)\left<M\right>\big{)}\bm{\alpha}_{0}$
$\displaystyle+\varepsilon^{-1}\left(\mathrm{i}\frac{\partial\bm{\alpha}_{0}}{\partial
z}+\left<M\right>\bm{\alpha}_{0}\right)=0$ (23)

The secularity condition now gives us a final (for our purposes) higher order
equation for $\bm{\alpha}_{0}$ that will be used to approximate the spectral
properties of the lattice,

$\mathrm{i}\frac{\partial\bm{\alpha}_{0}}{\partial z}+\left\langle
M\right\rangle\bm{\alpha}_{0}+\varepsilon\left\langle
N\right\rangle\bm{\alpha}_{0}=0$ (24)

where $\left\langle
N\right\rangle=\mathrm{i}\left[\left<MM_{1}\right>-\left<M_{1}\right>\left<M\right>\right]$.
The average equation (24) has constant coefficients and hence allows for
significant analytic insights; e.g. the computation of the Chern number.

For $M(z,\mathrm{\mathbf{k}})$ of the form derived in the tight-binding
analysis from Section 2 and determined by the nearest neighbor interactions in
the lattice, the matrix $\left\langle N\right\rangle$ will be determined by
all the next-nearest neighbors. This means, the averaged equation at first
order is

$\left\langle
M\right\rangle_{ij}=e^{-\mathrm{i}\mathrm{\mathbf{k}}\cdot(\mathrm{\mathbf{d}}_{j}-\mathrm{\mathbf{d}}_{i})}\sum_{\mathrm{\mathbf{w}}}q_{0}e^{\mathrm{i}\mathrm{\mathbf{k}}\cdot\mathrm{\mathbf{w}}}$
(25)

where $\left\langle q(z,\mathrm{\mathbf{k}})\right\rangle=q_{0}$ and the
elements of the first correct matrix, $\left\langle N\right\rangle$, take the
form

$\varepsilon\left\langle
N\right\rangle_{ij}=e^{-\mathrm{i}\mathrm{\mathbf{k}}\cdot(\mathrm{\mathbf{d}}_{j}-\mathrm{\mathbf{d}}_{i})}\displaystyle\sum_{\mathrm{\mathbf{w}}_{1},\mathrm{\mathbf{w}}_{2}}Q(\mathrm{\mathbf{w}}_{1},\mathrm{\mathbf{w}}_{2})\mathrm{e}^{\mathrm{i}\mathrm{\mathbf{k}}\cdot(\mathrm{\mathbf{w}}_{1}+\mathrm{\mathbf{w}}_{2})}$
(26)

where the vectors $\mathrm{\mathbf{w}}_{1}$ and $\mathrm{\mathbf{w}}_{2}$ are
such that $\mathrm{\mathbf{w}}_{1}$ points from an element of sublattice $i$
to any nearest neighbor and $\mathrm{\mathbf{w}}_{2}$ that starts at the tip
of $\mathrm{\mathbf{w}}_{1}$ and ends at an element of the sublattice $j$. So,
the averaged system (24) depends on next-nearest neighbors, even though the
original Floquet system (13) only depends on nearest neighbor interactions.
The coefficients for the next-nearest neighbor terms are defined by the
function

$\displaystyle
Q(\mathrm{\mathbf{w}}_{1},\mathrm{\mathbf{w}}_{2})=\mathrm{i}\Bigg{(}$
$\displaystyle\left\langle
q(z,\mathrm{\mathbf{w}}_{1})\int_{0}^{z}q(\tilde{z},\mathrm{\mathbf{w}}_{2})-\left\langle
q(\tilde{z},\mathrm{\mathbf{w}}_{2})\right\rangle\mathrm{d}\tilde{z}\right\rangle$
$\displaystyle-\left\langle\int_{0}^{z}q(\tilde{z},\mathrm{\mathbf{w}}_{1})-\left\langle
q(\tilde{z},\mathrm{\mathbf{w}}_{1})\right\rangle\mathrm{d}\tilde{z}\right\rangle\left\langle
q(z,\mathrm{\mathbf{w}}_{2})\right\rangle\Bigg{)}$

which can be found explicitly for the driving function (5) as

$Q(\mathrm{\mathbf{w}}_{1},\mathrm{\mathbf{w}}_{2})=\Big{(}\frac{\mathrm{i}}{2}s_{0}^{2}\kappa^{2}(\mathrm{\mathbf{w}}_{1}\wedge\mathrm{\mathbf{w}}_{2})+q_{0}s_{0}\kappa(w_{1x}-w_{2x})\Big{)}.$
(27)

Here $\wedge$ is the wedge product
$(w_{x},w_{y})\wedge(v_{x},v_{y})=w_{x}v_{y}-v_{x}w_{y}$.

##  4 Chern Number

We are now situated to analytically calculate the Chern number from the
related eigenvalue problem
$\bm{\alpha}_{0}(z,\mathrm{\mathbf{k}})\rightarrow\mathrm{e}^{-\mathrm{i}\lambda
z}\bm{\alpha}_{0}(\mathrm{\mathbf{k}})$

$\lambda\bm{\alpha}_{0}+\left(\left\langle
M\right\rangle+\varepsilon\left\langle
N\right\rangle\right)\bm{\alpha}_{0}=0.$ (28)

The number of eigenvalues/eigenvectors depends on the number of sublattices
$L$; e.g, for a honeycomb lattice $M$ is a $2\times 2$ matrix. For a
particular band, $\lambda^{\ell}(\mathrm{\mathbf{k}})$, the eigenvector,
$\bm{\alpha}_{0}^{\ell}(\mathrm{\mathbf{k}})$, parameterized by the
wavenumber, $\mathrm{\mathbf{k}}$, forms a natural vector bundle over the
torus of the Brillouin Zone [7]. The curvature $\Omega(\mathrm{\mathbf{k}})$
of the band (defined below) determines a classification within the de Rham
cohomology of the 2D torus [23]. This is called the Chern class and the Chern
number gives an integer indexing of the possible classes.

In concrete terms, for a particular eigenfunction the Berry curvature is given
by

$\Omega(\bm{k})=\Big{(}\nabla_{\mathrm{\mathbf{k}}}\times\bm{\mathcal{A}}\Big{)}\cdot\widehat{\mathrm{\mathbf{k}}}$
(29)

where

$\bm{\mathcal{A}}(\mathrm{\mathbf{k}})=\bm{\alpha}_{0}^{*}\cdot\nabla_{\mathrm{\mathbf{k}}}\bm{\alpha}_{0},$
(30)

$*$ represents complex conjugate and
$\nabla_{\mathrm{\mathbf{k}}}=\hat{\mathrm{\mathbf{i}}}\partial_{k_{x}}+\hat{\mathrm{\mathbf{j}}}\partial_{k_{y}}$.
$\bm{\mathcal{A}}(\mathrm{\mathbf{k}})$ is the Berry connection and
$\bm{\alpha}_{0}$ is chosen such that $\lVert\bm{\alpha}_{0}\rVert=1$. The
Chern class of the vector bundle is in fact independent of the choice of
connection. For any connection on the vector bundle, the associated curvature
is determined from the differential of the connection. While the curvatures
vary depending on the choice of connection, they all belong to the same Chern
class. It is traditional, however, in the field of Topological Photonics to
express the curvature in terms of the adiabatic theory developed by Berry [8].

The Chern number is computed from the curvature as

$C=\frac{1}{2\pi\mathrm{i}}\int_{\mathcal{B}}\Omega(\mathrm{\mathbf{k}})~{}\mathrm{d}k_{x}\mathrm{d}k_{y}$
(31)

where $\mathcal{B}$ is the Brillouin zone (reciprocal unit cell). When the
Berry connection is a continuously differentiable function of the wavenumber
$\mathrm{\mathbf{k}}$ over the entire Brillouin zone, Stokes’ theorem may be
directly applied to find a zero Chern number,

$C=\frac{1}{2\pi\mathrm{i}}\oint_{\partial\mathcal{B}}\bm{\mathcal{A}}\cdot\mathrm{d}\mathrm{\mathbf{k}}=0,$
(32)

since the connection, $\bm{\mathcal{A}}$, has period matching that of
$\mathcal{B}$. This means nonzero Chern numbers are associated with
discontinuities in the connection which in turn are associated with
discontinuities which appear in the phase of of the eigenfunction
$\bm{\alpha}_{0}$.

Global information about the curvature is found from local behavior of the
Berry connection by isolating and removing discontinuities of the spectrum
from the region of integration, $\mathcal{B}$. This can be done without
effecting the Chern number since the curvature is bounded and vanishing small
regions in the domain of integration correspond to vanishing small
contributions to the integral. The Berry curvature is gauge invariant,
however, this not true of the connection where a transformation in the
eigenfunction
$\bm{\alpha}_{0}(\mathrm{\mathbf{k}})\rightarrow\mathrm{e}^{\mathrm{i}\chi(\mathrm{\mathbf{k}})}\bm{\alpha}_{0}(\mathrm{\mathbf{k}})$
results in a corresponding transform in the connection,
$\mathcal{A}\rightarrow\mathcal{A}+i\nabla_{\bm{k}}\chi(\bm{k})$. As a
consequence, the discontinuities in the connection may be arbitrarily
positioned within the Brillouin zone. To remove this ambiguity the
eigenfunction can be chosen to have a constant phase in one component. The
location of discontinuities will now be uniquely defined and may be dealt with
analytically.

Alternatively, the Chern number may be obtained numerically using the method
presented in [24]. This method similarly relies on a fixed choice of gauge,
however it does not require knowledge of where the discontinuities of the
Berry phase lie. Thus, for higher order systems where a closed form for the
eigenvectors can’t be easily obtained, a numerical calculation for the Chern
number will suffice. Nonetheless, we note that in principal the calculation of
the Chern number can be carried out analytically as demonstrated in Section
4.2.

###  4.1 Honeycomb Lattice

The honeycomb lattice, named for its hexagonal tiling, is well known for its
connections to graphene [25] and for being the simplest lattice model giving
rise to Dirac points [26]. It can be thought of as two interlaced triangular
sublattices. (See Figure 1). For the honeycomb lattice, there are two Dirac
points in the Brillouin zone. Thus, there are two possible locations for
discontinuities to naturally occur in the eigenvectors. Although this doesn’t
effect the Chern number, the expression for the eigenvectors does have
consequences for the Berry connection and the local behavior being used to
arrive at the Chern number analytically.

Figure 1: The honeycomb lattice consists of two triangular sublattices. Blue
spheres represent the first sublattice and purple squares represent the second
sublattice. Dotted lines outline the lattice cells.

(a) Spectrum for the honeycomb lattice with $\rho=1$.

(b) Spectrum for the honeycomb lattice with $\rho=0.5$.

Figure 2: Bulk spectral surfaces for the honeycomb lattice over one Brillouin
zone.

In terms of the framework of Section 1.1, the Honeycomb lattice is defined by
$\mathrm{\mathbf{v}}_{1}=\left(\frac{3}{2},~{}\frac{\sqrt{3}}{2}\right)$,
$\mathrm{\mathbf{v}}_{2}=\left(\frac{3}{2},-\frac{\sqrt{3}}{2}\right)$ and
$\mathrm{\mathbf{d}}_{1}=(1,0)$, $\mathrm{\mathbf{d}}_{2}=(2,0)$ and
$\mathrm{\mathbf{w}}_{1}=\left(-\frac{1}{2},-\frac{\sqrt{3}}{2}\right)$,
$\mathrm{\mathbf{w}}_{2}=\left(-\frac{1}{2},\frac{\sqrt{3}}{2}\right)$,
$\mathrm{\mathbf{w}}_{3}=\left(1,0\right)$, where the distance between nearest
neighbors is taken to be one. In the wavenumber domain of equation (13), the
system for our honeycomb lattice is

$\mathrm{i}\frac{\mathrm{d}\bm{\alpha}}{\mathrm{d}z}+\begin{bmatrix}0&\gamma^{*}(z,\mathrm{\mathbf{k}})\\\
\gamma(z,\mathrm{\mathbf{k}})&0\end{bmatrix}\bm{\alpha}=0$ (33)

where

$\gamma(z,\mathrm{\mathbf{k}})=q(z,\mathrm{\mathbf{w_{3}}})+\rho\Big{(}q(z,\mathrm{\mathbf{w_{1}}})\mathrm{e}^{\mathrm{i}\mathrm{\mathbf{v}}_{1}\cdot\mathrm{\mathbf{k}}}+q(z,\mathrm{\mathbf{w_{2}}})\mathrm{e}^{\mathrm{i}\mathrm{\mathbf{v}}_{2}\cdot\mathrm{\mathbf{k}}}\Big{)}$
(34)

and $q(z,\mathrm{\mathbf{w}})$ is the coefficient function (15) defined
earlier. The parameter $\rho$ represents a stretching/contracting of the
lattice in the $y$-direction, in which case the strength of the connection
along the vectors $\mathrm{\mathbf{w}}_{1}$ and $\mathrm{\mathbf{w}}_{2}$
differs from the strength of the connection along the vector
$\mathrm{\mathbf{w}}_{3}$. This amounts to a breaking of the $120^{\circ}$
rotational symmetry of the lattice.

Substituting the matrix from (33) into the general averaged system (24), we
have the averaged matrix

$\left\langle
M\right\rangle=\begin{bmatrix}0&\langle\gamma\rangle^{*}(\mathrm{\mathbf{k}})\\\
\langle\gamma\rangle(\mathrm{\mathbf{k}})&0\end{bmatrix}$

where

$\left\langle\gamma\right\rangle(\mathrm{\mathbf{k}})=q_{0}\left(1+\rho(\mathrm{e}^{\mathrm{i}\mathrm{\mathbf{v}}_{1}\cdot\mathrm{\mathbf{k}}}+\mathrm{e}^{\mathrm{i}\mathrm{\mathbf{v}}_{2}\cdot\mathrm{\mathbf{k}}})\right)$
(35)

and the constant $q_{0}=\left\langle q(z,\mathrm{\mathbf{w}}\right\rangle)$
for all $\mathrm{\mathbf{w}}$ is the same as defined in Section 2 and can be
factored out. The order $\varepsilon$ correction term is given by

$\left\langle
N\right\rangle(\mathrm{\mathbf{k}})=\begin{bmatrix}\nu(\mathrm{\mathbf{k}})&0\\\
0&-\nu(\mathrm{\mathbf{k}})\end{bmatrix}$ (36)

where

$\displaystyle\nu(\mathrm{\mathbf{k}})=\rho
Q(\mathrm{\mathbf{w}}_{1},-\mathrm{\mathbf{w}}_{3})$
$\displaystyle\mathrm{e}^{\mathrm{i}\mathrm{\mathbf{k}}\cdot\mathrm{\mathbf{v}}_{1}}+\rho^{2}Q(\mathrm{\mathbf{w}}_{2},-\mathrm{\mathbf{w}}_{1})\mathrm{e}^{\mathrm{i}\mathrm{\mathbf{k}}\cdot(\mathrm{\mathbf{v}}_{2}-\mathrm{\mathbf{v}}_{1})}$
(37) $\displaystyle+\rho
Q(\mathrm{\mathbf{w}}_{3},-\mathrm{\mathbf{w}}_{2})\mathrm{e}^{-\mathrm{i}\mathrm{\mathbf{k}}\cdot\mathrm{\mathbf{v}}_{2}}+c.c.$

and “c.c.” stands for the complex conjugate, so $\nu(\mathrm{\mathbf{k}})$ is
a real valued function. Combining equations (35) and (36) we arrive at the
eigenvalue problem of the averaged system from equation (28)

$\lambda\bm{\alpha}_{0}+\begin{bmatrix}\varepsilon\nu(\mathrm{\mathbf{k}})&\left\langle\gamma\right\rangle^{*}(\mathrm{\mathbf{k}})\\\
\left\langle\gamma\right\rangle(\mathrm{\mathbf{k}})&-\varepsilon\nu(\mathrm{\mathbf{k}})\end{bmatrix}\bm{\alpha}_{0}=0.$
(38)

with eigenvalues

$\lambda(\mathrm{\mathbf{k}})=\pm\sqrt{\lvert\left\langle\gamma\right\rangle\rvert^{2}+\varepsilon^{2}\lvert\nu\rvert^{2}}$
(39)

This averaged equation for the honeycomb lattice is equivalent to the well-
known Haldane model [3] in Fourier space; this is due to the fact that the
first correction in the averaging introduces the next-nearest neighbor terms
that are used to break time reversal symmetry in Haldane’s model. In terms of
the Pauli matrices, equation (38) can be written as

$\displaystyle\lambda I\bm{\alpha}$
$\displaystyle+\sigma_{x}q_{0}\Big{(}1+\rho\cos(\mathrm{\mathbf{k}}\cdot\mathrm{\mathbf{v}}_{1})+\rho\cos(\mathrm{\mathbf{k}}\cdot\mathrm{\mathbf{v}}_{2})\Big{)}$
(40)
$\displaystyle+\sigma_{y}q_{0}\Big{(}1+\rho\sin(\mathrm{\mathbf{k}}\cdot\mathrm{\mathbf{v}}_{1})+\rho\sin(\mathrm{\mathbf{k}}\cdot\mathrm{\mathbf{v}}_{2})\Big{)}$
$\displaystyle-\frac{1}{\omega}\sigma_{z}s_{0}^{2}\kappa^{2}\frac{\sqrt{3}}{2}\Big{(}\rho\sin(\mathrm{\mathbf{k}}\cdot\mathrm{\mathbf{v}}_{1})+\rho^{2}\sin(\mathrm{\mathbf{k}}\cdot(\mathrm{\mathbf{v}}_{2}-\mathrm{\mathbf{v}}_{1}))+\rho\sin(-\mathrm{\mathbf{k}}\cdot\mathrm{\mathbf{v}}_{2})\Big{)}$
$\displaystyle+\frac{1}{\omega}\sigma_{z}q_{0}s_{0}\kappa\Big{(}\rho\cos(\mathrm{\mathbf{k}}\cdot\mathrm{\mathbf{v}}_{1})-2\rho^{2}\cos(\mathrm{\mathbf{k}}\cdot(\mathrm{\mathbf{v}}_{2}-\mathrm{\mathbf{v}}_{1}))+\rho\cos(-\mathrm{\mathbf{k}}\cdot\mathrm{\mathbf{v}}_{2})\Big{)}=0.$

This corresponds to the Haldane model with zero mass, i.e. inversion symmetry.

Below we give the analytic step for calculating the Chern number within the
framework of Section 2 and this paper in general. The eigenvectors of equation
(38) are

$\bm{\alpha}_{0}=\frac{1}{D(\mathrm{\mathbf{k}})}\begin{bmatrix}\left\langle\gamma\right\rangle^{*}\\\
-\left(\lambda+\varepsilon\nu\right)\end{bmatrix}$ (41)

where $D(\mathrm{\mathbf{k}})$ is a normalizing term such that
$\bm{\alpha}_{0}^{*}\cdot\bm{\alpha}_{0}=1$; it is given by

$\displaystyle D(\mathrm{\mathbf{k}})^{2}$
$\displaystyle=\left|\left\langle\gamma\right\rangle\right|^{2}+\left(\lambda+\varepsilon\nu\right)^{2}$
$\displaystyle=\left|\left\langle\gamma\right\rangle\right|^{2}+\lambda^{2}+2\varepsilon\lambda\nu+\varepsilon^{2}\nu$
$\displaystyle=2\lambda\left(\lambda+\varepsilon\nu\right)$ (42)

Note: As a result of the longitudinal driving ($\varepsilon\not=0$), the
honeycomb spectrum in equation (38) is no longer conical at the Dirac point.
Furthermore, at the Dirac point, where
$\left\langle\gamma\right\rangle(\mathrm{\mathbf{K}}_{D})=0$, there is a gap
width of $2\lambda(\mathrm{\mathbf{K}}_{D})=2\varepsilon\nu$.

(a) Magnitude and phase for
the first component of ${\bm{\alpha}}_{0}$.

(b) Magnitude for the second component of ${\bm{\alpha}}_{0}$.

Figure 3: Contour of integration that removes the critical points at
$\mathrm{\mathbf{K}}_{D}=\left(0,\frac{4\pi}{3\sqrt{3}}\right)$ and
$\mathrm{\mathbf{K}}^{\prime}_{D}=\left(0,-\frac{4\pi}{3\sqrt{3}}\right)$. The
magnitudes and phases are shown for the components of $\bm{\alpha}_{0}$ (41)
on the upper band. The magnitudes are represented by color saturation (with
white representing a magnitude of $0$) and the phases are represented by hue
(with the second component being strictly real and shown in grey scale).

At any point where
$\left\langle\gamma\right\rangle(\mathrm{\mathbf{K}}_{D})=0$, and therefore
$\lambda=\pm\varepsilon\left|\nu\right|$, the vanishing of the normalizing
function leads to a discontinuity in the eigenvector. So, at a critical point
$\mathrm{\mathbf{K}}_{D}$,

$D(\mathrm{\mathbf{K}}_{D})^{2}=2(\varepsilon\nu)^{2}\left(1+\mathrm{sgn}(\lambda\nu)\right)$
(43)

and thus $\bm{\alpha}_{0}$ has an integrable singularity if
$\mathrm{sgn}(\nu)=-\mathrm{sgn}(\lambda)$.

As an example, we take $\rho=1$. The critical points now occur at
$\mathrm{\mathbf{K}}_{D}=\left(0,\,\frac{4\pi}{\sqrt{3}}\right)$ and
$\mathrm{\mathbf{K}}^{\prime}_{D}=\left(0,-\frac{4\pi}{\sqrt{3}}\right)$. We
look for the local behavior of the eigenvectors near one of these points, as
well as the Berry connection and curvature which are derived from them. For
illustrative proposes, we take the expansion near

$\mathrm{\mathbf{k}}=\mathrm{\mathbf{K}}_{D}+\frac{1}{q_{0}}\delta\mathrm{\mathbf{k}}$
(44)

where $\left\lVert\delta\mathrm{\mathbf{k}}\right\rVert\ll 1$ and the
expansion around $\bm{K}_{D}^{\prime}$ follows similarly. Approximations for
the key terms are

$\displaystyle\left\langle\gamma\right\rangle(\mathrm{\mathbf{k}})$
$\displaystyle\approx-\frac{3}{2}(\mathrm{i}\delta k_{x}+\delta k_{y})$ (45a)
$\displaystyle\nu(\mathrm{\mathbf{k}})$
$\displaystyle\approx\nu(\mathrm{\mathbf{K}}_{D})\equiv\nu_{0}$ (45b)
$\displaystyle\lambda(\mathrm{\mathbf{k}})$
$\displaystyle\approx\pm\sqrt{\frac{9}{4}(\delta k_{x}^{2}+\delta
k_{y}^{2})+\varepsilon^{2}\nu_{0}^{2}}$ (45c)

where $\nu_{0}\in\mathbb{R}$. Substituting these approximations into
$\bm{\alpha}_{0}$ gives

$\bm{\alpha}_{0}(\mathrm{\mathbf{k}})\approx-\frac{1}{\sqrt{2\lambda}}\begin{bmatrix}\displaystyle\frac{\frac{3}{2}(\delta
k_{y}-\mathrm{i}\delta k_{x})}{\sqrt{\lambda+\varepsilon\nu_{0}}}\\\
\sqrt{\lambda+\varepsilon\nu_{0}}\end{bmatrix}.$ (46)

The local approximation of the Berry connection is

$\bm{\mathcal{A}}(\mathrm{\mathbf{k}})\approx\displaystyle\frac{9\mathrm{i}}{8\lambda(\lambda+\varepsilon\nu_{0})}(-\delta
k_{y}\hat{\mathrm{\mathbf{i}}}+\delta k_{x}\hat{\mathrm{\mathbf{j}}}).$ (47)

This shows that the discontinuity from the eigenvector carries through to the
connection in the form of a singularity that behaves like
$\frac{1}{\lVert\delta\mathrm{\mathbf{k}}\rVert}$ when
$\mathrm{sgn}(\lambda\nu_{0})=-1$. By a similar set of calculation, we find
that a potential singularity exists at $\mathrm{\mathbf{K}}^{\prime}_{D}$ as
well. The approximation of the Berry curvature is

$\Omega(\mathrm{\mathbf{k}})\approx-\varepsilon\displaystyle\frac{9\nu_{0}\mathrm{i}}{8\lambda^{3}}{\bm{\hat{\textnormal{k}}}}.$
(48)

Note that $|\lambda|>0$ and so the curvature is bounded. Due to the radial
symmetry of $\bm{\mathcal{A}}\cdot\,\mathrm{d}\mathrm{\mathbf{k}}$ around the
critical points, it’s convenient to define the disk regions

$\mathcal{D}(\bm{\delta})=\left\\{\bm{k}\middle|\left\lVert\bm{K}_{D}-\bm{k}\right\rVert\leq\delta\right\\},\mathcal{D}^{\prime}(\bm{\delta})=\left\\{\bm{k}\middle|\left\lVert\bm{K}^{\prime}_{D}-\bm{k}\right\rVert\leq\delta\right\\}$
(49)

with the punctured Brillouin zone

$\iint_{\mathcal{B}}=\displaystyle\lim_{\delta\rightarrow
0}\iint_{\mathcal{B}-\mathcal{D}(\bm{\delta})-\mathcal{D}^{\prime}(\bm{\delta})}.$
(50)

In this way, we may carve out a small region around any potential
singularities without effecting the the calculation of the curvature or the
Chern number. Once all singularities have been excised from the Brillouin
zone, by removing infinitesimal regions around $\mathrm{\mathbf{K}}_{D}$ and
$\mathrm{\mathbf{K}}^{\prime}_{D}$, we can apply Stokes theorem to represent
global information about the curvature in terms of the local behavior of
connection near the critical points

$\iint_{\mathcal{B}}\Omega(\mathrm{\mathbf{k}})\cdot{\bm{\hat{\textnormal{k}}}}\,\mathrm{d}k_{x}\mathrm{d}k_{y}=\oint_{\partial\mathcal{B}}\mathcal{A}(\mathrm{\mathbf{k}})\,\mathrm{d}\mathrm{\mathbf{k}}-C_{\mathrm{\mathbf{K}}_{D}}-C_{\mathrm{\mathbf{K}}_{D}^{\prime}}$
(51)

where

$\displaystyle C_{\mathrm{\mathbf{K}}_{D}}$
$\displaystyle=\lim_{\delta\rightarrow
0}\oint_{\partial\mathcal{D}(\delta)}\mathcal{A}(\mathrm{\mathbf{k}})\,\mathrm{d}\mathrm{\mathbf{k}}$
(52a) $\displaystyle C_{\mathrm{\mathbf{K}}^{\prime}_{D}}$
$\displaystyle=\lim_{\delta\rightarrow
0}\oint_{\partial\mathcal{D}^{\prime}(\delta)}\mathcal{A}(\mathrm{\mathbf{k}})\,\mathrm{d}\mathrm{\mathbf{k}}.$
(52b)

An illustration of the contours taken can be found in figure 3. Here we take
$\bm{\alpha}_{0}$ on the upper band. At $\mathrm{\mathbf{K}}_{D}$ the first
component of the eigenvector approaches a magnitude of 1 and exhibits a $2\pi$
change in phase around this point. This is the discontinuity that must be
integrated around. By contrast, the first component of the eigenvector is zero
at $\mathrm{\mathbf{K}}^{\prime}_{D}$ and same $2\pi$ phase change no longer
corresponds to a discontinuity. On the lower band this is reversed, which is
why we have chosen to remove both points from the domain. From the periodicity
of the eigenvectors in $\mathrm{\mathbf{k}}$, the boundary integral around
$\partial\mathcal{B}$ is zero.

Define
$\nu_{0}=\nu(\mathrm{\mathbf{K}}_{D})=-\nu(\mathrm{\mathbf{K}}^{\prime}_{D})$.
The Chern number is found to be

$\displaystyle C=$
$\displaystyle\frac{1}{2\pi\mathrm{i}}\left(-C_{\mathrm{\mathbf{K}}_{D}}-C_{\mathrm{\mathbf{K}}^{\prime}_{D}}\right)$
$\displaystyle=$ $\displaystyle\displaystyle\lim_{\bm{\delta}\rightarrow
0}\Bigg{[}-\left(\frac{1}{2}-\mathrm{sgn}(\lambda)\displaystyle\frac{\varepsilon\nu_{0}}{\sqrt{9\delta^{2}+4(\varepsilon\nu_{0})^{2}}}\right)$
$\displaystyle\quad+\left(\frac{1}{2}+\mathrm{sgn}(\lambda)\displaystyle\frac{\varepsilon\nu_{0}}{\sqrt{9\delta^{2}+4(\varepsilon\nu_{0})^{2}}}\right)\Bigg{]}$
$\displaystyle=$ $\displaystyle\,\mathrm{sgn}(\lambda\nu_{0}).$

For $\rho=1$, the Chern number is $-1$ for the upper band, $\lambda>0$, and
$+1$ for the lower band, $\lambda<0$. This is consistent with results from
[16, 27]. This topology persists for $\rho>1/2$. When $\rho=1/2$ the two
critical points merge at the corner of the Brillouin zone, i.e.,
$\mathrm{\mathbf{K}}_{D}=\mathrm{\mathbf{K}}^{\prime}_{D}=\left(0,\pm\frac{4\pi}{\sqrt{3}}\right)$,
as seen in Figure (2(b)). At this point $\nu(\mathrm{\mathbf{K}}_{D})=0$ and
so the gap in the spectrum caused by the driving of the lattice closes and a
Dirac point forms. For $\rho<1/2$, a gap opens that is caused by the
stretching of the lattice instead of the driving of the lattice. In the
analysis, this corresponds to $|\gamma|>0$ for all $\mathrm{\mathbf{k}}$ in
the Brillouin zone, so that no discontinuities exist in the phase of the
eigenvectors and the Chern Number on both bands is zero.

###  4.2 Lieb Lattice

Here we consider the Lieb lattice [28, 29] as an example. The Lieb lattice is
a 1/4-depleted square lattice, i.e., one in four lattice sites have been
removed, consisting of three sublattice sites in each unit cell. The Lieb
system features a single spectral degeneracy within the Brillouin zone where
upper and lower bands intersect with a middle flat-band at a Dirac point (see
Figure 4(b)). While the three distinct sub-lattices result in a 3rd order
system in $\mathrm{\mathbf{k}}$-space, the lone Dirac point streamlines the
calculations of the Chern number. As we will show, the Dirac point lines up
with a singularity in the Berry connection that occurs on both the upper and
lower bands.

(a) Diagram for the Lieb lattice with three lattice sites per unit cell.

(b) Spectrum for the Lieb lattice.

Figure 4: The Lieb Lattice consists of three square sublattices. The green
triangles represent the first sublattice, the purple squares represent the
second sublattice and the teal spheres represent the third sublattice.

In terms of the framework laid out in Section 1.1, the Lieb lattice is defined
by $\mathrm{\mathbf{v}}_{1}=(2,0)$, $\mathrm{\mathbf{v}}_{2}=(0,2)$ and
$\mathrm{\mathbf{d}}_{1}=(0,1)$, $\mathrm{\mathbf{d}}_{2}=(1,1)$,
$\mathrm{\mathbf{d}}_{3}=(1,0)$ and $\mathrm{\mathbf{w}}_{1}=(1,0)$,
$\mathrm{\mathbf{w}}_{2}=(0,1)$. The system of differential equations, in the
form (13), derived for this lattice is

$\mathrm{i}\frac{\mathrm{d}\bm{\alpha}}{\mathrm{d}z}+\begin{bmatrix}0&\gamma_{1}^{*}(z,\mathrm{\mathbf{k}})&0\\\
\gamma_{1}(z,\mathrm{\mathbf{k}})&0&\gamma_{2}^{*}(z,\mathrm{\mathbf{k}})\\\
0&\gamma_{2}(z,\mathrm{\mathbf{k}})&0\end{bmatrix}\bm{\alpha}=0$ (53)

where the $\gamma$ functions are given by

$\displaystyle\gamma_{1}(z,\mathrm{\mathbf{k}})$
$\displaystyle=q(z,\mathrm{\mathbf{w}}_{1})\mathrm{e}^{\mathrm{i}2k_{x}}+q(z,-\mathrm{\mathbf{w}}_{1})$
(54a) $\displaystyle\gamma_{2}(z,\mathrm{\mathbf{k}})$
$\displaystyle=q(z,\mathrm{\mathbf{w}}_{2})\mathrm{e}^{\mathrm{i}2k_{y}}+q(z,-\mathrm{\mathbf{w}}_{2})$
(54b)

and $q$ as defined in equation (15).

Looking at Figure 4, we see that there is no interaction between the first and
third sublattice, which correspond to lattice sites positioned between the
lattice sites of the square lattice created by the second sublattice. This
lack of interaction is reflected in the zeros found in the upper right and
lower left corners of the matrix from equation (53). A detailed derivation of
the specific $\gamma$ functions given here as well as a more generalized
version of this system appears in [16].

Substituting the matrix from the Lieb system (53) into the general averaged
equation (24) gives the averaged matrix

$\left\langle
M\right\rangle(\mathrm{\mathbf{k}})=\begin{bmatrix}0&\left\langle\gamma_{1}\right\rangle^{*}(\mathrm{\mathbf{k}})&0\\\
\left\langle\gamma_{1}\right\rangle(\mathrm{\mathbf{k}})&0&\left\langle\gamma_{2}\right\rangle^{*}(\mathrm{\mathbf{k}})\\\
0&\left\langle\gamma_{2}\right\rangle(\mathrm{\mathbf{k}})&0\end{bmatrix}.$
(55)

For the averaged $\gamma_{i}$ functions, we again use the fact that
$\left\langle q(\zeta,\mathrm{\mathbf{w}})\right\rangle=q_{0}$. This gives

$\displaystyle\left\langle\gamma_{1}\right\rangle(\mathrm{\mathbf{k}})$
$\displaystyle=q_{0}\left(\mathrm{e}^{\mathrm{i}2k_{x}}+1\right)$ (56a)
$\displaystyle\left\langle\gamma_{2}\right\rangle(\mathrm{\mathbf{k}})$
$\displaystyle=q_{0}\left(\mathrm{e}^{\mathrm{i}2k_{y}}+1\right)$ (56b)

At order $\varepsilon$, the correction term is

$\left\langle
N\right\rangle(\mathrm{\mathbf{k}})=\begin{bmatrix}0&0&-\mathrm{i}\varepsilon\,\nu^{*}(\mathrm{\mathbf{k}})\\\
0&0&0\\\ \,\mathrm{i}\varepsilon\,\nu(\mathrm{\mathbf{k}})&0&0\end{bmatrix}$
(57)

where

$\displaystyle\mathrm{i}\nu(\mathrm{\mathbf{k}})=$
$\displaystyle~{}Q(\mathrm{\mathbf{w}}_{1},-\mathrm{\mathbf{w}}_{2})+Q(-\mathrm{\mathbf{w}}_{1},\mathrm{\mathbf{w}}_{2})\mathrm{e}^{\mathrm{i}\mathrm{\mathbf{k}}\cdot(\mathrm{\mathbf{v}}_{2}-\mathrm{\mathbf{v}}_{1})}$
(58)
$\displaystyle+Q(\mathrm{\mathbf{w}}_{1},\mathrm{\mathbf{w}}_{2})\mathrm{e}^{\mathrm{i}\mathrm{\mathbf{k}}\cdot\mathrm{\mathbf{v}}_{2}}+Q(-\mathrm{\mathbf{w}}_{1},-\mathrm{\mathbf{w}}_{2})\mathrm{e}^{-\mathrm{i}\mathrm{\mathbf{k}}\cdot\mathrm{\mathbf{v}}_{1}}$

and $Q$ is defined by equation (27). The nonzero terms in (57) correspond to
next-nearest neighbor interactions between the first and third lattice sites.
The eigenvalue problem for approximating the spectrum of the Lieb lattice is
now found to be

$\lambda\bm{\alpha}_{0}+\begin{bmatrix}0&\left\langle\gamma_{1}\right\rangle^{*}&-\mathrm{i}\varepsilon\nu^{*}\\\
\left\langle\gamma_{1}\right\rangle&0&\left\langle\gamma_{2}\right\rangle^{*}\\\
\mathrm{i}\varepsilon\nu&\left\langle\gamma_{2}\right\rangle&0\end{bmatrix}\bm{\alpha}_{0}=0.$
(59)

This systems has the eigenvalues

$\lambda=0,\pm\sqrt{\left|\left\langle\gamma_{1}\right\rangle\right|^{2}+\left|\left\langle\gamma_{2}\right\rangle\right|^{2}+\varepsilon^{2}\left|\nu\right|^{2}}.$
(60)

Note: The constant term in the characteristic polynomial,
$\left\langle\gamma_{1}\right\rangle\left\langle\gamma_{2}\right\rangle\nu^{*}-\left\langle\gamma_{1}\right\rangle^{*}\left\langle\gamma_{2}\right\rangle^{*}\nu$,
is zero and we find that
$\left\langle\gamma_{1}\right\rangle\left\langle\gamma_{2}\right\rangle\nu^{*}$
is strictly real. The eigenvectors can be expressed in terms of
$\lambda(\mathrm{\mathbf{k}})$ as

$\bm{\alpha}_{0}=\displaystyle\frac{1}{D(\mathrm{\mathbf{k}})}\begin{bmatrix}\lambda\left\langle\gamma_{1}\right\rangle^{*}+\mathrm{i}\varepsilon\nu^{*}\left\langle\gamma_{2}\right\rangle\\\
-\lambda^{2}+\varepsilon^{2}\left|\nu\right|^{2}\\\
\lambda\left\langle\gamma_{2}\right\rangle-\mathrm{i}\varepsilon\nu\left\langle\gamma_{1}\right\rangle^{*}\end{bmatrix}$
(61)

where $D(\mathrm{\mathbf{k}})$ is a normalizing function

$\displaystyle D(\mathrm{\mathbf{k}})^{2}=$
$\displaystyle\lambda^{2}(2\left|\left\langle\gamma_{1}\right\rangle\right|^{2}+2\left|\left\langle\gamma_{2}\right\rangle\right|^{2}-\varepsilon^{2}\left|\nu\right|^{2})$
$\displaystyle+\varepsilon^{2}\left|\nu\right|^{2}(\left|\left\langle\gamma_{1}\right\rangle\right|^{2}+\left|\left\langle\gamma_{2}\right\rangle\right|^{2}+\varepsilon^{2}\left|\nu\right|^{2})$
(62)

which makes $\lVert\bm{\alpha}_{0}\rVert=1$. In this case, we have chosen to
fix the second component of $\bm{\alpha}_{0}$ as strictly real. When not on
the flat band, i.e., $\lambda=0$, $D$ further reduces to
$D(\mathrm{\mathbf{k}})^{2}=2\lambda^{2}(\left|\left\langle\gamma_{1}\right\rangle\right|^{2}+\left|\left\langle\gamma_{2}\right\rangle\right|^{2})$.
Notice that $D(\mathrm{\mathbf{k}})$ vanishes when
$\left\langle\gamma_{1}\right\rangle=\left\langle\gamma_{2}\right\rangle=0$,
which corresponds to the Dirac point of the unperturbed spectrum, i.e.,
$\varepsilon=0$. This has the effect of creating a discontinuity from the
twist in the phase around the zeros of the first and third components of the
normalized eigenvector (61) (see Figure 5). In terms of the magnitude, the
zeros of $D$ are removable discontinuities and the important discontinuity
appears only in the phase of the eigenvector. As such, there is a $2\pi$
change in the eigenfunction or so called Berry phase for a loop around the
discontinuity.

(a) Magnitude and phase for the first component of $\bm{\alpha}_{0}$.

(b) Magnitude and phase for the second component of $\bm{\alpha}_{0}$.

(c) Magnitude and phase for third component of $\bm{\alpha}_{0}$.

Figure 5: Contour of integration that removes the critical point at
$\bm{K}_{D}=\left(\frac{\pi}{2},\frac{\pi}{2}\right)$. The magnitudes and
phases are shown for the components of $\bm{\alpha}_{0}(\bf k)$ on the upper
band. The magnitudes are represented by color saturation (with white
representing a magnitude of $0$) and the phases are represented by hue (with
the second component of $\bm{\alpha}_{0}$ being strictly real and shown in
grey scale).

The discontinuity in the spectrum occurs at
$\bm{K}_{D}=\left(\frac{\pi}{2},\frac{\pi}{2}\right)$, which means our
calculation of the Chern number will be determined by the local behavior near
this point. For
$\mathrm{\mathbf{k}}=\mathrm{\mathbf{K}}_{D}+\delta\mathrm{\mathbf{k}}/q_{0}$
we have the approximations

$\displaystyle\left\langle\gamma_{1}\right\rangle(\mathrm{\mathbf{k}})$
$\displaystyle\approx-2\mathrm{i}\,\delta k_{x}$ (63a)
$\displaystyle\left\langle\gamma_{2}\right\rangle(\mathrm{\mathbf{k}})$
$\displaystyle\approx-2\mathrm{i}\,\delta k_{y}$ (63b)
$\displaystyle\nu(\mathrm{\mathbf{k}})$
$\displaystyle\approx\nu(\mathrm{\mathbf{K}}_{D})\equiv\nu_{0}$ (63c)

where $\nu_{0}$ is a real nonzero constant. In this region the set of
eigenvalues are approximated by
$\lambda(\mathrm{\mathbf{k}})\approx\left\\{0,\pm\sqrt{4\delta
k_{x}^{2}+4\delta k_{y}^{2}+\varepsilon^{2}\nu_{0}^{2}}\right\\}$. In terms of
their graphs, these spectral surfaces correspond to a flat plane and a
hyperboloid of two sheets. On the upper and lower bands, the related
eigenvectors are given by

$\bm{\alpha}_{0}(\mathrm{\mathbf{k}})\approx\displaystyle\frac{2}{D(\mathrm{\mathbf{k}})}\begin{bmatrix}-\mathrm{i}\delta
k_{x}\lambda+\delta k_{y}\varepsilon\nu_{0}\\\ -2\left(\delta k_{x}^{2}+\delta
k_{y}^{2}\right)\\\ ~{}~{}\mathrm{i}\delta k_{y}\lambda+\delta
k_{x}\varepsilon\nu_{0}\end{bmatrix}$ (64)

where $D(\mathrm{\mathbf{k}})^{2}\approx 8\left(\delta k_{x}^{2}+\delta
k_{y}^{2}\right)\lambda^{2}$.

From here, the approximation of the Berry connection is found to be

$\mathcal{A}=\mathcal{A}_{1}\hat{\mathrm{\mathbf{i}}}+\mathcal{A}_{2}\hat{\mathrm{\mathbf{j}}}\approx\displaystyle\frac{\varepsilon\nu_{0}\mathrm{i}}{\lambda\left(\delta
k_{x}^{2}+\delta k_{y}^{2}\right)}\left(\delta
k_{y}\hat{\mathrm{\mathbf{i}}}-\delta k_{x}\hat{\mathrm{\mathbf{j}}}\right)$
(65)

and the discontinuity in the eigenvectors is seen to produce a singularity in
the Berry connection. Note: On the flat band where $\lambda=0$ the connection
is zero. The approximation for the Berry curvature is then

$\Omega=(\partial_{k_{x}}\mathcal{A}_{2}-\partial_{k_{y}}\mathcal{A}_{1}){\bm{\hat{\textnormal{k}}}}\approx\displaystyle\frac{\mathrm{i}4\varepsilon\nu_{0}}{\lambda^{3}}{\bm{\hat{\textnormal{k}}}}.$
(66)

Unlike the eigenvectors and the Berry connection, the Berry curvature for the
upper and lower bands are both smooth and bounded since $|\lambda|>0$. As a
result, analogous to what was done in the honeycomb lattice calculation, we
can remove an infinitesimal region around
$\mathrm{\mathbf{k}}=\mathrm{\mathbf{K}}_{D}$ from the Brillouin zone without
effecting the total curvature or Chern number. Once the singular point has
been removed in this manner, Stokes’ theorem can be employed to calculate the
Chern number using the local behavior of the connection instead of the global
behavior of the curvature. The radial symmetry of the differential
$\bm{\mathcal{A}}\cdot\mathrm{d}\mathrm{\mathbf{k}}$ implies that removing a
circular disk
$\mathcal{D}(\delta)=\left\\{\bm{k}\middle|\lVert\bm{k}-\bm{K}_{D}\rVert\leq\delta\right\\}$
will facilitate the direct calculation of all relevant integrals.

The particular steps of this reasoning are outlined below

$\displaystyle C$
$\displaystyle=\displaystyle\frac{1}{2\pi\mathrm{i}}\displaystyle\iint_{\mathcal{B}}\Omega(\mathrm{\mathbf{k}})\,\,\mathrm{d}k_{x}\mathrm{d}k_{y}$
(67a) $\displaystyle=\displaystyle\lim_{\delta\rightarrow
0}\displaystyle\frac{1}{2\pi\mathrm{i}}\displaystyle\iint_{\mathcal{B}-\mathcal{D}(\delta)}\Big{(}\nabla_{\mathrm{\mathbf{k}}}\times\bm{\mathcal{A}}\Big{)}\cdot{\bm{\hat{\textnormal{k}}}}\,\mathrm{d}k_{x}\mathrm{d}k_{y}$
(67b) $\displaystyle=\displaystyle\lim_{\delta\rightarrow
0}\displaystyle\frac{1}{2\pi\mathrm{i}}\left(\oint_{\partial\mathcal{B}}\mathcal{A}(\mathrm{\mathbf{k}})\cdot\mathrm{d}\bm{k}-\displaystyle\oint_{\partial\mathcal{D}(\delta)}\mathcal{A}(\mathrm{\mathbf{k}})\cdot\mathrm{d}\bm{k}\right)$
(67c)
$\displaystyle=\displaystyle\frac{1}{2\pi\mathrm{i}}\displaystyle\lim_{\delta\rightarrow
0}\left(-\displaystyle\oint_{\partial\mathcal{D}(\delta)}\mathcal{A}(\bm{k})\cdot\mathrm{d}\bm{k}\right)$
(67d)
$\displaystyle=\displaystyle\frac{1}{2\pi}\displaystyle\lim_{\delta\rightarrow
0}\left(-\mathrm{sgn}(\lambda)\int_{0}^{2\pi}\displaystyle\frac{\varepsilon\nu_{0}}{\sqrt{4\delta^{2}+(\varepsilon\nu_{0})^{2}}}\mathrm{d}\,\theta\right)$
(67e) $\displaystyle=\displaystyle\lim_{\delta\rightarrow
0}\left(-\mathrm{sgn}(\lambda)\displaystyle\frac{\varepsilon\nu_{0}}{\sqrt{4\delta^{2}+(\varepsilon\nu_{0})^{2}}}\right)$
(67f) $\displaystyle=-\mathrm{sgn}(\lambda\nu_{0}),$ (67g)

where we have used the periodicity of the Berry connection along the boundary
of the Brillouin zone. The end result is a pair of opposite Chern numbers of
$\pm 1$ for the top and bottom bands. Numerical calculation of coefficients in
the function $\nu$ shows that $\nu_{0}>0$ for a wide range of parameters and
the upper band has a Chern number $-1$ while the bottom band has a Chern
number of $+1$. The flat middle band has a Chern number of $0$, which can be
found by direct application of Stokes’ theorem since there are no
discontinuities present in the Berry connection. This is in agreement with
numerical calculations of the spectrum in [16].

###  4.3 Kagome Lattice

Here we consider the kagome lattice (see Figure 6) as an example [30]. In a
similar way to how the Lieb lattice is a 1/4-depleted square lattice, the
kagome lattice is a triangular lattice where one in every four lattice sites
has been removed. In Figure 6, the lattice cell is drawn so that the corners
are located where the lattice sites would have been removed from a triangular
lattice. When considering the behavior in the Brillouin zone the, the kagome
system features two Dirac points between the middle and top bands as well as a
third parabolic degeneracy between the nearly flat bottom band and the middle
band (see Figure 6(b)). While the discontinuities of the eigenfunctions
bifurcate from the Dirac points for small oscillations ($\kappa\ll 1$), they
do not remain stationary as is the case for honeycomb and Lieb lattices.

(a) Diagram for the Kagome lattice and except for a single lattice cell.

(b) Spectrum for the kagome lattice.

Figure 6: The kagome kattice consists of three sublattices. The green
triangles represent the first sublattice, the blue spheres represent the
second sublattice and the purple squares represent the third sublattice.

In terms of the framework laid out in Section 1.1, the kagome lattice is
defined by $\mathrm{\mathbf{v}}_{1}=(2,0)$,
$\mathrm{\mathbf{v}}_{2}=(-1,\sqrt{3})$ (though not used to define the lattice
cell we define $\mathrm{\mathbf{v}}_{3}=(-1,-\sqrt{3})$) and
$\mathrm{\mathbf{d}}_{1}=\left(0,\sqrt{3}\right)$,
$\mathrm{\mathbf{d}}_{2}=\left(-\frac{1}{2},\frac{\sqrt{3}}{2}\right)$,
$\mathrm{\mathbf{d}}_{3}=(\frac{1}{2},\frac{\sqrt{3}}{2})$. From here we can
derive the nearest neighbor vectors $\mathrm{\mathbf{w}}_{1}=(1,0)$,
$\mathrm{\mathbf{w}}_{2}=\left(-\frac{1}{2},\frac{\sqrt{3}}{2}\right)$, and
$\mathrm{\mathbf{w}}_{3}=\left(-\frac{1}{2},-\frac{\sqrt{3}}{2}\right)$. The
system of differential equations, in the form (13), derived for this lattice
is

$\mathrm{i}\frac{\mathrm{d}\bm{\alpha}}{\mathrm{d}z}+\begin{bmatrix}0&\gamma_{3}^{*}(z,\mathrm{\mathbf{k}})&\gamma_{2}(z,\mathrm{\mathbf{k}})\\\
\gamma_{3}(z,\mathrm{\mathbf{k}})&0&\gamma_{1}^{*}(z,\mathrm{\mathbf{k}})\\\
\gamma_{2}^{*}(z,\mathrm{\mathbf{k}})&\gamma_{1}(z,\mathrm{\mathbf{k}})&0\end{bmatrix}\bm{\alpha}=0$
(68)

where the $\gamma_{j}$ and their averages
$\left\langle\gamma_{j}\right\rangle$ functions are given by

$\displaystyle\gamma_{j}(z,\mathrm{\mathbf{k}})$
$\displaystyle=q(z,\mathrm{\mathbf{w}}_{j})\mathrm{e}^{\mathrm{i}\mathrm{\mathbf{k}}\cdot\mathrm{\mathbf{v}}_{j}}+q^{*}(z,\mathrm{\mathbf{w}}_{j})$
(69a) $\displaystyle\left\langle\gamma_{j}\right\rangle(\mathrm{\mathbf{k}})$
$\displaystyle=q_{0}\Big{(}\mathrm{e}^{\mathrm{i}\mathrm{\mathbf{k}}\cdot\mathrm{\mathbf{v}}_{j}}+1\Big{)}$
(69b)

for the function $q$ previously defined in (15). Looking at Figure 6, we see
that all interactions come in pairs pointing in opposite directions.

Plugging (68) into the general averaged equation (24) yields

$\left\langle
M\right\rangle(\mathrm{\mathbf{k}})=\begin{bmatrix}0&\left\langle\gamma_{3}\right\rangle^{*}&\left\langle\gamma_{2}\right\rangle~{}\\\
\left\langle\gamma_{3}\right\rangle~{}&0&\left\langle\gamma_{1}\right\rangle^{*}\\\
\left\langle\gamma_{2}\right\rangle^{*}&\left\langle\gamma_{1}\right\rangle~{}&0\end{bmatrix}.$
(70)

and at order $\varepsilon$, the correction term is

$\left\langle
N\right\rangle(\mathrm{\mathbf{k}})=\begin{bmatrix}0&\nu_{21}(\mathrm{\mathbf{k}})&\nu_{31}(\mathrm{\mathbf{k}})\\\
\nu_{12}(\mathrm{\mathbf{k}})&0&\nu_{32}(\mathrm{\mathbf{k}})\\\
\nu_{13}(\mathrm{\mathbf{k}})&\nu_{31}(\mathrm{\mathbf{k}})&0\end{bmatrix}$
(71)

where

$\nu_{ij}(\mathrm{\mathbf{k}})=\mathrm{e}^{-\mathrm{i}\mathrm{\mathbf{k}}\cdot(\mathrm{\mathbf{d}}_{j}+\mathrm{\mathbf{d}}_{i})}\Bigg{(}Q(\mathrm{\mathbf{w}}_{i},\mathrm{\mathbf{w}}_{j})\mathrm{e}^{\mathrm{i}\mathrm{\mathbf{k}}\cdot(\mathrm{\mathbf{w}}_{i}+\mathrm{\mathbf{w}}_{j})}+Q(-\mathrm{\mathbf{w}}_{i},\mathrm{\mathbf{w}}_{j})\mathrm{e}^{\mathrm{i}\mathrm{\mathbf{k}}\cdot(-\mathrm{\mathbf{w}}_{i}+\mathrm{\mathbf{w}}_{j})}-c.c.\Bigg{)}$

and “c.c.” stands for complex conjugate. We can formulate the averaged
eigenfunctions and the corresponding connection, curvature and the Chern
numbers using the same framework from Section 4.2. The Chern numbers are found
numerically (using the algorithm in [31]) to be $-1$, $0$, and $1$ for the
bottom, middle, and top band respectively ; this agrees with the results in
[16].

###  4.4 1/5 Depleted Lattice

Here we consider the 1/5 depleted lattice [32] as an example. As a 1/5
depleted lattice, one in every five lattice sites is being removed (see e.g.
Figure 7) from a square lattice. The choice of lattice cell is chosen such
that the corners of the cell lie at the location of the removed lattice sites.
The spectrum for the 1/5 depleted lattice has two Dirac points similar to the
kind found for the Lieb lattice in where a flat band sits between upper and
lower sections of the conical spectrum. In this case, the middle band is only
locally flat near the Dirac point, but not across the entire Brouillon zone.
Though they do not touch, there is no gap in $\lambda$ between the lower
middle and upper middle bands.

(a) Diagram for the 1/5 depleted lattice with four lattice sites per unit
cell.

(b) Spectrum for the 1/5 depleted lattice.

Figure 7: The 1/5-depleted Lattice consists of four sublattices. The green
triangles represent the first sublattice, the blue squares represent the
second sublattice, the purple squares represent the third sublattice, and the
red pentagrams represent the fourth sublattice.

In terms of the framework laid out in Section 1.1, the 1/5 depleted lattice is
defined by $\mathrm{\mathbf{v}}_{1}=(1,-2)$, $\mathrm{\mathbf{v}}_{2}=(2,1)$
and $\mathrm{\mathbf{d}}_{1}=(1,-1)$, $\mathrm{\mathbf{d}}_{2}=(2,-1)$,
$\mathrm{\mathbf{d}}_{3}=(2,0)$, $\mathrm{\mathbf{d}}_{4}=(1,0)$. From here we
can derive the nearest neighbor vectors $\mathrm{\mathbf{w}}_{1}=(1,0)$ and
$\mathrm{\mathbf{w}}_{2}=(0,1)$. The system of differential equations, in the
form (13), derived for this lattice is

$\mathrm{i}\frac{\mathrm{d}\bm{\alpha}}{\mathrm{d}z}+\begin{bmatrix}0&q(z,\mathrm{\mathbf{w}}_{1})&\gamma_{1}^{*}(z,\mathrm{\mathbf{k}})&q(z,\mathrm{\mathbf{w}}_{2})\\\
q^{*}(z,\mathrm{\mathbf{w}}_{1})&0&q(z,\mathrm{\mathbf{w}}_{2})&\gamma_{2}^{*}(z,\mathrm{\mathbf{k}})\\\
\gamma_{1}(z,\mathrm{\mathbf{k}})&q^{*}(z,\mathrm{\mathbf{w}}_{2})&0&q^{*}(z,\mathrm{\mathbf{w}}_{1})\\\
q^{*}(z,\mathrm{\mathbf{w}}_{2})&\gamma_{2}(z,\mathrm{\mathbf{k}})&q(z,\mathrm{\mathbf{w}}_{1})&0\end{bmatrix}\bm{\alpha}=0$
(72)

where the $\gamma_{j}$ functions are given explicitly by

$\displaystyle\gamma_{1}(z,\mathrm{\mathbf{k}})$
$\displaystyle=q(z,\mathrm{\mathbf{w}}_{1})\mathrm{e}^{\mathrm{i}\mathrm{\mathbf{k}}\cdot\mathrm{\mathbf{v}}_{2}}$
(73a) $\displaystyle\gamma_{2}(z,\mathrm{\mathbf{k}})$
$\displaystyle=q(z,\mathrm{\mathbf{w}}_{2})\mathrm{e}^{-\mathrm{i}\mathrm{\mathbf{k}}\cdot\mathrm{\mathbf{v}}_{1}}$
(73b)

and we note that $q^{*}(z,\mathrm{\mathbf{w}})=q(z,-\mathrm{\mathbf{w}})$.

It is often a feature of more complicated lattices that there are few
interactions between sublattices. Here for example, each of the four
sublattices of the 1/5 depleted lattice has only one nearest neighbor on each
of the other sublattices, but there were two interactions in the cases of Lieb
and kagome lattices and three for the hexagonal lattice. This means the
individual elements of $M$ become simpler as the order of $M$ increases.

Substituting (72) into the general averaged equation (24) gives

$\left\langle
M\right\rangle(\mathrm{\mathbf{k}})=q_{0}\begin{bmatrix}0&1&\mathrm{e}^{-\mathrm{i}\mathrm{\mathbf{k}}\cdot\mathrm{\mathbf{v}}_{2}}&1\\\
1&0&1&\mathrm{e}^{\mathrm{i}\mathrm{\mathbf{k}}\cdot\mathrm{\mathbf{v}}_{1}}\\\
\mathrm{e}^{\mathrm{i}\mathrm{\mathbf{k}}\cdot\mathrm{\mathbf{v}}_{2}}&1&0&1\\\
1&\mathrm{e}^{-\mathrm{i}\mathrm{\mathbf{k}}\cdot\mathrm{\mathbf{v}}_{1}}&1&0\\\
\end{bmatrix}.$ (74)

At order $\varepsilon$, the correction term is

$\left\langle
N\right\rangle(\mathrm{\mathbf{k}})=\begin{bmatrix}\nu_{0}&-\nu_{1}^{*}~{}\mathrm{e}^{-\mathrm{i}\mathrm{\mathbf{k}}\cdot\mathrm{\mathbf{v}}_{2}}&0&\nu_{1}^{*}~{}\mathrm{e}^{\mathrm{i}\mathrm{\mathbf{k}}\cdot\mathrm{\mathbf{v}}_{1}}\\\
-\nu_{1}\mathrm{e}^{\mathrm{i}\mathrm{\mathbf{k}}\cdot\mathrm{\mathbf{w}}_{2}}&0&-\nu_{1}^{*}\mathrm{e}^{\mathrm{i}\mathrm{\mathbf{k}}\cdot\mathrm{\mathbf{v}}_{1}}&0\\\
0&-\nu_{1}\mathrm{e}^{-\mathrm{i}\mathrm{\mathbf{k}}\cdot\mathrm{\mathbf{v}}_{1}}&-\nu_{0}&\nu_{1}\mathrm{e}^{\mathrm{i}\mathrm{\mathbf{k}}\cdot\mathrm{\mathbf{v}}_{2}}\\\
\nu_{1}\mathrm{e}^{-\mathrm{i}\mathrm{\mathbf{k}}\cdot\mathrm{\mathbf{v}}_{1}}&0&\nu_{1}^{*}\mathrm{e}^{-\mathrm{i}\mathrm{\mathbf{k}}\mathrm{\mathbf{v}}_{2}}&0\end{bmatrix}$
(75)

where $\nu_{0}=2q_{0}s_{0}\kappa$ and
$\nu_{1}=\frac{1}{2}\big{(}\nu_{0}+\mathrm{i}s_{0}^{2}\kappa^{2}\big{)}$. As
before, using the averaging framework and the algorithm in [31] the Chern
numbers are found as $1$, $1$, $-1$ and $-1$ for the bands running from bottom
to top. Due to space limitations we omit intermediate details.

##  5 Nonlinear Envelope Equations

From the linear analysis, we have found that even though the Chern Number is a
global property of the spectrum, it can be successfully calculated from the
local behavior near critical points in the spectrum under suitable
restrictions on the eigenvectors. We now consider the nonlinear dynamics of
wave packets in the gap near the critical points of the linear spectrum. For
this, we assume weak nonlinearity.

$\mathrm{i}\displaystyle\frac{\partial\bm{\alpha}}{\partial
z}+M\left(z,\mathrm{\mathbf{k}}\right)\bm{\alpha}+\varepsilon\widetilde{\sigma}\mathbb{N}(\bm{\alpha})\bm{\alpha}=0.$
(76)

where the nonlinear term once again is defined as
$\mathbb{N}(\bm{\alpha})=\mathrm{diag}(|\alpha^{1}|^{2},\ldots,|\alpha^{\ell}|^{2})$.

Following the steps from Section 3 the averaged nonlinear system when the
nonlinear term first appears is found to be

$\mathrm{i}\frac{\partial\bm{\alpha}}{\partial
z}+\left(\left<M\right>(\mathrm{\mathbf{k}})+\varepsilon\left\langle
N\right\rangle(\mathrm{\mathbf{k}})\right)\bm{\alpha}+\varepsilon\widetilde{\sigma}\mathbb{N}(\bm{\alpha})\bm{\alpha}=0.$
(77)

Note that we are assuming a balance between the nonlinearity and the higher
order term due to the averaging.

Near a critical point of the linear spectrum,
$\mathrm{\mathbf{k}}=\mathrm{\mathbf{K}}_{D}+\varepsilon\mathrm{\mathbf{k}}$,
we can expand the linear terms as

$\left\langle M\right\rangle+\varepsilon\left\langle
N\right\rangle\approx\left\langle
M\right\rangle(\mathrm{\mathbf{K}}_{D})+\varepsilon\mathrm{\mathbf{k}}\cdot\nabla_{\mathrm{\mathbf{k}}}\left\langle
M\right\rangle+\varepsilon\left\langle
N\right\rangle(\mathrm{\mathbf{K}}_{D})+\ldots$ (78)

where the gradient of the matrix $\left\langle M\right\rangle$ is applied
element-wise.

To find an equation for the slowly-varying envelope of a wave packet, we
expand the solution in form

$\bm{\alpha}=\left(\bm{\alpha}_{0}\left(Z,\varepsilon\mathrm{\mathbf{k}}\right)+\varepsilon\bm{\alpha}_{1}\left(Z,\varepsilon\mathrm{\mathbf{k}}\right)+\varepsilon^{2}\bm{\alpha}_{2}\left(Z,\varepsilon\mathrm{\mathbf{k}}\right)+\ldots\right)\mathrm{e}^{-\mathrm{i}\lambda_{0}z}.$
(79)

where $Z=\varepsilon z$ is a slow propagation variable. At leading order,
$\bm{\alpha}_{0}$ satisfies the eigenvalue problem

$\lambda_{0}\bm{\alpha}_{0}+\left\langle
M\right\rangle(\mathrm{\mathbf{K}}_{D})\bm{\alpha}_{0}=0.$ (80)

Equation (80) represents the unperturbed, $\varepsilon=0$, lattice at the
critical point, $\mathrm{\mathbf{K}}_{D}$, corresponds to a Dirac point of the
linear spectrum. We assume linearly independent solutions (two in the case of
Dirac , three in the case of the Lieb and kagome lattices). Superposition of
$L$ slowly-varying envelopes, one for each eigenvector, yields

$\bm{\alpha}_{0}\left(Z,\varepsilon\mathrm{\mathbf{k}}\right)=\sum_{i=1}^{L}\widehat{E}_{i}\left(Z,\varepsilon\mathrm{\mathbf{k}}\right)\bm{\alpha}^{i}_{0}\left(\varepsilon\mathrm{\mathbf{k}}\right).$
(81)

At $\mathrm{O}(\varepsilon)$, the nonlinear term balances against the linear
terms of the band approximation

$\displaystyle\lambda_{0}\bm{\alpha}_{1}+\left\langle M\right\rangle$
$\displaystyle(K_{D})\bm{\alpha}_{1}=$ $\displaystyle-\sum_{i=1}^{L}\Bigg{(}$
$\displaystyle\mathrm{i}\frac{\partial\widehat{E}_{i}}{\partial
Z}\bm{\alpha}_{0}^{i}+\mathrm{\mathbf{k}}\cdot\nabla\left\langle
M\right\rangle(\mathrm{\mathbf{K}}_{D})\widehat{E}_{i}\bm{\alpha}_{0}^{i}+\left\langle
N\right\rangle(\mathrm{\mathbf{K}}_{D})\widehat{E}_{i}\bm{\alpha}_{0}^{i}\Bigg{)}$
(82)
$\displaystyle-\widetilde{\sigma}\mathbb{N}(\bm{\alpha}_{0})\bm{\alpha}_{0}.$

This represents a maximal balance between not only the slowly varying
envelope, the slow local dispersion of the unperturbed (non-driven) lattice,
and the strength of the nonlinearity, but also the warping of the dispersion
relation due to the driving of the lattice. Hence, we are either finding
solutions that exist in the gap, as is the case for the honeycomb and kagome
lattices, or on the flat band, as is the case for the Lieb and 1/5 depleted
lattices. As mentioned previously, $\left\langle M\right\rangle$ is a
Hermitian matrix and so the Fredholm alternative says the right hand side must
be orthogonal to the $L$ eigenvectors.

In the case of the honeycomb and Lieb lattices, things are further reduced
since the eigenvectors span the whole space at the point of degeneracy. We
assume $\bm{\alpha}_{0}$ can be written as a vector of envelopes,
$\left(\widehat{E}_{1}(Z),\ldots,\widehat{E}_{L}(Z)\right)^{T}$, and after
applying orthogonality the right hand side yields the following system of
nonlinear system of equations for
$\widehat{\mathrm{\mathbf{E}}}=\left(\widehat{E}_{1},\ldots,\widehat{E}_{L}\right)^{T}$

$\mathrm{i}\frac{\partial\widehat{\mathrm{\mathbf{E}}}}{\partial
Z}+\mathbb{L}\left(\mathrm{\mathbf{k}}\right)\widehat{\mathrm{\mathbf{E}}}+\widetilde{\sigma}\mathbb{N}\left(\widehat{\mathrm{\mathbf{E}}}\right)\widehat{\mathrm{\mathbf{E}}}=0$
(83)

where

$\mathbb{L}(\mathrm{\mathbf{k}})=\nabla_{\mathrm{\mathbf{k}}}\left\langle
M\right\rangle(\mathrm{\mathbf{K}}_{D})\cdot{\bm{k}}+\left\langle
N\right\rangle(\mathrm{\mathbf{K}}_{D}).$ (84)

Finally, moving back to physical space, in the neighborhood of $K_{D}$ we
replace $\mathrm{\mathbf{k}}$ with $-\mathrm{i}\nabla_{\mathrm{\mathbf{R}}}$
for the slow spatial variable
$\mathrm{\mathbf{R}}=\varepsilon\mathrm{\mathbf{r}}$.

$\mathrm{i}\frac{\partial\mathrm{\mathbf{E}}}{\partial
Z}+\mathbb{L}(-\mathrm{i}\nabla_{\mathrm{\mathbf{R}}})\mathrm{\mathbf{E}}+\widetilde{\sigma}\mathbb{N}(\mathrm{\mathbf{E}})\mathrm{\mathbf{E}}=0.$
(85)

It should be noted that we have treated the nonlinear terms independently of
the liner terms when transforming between the physical and wave number
domains. This is a natural approach when the dispersive and nonlinear terms do
not interact and results in the same final envelope equations as a more
rigorous accounting of the analysis; see [20] for more details on this type of
argument.

For the honeycomb lattice, the envelope equations are

$\displaystyle\mathrm{i}\partial_{Z}E_{1}-\frac{3}{2}(-\partial_{X}E_{2}$
$\displaystyle+\mathrm{i}\partial_{Y}E_{2})+\nu_{0}E_{1}+\widetilde{\sigma}\left|E_{1}\right|^{2}E_{1}=0$
(86a) $\displaystyle\mathrm{i}\partial_{Z}E_{2}-\frac{3}{2}(\partial_{X}E_{1}$
$\displaystyle+\mathrm{i}\partial_{Y}E_{1})-\nu_{0}E_{2}+\widetilde{\sigma}\left|E_{2}\right|^{2}E_{2}=0$
(86b)

where $\nu_{0}=\nu(K_{D})$ and we assume a perfect honeycomb lattice $\rho=1$.
This is the nonlinear Dirac system studied earlier cf. [18] but now with an
important additional linear term proportional to $\nu_{0}$ that opens a gap in
the dispersion relation sometimes called the Dirac equation with varying mass
[33]. The corresponding dispersion relations for plane waves
$\exp\left[i(KX+LY-\Omega Z)\right]$ are

$\Omega(K,L)=\pm\sqrt{\nu_{0}^{2}+\frac{9}{4}\left(K^{2}+L^{2}\right)}.$ (87)

Note: $\Omega(K,L)$ has a dispersion relation similar to that of two-
dimensional Klein-Gordon equations. As such, there is dispersion and hence
localized initial data will disperse for large $z$ [20]. This spectrum is the
continuous analog of the discrete eigenvalues in equation (2), with a spectral
gap at $(K,L)=(0,0)$ of width $2|\nu_{0}|$.

For the system (86) we consider Gaussian initial envelopes over both carrier
eigenfunctions. The effect of the lattice rotation is to turn the conical
diffraction usually associated with Dirac points [19], into the spiral
diffraction appearing in Figure 8. Here the nonlinear effects are relatively
weak, $\widetilde{\sigma}=0.10$, in comparison to the effects of lattice
rotation, $\nu_{0}=2$. As the strength of the nonlinearity is increased, the
nonlinear Dirac system transitions to a regime with wave collapse. This
suggests the existence of soliton-like solution which persist over a long
distance at some critical strength of the nonlinearity. Figure 9 shows an
example of such a soliton-like solution and we see that the light focus in the
second second component of the system. Figure 10 shows that as nonlinearity
further increases both components collapse.

Figure 8: Evolution of an initial Gaussian pulse in the nonlinear Dirac system
(86) for $\nu_{0}=2$ and $\widetilde{\sigma}=0.1$. The inclusion of the term
$\nu_{0}$ term creates a spiral diffraction pattern. Figure 9: Evolution of
the nonlinear Dirac system (86) for $\nu_{0}=2$ and $\widetilde{\sigma}=0.75$.
A persistent, soliton-like, pulse develops in the second component. Figure 10:
Evolution of an initial Gaussian pulse in the nonlinear Dirac system (86) for
$\nu_{0}=2$ and $\widetilde{\sigma}=1$. For nonlinear above the critical
threshold the pulse will undergo collapse in finite time.

For the Lieb lattice, the envelope equations are

$\displaystyle\mathrm{i}\partial_{Z}E_{1}$
$\displaystyle-\partial_{X}E_{2}-\mathrm{i}\nu_{0}E_{3}+\widetilde{\sigma}\left|E_{1}\right|^{2}E_{1}=0$
(88a) $\displaystyle\mathrm{i}\partial_{Z}E_{2}$
$\displaystyle+\partial_{X}E_{1}-\partial_{Y}E_{3}+\widetilde{\sigma}\left|E_{2}\right|^{2}E_{2}=0$
(88b) $\displaystyle\mathrm{i}\partial_{Z}E_{3}$
$\displaystyle+\partial_{Y}E_{2}+\mathrm{i}\nu_{0}E_{1}+\widetilde{\sigma}\left|E_{3}\right|^{2}E_{3}=0.$
(88c)

They appear to be a novel nonlinear system. The corresponding dispersion
relations are

$\mu(K,L)=0,~{}~{}\pm\sqrt{\nu_{0}^{2}+K^{2}+L^{2}}.$ (89)

with a gap width of $2|\nu_{0}|$ at the origin. Comparing these continuous
eigenvalues to the discrete versions in equation (60), we observe a natural
extension including a zero plane separating two hyperboloid sheets. Moreover,
due to the vanishing of one of the dispersion relation functions localized
initial values will generally not decay for large $z$.

Figure 11: Evolution of an initial pulse in the nonlinear Lieb system (88)
with $\nu_{0}=2$ and $\widetilde{\sigma}=0.6$. Energy oscillates between the
$E_{1}$ and $E_{3}$ modes, mediated by $E_{2}$.

Again we consider initial conditions where a Gaussian envelope is taken on all
(three) eigenfunctions. In Figure 11, we see an example nonlinear evolution
for a pulse in which the energy of the system oscillates between the first and
third components of the system. These oscillations are a feature of the
related linear system and the focusing nonlinearity effect. Similar to the
nonlinear Dirac system, the nonlinear Lieb system transitions to a regime of
wave collapse as the nonlinear coefficient is increased.

##  6 Conclusions

Tight-binding approximations are constructed for a longitudinally driven two-
dimensional lattices with uniform lattice sites. The general lattice is doubly
periodic in the transverse directions and consists of an arbitrary number of
sublattice sites within each unit cell. An orbital expansion results in a
system of ODE’s with order equal to the number of sublattice sites. For fast
oscillations in the longitudinal direction, an averaged equation is derived
with constant coefficients in terms of $z$ and parameterized by the transverse
wavenumber, $\bm{k}$.

These averaged models have the effect of introducing next-nearest neighbor
interactions into the system. For the honeycomb lattice the averaging process
leads to the well known Haldane Model in Fourier space. Analytic calculation
of the Berry connection and curvature means that global information about the
topology of the spectral bands can be obtained. Only a local approximation of
the eigenvectors is required for the connection and curvature since the Chern
number is determined by the local behavior near discontinuities. Detailed
analysis is also carried out for the Lieb lattice.

Nonlinear wave equations are derived that govern the envelopes of the
underlying eigenvectors. In the case of a honeycomb lattice the resulting
nonlinear Dirac system contains an extra ”mass” term, in comparison to
previous models [18, 19]. Numerical computation shows spiral diffraction
patterns for relatively small nonlinear effects. As the effect of the
nonlinearity is increased the wave envelopes tend towards a soliton type
structure and at larger nonlinear coefficent the envelopes undergo collapse.
For Dirac points of the spectrum with a flatband all crossing the point, as
appears in the Lieb and $1/5$-depleted lattices, a novel nonlinear system is
derived for the envelope dynamics. In this system simulations show energy
oscillates between the eigenfunctions.

The authors have no competing interests to declare.

This article has no additional data.

## Acknowledge

This work was partially supported by AFOSR under Grant No. FA9550-19-1-0084
and NSF under Grant No. DMS-2005343.

## References

  * [1] K. V. Klitzing, G. Dorda, and M. Pepper.  New method for high-accuracy determination of the fine-structure constant based on quantized hall resistance.  Phys. Rev. Lett., 45:494–497, Aug 1980.
  * [2] D. C. Tsui, H. L. Stormer, and A. C. Gossard.  Two-dimensional magnetotransport in the extreme quantum limit.  Phys. Rev. Lett., 48:1559–1562, May 1982.
  * [3] F. D. M. Haldane.  Model for a quantum hall effect without landau levels: Condensed-matter realization of the ”parity anomaly”.  Phys. Rev. Lett., 61:2015–2018, Oct 1988.
  * [4] F. D. M. Haldane and S. Raghu.  Possible realization of directional optical waveguides in photonic crystals with broken time-reversal symmetry.  Phys. Rev. Lett., 100:013904, Jan 2008.
  * [5] D. J. Thouless, M. Kohmoto, M. P. Nightingale, and M. den Nijs.  Quantized hall conductance in a two-dimensional periodic potential.  Phys. Rev. Lett., 49:405–408, Aug 1982.
  * [6] J. E. Avron, R. Seiler, and B. Simon.  Homotopy and quantization in condensed matter physics.  Phys. Rev. Lett., 51:51–53, Jul 1983.
  * [7] Barry Simon.  Holonomy, the quantum adiabatic theorem, and berry’s phase.  Phys. Rev. Lett., 51:2167–2170, Dec 1983.
  * [8] Michael Victor Berry.  Quantal phase factors accompanying adiabatic changes.  Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences, 392(1802):45–57, 1984.
  * [9] Mahito Kohmoto.  Topological invariant and the quantization of the hall conductance.  Annals of Physics, 160(2):343–354, 1985.
  * [10] Zheng Wang, Y. D. Chong, John D. Joannopoulos, and Marin Soljačić.  Reflection-free one-way edge modes in a gyromagnetic photonic crystal.  Phys. Rev. Lett., 100:013905, Jan 2008.
  * [11] Mikael C. Rechtsman, Julia M. Zeuner, Yonatan Plotnik, Yaakov Lumer, Daniel Podolsky, Felix Dreisow, Stefan Nolte, Mordechai Segev, and Alexander Szameit.  Photonic floquet topological insulators.  Nature, 496(7444):196–200, 2013.
  * [12] Mark J. Ablowitz, Christopher W. Curtis, and Yi-Ping Ma.  Linear and nonlinear traveling edge waves in optical honeycomb lattices.  Phys. Rev. A, 90:023813, Aug 2014.
  * [13] Ling Lu, John D. Joannopoulos, and Marin Soljačić.  Topological photonics.  Nature Photonics, 8(11):821–829, 2014.
  * [14] Tomoki Ozawa, Hannah M. Price, Alberto Amo, Nathan Goldman, Mohammad Hafezi, Ling Lu, Mikael C. Rechtsman, David Schuster, Jonathan Simon, Oded Zilberberg, and Iacopo Carusotto.  Topological photonics.  Rev. Mod. Phys., 91:015006, Mar 2019.
  * [15] Mark J. Ablowitz and Justin T. Cole.  Tight-binding methods for general longitudinally driven photonic lattices: Edge states and solitons.  Phys. Rev. A, 96:043868, Oct 2017.
  * [16] Mark J. Ablowitz and Justin T. Cole.  Topological insulators in longitudinally driven waveguides: Lieb and kagome lattices.  Phys. Rev. A, 99:033821, Mar 2019.
  * [17] Mark J. Ablowitz and Justin T. Cole.  Discrete approximation of topologically protected modes in magneto-optical media.  Phys. Rev. A, 101:023811, Feb 2020.
  * [18] Mark J. Ablowitz, Sean D. Nixon, and Yi Zhu.  Conical diffraction in honeycomb lattices.  Phys. Rev. A, 79:053830, May 2009.
  * [19] Mark J. Ablowitz and Yi Zhu.  Evolution of bloch-mode envelopes in two-dimensional generalized honeycomb lattices.  Phys. Rev. A, 82:013840, Jul 2010.
  * [20] M.J. Ablowitz.  Nonlinear Dispersive Waves: Asymptotic Analysis and Solitons.  Cambridge Texts in Applied Mathematics. Cambridge University Press, 2011\.
  * [21] Alexander Szameit and Stefan Nolte.  Discrete optics in femtosecond-laser-written photonic structures.  Journal of Physics B: Atomic, Molecular and Optical Physics, 43(16):163001, jul 2010.
  * [22] Amir Sagiv and Michael I. Weinstein.  Effective gaps in continuous floquet hamiltonians.  SIAM Journal on Mathematical Analysis, 54(1):986–1021, 2022.
  * [23] S S Chern, W H Chen, and K S Lam.  Lectures on Differential Geometry.  World Scientific, 1999.
  * [24] Takahiro Fukui, Yasuhiro Hatsugai, and Hiroshi Suzuki.  Chern numbers in discretized brillouin zone: Efficient method of computing (spin) hall conductances.  Journal of the Physical Society of Japan, 74(6):1674–1677, 2005\.
  * [25] Matthew J. Allen, Vincent C. Tung, and Richard B. Kaner.  Honeycomb carbon: A review of graphene.  Chemical Reviews, 110(1):132–145, 2010.
  * [26] Charles L. Fefferman and Michael I. Weinstein.  Honeycomb lattice potentials and dirac points.  J. Amer. Math. Soc., 25:1169–1220, 2012.
  * [27] Giovanna Marcelli, Domenico Monaco, Massimo Moscolari, and Gianluca Panati.  The haldane model and its localization dichotomy.  Rendiconti di Mat. e delle Sue Appl., 39:307–327, 2018.
  * [28] Elliott H. Lieb.  Two theorems on the hubbard model.  Phys. Rev. Lett., 62:1201–1204, Mar 1989.
  * [29] C. Weeks and M. Franz.  Topological insulators on the lieb and perovskite lattices.  Phys. Rev. B, 82:085310, Aug 2010.
  * [30] H.-M. Guo and M. Franz.  Topological insulator on the kagome lattice.  Phys. Rev. B, 80:113102, Sep 2009.
  * [31] Takahiro Fukui, Yasuhiro Hatsugai, and Hiroshi Suzuki.  Chern numbers in discretized brillouin zone: Efficient method of computing (spin) hall conductances.  Journal of the Physical Society of Japan, 74(6):1674–1677, 2005\.
  * [32] L. O. Manuel, M. I. Micheletti, A. E. Trumper, and H. A. Ceccatto.  Heisenberg model on the $\frac{1}{5}$-depleted square lattice and the ${\mathrm{cav}}_{4}{\mathrm{o}}_{9}$ compound.  Phys. Rev. B, 58:8490–8494, Oct 1998.
  * [33] Pipi Hu, Peng Xie, and Yi Zhu.  Traveling edge states in massive dirac equations along slowly varying edges.  2022\.

##  Appendix A Tight binding details

A heuristic understanding of the $\mathcal{L}$ operators in (9) can be reached
by considering a lattice with a single pair of lattice sites centered at
$\mathrm{\mathbf{r}}_{1}$ and $\mathrm{\mathbf{r}}_{2}$, i.e.,
$V(\mathrm{\mathbf{r}})=V_{0}^{2}-\widetilde{V}(\mathrm{\mathbf{r}}-\mathrm{\mathbf{r}}_{1})-\widetilde{V}(\mathrm{\mathbf{r}}-\mathrm{\mathbf{r}}_{2})$.
This corresponds to an approximation with scalar coefficient functions
$a_{1}(z)$ and $a_{2}(z)$,

$\psi(\mathrm{\mathbf{r}},z)\approx\Big{(}a_{1}(z)\phi_{1}(\mathrm{\mathbf{r}})+a_{2}(z)\phi_{2}(\mathrm{\mathbf{r}})\Big{)}\mathrm{e}^{-\mathrm{i}\left(\mu+V_{0}^{2}\right)z},$
(90)

where
$\phi_{j}(\mathrm{\mathbf{r}})=\phi_{0}(\mathrm{\mathbf{r}}-\mathrm{\mathbf{r}}_{j})$.
Substituting this directly into equation (6) gives

$\displaystyle 0=$
$\displaystyle~{}~{}\mathrm{i}a_{1}^{\prime}(z)\phi_{1}-\mathrm{i}a_{1}(z)\nabla\phi_{1}\cdot\mathrm{\mathbf{h}}^{\prime}(z)+a_{1}(z)\widetilde{V}_{2}\phi_{1}$
(91)
$\displaystyle+\mathrm{i}a_{2}^{\prime}(z)\phi_{1}-\mathrm{i}a_{2}(z)\nabla\phi_{2}\cdot\mathrm{\mathbf{h}}^{\prime}(z)+a_{2}(z)\widetilde{V}_{1}\phi_{2}$
$\displaystyle+\sigma\Big{(}\lvert a_{1}\rvert^{2}a_{1}\phi_{1}^{3}+2\lvert
a_{1}\rvert^{2}a_{2}\phi_{1}^{2}\phi_{2}+a_{1}^{2}a_{2}^{*}\phi_{1}^{2}\phi_{2}$
$\displaystyle~{}~{}~{}~{}~{}~{}+\lvert
a_{1}\rvert^{2}a_{2}\phi_{2}^{3}+2\lvert
a_{2}\rvert^{2}a_{1}\phi_{2}^{2}\phi_{1}+a_{2}^{2}a_{1}^{*}\phi_{2}^{2}\phi_{1}\Big{)}$

where
$\widetilde{V}_{j}(\mathrm{\mathbf{r}})=\widetilde{V}(\mathrm{\mathbf{r}}-\mathrm{\mathbf{r}}_{j})$.
Taking the inner product of (91) with $\phi_{j}$, we keep terms that scale
with $V_{0}$ and are large due to the deep lattice assumption, $V_{0}\gg 1$.
This includes terms with $\widetilde{V}$ and terms with $\nabla\phi_{j}$ that
scale with $V_{0}$ due to the rapid decay rate of the tails. Otherwise, tail
interactions are exponentially small and we drop all but the onsite nonlinear
effects. This results in the system

$\displaystyle\mathrm{i}\frac{\mathrm{d}a_{1}}{\mathrm{d}z}$
$\displaystyle+p_{0}a_{1}+\Big{(}q_{0}-\mathrm{i}s_{0}\mathrm{\mathbf{h}}^{\prime}(z)\cdot(\mathrm{\mathbf{r}}_{2}-\mathrm{\mathbf{r}}_{1})\Big{)}a_{2}+\widetilde{\sigma}\left|a_{1}\right|^{2}a_{1}=0$
(92a) $\displaystyle\mathrm{i}\frac{\mathrm{d}a_{2}}{\mathrm{d}z}$
$\displaystyle+p_{0}a_{2}+\Big{(}q_{0}-\mathrm{i}s_{0}\mathrm{\mathbf{h}}^{\prime}(z)\cdot(\mathrm{\mathbf{r}}_{1}-\mathrm{\mathbf{r}}_{2})\Big{)}a_{1}+\widetilde{\sigma}\left|a_{2}\right|^{2}a_{2}=0$
(92b)

where all coefficients are determined by the orbital form $\phi_{0}(\bm{r})$
and distance between the nearest neighbor lattice sites.

$\displaystyle
p_{0}=\iint\widetilde{V}(\mathrm{\mathbf{r}}-\Delta\mathrm{\mathbf{r}})\phi_{0}(\mathrm{\mathbf{r}})^{2}\mathrm{d}r_{x}\mathrm{d}r_{y},\quad
q_{0}=\iint\widetilde{V}(\mathrm{\mathbf{r}})\phi_{0}(\mathrm{\mathbf{r}}-\Delta\mathrm{\mathbf{r}})\phi_{0}(\mathrm{\mathbf{r}})\mathrm{d}r_{x}\mathrm{d}r_{y}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}s_{0}\Delta\mathrm{\mathbf{r}}\approx\iint\nabla\phi_{0}(\mathrm{\mathbf{r}}-\Delta\mathrm{\mathbf{r}})\phi_{0}(\mathrm{\mathbf{r}})\mathrm{d}r_{x}\mathrm{d}r_{y}$

In the case of the last integral, this is an asymptotic approximation that has
been verified against numerical results. This suggests an integrating factor
$a_{j}=\widetilde{a}_{j}e^{\mathrm{i}\,p_{0}z}$, to remove onsite interactions
from the system. In term of our binary lattice example this gives

$\displaystyle\mathrm{i}\frac{\mathrm{d}\widetilde{a}_{1}}{\mathrm{d}z}+\Big{(}q_{0}-\mathrm{i}s_{0}\mathrm{\mathbf{h}}^{\prime}(z)\cdot(\mathrm{\mathbf{r}}_{2}-\mathrm{\mathbf{r}}_{1})\Big{)}\widetilde{a}_{2}+\widetilde{\sigma}\left|\widetilde{a}_{1}\right|^{2}\widetilde{a}_{1}$
$\displaystyle=0$ (93a)
$\displaystyle\mathrm{i}\frac{\mathrm{d}\widetilde{a}_{2}}{\mathrm{d}z}+\Big{(}q_{0}-\mathrm{i}s_{0}\mathrm{\mathbf{h}}^{\prime}(z)\cdot(\mathrm{\mathbf{r}}_{1}-\mathrm{\mathbf{r}}_{2})\Big{)}\widetilde{a}_{1}+\widetilde{\sigma}\left|\widetilde{a}_{2}\right|^{2}\widetilde{a}_{2}$
$\displaystyle=0.$ (93b)

We now see that for $\widetilde{a}_{1}$ the interaction depends on the vector
$\mathrm{\mathbf{r}}_{2}-\mathrm{\mathbf{r}}_{1}$, which points to the
neighboring lattice site; and vice versa for $\widetilde{a}_{2}$. We define
the coefficient function given earlier in equation (15),

$q(z,\mathrm{\mathbf{w}})=q_{0}-\mathrm{i}s_{0}\mathrm{\mathbf{h}}^{\prime}(z)\cdot\mathrm{\mathbf{w}},$
(94)

and the system of equations for a full lattice of the form (3) can now be seen
as a superposition of all pairwise interactions between nearest neighbor
lattice sites.