# Mirror Chern numbers in the hybrid Wannier representation

Tomáš Rauch  Friedrich-Schiller-University Jena, 07743 Jena, Germany  Thomas
Olsen  Computational Atomic-scale Materials Design, Department of Physics,
Technical University of Denmark, 2800 Kgs. Lyngby Denmark  David Vanderbilt
Department of Physics and Astronomy, Rutgers University, Piscataway, New
Jersey 08854-8019, USA  Ivo Souza  Centro de Física de Materiales, Universidad
del País Vasco, 20018 San Sebastián, Spain  Ikerbasque Foundation, 48013
Bilbao, Spain

###### Abstract

The topology of electronic states in band insulators with mirror symmetry can
be classified in two different ways. One is in terms of the mirror Chern
number, an integer that counts the number of protected Dirac cones in the
Brillouin zone of high-symmetry surfaces. The other is via a $\mathbbm{Z}_{2}$
index that distinguishes between systems that have a nonzero quantized orbital
magnetoelectric coupling (“axion-odd”), and those that do not (“axion-even”);
this classification can also be induced by other symmetries in the magnetic
point group, including time reversal and inversion. A systematic
characterization of the axion $\mathbbm{Z}_{2}$ topology has previously been
obtained by representing the valence states in terms of hybrid Wannier
functions localized along one chosen crystallographic direction, and
inspecting the associated Wannier band structure. Here we focus on mirror
symmetry, and extend that characterization to the mirror Chern number. We
choose the direction orthogonal to the mirror plane as the Wannierization
direction, and show that the mirror Chern number can be determined from the
winding numbers of the touching points between Wannier bands on mirror-
invariant planes, and from the Chern numbers of flat bands pinned to those
planes. In this representation, the relation between the mirror Chern number
and the axion $\mathbbm{Z}_{2}$ index is readily established. The formalism is
illustrated by means of ab initio calculations for SnTe in the monolayer and
bulk forms, complemented by tight-binding calculations for a toy model.

##  I Introduction

The band theory of solids has been enriched in recent years by a vigorous
study of its topological aspects. That effort resulted in a systematic
topological classification of insulators on the basis of symmetry, and in the
identification of a large number of topological materials. After an initial
focus on the role of time-reversal symmetry, it was realized that
crystallographic symmetries could also protect topological behaviors, leading
to the notion of “topological crystalline insulators.”

The assignment of an insulator to a particular topological class can be made
by evaluating the corresponding topological invariant. Depending on the
protecting symmetry, that invariant may assume one of two values
($\mathbbm{Z}_{2}$ classification), or it may assume any integer value
($\mathbbm{Z}$ classification). Other types of classifications such as
$\mathbbm{Z}_{4}$ also occur, but they do not concern us here. When the
invariant vanishes the system is classified as trivial, and otherwise it is
classified as nontrivial or topological. Topological insulators typically
display robust gapless states at the boundary, which provide an experimental
signature of topological behavior.

In some cases, the same symmetry may induce two different topological
classifications. This happens for example with mirror symmetry, where a
$\mathbbm{Z}$ classification in terms of the mirror Chern number (MCN) Teo
_et al._ (2008); Ando and Fu (2015) coexists with a $\mathbbm{Z}_{2}$
classification based on the quantized axion angle. The two classifications are
not independent, and elucidating the relation between them is one goal of the
present work.

The axion $\mathbbm{Z}_{2}$ classification of three-dimensional (3D)
insulators is based on the orbital magnetoelectric effect. In brief, the
isotropic part of the linear orbital magnetoelectric tensor is conveniently
expressed in terms of the axion angle $\theta$, which is only defined modulo
$2\pi$ as a bulk property. In the presence of “axion-odd” symmetries that flip
its sign, the axion angle can only assume two values: $\theta=0$ (trivial),
and $\theta=\pi$ (topological) Qi  _et al._ (2008); Essin  _et al._ (2009);
Vanderbilt (2018); Armitage and Wu (2019); Nenno  _et al._ (2020); Sekine and
Nomura (2021).

The axion $\mathbbm{Z}_{2}$ index was originally introduced for time-reversal
invariant insulators, where it was shown to be equivalent to the “strong”
$\mathbbm{Z}_{2}$ index $\nu_{0}=0$ or $1$, that is, $\theta=\pi\nu_{0}$. More
generally, axion-odd symmetries can be classified as proper rotations combined
with time reversal (including time reversal itself), and improper rotations
(including inversion and reflection) not combined with time reversal; in both
cases, the associated symmetry operation in the magnetic space group may
include a fractional translation. This results in a large number of magnetic
space groups that can host axion-odd topological insulators. A recent
realization is the MnBi2Te4 family of antiferromagnetic materials Otrokov  _et
al._ (2019); Nenno  _et al._ (2020); Sekine and Nomura (2021), whose axion
topology is protected by the time reversal operation combined with a half-
lattice translation as envisioned in Ref. Mong  _et al._ (2010).

To aid the computational search for axionic topological insulators, it is
useful to devise simple procedures for determining the (quantized) axion angle
$\theta$. Unfortunately, subtle gauge issues make its direct evaluation from
the valence Bloch states rather challenging in general Vanderbilt (2018).
Notable exceptions are centrosymmetric insulators, both nonmagnetic and
magnetic. For such systems, the axion $\mathbbm{Z}_{2}$ index can be obtained
by counting the number of odd-parity states at high-symmetry points in the
Brillouin zone (BZ) Fu and Kane (2007); Turner  _et al._ (2012).

Recently, an alternative procedure was introduced based on representing the
valence states in terms of hybrid Wannier (HW) functions that are maximally
localized along a chosen crystallographic direction $z$. The HW centers along
$z$, also known as “Wilson-loop eigenvalues,” form a band structure when
plotted as a function of $k_{x}$ and $k_{y}$; in the presence of one or more
axion-odd symmetries, the quantized $\theta$ value can be determined from this
“Wannier band structure,” often by mere visual inspection Varnava  _et al._
(2020).

In the HW representation, axion-odd symmetries are naturally classified as
“$z$-preserving” or “$z$-reversing,” and the rules for deducing the axion
$\mathbbm{Z}_{2}$ index are different in each case (they also depend on
whether or not the symmetry operation involves a fractional translation along
$z$) Varnava  _et al._ (2020). Time reversal is an example of a $z$-preserving
operation, while inversion is $z$ reversing. Mirror operations may be placed
in one group or the other, depending on whether the Wannierization direction
$z$ lies in the reflection plane (vertical mirror) or is orthogonal to it
(horizontal mirror). In this work we make the latter choice, so that the
mirror operation of interest becomes

$M_{z}:z\rightarrow-z\,,$ (1)

which is manifestly $z$ reversing.

A simple mirror symmetry without a glide component protects not only the axion
$\mathbbm{Z}_{2}$ classification, but also a $\mathbbm{Z}$ or
$\mathbbm{Z}\times\mathbbm{Z}$ classification based on one or two MCNs,
depending on the type of mirror. This raises the question of whether the HW
representation might also be convenient for determining the MCNs, and for
illuminating their relationship to the quantized axion angle.

In this work, we address the above questions by investigating in detail the
Wannier bands in the presence of $M_{z}$ symmetry. We clarify the generic
behaviors that are expected, and discuss the rules for deducing the MCNs. By
comparing those rules with the ones obtained in Ref. Varnava  _et al._ (2020)
for the axion $\mathbbm{Z}_{2}$ index, we establish the relation between the
two classifications.

The paper is organized as follows. In Sec. II we first distinguish between
“type-1” and “type-2” crystallographic mirror operations; we then review the
definitions of Chern invariants and MCNs in terms of the Bloch states in the
filled bands; finally, we introduce maximally localized HW functions spanning
the valence states, and assign Chern numbers to isolated groups of Wannier
bands. This background material sets the stage for the developments in the
remainder of the paper. In Sec. III we discuss the generic features of the
Wannier band structure in the presence of $M_{z}$ symmetry, and obtain a
relation between Chern numbers and winding numbers in groups of bands touching
on a mirror plane. The rules for deducing the MCNs from the Chern numbers and
winding numbers on the mirror planes are given in Sec. IV, where their
relation to the quantized axion angle is also established. In Sec. V we
describe the numerical methods that are used in Sec. VI to apply the formalism
to several prototypical systems. We summarize and conclude in Sec. VII, and
present in three Appendices some derivations that were left out of the main
text.

##  II Preliminaries

###  II.1 Two types of crystallographic mirrors

Figure 1: The upper panel shows schematically a pair of 2D crystals lying on
the $(x,z$) plane; each has one atom per primitive cell (black dots), and
lattice constant $c$ along $z$. The crystal on the left has a rectangular
lattice and a type-1 horizontal mirror, with inequivalent mirror lines
$z=0\text{ mod $c$}$ (A) and $z=c/2\text{ mod $c$}$ (B), shown as dashed
lines; the one on the right has a centered rectangular lattice and a type-2
mirror, with equivalent mirror lines A and B. The lattice vectors ${\bf
a}_{3}$ and $\widetilde{{\bf a}}_{3}$ are defined in the main text. The lower
panel shows the reciprocal lattices, with a separation of $2\pi/c$ between
horizontal lattice lines G. On the left the periodicity along $k_{z}$ is
$2\pi/c$, and hence both $k_{z}=0\text{ mod $2\pi/c$}$ (G) and
$k_{z}=\pi/c\text{ mod $2\pi/c$}$ (X) are pointwise-invariant mirror lines, as
indicated by the dashed lines. On the right, where the periodicity along
$k_{z}$ is $4\pi/c$, G is a mirror-invariant line but X is not. The associated
Brillouin zones are indicated by the shaded green areas.

We begin by observing that if a crystal is left invariant under an $M_{z}$
reflection operation, then its Bravais lattice must contain vectors pointing
along $z$. To construct the shortest such vector ${\bf a}_{3}=c\hat{\bf z}$,
we pick the shortest vector $\widetilde{{\bf a}}_{3}$ connecting lattice
points on adjacent horizontal lattice planes. If $\widetilde{{\bf a}}_{3}$
points along $z$ then we take it as ${\bf a}_{3}$, and we say that the mirror
is of type 1. Otherwise we choose the vector ${\bf a}_{3}=\widetilde{{\bf
a}}_{3}-M_{z}\widetilde{{\bf a}}_{3}$ connecting second-neighbor lattice
planes, and the mirror is of type 2.

The two types of crystallographic mirrors are exemplified in 2D in Fig. 1,
where the mirror lines $z=0$ and $c/2$ are labeled A and B, and the
reciprocal-space lines $k_{z}=0$ and $k_{z}=\pi/c$ are labeled G and X. The
same notation will be used in 3D, where A and B (G and X) become planes in
real (reciprocal) space.

The distinction between mirror operations that leave pointwise invariant two
inequivalent planes in the BZ, and those that leave invariant only one BZ
plane, was made in Refs. Varjas  _et al._ (2015); Fulga  _et al._ (2016).
Since MCNs are defined on such planes Teo  _et al._ (2008); Ando and Fu
(2015), a 3D insulator with a type-1 mirror is characterized by two separate
MCNs $\mu_{\rm G}$ and $\mu_{\rm X}$, while for a type-2 mirror there is a
single MCN $\mu_{\rm G}$. If the crystallographic space group contains
additional mirror operations, there will be additional MCNs associated with
them.

###  II.2 Chern invariants in band insulators

####  II.2.1 Generic insulators

Before introducing MCNs for insulators with reflection symmetry, let us define
Chern invariants for generic 2D and 3D band insulators in terms of the ${\bf
k}$-space Berry curvature of the valence states Vanderbilt (2018).

In 2D, the Berry curvature of a Bloch state $|\psi_{n{\bf k}}\rangle$ with
cell-periodic part $|u_{n{\bf k}}\rangle$ is a scalar defined as

$\Omega_{n{\bf k}}=-2{\rm Im\,}\langle\partial_{k_{x}}u_{n{\bf
k}}|\partial_{k_{y}}u_{n{\bf k}}\rangle$ (2)

where ${\bf k}=(k_{x},k_{y})$, and the Chern number is given by

$C=\frac{1}{2\pi}\int_{\rm 2DBZ}\sum_{n=1}^{J}\,\Omega_{n{\bf k}}\,d^{2}k$ (3)

where the summation is over the $J$ filled energy bands. Since the Berry
curvature has units of length squared, $C$ is a dimensionless number, and for
topological reasons it must be an integer. The Chern number is a global
property of the manifold of occupied states, remaining invariant under
multiband gauge transformations described by $J\times J$ unitary matrices at
each ${\bf k}$, and it vanishes when the crystal has time-reversal symmetry.
If a 2D magnetic crystal has a nonzero Chern number $C$, when that crystal is
terminated at an edge there will be $|C|$ edge modes crossing the bulk gap,
whose chirality will depend on the sign of $C$.

3D insulators are characterized by a Chern vector

${\bf K}=\frac{1}{2\pi}\int_{\rm 3DBZ}\sum_{n=1}^{J}\,{\bm{\Omega}}_{n{\bf
k}}\,d^{3}k\,,$ (4)

where now ${\bf k}=(k_{x},k_{y},k_{z})$ and the Berry curvature has become a
vector field, ${\bm{\Omega}}_{n{\bf k}}=-{\rm Im\,}\langle\partial_{\bf
k}u_{n{\bf k}}|\times|\partial_{\bf k}u_{n{\bf k}}\rangle$. The Chern vector
has units of inverse length, and is quantized to be a reciprocal-lattice
vector. Like the Chern number in 2D, the Chern vector always vanishes in
nonmagnetic crystals.

Given a set of lattice vectors ${\bf a}_{j}$ and dual reciprocal-lattice
vectors ${\bf b}_{j}$, the expansion ${\bf K}=\sum_{j}\,C_{j}{\bf b}_{j}$
defines a triad of integer Chern indices $C_{j}$. Let us orient the Cartesian
axes such that ${\bf a}_{3}=c\hat{\bf z}$. The vectors ${\bf b}_{1}$ and ${\bf
b}_{2}$ then lie on the $(x,y)$ plane, and the third Chern index can be
expressed as

$C_{3}=\frac{c}{2\pi}\int_{0}^{2\pi/c}\,C(k_{z})\,dk_{z}\,,$ (5)

where

$C(k_{z})=\frac{1}{2\pi}\int_{\rm
2DBZ}\sum_{n=1}^{J}\,\Omega^{z}_{n}(k_{x},k_{y},k_{z})\,dk_{x}dk_{y}\,.$ (6)

The integral in Eq. (6) is over a slice of the 3D BZ spanned by ${\bf b}_{1}$
and ${\bf b}_{2}$ at fixed $k_{z}$. By viewing it as an effective 2D BZ and
comparing with Eq. (3), it becomes clear that $C(k_{z})$ is a Chern number;
and since in a gapped system its integer value cannot change with the
continuous parameter $k_{z}$, Eq. (5) reduces to $C_{3}=C(k_{z})$ evaluated at
any $k_{z}$. The Chern indices of 3D insulators can therefore be evaluated as
Chern numbers defined over individual BZ slices.

####  II.2.2 Mirror-symmetric insulators

We now consider a 3D crystalline insulator with mirror symmetry $M_{z}$, and
assume that its Chern vector ${\bf K}$ vanishes. A new integer-valued
topological index, the MCN, can be defined for such a system as follows Teo
_et al._ (2008); Ando and Fu (2015).

On the mirror-invariant BZ planes, G and possibly X, the energy eigenstates
are also eigenstates of $M_{z}$. The eigenvalues are $i^{F}p$, where $p=\pm 1$
is the “mirror parity” and $F=0$ or $1$ when the electrons are treated as
spinless or spinful particles, respectively. The occupied Bloch states on
those planes can therefore be grouped into “even” ($p=+1$) and “odd” ($p=-1$)
sectors under reflection about the A plane $z=0$, each carrying its own Chern
number. The Chern numbers of the two sectors on the G plane $k_{z}=0$ are
given by

$C_{\rm G}^{\pm}=\frac{1}{2\pi}\int_{\rm 2DBZ}\sum_{n=1}^{J}\,f_{n{\bf
k}}^{\pm}\Omega^{z}_{n}(k_{x},k_{y},k_{z}=0)\,dk_{x}dk_{y}\,,$ (7)

where $f_{n{\bf k}}^{+}=1-f_{n{\bf k}}^{-}$ equals one or zero for a state
with $p=\pm 1$, respectively. The MCN is defined as

$\mu_{\rm G}=\frac{1}{2}\left(C_{\rm G}^{+}-C_{\rm G}^{-}\right)\,,$ (8)

and it is guaranteed to be an integer since $C_{\rm G}^{+}+C_{\rm
G}^{-}=C_{3}$ vanishes by assumption. If the mirror is of type 1, the plane X
carries a second MCN

$\mu_{\rm X}=\frac{1}{2}\left(C_{\rm X}^{+}-C_{\rm X}^{-}\right)\,,$ (9)

where $C_{\rm X}^{\pm}$ is obtained by replacing $k_{z}=0$ with $k_{z}=\pi/c$
in Eq. (7). The MCNs remain invariant under multiband gauge transformations
that do not mix the two mirror-parity sectors. When they are nonzero,
protected gapless modes appear on surfaces that retain the mirror symmetry
$M_{z}$, with $|\mu_{\rm G}|$ and $|\mu_{\rm X}|$ counting the number of Dirac
cones on the two $M_{z}$-invariant lines in the surface BZ Hsieh  _et al._
(2012).

In the case of a 2D or quasi-2D insulator with reflection symmetry $M_{z}$
about its own plane, the entire 2D BZ is left invariant under $M_{z}$. Such a
system has a unique MCN

$\mu_{\rm 2D}=\frac{1}{2}\left(C_{+}-C_{-}\right)\,,$ (10)

where $C_{+}$ and $C_{-}$ are obtained by inserting the 2D Berry curvature of
Eq. (2) in Eq. (7). When the net Chern number $C=C_{+}+C_{-}$ vanishes,
$|\mu_{\rm 2D}|$ becomes an integer that counts the number of pairs of
counterpropagating chiral edge modes Liu  _et al._ (2014).

We note in passing that spin-orbit coupling is required to obtain non-
vanishing MCNs in systems that are either non-magnetic or whose magnetic order
is collinear.

###  II.3 The hybrid Wannier representation

####  II.3.1 Hybrid Wannier functions and Wannier bands

HW functions are obtained from the valence Bloch states of a 2D or 3D
crystalline insulator by carrying out the Wannier construction along a chosen
reciprocal-lattice direction. They are therefore localized along one direction
only, in contrast to ordinary Wannier functions which are localized in all
spatial directions.

Let us momentarily return to a generic 3D insulating crystal, not necessarily
mirror-symmetric. We denote by $z$ the chosen localization direction and let
${\bm{\kappa}}=(k_{x},k_{y})$, so that the wavevector in the 3D BZ becomes
${\bf k}=({\bm{\kappa}},k_{z}$). Given a gauge for the Bloch states that is
periodic in $k_{z}$,
$|\psi_{n{\bm{\kappa}},k_{z}+2\pi/c}\rangle=|\psi_{n{\bm{\kappa}}k_{z}}\rangle$,
the corresponding HW functions are defined as

$|h_{ln{\bm{\kappa}}}\rangle=\frac{1}{2\pi}\int_{-\pi/c}^{\pi/c}\,e^{-ik_{z}lc}e^{-i{\bm{\kappa}}\cdot{\bf
r}}|\psi_{n{\bm{\kappa}}k_{z}}\rangle\,dk_{z},$ (11)

where the index $l$ runs over unit cells along $z$, and $n$ runs over the $J$
HW functions in one unit cell. By factoring out $e^{-i{\bm{\kappa}}\cdot{\bf
r}}$, we have made the HW functions cell periodic in the in-plane directions,
$h_{ln{\bm{\kappa}}}({\bf r}+{\bf R})=h_{ln{\bm{\kappa}}}({\bf r})$ for any
in-plane lattice vector ${\bf R}$. This will be convenient later on when we
define Berry curvatures and Chern numbers in the HW representation.

For each ${\bm{\kappa}}$ in the projected 2D BZ, we choose the multiband gauge
for the Bloch states in such a way that the HW functions have the smallest
possible quadratic spread along $z$. Such maximally-localized HW functions
satisfy the eigenvalue equation Marzari and Vanderbilt (1997)

$P_{\bm{\kappa}}zP_{\bm{\kappa}}|h_{ln{\bm{\kappa}}}\rangle=z_{ln{\bm{\kappa}}}|h_{ln{\bm{\kappa}}}\rangle,$
(12)

where $P_{\bm{\kappa}}$ is the projection operator onto the space of valence
states with in-plane wave vector ${\bm{\kappa}}$. The eigenvalues in Eq. (12)
are the HW centers

$z_{ln{\bm{\kappa}}}=\langle
h_{ln{\bm{\kappa}}}|z|h_{ln{\bm{\kappa}}}\rangle\,,$ (13)

which form Wannier bands. These are periodic in real space along $z$, as well
as in the in-plane reciprocal space,

$z_{ln{\bm{\kappa}}}=z_{0n{\bm{\kappa}}}+lc\,,\quad z_{ln,{\bm{\kappa}}+{\bf
G}}=z_{ln{\bm{\kappa}}}\,,$ (14)

where ${\bf G}$ is an in-plane reciprocal lattice vector.

A Wannier band structure is said to be gapped if it contains at least one
Wannier band per vertical cell that is separated from the band below by a
finite gap at all ${\bm{\kappa}}$. When that is the case, we choose the cell
contents in such a way that the first band, $n=1$, has a gap below it.

####  II.3.2 Chern numbers of Wannier bands

The Berry curvature of a HW state is defined as

$\Omega_{ln}=-2\,{\rm
Im\,}\langle\partial_{k_{x}}h_{ln}|\partial_{k_{y}}h_{ln}\rangle\,,$ (15)

and periodicity along $z$ implies that $\Omega_{ln}=\Omega_{0n}$. (Here and in
the following, we will frequently drop the index ${\bm{\kappa}}$.) When the
Wannier spectrum is gapped, it becomes possible to associate a Chern number
with each isolated group $a$ of bands within a vertical cell,

$C_{la}=\frac{1}{2\pi}\int_{\rm 2DBZ}\sum_{n\in
a}\,\Omega_{ln}\,d^{2}k=C_{0a}\,.$ (16)

From the HW states in a given group, one can construct Bloch-like states at
any ${\bf k}=(k_{x},k_{y},k_{z})$ by inverting Eq. (11). In general these are
not energy eigenstates, and their band indices label Wannier bands rather than
energy bands. Their Berry curvatures along $z$ are given by

$\Omega^{z}_{n}(k_{x},k_{y},k_{z})=\sum_{l}\,e^{ik_{z}lc}\Omega_{0n,ln}(k_{x},k_{y})\,,$
(17)

where

$\Omega_{0n,ln}=i\langle\partial_{k_{x}}h_{0n}|\partial_{k_{y}}h_{ln}\rangle-i\langle\partial_{k_{y}}h_{0n}|\partial_{k_{x}}h_{ln}\rangle$
(18)

is a matrix generalization of Eq. (15) Taherinejad and Vanderbilt (2015). To
evaluate the net Chern number $C_{a}(k_{z})$ of that group of Bloch-like
states on a slice of the 3D BZ, we insert Eq. (17) in Eq. (6) and restrict the
summation over $n$ to $n\in a$. The contributions from the $l\not=0$ terms
drop out,111The expression for $C_{a}(k_{z})$ involves
$\int_{0}^{2\pi/a}\partial_{k_{x}}Y_{0n,ln}(k_{x})\,dk_{x}$ where
$Y_{0n,ln}(k_{x})=\int_{0}^{2\pi/b}A^{y}_{0n,ln}(k_{x},k_{y})\,dk_{y}$, and
another similar integral
$\int_{0}^{2\pi/b}\partial_{k_{y}}X_{0n,ln}(k_{y})\,dk_{y}$. When $l\not=0$
the quantity $Y_{0n,ln}(k_{x})$ becomes fully invariant under band-diagonal
gauge transformations of the HW states. Hence its value at $k_{x}=2\pi/a$ must
be the same as at $k_{x}=0$, and the integral vanishes. yielding

$C_{a}(k_{z})=C_{0a}\,.$ (19)

Hence the Chern numbers are the same in the Bloch-like and HW representations,
as expected since the two representations are related by a unitary
transformation. When the group $a$ comprises all $J$ Wannier bands in one
vertical cell, its Chern number becomes equal to the Chern index $C_{3}$ of
Eq. (5), which vanishes by assumption.

##  III Mirror-symmetric Wannier bands

With the above background material in hand, we now return to our system of
interest – a 3D insulator with $M_{z}$ symmetry – and construct HW functions
localized along the direction $z$ orthogonal to the mirror plane. We begin
this section by discussing the generic features of Wannier band structures
with $M_{z}$ symmetry.

###  III.1 Flat vs dispersive bands, and the uniform parity assumption

If $M_{z}$ is a symmetry of the system, the operator $PzP$ anticommutes with
$M_{z}$. It follows that if a HW function $|h_{ln}\rangle$ satisfies Eq. (12)
with eigenvalue $z_{ln}$, $M_{z}|h_{ln}\rangle$ satisfies it with eigenvalue
$-z_{ln}$. Since $z_{ln}$ is only defined modulo $c$, two situations may
occur. (i) $|h_{ln}\rangle$ and $M_{z}|h_{ln}\rangle$ are orthogonal, in which
case a pair of dispersive bands appear at $\pm z_{ln}$. (ii) $|h_{ln}\rangle$
and $M_{z}|h_{ln}\rangle$ are the same up to a phase, in which case
$|h_{ln}\rangle$ is an eigenstate of $M_{z}$, and a single flat band appears
at either $z=0$ (A plane) or $z=c/2$ (B plane). The Wannier bands of the
system can therefore be classified into flat bands of even or odd mirror
parity at A; flat bands of even or odd mirror parity at B; and dispersive
pairs appearing at $\pm z$.

If there are several flat bands on a given mirror plane and not all of them
have the same parity, those of opposite parity will generally have a nonzero
$PzP$ matrix element between them, and will tend to hybridize and split to
form dispersive pairs. Thus, all flat bands pinned at A are expected to have
the same parity $p_{\rm A}$, and all flat bands pinned at B are expected to
have the same parity $p_{\rm B}$. Following Ref. Varnava  _et al._ (2020), we
call this the “uniform parity” assumption. As discussed in Ref. Varnava  _et
al._ (2020), this assumption is closely related to a well-known theorem on the
minimum number of zero-energy modes in bipartite lattices Sutherland (1986);
Lieb (1989); Ramachandran  _et al._ (2017).

Under the uniform parity assumption, the numbers $\overline{N}_{\rm A}$ and
$\overline{N}_{\rm B}$ of flat bands at A and B can be expressed in terms of
the imbalance between even- and odd-parity valence Bloch states at the mirror-
invariant plane(s) in the BZ. For a type-1 mirror we have

$\overline{N}_{\rm A}=\frac{1}{2}\left|\Delta N_{\rm G}+\Delta N_{\rm
X}\right|$ (20)

and

$\overline{N}_{\rm B}=\frac{1}{2}\left|\Delta N_{\rm G}-\Delta N_{\rm
X}\right|\,,$ (21)

where $\Delta N_{\rm G}$ and $\Delta N_{\rm X}$ denote the excess of even over
odd valence states at G and X, respectively. Hence if the mirror-parity
content is balanced at both G and X, flat Wannier bands are absent from both A
and B; if it is balanced only at G but not at X or vice-versa, the same number
of flat bands is present at A and at B; and if it is unbalanced at both G and
X, the number of flat bands at B can differ from the number at A. The
corresponding relation for a type-2 mirror is

$\overline{N}_{\rm A}=\overline{N}_{\rm B}=\frac{1}{2}|\Delta N_{\rm G}|\,.$
(22)

Equations (20-22) are derived in Appendix A.

###  III.2 Types of generic degeneracies

In this section, we consider the types of degeneracies that are typical of the
Wannier spectra of insulators with $M_{z}$ symmetry. We call a degeneracy
generic when it occurs without the assistance of any symmetries other than
$M_{z}$. If in addition the degeneracy is codimension protected, we call it
accidental.

Accidental degeneracies away from the A and B planes have codimension three,
and hence they require fine tuning. On the mirror planes, there are two types
of generic degeneracies: multiple flat bands pinned to the same plane, and
accidental touchings, at isolated points in the 2D BZ, between one or more
pairs of dispersive bands. Other possibilities such as nodal lines are non-
generic and will not be considered further. In the following we focus on the A
plane $z=0$, but the discussion would be identical for the B plane $z=c/2$.

####  III.2.1 Point nodes between pairs of dispersive bands

If there are no flat bands pinned at $z=0$, any bands near $z=0$ must come in
dispersive pairs at $\pm z$. If there is a single pair, we construct from the
two HW functions at each ${\bm{\kappa}}$ a pair of orthogonal states with
opposite parities about $z=0$. In this basis, the $z$ operator is represented
by a matrix of the form

$\begin{pmatrix}0&f_{\bm{\kappa}}\\\ f^{*}_{\bm{\kappa}}&0\end{pmatrix}\,,$
(23)

with eigenvalues $z_{\bm{\kappa}}=\pm|f_{\bm{\kappa}}|$. The two bands touch
at $z=0$ when $|f_{\bm{\kappa}}|=0$, and for that to happen both the real and
imaginary parts of $f_{\bm{\kappa}}$ must vanish; this means that such
degeneracies have codimension two, and hence they occur at isolated points in
the 2D BZ. (When the bands disperse linearly close to the nodal point, the
degeneracy is called a “Dirac node.”) If more than one dispersive band pair is
involved, $f_{\bm{\kappa}}$ becomes a matrix. The degeneracy condition
$\det(f_{\bm{\kappa}})=0$ again leads to point nodes on the $z=0$ plane.
Generically, these are simple nodes where only two bands meet. However, with
additional symmetries or fine tuning, more than one pair of bands may become
degenerate at a given node.

In summary, pairs of dispersive Wannier bands can touch accidentally at
isolated points on a mirror plane free of flat bands. We note that the same
happens, and for the same mathematical reasons, with the energy bands of
models with sublattice symmetry Ramachandran  _et al._ (2017).

####  III.2.2 Flat bands repel point nodes

When one or more flat bands are present at $z=0$, they gap out the point
nodes. Let us show this for the simplest case of one flat band surrounded by a
dispersive pair. Choosing a basis of $M_{z}$ eigenstates within this three-
band space, the matrix representation of the $z$ operator takes the form

$\begin{pmatrix}0&f_{\bm{\kappa}}&g_{\bm{\kappa}}\\\
f^{*}_{\bm{\kappa}}&0&0\\\ g^{*}_{\bm{\kappa}}&0&0\end{pmatrix},$ (24)

where we have chosen the first basis state to have the opposite mirror parity
from the other two. The eigenvalues are $z_{\bm{\kappa}}=0$ (flat band) and
$z_{\bm{\kappa}}=\pm\sqrt{|f_{\bm{\kappa}}|^{2}+|g_{\bm{\kappa}}|^{2}}$
(dispersive pair). An accidental degeneracy between the pair requires the real
and imaginary parts of both $f_{\bm{\kappa}}$ and $g_{\bm{\kappa}}$ to vanish
(codimension four). In general this cannot be achieved by adjusting
${\bm{\kappa}}$ alone; it also requires fine tuning the parameters
$f_{\bm{\kappa}}$ and $g_{\bm{\kappa}}$.

In conclusion, flat bands and point nodes do not generally coexist on a mirror
plane. Although we have only shown this for the case of one flat band plus one
dispersive pair, the same result is expected to hold when several flat bands
and/or dispersive pairs are present. That scenario has in fact been considered
for the analogous problem of energy bands in models with sublattice symmetry
Ramachandran  _et al._ (2017).

####  III.2.3 Spinful time-reversal symmetry excludes flat bands

The presence of flat bands on the mirror planes can sometimes be ruled out on
the basis of symmetry. This is the case for a crystal that has both $M_{z}$
symmetry and spinful time-reversal symmetry $\mathcal{T}$. Since
$[P_{\bm{\kappa}}zP_{\bm{\kappa}},\mathcal{T}]=0$, the standard Kramers-
degeneracy argument applies to the Wannier bands: if $|h_{\bm{\kappa}}\rangle$
is an eigenstate of $P_{\bm{\kappa}}zP_{\bm{\kappa}}$ with eigenvalue
$z_{\bm{\kappa}}$, then
$|h^{\prime}_{-{\bm{\kappa}}}\rangle=\mathcal{T}|h_{\bm{\kappa}}\rangle$ is an
orthogonal eigenstate with the same eigenvalue. Now suppose that
$|h_{\bm{\kappa}}\rangle$ is a flat-band state at A, with $M_{z}$ eigenvalue
$\lambda=\pm i$. Then $|h^{\prime}_{-{\bm{\kappa}}}\rangle$ is also a flat-
band state, and using $[M_{z},\mathcal{T}]=0$ we find that its mirror
eigenvalue is $\lambda^{*}=-\lambda$. Since the two flat bands have opposite
mirror eigenvalues, they will generally hybridize to form a dispersive pair.

Another example is a crystal that has both $M_{z}$ symmetry, and spinful
$\mathcal{T}$ combined with inversion $\mathcal{I}$. The combined symmetry
$\mathcal{I}*\mathcal{T}$ renders the energy bands Kramers-degenerate at every
${\bf k}$, and since $[M_{z},\mathcal{I}*\mathcal{T}]=0$ and $M_{z}$ has
purely imaginary eigenvalues, Kramers pairs of Hamiltonian eigenstates on the
invariant BZ planes have opposite $M_{z}$ eigenvalues. The mirror-parity
content is therefore balanced on those planes, and from Eqs. (20-22) we
conclude that both $\overline{N}_{\rm A}$ and $\overline{N}_{\rm B}$ vanish.
(Note that while the energy bands are Kramers degenerate in the presence of
$\mathcal{I}*\mathcal{T}$ symmetry, the Wannier bands are not. The difference
is that $\mathcal{I}*\mathcal{T}$ commutes with the Hamiltonian, but it
anticommutes with $PzP$.)

In summary, spinful time-reversal symmetry, either by itself or in combination
with inversion, rules out the presence of flat Wannier bands on the mirror
planes (under the uniform parity assumption).

###  III.3 Chern numbers in gapped band structures

When an $M_{z}$-symmetric Wannier band structure is gapped, the $J$ bands per
cell can be grouped into three internally connected collections Varnava  _et
al._ (2020): one containing bands that are pinned at A (over the entire 2D BZ
or at isolated ${\bm{\kappa}}$ points), another containing bands that are
pinned at B, and a third containing “unpinned” bands, in the sense that they
do not touch the mirror planes anywhere in the 2D BZ. In Ref. Varnava  _et
al._ (2020) these three collections were called origin-centered, boundary-
centered, and uncentered, respectively.

Letting $\alpha=\text{A or B}$, in each vertical cell $l$ there are in general

  * •

$\overline{N}_{\alpha^{+}}$ flat bands at $\alpha$ of even parity,

  * •

$\overline{N}_{\alpha^{-}}$ flat bands at $\alpha$ of odd parity,

  * •

$\widetilde{N}_{\alpha}$ dispersive bands touching at $\alpha$

in the $\alpha$-pinned collection, and $\widetilde{N}_{\rm UC}$ dispersive
bands in the unpinned collection. (At this stage we do not yet assume uniform
parity for the flat bands, nor do we invoke the fact that flat bands repel
point nodes.) In the home cell $l=0$, the dispersive bands in the A-pinned
collection come in pairs at $\pm z$, and those in the B-pinned collection come
in pairs at $z$ and $c-z$. In the case of the unpinned collection we have a
choice, since the mirror-symmetric partners never become degenerate; for
definiteness, we choose the contents of the home cell so that the bands in the
unpinned collection come in pairs at $\pm z$.

For each of the seven groups listed above, we can add up the Chern numbers in
that group to get $\overline{C}_{\alpha^{\pm}}$, $\widetilde{C}_{\alpha}$, and
$\widetilde{C}_{\rm UC}$, keeping in mind that their sum $C_{3}$ vanishes by
assumption,

$C_{\rm A}+C_{\rm B}+\widetilde{C}_{\rm UC}=0\,,$ (25)

where
$C_{\alpha}=\overline{C}_{\alpha^{+}}+\overline{C}_{\alpha^{-}}+\widetilde{C}_{\alpha}$
is the net Chern number of the $\alpha$-pinned collection. We further
decompose each of the three dispersive band subspaces into even and odd
sectors under reflection about their centers, and assign separate Chern
numbers to them,

$\displaystyle\widetilde{C}_{\alpha}$
$\displaystyle=\widetilde{C}_{\alpha^{+}}+\widetilde{C}_{\alpha^{-}}\,,$ (26a)
$\displaystyle\widetilde{C}_{\rm UC}$ $\displaystyle=\widetilde{C}_{{\rm
UC}^{+}}+\widetilde{C}_{{\rm UC}^{-}}\,.$ (26b)

In Appendix B we show that

$\widetilde{C}_{\alpha^{+}}-\widetilde{C}_{\alpha^{-}}=W_{\alpha}\,,$ (27)

where $W_{\alpha}$ is the sum of the winding numbers of all the nodal points
in the projected 2D BZ on the $\alpha$ mirror plane.

The winding number of a nodal point ${\bm{\kappa}}_{j}$ is defined as Asbóth
_et al._ (2016)

$W_{j}=\frac{1}{2\pi}\oint_{c_{j}}\partial_{\bm{\kappa}}\gamma_{\bm{\kappa}}\cdot
d{\bm{\kappa}}\,,$ (28)

where the integral is over a small circle around the node. $W_{j}$ is an
integer, typically taking values $\pm 1$ according to how the phase
$\gamma_{\bm{\kappa}}$ changes going around the node. In the simplest case
where a single pair of bands meet at the node, $\gamma_{\bm{\kappa}}$ is the
phase angle of the complex matrix element $f_{\bf k}$ appearing in Eq. (23).
If two or more pairs of bands meet at a node, $f_{\bf k}$ becomes a matrix and
$\gamma_{\bm{\kappa}}$ becomes the phase angle of its determinant (see Sec.
V.3).

Combining Eqs. (26a) and (27) we obtain

$W_{\alpha}=\widetilde{C}_{\alpha}-2\widetilde{C}_{\alpha^{-}}\,,$ (29)

which shows that $\widetilde{C}_{\alpha}$ has the same even or odd parity as
$W_{\alpha}$. Since band pairs in the unpinned collection do not touch on the
special planes, by applying the same argument in Appendix B that leads to Eq.
(27) we obtain

$\widetilde{C}_{{\rm UC}^{+}}=\widetilde{C}_{{\rm UC}^{-}}\,,$ (30)

which implies that their sum $\widetilde{C}_{\rm UC}$ is always an even
number.222The fact that $\widetilde{C}_{\rm UC}$ is even can also be seen as
follows Varnava  _et al._ (2020). The unpinned collection is formed by two
disconnected groups of bands related by $M_{z}$ symmetry, which imposes the
same Berry curvature at every ${\bm{\kappa}}$ in the two groups, and hence the
same Chern number.

##  IV Mirror Chern numbers in the hybrid Wannier representation

We are finally ready to evaluate the MCNs in the HW representation, and then
relate them to the axion $\mathbbm{Z}_{2}$ index. In Sec. IV.1 we consider the
case of a gapped Wannier spectrum, and in Sec. IV.2 we treat the gapless case.

###  IV.1 Gapped Wannier band structure

To recap, a generic gapped Wannier band structure with $M_{z}$ symmetry
consists of seven band collections per cell. The four that are flat have well-
defined mirror parities, and the three that are dispersive can be decomposed
into even and odd sectors. This yields a total of ten HW groups with well-
defined parities, each carrying its own Chern number.

####  IV.1.1 Type-1 mirrors

Table 1:  Parities under a type-1 mirror $M_{z}$ of Bloch-like states
constructed from HW functions that are maximally localized along $z$. For
spinful electrons, the parity is said to be “even” or “odd” when the $M_{z}$
eigenvalue is $+i$ or $-i$. Bloch representation
---
${\rm G}^{+}=\text{even about A (and even about B)}$
${\rm G}^{-}=\text{odd\,\, about A (and odd\,\, about B)}$
${\rm X}^{+}\,=\text{even about A (and odd\,\, about B)}$
${\rm X}^{-}\,=\text{odd\,\, about A (and even about B)}$
Hybrid Wannier representation
${\rm A}^{+}$ = even about A, generates ${\rm G}^{+}$ and ${\rm X}^{+}$
${\rm A}^{-}$ = odd about A, generates ${\rm G}^{-}$ and ${\rm X}^{-}$
${\rm B}^{+}$ = even about B, generates ${\rm G}^{+}$ and ${\rm X}^{-}$
${\rm B}^{-}$ = odd about B, generates ${\rm G}^{-}$ and ${\rm X}^{+}$
pairs C and ${\rm C}^{\prime}$, generates ${\rm G}^{+}{\rm G}^{-}$ and ${\rm
X}^{+}{\rm X}^{-}$

To evaluate the MCNs $\mu_{\rm G}$ and $\mu_{\rm X}$, we construct from each
of the ten HW groups a group of Bloch-like states by performing Bloch sums
along $z$, and recall from Eq. (19) that their Chern numbers on any
constant-$k_{z}$ BZ slice (and, in particular, at G and X) are the same as the
Chern numbers of the parent HW groups. The final needed ingredient is Table 1,
which tells the mirror parities at G and X of the Bloch groups coming from
each of the HW groups. That table is valid for both spinless and spinful
mirror symmetry $M_{z}$, and it agrees with the parity rules for inversion
symmetry ${\cal I}$ in 1D Varnava  _et al._ (2020); this is consistent with
the fact that $M_{z}={\cal I}*C_{2}^{z}$ acts along $z$ in the same way as
${\cal I}$.

To evaluate $\mu_{\rm G}$, we need to split the occupied Bloch space at G into
even- and odd-parity sectors about A. According to Table 1, their Chern
numbers are

$C_{\rm G}^{\pm}=\left(\overline{C}_{{\rm A}^{\pm}}+\widetilde{C}_{{\rm
A}^{\pm}}+\widetilde{C}_{\rm UC}^{\pm}\right)+\left(\overline{C}_{{\rm
B}^{\pm}}+\widetilde{C}_{{\rm B}^{\pm}}\right)\,,$ (31)

where the first and second groups of terms correspond to Wannier groups that
are even or odd about A and B, respectively. Inserting this expression into
Eq. (8) for $\mu_{\rm G}$ and then using Eqs. (27) and (30), we find

$2\mu_{\rm G}=\left(\overline{C}_{{\rm A}^{+}}-\overline{C}_{{\rm
A}^{-}}\right)+\left(\overline{C}_{{\rm B}^{+}}-\overline{C}_{{\rm
B}^{-}}\right)+W_{\rm A}+W_{\rm B}\,.$ (32)

Under the uniform parity assumption the first group of terms becomes $p_{\rm
A}\overline{C}_{\rm A}$, where $\overline{C}_{\rm A}$ is the total Chern
number of the flat bands at A, all of the same parity $p_{\rm A}=\pm 1$;
similarly, the second group becomes $p_{\rm B}\overline{C}_{\rm B}$. Thus we
arrive at

$\mu_{\rm G}=\frac{1}{2}\left(p_{\rm A}\overline{C}_{\rm A}+W_{\rm
A}\right)+\frac{1}{2}\left(p_{\rm B}\overline{C}_{\rm B}+W_{\rm B}\right)\,,$
(33)

and via similar steps Eq. (9) for $\mu_{\rm X}$ turns into

$\mu_{\rm X}=\frac{1}{2}\left(p_{\rm A}\overline{C}_{\rm A}+W_{\rm
A}\right)-\frac{1}{2}\left(p_{\rm B}\overline{C}_{\rm B}+W_{\rm B}\right)\,.$
(34)

Out of the three collections in a type-1 disconnected band structure, the
uncentered collection does not contribute to the MCNs; and the A-centered and
B-centered ones contribute as in Eqs. (33) and (34).

Equations (33) and (34) are a central result of this work, and in the
following sections we will extract several conclusions from them. In practical
applications, those equations can often be simplified: since flat bands and
point nodes do not generically coexist on the mirror planes, at least one of
the two terms inside each parenthesis will typically vanish.

Before proceeding, let us verify that Eq. (33) correctly yields an integer
value for $\mu_{\rm G}$ when $C_{3}=0$. First we eliminate the winding numbers
from Eq. (33) with the help of Eq. (29), and then we take mod 2 on both sides
of the resulting equation to find

$\displaystyle 2\mu_{\rm G}\text{ mod 2}$
$\displaystyle=\left(\overline{C}_{\rm A}+\widetilde{C}_{\rm
A}+\overline{C}_{\rm B}+\widetilde{C}_{\rm B}\right)\text{ mod 2}$
$\displaystyle=-\widetilde{C}_{\rm UC}\text{ mod 2}\,,$ (35)

where Eq. (25) was used to go from the first to the second line. Given that
$\widetilde{C}_{\rm UC}$ is an even number, we conclude that $\mu_{\rm G}$ is
an integer. The proof is identical for Eq. (34).

We emphasize that the separate contributions from the A- and B-centered
collection to Eqs. (33) and (34) are not always integer-valued. As can be seen
from Eq. (37) below, those contributions assume half-integer values when the
axion angle is quantized to $\theta=\pi$ by mirror symmetry; a concrete
example where this happens will be given in Sec. VI.3.

####  IV.1.2 Relation to the quantized axion coupling

As mentioned in the Introduction, mirror symmetry belongs to the group of
“axion-odd” symmetries that reverse the sign of the axion angle $\theta$. When
one or more such symmetries are present in a 3D insulator with a vanishing
Chern vector, $\theta$ is restricted to be zero or $\pi$ mod $2\pi$, becoming
a $\mathbbm{Z}_{2}$ topological index.

In the case of mirror symmetry, where the band topology is already
characterized by the MCNs, there should be a relation between them and the
quantized $\theta$ value. Below we derive that relation for an insulator with
a type-1 mirror and a gapped Wannier spectrum. To that end, we make use of the
formalism of Ref. Varnava  _et al._ (2020) for expressing $\theta$ in the HW
representation.

First we write $\mu_{\rm G}+\mu_{\rm X}$ by combining Eqs. (33) and (34), and
eliminate the winding numbers using Eq. (29). Then we take mod 2 on both sides
to find

$\left(\mu_{\rm G}+\mu_{\rm X}\right)\text{ mod 2}=C_{\rm A}\text{ mod 2}\,.$
(36)

Comparing with the relation $\theta/\pi=C_{\rm A}\text{ mod 2}$ Varnava  _et
al._ (2020), valid for a gapped spectrum in the presence of a $z$-reversing
axion-odd symmetry such as $M_{z}$, we conclude that

$\frac{\theta}{\pi}=\left(\mu_{\rm G}+\mu_{\rm X}\right)\text{ mod 2}\,.$ (37)

Thus, the system is axion-even ($\theta=0$) or axion-odd ($\theta=\pi$)
depending on whether the sum of the two MCNs associated with $M_{z}$ is even
or odd. Previously, this result had been inferred from an argument based on
counting Dirac cones in the surface BZ Varjas  _et al._ (2015); Fulga  _et
al._ (2016). Here, we have obtained it directly as a formal relation between
bulk quantities expressed in the HW representation. As we will see shortly,
the same relation holds when the Wannier spectrum is gapless.

####  IV.1.3 Type-2 mirrors

In a crystal with a type-2 mirror, where the planes A and B are equivalent and
G is the only mirror-invariant plane in reciprocal space, the unique MCN
$\mu_{\rm G}$ is obtained by setting $p_{\rm B}=p_{\rm A}$, $\overline{C}_{\rm
B}=\overline{C}_{\rm A}$, and $W_{\rm B}=W_{\rm A}$ in Eq. (33),

$\mu_{\rm G}=p_{\rm A}\overline{C}_{\rm A}+W_{\rm A}\,.$ (38)

If flat bands are present at A, they repel the point nodes. Hence $W_{\rm
A}=0$, and therefore $|\mu_{\rm G}|=|\overline{C}_{\rm A}|$. Interestingly, in
this case the magnitude of the MCN does not depend on the parity of the flat-
band states; this simplifies considerably its numerical evaluation, since one
does not need to know how the basis orbitals transform under $M_{z}$. Given
that only the magnitude (not the sign) of the MCN is needed to establish the
bulk-boundary correspondence, this is a potentially useful result.

Inserting Eq. (29) for $W_{\rm A}$ in Eq. (38), taking mod 2 on both sides,
and again comparing with $\theta/\pi=C_{\rm A}\text{ mod 2}$, we conclude that
in this case the relation between the axion $\mathbbm{Z}_{2}$ index and the
MCN reads

$\frac{\theta}{\pi}=\mu_{\rm G}\text{ mod 2}\,,$ (39)

as stated in Ref. Fulga  _et al._ (2016).

####  IV.1.4 Weakly coupled layered crystals

Consider a crystal composed of weakly coupled identical layers that remain
invariant under reflection about their own planes. Following Ref. Kim  _et
al._ (2015), we assume that the layers are stacked exactly vertically. In this
case the reflection symmetry about the individual layers becomes a type-1
mirror of the 3D structure, with two separate MCNs $\mu_{\rm G}$ and $\mu_{\rm
X}$. In the fully decoupled limit where there is no $k_{z}$ dependence the G
and X reciprocal planes become equivalent, so that $\mu_{\rm X}=\mu_{\rm
G}\equiv\mu_{\rm 2D}$ where $\mu_{\rm 2D}$ is the MCN of an isolated layer
[Eq. (10)]. But since the MCNs are integers, they cannot change if a weak
interlayer coupling is introduced, and from Eqs. (33) and (34) we obtain

$\mu_{\rm 2D}=\frac{1}{2}\left(p_{\rm A}\overline{C}_{\rm A}+W_{\rm A}\right)$
(40)

for the unique MCN of a weakly-coupled layered crystal.

If flat bands are present at A (the plane of a layer), then $W_{\rm A}=0$ and
the net Chern number of the valence bands becomes $\overline{C}_{\rm
A}+\widetilde{C}_{\rm UC}$; since the net Chern number vanishes by assumption
and $\widetilde{C}_{\rm UC}$ is even, $\mu_{\rm 2D}=p_{\rm A}\overline{C}_{\rm
A}/2$ is clearly an integer. In this case $|\mu_{\rm 2D}|$ can be determined
without knowing the parity of the flat-band states, as in the case of a type-2
mirror with flat bands.

Let us now evaluate the axion $\mathbbm{Z}_{2}$ index. Since $\mu_{\rm
G}+\mu_{\rm X}=2\mu_{\rm 2D}$ is an even number, Eq. (37) yields

$\theta=0\text{ mod $2\pi$}\,.$ (41)

This is consistent with the assertion made in Ref. Kim  _et al._ (2015) that
weakly-coupled layered topological crystalline insulators are analogous to
“weak topological insulators” with a vanishing strong $\mathbbm{Z}_{2}$
invariant $\nu_{0}$.

###  IV.2 Gapless Wannier band structure

Let us now apply our formalism to a $M_{z}$-symmetric system with a gapless
Wannier spectrum. We start out by noting that such a spectrum must have
degeneracies at both A and B. On those special planes the codimension is two,
so point nodes are allowed. Flat bands can be ruled out since they would repel
any nodes and generate a gap, and we assume that nodal lines are absent as
well.

We are left with a scenario where there are point nodes at both A and B, and
these are connected by Wannier bands. The only way this can happen without the
assistance of other symmetries is if there are only two Wannier bands, one in
each half unit cell, since otherwise there is generically a gap somewhere in
each half cell (accidental degeneracies away from A and B are not protected,
since the codimension is three). With the assistance of other symmetries, the
gapless spectrum may contain more than two bands per cell.

To treat the above scenario, we temporarily add a symmetric pair of occupied
orbitals at degeneracy-free planes $\pm z_{0}$, and initially do not let them
hop at all (completely isolated). This will introduce flat bands on those
planes. Now let the added orbitals hybridize with other orbitals. Since
accidental degeneracies away from the mirror planes are not protected, gaps
will generally open up between the new and the old Wannier bands (the only
exceptions to this rule are treated in the next paragraph). And since the
added orbitals are topologically trivial, they have no effect on the MCNs,
which can now be evaluated using the formalism of Sec. IV.1 for gapped
spectra. Setting $\overline{C}_{\rm A}=\overline{C}_{\rm B}=0$ in Eqs. (33)
and (34) therein, we obtain

$\mu_{\rm G}=\frac{1}{2}\left(W_{\rm A}+W_{\rm B}\right)$ (42)

and

$\mu_{\rm X}=\frac{1}{2}\left(W_{\rm A}-W_{\rm B}\right)\,.$ (43)

But since $W_{\rm A}$ and $W_{\rm B}$ cannot be affected by orbitals inserted
far from the A and B planes, we conclude that Eqs. (42) and (43) can be
directly applied to the original system with a gapless Wannier spectrum.

The above argument needs to be refined if the system is an axion-odd insulator
that has, in addition to $M_{z}$ symmetry, one or more axion-odd symmetries
that are $z$ preserving and symmorphic (e.g., spinful time reversal or
vertical mirrors). The Wannier spectrum is then guaranteed to be gapless, with
adjacent bands touching at an odd number of Dirac nodes Varnava  _et al._
(2020). The solution is to weakly break all such symmetries via some low-
symmetry perturbation; the band connectivity then becomes “fragile,” allowing
gaps to open up once the added orbitals hybridize with the original ones
Wieder and Bernevig (2018); Varnava  _et al._ (2020). The rest of the argument
proceeds as before, again with the conclusion that Eqs. (42) and (43) can be
directly applied to the original system with a gapless spectrum. This scenario
is illustrated in Sec. VI.3.2, where the orbital insertion itself acts as the
symmetry-lowering perturbation.

To conclude, let us show that the relation (37) between the MCNs and the axion
angle remains valid for gapless spectra. Equations (42) and (43) give
$\mu_{\rm G}+\mu_{\rm X}=W_{\rm A}$, while $\theta$ is equal to the sum of
Berry phases of vanishingly small loops around the nodes at A Varnava  _et
al._ (2020). Since those Berry phases divided by $\pi$ are equal to the node
winding numbers modulo 2 Park and Marzari (2011), Eq. (37) is immediately
recovered.

##  V Methods

###  V.1 Tight-binding, ab initio, and Wannier methods

In this work, the formalism for evaluating MCNs in the HW representation is
implemented in the tight-binding (TB) framework, using a modified version of
the PythTB code pyt . Illustrative calculations are carried out for 2D and 3D
models with mirror symmetry; some are simple toy models, while others are
obtained from ab initio calculations as described below. Each model is
specified by providing the on-site energies, the hopping amplitudes, and the
matrix elements of the position and mirror operators.

In the TB literature, it is common to assume that the position operator is
represented by a diagonal matrix in the TB basis,

$\langle\varphi_{{\bf R}i}|{\bf r}|\varphi_{{\bf R}^{\prime}j}\rangle=({\bf
R}+\bm{\tau}_{i})\delta_{{\bf R},{\bf R}^{\prime}}\delta_{ij}$ (44)

where $\bm{\tau}_{i}$ is the location of the $i$th basis orbital in the home
cell ${\bf R}=\mathbf{0}$. This approximation is problematic for calculating
the Wannier bands of unbuckled monolayers, since it forces all bands to lie
flat on the $z=0$ plane: when all basis orbitals lie on the $z=0$ plane and
all off-diagonal matrix elements $\langle\varphi_{{\bf R}i}|z|\varphi_{{\bf
R}^{\prime}j}\rangle$ vanish, the matrix $Z_{\bm{\kappa}}$ that is
diagonalized to obtain the HW centers [see Eqs. (45) and (46)] is the null
matrix.

To apply our formalism to flat monolayers, any flat Wannier bands that may be
present must be robust and satisfy the uniform parity assumption, while all
other bands must be dispersive. To ensure that this is so, one should retain
some off-diagonal $z$ matrix elements. For models based on ab initio Wannier
functions this occurs naturally, since the position matrix elements between
the Wannier functions are explicitly calculated, and they are generally
nonzero for nearby Wannier functions. In the case of toy models, one needs to
assign nonzero values to some of the off-diagonal $z$ matrix elements under
reasonable assumptions.

The material chosen for the ab initio calculations is SnTe, which we study as
a flat monolayer in Sec. VI.1 and as a bulk phase in Sec. VI.2. We first
calculate the electronic structure from density-functional theory (DFT) using
the GPAW code Enkovaara  _et al._ (2010), and then use the Wannier90 code
Mostofi  _et al._ (2014) to construct well-localized Wannier functions.
Lastly, TB models are generated by tabulating the matrix elements of the Kohn-
Sham Hamiltonian and of the position operator between those Wannier functions.

The self-consistent DFT calculations are performed without including spin-
orbit coupling, which is added afterwards non-selfconsistently Olsen (2016).
We use the Perdew-Burke-Ernzerhof exchange-correlation functional Perdew  _et
al._ (1996, 1997), and describe the valence-core interaction via the projector
augmented wave method Blöchl (1994). The valence states are expanded in a
plane-wave basis with an energy cutoff of 600 eV, and the BZ is sampled on
$\Gamma$-centered uniform grids containing $6\times 6\times 1$ and $6\times
6\times 6$ points for monolayer and bulk SnTe, respectively. The projector
augmented wave setup includes the 4$d$ semicore states of Sn in addition to
the 5$s$ and 5$p$ states of Sn and Te, yielding a total of 20 valence
electrons for each SnTe formula unit (one per cell for the monolayer, and two
for the bulk).

For each formula unit, we construct 16 spinor Wannier functions of $s$ and $p$
character spanning the upper-valence and low-lying conduction band states. The
Sn 4$d$ states, which give rise to flat bands lying 22 eV below the Fermi
level, are excluded from the Wannier construction.

As a first step towards obtaining well-localized Wannier functions, we extract
from the space of ab initio Bloch eigenstates at each grid point ${\bf k}$ an
$N$-dimensional subspace with the desired orbital character ($N=16$ for the
monolayer, and $N=32$ for the bulk). This is achieved via the “band
disentanglement” procedure of Ref. Souza  _et al._ (2001), which involves
specifying two energy windows, known as the inner and the outer window, and a
set of trial orbitals. The outer window encloses all the valence bands except
for the 4$d$ semicore states, as well as all the low-lying conduction states
of $5s$ and $5p$ character. To ensure that the valence states are exactly
preserved in the disentangled subspace, we “freeze” them inside an inner
window. An initial guess for the target subspace is obtained by projecting
atom-centered $s$ and $p$ trial orbitals onto the outer-window states. This is
followed by an iterative procedure that yields an optimally-smooth
disentangled subspace across the BZ Souza  _et al._ (2001).

Having extracted a suitable Bloch subspace, we proceed to construct well-
localized $s$\- and $p$-like Wannier functions spanning that subspace. This is
done by projecting onto it the same $s$ and $p$ trial orbitals that were used
in the disentanglement step, and then orthogonalizing the resulting orbitals
via the Löwdin scheme Marzari and Vanderbilt (1997). This one-shot procedure,
without additional maximal-localization steps Marzari and Vanderbilt (1997),
ensures that the Wannier functions retain the orbital character of the trial
orbitals.

To assess the quality of the Wannier basis we calculate the energy bands from
the Hamiltonian matrix elements in that basis Souza  _et al._ (2001), and find
that they are in excellent agreement with the ab initio bands obtained using
the GPAW code Olsen  _et al._ (2019).

In addition to the Hamiltonian and position matrix elements, we also require
the matrix elements of the mirror operator $M_{z}$ in the Wannier basis. These
are needed to determine the winding numbers of the nodal touchings between
Wannier bands on the mirror planes (see Sec. V.3), as well as the mirror
parities $p_{\rm A}$ and $p_{\rm B}$ of the flat-band states. To set up the
matrix representation of $M_{z}$, we assume that the Wannier functions
transform under $M_{z}$ in the same way as pure $s$ and $p$ orbitals. We find
that the eigenstates of the Wannier Hamiltonian on the mirror-invariant BZ
planes are, to a good approximation, eigenstates of this approximate $M_{z}$
operator, which validates that assumption.

###  V.2 Construction of hybrid Wannier functions and Wannier bands

Formally, maximally-localized HW functions satisfy the eigenvalue equation
(12). For a 2D or quasi-2D system extended along $x$ and $y$, the matrix
elements of the $z$ operator appearing in that equation are well defined. It
is therefore straightforward to set up the matrix

$Z_{mn{\bf k}}=\langle\psi_{m{\bf k}}|z|\psi_{n{\bf k}}\rangle\,,$ (45)

where ${\bf k}=(k_{x},k_{y})$ and $m$ and $n$ run over the $J$ occupied energy
bands, and to diagonalize it,

$\left[U^{\dagger}_{\bf k}Z_{\bf k}U_{\bf k}\right]_{mn}=z_{m{\bf
k}}\delta_{mn}\,.$ (46)

The eigenvalues are the HW centers, and from the eigenvectors (the columns of
the $U_{\bf k}$ matrix) we can construct the maximally-localized HW functions
according to

$|h_{n{\bf k}}\rangle=\sum_{m}\,e^{-i{\bf k}\cdot{\bf r}}|\psi_{m{\bf
k}}\rangle U_{mn{\bf k}}\,,$ (47)

where the phase factor has been included to render them in-plane periodic.

For bulk systems, which are extended in all directions including the
wannierization direction $z$, the above procedure fails because the matrix
elements in Eq. (45) become ill defined. In such cases, it is still possible
to construct maximally-localized HW functions by working in reciprocal space.
We now write ${\bf k}=({\bm{\kappa}},k_{z})$, and choose a uniform grid; for
each point ${\bm{\kappa}}$ in the projected 2D BZ, the problem reduces to the
construction of 1D maximally-localized Wannier functions along $z$. The
procedure is detailed in Refs. Vanderbilt (2018); Marzari and Vanderbilt
(1997). Briefly, the first step is to establish a “twisted parallel transport
gauge” for the valence Bloch states along the string of $k_{z}$ points at each
${\bm{\kappa}}$, obtaining as a byproduct the HW centers
$z_{ln{\bm{\kappa}}}$. The maximally-localized HW functions
$|h_{ln{\bm{\kappa}}}\rangle$ are then constructed in this gauge using Eq.
(11), with the integral over $k_{z}$ replaced by a summation over the string
of $k_{z}$ points.

###  V.3 Winding number of a point node of order $N$

####  V.3.1 Definition

Earlier, we defined the winding number of a point node where two Wannier bands
meet on a mirror plane. Since there are situations where $N>1$ pairs of bands
meet at a node, we need to generalize that definition to handle such “higher-
order” nodes.

Given a point node ${\bm{\kappa}}_{j}$ of order $N\geq 1$, we introduce the
$2N\times 2N$ matrix representation of $M_{z}$ at a nearby point
${\bm{\kappa}}$,

${\cal M}^{z}_{mn{\bm{\kappa}}}=\langle
h_{m{\bm{\kappa}}}|M_{z}|h_{n{\bm{\kappa}}}\rangle\,.$ (48)

Here, $m$ and $n$ run over the $2N$ Wannier bands that meet at
${\bm{\kappa}}_{j}$. By diagonalizing ${\cal M}^{z}_{\bm{\kappa}}$ and then
transforming the $|h_{n{\bm{\kappa}}}\rangle$ states accordingly [see Eqs.
(46) and (47)], we obtain a new set of $2N$ states
$|\tilde{h}_{n{\bm{\kappa}}}\rangle$. Like the original ones they are cell-
periodic in plane and localized along $z$, but they have definite mirror
parities. We choose the first $N$ to be even under $M_{z}$, and denote them as
$|\tilde{h}^{+}_{l{\bm{\kappa}}}\rangle$; the remaining $N$ are odd under
$M_{z}$, and we denote them as $|\tilde{h}^{-}_{l{\bm{\kappa}}}\rangle$. In
both cases, $l$ goes from 1 to $N$. The matrix representation of $z$ in the
new basis takes the form of Eq. (23), where $f_{\bm{\kappa}}$ is the $N\times
N$ matrix with elements

$f_{ll^{\prime}{\bm{\kappa}}}=\langle\tilde{h}^{+}_{l{\bm{\kappa}}}|z|\tilde{h}^{-}_{l^{\prime}{\bm{\kappa}}}\rangle\,.$
(49)

Letting

$\gamma_{\bm{\kappa}}=\arg(\det f_{\bm{\kappa}})\,,$ (50)

the winding number can be evaluated from Eq. (28) irrespective of the order
$N$ of the node.

####  V.3.2 Numerical evaluation

Suppose a single pair of Wannier bands meet at a point node
${\bm{\kappa}}_{j}$. To evaluate the winding number (28), the phase
$\gamma_{\bm{\kappa}}$ must be smooth on $c_{j}$. In practice, we establish a
smooth gauge for the states $|\tilde{h}^{\pm}_{\bm{\kappa}}\rangle$ as
follows. We pick a representation of the two states at a reference point
${\bm{\kappa}}^{\prime}_{j}$ in the vicinity of the node. Then at any point
${\bm{\kappa}}^{\prime}_{j}+\Delta{\bm{\kappa}}$ on the circle $c_{j}$ we
choose the gauge by enforcing maximal phase alignment with the states at
${\bm{\kappa}}^{\prime}_{j}$, i.e., by requiring that the overlaps
$\langle\tilde{h}^{+}_{{\bm{\kappa}}^{\prime}_{j}}|\tilde{h}^{+}_{{\bm{\kappa}}^{\prime}_{j}+\Delta{\bm{\kappa}}}\rangle$
and
$\langle\tilde{h}^{-}_{{\bm{\kappa}}^{\prime}_{j}}|\tilde{h}^{-}_{{\bm{\kappa}}^{\prime}_{j}+\Delta{\bm{\kappa}}}\rangle$
are real and positive. In other words, we carry out a one-step parallel
transport from ${\bm{\kappa}}^{\prime}_{j}$ to each circumference point.

If several pairs of bands meet at a node, the strategy is basically the same.
The only difference is that one must now use the multiband version of the
parallel-transport procedure Vanderbilt (2018); Marzari and Vanderbilt (1997).

##  VI Numerical results

In this section, we use our formalism to calculate the MCNs of three different
systems. The first is an unbuckled monolayer of SnTe, a topological
crystalline insulator protected by reflection symmetry about its plane. The
second is rocksalt SnTe, a 3D topological crystalline insulator protected by a
type-2 mirror. Our last example is a 3D toy model based on a modified Dirac
equation. It is both a strong topological insulator protected by time-reversal
symmetry, and a topological crystalline insulator with a type-1 mirror. In the
first example the Wannier spectrum is trivially gapped, while in the other two
it is gapless.

###  VI.1 Unbuckled monolayer of SnTe

Figure 2: (a) Atomic structure of monolayer SnTe. The black square is the
conventional unit cell with lattice constant $a$, and the red square is the
primitive cell with lattice constant $a^{\prime}=a/\sqrt{2}$. (b) Brillouin
zone and high-symmetry points.

Figure 3: (a) Energy bands of monolayer SnTe, with the $s$-type lower valence
bands that are exluded from the Wannierization shown in grey. All bands are
doubly degenerate, and the Fermi level is indicated by the dashed line. (b)
Wannier bands obtained from the Bloch states in the six $p$-type upper valence
bands. (c) Heatmap plot of the gap function of Eq. (52) for the central pair
of Wannier bands, where zero-gap points (nodal points) appear as dark spots.
Those with winding numbers $W_{j}=\pm 1$ are indicated by red or blue circles,
while the one with $W_{j}=-3$ at the $\Gamma$ point is indicated by a blue
triangle. Dashed circles denote pairs of nearby nodes with equal and opposite
winding numbers. When a node falls on the BZ boundary, only one of the
periodic images is shown.

The structure we consider is shown in Fig. 2(a). It consist of a single
unbuckled layer of Sn and Te atoms arranged in a checkerboard pattern, which
can be viewed as a single (001) layer of the bulk rocksalt structure.

DFT calculations reveal that the system with an optimized lattice constant of
$a=6.16$ Å is situated 0.4 eV above the convex hull and is dynamically
unstable Haastrup  _et al._ (2018), and that a buckled structure that breaks
mirror symmetry is energetically favored Kobayashi (2015). These results imply
that a flat SnTe monolayer is not likely to be experimentally relevant. This
system is nevertheless ideally suited for illustrating our methodology, since
it has reflection symmetry about its own plane and the associated MCN is
nonzero Liu  _et al._ (2015).

We carry out calculations using the primitive cell containing one formula
unit. The Wannier-interpolated energy bands are shown in Fig. 3(a), where all
bands are doubly degenerate due to time-reversal and inversion symmetry. There
is a robust inverted gap ($0.3$ eV) at the X point, and a tiny indirect gap
($0.17$ meV) around the X point; when the lattice expands the indirect gap
increases, and when it shrinks the system turns into a band overlap semimetal
Liu  _et al._ (2015); Kobayashi (2015). The lowest four valence bands are
predominantly $s$-type, and the remaining six (plotted in red) are
predominantly $p$-type.

Figure 3(b) shows the Wannier bands calculated from the Bloch states in the
$p$-type upper valence bands. The spectrum consists of three mirror-symmetric
band pairs that touch on the A plane $z=0$ at isolated points in the 2D BZ.
There are no flat bands on that plane, as expected from the presence of time-
reversal symmetry (Sec. III.2.3). Equation (40) therefore reduces to

$\mu_{\rm 2D}=\frac{1}{2}W_{\rm A}\,,$ (51)

and the MCN can be determined by evaluating the winding numbers of the nodal
points on the A plane.

To locate those nodal points, we plot in Fig. 3(c) the “gap function”

$g_{\bf k}=-\log(\Delta z_{\bf k}/c)\,,$ (52)

where $\Delta z({\bf k})$ is the separation between the central pair of bands.
Regions with a small gap appear in dark gray, and nodal points as dark spots.
The positions and winding numbers of all the nodal points are indicated in the
figure, where we have included only one of the periodic images when a node
falls on the BZ boundary. At $\Gamma$ and M there are nodes where three pairs
of Wannier bands touch, with winding numbers $W_{j}=-3$ and $W_{j}=+1$,
respectively. All other nodes on the $z=0$ plane are simple Dirac nodes where
only the two central bands meet, and they have $W_{j}=\pm 1$. Adding up the
winding numbers of the 36 nodal points in the BZ we obtain $W_{\rm A}=-4$, and
from Eq. (51) we conclude that the group of six $p$-type valence bands has a
MCN of $-2$.

We repeat the calculation for the four $s$-type lower valence bands, and find
that their net winding number vanishes. The net MCN of the occupied states is
therefore $\mu_{\rm 2D}=-2$, with the nontrivial topology coming from the $p$
states. This result agrees with the value $|\mu_{\rm 2D}|=2$ inferred from a
$k\cdot p$ analysis of the simultaneous band inversions at the two X points in
the BZ Liu  _et al._ (2014, 2015).

###  VI.2 Bulk SnTe

Bulk SnTe, which crystallizes in the rocksalt structure, is known both from
theory Hsieh  _et al._ (2012) and experiment Tanaka  _et al._ (2012) to be a
topological crystalline insulator. The symmetry protecting its nontrivial band
topology is reflection about the $\\{110\\}$ family of planes. (Instead, the
(001) mirror symmetry responsible for the topological state of the monolayer
is topologically trivial in the bulk crystal.)

The lattice is face-centered cubic lattice, so that the shortest lattice
vector perpendicular to the (110) planes is ${\bf a}_{3}=a\hat{\bf
x}/2+a\hat{\bf y}/2$. Since its length is twice the separation between
adjacent planes, the (110) mirror operation is of type 2, as is typical of
centered lattices (see Fig. 1).

For our simulations we pick a tetragonal cell subtended by ${\bf
a}_{1}=-a\hat{\bf x}/2+a\hat{\bf y}/2$, ${\bf a}_{2}=a\hat{\bf z}$, and ${\bf
a}_{3}$, and reorient the axes such that those vectors point along $\hat{\bf
x}$, $\hat{\bf y}$, and $\hat{\bf z}$, respectively. In this new frame, the
(110) mirror operation of interest becomes $M_{z}$. The simulation cell with
two formula units is shown in Fig. 4(a), and the associated BZ in Fig. 4(b).

Figure 4: (a) Rocksalt structure of bulk SnTe in a tetragonal conventional
cell. $a$ is the lattice constant of the conventional cubic cell, and
$b=c=a/\sqrt{2}$. Green planes are equivalent mirror planes. (b) Brillouin
zone associated with the tetragonal cell, with its high-symmetry points
indicated in red and the unique $M_{z}$-invariant plane in green. The
projected 2D Brillouin zone with its high-symmetry points is shown on top.

Figure 5: (a) Energy bands of bulk SnTe along high-symmetry lines of the
folded tetragonal BZ. The Fermi level is indicated by the dashed line. (b)
Wannier band structure obtained from the full set of valence states. (c)
Detail of the Wannier bands around the $z=0$ mirror plane. (d) Heatmap plot of
the gap function of Eq. (52) for the central pair of Wannier bands around
$z=0$, with the nodal points color-coded as in Fig. 3(c).

In Fig. 5(a) we present the energy bands calculated along the high-symmetry
lines of the folded BZ. The nontrivial topology arises from simultaneous band
inversions at the two L points in the unfolded BZ Hsieh  _et al._ (2012),
which map onto the two R points in Fig. 4(b). The inverted band gap at R and
the global indirect band gap amount to $0.3$ and $0.1$ eV, respectively.

From the full set of valence band states, we construct HW functions localized
along $z$. The Wannier spectrum is shown in Fig. 5(b). Its periodicity is
$c/2$ because the cell is doubled along $z$, and only one period is shown. The
spectrum is gapless, with two pairs of bands crossing in opposite directions,
between $\overline{\rm X}$ and $\overline{\Gamma}$, the gap centered at
$z=c/4$ (only one of the two crossings is shown). This spectral flow arises
from the nonzero MCN associated with $M_{y}$ symmetry (equivalent to $M_{z}$),
which leaves invariant the BZ plane containing the $\Gamma$, X, ${\rm R}_{2}$,
and ${\rm Y}_{2}$ points. For a discussion of such “in-plane” Wannier flow
associated with a nonzero MCN, see Ref. Gresch  _et al._ (2017).

Since $M_{z}$ is a type-2 mirror, we evaluate its unique MCN using Eq. (38).
And since the Wannier spectrum is gapless, and hence devoid of flat bands, we
set $\overline{C}_{\rm A}=0$ in that equation to obtain

$\mu_{\rm G}=W_{\rm A}\,,$ (53)

which says that the MCN equals the sum of the winding numbers of all the point
nodes on the $z=0$ plane.

As indicated in Fig. 5(d), there are 16 independent point nodes in total on
that plane, all of them simple nodes where only two bands meet. Seven have
winding numbers $+1$ and the other nine have winding numbers $-1$, yielding
$\mu_{\rm G}=-2$ for the MCN. This value is in agreement with that originally
obtained in Ref. Hsieh  _et al._ (2012) from a $k\cdot p$ analysis of the band
inversions. Using Eq. (39), we confirm that the system is axion-trivial.

###  VI.3  Modified Dirac model on a cubic lattice

In this section we study a 3D toy model constructed by first modifying the
free Dirac equation to enable topological phases for certain parameter values,
and then placing it on a cubic lattice. The 4$\times$4 Hamiltonian matrix in
reciprocal space reads Shen  _et al._ (2011); Rauch  _et al._ (2017)

$H({\bf k})=\left(\begin{matrix}m-2MK(\mathbf{k})&0&c\sin k_{z}&c(\sin
k_{x}-i\sin k_{y})\\\ 0&m-2MK(\mathbf{k})&c(\sin k_{x}+i\sin k_{y})&-c\sin
k_{z}\\\ c\sin k_{z}&c(\sin k_{x}-i\sin k_{y})&-m+2MK(\mathbf{k})&0\\\ c(\sin
k_{x}+i\sin k_{y})&-c\sin k_{z}&0&-m+2MK(\mathbf{k})\end{matrix}\right)\,,$
(54)

where $K(\mathbf{k})=3-\cos k_{x}-\cos k_{y}-\cos k_{z}$, and $c$, $m$, and
$M$ are dimensionless parameters inherited from the original isotropic
modified Dirac equation Shen  _et al._ (2011) by setting the rest mass
$m_{0}c^{2}$ to be the energy scale of the model Rauch  _et al._ (2017).

Figure 6: Topological phase diagram of the model of Eq. (54) for $c=1.0$.
Orange and blue regions denote axion-even ($\theta=0$) and axion-odd
($\theta=\pi$) phases, respectively.

The topological phase diagram of the half-filled model is shown in Fig. 6 for
$c=1.0$. The system is gapped except on the $m=0,4M,8M,12M$ lines, where the
gap closes at $\Gamma=(0,0,0)$, ${\rm X}=(\pi,0,0)$, ${\rm M}=(\pi,\pi,0)$,
and ${\rm A}=(\pi,\pi,\pi)$, respectively. As shown in Appendix C, those
metallic lines separate axion-trivial from axion-odd insulating phases.

The axion angle is quantized by several axion-odd symmetries. Some are
$z$-reversing (inversion and horizontal mirror $M_{z}$), and others are
$z$-preserving (spinful time reversal and vertical mirrrors). As $M_{z}$ is a
type-1 mirror, it protects two MCNs that are related to the axion angle by Eq.
(37).

####  VI.3.1 Axion-odd phase with protected Wannier flow

For our numerical tests we set $c=m=1.0$ and $M=0.5$ to put the model in the
axion-odd phase. The energy band structure is shown in Fig. 7(a). The bands
are pairwise degenerate due to the presence of time-reversal and inversion
symmetry, with a finite gap between the two pairs over the entire BZ. The
Fermi level is placed at midgap.

Figure 7: (a) Energy bands of the model described by Eq. (54) with $c=m=1.0$
and $M=0.5$. The bands are doubly degenerate, and the Fermi level (dashed
line) has been placed at midgap. (b) Wannier band structure obtained from the
valence states. (c) and (d) Heatmap plots of the gap function of Eq. (52)
about the $z=0$ and $z=c/2$ planes, respectively, with the nodal points color-
coded as in Fig. 3(c).

Since the system is axion-odd and has $z$-preserving axion-odd symmetries, the
connectivity (or “flow”) of the Wannier bands is topologically protected
Varnava  _et al._ (2020). In particular, spinful time reversal symmetry
requires that the two bands per vertical cell are glued together as follows:
one band touches the band above at one of the four time-reversal invariant
momenta (TRIM), and it touches the periodic image below at the other three. As
for the $z$-reversing axion-odd symmetries, the effect of $M_{z}$ is to pin
the up-touching to one of the mirror planes and the three down-touchings to
the other, while inversion further constrains the four touchings to occur at
TRIM on those planes, as already mandated by time reversal.

The pattern of band touchings described above is confirmed by Fig. 7(b), where
we plot the Wannier bands. They were obtained by placing at the origin the
four basis orbitals that belong to the home unit cell, and making the diagonal
approximation of Eq. (44) for the position matrix. There is one band touching
at $\overline{\Gamma}$ on the B plane, and three more on the A plane: one at
$\overline{\rm M}$, and the others at the two $\overline{\rm X}$ points.

Since the Wannier spectrum is gapless, the MCNs $\mu_{\rm G}$ and $\mu_{\rm
X}$ are given respectively by the half-sum and the half-difference of the net
winding numbers on the A and B planes [Eqs. (42) and (43)]. As indicated in
the gap-function plots of Figs. 7(c,d), the three nodes at A give $W_{\rm
A}=-1$ and the single node at B gives $W_{\rm B}=-1$, so that $\mu_{\rm G}=-1$
and $\mu_{\rm X}=0$. Note that $\mu_{\rm G}+\mu_{\rm X}$ is an odd number, as
required by Eq. (37) for an axion-odd system.

####  VI.3.2 Axion-odd phase with fragile Wannier flow

If the $z$-preserving axion-odd symmetries of the model (time reversal and
vertical mirrors) are weakly broken, the system will remain in an axion-odd
phase protected by $M_{z}$ and inversion. But since these are $z$-reversing
operations, the Wannier spectrum is no longer topologically required to be
gapless. The Wannier flow is only protected in a “fragile” sense, and it can
be destroyed, while preserving $M_{z}$, by adding some weakly-coupled trivial
bands to the valence manifold Varnava  _et al._ (2020); Wieder and Bernevig
(2018). Below we carry out this procedure in two different ways, and confirm
that the MCNs remain the same as in the original model.

##### Insertion of a symmetric pair of occupied orbitals

Figure 8: (a) Energy bands of the same model as in Fig. 7, after adding an
extra pair of occupied orbitals with $E=-4.0$ at $z=\pm 0.2c$ and coupling
them to the other orbitals. The bands are doubly degenerate, and the Fermi
level (dashed line) has been placed at midgap. (b) Wannier band structure
obtained from the valence states, with small gaps around $z=\pm 0.2c$ due to
the added orbitals.

Here we implement the strategy outlined in Sec. IV.2. We insert in the unit
cell two more orbitals, denoted as $|5\rangle$ and $|6\rangle$, that have
opposite spins and the same on-site energy $E=-4.0$. To break time reversal
and the vertical mirrors while preserving $M_{z}$ and inversion, we place the
spin-up orbital $|5\rangle$ at $(x,y,z)=(0.0,0.0,0.2c)$, and the spin-down
orbital $|6\rangle$ at $(x,y,z)=(0.0,0.0,-0.2c)$, keeping the original
orbitals $|1\rangle$ to $|4\rangle$ at the origin. Finally, we couple the new
orbitals to the old via the matrix elements $\langle 5|H|1\rangle=\langle
6|H|2\rangle=0.5$. The resulting model retains the $M_{z}$ and inversion
symmetries of the original model, and it breaks the time-reversal and vertical
mirror symmetries in the $Z$ matrix of Eq. (45) (but not in the Hamiltonian).

The energy and Wannier band structures are plotted in Figs. 8(a,b). Because
the Hamiltonian has both inversion and time-reveral symmetry, the energy bands
remain doubly degenerate as in Fig. 7(a). The breaking of the $z$-preserving
symmetries in the $Z$ matrix is reflected in the Wannier spectrum which is no
longer connected as in Fig. 7(b), with small gaps opening up near $z=\pm
0.2c$. The node at $\overline{\Gamma}$ on the B plane and those at
$\overline{\rm X}_{1}$, $\overline{\rm X}_{2}$, and $\overline{\rm M}$ on the
A plane remain intact, protected by $M_{z}$ and inversion. Their winding
numbers are also unchanged, leading to the same MCNs as in the original model.

##### Insertion of a single occupied orbital at $z=0$.

Figure 9: (a) Energy bands of the same model as in Fig. 7, after adding an
extra occupied orbital at $z=0$ and coupling it to the other orbitals. The
Fermi level (dashed line) has been placed in the gap. (b) Wannier band
structure obtained from the valence states. The added orbital generates a flat
band at $z=0$, which repels the nodal points on that plane (lower panel).

An alternative way of opening up a gap in the Wannier spectrum is to insert a
flat band on a mirror plane. To illustrate this procedure, we add at the
origin a single spin-up orbital $|5\rangle$ with on-site energy $E=-4.0$ and
odd parity about that plane, and couple it to the model via $\langle
5|H|1\rangle=\langle 5|H|4\rangle=2.0$. Because the orbital is spin-polarized,
it breaks time reversal; and because the spin points in the vertical
direction, it also breaks all vertical mirrors while preserving $M_{z}$. In
addition, the coupling terms break inversion symmetry, leaving $M_{z}$ as the
only axion-odd symmetry. The energy bands of the modified model are shown in
Fig. 9(a). A new band has appeared below the other four, so that there are now
three valence bands in total, leading to three Wannier bands.

The added orbital, which belongs to the ${\rm A}^{+}$ class in Table 1,
generates an extra even-parity state at both G and X. This creates an
imbalance $\Delta N_{\rm G}=\Delta N_{\rm X}=1$ between even- and odd-parity
states on the two mirror-invariant BZ planes, which according to Eq. (20)
results in a flat band at A. We emphasize that this extra band remains flat
even after the added orbital is coupled to the model, as long as the coupling
terms respect $M_{z}$ symmetry. As already mentioned, those terms are chosen
to break inversion symmetry. This is needed to ensure that the three point
nodes on the A plane are repelled by the flat band in the manner described in
Sec. III.2.2, since inversion symmetry would otherwise protect them.

The resulting Wannier bands are displayed in the upper panel of Fig. 9(b);
because of the lowered symmetry, the node at $z=c/2$ is no longer pinned to
$\overline{\Gamma}$ as in Fig. 7(b). The lower panel reveals a perfectly flat
band at $z=0$, well separated from a pair of dispersive bands whose three
touchings on the $z=0$ plane in Fig. 7(c) have been gapped out. Under these
circumstances, Eqs. (33) and (34) for the MCNs reduce to

$\mu_{\rm G}=\tfrac{1}{2}(p_{\rm A}\overline{C}_{\rm A}+W_{\rm B})$ (55)

and

$\mu_{\rm X}=\tfrac{1}{2}(p_{\rm A}\overline{C}_{\rm A}-W_{\rm B})\,.$ (56)

The single node at B has the same winding number $W_{\rm B}=-1$ as in the
original model, while the net winding number $W_{\rm A}=-1$ of the gapped-out
nodes at A has been transferred to the index $p_{\rm A}\overline{C}_{\rm A}$
of the flat band ($p_{\rm A}=-1$, and $\overline{C}_{\rm A}=+1$). Overall, the
MCNs remain unchanged.

##  VII Summary

In summary, we have investigated the topological properties of mirror-
symmetric insulating crystals from the viewpoint of HW functions localized
along the direction orthogonal to the mirror plane. We first clarified the
generic behaviors of the associated Wannier bands, and then derived a set of
rules for deducing the MCNs. To validate and illustrate the formalism, we
applied it to SnTe in the monolayer and bulk forms, and to a toy model of an
axion-odd insulator.

In the HW representation, the MCNs are expressed in terms of a set of integer-
valued properties of the Wannier bands on the mirror planes: the Chern numbers
and mirror parities of flat bands lying on those planes, and the winding
numbers of the touching points on those planes between symmetric pairs of
dispersive bands. One advantage of this representation is that it reveals the
relation between the MCNs and the axion $\mathbbm{Z}_{2}$ index from purely
bulk considerations. That relation is far from obvious in the standard Bloch
representation, and previously it had only been obtained via an indirect
argument involving surface states.

In some cases the axion $\mathbbm{Z}_{2}$ index can be determined by visual
inspection of the Wannier band structure, e.g., by counting the number of
nodal points between certain bands Varnava  _et al._ (2020). We have found
that mere visual inspection does not suffice for obtaining the MCNs since it
does not reveal, for example, the relative signs of the winding numbers of
different nodes.

Interestingly, in certain cases where flat Wannier bands are present the
magnitudes of the MCN can be determined without having to divide the occupied
manifold into two mirror sectors. This follows from the uniform-parity
assumption for the flat bands, which has no counterpart in the Bloch
representation. Since the determination of the mirror parities is the most
cumbersome step in the calculation of MCNs, this feature of the HW formalism
could lead to a more automated algorithm for computing MCNs. Even without such
further developments, the formalism has already proven useful for discussing
the topological classification of mirror-symmetric insulators.

###### Acknowledgements.

Work by T.R. was supported by the Deutsche Forschungsgemeinschaft Grant No. Ra
3025/1-1 from the Deutsche Forschungsgemeinschaft. Work by D.V. was supported
by National Science Foundation Grant DMR-1954856. Work by I.S. was supported
by Grant No. FIS2016-77188-P from the Spanish Ministerio de Economía y
Competitividad.

##  Appendix A Derivation of Eqs. (20-22)

According to Table 1, the numbers of occupied states with each mirror parity
at G and X are

$\displaystyle N_{{\rm G}^{\pm}}$ $\displaystyle=\overline{N}_{{\rm
A}^{\pm}}+\overline{N}_{{\rm B}^{\pm}}+\frac{1}{2}\widetilde{N}\,,$ (57a)
$\displaystyle N_{{\rm X}^{\pm}}$ $\displaystyle=\overline{N}_{{\rm
A}^{\pm}}+\overline{N}_{{\rm B}^{\mp}}+\frac{1}{2}\widetilde{N}\,,$ (57b)

where $\widetilde{N}=\widetilde{N}_{\rm A}+\widetilde{N}_{\rm
B}+\widetilde{N}_{\rm UC}$ is the total number of dispersive Wannier bands per
cell. Letting $\Delta N_{{\rm G}}=N_{{\rm G}^{+}}-N_{{\rm G}^{-}}$ and
$\Delta\overline{N}_{{\rm A}}=\overline{N}_{{\rm A}^{+}}-\overline{N}_{{\rm
A}^{-}}$, and defining $\Delta N_{{\rm X}}$ and $\Delta N_{{\rm B}}$ in the
same way, we find

$\displaystyle\Delta\overline{N}_{{\rm A}}$
$\displaystyle=\frac{1}{2}\left(\Delta N_{{\rm G}}+\Delta N_{{\rm
X}}\right)\,,$ (58a) $\displaystyle\Delta\overline{N}_{{\rm B}}$
$\displaystyle=\frac{1}{2}\left(\Delta N_{{\rm G}}-\Delta N_{{\rm
X}}\right)\,.$ (58b)

Under the uniform parity assumption $|\Delta\overline{N}_{{\rm
A}}|=\overline{N}_{\rm A}$ and $|\Delta\overline{N}_{{\rm
B}}|=\overline{N}_{\rm B}$, resulting in Eqs. (20) and (21). In the case of a
type-2 mirror A and B are equivalent, and from Eq. (57a)
$\Delta\overline{N}_{\rm A}+\Delta\overline{N}_{\rm B}=\Delta N_{\rm G}$.
Hence $\Delta\overline{N}_{\rm A}=\Delta\overline{N}_{\rm B}=\Delta N_{\rm
G}/2$, yielding Eq. (22) under the same assumption.

##  Appendix B Derivation of Eq. (27)

Let us prove Eq. (27) for the case of a single pair of dispersive Wannier
bands connected by point nodes on the A plane. In this case the matrix
$f_{\bm{\kappa}}$ of Eq. (49) reduces to the scalar

$f_{\bm{\kappa}}\equiv\langle\widetilde{h}_{\bm{\kappa}}^{+}|z|\widetilde{h}_{\bm{\kappa}}^{-}\rangle=|f_{\bm{\kappa}}|e^{i\gamma_{\bm{\kappa}}}\,,$
(59)

where $|\widetilde{h}_{\bm{\kappa}}^{\pm}\rangle$ are states of even or odd
mirror parity constructed from the pair of HW functions as described in Sec.
V.3.1. These states are cell-periodic in plane and localized along $z$, and we
also define new states
$|\psi_{\bm{\kappa}}^{\pm}\rangle=e^{i{\bm{\kappa}}\cdot{\bf
r}}|\widetilde{h}_{\bm{\kappa}}^{\pm}\rangle$ that are Wannier-like along $z$
and Bloch-like in plane.

When the Chern numbers $\widetilde{C}_{{\rm A}^{\pm}}$ are nonzero, it becomes
impossible to choose a gauge for the states $|\psi_{\bm{\kappa}}^{\pm}\rangle$
that is both smooth and periodic in the projected 2D BZ Vanderbilt (2018). We
assume a square BZ with $k_{x},k_{y}\in[0,2\pi]$, and choose a smooth but
nonperiodic gauge for the $|\psi_{\bm{\kappa}}^{-}\rangle$ states. To
characterize the lack of periodicity, let the phase relations between the
edges of the BZ be

$|\psi^{-}_{\rm R}\rangle=e^{-i\mu}|\psi^{-}_{\rm
L}\rangle\,,\quad|\psi^{-}_{\rm T}\rangle=e^{-i\nu}|\psi^{-}_{\rm
B}\rangle\,,$ (60)

where $\\{\text{L,R,T,B}\\}=\\{\text{left,right,top,bottom}\\}$,
$\mu=\mu(k_{y})$, and $\nu=\nu(k_{x})$. Also let

$\Delta\mu=\mu(2\pi)-\mu(0)\,,\quad\Delta\nu=\nu(2\pi)-\nu(0)\,.$ (61)

When computing the Berry phase around the BZ boundary as an integral of the
connection ${\bf
A}_{\bm{\kappa}}^{-}=i\langle\widetilde{h}_{\bm{\kappa}}^{-}|\partial_{\bm{\kappa}}\widetilde{h}_{\bm{\kappa}}^{-}\rangle$,

$\phi_{-}=\oint_{\partial\text{BZ}}{\bf A}_{\bm{\kappa}}^{-}\cdot
d{\bm{\kappa}}\,,$ (62)

the contribution from the L and R segments cancel except for terms coming from
$\mu$, and similarly for the top and bottom segments. It follows that

$\phi_{-}=\Delta\mu-\Delta\nu\,.$ (63)

We assume a smooth but nonperiodic gauge for the
$|\psi^{+}_{\bm{\kappa}}\rangle$ states as well, so that the phase
$\gamma_{\bm{\kappa}}$ in Eq. (59) becomes a smooth function of
${\bm{\kappa}}$ (except at the nodes, where $f_{\bm{\kappa}}$ vanishes and
$\gamma_{\bm{\kappa}}$ becomes ill defined). Now we phase-align
$|\psi^{+}_{\bm{\kappa}}\rangle$ with $|\psi^{-}_{\bm{\kappa}}\rangle$ by re-
gauging as follows,

$|\psi^{+}_{\bm{\kappa}}\rangle^{\prime}=e^{i\gamma_{\bm{\kappa}}}|\psi^{+}_{\bm{\kappa}}\rangle\,.$
(64)

(In this new gauge $f^{\prime}_{\bm{\kappa}}$ is real, and
$\gamma^{\prime}_{\bm{\kappa}}$ is zero everywhere.) This will make a gauge
for $|\psi^{+}_{\bm{\kappa}}\rangle^{\prime}$ that is also nonperiodic. For
the moment we only assume that this gauge is smooth in a neighborhood
extending some small distance inside the boundary; we ignore what is going on
deeper inside. It is not hard to see that the same relations as in Eq. (60),
with the same functions $\mu$ and $\nu$, apply to the
$|\psi^{+}_{\bm{\kappa}}\rangle^{\prime}$ states, and it follows that

$\phi^{\prime}_{+}=\phi_{-}\quad\text{(call it $\phi$)\,.}$ (65)

Now, in the case of the $|\psi^{-}_{\bm{\kappa}}\rangle$ states the interior
was smooth, so by applying Stokes’ theorem to

$2\pi\widetilde{C}_{{\rm A}^{-}}=\int_{\rm
BZ}\Omega_{\bm{\kappa}}^{-}\,d^{2}k$ (66)

where
$\Omega^{-}_{\bm{\kappa}}=\partial_{k_{x}}A^{-}_{{\bm{\kappa}},y}-\partial_{k_{y}}A^{-}_{{\bm{\kappa}},x}$
is the Berry curvature of state $|u^{-}_{\bm{\kappa}}\rangle$, we get

$2\pi\widetilde{C}_{{\rm A}^{-}}=\phi\,.$ (67)

If the interior of $|\psi^{+}_{\bm{\kappa}}\rangle^{\prime}$ were also smooth,
we would conclude that $\widetilde{C}_{{\rm A}^{+}}=\widetilde{C}_{{\rm
A}^{-}}$. Conversely, when the MCN is nonzero there must exist nonanalytic
points where the phase of $|u^{+}_{\bm{\kappa}}\rangle^{\prime}$ changes
discontinuously. Those points are precisely the nodes of $f_{\bm{\kappa}}$,
which we label by $j$; they act as vortex singularities of the Berry
connection

$\left({\bf A}_{\bm{\kappa}}^{+}\right)^{\prime}={\bf
A}_{\bm{\kappa}}^{+}-\partial_{\bm{\kappa}}\gamma_{\bm{\kappa}}\,,$ (68)

and we extract their winding numbers $W_{j}$ using Eq. (28). Let $S$ be the
interior of the projected BZ with a small circle $c_{j}$ cut around each node,
and apply Stokes’ theorem over the region $S$ to find

$\int_{S}\Omega_{\bm{\kappa}}^{+}\,d^{2}k=\int_{\partial\text{BZ}}\left({\bf
A}_{\bm{\kappa}}^{+}\right)^{\prime}\cdot
d{\bm{\kappa}}-\sum_{j}\oint_{c_{j}}\left({\bf
A}_{\bm{\kappa}}^{+}\right)^{\prime}\cdot d{\bm{\kappa}}\,.$ (69)

The first term on the right-hand side is equal to
$\phi^{\prime}_{+}=\phi=2\pi\widetilde{C}_{{\rm A}^{-}}$. In the limit of
small circles the left-hand side becomes $2\pi\widetilde{C}_{{\rm A}^{+}}$,
and the second term on the right-hand side reduces to $2\pi\sum_{j}\,W_{j}$
(this follows from Eq. (68) by noting that ${\bf A}^{+}_{\bm{\kappa}}$ is
smooth everywhere). Thus $\widetilde{C}_{{\rm A}^{+}}-\widetilde{C}_{{\rm
A}^{-}}$ equals $W_{\rm A}=\sum_{j\in{\rm A}}\,W_{j}$, which is what we set
out to prove. The same result holds if more than one pair of bands meet at
some of the point nodes, in which case $\gamma_{\bm{\kappa}}$ is given by the
more general expression in Eq. (50).

##  Appendix C Phase diagram of the modified Dirac model on a cubic lattice

Figure 10: Wannier bands of the modified Dirac model on a cubic lattice [Eq.
(54)], for $m=1.0$ and varying $M$.

In this Appendix, we map out the topological phase diagram of the model of Eq.
(54) as a function of the parameters $m$ and $M$, for $c=1.0$. The band gap
closes for $m=0,4M,8M,12M$ at the points $\Gamma$, $\rm X$, $\rm M$, and $\rm
A$, respectively Shen (2012). Those lines in the phase diagram mark the
topological phase transitions between axion-even and axion-odd phases.

To decide which phases are trivial and which are topological, it is sufficient
to inspect the Wannier band structures in Fig. 10, obtained for representative
states in each of the four phases along the $m=1.0$ line. Since the model has
several axion-odd symmetries (time reversal, inversion, and multiple mirrors),
we can base our analysis on either of them, applying in each case the rules
given in Ref. Varnava  _et al._ (2020) to determine the axion
$\mathbbm{Z}_{2}$ index. In the following, we choose to focus on time-reversal
symmetry.

The Wannier spectrum of an axion-odd phase with spinful time-reversal symmetry
must be gapless, with each band touching the band above at one of the four
TRIM and the band below at the other three (or vice-versa). From this
criterion we conclude that Figs. 10(a,c) correspond to axion-trivial phases,
and Figs. 10(b,d) to axion-odd topological phases. Hence the system is
topological for $0<m/M<4$ and $8<m/M<12$, producing the phase diagram in Fig.
6. This is in agreement with Ref. Shen, 2012, where the strong topological
index $\nu_{0}=\theta/\pi$ of each phase was determined from the parity
eigenvalues of the Bloch states at the eight TRIM in the 3D BZ Fu and Kane
(2007).

## References

  * Teo  _et al._ (2008) J. C. Y. Teo, L. Fu, and C. L. Kane, “Surface states and topological invariants in three-dimensional topological insulators: Application to ${\text{Bi}}_{1-x}{\text{Sb}}_{x}$,” Phys. Rev. B 78, 045426 (2008).
  * Ando and Fu (2015) Y. Ando and L. Fu, “Topological Crystalline Insulators and Topological Superconductors: From Concepts to Materials,” Annu. Rev. Condens. Matter Phys. 6, 361 (2015).
  * Qi  _et al._ (2008) X.-L. Qi, T. L. Hughes, and S.-C. Zhang, “Topological field theory of time-reversal invariant insulators,” Phys. Rev. B 78, 195424 (2008).
  * Essin  _et al._ (2009) A. M. Essin, J. E. Moore, and D. Vanderbilt, “Magnetoelectric Polarizability and Axion Electrodynamics in Crystalline Insulators,” Phys. Rev. Lett. 102, 146805 (2009).
  * Vanderbilt (2018) D. Vanderbilt,  _Berry Phases in Electronic Structure Theory: Electric Polarization, Orbital Magnetization and Topological Insulators_ (Cambridge University Press, Cambridge (United Kingdom), 2018).
  * Armitage and Wu (2019) N. P. Armitage and Liang Wu, “On the matter of topological insulators as magnetoelectrics,” SciPost Phys. 6, 46 (2019).
  * Nenno  _et al._ (2020) D. M. Nenno, C. A. C. Garcia, J. Gooth, C. Felser, and P. Narang, “Axion physics in condensed-matter systems,” Nat. Rev. Phys. 2, 682 (2020).
  * Sekine and Nomura (2021) A. Sekine and K. Nomura, “Axion electrodynamics in topological materials,” J. Appl. Phys. 129, 141101 (2021).
  * Otrokov  _et al._ (2019) M. M. Otrokov _et al._ , “Prediction and observation of an antiferromagnetic topological insulator,” Nature 576, 416 (2019).
  * Mong  _et al._ (2010) R. S. K. Mong, A. M. Essin, and J. E. Moore, “Antiferromagnetic topological insulators,” Phys. Rev. B 81, 245209 (2010).
  * Fu and Kane (2007) L. Fu and C. L. Kane, “Topological insulators with inversion symmetry,” Phys. Rev. B 76, 045302 (2007).
  * Turner  _et al._ (2012) A. M. Turner, Y. Zhang, R. S. K. Mong, and A. Vishwanath, “Quantized response and topology of magnetic insulators with inversion symmetry,” Phys. Rev. B 85, 165120 (2012).
  * Varnava  _et al._ (2020) N. Varnava, I. Souza, and D. Vanderbilt, “Axion coupling in the hybrid Wannier representation,” Phys. Rev. B 101, 155130 (2020).
  * Varjas  _et al._ (2015) D. Varjas, F. de Juan, and Y.-M. Lu, “Bulk invariants and topological response in insulators and superconductors with nonsymmorphic symmetries,” Phys. Rev. B 92, 195116 (2015).
  * Fulga  _et al._ (2016) I. C. Fulga, N. Avraham, H. Beidenkopf, and A. Stern, “Coupled-layer description of topological crystalline insulators,” Phys. Rev. B 94, 125405 (2016).
  * Hsieh  _et al._ (2012) T. H. Hsieh, H. Lin, J. Liu, W. Duan, A. Bansil, and L. Fu, “Topological crystalline insulators in the SnTe material class,” Nat. Commun. 3, 982 (2012).
  * Liu  _et al._ (2014) J. Liu, T. H. Hsieh, P. Wei, W. Duan, J. Moodera, and L. Fu, “Spin-filtered edge states with an electrically tunable gap in a two-dimensional topological crystalline insulator,” Nature Mater. 13, 178 (2014).
  * Marzari and Vanderbilt (1997) N. Marzari and D. Vanderbilt, “Maximally localized generalized Wannier functions for composite energy bands,” Phys. Rev. B 56, 12847 (1997).
  * Taherinejad and Vanderbilt (2015) M. Taherinejad and D. Vanderbilt, “Adiabatic Pumping of Chern-Simons Axion Coupling,” Phys. Rev. Lett. 114, 096401 (2015).
  * Sutherland (1986) B. Sutherland, “Localization of electronic wave functions due to local topology,” Phys. Rev. B 34, 5208 (1986).
  * Lieb (1989) E. H. Lieb, “Two theorems on the Hubbard model,” Phys. Rev. Lett. 62, 1201 (1989).
  * Ramachandran  _et al._ (2017) A. Ramachandran, A. Andreanov, and S. Flach, “Chiral flat bands: Existence, engineering, and stability,” Phys. Rev. B 96, 161104(R) (2017).
  * Asbóth  _et al._ (2016) J. A. Asbóth, A. Pályi, and L. Oroszlány,  _A Short Course on Topological Insulators_ (Springer, Cham, 2016).
  * Kim  _et al._ (2015) Y. Kim, C. L. Kane, E. J. Mele, and A. M. Rappe, “Layered Topological Crystalline Insulators,” Phys. Rev. Lett. 115, 086802 (2015).
  * Wieder and Bernevig (2018) B. J. Wieder and B. A. Bernevig, “The Axion Insulator as a Pump of Fragile Topology,”  (2018), arXiv:1810.02373 .
  * Park and Marzari (2011) C.-H. Park and N. Marzari, “Berry phase and pseudospin winding number in bilayer graphene,” Phys. Rev. B 84, 205440 (2011).
  * (27) The PythTB code package is available at http://www.physics.rutgers.edu/pythtb/about.html.
  * Enkovaara  _et al._ (2010) J. Enkovaara, C. Rostgaard, J. J. Mortensen, J. Chen, M. Dułak, L. Ferrighi, J. Gavnholt, C. Glinsvad, V. Haikola, H. A. Hansen, H. H. Kristoffersen, M. Kuisma, A. H. Larsen, L. Lehtovaara, M. Ljungberg, O. Lopez-Acevedo, P. G. Moses, J. Ojanen, T. Olsen, V. Petzold, N. A. Romero, J. Stausholm-Møller, M. Strange, G. A. Tritsaris, M. Vanin, M. Walter, B. Hammer, H. Häkkinen, G. K. H. Madsen, R. M. Nieminen, J. K. Nørskov, M. Puska, T. T. Rantala, J. Schiøtz, K. S. Thygesen, and K. W. Jacobsen, “Electronic structure calculations with GPAW: a real-space implementation of the projector augmented-wave method,” J. Phys. Condens. Matter 22, 253202 (2010).
  * Mostofi  _et al._ (2014) A. A. Mostofi, J. R. Yates, G. Pizzi, Y.-S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, “An updated version of wannier90: A tool for obtaining maximally-localised Wannier functions,” Comput. Phys. Commun. 185, 2309 (2014).
  * Olsen (2016) T. Olsen, “Designing in-plane heterostructures of quantum spin Hall insulators from first principles: $1\text{T}^{\prime}$-MoS2 with adsorbates,” Phys. Rev. B 94, 235106 (2016).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, “Generalized Gradient Approximation Made Simple,” Phys. Rev. Lett. 77, 3865 (1996).
  * Perdew  _et al._ (1997) J. P. Perdew, K. Burke, and M. Ernzerhof, “Generalized Gradient Approximation Made Simple [Phys. Rev. Lett. 77, 3865 (1996)],” Phys. Rev. Lett. 78, 1396(E) (1997).
  * Blöchl (1994) P. E. Blöchl, “Projector augmented-wave method,” Phys. Rev. B 50, 17953 (1994).
  * Souza  _et al._ (2001) I. Souza, N. Marzari, and D. Vanderbilt, “Maximally localized Wannier functions for entangled energy bands,” Phys. Rev. B 65, 035109 (2001).
  * Olsen  _et al._ (2019) T. Olsen, E. Andersen, T. Okugawa, D. Torelli, T. Deilmann, and K. S. Thygesen, “Discovering two-dimensional topological insulators from high-throughput computations,” Phys. Rev. Mater. 3, 024005 (2019).
  * Haastrup  _et al._ (2018) Sten Haastrup, Mikkel Strange, Mohnish Pandey, Thorsten Deilmann, Per S. Schmidt, Nicki F. Hinsche, Morten N. Gjerding, Daniele Torelli, Peter M. Larsen, Anders C. Riis-Jensen, Jakob Gath, Karsten W. Jacobsen, Jens Jørgen Mortensen, Thomas Olsen, and Kristian S. Thygesen, “The Computational 2D Materials Database: High-throughput modeling and discovery of atomically thin crystals,” 2D Mater. 5, 042002 (2018).
  * Kobayashi (2015) K. Kobayashi, “Electronic states of SnTe and PbTe (001) monolayers with supports,” Surf. Sci. 639, 54 (2015).
  * Liu  _et al._ (2015) J. Liu, X. Qian, and L. Fu, “Crystal Field Effect Induced Topological Crystalline Insulators In Monolayer IV–VI Semiconductors,” Nano Lett. 15, 2657 (2015).
  * Tanaka  _et al._ (2012) Y. Tanaka, Z. Ren, T. Sato, K. Nakayama, S. Souma, T. Takahashi, K. Segawa, and Y. Ando, “Experimental realization of a topological crystalline insulator in SnTe,” Nat. Phys. 8, 800 (2012).
  * Gresch  _et al._ (2017) D. Gresch, G. Autès, O. V. Yazyev, M. Troyer, D. Vanderbilt, B. A. Bernevig, and A. A. Soluyanov, “Z2Pack: Numerical implementation of hybrid Wannier centers for identifying topological materials,” Phys. Rev. B 95, 075146 (2017).
  * Shen  _et al._ (2011) S.-Q. Shen, W.-Y. Shan, and H.-Z. Lu, “Topological Insulator and the Dirac Equation,” SPIN 01, 33 (2011).
  * Rauch  _et al._ (2017) T. Rauch, H. Nguyen Minh, J. Henk, and I. Mertig, “Model for ferromagnetic Weyl and nodal line semimetals: Topological invariants, surface states, anomalous and spin Hall effect,” Phys. Rev. B 96, 235103 (2017).
  * Shen (2012) S.-Q. Shen,  _Topological Insulators – Dirac Equation in Condensed Matters_ (Springer, Berlin, Heidelberg, 2012).