#  Quantized surface magnetism and higher-order topology:
Application to the Hopf insulator

Penghao Zhu  Department of Physics and Institute for Condensed Matter Theory,
University of Illinois at Urbana-Champaign, Urbana, Illinois 61801, USA
Taylor L. Hughes  Department of Physics and Institute for Condensed Matter
Theory, University of Illinois at Urbana-Champaign, Urbana, Illinois 61801,
USA  A. Alexandradinata  Corresponding author<EMAIL_ADDRESS>Department of Physics and Institute for Condensed Matter Theory, University of
Illinois at Urbana-Champaign, Urbana, Illinois 61801, USA

###### Abstract

We identify topological aspects of the subextensive magnetic moment
contributed by the surfaces of a three-dimensional crystallite – assumed to be
insulating in the bulk as well as on all surface facets, with trivial Chern
invariants in the bulk. The geometric component of this subextensive moment is
given by its derivative with respect to the chemical potential, at zero
temperature and zero field, per unit surface area, and hence corresponds to
the surface magnetic compressibility. The sum of the surface compressibilities
contributed by two opposite facets of a cube-shaped crystallite is quantized
to an integer multiple of the fundamental constant $e/hc$; this integer is in
one-to-one correspondence with the net chirality of hinge modes on the surface
of the crystallite, manifesting a link with higher-order topology. The
contribution by a single facet to the magnetic compressibility need not be
quantized to integers; however, symmetry and/or Hilbert-space constraints can
fix the single-facet compressibility to half-integer multiples of $e/hc$, as
will be exemplified by the Hopf insulator.

††preprint: APS/123-QED

Disentangling bulk and boundary effects in orbital magnetism [1, 2, 3] has
long been a subtle affair. Boundary currents contribute extensively to the
orbital magnetic moment, yet the boundary current cannot be completely
attributed to single-particle states localized to the boundary [4, 5]. The
recent formulation of a local marker for the orbital magnetization allows one,
for the first time, to unambiguously distinguish boundary contributions to the
orbital magnetic moment from bulk contributions [6]. For band insulators, it
was subsequently realized that the boundary contribution can be extensive only
if the insulator has a nonzero Hall conductivity[5].

What is missed in previous works is the possibility of boundary contributions
to the magnetic moment which scale subextensively with the area of the surface
in the thermodynamic limit. This work aims to identify topological aspects of
the subextensive magnetic moment contributed by surfaces of a three-
dimensional crystallite – with trivial Chern invariants in the bulk. Such
aspects arise when there is a spectral gap to excitations not just in the
bulk, but also on all surface facets, as is assumed throughout this work.

For illustration, consider the magnetic moment contributed by the top and
bottom facets of a cubic slab depicted in Fig. 1. The geometric component of
this subextensive moment, as we will show, is given by its derivative with
respect to the chemical potential, at zero temperature and zero field, per
unit surface area. This quantity, henceforth referred to as the magnetic
compressibility, is known to be quantized to integer multiples of the
fundamental constant $e/hc$ for strictly two-dimensional, bulk insulators, and
the integer invariant has been identified with the net chirality of the edge
modes [2, 6, 7, 5]. Generalizing this statement to a three-dimensional cubic
slab, the analogous integer invariant becomes the net chirality of hinge modes
localized to the one-dimensional intersections of the four side facets, as
illustrated in Fig. 2. This means that a three-dimensional insulator with a
nontrivial magnetic compressibility is also a higher-order topological
insulator [8, 9, 10, 11, 12, 13, 14], as we will establish.

Figure 1: Schematic illustration of a cubic slab with surface magnetism
(indicated by current loops) and hinge modes (purple arrows). On the right, a
cross section of the left slab illustrates the exponential decay of the
surface magnetization.

This is one of several links between surface magnetism and higher-order
topology that we will explore in Sec. I. Unlike strictly two-dimensional
systems, the contribution of a single facet to the magnetic compressibility
need not be quantized to integers. However, we find that certain symmetries,
and/or Hilbert-space constraints, can fix the single-facet compressibility to
half-integer multiples of $e/hc$. Our main case study [cf. Sec. II] for this
half-integer quantization is the Hopf insulator [15, 16], a three-dimensional
magnetic topological insulator inspired by the Hopf map in differential
topology [17].

Beyond the Hopf insulator, we also consider a similar (but strictly distinct)
class of insulators with zero bulk magnetization at zero field, but with
nontrivially-quantized, single-facet magnetic compressibility. Such a phase of
matter does not fit into standard classifications under the family tree of
magnetism [18] (paramagnetism, ferromagnetism, ferrimagnetism, etc.), and we
refer to it as a topomagnetic insulator. In Sec. III, we show that some (but
not all) higher-order topological insulators known in the literature[19, 20,
21, 10] are topomagnetic insulators. The concluding Sec. IV contains a partial
summary of our results and suggestions on how the subextensive magnetic moment
can be measured.

##  I Faceted orbital magnetization and higher-order topology

After a brief review of the local marker for orbital magnetization in Sec.
I.1, we present a formula for the layer-resolved magnetization in Sec. I.2,
which is used to isolate a single facet’s contribution to the orbital magnetic
moment. The associated compressibility of this moment is shown to relate to
the geometric component of the surface anomalous Hall conductivity in Sec.
I.3. We also derive therein the correspondence between compressibility of
facets and chiral hinge modes, by arguments involving Amperean loops and the
Streda formula.

###  I.1 Review of local marker for orbital magnetization

For an insulating crystallite in contact with an electron reservoir (at
chemical potential $\mu$), the zero-temperature magnetization is related to
the grand thermodynamic potential $\Omega$ as

$\displaystyle{M=-\frac{1}{V}\frac{\partial\Omega}{\partial
B}\bigg{|}_{\mu}=-\frac{1}{V}\frac{\partial E}{\partial B}+\mu\frac{\partial
n}{\partial B}:=m_{1}+m_{2},}$ (1)

where $n=\frac{N}{V}$ is the electron density, $E$ is the internal energy, and
$B$ is the external magnetic field. The $m_{1}$ term contains an orbital
contribution originating from the minimal coupling to the magnetic field, as
well as a spin contribution originating from the Zeeman coupling to the
intrinsic spin magnetic moment. We refer to the orbital component of $m_{1}$
plus the entirety of $m_{2}$ as the orbital magnetization $m_{orb}$.

We shall primarily be concerned with the zero-$B$ orbital magnetization of
band insulators with a spectral gap in both the bulk and all surface facets of
a three-dimensional crystallite. Since we are not considering a fully periodic
system with translation symmetry it is important that we introduce a method to
calculate the magnetization in position-space. Fortunately, it was already
shown in Ref 6 that $m_{orb}$ can be expressed as an integral over continuous
position-space of a local marker $\mathfrak{m}$ for orbital magnetization. For
convenience, we present here the tight-binding analog of their formula:

$\displaystyle
M_{\gamma}=\frac{1}{V}\sum_{\bm{R}}\mathfrak{m}_{\gamma}(\bm{R}),$ (2)
$\displaystyle\mathfrak{m}_{\gamma}(\bm{R})$ $\displaystyle=\frac{e}{2\hbar
c}\epsilon_{\alpha\beta\gamma}\operatorname{Im}\sum_{a}\langle\bm{R},a|\bar{P}\hat{r}_{\alpha}\bar{Q}\bar{H}\bar{Q}\hat{r}_{\beta}\bar{P}$
$\displaystyle\quad-\bar{Q}\hat{r}_{\alpha}\bar{P}(\bar{H}-2\mu)\bar{P}\hat{r}_{\beta}\bar{Q}|\bm{R},a\rangle,$

where $\bar{H}$ is the Hamiltonian, and $\bar{P}$ and $\bar{Q}$ are the
projectors onto the occupied and unoccupied subspaces of $\bar{H}$,
respectively. $\bar{H},\bar{P}$ and $\bar{Q}$ act in a Hilbert space spanned
by the orthonormal basis vectors $\\{|\bm{R},a\rangle\\}_{\bm{R},a}$, where
$|\bm{R},a\rangle$ is a ket state of an electron localized (as a Wannier
function) to a Bravais-lattice vector $\bm{R}=(R_{x},R_{y},R_{z})$, and the
index $a$ labels the basis vectors of the reduced Hilbert space of one
primitive unit cell. The action of the position operator $\hat{r}_{\alpha}$ on
these states satisfies
$\hat{r}_{\alpha}|\bm{R},a\rangle=R_{\alpha}|\bm{R},a\rangle$. By assumption,
an energy gap separates the occupied and unoccupied subspaces of $\bar{H}$, so
that the matrix elements of $\bar{P}$ decay exponentially with increasing
position-space separation [22, 23]. This gives an exponentially sharp
distinction between the local marker $\mathfrak{m}_{\gamma}(\bm{R})$ evaluated
on the boundary versus in the bulk [6]. Using this formalism we immediately
observe that the $\mu$-dependent term in Eq. (2) is geometric, i.e., it is
simply proportional to the local Chern marker [24]. Alternatively, this
proportionality can be derived by applying the Streda formula [25, 26] to
$m_{2}=\mu\partial n/\partial B$.

###  I.2 Faceted orbital magnetization

To model a facet on a three-dimensional crystallite, let us consider a half-
infinite slab with surface normal vector $+\hat{z}$. The orbital magnetic
moment of the slab can be decomposed into a sum over layers indexed by
$l=1,2,\ldots$, where $l=1$ lies closest to the surface termination. By
summing the local marker over each layer, we show in Appendix A that the
orbital magnetization (contributed by each layer) is expressible as an
integral over the reduced Brillouin zone
($rBZ\ni\bm{k}_{\perp}=(k_{x},k_{y})$):

$\displaystyle
M_{z}(l):=\frac{e}{c\hbar}\operatorname{Im}\int_{rBZ}\frac{d\mathbf{k}_{\perp}}{(2\pi)^{2}}\bigg{[}(g_{\mathbf{k}_{\perp},xy}(l)$
(3) $\displaystyle+h_{\mathbf{k}_{\perp},xy}(l))-2\mu
f_{\mathbf{k}_{\perp},xy}(l))\bigg{]},$
$\begin{array}[]{c}g_{\mathbf{k}_{\perp},xy}(l)=\operatorname{Tr}_{\mathrm{cell},z}\left[\mathcal{P}_{l}\partial_{k_{x}}\widetilde{P}\widetilde{Q}\widetilde{H}\widetilde{Q}\partial_{k_{y}}\widetilde{P}\right],\\\
h_{\mathbf{k}_{\perp},xy}(l)=\operatorname{Tr}_{\mathrm{cell},z}\left[\mathcal{P}_{l}\widetilde{Q}\partial_{k_{y}}\widetilde{P}\widetilde{H}\partial_{k_{x}}\widetilde{P}\widetilde{Q}\right],\\\
f_{\mathbf{k}_{\perp},xy}(l)=\operatorname{Tr}_{\mathrm{cell},z}\left[\mathcal{P}_{l}\widetilde{Q}\partial_{k_{y}}\widetilde{P}\partial_{k_{x}}\widetilde{P}\widetilde{Q}\right],\end{array}$
(4)

where $g,h,f$ correspond, in order, to the three terms in Eq. (2). The
quantity $\widetilde{P}(\mathbf{k}_{\perp})$
($\widetilde{Q}(\mathbf{k}_{\perp})$) is the projector to the occupied
(unoccupied) subspace of the Hamiltonian $\widetilde{H}(\mathbf{k}_{\perp})$;
note the arguments have been suppressed in the above equations above. In a
slab geometry, these operators have matrix elements labelled as, for example,
$\widetilde{H}(\mathbf{k}_{\perp})_{R_{z}R_{z^{\prime}},ab}$, and
$``\operatorname{Tr}_{\mathrm{cell},z}$” is the trace with respect to the
composite index $(R_{z},a)$. Finally, the operator
$\mathcal{P}_{l}(\mathbf{k}_{\perp})$ is the projector onto the $R_{z}=l$
layer in the $\mathbf{k}_{\perp}$ momentum sector. In direct correspondence
with the local Chern marker [24], the imaginary part of $f$ in Eq. (3) is the
layer-resolved Berry curvature [27], and its integral over the $rBZ$ gives the
contribution (of each layer) to the Hall conductivity:

$C_{+\hat{z}}(l)=\frac{1}{2\pi}\int_{rBZ}d\bm{k}_{\perp}\operatorname{Im}f_{\mathbf{k}_{\perp},xy}(l).$
(5)

Now we can define the _faceted orbital magnetization_ as the orbital
magnetization contributed by the slab surface (or more generally a facet of a
three-dimensional crystallite). This quantity is obtained by summing the
layer-resolved magnetization in the vicinity of the surface, and subtracting
the bulk contribution:

$M_{+\hat{z}}=\sum_{l=1}^{bulk}\bigg{(}M_{z}(l)-M_{z,bulk}\bigg{)},$ (6)

with $M_{z,bulk}$ defined as the layer-resolved magnetization in the bulk
region; in a slab geometry, $M_{z,bulk}$ is the asymptotic value of $M_{z}(l)$
for large $l$. By $\sum_{l=1}^{bulk}$, we mean that the summation should be
carried out into the bulk region. The exponential decay of
$\tilde{P}_{R_{z}R_{z}^{\prime};ab}$ in $|R_{z}-R_{z}^{\prime}|$ guarantees
that the sum reaches its asymptotic value exponentially. The magnetic moment
contributed by a facet of area $A$ is then $M_{+\hat{z}}A$ – this form
manifests that the faceted magnetic moment is subextensive.

###  I.3 Faceted magnetic compressibility and chiral hinge modes

Differentiating Eq. (6) with respect to the chemical potential, we obtain the
magnetic compressibility contributed by a single facet (in short, the faceted
compressibility):

$\frac{dM_{+\hat{z}}}{d\mu}=-\frac{e}{hc}\sum_{l=1}^{bulk}C_{+\hat{z}}(l).$
(7)

The above equation can be derived straightforwardly from Eqs. (3)-(5). Note
that the derivative of $M_{z,bulk}$ vanishes by our initial assumption of the
triviality of the bulk Chern invariants. Hence, we can then identify
$\sum_{l=1}^{bulk}C_{+\hat{z}}(l)e^{2}/h$ as the geometric component of the
surface anomalous Hall conductivity (in short, the geometric SAHC, which we
denote by $\sigma_{+\hat{z}}^{g}$).[27, 28]. More generally, for any facet
with normal vector $\hat{\bm{n}}$,

$\frac{dM_{\hat{\bm{n}}}}{d\mu}=-\frac{1}{ec}\sigma_{\hat{\bm{n}}}^{g},$ (8)

where $\sigma_{\hat{\bm{n}}}^{g}$ is the geometric SAHC on the surface facet
with normal vector $\hat{\bm{n}}$. Because the spin contribution to the
magnetization does not depend on the chemical potential [cf. Sec. I.1], the
above equation represents the complete compressibility contributed by a single
facet, and not just its orbital component.

Figure 2: Geometry of facets (blue), hinges (green and orange), and Amperean
loops (red and black). The arrows on facets indicate the direction of the
faceted magnetic moment, assuming $M_{\hat{\bm{n}}}>0$ for each of
$\hat{\bm{n}}\in\\{+\hat{z},+\hat{y},-\hat{z}\\}$.

Let us emphasize a distinction between the faceted magnetization and the
magnetization of a strictly two-dimensional system $M_{2D}$. The latter
satisfies $dM_{2D}/d\mu=-(e/\hbar c)C$ where $C$ is the integer-valued Chern
number [6, 2], which remains quantized even in disordered systems[7]. In
contrast, $M_{\hat{\bm{n}}}$, as defined on a surface facet of a three-
dimensional crystallite, satisfies Eq. (7). Here the quantity
$\sum_{l=1}^{bulk}C_{\hat{\bm{n}}}(l)$ is not necessarily fixed to discrete
values. However, we will now show that the difference of
$\sum_{l=1}^{bulk}C_{\hat{\bm{n}}}(l)$ on adjacent facets must be quantized to
an integer, with this integer equal to (minus) the net chirality of the hinge
modes at the interface of the two facets.

Let us derive this equality from the classical electromagnetic relation
between magnetization and current: $\oint\mathbf{M}\cdot d\mathbf{l}=I/c$
[29], with $\bm{M}$ the total magnetization defined per unit volume, as in Eq.
(2). We construct an Amperean loop to encircle a hinge between two facets, as
illustrated by the red loop in Fig. 2. By assumption, the local Chern marker
in the bulk is trivial, hence $d\bm{M}/d\mu$ receives contributions only from
the faceted magnetization: $M_{+\hat{z}}\delta(z-z_{0})\hat{z}$ for the top
facet (colored blue in Fig. 2), and $M_{+\hat{y}}\delta(y-y_{0})\hat{y}$ for
the right facet (also colored blue). This means that

$\displaystyle\frac{dM_{+\hat{z}}}{d\mu}-\frac{dM_{+\hat{y}}}{d\mu}=\frac{1}{c}\frac{dI}{d\mu},$
(9)

for a current $I$ flowing through the Amperean loop. By assumption, there is
an energy gap to excitations for both bulk and surface states. Thus, a
variation of $\mu$ within this gap can only change the filling of hinge states
localized to the interface colored green in Fig. 2, and $dI/d\mu$ is
completely determined by the hinge state current. Each hinge state with wave
number $k$ and velocity $v(k)$ contributes $-ev(k)/L$ to the current $I$. Let
$s$ be a label for each hinge mode crossing the Fermi level with velocity
$v_{F}^{(s)}$ and wave number $k_{F}^{(s)}$. For an infinitesimal increase
$\delta\mu$, the resultant change in current is

$\displaystyle{\delta I=-\frac{e}{L}\sum_{s}\frac{L|\delta
k_{F}^{(s)}|}{2\pi}v_{F}^{(s)},}$ (10)

where $\delta k_{F}^{(s)}=\delta\mu/v^{(s)}\hbar$ is the change in Fermi wave
number for mode $s$. It follows that $\delta
I/\delta\mu=-(e/h)\sum_{s}\text{sgn}[v_{F}^{(s)}]$, hence

$\displaystyle{\frac{dM_{+\hat{z}}}{d\mu}-\frac{dM_{+\hat{y}}}{d\mu}=-\frac{e}{hc}\sum_{s}\text{sgn}[v_{F}^{(s)}],}$
(11)

which proves our claim. We remark that Eq. (11) is not unexpected given that
the faceted compressibility is related to the geometric SAHC through Eq.(8);
it is well-known that the differential geometric SAHC between two adjacent
facets determines the chirality of states localized at their interface [21,
30, 31]. Our present derivation, however, emphasizes the heretofore unexplored
connection between chiral hinge states and the faceted compressibility.

By a simple generalization of the above argument to the black Amperean loop in
Fig. 2, we derive

$\displaystyle{\frac{dM_{+\hat{z}}}{d\mu}-\frac{dM_{-\hat{z}}}{d\mu}=-\frac{e}{hc}\bigg{(}\sum_{s}\text{sgn}[v_{F}^{(s)}]+\sum_{s^{\prime}}\text{sgn}[v_{F}^{(s^{\prime})}]\bigg{)},}$
(12)

where $M_{-\hat{z}}$ is the magnetization of the bottom facet, and
$s^{\prime}$ is an index for hinge modes localized to the bottom hinge colored
orange in Fig. 2. The relative minus sign on the left-hand side of Eq. (12) is
because the faceted magnetization is defined to point from the interior to the
exterior.

Let us present an alternative and instructive derivation of Eq. (12). Because
the bulk magnetization has vanishing compressibility (since bulk Chern numbers
vanish), the right-hand side of Eq. (12) may be identified with the
compressibility of the total magnetization: $M_{z}=\sum_{l}M_{z}(l)$, where
the sum is over all layers in the slab; $M_{z}$ is defined with respect to a
fixed Euclidean coordinate system indicated on the right-bottom corner of Fig.
2. Applying the thermodynamic relation [in Eq. (1)] to the case of a
magnetization defined _per unit area_ ,

$\frac{\partial M_{z}}{\partial\mu}=\frac{1}{A_{xy}}\frac{\partial N}{\partial
B}=\frac{e}{hc}C_{\mathrm{slab},z},$ (13)

where $N$ is the total electron number and $A_{xy}$ is the area of the slab.
The second equality above relates $\partial N/\partial B$ (at fixed $\mu$ and
zero temperature) to the first Chern number ($C_{\mathrm{slab},z}$) of the
slab, in accordance with Streda’s formula [25]. By the bulk-boundary
correspondence for Chern insulators, there must be chiral modes localized to
the ‘edge’ of the slab. In this case, the ‘edge’ is a macroscopically large
side-surface facet. However, we have assumed that any surface states are
gapped, and thus the chiral ‘edge’ modes are identified with hinge modes
localized to either hinge.

The formulation of a faceted magnetization [cf. Eqs. (3), (4), (6)], and the
relation – between the faceted compressibility and the geometric SAHC [Eqs.
(8)] and chiral hinge states [Eqs. (11),(12)] – are the main results of this
section. With these formal preparations out of the way, we proceed to an
illustrative example.

##  II Case study of the Hopf insulator

The Hopf insulator is a three-dimensional bulk insulator with a two-band, bulk
Hamiltonian $H(\bm{k})$ that cannot be continuously deformed to a trivial
$\bm{k}$-independent Hamiltonian while the bulk energy gap and bulk
translation symmetry are preserved [15, 16, 32]. This topological obstruction
exists despite the Hopf insulator having trivial first Chern class – meaning
the Chern number of any two-dimensional slice in the three-dimensional
Brillouin zone vanishes. (We leave to future investigations the related notion
of Hopf-Chern insulators [33], which have a nontrivial first Chern class in
the bulk.) Instead of a bulk Chern invariant, the bulk topology of the Hopf
insulator is diagnosed by the Hopf invariant, which is equivalent to a BZ-
integral of the Abelian Chern-Simons three-form [15]:

$\chi=-\frac{1}{4\pi^{2}}\int_{BZ}d^{3}k\bm{A}(\bm{k})\cdot\bm{\mathcal{F}}(\bm{k}),$
(14)

where $\bm{A}(\bm{k})=i\langle u|\nabla_{\bm{k}}u\rangle$ is the Berry
connection of the valence band, and
$\bm{\mathcal{F}}(\bm{k})=\nabla\times\bm{A}(\bm{k})$ is the associated Berry
curvature. This Chern-Simons formula is gauge-invariant if the valence band
has unit rank, and is furthermore quantized to integer values if the
conduction band also has unit rank; both conditions are automatically
satisfied because the model we use for $H(\bm{k})$ is a sum of Pauli matrices.
(The rank of a band defined over the BZ or rBZ is defined as the number of
independent wave functions that span the band at each $\bm{k}\in BZ$, resp.
$\bm{k}_{\perp}\in rBZ$.)

Being an axion insulator with quantized magneto-electric polarizability [32],
the Hopf insulator is a prime candidate for the quantized magnetic
compressibility introduced in Sec. I.1. Additionally, the Hopf insulator
presents a unique opportunity to extend the Chern-Simons/axion theory of
magneto-electric polarizability to the case where the Chern-Simons action is
gauge-invariant, unlike other known axion insulators [27, 34, 35, 36, 4, 21,
37, 19, 20, 38, 28].

Let us summarize the main results of our case study:

(i) The faceted magnetic compressibility of the Hopf insulator is a half-
integer multiple of $e/hc$:

$\frac{dM_{\hat{\bm{n}}}}{d\mu}=\frac{e}{hc}\left(\frac{\chi}{2}-C_{v}(\hat{\bm{n}})\right).$
(15)

The half-integer includes a contribution from the bulk integer invariant
$\chi\in\mathbb{Z}$, and also a surface-dependent contribution from the
integer invariant $C_{v}(\hat{\bm{n}})\in\mathbb{Z}$. Roughly speaking,
$C_{v}$ is the first Chern number of the occupied bands localized to a surface
facet; the precise formulation of this quantity requires a Bloch-Wannier
representation that is elaborated in Sec. II.1.

(ii) Eq. (15) is derived by relating [cf. Eq. (8)] the facet compressibility
to the geometric component of the surface anomalous Hall conductivity (SAHC);
the latter quantity is related to the aforementioned topological invariants by

$\displaystyle{\sigma^{g}_{\hat{\bm{n}}}=\left(-\frac{\chi}{2}+C_{v}(\hat{\bm{n}})\right)\frac{e^{2}}{h},}$
(16)

as we prove in Sec. II.2.

(iii) From Eq. (16), we argue that if two neighboring facets differ in their
values for $C_{v}(\hat{\bm{n}})$, then there must be chiral hinge mode(s)
localized to the interface of the two facets – this follows from a generalized
Laughlin argument (for surface facets) that is described in Sec. II.2.2.
Furthermore, we show for the Hopf insulator that the existence of hinge modes
is guaranteed by a symmetry that combines particle-hole conjugation with
spatial inversion, as demonstrated in Sec. II.3.

(iv) Finally in Sec. II.4, we numerically simulate the quantization of Eq.
(15) for a model Hamiltonian.

###  II.1 Formulation of a Chern number for states localized to a facet

We have proposed, roughly speaking, that the integer invariant $C_{v}$ in Eq.
(15) is the Chern number for occupied states localized to a facet. To have a
properly defined Chern number (as an integral of the Berry curvature over the
reduced Brillouin zone $rBZ$), it is necessary to identify an entire band of
occupied states that is localized to the facet, with a projector that is
continuously defined over the _entire_ $rBZ$. This does not generally hold
when one evaluates energy eigenstates of topological band systems with a
surface termination. The more likely scenario is that of a ‘partial surface
band’, whose energy spectrum is attached to a bulk energy band in some region
of the $rBZ,$ such that the states localized to a facet (in position space)
are confined to a subregion $\mathfrak{D}$ of the $rBZ$ (in momentum space),
as exemplified in Fig. 3(a). In general, the integral of the Berry curvature
over $\mathfrak{D}$ is not quantized to integer values.

To obtain an integer-valued Chern invariant for surface-localized states in
the more general scenario of a partial surface band, one has to relax our
previous assumption that said states are energy eigenstates of the single-
particle Hamiltonian. Instead, we propose to define the Chern invariant
$C_{v}$ for surface-localized states that are eigenstates of the projected
position operator on a half-infinite geometry, as explained in Sec II.1.1.
(Note that eigenstates of the projected position operator are simply linear
combinations of the occupied energy eigenstates.) Then, in Sec II.1.2, we
describe the conditions that make $C_{v}$ a well-defined topological
invariant. One implication is that if a three-dimensional, bounded crystallite
is insulating in the bulk and on all surface facets, and has a bulk Pauli-
matrix Hamiltonian in the trivial first Chern class, then bulk translational
symmetry alone gives a topological classification specified the bulk invariant
$\chi$ [cf. Eq. (14)] and the set of $C_{v}$ for all facets.

Figure 3: (a) Partial surface bands in the open-boundary spectrum of the model
Hamiltonian in Eq. (19), with $u=2$. The states colored red are localized to
surface facets. $\Delta E_{b,m}$ is minimal energy gap in the bulk. $\epsilon$
is set to be the natural energy scale in this paper. $m_{D}$ is the surface
Dirac mass. (b) Schematic illustration of Bloch-Wannier bands for a
superlattice model described in Sec. II.1.1.

####  II.1.1 Chern invariant $C_{v}$ from the projected position operator

Given a bulk insulator with an occupied subspace of states (represented by the
projector $P$) and an unoccupied subspace (with a projector $Q$), we consider
the projected position operators $P\hat{Z}P$ and $Q\hat{Z}Q$, with $\hat{Z}$
the position operator in the $z$ direction. Throughout this subsection we
impose periodic boundary conditions in the $x$ and $y$ directions.

While the eventual goal of this section is to define the Chern invariant
$C_{v}$ through the projected position operator defined on a half-infinite
slab geometry (with $\hat{Z}$ taking only positive values), it is useful to
review some properties of $P\hat{Z}P$ on an infinite geometry (with $\hat{Z}$
taking also negative values). These properties hold equally well for
$Q\hat{Z}Q$, with the semantic replacement $P\rightarrow Q$. Since
translations perpendicular to $z$ are symmetries of $P\hat{Z}P$, each
eigenvalue of $P\hat{Z}P$ can be labelled by a reduced wavevector
$\bm{k}_{\perp}{=}(k_{x},k_{y})$ in a reduced Brillouin zone rBZ, and
eigenvalues that are continuously connected form a band. The corresponding
eigenstates, known as Bloch-Wannier bands [39, 40], are extended (in the
$xy$-plane) as Bloch functions with crystal wavevector $\bm{k}_{\perp}$, but
exponentially localized as a Wannier function in the $z$ direction. Such
exponential localization makes the expectation value
$\langle{\hat{Z}}\rangle,$ for a Bloch-Wannier eigenstate, well-defined, with
$\langle{\hat{Z}}\rangle$ equalling an eigenvalue of $P\hat{Z}P$. Under
translation in the $z$-direction by a lattice period (set to unity),
$P\hat{Z}P\rightarrow P(\hat{Z}+1)P$, hence each eigenvalue in a given
$\bm{k}_{\perp}$ sector belongs to an infinitely-extended Wannier-Stark ladder
where adjacent rungs having unit spacing [41]. Generally, if the bulk valence
energy band has rank $N$, the spectrum of $P\hat{Z}P$ (at each
$\bm{k}_{\perp}$) comprises $N$ Wannier-Stark ladders – each of them infinite
in extent. For the Hopf insulator, $N=1$.

Let us now consider a half-infinite slab geometry with $\hat{Z}$ taking only
positive values. We assume that a gap exists in the energy spectrum over the
reduced Brillouin zone (rBZ) of the surface facet, and separates two
orthogonal subspaces projected by $P$ (for occupied states) and $Q$ (for
unoccupied states). (Such a gap is consistent with the possible existence of
chiral hinge modes, because the surface geometry considered here is periodic
without hinges.) Let $j=1,2,3\ldots$ label Bloch-Wannier bands with $j=1$
being closest to the surface, $j=2$ the next closest, and so on. Owing to the
exponential decay of $P$ in the coordinate representation [22, 23], Bloch-
Wannier bands sufficiently far from the surface are indistinguishable (up to
exponentially small corrections) from Bloch-Wannier bands obtained in the
above-described infinite geometry. We describe the ‘sufficiently far’ Bloch-
Wannier bands as bulk-like, while Bloch-Wannier bands in the infinite geometry
are bulk Bloch-Wannier bands; every other Bloch-Wannier band is a surface
band. In the usual convention, we attach the additional qualifier ‘valence’
and ‘conduction’ to objects associated to the occupied and unoccupied
subspaces respectively, e.g., valence bulk-like Bloch-Wannier bands.

Let us denote the Chern number of the $j$-th Bloch-Wannier band in the
occupied subspace as $C_{v}(+\hat{z},j)$. The faceted valence Chern number
$C_{v}(+\hat{z})$ is defined as the net Chern number of all surface valence
bands: $\sum_{j}C_{v}(+\hat{z},j)$, with the sum extending to the bulk-like
region. This sum is uniquely defined because bulk-like Bloch-Wannier bands
have zero Chern number, by assumption. Analogously, we can define
$C_{c}(+\hat{z},j)$ and $C_{c}(+\hat{z})$ for the conduction band, $i.e.$,

$\begin{array}[]{c}C_{v}(+\hat{z})=\sum_{j=1}^{bulk}C_{v}(+\hat{z},j),\\\
C_{c}(+\hat{z})=\sum_{j=1}^{bulk}C_{c}(+\hat{z},j).\end{array}$ (17)

In numerics, one can calculate the faceted valence/conduction Chern numbers by
summing $j$ over half a slab, assuming the slab height is much greater than
the exponential decay length of $P$ and $Q$. The definitions in Eq. (17) are
straightforwardly generalized to $C_{v,c}(\hat{\bm{n}})$ any insulating
surface facet with normal vector $\hat{\bm{n}}$. Then, we define the faceted
Chern number of a surface as the sum of the faceted valence Chern number and
the faceted conduction Chern number:

$\displaystyle C(\hat{\bm{n}})=C_{v}(\hat{\bm{n}})+C_{c}(\hat{\bm{n}})=\chi,$
(18)

where the second equality relates the faceted Chern number to the bulk Hopf
invariant $\chi$ [cf. Eq.(14)] – a bulk-boundary correspondence proven in Ref.
32.

This correspondence suggests that the faceted valence Chern number
$C_{v}(+\hat{z})$, just like the Hopf invariant, is only well-defined if the
bulk valence energy band is unit-rank with trivial first Chern class. To
appreciate this point, we offer an example with non-unit rank where $C_{v}$ is
not well-defined. Our example is constructed by stacking decoupled, two-
dimensional Chern insulators (in the $z$-direction) to form a three-
dimensional superlattice with two layers per superlattice period, and with the
Chern number alternating in sign between adjacent layers. The result is a
topologically trivial insulator with two Bloch-Wannier bands per lattice
period in the z-direction, with the Chern number alternating in sign between
adjacent Bloch-Wannier bands, as schematically shown in Fig. 3(b). It follows
that $C_{v}$ is not uniquely defined and depends on how the sum in Eq. (17) is
truncated: if the number ($N$) of Bloch-Wannier bands included in the
summation is even (see dotted square in Fig. 3 (b)), then $C_{v}=0$, while if
$N$ is odd, $C_{v}=1.$

To close this subsection, let us illustrate the notions of Bloch-Wannier
bands, faceted Chern numbers, and the bulk-boundary correspondence for a
tight-binding model of the Hopf insulator constructed by Moore, Ran, and Wen
(MRW) [15]. The Hamiltonian of the MRW model is constructed from the well-
known Hopf map:

$\displaystyle z=(z_{1}+iz_{2},z_{3}+iz_{4})^{T},$ (19)
$\displaystyle\mathbf{d}=z^{{\dagger}}\bm{\tau}z,\
\bm{\tau}=(\tau_{x},\tau_{y},\tau_{z}),$ $\displaystyle
H_{MRW}(\mathbf{k})=\mathbf{d}\cdot\bm{\tau},$

where $\tau_{x},\tau_{y},\tau_{z}$ are Pauli matrices, and

$\displaystyle z_{1}=\sin k_{x},z_{2}=\sin k_{y},z_{3}=\sin k_{z},$ (20)
$\displaystyle z_{4}=u-\cos k_{x}-\cos k_{y}-\cos k_{z}.$

If we take $1<u<3$, then the above model is a Hopf insulator with Hopf
invariant $\chi=1$. In numerical calculations, we always use the MRW model
with $u=2$ and $\chi=1$ unless otherwise specified.

To calculate the Bloch-Wannier bands in the presence of a surface, we Fourier
transform the Bloch Hamiltonian of the MRW model to obtain the position-space
Hamiltonian. Then we terminate the hopping matrix elements across the surfaces
of a finite slab that is open in the $z$-direction and periodic in $x$ and
$y$. After modifying the surface Hamiltonian to gap out both the top and
bottom surfaces, we diagonalize $P\hat{Z}P$ and $Q\hat{Z}Q$ to derive Bloch-
Wannier bands as shown in Fig. 4(a) and (b). If we sum the Chern number of the
top five bands of $P\hat{z}P$ in Fig. 4(a), we get $C_{v}(+\hat{z})=1$. If we
sum the Chern number of the top five bands of $Q\hat{z}Q$ in Fig. 4(b), we get
$C_{c}(+\hat{z})=0$. Then, we have
$C(+\hat{z})=C_{v}(+\hat{z})+C_{c}(+\hat{z})=1$, which is consistent with
$\chi=1$, according to the bulk-boundary correspondence in Eq. (18).

Figure 4: Spectrum of (a) $P\hat{z}P$ and (b) $Q\hat{Z}Q$ for the MRW model
with $u=2$ and $\chi=1$, for a finite slab with a width of ten unit cells in
the $z$ direction. Except for Bloch-Wannier bands in the vicinity of the top
and bottom facets, all other bands are bulk-like with near-unit separation.

####  II.1.2 $C_{v}$ as a topological invariant

Let us prove that $C_{v}$ is invariant under continuous deformations of the
Hamiltonian that preserve the bulk translation symmetry as well as the energy
gap separating valence and conduction subspaces (on a half-infinite geometry).

The only way for $C_{v}$ to change is by the transference of Berry-curvature
quanta to (or away from) the facet. If the energy gap is preserved, then there
can be no transference between occupied and unoccupied states. What remains is
to rule out transference between occupied states. Because Bloch-Wannier bands
have differing spatial centers $\langle\hat{Z}\rangle$, the transference of a
Berry-curvature quantum between neighboring Bloch-Wannier bands effectively
moves the quantum in real space. Such transference can only occur if the
spectral gap separating the two Bloch-Wannier bands closes then reopens;
generically, the gap closes at a Dirac point. The transference of Berry quanta
by Dirac-point touchings of Bloch-Wannier bands is sometimes referred to as an
axion pump, and has been explored in Refs. 42, 43, 31 in the context of
topological insulators whose bulk valence bands have rank larger than one.
Yet, in our context, if bulk translation symmetry is preserved, it is not
possible for bulk/bulk-like Bloch-Wannier bands (with unit rank) to touch,
because they are always separated by the lattice period in the $z$-direction;
by implication, the axion pump is forbidden. An analogous argument proves that
$C_{c}$ (for the unoccupied subspace) is also a topological invariant.

We can provide an alternative derivation of this result by assuming
translation invariance in all three spatial directions. An adiabatic cycle for
a bulk pump of Berry-curvature quanta corresponds to a nontrivial second Chern
number for the Hamiltonian $H(\bm{k},\lambda)$[42, 32], where $\lambda$ is a
cyclic parameter for a family of Bloch Hamiltonians $H(\bm{k})$ defined over
the 3D Brillouin zone. The second Chern number is defined as a four-
dimensional integral of the second Chern character. Generally, the $n$’th
Chern character is given by

$\displaystyle{ch_{n}(F)=\text{Tr}\left[\frac{1}{n!}\bigg{(}\frac{iF}{2\pi}\bigg{)}^{n}\right],}$
(21)

with $F^{2}=F\wedge F$ the exterior product of the Berry field strength.
However, the second Chern number is always trivial for a non-degenerate band
spanned by a single Bloch function [of $(\bm{k},\lambda)$] that is continuous
with respect to $\bm{k}$. This follows from a relation between the first two
Chern characters and the second Chern class $c_{2}$:[44]

$\displaystyle{c_{2}(F)=\frac{1}{2}ch_{1}(F)\wedge ch_{1}(F)-ch_{2}(F),}$ (22)

and the triviality of the second Chern class (i.e., $\int
d^{3}kd\lambda\,c_{2}{=}0$) for line bundles [45]. Note $\int
d^{3}kd\lambda\,ch_{1}{\wedge}ch_{1}{=}0$ because, by assumption, the first
Chern numbers vanish in the $k_{x}k_{y}$, $k_{y}k_{z}$ and $k_{z}k_{x}$
planes. This illustrates that adiabatic deformations of the Bloch Hamiltonian
will not transfer any Berry curvature quanta away from the facet. A similar
argument was used in Ref. 32 to demonstrate the topological invariance of the
faceted Chern number $C$ [cf. Eq. (17)].

###  II.2 Surface anomalous Hall conductivity

The half-integer quantization of the faceted compressibility of the Hopf
insulator [cf. Eq. (15)] can be proven once we relate the geometric SAHC to
the bulk topological invariant $\chi$ [cf. Eq. (14)] and the faceted valence
Chern number $C_{v}$ [cf. Eq. (17)] through Eq. (16). As preliminary steps, we
will first review some salient aspects of the general theory of surface Hall
conductivity in Sec. II.2.1, then generalize Laughlin’s argument to surface
facets in Sec. II.2.2. This general theory is subsequently applied to the Hopf
insulator in Sec. II.2.3.

####  II.2.1 Review of the general theory of magnetoelectric polarizability
and SAHC

For an insulator with a gapped surface facet (with normal vector
$\hat{\bm{n}}$), the application of an electric field results in a surface
Hall current:
$\bm{j}_{\hat{\bm{n}}}=\sigma_{\hat{\bm{n}}}\hat{\bm{n}}\times\bm{E}$. The
total SAHC $\sigma_{\hat{\bm{n}}}$ can be separated into geometric
($\sigma_{\hat{\bm{n}}}^{g}$) and non-geometric contributions[20, 27, 34, 35,
4, 21, 28, 38, 30, 46]:

The geometric contribution comprises two terms, the first of which is
proportional to the axion angle $\theta.$ The value of $\theta$ is expressible
as a Brillouin-zone integral of the Chern-Simons three-form of the Berry gauge
field [20]. For a unit-rank valence energy band, the Chern-Simons three-form
is Abelian, and the axion angle reduces to $\pi$ multiplied by the right-hand
side of Eq. (14), i.e., it is $\pi\chi.$ The second term in
$\sigma_{\hat{\bm{n}}}^{g}$ is an integer multiple of $e^{2}/h$ that depends
on the surface preparation. The non-geometric contribution to
$\sigma_{\hat{\bm{n}}}$ also comprises two terms, one being facet-independent
(the “cross-gap” term $\alpha_{iso}^{cg}$), and the other facet-dependent.
Analytic expressions for both these terms can be found in Refs. 34, 35, but
will not be the focus of this work.

####  II.2.2 SAHC of insulators with trivial bulk Chern invariants: a
generalized Laughlin’s argument

The differential SAHC between any two adjacent facets is an integer multiple
of $e^{2}/h$, with this integer equal to the net chirality of hinge modes
localized to the interface of the two facets. This follows from a
generalization of Laughlin’s argument to surface facets that we now develop.

Laughlin’s original argument was applied to prove the integer quantization of
the Hall conductivity for two-dimensional electron gases [47], using flux
threaded in a cylindrical geometry. Our generalized argument is based on the
geometry illustrated in Fig. 5(a): two adjacent surface facets with normal
vectors $\hat{\bm{n}}_{1}$ and $\hat{\bm{n}}_{2}$, having an interfacial
“hinge” of length $L_{x}$ parallel to $+\hat{x}$. Periodic boundary conditions
are imposed in the $x$ direction, while each facet is assumed to be semi-
infinite in the direction orthogonal to $x$.

Let us introduce a time-dependent flux $\Phi(t)$ into the Hamiltonian by the
minimal coupling $p_{x}\rightarrow p_{x}-(e/c)\Phi/L$, with $p_{x}$ the
canonical momentum operator. One can remove the flux from the Hamiltonian via
a unitary transformation which twists the boundary conditions of single-body
wave functions as:
$\psi(x,y,z)=\exp\left(-i\Phi/\Phi_{0}\right)\psi(x+L_{x},y,z)$, with
$\Phi_{0}=hc/e$ the flux quantum. The spectrum of the Hamiltonian is identical
for any $\Phi$ that is an integer multiple of $\Phi_{0}$. However, it is
possible that a filled Fermi sea ground state (at $\Phi=0$) is adiabatically
connected to an excited many-body state after a flux quantum is threaded (at
$\Phi=\Phi_{0}$). This happens only if there are states (near the Fermi level)
that are extended in the $x$-direction; such extension makes them sensitive to
twisting of the boundary condition along $x$. Since we have assumed the bulk
and surfaces are gapped, any extended state (near the Fermi level) will be a
quasi-1D Bloch wave localized to the hinge. If hinge modes crossing the Fermi
level have a net chirality for propagation in the $x$-direction, then the
many-body state at $\Phi=\Phi_{0}$ is excited when compared to $\Phi=0$. For
illustration, if the spectrum contains a single chiral hinge mode (illustrated
in Fig. 5 (b)), then all single-body states are advanced in quasimomentum by
$k_{x}\rightarrow k_{x}+2\pi/L_{x}$, resulting in one additional filled state
just above the Fermi level.

Figure 5: (a) Two adjacent facets with normal vectors $\hat{\bm{n}}_{1}$ and
$\hat{\bm{n}}_{2}$. Their interface is emphasized by a purple colored line.
(b) Illustration of spectral flow for a single chiral hinge mode. The purple
dashes indicate hinge states with wavenumber $k_{x}$ that takes discrete
values owing to finite-size quantization.

More generally, let the net chirality of the hinge modes be $\Delta
N=N_{R}-N_{L}$, implying that the number of occupied hinge states increases by
$\Delta N$ upon inserting of a flux quantum. By particle conservation, these
states must originate from a net flow of current transverse to the hinge:
$(-e)\Delta N=L\int_{0}^{T_{0}}(j^{(1)}-j^{(2)})dt$, with $j^{(i)}$ the areal
current density on the $i$-th facet. We note that there is no current
contributed from the bulk, because the bulk is insulating, and the bulk Chern
invariants are assumed to vanish. A time-dependent flux $\Phi(t)=LA_{x}(t)$
can be viewed as proportional to a time-dependent electromagnetic vector
potential, which implies a non-vanishing electric field: $E_{x}=(-1/c)\partial
A_{x}/\partial t$. Assuming $E_{x}$ is constant in time, then
$\Phi(T_{0})=\Phi_{0}$ at time $T_{0}=-h/eE_{x}L$. Substituting this time
scale into the above expression for $\Delta N$, and further identifying
$j^{(i)}=\sigma_{\hat{\bm{n}}_{i}}E_{x}$ as a Hall current (since both
surfaces are gapped), we derive that the difference in Hall conductivity is
determined by the net chirality:
$\sigma_{\hat{\bm{n}}_{1}}-\sigma_{\hat{\bm{n}}_{2}}=(e^{2}/h)\Delta N$. This
completes our argument.

We can now use this result to interpret the general expression for the SAHC in
Eq. LABEL:eq:_surfhall_general. First, we argue that in three-dimensional
insulators whose bulk valence band has trivial Chern class, the anisotropic
term $\frac{1}{2}\widetilde{\alpha}_{ab}\hat{n}_{a}\hat{n}_{b}$ in Eq.
(LABEL:eq:_surfhall_general) must vanish. This follows because on one hand,
(i) if $\widetilde{\alpha}$ were nonzero,
$\frac{1}{2}\widetilde{\alpha}_{ab}\hat{n}_{a}\hat{n}_{b}$ can in principle
take arbitrary real values depending on the facet orientation $\hat{\bm{n}}$,
implying that the differential
$\frac{1}{2}\widetilde{\alpha}_{ab}\hat{n}_{a}\hat{n}_{b}$ between two facets
can take arbitrary real values. On the other hand, (ii) we just proved that
the differential SAHC between two adjacent facets must be an integer multiple
of $e^{2}/h$. Since the real-valued arbitrariness of (i) is incompatible with
the integer-valued quantization of (ii), our previous assumption of
$\widetilde{\alpha}\neq 0$ cannot be true.

Secondly, we note that $\Delta N$ can be identified with the differential
$m_{s}$ between two adjacent facets. After all, having just argued that
$\widetilde{\alpha}=0$, $m_{s}$ is the only remaining nonzero quantity in Eq.
(LABEL:eq:_surfhall_general) that depends on the facet orientation.

####  II.2.3 SAHC of the Hopf insulator

Specializing the general expression of the SAHC in Eq.
(LABEL:eq:_surfhall_general) to the case of the Hopf insulator, we obtain

$\begin{array}[]{c}\sigma_{\hat{\bm{n}}}=\sigma_{\hat{\bm{n}}}^{g}-\alpha_{iso}^{cg},\;\;\;\;\sigma_{\hat{\bm{n}}}^{g}=\left(-\frac{\chi}{2}+m_{s}(\hat{\bm{n}})\right)\frac{e^{2}}{h}.\end{array}$
(23)

We have replaced the (Abelian) axion theta angle by the Hopf invariant:
$\theta_{Abelian}=\pi\chi$, in accordance with Eq. (14). (Note that $\chi$ is
gauge-invariant, unlike the non-Abelian version where
$\theta_{NA}\equiv\theta_{NA}+2\pi$.) Moreover, the term involving
$\tilde{\alpha}$ vanishes as an implication of the generalized Laughlin
argument given above.

Next we will argue that we can identify the surface-dependent integer
invariant $m_{s}(\hat{\bm{n}})$ with the faceted valence Chern number (c.f.
Eq. (17)):

$\displaystyle{m_{s}(\hat{\bm{n}})=C_{v}(\hat{\bm{n}}).}$ (25)

To prove this identity, consider a half-infinite slab of the Hopf insulator
with surface-normal vector $\hat{\bm{n}},$ and total SAHC
$\sigma_{\hat{\bm{n}}}$. The surface valence bands have the net Chern number
$C_{v}(\hat{\bm{n}}),$ which depends on details of the surface termination. As
a thought experiment, let us absorb a two-dimensional Chern-insulating layer
onto the surface which has the opposite Chern number $-C_{v}(\hat{\bm{n}})$ in
its valence subspace. We do not, however, introduce any Hamiltonian matrix
elements between the absorbed layer and the original system. The combined
system is then a half-infinite, Hopf-insulating slab whose occupied subspace
can be spanned by a set of exponentially-localized Wannier functions. Then by
following the derivations in Refs. 35, 4 (which apply only to systems with a
Wannier representation), we would deduce that the SAHC of the combined system
equals $-\chi/2(e^{2}/h)-\alpha_{iso}^{cg}$. In the assumed absence of
coupling between the Chern-insulating layer and the original slab, the SAHC of
the combined system must equal the sum of the SAHC of the constituent systems,
thus

$\sigma_{\hat{\bm{n}}}-C_{v}(\hat{\mathbf{n}})\frac{e^{2}}{h}=-\frac{\chi}{2}\frac{e^{2}}{h}-\alpha_{iso}^{cg}.$
(26)

Comparing this to the first equation in Eq.(23), we deduce
$m_{s}(\hat{\bm{n}})=C_{v}(\hat{\bm{n}})$, and hence, for the Hopf insulator

$\sigma_{\hat{\bm{n}}}^{g}=\left(-\frac{\chi}{2}+C_{v}(\hat{\bm{n}})\right)\frac{e^{2}}{h}.$
(27)

We remark that for Hopf insulators whose bulk Hamiltonian $H(\bm{k})$ is
traceless, as exemplified by the MRW model in Eq. (19), the SAHC is purely
geometric. This is because the traceless condition ensures that the
‘degeneracy’ and ‘reflection’ conditions [discussed in Ref 34] are satisfied,
which imply the vanishing of the non-geometric, ‘cross-gap’ SAHC.

###  II.3 Chiral hinge modes and their symmetry protection

We have just determined that the net chirality of hinge modes separating two
surface facets equals the differential $C_{v}$ between the two adjacent
facets. Given a Hopf-insulating crystallite with a fixed bulk invariant
$\chi,$ it is possible to manufacture a surface potential that leads to many
choices of $\\{C_{v}(\hat{\bm{n}})\\}_{\hat{\bm{n}}}$, and by implication many
possible networks of hinge modes. Three representative examples consistent
with half filling of the MRW model are illustrated in Fig. 6(a-c), with blue
indicating a facet with $C_{v}=+1$, light orange indicating a facet with
$C_{v}=0$, and red indicating a unidirectional hinge mode. These three
simulations differ only in the surface preparations, while their bulk is
identically given by the MRW model [cf. Eq. (19)] with parameter $u=2,$ and
bulk invariant $\chi=1$. For the two examples in Fig. 6(a-b), the dispersion
of the hinge modes propagating in the $\pm x$ direction are plotted in Fig.
6(d-e) respectively.

Figure 6:  (a-c) shows the possible networks of hinge modes for a cubic Hopf-
insulating crystalline with bulk Hopf invariant $\chi=1$. Facets with
$C_{v}=1$ ($C_{v}=0$) are colored light orange (blue, respectively). Black
arrows represent the the group velocity of chiral hinge modes at the Fermi
level. The particle-hole inversion symmetry of Eq. (28) is imposed in panels
(b) and (c), but not in (a). The dispersion of hinge modes localized to the
red-highlighted hinges in panel (a) [resp. (b)] are plotted in panel (d)
[resp. (e)]. Green-colored hinges also hosts hinge modes which are, however,
not shown in an energy-momentum plot. Simulations for (d-e) were performed for
the MRW model ($u=2$) on a finite slab that is periodic in the $x$ direction,
with a width of 10 unit cells in the $z$-direction, and a width of 40 unit
cells in the $y$-direction. Figure 7:  Simulations on a finite slab that is
open in the $y$-direction, with 20 unit cells. Figure (a) plots
$\ln\left|M_{y}(l)/(e\epsilon/hc)\right|$ _v.s._ the layer index $l$, with
$l=1$ closest to the surface termination; we have chosen the surface Dirac
mass to be $m_{D}=0.1\epsilon$. Exponential decay of the magnetization away
from the surface is observed. (b) plots $M_{+\hat{y}}$ (at $\mu=0$) as a
function of $m_{D}$. (c) plots $M_{+\hat{y}}$ _vs_ $\mu$ with
$m_{D}=0.5\epsilon$. The constant slope as a function of $\mu$ indicates a
quantized compressibility.

Interestingly, it is possible that all facets have identical $C_{v}$ such that
there are no chiral hinge modes at all. Hence, one might consider if a
symmetry exists that might prevents this possibility and thus guarantee
higher-order topology with robust hinge modes. As a candidate for such a
symmetry, let us consider the following transformation of annihilation to
creation operators

$\displaystyle{\mathcal{C}^{\prime}:c_{\bm{R},a}\rightarrow\epsilon_{ba}c_{-\bm{R},b}^{{\dagger}}.}$
(28)

Here, $c_{\bm{R},a}^{{\dagger}}$ creates an electron of orbital character
$a\in\\{1,2\\}$ at lattice site $\bm{R}$, and $\epsilon_{ab}$ is the two-
dimensional Levi-Civita tensor. Since $\mathcal{C}^{\prime}$ is the
composition of a spatial inversion ${\cal I}:c_{\bm{R},a}\rightarrow
c_{-\bm{R},a}$ and a particle-hole conjugation ${\cal
P}:c_{\bm{R},a}\rightarrow\epsilon_{ba}c_{\bm{R},b}^{{\dagger}}$ (with pseudo-
spin flip), we refer to $\mathcal{C}^{\prime}={\cal P}{\cal I}$ as a particle-
hole inversion symmetry. Note that neither ${\cal P}$ nor ${\cal I}$ symmetry
is individually imposed in the following discussion. The
$\mathcal{C}^{\prime}$ symmetry of the second-quantized Hamiltonian translates
to the following antisymmetry of the first-quantized Hamiltonian in the
momentum representation:

$\displaystyle{\tau_{y}H^{\star}(\bm{k})\tau_{y}=-H(\bm{k}),\;\;\;\;\bm{k}\in
BZ,}$ (29)

where $H^{\star}(\bm{k})$ means the complex conjugate of the matrix
$H(\bm{k})$. This constraint is satisfied by any traceless, two-by-two matrix
Hamiltonian.

We propose that $\mathcal{C}^{\prime}$ symmetry guarantees hinge modes for
Hopf insulators with odd bulk invariant $\chi$, by ensuring that $C_{v}$ on
inversion-related facets are distinct, as illustrated in Fig. 6(b) and (e). To
derive the necessity of hinge modes, let us express $C=C_{v}+C_{c}$ as an
integral (over the rBZ) of the Abelian Berry curvature of the surface valence
(and conduction) bands:

$\displaystyle C(\hat{\bm{n}})$
$\displaystyle=\frac{1}{2\pi}\int_{rBZ}d^{2}k\big{\\{}Tr\left[\bm{\mathcal{F}}(\bm{k})_{\hat{\bm{n}},v}\cdot\hat{\bm{n}}\right]$
(30)
$\displaystyle+Tr\left[\bm{\mathcal{F}}(\bm{k})_{\hat{\bm{n}},c}\cdot\hat{\bm{n}}\right]\big{\\}}=C_{v}(\hat{\bm{n}})+C_{c}(\hat{\bm{n}}).$

We can understand how $\mathcal{C}^{\prime}$ constrains the Berry curvature by
decomposing its action as described above. Since the Berry curvature is a
pseudo-vector, the spatial inversion transformation sends
$\bm{\mathcal{F}}(\bm{k})_{\hat{\bm{n}},v}$ to
$\bm{\mathcal{F}}(-\bm{k})_{-\hat{\bm{n}},v}$. Next, the particle-hole
conjugation transformation subsequently sends
$\bm{\mathcal{F}}(-\bm{k})_{-\hat{\bm{n}},v}$ to
$-\bm{\mathcal{F}}(\bm{k})_{-\hat{\bm{n}},c}$. The minus sign is due to the
anti-unitary action of particle-hole conjugation in $\bm{k}$-space. Thus,
$\mathcal{C}^{\prime}$ ultimately constrains
$\bm{\mathcal{F}}(\bm{k})_{\hat{\bm{n}},v}=-\bm{\mathcal{F}}(\bm{k})_{-\hat{\bm{n}},c}$,
and therefore

$\displaystyle
C_{v}(\hat{\bm{n}})=\frac{1}{2\pi}\int_{rBZ}d^{2}kTr\left[\bm{\mathcal{F}}(\bm{k})_{\hat{\bm{n}},v}\cdot\hat{\bm{n}}\right]$
(31)
$\displaystyle=\frac{1}{2\pi}\int_{rBZ}d^{2}kTr\left[\bm{\mathcal{F}}(\bm{k})_{-\hat{\bm{n}},c}\cdot\left(-\hat{\bm{n}}\right)\right]=C_{c}(-\hat{\bm{n}}).$

Combining this with the bulk-boundary correspondence [cf. Eq. (18)], we find
that if $\chi$ is odd, the parities of $C_{v}(\hat{\bm{n}})$ and
$C_{v}(-\hat{\bm{n}})$ are non-identical. This would indicate that some change
of the SAHC is required over neighboring facets, and hence there will be
localized hinge modes. However, for $\chi$ even, it is possible that
$C_{v}(\hat{\bm{n}})=C_{v}(-\hat{\bm{n}})=\chi/2$ for all $\hat{\bm{n}}$,
implying no generic, robust hinge modes.

Upon enlargement of the Hilbert space by addition of trivial bands, it was
discovered that Hopf insulators (with odd Hopf invariants) remain
homotopically inequivalent to a trivial insulator – assuming an antisymmetry
constraint that is a higher-rank generalization of Eq. (29) [48]. Whether or
not these higher-rank insulators are also higher-order topological insulators
is a question we leave for future investigation.

###  II.4 Numerical simulation of the magnetic compressibility

For the MRW model, we will numerically demonstrate that the layer-resolved
magnetization is indeed exponentially-localized to a surface facet, and
sensitive to details of the surface potential. On the other hand, the faceted
magnetic compressibility $dM_{\hat{\bm{n}}}/d\mu$ is a half-integer-valued
topological invariant if the Hopf invariant $\chi$ is odd, no matter the
details of the surface termination.

Beginning with the position-space representation of the MRW model with
parameters ($u=2,\chi=1$), we diagonalize the Hamiltonian with a surface
termination in the $+\hat{y}$ direction, by abruptly truncating all
Hamiltonian matrix elements that extend across the surface. With this
prescription, the surface energy bands touch at one Dirac point owing to an
artifact of the simplified model. To gap the surface, we split this Dirac-
point degeneracy by adding the Hamiltonian term:
$m_{D}(\mathcal{P}_{l=2}-\mathcal{P}_{l=1})$, where $m_{D}$ can be interpreted
as a Dirac mass, and $\mathcal{P}_{l=1,2}$ projects to the first and second
layers (counting from the surface termination). We then fix the Fermi level in
the resultant energy gap, and compute $M_{y}(l)$ on the surface with normal
vector $+\hat{y}$ through Eq.(3) (with $z$ replaced by $y$). Fig. 7(a) plots
the absolute value $|M_{y}(l)|$ vs $l$, for the Dirac mass
$m_{D}=0.1\epsilon$, where $\epsilon$ is half the bulk energy gap [see Fig.
3(a)]. This simulation demonstrates that the faceted orbital magnetization is
indeed exponentially localized at the surface. What is more, as shown in Fig.
7(b), the faceted orbital magnetization on the surface with normal vector
$+\hat{y}$ can be tuned continuously by varying the Dirac mass $m_{D}$. Since
the Dirac mass is experimentally tunable (at least on the surface of a
$\mathbb{Z}_{2}$ topological insulator [49]), we see that, in principle, the
faceted orbital magnetization is tunable by surface manipulations.

In contrast, Fig. 7(c) shows that $M_{+\hat{y}}$ depends linearly on $\mu$
with slope equal to $\sigma_{+\hat{y}}^{g}=(1/2)(e^{2}/h)$. This is consistent
with Eq. (23) which states that the slope equals $(C_{v}-\chi/2)e^{2}/h$;
here, $\chi=1$ follows from our choice of the bulk Hamiltonian parameter
$u=2$, and $C_{v}=1$ on the surface facet with normal vector $+\hat{y}$, due
to the specific nature of our surface termination.

##  III Topomagnetic insulators (and beyond)

The half-integer-quantized faceted compressibility is not a unique property of
the Hopf insulator; rather, it applies to any insulator whose geometric SAHC
is a half-integer multiple of $e^{2}/h$, according to the general relation
between compressibility and surface anomalous Hall conductivity in Eq. (8).
This class of insulators, sometimes referred to as axion insulators, are
exemplified by (i) the three-dimensional $\mathbb{Z}_{2}$ topological
insulators [19, 20] with gapped surface states [20, 21], (ii) the 3D, second-
order topological insulators protected by a composition of time reversal with
four-fold rotation ($TC_{4}$) [10], and (iii) certain anti-ferromagnetic
topological insulators with bulk axion angle $\theta=\pi$ and with gapped
surfaces [50, 51].

In examples (i-ii), the bulk magnetization vanishes owing to symmetry – time
reversal in the first case, and $TC_{4}$ in the second. This is because both
spin and orbital magnetization [cf. Sec. I.1] transform as a pseudo-vector
under crystallographic spatial transformations, and are moreover odd under
time reversal, as shown in Appendix B. In example (iii), the bulk
magnetization vanishes because the magnetic order is anti-ferromagnetic.
(i-iii) represent a counter-intuitive phase of matter which is insulating (in
the bulk and on surface facets) with vanishing bulk magnetization under zero
field, but has a half-integer-quantized magnetic compressibility on its
surface – we call this a topomagnetic insulator. Having zero bulk
magnetization renders the faceted magnetization more readily measurable, with
techniques that are tentatively proposed in Sec. IV.

###  III.1 The Hopf insulator as an anisotropic analog of the topomagnetic
insulator

Not all axion insulators are topomagnetic, e.g., time-reversal-broken axion
insulators protected by spatial inversion symmetry [52, 20, 30, 27, 53]
generally have nonvanishing bulk magnetization. The Hopf insulator is another
axion insulator that is not topomagnetic, however certain symmetries render
its bulk magnetization to vanish in two of three independent spatial
directions. In this sense, the Hopf insulator realizes an anisotropic analog
of a topomagnetic insulator.

The simplest symmetry that would constrain the bulk magnetization to vanish is
time reversal. However, this symmetry is incompatible with the Hopf invariant
$\chi\neq 0$. This is because $\chi$, being an integral of the Chern-Simons
three-form [cf. Eq. (14)], transforms like a pseudoscalar under
crystallographic spatial transformations, and is odd under time reversal.

The best alternative one can achieve for the Hopf insulator is to impose
symmetry such that the bulk magnetization vanishes in two independent
directions. For example, time reversal $\hat{T}$ composed with the reflection
$\hat{m}_{x}:x\rightarrow-x,$ constrains the $x$-component of the
magnetization to vanish, and is also compatible with $\chi\neq 0$. It is
possible to simultaneously have $\hat{T}\hat{m}_{x}$ and $\hat{T}\hat{m}_{y}$
symmetries, hence forbidding magnetization in the $x$ and $y$-directions.
However, if $\hat{T}\hat{m}_{z}$ is further imposed, the product of these
three operations inverts $(x,y,z,t)\rightarrow(-x,-y,-z,-t)$, which ensures
the Berry curvature vanishes at each $\bm{k}$, hence $\chi=0$. Incidentally,
the MRW model [cf. Eq. (19)] is manifestly symmetric under
$\hat{T}\hat{m}_{x}$ and $\hat{T}\hat{m}_{y}$:

$\begin{array}[]{c}\tau_{z}H_{MRW}^{\star}(\bm{k})\tau_{z}=H_{MRW}(k_{x},-k_{y},-k_{z}),\\\
H_{MRW}^{\star}(\bm{k})=H_{MRW}(-k_{x},k_{y},-k_{z}),\end{array}$ (32)

making the MRW phase an anisotropic topomagnetic insulator.

As a curious sidenote, because the MRW model Hamiltonian is traceless, it
possesses an antisymmetry [cf. Eq. (29)] that ensures that the orbital (but
not the spin) magnetization vanishes in all directions. This follows from a
more general result: any three-dimensional insulator with trivial bulk Chern
invariants and the following antisymmetry

$\displaystyle{UH^{\star}(\mathbf{k})U^{{\dagger}}=-H(\mathbf{k}),\;\;\;\;U^{\scriptstyle{\dagger}}=U^{-1},}$
(35)

must have vanishing extensive orbital magnetization in the bulk. This
statement holds for Hamiltonians of any rank, and in particular it applies to
Pauli-matrix Hamiltonians, as proven in Appendix C.

##  IV Summary and outlook

We have investigated the subextensive magnetic moment $m$ contributed by
crystalline facets of a three-dimensional insulator in the trivial first Chern
class, assuming the existence of a bulk and surface energy gap. This magnetic
moment has both dynamical and geometric contributions, where the latter is a
product of the chemical potential $\mu$ with a momentum-space integral of the
Berry curvature (for states localized to the facet). Consequently, the
magnetic compressibility $A^{-1}dm/d\mu$ (contributed by one facet, per unit
facet area $A$) is $e/hc$ times the geometric part of the surface anomalous
Hall conductivity. In some cases the geometric contribution is restricted to a
half-integer multiple of $e/hc$, owing either to symmetry or to Hilbert-space
restrictions.

Geometry is ever linked with topology: two adjacent facets with differing
geometrical conductivity must host chiral hinge mode(s) at their interface,
where the net chirality is an integer-valued topological invariant. Thus we
establish a connection between quantized surface magnetism and higher-order
topology. When half-integer-quantization of $A^{-1}dm/d\mu$ is combined with a
vanishing bulk magnetization, such a phase of matter is referred to as a
topomagnetic insulator.

An advantage of topomagnetic insulator is that the faceted magnetization is
more readily measurable if the bulk contribution vanishes. Since the
magnetization contributed by a single facet has never been measured (to our
knowledge), we speculate on three existing experimental techniques that might
be augmented to do the job:

(i) Torque magnetometry couples to the total moment of a crystallite [54],
which is a sum of contributions from the bulk and from all facets. The former
contribution vanishes by assumption. Note that the contribution from opposite
facets of a crystallite generically do not cancel out, as exemplified by our
Hopf case study.

(ii) Magnetic neutron diffraction probes the magnetization [55] of a single
facet, assuming the cross-sectional area of the neutron beam is smaller than
the faceted area. While the penetration depth of a neutron beam is typically
much larger than the penetration depth of surface states, there is no need to
isolate the faceted contribution if the bulk magnetization vanishes.

(iii) Measuring second harmonic generation is also a surface-sensitive
diagnostic of magnetization in centrosymmetric materials [56].

The above three techniques potentially measure the faceted magnetization, but
to further measure the faceted compressibility, a means must be found to vary
the chemical potential, e.g., by doping [49], or gating.

## Acknowledgements

P.Z and T.L.H. were supported by ARO MURI W911NF2020166. A.A. was supported by
the Gordon and Betty Moore Foundation EPiQS Initiative through Grant No. GBMF
4305 and GBMF 8691 at the University of Illinois.

##  Appendix A Derivation of Eq.(3) from local orbital magnetization

###  A.1 Derivation of Eq.(3) under OBC

To derive Eq.(3) from Eq.(2) under OBC in a slab geometry where $z$ direction
is open and $x$ and $y$ directions are periodic, we introduce
$|\bm{k}_{\perp},R_{z},a\rangle=\sum_{\bm{R}_{\perp}}\sqrt{A_{xy,\mathrm{cell}}}e^{i\bm{k}_{\perp}\cdot\bm{R}_{\perp}}|\bm{R},a\rangle$,
where $A_{xy,\mathrm{cell}}$ is the area in the $x$-$y$ plane of a primitive
unit cell, and $\bm{k}_{\perp}=(k_{x},k_{y})$, and
$\bm{R}_{\perp}=(R_{x},R_{y})$. The normalization of
$|\bm{k}_{\perp},R_{z},a\rangle$ is
$\langle\bm{k}_{\perp},R_{z},a|\bm{k}_{\perp}^{\prime},R^{\prime}_{z},b\rangle=(2\pi)^{2}\delta^{(2)}(\bm{k}_{\perp}-\bm{k}^{\prime}_{\perp})\delta_{R_{z}R_{z}^{\prime}}\delta_{ab}$,
where $\delta^{(2)}(\bm{k}_{\perp}-\bm{k}^{\prime}_{\perp})$ is the two-
dimensional Dirac delta function. Note that

$\langle\bm{k}_{\perp},R_{z},a|\bar{O}|\bm{k}_{\perp}^{\prime},R_{z}^{\prime},b\rangle=(2\pi)^{2}\delta^{(2)}(\bm{k}_{\perp}-\bm{k}^{\prime}_{\perp})\widetilde{O}_{R_{z}R_{z}^{\prime},ab},\quad\mathrm{for}\quad
O=P,Q,H.$ (36)

Insert the identity
$(1/(2\pi)^{2})\sum_{R_{z},a}\int_{rBZ}d\bm{k}_{\perp}|\bm{k}_{\perp},R_{z},a\rangle\langle\bm{k}_{\perp},R_{z},a|=\mathbbm{1}$
into Eq.(2), and sum over $\bm{R}_{\perp}$ fixing $R_{z}=l$, we have

$\displaystyle\frac{1}{A}\sum_{\bm{R}_{\perp}}\mathfrak{m}_{z}(\bm{R}_{\perp},R_{z}=l)=\frac{1}{A}\sum_{\bm{R}_{\perp}}\frac{e}{\hbar
c}\operatorname{Im}\sum_{a,b,c}\sum_{R_{z}^{\prime},R_{z}^{\prime\prime}}\int_{rBZ}\frac{d\bm{k}_{\perp}}{(2\pi)^{2}}\frac{d\bm{k}^{\prime}_{\perp}}{(2\pi)^{2}}\langle\bm{R}_{\perp},l,a|\bm{k}_{\perp},R_{z}^{\prime},b\rangle\langle\bm{k}_{\perp},R_{z}^{\prime},b|$
(37)
$\displaystyle\times\bar{P}\hat{r}_{x}\bar{Q}\bar{H}\bar{Q}\hat{r}_{y}\bar{P}-\bar{Q}\hat{r}_{x}\bar{P}(\bar{H}-2\mu)\bar{P}\hat{r}_{y}\bar{Q}|\bm{k}^{\prime}_{\perp},R_{z}^{\prime\prime},c\rangle\langle\bm{k}^{\prime}_{\perp},R_{z}^{\prime\prime},c|\bm{R}_{\perp},l,a\rangle$
$\displaystyle=\frac{1}{A}\sum_{\bm{R}_{\perp}}A_{xy,\mathrm{cell}}\frac{e}{\hbar
c}\operatorname{Im}\int_{rBZ}\frac{d\bm{k}_{\perp}}{(2\pi)^{2}}\frac{d\bm{k}^{\prime}_{\perp}}{(2\pi)^{2}}e^{i(\bm{k}_{\perp}-\bm{k}^{\prime}_{\perp})\cdot\bm{R}_{\perp}}\langle\bm{k}_{\perp},l,a|\bar{P}\hat{r}_{x}\bar{Q}\bar{H}\bar{Q}\hat{r}_{y}\bar{P}-\bar{Q}\hat{r}_{x}\bar{P}(\bar{H}-2\mu)\bar{P}\hat{r}_{y}\bar{Q}|\bm{k}^{\prime}_{\perp},l,a\rangle,$

where $A$ is the area of the $x$-$y$ plane of the system, and
$A=\sum_{\bm{R}_{\perp}}A_{xy,\mathrm{cell}}$. With Eq.(36), follow steps
closely analogous to Eq.(43), we can derive

$\displaystyle\sum_{a}\langle\bm{k}_{\perp},l,a|\bar{P}\hat{r}_{x}\bar{Q}\bar{H}\bar{Q}\hat{r}_{y}\bar{P}|\bm{k}^{\prime}_{\perp},l,a\rangle$
(38)
$\displaystyle=(2\pi)^{2}\delta^{(2)}(\bm{k}_{\perp}-\bm{k}^{\prime}_{\perp})\sum_{R_{z}^{1},R_{z}^{2},R_{z}^{3},R_{z}^{4}}\operatorname{Tr}_{\mathrm{cell}}\left[\partial_{k_{x}}\widetilde{P}_{lR_{z}^{1}}\widetilde{Q}_{R_{z}^{1}R_{z}^{2}}\widetilde{H}_{R_{z}^{2}R_{z}^{3}}\widetilde{Q}_{R_{z}^{3}R_{z}^{4}}\partial_{k_{y}}\widetilde{P}_{R_{z}^{4}l}\right]$
$\displaystyle=(2\pi)^{2}\delta^{(2)}(\bm{k}_{\perp}-\bm{k}^{\prime}_{\perp})\operatorname{Tr}_{\mathrm{cell},z}\left[\mathcal{P}_{l}\partial_{k_{x}}\widetilde{P}\widetilde{Q}\widetilde{H}\widetilde{Q}\partial_{k_{y}}\widetilde{P}\right]=(2\pi)^{2}\delta^{(2)}(\bm{k}_{\perp}-\bm{k}^{\prime}_{\perp})g_{\bm{k}_{\perp},xy}(l),$

where $\mathcal{P}_{l}$ is the projector on the layer with $R_{z}=l$, and we
remind that $\operatorname{Tr}_{\mathrm{cell},z}$ is the trace respective to
the index $(R_{z},a)$. The last equality used the definition of
$g_{\bm{k}_{\perp},xy}(l)$ in Eq.(4). Similarly, we can also derive

$\displaystyle\langle\bm{k}_{\perp},l,a|\bar{Q}\hat{r}_{x}\bar{P}\bar{H}\bar{P}\hat{r}_{y}\bar{Q}|\bm{k}^{\prime}_{\perp},l,a\rangle=-(2\pi)^{2}\delta^{(2)}(\bm{k}_{\perp}-\bm{k}^{\prime}_{\perp})\operatorname{Tr}_{\mathrm{cell},z}\left[\mathcal{P}_{l}\widetilde{Q}\partial_{k_{x}}\widetilde{P}\widetilde{H}\widetilde{P}\partial_{k_{y}}\widetilde{Q}\right]$
(39)
$\displaystyle=(2\pi)^{2}\delta^{(2)}(\bm{k}_{\perp}-\bm{k}^{\prime}_{\perp})\operatorname{Tr}_{\mathrm{cell},z}\left[\mathcal{P}_{l}\widetilde{Q}\partial_{k_{x}}\widetilde{P}\widetilde{H}\partial_{k_{y}}\widetilde{P}\widetilde{Q}\right]=(2\pi)^{2}\delta^{(2)}(\bm{k}_{\perp}-\bm{k}^{\prime}_{\perp})h_{\bm{k}_{\perp},yx}(l),$

and

$\displaystyle\langle\bm{k}_{\perp},l,a|\bar{Q}\hat{r}_{x}\bar{P}\bar{P}\hat{r}_{y}\bar{Q}|\bm{k}^{\prime}_{\perp},l,a\rangle=-(2\pi)^{2}\delta^{(2)}(\bm{k}_{\perp}-\bm{k}^{\prime}_{\perp})\operatorname{Tr}_{\mathrm{cell},z}\left[\mathcal{P}_{l}\widetilde{Q}\partial_{k_{x}}\widetilde{P}P\partial_{k_{y}}\widetilde{Q}\right]$
(40)
$\displaystyle=(2\pi)^{2}\delta^{(2)}(\bm{k}_{\perp}-\bm{k}^{\prime}_{\perp})\operatorname{Tr}_{\mathrm{cell},z}\left[\mathcal{P}_{l}\widetilde{Q}\partial_{k_{x}}\widetilde{P}\partial_{k_{y}}\widetilde{P}\widetilde{Q}\right]=(2\pi)^{2}\delta^{(2)}(\bm{k}_{\perp}-\bm{k}^{\prime}_{\perp})f_{\bm{k}_{\perp},yx}(l).$

In Eq.(39) and Eq.(40), we used
$\widetilde{P}\partial\widetilde{Q}=-\partial\widetilde{P}\widetilde{Q}$; the
last steps of both Eq.(39) and Eq.(40) again used the definition of
$h_{\bm{k}_{\perp},yx}(l)$ and $f_{\bm{k}_{\perp},yx}(l)$ in Eq.(4) . Note
that from Eq.(4), it is straightforward to see
$\operatorname{Im}(\widetilde{o}_{\bm{k}_{\perp},xy})=-\operatorname{Im}(\widetilde{o}_{\bm{k}_{\perp},yx})$
for $o=g,h,f$. With this fact, we substitute Eq.(38), Eq.(39) and Eq.(40) into
Eq.(37), we can finally derive Eq.(3).

###  A.2 Equivalence of Eq.(2) and Eq.(56) under PBC

To prove Eq.(2) and Eq .(56) with $\gamma=z$ are equivalent under PBC, we
define $|\bm{k},a\rangle$ as the Fourier transformation of $|\bm{R},a\rangle$
such that
$|\bm{k},a\rangle=\sum_{\bm{R}}\sqrt{V_{\mathrm{cell}}}e^{i\bm{k}\cdot\bm{R}}|\bm{R},a\rangle$,
where $V_{\mathrm{cell}}$ is the volume of a primitive unit cell and satisfies
$\sum_{\bm{R}}V_{\mathrm{cell}}=V$. Then normalization of $|\bm{k},a\rangle$
is
$\langle\bm{k},a|\bm{k}^{\prime},b\rangle=(2\pi)^{3}\delta^{(3)}(\bm{k}-\bm{k}^{\prime})\delta_{ab}$,
where $\delta^{(3)}(\bm{k}-\bm{k}^{\prime})$ is the three-dimensional Dirac
delta function. Note that

$\langle\bm{k},a|\bar{O}|\bm{k}^{\prime},b\rangle=(2\pi)^{3}\delta^{(3)}(\bm{k}-\bm{k}^{\prime})O_{ab}(\bm{k}),\quad\mathrm{for}\quad
O=P,Q,H.$ (41)

We will show that

$\displaystyle(2\pi)^{3}\delta^{(3)}(\bm{k}-\bm{k}^{\prime})g_{\mathbf{k},xy}=(2\pi)^{3}\delta^{(3)}(\bm{k}-\bm{k}^{\prime})\operatorname{Tr}_{\mathrm{cell}}\left[\partial_{k_{x}}PQHQ\partial_{k_{y}}P\right]=\sum_{a}\langle\bm{k},a|\bar{P}\hat{r}_{x}\bar{Q}\bar{H}\bar{Q}\hat{r}_{y}\bar{P}|\bm{k}^{\prime},a\rangle.$
(42)

We remind that $\operatorname{Tr}_{\mathrm{cell}}$ is the trace respective to
the unit-cell index $a$. We prove the above equation by inserting the identity
$(1/(2\pi)^{3})\sum_{a}\int_{BZ}d\bm{k}|\bm{k},a\rangle\langle\bm{k},a|=\mathbbm{1}$
into the RHS of the above equation:

$\displaystyle\sum_{a}\langle\bm{k},a|\bar{P}\hat{r}_{x}\bar{Q}\bar{H}\bar{Q}\hat{r}_{y}\bar{P}|\bm{k}^{\prime},a\rangle=\int_{BZ}\frac{d\mathbf{k}^{1}}{(2\pi)^{3}}\frac{d\mathbf{k}^{3}}{(2\pi)^{3}}\frac{d\mathbf{k}^{3}}{(2\pi)^{3}}\frac{d\mathbf{k}^{4}}{(2\pi)^{3}}\sum_{a,b,c,d,e}\langle\bm{k},a|\bar{P}|\bm{k}^{1},b\rangle\langle\bm{k}^{1},b|\hat{r}_{x}\bar{Q}|\bm{k}^{2},c\rangle\langle\bm{k}^{2},c|$
(43)
$\displaystyle\quad\times\bar{H}|\bm{k}^{3},d\rangle\langle\bm{k}^{3},d|\bar{Q}|\bm{k}^{4},e\rangle\langle\bm{k}^{4},e|\hat{r}_{y}\bar{P}|\bm{k}^{\prime},a\rangle$
$\displaystyle=(2\pi)^{3}\int_{BZ}d\mathbf{k}^{1}d\mathbf{k}^{2}d\mathbf{k}^{3}d\mathbf{k}^{4}\operatorname{Tr}_{\mathrm{cell}}\big{\\{}P(\bm{k})\delta^{(3)}(\bm{k}-\bm{k}^{1})i\partial_{k^{1}_{x}}\left[Q(\bm{k}^{1})\delta^{(3)}(\bm{k}^{1}-\bm{k}^{2})\right]H(\bm{k}^{2})\delta^{(3)}(\bm{k}^{2}-\bm{k}^{3})Q(\bm{k}^{3})$
$\displaystyle\quad\times\delta^{(3)}(\bm{k}^{3}-\bm{k}^{4})i\partial_{k^{4}_{y}}\left[P(\bm{k}^{4})\delta^{(3)}(\bm{k}^{4}-\bm{k}^{\prime})\right]\big{\\}}$
$\displaystyle=(2\pi)^{3}\int_{BZ}d\mathbf{k}^{1}d\mathbf{k}^{2}d\mathbf{k}^{3}d\mathbf{k}^{4}\delta^{(3)}(\bm{k}-\bm{k}^{1})\delta^{(3)}(\bm{k}^{1}-\bm{k}^{2})\delta^{(3)}(\bm{k}^{2}-\bm{k}^{3})\delta^{(3)}(\bm{k}^{3}-\bm{k}^{4})\delta^{(3)}(\bm{k}^{4}-\bm{k}^{\prime})$
$\displaystyle\quad\times\operatorname{Tr}_{\mathrm{cell}}\left[P(\bm{k})i\partial_{k^{1}_{x}}Q(\bm{k}^{1})H(\bm{k}^{2})Q(\bm{k}^{3})i\partial_{k^{4}_{y}}P(\bm{k}^{4})\right]$
$\displaystyle=-(2\pi)^{3}\delta^{(3)}(\bm{k}-\bm{k}^{\prime})\operatorname{Tr}_{\mathrm{cell}}\left[P(\bm{k})\partial_{k_{x}}Q(\bm{k})H(\bm{k})Q(\bm{k})\partial_{k_{y}}P(\bm{k})\right]$
$\displaystyle=(2\pi)^{3}\delta^{(3)}(\bm{k}-\bm{k}^{\prime})\operatorname{Tr}_{\mathrm{cell}}\left[\partial_{k_{x}}P(\bm{k})Q(\bm{k})H(\bm{k})Q(\bm{k})\partial_{k_{y}}P(\bm{k})\right]=(2\pi)^{3}\delta^{(3)}(\bm{k}-\bm{k}^{\prime})g_{\bm{k},xy}$

From the second equality to the third equality in the above equation, we used
the fact that $P(\bm{k})\delta^{(3)}(\bm{k}-\bm{k}^{1})Q(\bm{k}^{1})=0$, so
that we have

$\displaystyle
P(\bm{k})\delta^{(3)}(\bm{k}-\bm{k}^{1})i\partial_{k^{1}_{x}}\left[Q(\bm{k}^{1})\delta^{(3)}(\bm{k}^{1}-\bm{k}^{2})\right]$
(44)
$\displaystyle=P(\bm{k})\delta^{(3)}(\bm{k}-\bm{k}^{1})\left\\{\left[i\partial_{k^{1}_{x}}Q(\bm{k}^{1})\right]\delta^{(3)}(\bm{k}^{1}-\bm{k}^{2})+Q(\bm{k}^{1})\left[i\partial_{k^{1}_{x}}\delta^{(3)}(\bm{k}^{1}-\bm{k}^{2})\right]\right\\}$
$\displaystyle=P(\bm{k})\delta_{k_{z},k_{z}^{1}}\delta^{(3)}(\bm{k}-\bm{k}^{1})\left[i\partial_{k^{1}_{x}}Q(\bm{k}_{1})\right]\delta^{(3)}(\bm{k}^{1}-\bm{k}^{2}),$

and

$\displaystyle
Q(\bm{k}^{3})\delta^{(3)}(\bm{k}^{3}-\bm{k}^{4})i\partial_{k^{4}_{y}}\left[P(\bm{k}^{4})\delta^{(3)}(\bm{k}^{4}-\bm{k}^{\prime})\right]$
(45)
$\displaystyle=Q(\bm{k}^{3})\delta^{(3)}(\bm{k}^{3}-\bm{k}^{4})\left\\{\left[i\partial_{k^{4}_{y}}P(\bm{k}^{4})\right]\delta^{(3)}(\bm{k}^{4}-\bm{k}^{\prime})+P(\bm{k}^{4})\left[i\partial_{k^{4}_{y}}\delta^{(3)}(\bm{k}^{4}-\bm{k}^{\prime})\right]\right\\}$
$\displaystyle=Q(\bm{k}^{3})\delta^{(3)}(\bm{k}^{3}-\bm{k}^{4})\left[i\partial_{k^{4}_{y}}P(\bm{k}^{4})\right]\delta^{(3)}(\bm{k}^{4}-\bm{k}^{\prime}).$

Following the same procedure, we can also get

$\begin{array}[]{c}(2\pi)^{3}\delta^{(3)}(\bm{k}-\bm{k}^{\prime})h_{\mathbf{k},yx}=\sum_{a}\langle\bm{k},a|\bar{Q}\hat{r}_{x}\bar{P}\bar{H}\bar{P}\hat{r}_{y}\bar{Q}|\bm{k}^{\prime},a\rangle,\\\
(2\pi)^{3}\delta^{(3)}(\bm{k}-\bm{k}^{\prime})f_{\mathbf{k},yx}=\sum_{a}\langle\bm{k},a|\bar{Q}\hat{r}_{x}\bar{P}\bar{P}\hat{r}_{y}\bar{Q}|\bm{k}^{\prime},a\rangle.\end{array}$
(46)

Then, by inserting the identity
$(1/(2\pi)^{3})\sum_{a}\int_{BZ}d\bm{k}|\bm{k},a\rangle\langle\bm{k},a|=\mathbbm{1}$
into Eq.(2), we have

$\displaystyle
M_{z}=\frac{1}{V}\sum_{\bm{R}}\mathfrak{m}_{z}(\bm{R})=\frac{1}{V}\sum_{\bm{R}}\frac{e}{\hbar
c}\operatorname{Im}\sum_{a,b,c}\int_{BZ}\frac{d\bm{k}}{(2\pi)^{3}}\frac{d\bm{k}^{\prime}}{(2\pi)^{3}}\langle\bm{R},a|\bm{k},b\rangle\langle\bm{k},b|\bar{P}\hat{r}_{x}\bar{Q}\bar{H}\bar{Q}\hat{r}_{y}\bar{P}$
(49)
$\displaystyle\quad-\bar{Q}\hat{r}_{x}\bar{P}(\bar{H}-2\mu)\bar{P}\hat{r}_{y}\bar{Q}|\bm{k}^{\prime},c\rangle\langle\bm{k}^{\prime},c|\bm{R},a\rangle$
$\displaystyle=\frac{1}{V}\sum_{\bm{R}}\frac{e}{\hbar
c}\operatorname{Im}\sum_{a,b,c}V_{\mathrm{cell}}\delta_{ab}\delta_{ca}\int_{BZ}\frac{d\bm{k}}{(2\pi)^{3}}\frac{d\bm{k}^{\prime}}{(2\pi)^{3}}e^{i(\bm{k}-\bm{k}^{\prime})\cdot\bm{R}}\langle\bm{k},b|\bar{P}\hat{r}_{x}\bar{Q}\bar{H}\bar{Q}\hat{r}_{y}\bar{P}-\bar{Q}\hat{r}_{x}\bar{P}(\bar{H}-2\mu)\bar{P}\hat{r}_{y}\bar{Q}|\bm{k}^{\prime},c\rangle$
$\displaystyle=\frac{1}{V}\sum_{\bm{R}}\frac{e}{\hbar
c}\operatorname{Im}\sum_{a}V_{\mathrm{cell}}\int_{BZ}\frac{d\bm{k}}{(2\pi)^{3}}\frac{d\bm{k}^{\prime}}{(2\pi)^{3}}e^{i(\bm{k}-\bm{k}^{\prime})\cdot\bm{R}}\langle\bm{k},a|\bar{P}\hat{r}_{x}\bar{Q}\bar{H}\bar{Q}\hat{r}_{y}\bar{P}-\bar{Q}\hat{r}_{x}\bar{P}(\bar{H}-2\mu)\bar{P}\hat{r}_{y}\bar{Q}|\bm{k}^{\prime},a\rangle$
$\displaystyle=\frac{1}{V}\sum_{\bm{R}}\frac{e}{\hbar
c}\operatorname{Im}V_{\mathrm{cell}}\int_{BZ}\frac{d\bm{k}}{(2\pi)^{3}}d\bm{k}^{\prime}e^{i(\bm{k}-\bm{k}^{\prime})\cdot\bm{R}}\delta^{(3)}(\bm{k}-\bm{k}^{\prime})\left(g_{\bm{k},xy}-h_{\bm{k},yx}+2\mu
f_{\bm{k},yx}\right)$ $\displaystyle=\frac{e}{\hbar
c}\frac{1}{V}\sum_{\bm{R}}V_{\mathrm{cell}}\operatorname{Im}\int_{BZ}\frac{d\bm{k}}{(2\pi)^{3}}\left(g_{\bm{k},xy}-h_{\bm{k},yx}+2\mu
f_{\bm{k},yx}\right)$ $\displaystyle=\frac{e}{\hbar
c}\operatorname{Im}\int_{BZ}\frac{d\bm{k}}{(2\pi)^{3}}\left(g_{\bm{k},xy}-h_{\bm{k},yx}+2\mu
f_{\bm{k},yx}\right)$

Note that from Eq.(57), it is straightforward to see
$o_{\bm{k},xy}=o^{\star}_{\bm{k},yx}$ for $o=g,h,f$, which leads to
$\operatorname{Im}(o_{\bm{k},xy})=-\operatorname{Im}(o^{\star}_{\bm{k},yx})$.
With this fact, the above equation becomes

$\displaystyle
M_{z}=\frac{1}{V}\sum_{\bm{R}}\mathfrak{m}_{z}(\bm{R})=\frac{e}{\hbar
c}\operatorname{Im}\int_{BZ}\frac{d\bm{k}}{(2\pi)^{3}}\left(g_{\bm{k},xy}+h_{\bm{k},xy}-2\mu
f_{\bm{k},xy}\right).$ (50)

If we fix $\gamma=z$ and $d=3$ in Eq.(56), we have

$\displaystyle M_{z}$
$\displaystyle=\frac{e}{c\hbar}\operatorname{Im}\int_{BZ}\frac{d\mathbf{k}}{(2\pi)^{3}}\left\\{\epsilon_{zxy}\bigg{[}\frac{1}{2}(g_{\mathbf{k},xy}+h_{\mathbf{k},xy})-\mu
f_{\mathbf{k},xy}\bigg{]}+\epsilon_{zyx}\bigg{[}\frac{1}{2}(g_{\mathbf{k},yx}+h_{\mathbf{k},yx})-\mu
f_{\mathbf{k},yx}\bigg{]}\right\\}$ (51)
$\displaystyle=\frac{e}{c\hbar}\operatorname{Im}\int_{BZ}\frac{d\mathbf{k}}{(2\pi)^{3}}\left\\{\bigg{[}\frac{1}{2}(g_{\mathbf{k},xy}+h_{\mathbf{k},xy})-\mu
f_{\mathbf{k},xy}\bigg{]}-\bigg{[}\frac{1}{2}(g_{\mathbf{k},yx}+h_{\mathbf{k},yx})-\mu
f_{\mathbf{k},yx}\bigg{]}\right\\}$
$\displaystyle=\frac{e}{c\hbar}\operatorname{Im}\int_{BZ}\frac{d\mathbf{k}}{(2\pi)^{3}}(g_{\mathbf{k},xy}+h_{\mathbf{k},xy}-2\mu
f_{\mathbf{k},xy}),$

which is clearly equivalent to Eq.(50).

##  Appendix B Transformation of the magnetization under crystallographic
spacetime group

If the bulk Chern class is trivial, the bulk orbital magnetization has
vanishing compressibility, and is simply related to a trace of the orbital
angular momentum operator:
$\mathbf{M}=(-e/2Vc)\textit{tr}(\bar{P}\hat{\bm{r}}\times\hat{\bm{v}})$.
Note“tr” is the trace over the whole Hilbert space, $V$ is the volume of the
system, and $\bar{P}$ is the projector onto the occupied subspace of a
Hamiltonian $\bar{H}$ in the real space basis. The $\hat{\bm{r}}$ and
$\hat{\bm{v}}$ are position operator and velocity operator respectively.

We show here that $\mathbf{M}$ transforms as a pseudo-vector under the
crystallographic spatial transformations and is moreover odd under time
reversal. Consider an element $g$ in the spacetime symmetry group such that

$\hat{g}\hat{r}_{i}\hat{g}^{-1}=\breve{g}_{ij}\hat{r}_{j}+\delta_{i},\breve{g}^{-1}=\breve{g}^{t}\in\mathbb{R}.$
(52)

If a Hamiltonian has the symmetry $\hat{g}$, $i.e.$,
$\hat{g}\bar{H}\hat{g}^{-1}=\bar{H}$ and thus
$\hat{g}\bar{P}\hat{g}^{-1}=\bar{P}$, given
$\hat{v}_{l}=i[\bar{H},\hat{r}_{l}]/\hbar$, then we have

$\displaystyle\hat{g}\hat{v}_{l}\hat{g}^{-1}$
$\displaystyle=(-1)^{s(g)}(\breve{g}_{lj}\frac{i}{\hbar}[\bar{H},\hat{r}_{j}])$
(53) $\displaystyle=(-1)^{s(g)}(\breve{g}_{lj}\hat{v}_{j}),$

where $s(g)=1$ if there is time-reversal in $g$, and otherwise $s(g)=0$.
Substitute Eq.(52) and Eq.(53) into
$\mathbf{M}=(-e/2Vc)\textit{tr}(\bar{P}\hat{\bm{r}}\times\hat{\bm{v}})$, and
note that $\textit{tr}(\bar{P}\hat{v}_{l})=0$ because the total current is
zero in insulators under zero external field, we have

$\displaystyle\mathbf{M}_{i}$
$\displaystyle=\frac{e}{2Vc}\epsilon_{ijl}\textit{tr}(\hat{g}\bar{P}\hat{g}^{-1}\hat{g}\hat{r}_{j}\hat{g}^{-1}\hat{g}\hat{v}_{l}\hat{g}^{-1})$
(54)
$\displaystyle=-\frac{e}{2Vc}(-1)^{s(g)}\epsilon_{ijl}\breve{g}_{jm}\breve{g}_{ln}\textit{tr}(\bar{P}\hat{r}_{m}\hat{v}_{n}).$

Multiply
$\delta_{ip}=\breve{g}_{iq}\breve{g}^{-1}_{qp}=\breve{g}_{iq}\breve{g}_{pq}$
on both side of the above equation, given
$\epsilon_{ijl}\breve{g}_{iq}\breve{g}_{jm}\breve{g}_{ln}=\det{\breve{g}}\epsilon_{qmn}$,
we can finally derive

$\mathbf{M}=(-1)^{s(g)}\det{\breve{g}}(\breve{g}\mathbf{M}).$ (55)

Note $\det{\breve{g}}=-1$ for improper rotations, hence $\mathbf{M}$ is
invariant under spatial inversion.

##  Appendix C Vanishing bulk orbital magnetization due to the antisymmetry

To prove that the bulk orbital magnetization can be ruled out by the
antisymmetry in Eq. (29), let us first rewrite the orbital magnetization in
the momentum representation:[2, 3]

$\displaystyle M_{\gamma}$
$\displaystyle=\frac{e}{c\hbar}\operatorname{Im}\int_{BZ}\frac{d\mathbf{k}}{(2\pi)^{d}}\epsilon_{\alpha\beta\gamma}\bigg{[}\frac{1}{2}(g_{\mathbf{k},\alpha\beta}+h_{\mathbf{k},\alpha\beta})-\mu
f_{\mathbf{k},\alpha\beta}\bigg{]},$ (56)

where

$\begin{array}[]{c}g_{\mathbf{k},\alpha\beta}=\operatorname{Tr}_{\mathrm{cell}}\left[\partial_{\alpha}PQHQ\partial_{\beta}P\right],\\\
h_{\mathbf{k},\alpha\beta}=\operatorname{Tr}_{\mathrm{cell}}\left[Q\partial_{\beta}PH\partial_{\alpha}PQ\right],\\\
f_{\mathbf{k},\alpha\beta}=\operatorname{Tr}_{\mathrm{cell}}\left[Q\partial_{\beta}P\partial_{\alpha}PQ\right],\end{array}$
(57)

$H(\mathbf{k})$ is a Bloch Hamiltonian with an energy gap separating occupied
and unoccupied bands, which are respectively projected by $P(\bm{k})$ and
$Q(\bm{k})$. (In the above equation, the arguments of $P,Q,H$ have been
omitted for simplicity.) All of $P(\bm{k}),Q(\bm{k})$ and $H(\bm{k})$ are
matrices with row index $a$ labelling the basis vectors of the Hilbert space
of one primitive unit cell, and “$\operatorname{Tr}_{\mathrm{cell}}$” is the
trace with respect to the unit-cell index $a$. If $H(\mathbf{k})$ is an
insulating Bloch Hamiltonian satisfying Eq. (35), then
$UP^{\star}(\mathbf{k})U^{{\dagger}}=Q(\mathbf{k})$. Substituting this into
Eq.(57), we have

$\displaystyle
h_{\mathbf{k},\alpha\beta}=\operatorname{Tr}_{\mathrm{cell}}\left[\partial_{\beta}PH\partial_{\alpha}PQ\right]=\operatorname{Tr}_{\mathrm{cell}}\left[\partial_{\beta}PH\partial_{\alpha}PUP^{\star}U^{{\dagger}}\right]$
(61)
$\displaystyle=\operatorname{Tr}_{\mathrm{cell}}[\partial_{\beta}PUU^{{\dagger}}HUU^{{\dagger}}\partial_{\alpha}PUP^{\star}U^{{\dagger}}]$
$\displaystyle=\operatorname{Tr}_{\mathrm{cell}}[(U^{{\dagger}}\partial_{\beta}PU)(U^{{\dagger}}HU)(U^{{\dagger}}\partial_{\alpha}PU)P^{\star}]$
$\displaystyle=-\operatorname{Tr}_{\mathrm{cell}}\left[\partial_{\beta}Q^{\star}H^{\star}\partial_{\alpha}Q^{\star}P^{\star}\right]$
$\displaystyle=-\operatorname{Tr}_{\mathrm{cell}}\left[\partial_{\beta}P^{\star}H^{\star}\partial_{\alpha}P^{\star}P^{\star}\right]$
$\displaystyle=-\operatorname{Tr}_{\mathrm{cell}}\left[P\partial_{\alpha}PH\partial_{\beta}P\right].$

For $g_{\mathbf{k},\alpha\beta}$ in Eq.(57), using
$\partial_{\alpha}PQ=P\partial_{\alpha}P$ and
$Q\partial_{\beta}P=\partial_{\beta}PP$ we can do some manipulations:

$\displaystyle
g_{\mathbf{k},\alpha\beta}=\operatorname{Tr}_{\mathrm{cell}}\left[\partial_{\alpha}PQHQ\partial_{\beta}P\right]=\operatorname{Tr}_{\mathrm{cell}}\left[P\partial_{\alpha}PH\partial_{\beta}PP\right]$
(62)
$\displaystyle=\operatorname{Tr}_{\mathrm{cell}}\left[P\partial_{\alpha}PH\partial_{\beta}P\right].$

Compare the above two equations, we can clearly see
$h_{\bm{k},\alpha\beta}=-g_{\bm{k},\alpha\beta}$. Thus, the first two terms of
Eq.(56) cancel each other. The third and last term is proportional to the
bulk-BZ-integral of $f_{\mathbf{k},\alpha\beta}$, and therefore vanishes along
with the bulk Chern invariants, completing our proof.

## References

  * Thonhauser  _et al._ [2005] T. Thonhauser, D. Ceresoli, D. Vanderbilt, and R. Resta, Orbital magnetization in periodic insulators, Phys. Rev. Lett. 95, 137205 (2005).
  * Ceresoli  _et al._ [2006] D. Ceresoli, T. Thonhauser, D. Vanderbilt, and R. Resta, Orbital magnetization in crystalline solids: Multi-band insulators, chern insulators, and metals, Phys. Rev. B 74, 024408 (2006).
  * Shi  _et al._ [2007] J. Shi, G. Vignale, D. Xiao, and Q. Niu, Quantum theory of orbital magnetization and its generalization to interacting systems, Phys. Rev. Lett. 99, 197202 (2007).
  * Thonhauser [2011] T. Thonhauser, Theory of orbital magnetization in solids, International Journal of Modern Physics B 25, 1429 (2011).
  * Bianco and Resta [2016] R. Bianco and R. Resta, Orbital magnetization in insulators: Bulk versus surface, Phys. Rev. B 93, 174417 (2016).
  * Bianco and Resta [2013] R. Bianco and R. Resta, Orbital magnetization as a local property, Phys. Rev. Lett. 110, 087202 (2013).
  * Schulz-Baldes and Teufel [2013] H. Schulz-Baldes and S. Teufel, Orbital polarization and magnetization for independent particles in disordered media, Communications in Mathematical Physics 319, 649 (2013).
  * Langbehn  _et al._ [2017] J. Langbehn, Y. Peng, L. Trifunovic, F. von Oppen, and P. W. Brouwer, Reflection-symmetric second-order topological insulators and superconductors, Phys. Rev. Lett. 119, 246401 (2017).
  * Benalcazar  _et al._ [2017] W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, Electric multipole moments, topological multipole moment pumping, and chiral hinge states in crystalline insulators, Physical Review B 96, 245115 (2017).
  * Schindler  _et al._ [2018a] F. Schindler, A. M. Cook, M. G. Vergniory, Z. Wang, S. S. Parkin, B. A. Bernevig, and T. Neupert, Higher-order topological insulators, Science advances 4, eaat0346 (2018a).
  * Schindler  _et al._ [2018b] F. Schindler, Z. Wang, M. G. Vergniory, A. M. Cook, A. Murani, S. Sengupta, A. Y. Kasumov, R. Deblock, S. Jeon, I. Drozdov, _et al._ , Higher-order topology in bismuth, Nature physics 14, 918 (2018b).
  * Kooi  _et al._ [2018] S. H. Kooi, G. van Miert, and C. Ortix, Inversion-symmetry protected chiral hinge states in stacks of doped quantum hall layers, Phys. Rev. B 98, 245102 (2018).
  * Geier  _et al._ [2018] M. Geier, L. Trifunovic, M. Hoskam, and P. W. Brouwer, Second-order topological insulators and superconductors with an order-two crystalline symmetry, Phys. Rev. B 97, 205135 (2018).
  * Trifunovic and Brouwer [2019] L. Trifunovic and P. W. Brouwer, Higher-order bulk-boundary correspondence for topological crystalline phases, Phys. Rev. X 9, 011012 (2019).
  * Moore  _et al._ [2008] J. E. Moore, Y. Ran, and X.-G. Wen, Topological surface states in three-dimensional magnetic insulators, Phys. Rev. Lett. 101, 186805 (2008).
  * Deng  _et al._ [2013] D.-L. Deng, S.-T. Wang, C. Shen, and L.-M. Duan, Hopf insulators and their topologically protected surface states, Phys. Rev. B 88, 201105 (2013).
  * Hopf [1964] H. Hopf, Über die abbildungen der dreidimensionalen sphäre auf die kugelfläche, in  _Selecta Heinz Hopf_ (Springer, 1964) pp. 38–63.
  * Mattis [2006] D. C. Mattis,  _The theory of magnetism made simple: an introduction to physical concepts and to some useful mathematical methods_ (World Scientific Publishing Company, 2006).
  * Fu  _et al._ [2007] L. Fu, C. L. Kane, and E. J. Mele, Topological insulators in three dimensions, Phys. Rev. Lett. 98, 106803 (2007).
  * 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).
  * Sitte  _et al._ [2012] M. Sitte, A. Rosch, E. Altman, and L. Fritz, Topological insulators in magnetic fields: Quantum hall effect and edge channels with a nonquantized $\theta$ term, Phys. Rev. Lett. 108, 126807 (2012).
  * Cloizeaux [1964] J. D. Cloizeaux, Energy bands and projection operators in a crystal: Analytic and asymptotic properties, Phys. Rev. 135, A685 (1964).
  * Thonhauser and Vanderbilt [2006] T. Thonhauser and D. Vanderbilt, Insulator/chern-insulator transition in the haldane model, Phys. Rev. B 74, 235111 (2006).
  * Bianco and Resta [2011] R. Bianco and R. Resta, Mapping topological order in coordinate space, Phys. Rev. B 84, 241106 (2011).
  * Streda [1982] P. Streda, Theory of quantised hall conductivity in two dimensions, Journal of Physics C: Solid State Physics 15, L717 (1982).
  * Xiao  _et al._ [2010] D. Xiao, M.-C. Chang, and Q. Niu, Berry phase effects on electronic properties, Rev. Mod. Phys. 82, 1959 (2010).
  * 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).
  * Rauch  _et al._ [2018] T. c. v. Rauch, T. Olsen, D. Vanderbilt, and I. Souza, Geometric and nongeometric contributions to the surface anomalous hall conductivity, Phys. Rev. B 98, 115108 (2018).
  * Jackson [2007] J. D. Jackson,  _Classical electrodynamics_ (John Wiley & Sons, 2007).
  * Varnava and Vanderbilt [2018] N. Varnava and D. Vanderbilt, Surfaces of axion insulators, Phys. Rev. B 98, 245117 (2018).
  * Olsen  _et al._ [2020] T. Olsen, T. c. v. Rauch, D. Vanderbilt, and I. Souza, Gapless hinge states from adiabatic pumping of axion coupling, Phys. Rev. B 102, 035166 (2020).
  * Alexandradinata  _et al._ [2019] A. Alexandradinata, A. Nelson, and A. A. Soluyanov, The actually robust surface signature of a hopf insulator: Bulk-to-boundary flow of berry curvature beyond the anomaly inflow paradigm, arXiv preprint arXiv:1910.10717 (2019).
  * Kennedy [2016] R. Kennedy, Topological hopf-chern insulators and the hopf superconductor, Phys. Rev. B 94, 035137 (2016).
  * Essin  _et al._ [2010] A. M. Essin, A. M. Turner, J. E. Moore, and D. Vanderbilt, Orbital magnetoelectric coupling in band insulators, Phys. Rev. B 81, 205104 (2010).
  * Malashevich  _et al._ [2010] A. Malashevich, I. Souza, S. Coh, and D. Vanderbilt, Theory of orbital magnetoelectric response, New Journal of Physics 12, 053032 (2010).
  * Coh  _et al._ [2011] S. Coh, D. Vanderbilt, A. Malashevich, and I. Souza, Chern-simons orbital magnetoelectric coupling in generic insulators, Phys. Rev. B 83, 085108 (2011).
  * Yu  _et al._ [2011] R. Yu, X. L. Qi, A. Bernevig, Z. Fang, and X. Dai, Equivalent expression of $\mathbbm{Z}_{2}$ topological invariant for band insulators using the non-abelian berry connection, Phys. Rev. B 84, 075119 (2011).
  * Vanderbilt [2018] D. Vanderbilt,  _Berry Phases in Electronic Structure Theory: Electric Polarization, Orbital Magnetization and Topological Insulators_ (Cambridge University Press, 2018).
  * Taherinejad  _et al._ [2014] M. Taherinejad, K. F. Garrity, and D. Vanderbilt, Wannier center sheets in topological insulators, Phys. Rev. B 89, 115102 (2014).
  * Soluyanov and Vanderbilt [2011] A. A. Soluyanov and D. Vanderbilt, Wannier representation of $\mathbbm{Z}_{2}$ topological insulators, Phys. Rev. B 83, 035108 (2011).
  * Wannier [1960] G. H. Wannier, Wave functions and effective Hamiltonian for Bloch electrons in an electric field, Phys. Rev. 117, 432 (1960).
  * Taherinejad and Vanderbilt [2015] M. Taherinejad and D. Vanderbilt, Adiabatic pumping of chern-simons axion coupling, Phys. Rev. Lett. 114, 096401 (2015).
  * Olsen  _et al._ [2017] T. Olsen, M. Taherinejad, D. Vanderbilt, and I. Souza, Surface theorem for the chern-simons axion coupling, Phys. Rev. B 95, 075137 (2017).
  * Stone and Goldbart [2008] M. Stone and P. Goldbart,  _Mathematics for Physics_ (Pimander-Casaubon, Alexandria, 2008).
  * Hatcher [2009] A. Hatcher,  _Vector Bundles and K-Theory_ (Autoedición, 2009).
  * Armitage and Wu [2019] N. P. Armitage and L. Wu, On the matter of topological insulators as magnetoelectrics, SciPost Phys. 6, 46 (2019).
  * Laughlin [1981] R. B. Laughlin, Quantized hall conductivity in two dimensions, Phys. Rev. B 23, 5632 (1981).
  * Liu  _et al._ [2017] C. Liu, F. Vafa, and C. Xu, Symmetry-protected topological hopf insulator and its generalizations, Phys. Rev. B 95, 161116 (2017).
  * Chen  _et al._ [2010] Y. Chen, J.-H. Chu, J. Analytis, Z. Liu, K. Igarashi, H.-H. Kuo, X. Qi, S.-K. Mo, R. Moore, D. Lu, _et al._ , Massive dirac fermion on the surface of a magnetically doped topological insulator, Science 329, 659 (2010).
  * Mong  _et al._ [2010] R. S. K. Mong, A. M. Essin, and J. E. Moore, Antiferromagnetic topological insulators, Phys. Rev. B 81, 245209 (2010).
  * Xu  _et al._ [2019] Y. Xu, Z. Song, Z. Wang, H. Weng, and X. Dai, Higher-order topology of the axion insulator ${\mathrm{euin}}_{2}{\mathrm{as}}_{2}$, Phys. Rev. Lett. 122, 256402 (2019).
  * Hughes  _et al._ [2011] T. L. Hughes, E. Prodan, and B. A. Bernevig, Inversion-symmetric topological insulators, Phys. Rev. B 83, 245132 (2011).
  * Fang  _et al._ [2012] C. Fang, M. J. Gilbert, and B. A. Bernevig, Bulk topological invariants in noninteracting point group symmetric insulators, Phys. Rev. B 86, 115112 (2012).
  * Shoenberg [1984] D. Shoenberg,  _Magnetic oscillations in metals_ (Cambridge University Press, The Edinburgh Building, Cambridge CB2 8RU, UK, 1984).
  * Hirst [1997] L. L. Hirst, The microscopic magnetization: concept and application, Rev. Mod. Phys. 69, 607 (1997).
  * Kirilyuk and Rasing [2005] A. Kirilyuk and T. Rasing, Magnetization-induced-second-harmonic generation from surfaces and interfaces, J. Opt. Soc. Am. B 22, 148 (2005).