††thanks: On the leave from Institut für Theorie der Statistischen Physik,
RWTH Aachen, 52056 Aachen, Germany

# Topological charge quantization on localized imperfections in crystalline
insulators and the nearsightedness principle of Kohn

Kiryl Piasotski<EMAIL_ADDRESS>Institut für Theorie der
Kondensierten Materie, Karlsruher Institut für Technologie, 76131 Karlsruhe,
Germany  Institut für QuantenMaterialien und Technologien, Karlsruher Institut
für Technologie, 76021 Karlsruhe, Germany

###### Abstract

We study the quantization of the excess charge on $N$ localized (ultra-
screened) impurities in $d$-dimensional crystalline insulating systems.
Solving Dyson’s equation, we demonstrate that such charges are topological, by
expressing them as winding numbers of appropriate functionals of bulk position
space Green’s functions. We discuss the ties of our topological invariant with
the nearsightedness principle of W. Kohn, stating that the electronic charge
density at fixed chemical potential depends on the external field only
locally, meaning that localized perturbations by external fields may only
result in localized charge redistributions. We arrive at the same conclusion
by demonstrating that an adiabatic perturbation comprised of a variation of
impurities’ positions and/or strengths may only result in the change in the
occupancy of impurity-localized bound states sitting, energy-wise, close to
the Fermi level. Finally, we conclude by discussing the relations of the
nearsightedness principle with the topological invariants characterizing the
boundary charge.

##  I Introduction

With the discovery of the quantum Hall effect [Klitzing_1980, ] and its
topological origins [Thouless_1982, , Hatsugai_1993, ] the study of the
topological structures in condensed matter systems became a style rather than
a fashion. The possibly best-known contemporary example of that is the field
of topological insulators, where the boundary-localized edge states, with
topologically granted existence and robustness, are being studied (popular
reviews are [Hasan_2010, , Qi_2011, ]). Despite being a well-defined endeavor
with its very own periodic table [Kitaev_2009, ], this discipline leaves a
number of questions open. One of these regards the direct experimental
accessibility of these topological surface states. In particular, aside from
the subspace of such states, their basis is incomplete for a description of
the physical system accommodating them – a topological insulator, making it
highly questionable whether an actual physical observable may be expanded into
a basis of intra- and inter-surface state transition operators. This is, for
example, not true of the excess charge density that these insulators
accumulate at their boundaries as it, as such, also features the exponentially
localized contributions of all of the occupied extended states. Despite that,
it is clear that as the surface states also contribute to such an observable,
the change in their occupancy has to have an observable effect.

In a series of recent works [Pletyukhov_2020, , Muller_2021, , Miles_2021, ,
Pias_2022, ], the topological properties of the boundary-localized electronic
excess charges (the boundary charges) in unidimensional crystals were
examined. In particular, a pair of topological invariants characterizing the
boundary charge upon two bulk energy spectrum-preserving transformations of
crystal’s potential, translations and local inversions, were devised.
Specifically, it was demonstrated that upon local inversion (inversion of
coordinates within the unit cell), the boundary charge maps to its negative,
up to an integral topological quantum number known as the interface invariant.
Likewise, upon the lattice translation by $x_{\varphi}$, the boundary charge
was shown to grow linearly with the shift variable $x_{\varphi}$ (with the
slope being the unit cell-averaged average charge density in the bulk
$\bar{\rho}=\frac{\nu}{L}$, with $\nu$ – filling factor and $L$-system’s
period), whilst performing discontinuous downward jumps by a unit of the
electron charge, as quantified by another topological quantum number – the
boundary invariant. These topological invariants were shown to be generated by
the spectral flow of the energies corresponding to the edge states inside the
energy gap that hosts the chemical potential, in complete analogy with the
integer quantum Hall effect [Hatsugai_1993, ]. As opposed to the edge states
in topological insulators, the quantization of these invariants does not rely
on the internal symmetries of the bulk Bloch’s Hamiltonian (such as particle-
hole or time-reversal symmetries) and is instead guaranteed by a number of
fundamental physical principles, such as charge conservation, Pauli principle,
and the nearsightedness principle of W. Kohn [Kohn_1978, , Kohn_1996, ,
Prodan_2005, , Prodan_2006, ] (to be discussed further on). Moreover, these
invariants are directly linked with the properties of an experimental
observable, a privilege shared by both the quantum Hall effect and the
topological defects (see Ref. [Trebin_1982, ] for a review), while not being
entirely clear in the domain of the topological insulators.

Further, in a different paper [Pletyukhov_2020_prr, ], rational quantization
of boundary and interface charges was discussed. Particularly, with the aid of
the aforementioned physical principles, a general framework for studying
quantized charges in one dimension was laid down, allowing us to quantify all
possible quantization patterns of the boundary charge in terms of the non-
symmorphic symmetries of the crystal. The charges on the interfaces between
pairs of insulators sharing their bulk properties were demonstrated to follow
a lattice version of the Goldstone-Wilczek formula [Goldstone_1981, ],
relating the interface charge to the sum of the boundary charges right and
left to their septum, modulo an unknown integer generated by the local
coupling between the two subsystems.

The key feature of the method developed in Ref. [Pletyukhov_2020_prr, ] is
this “modulo an unknown integer” paradigm, arising from the nearsightedness
principle of the electronic matter. As such, the nearsightedness principle
tells us that (see Ref. [Prodan_2005, , Prodan_2006, ]), in insulators,
localized perturbations by external fields may result in localized charge
redistributions only. To be more specific, the corrections beyond the
characteristic length scale $\xi_{g}=\frac{v_{F}}{E_{g}}$ (where $v_{F}$ and
$E_{g}$ are the Fermi velocity and the gap opening up at the Fermi level, see
Ref. [Weber_2021, ] for example) are exponentially suppressed. An even further
refinement of this statement would be that such perturbations may only remove
or add an additional number of bound states whose wave functions are localized
around the corresponding perturbations. One of the key purposes of the present
paper is to substantiate this claim mathematically, which turns out to be
possible in pretty general $d$-dimensional models.

To be more specific, this paper concerns the topological properties of the
electronic excess charges accumulated around point-like defects in
$d$-dimensional insulators. Although we purposefully specify the Hamiltonian
of the crystal under consideration to make our exposition more transparent,
the derivations presented in this manuscript are shown to be independent of
its choice. What indeed matters is that the spectrum of the clean system
consists of the energy bands occasionally separated by the energy gaps, that
is, there exists at least one bulk energy gap into which we can put the
chemical potential to promote the resulting statistical system into an
insulator.

Furthermore, neither we specify the internal structure of the impurity
vertices, nor do we assume any particular arrangement of them, making our
analysis applicable to a wide range of experimental setups. In particular,
quite conventionally, we may assume that a number of randomly located point-
like impurities exerting an ultra-screened electrostatic force on the system’s
electrons are scattered through the charge sampling region of a crystal under
consideration. A slightly less familiar situation is inspired by the work of
Nomura and Nagaosa [Nomura_2010, ] and may be formulated as follows. Assuming
that a crystal is further magnetic, we know that, in an insulating regime, its
ground state may accurately be described by a Heisenberg model that, by
itself, features topological defects. A familiar example of such a defect
would be a magnetic skyrmion or a hedgehog texture. Assuming that the total
spin of atoms comprising our crystal is large, these textures may be seen as
an arrangement of classical magnetic moments nailed down to the atomic
positions. Their interaction with the electron’s spin degree of freedom may
then be written as a sum of the Zeeman-like terms, each weighted with the
Dirac $\delta$-function centered at the position of the corresponding atom.

Quite generically, we show that the total electronic excess charge accumulated
around these defects is an integer-valued topological invariant, which we
express as a contour integral winding number of an appropriate functional of
bulk position space Green’s functions. Further analysis of this topological
quantum number reveals that upon an adiabatic modification of positions and/or
vertex functions of the localized scattering centers, the value of the
invariant may only be affected by the change in the occupancy of the
imperfection-localized bound states in the process of the spectral flow of
their eigenenergies inside the chemical potential-accommodating energy gap.
This observation allows for an immediate interpretation in terms of the
nearsightedness principle, as well as for a direct read-off of the central
memo of Ref. [Pletyukhov_2020_prr, ]: “localized perturbations in insulators
result in localized charge redistribution, leading to an addition/removal of
the corresponding perturbation-localized bound states to/from the occupied
spectral region”. We conclude our analysis by commenting on the relation
between the nearsightedness principle and the topological invariants
characterizing the boundary charge.

In what follows, we set the reduced Plank’s constant $\hbar$ and the electron
charge $e$ equal to unity $\hbar=e=1$.

##  II Adiabatic response of the excess charge to localized perturbations in
an insulating state

###  II.1 A translationally invariant model

In the following, we shall specifically refer to an electronic system governed
by the following Hamiltonian

$\displaystyle
H^{(0)}_{\mathbb{x}}=\frac{\mathbb{p}^{2}}{2m}+\frac{1}{2m}\sum_{j=1}^{d}\\{\tilde{A}_{j}(\mathbb{x}),p_{j}\\}+V(\mathbb{x}),$
(1)

with $V(\mathbb{x})$ and $\tilde{A}_{j}(\mathbb{x}),\ j=1,\ \dots,\ d$ being
the lattice periodic $N_{c}\times N_{c}$ Hermitian matrices. More
specifically,

$\displaystyle\begin{Bmatrix}V(\mathbb{x})\\\
\tilde{\mathbb{A}}(\mathbb{x})\end{Bmatrix}=\begin{Bmatrix}V(\mathbb{x}+\mathbb{R}_{\mathbb{m}})\\\
\tilde{\mathbb{A}}(\mathbb{x}+\mathbb{R}_{\mathbb{m}})\end{Bmatrix},\quad\forall\mathbb{m}\in\mathbb{Z}^{d},$
(2)

where $\mathbb{R}_{\mathbb{m}}=\sum_{j=1}^{d}m_{j}\mathbb{a}_{j}$ is a lattice
vector characterized by a $d$-dimensional vector of integers
$\mathbb{m}=\begin{pmatrix}m_{1}&\ldots&m_{d}\end{pmatrix}^{T}$, specifying
its components in the basis of primitive vectors $\\{\mathbb{a}_{j}\\}_{j}$
spanning the unit cell of a Bravais lattice. Furthermore, $\mathbb{p}$ and
$\mathbb{x}$ are vectorial momentum and position operators comprised of the
individual components $p_{j}=-i\frac{\partial}{\partial x_{j}}$ and $x_{j}$.

This model naturally generalizes the one recently studied in Ref. [Pias_2022,
] in connection with the universal properties of one-dimensional boundary
charge, to higher dimensions. We remark that other models of multi-dimensional
periodic structures [other_models, ] are expected to share the same physics,
as the effects we are about to describe are rather generic to an insulating
state.

Translationally invariant systems are characterized by their band structure,
comprised of the individual energy bands dispersing as
$\epsilon_{\alpha,\mathbb{k}},\ \alpha=1,\ 2,\ \dots$, as a function of the
vectorial quasimomentum variable $\mathbb{k}$, confined to the first Brillouin
zone of the reciprocal space. The eigenstates of the Hamiltonian to which
$\epsilon_{\alpha,\mathbb{k}}$ are the corresponding eigenvalues are known as
Bloch functions $\psi_{\alpha,\mathbb{k}}(\mathbb{x})$, and may be generically
expressed as

$\displaystyle\psi_{\alpha,\mathbb{k}}(\mathbb{x})=e^{i\mathbb{k}\cdot\mathbb{x}}u_{\alpha,\mathbb{k}}(\mathbb{x}),$
(3)

where $u_{\alpha,\mathbb{k}}(\mathbb{x})$ in the $N_{c}$-component object and
is lattice periodic in the same sense as vector and scalar potentials are
$u_{\alpha,\mathbb{k}}(\mathbb{x})=u_{\alpha,\mathbb{k}}(\mathbb{x}+\mathbb{R}_{\mathbb{m}}),\
\forall\mathbb{m}\in\mathbb{Z}^{d}$. The completeness and identity resolution
relations may be written as

$\displaystyle\frac{V_{\text{UC}}}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}d^{(d)}\mathbb{x}\psi_{\alpha,\mathbb{k}}^{\dagger}(\mathbb{x})\psi_{\alpha^{\prime},\mathbb{k}^{\prime}}(\mathbb{x})$
$\displaystyle=\delta_{\alpha,\alpha^{\prime}}\delta^{(d)}(\mathbb{k}-\mathbb{k}^{\prime}),$
(4)
$\displaystyle\frac{V_{\text{UC}}}{(2\pi)^{d}}\sum_{\alpha=1}^{\infty}\int_{\text{BZ}}d^{(d)}\mathbb{k}\psi_{\alpha,\mathbb{k}}(\mathbb{x})\psi_{\alpha,\mathbb{k}}^{\dagger}(\mathbb{x}^{\prime})$
$\displaystyle=1_{N_{c}}\delta^{(d)}(\mathbb{x}-\mathbb{x}^{\prime}),$ (5)

where $V_{\text{UC}}$ is the volume of the unit cell, defined via

$\displaystyle
V_{\text{UC}}=\int_{\text{UC}}d^{(d)}\mathbb{x}=\det\begin{pmatrix}\mathbb{a}_{1}\Big{|}&\ldots&\Big{|}\mathbb{a}_{d}\end{pmatrix}.$
(6)

When studying charge, it is more convenient to introduce the retarded single-
particle Green’s function, containing the information on both the eigenstates
and the energy spectrum. In thermodynamic equilibrium, the Laplace image of
the latter is defined as the resolvent of the single-particle Hamiltonian (1)

$\displaystyle[z-H^{(0)}_{\mathbb{x}}]G^{(0)}(\mathbb{x},\mathbb{x}^{\prime})=1_{N_{c}}\delta^{(d)}(\mathbb{x}-\mathbb{x}^{\prime}),$
(7)

where $z$ is the complex energy variable, defined in terms of the physical
frequency variable $\omega$ as $z=\omega+i\eta$, where $\eta\rightarrow
0^{+}$. Owing to the identity resolution relation (5) we can establish the
conventional Lehmann representation

$\displaystyle
G^{(0)}(\mathbb{x},\mathbb{x}^{\prime})=\frac{V_{\text{UC}}}{(2\pi)^{d}}\sum_{\alpha=1}^{\infty}\int_{\text{BZ}}d^{(d)}\mathbb{k}\frac{\psi_{\alpha,\mathbb{k}}(\mathbb{x})\psi_{\alpha,\mathbb{k}}^{\dagger}(\mathbb{x}^{\prime})}{z-\epsilon_{\alpha,\mathbb{k}}}.$
(8)

Further, using the completeness of the basis (4), in Appendix A, we establish
the following important fusion rule for the bare propagators

$\displaystyle\int_{\mathbb{R}^{d}}d^{(d)}\mathbb{x}^{\prime}G^{(0)}(\mathbb{x},\mathbb{x}^{\prime})G^{(0)}(\mathbb{x}^{\prime},\mathbb{x}^{\prime\prime})=-\frac{\partial}{\partial\omega}G^{(0)}(\mathbb{x},\mathbb{x}^{\prime\prime}).$
(9)

As it is shown in Appendix A, this relation holds pretty generally, without
any reference to the Hamiltonian (1).

###  II.2 Localized perturbations and Dyson’s equation

Now we perturb the translationally invariant (on the scale of the unit cell)
system by a finite number of point-like impurities

$\displaystyle\tilde{V}(\mathbb{x})=\sum_{n=1}^{N}\tilde{V}^{(n)}_{0}\delta^{(d)}(\mathbb{x}-\mathbb{x}_{n}),$
(10)

where $\tilde{V}^{(n)}_{0}$ are $N_{c}\times N_{c}$ matrices describing the
action of the $n^{\text{th}}$ impurity on the channel space. This action is
further assumed to be local as prescribed by Dirac delta-function
$\delta^{(d)}(\mathbb{x}-\mathbb{x}_{n})$ centered at the impurity position
$\mathbb{x}_{n}$.

Let us remark that the problem of a Dirac delta-function potential is well-
known to be ill-defined in spatial dimensions higher than $d=1$. In our
analysis, this is manifested in the ill-definiteness of the bulk position
space Green’s function at equal spatial arguments
$G^{(0)}(\mathbb{x},\mathbb{x})$ due to the divergence of the defining
integrals (8) in the ultraviolet. Such a divergence is not physical and has to
be circumvented by an appropriate regularization scheme. In particular, in the
metallic case $\tilde{\mathbb{A}}(\mathbb{x})=0,\ V(\mathbb{x})=0$, in $d>1$
the problem of the delta-potential has been extensively studied in both
physical [Park_1995, , Atkinson_1975, , Jackiw_1991, ] and mathematical
[Friedman_1972, ] literature and several meaningful regularization techniques
were proposed and shown to produce physically sensible results. Since the
presence of the energy gaps is of no importance in the deep ultraviolet
regime, the same methods may be applied in our case.

The Dyson’s equation for the full Green’s function of the system is given by

$\displaystyle G(\mathbb{x},\mathbb{x}^{\prime})=$ $\displaystyle
G^{(0)}(\mathbb{x},\mathbb{x}^{\prime})$
$\displaystyle+\sum_{n=1}^{N}G^{(0)}(\mathbb{x},\mathbb{x}_{n})\tilde{V}^{(n)}_{0}G(\mathbb{x}_{n},\mathbb{x}^{\prime}).$
(11)

First we want to consistently solve for the functions
$G(\mathbb{x}_{n},\mathbb{x}^{\prime}),\ n=1,\ \dots,\ N$. This problem is
brought to the solution of the following matrix equation

$\displaystyle\mathcal{M}(z)\mathcal{D}(\mathbb{x}^{\prime})=\mathcal{D}^{(0)}(\mathbb{x}^{\prime}),$
(12)

where $\mathcal{M}(z)$ is the $N_{c}\cdot N\times N_{c}\cdot N$ block matrix
defined by

$\displaystyle\mathcal{M}(z)=$ $\displaystyle 1_{N_{c}\cdot
N}-\mathcal{G}^{(0)}(z)\tilde{\mathcal{V}}_{0},$ (13)
$\displaystyle(\mathcal{G}^{(0)}(z))_{n,n^{\prime}}=$ $\displaystyle
G^{(0)}(\mathbb{x}_{n},\mathbb{x}_{n^{\prime}}),\
(\tilde{\mathcal{V}}_{0})_{n,n^{\prime}}=\delta_{n,n^{\prime}}\tilde{V}_{0}^{(n)}.$
(14)

Likewise, $\mathcal{D}(\mathbb{x}^{\prime})$ and
$\mathcal{D}^{(0)}(\mathbb{x}^{\prime})$ are the $N_{c}\cdot N\times N_{c}$
matrices comprised of the full $G(\mathbb{x}_{n},\mathbb{x}^{\prime})$ and
bare $G^{(0)}(\mathbb{x}_{n},\mathbb{x}^{\prime})$ propagators, respectively.
With these notations we obtain

$\displaystyle G(\mathbb{x},\mathbb{x}^{\prime})=$ $\displaystyle
G^{(0)}(\mathbb{x},\mathbb{x}^{\prime})+\mathcal{D}^{(0)\dagger}(\mathbb{x})\tilde{\mathcal{V}}_{0}\mathcal{D}(\mathbb{x}^{\prime})$
$\displaystyle=$ $\displaystyle
G^{(0)}(\mathbb{x},\mathbb{x}^{\prime})+\mathcal{D}^{(0)\dagger}(\mathbb{x})\tilde{\mathcal{V}}_{0}\mathcal{M}^{-1}(z)\mathcal{D}^{(0)}(\mathbb{x}^{\prime}),$
(15)

where in our definition the Hermitian conjugate does not affect the
$z$-variable, i.e.

$\displaystyle(G^{(0)}(\mathbb{x},\mathbb{x}^{\prime}))^{\dagger}=G^{(0)}(\mathbb{x}^{\prime},\mathbb{x}).$
(16)

###  II.3 Measuring the excess charge

We define the excess charge density operator in the following manner:

$\displaystyle\delta\widehat{\rho}(\mathbb{x})=\widehat{\rho}(\mathbb{x})-\bar{\rho},$
(17)

where

$\displaystyle\widehat{\rho}(\mathbb{x})=\widehat{\psi}^{\dagger}(\mathbb{x})\widehat{\psi}(\mathbb{x}),$
(18)

is the density operator, expressed in terms of the $N_{c}$-component fermionic
field operators $\widehat{\psi}(\mathbb{x})$ and
$\widehat{\psi}^{\dagger}(\mathbb{x})$. The field operators
$\widehat{\psi}(\mathbb{x})$ and $\widehat{\psi}^{\dagger}(\mathbb{x})$ are
further assumed to destroy/create excitations of the full Hamiltonian
including the effect of localized scattering centers in Eq. (10). The constant
contribution $\bar{\rho}$ describes the unit cell-averaged average charge
density in the bulk:

$\displaystyle\bar{\rho}=$
$\displaystyle\frac{1}{V_{\text{UC}}}\int_{V_{\text{UC}}}d^{(d)}\mathbb{x}\rho^{(0)}(\mathbb{x}),$
(19) $\displaystyle\rho^{(0)}(\mathbb{x})=$
$\displaystyle\braket{\widehat{\psi}^{(0)\dagger}(\mathbb{x})\widehat{\psi}^{(0)}(\mathbb{x})}$
$\displaystyle=$
$\displaystyle-\frac{1}{\pi}\text{Im}\int^{\mu}_{-\infty}d\omega\text{tr}\left\\{G^{(0)}(\mathbb{x},\mathbb{x})\right\\},$
(20)

where the field operators $\widehat{\psi}^{(0)}(\mathbb{x})$ and
$\widehat{\psi}^{(0)\dagger}(\mathbb{x})$ describe the excitations of the
translationally invariant system, $\mu$ denotes the chemical potential, and
$G^{(0)}(\mathbb{x},\ \mathbb{x}^{\prime})$ is the bare Green’s function
defined by Eqs. (7) and (8).

We measure the excess charge with the help of the classical device, described
by the envelope function $f(\mathbb{x})$ (see Refs. [Pletyukhov_2020, ,
Muller_2021, , Miles_2021, , Pias_2022, ] and Ref. [Kivelson_1982, ] for
similar definitions). To be more specific, we define the excess charge
operator as

$\displaystyle\delta\widehat{Q}=\int_{\mathbb{R}^{d}}d^{(d)}\mathbb{x}f(\mathbb{x})\delta\widehat{\rho}(\mathbb{x}).$
(21)

It is sensible to define the function $f(\mathbb{x})$ relative to a certain
point $\mathbb{x}_{p}$, to which the charge probe is applied, and further
assume that the charge is sampled equivalently in all directions
$f(\mathbb{x})=f(|\mathbb{x}-\mathbb{x}_{p}|)$. Additionally, we assume that
all of the charge $f(|\mathbb{x}-\mathbb{x}_{p}|)\approx 1$ is sampled in
sufficiently large vicinity of the sampling point $\mathbb{x}_{p}$, while the
envelope function smoothly decays to zero
$f(|\mathbb{x}-\mathbb{x}_{p}|)\rightarrow 0$ far away from $\mathbb{x}_{p}$.
For that matter, it is convenient to choose

$\displaystyle
f(|\mathbb{x}-\mathbb{x}_{p}|)=1-\Theta_{l_{p}}(|\mathbb{x}-\mathbb{x}_{p}|-L_{p}),$
(22)

where $\Theta_{l_{p}}(|\mathbb{x}-\mathbb{x}_{p}|-L_{p})$ is some
representation of the Heaviside function broadened by $l_{p}$. The length
scales characteristic of the charge probe are assumed to satisfy

$\displaystyle L_{p}\gg l_{p}\gg\xi_{g},$ (23)

where $\xi_{g}\simeq\frac{v_{F}}{E_{g}}$ is the charge localization length in
an insulator (also it is the charge correlation length, defining the
exponential decay length of the density-density correlation function, see Ref.
[Weber_2021, ]), roughly defined as the ratio between the Fermi velocity
$v_{F}$ and size of the energy gap at the Fermi level $E_{g}$.

###  II.4 Topological invariant characterizing the excess charge

Let us assume that $N$ impurities, as characterized by the potential (10), are
placed in a region of a crystal falling into the sampling district of the
envelope function $|\mathbb{x}|\lesssim L_{p}$. We define the total excess
charge as the zero temperature expectation value of the excess charge operator
in the grandcanonical equilibrium density matrix, so that

$\displaystyle\delta{Q}=$
$\displaystyle\braket{\delta\widehat{Q}}=\int_{\mathbb{R}^{d}}d^{(d)}\mathbb{x}f(\mathbb{x})(\rho(\mathbb{x})-\bar{\rho}),$
(24) $\displaystyle\rho(\mathbb{x})=$
$\displaystyle-\frac{1}{\pi}\text{Im}\int_{-\infty}^{\mu}d\omega\text{tr}\left\\{G(\mathbb{x},\mathbb{x})\right\\}.$
(25)

With the help of the representation (15), we obtain

$\displaystyle\delta{Q}=Q^{\prime}+Q_{P},$ (26)

where $Q^{\prime}$ contains the Friedel charge as well as the charge due to
the impurity-localized bound states

$\displaystyle Q^{\prime}$
$\displaystyle=\int_{\mathbb{R}^{d}}d^{(d)}\mathbb{x}f(\mathbb{x})\rho^{\prime}(\mathbb{x}),$
(27) $\displaystyle\rho^{\prime}(\mathbb{x})$
$\displaystyle=-\frac{1}{\pi}\text{Im}\int_{-\infty}^{\mu}d\omega\text{tr}\left\\{\mathcal{D}^{(0)\dagger}(\mathbb{x})\tilde{\mathcal{V}}_{0}\mathcal{M}^{-1}(z)\mathcal{D}^{(0)}(\mathbb{x})\right\\},$
(28)

while $Q_{P}$ is the so-called polarization charge given by

$\displaystyle
Q_{P}=\int_{\mathbb{R}^{d}}d^{(d)}\mathbb{x}f(\mathbb{x})(\rho^{(0)}(\mathbb{x})-\bar{\rho}),$
(29)

and, with the help of the properties of the envelope function, is shown to be
zero $Q_{P}=0$ in Appendix B. It hence follows that

$\displaystyle\delta{Q}=$ $\displaystyle
Q^{\prime}=-\frac{1}{\pi}\text{Im}\int_{\mathbb{R}^{d}}d^{(d)}\mathbb{x}f(\mathbb{x})$
$\displaystyle\times\int_{-\infty}^{\mu}d\omega\text{tr}\left\\{\mathcal{D}^{(0)\dagger}(\mathbb{x})\tilde{\mathcal{V}}_{0}\mathcal{M}^{-1}(z)\mathcal{D}^{(0)}(\mathbb{x})\right\\}.$
(30)

Due to the branch cuts and poles of the $T$-matrix
$T(\mathbb{x},\mathbb{x}^{\prime})=\sum_{n,n^{\prime}}[\tilde{\mathcal{V}}_{0}\mathcal{M}^{-1}(z)]_{n,n^{\prime}}\delta(\mathbb{x}-\mathbb{x}_{n})\delta(\mathbb{x}^{\prime}-\mathbb{x}_{n^{\prime}})$,
the integrand of the outer integral is exponentially suppressed $\sim
e^{-|\mathbb{x}|/\xi_{g}}$ at large $\mathbb{x}$, allowing us to set
$f(\mathbb{x})=1$.

Interchanging the order of the integrals, we consider

$\displaystyle\int_{\mathbb{R}^{d}}d^{(d)}\mathbb{x}\text{tr}\left\\{\mathcal{D}^{(0)\dagger}(\mathbb{x})\tilde{\mathcal{V}}_{0}\mathcal{M}^{-1}(z)\mathcal{D}^{(0)}(\mathbb{x})\right\\}$
$\displaystyle=-\sum_{n,n^{\prime}=1}^{N}\text{tr}\left\\{[\mathcal{M}^{-1}(z)]_{n,n^{\prime}}\frac{\partial}{\partial\omega}G^{(0)}(\mathbb{x}_{n^{\prime}},\mathbb{x}_{n})\tilde{V}_{0}^{(n)}\right\\}$
$\displaystyle=\sum_{n,n^{\prime}=1}^{N}\text{tr}\left\\{[\mathcal{M}^{-1}(z)]_{n,n^{\prime}}\frac{\partial}{\partial\omega}[\mathcal{M}(z)]_{n^{\prime},n}\right\\}$
$\displaystyle=\sum_{n=1}^{N}\text{tr}\left\\{\left[\mathcal{M}^{-1}(z)\frac{\partial}{\partial\omega}\mathcal{M}(z)\right]_{n,n}\right\\}$
$\displaystyle=\frac{\partial}{\partial\omega}\text{tr}\left\\{\log\mathcal{M}(z)\right\\}=\frac{\partial}{\partial\omega}\log\det\left\\{\mathcal{M}(z)\right\\},$
(31)

where, in the last line, trace and determinant of the full $N_{c}\cdot N\times
N_{c}\cdot N$ block matrix $\mathcal{M}(z)$ are understood. Using the result
in Eq. (31), we arrive at the following compact formula for the total excess
charge

$\displaystyle\delta{Q}=$
$\displaystyle-\frac{1}{\pi}\text{Im}\int_{-\infty}^{\mu}d\omega\frac{\partial}{\partial\omega}\log\det\left\\{\mathcal{M}(z)\right\\}.$
(32)

To see why the integral in Eq. (32) may take on integral values only, in
Appendix C we find an alternative contour integral representation

$\displaystyle\delta{Q}=$ $\displaystyle-\oint_{C}\frac{dz}{2\pi
i}\frac{\partial}{\partial z}\log\det\left\\{\mathcal{M}(z)\right\\},$ (33)

where $C$ is an arbitrary non-self-intersecting curve that crosses the real
axis at two points only, below the lowest eigenvalue of the full Hamiltonian
and at the chemical potential $\mu$, and the direction of $C$ is assumed to be
clockwise.

In the representation (33), the excess charge $\delta{Q}$ is necessarily an
integer as it is expressed as a contour integral winding number and the
chemical potential is by definition inside one of the energy gaps (we focus on
the insulating systems solely). In other words, the integral in Eq. (33)
measures the degree of the mapping $S^{1}\rightarrow S^{1}$ and is thus a
member of the only non-trivial homotopy group of the unit circle
$\pi_{1}(S^{1})=\mathbb{Z}$.

In particular, the integral in Eq. (33), is a sum of two distinct
contributions: the contribution of the branch cuts corresponding to the
extended or scattering states, and the contribution of poles corresponding to
the imperfection-localized bound states.

The bands in multidimensional ($d>1$) and/or multichannel ($N_{c}>1$) systems
are typically composite, i.e. overlapping with one another along the frequency
axis. For that matter, it is convenient to choose the branch cuts to connect
the bottom of the lowest sub-band with the top of the highest one, within
every patch of the energy bands surrounded by a pair of energy gaps.

The bound state poles, determined as a solution of
$\det\left\\{\mathcal{M}(z)\right\\}\big{|}_{z\in\mathbb{R}}=0$, are located
on the complement of the bare Hamiltonian’s spectrum, i.e. inside the energy
gaps and, in some cases (e.g. an attractive scalar impurity), below the bottom
of the lowest energy band of the unperturbed Hamiltonian.

##  III Relation with the nearsightedness principle

Figure 1: A schematic illustration of how the spectral flow of the energies of
the imperfection-circumscribing bound states sitting inside the gap that
accommodates the Fermi level affects the total excess charge. The spectrum of
the system is visualized through the local spectral density as looked down on
the complex frequency plane. The occupied part of the spectrum is demonstrated
in blue, while the yellow color marks its complement (the states of the system
that are unoccupied). Panel a) shows a rectangular contour $C$ encircling the
occupied spectral region. Panels b) and c) show the zoomed-in vicinity of the
chemical potential before and after the perturbation. As is demonstrated in
panel c), the spectral flow results in the removal of a single bound state,
carrying away a unity of the electron charge from the system (an inverse
process is of course also possible).

###  III.1 Discussion

Now we would like to discuss the topological invariant (33) in greater detail.
In what follows, we specify the contour $C$ as a rectangle of length $\mu-B$
in the real direction and width $2\eta$ in the imaginary one. Here $B$ is by
definition an energy lying below the lowest eigenvalue of the total
Hamiltonian $H_{\mathbb{x}}=H_{\mathbb{x}}^{(0)}+\tilde{V}(\mathbb{x})$ (i.e.
$B\in(-\infty,\min\\{\text{spec}\\{H_{\mathbb{x}}\\}\\})$), and $\eta$ is not
necessarily an infinitesimal positive but is rather a finite positive number
(which is allowed as the integral is invariant under such contour deformations
(see Appendix C)). Furthermore, we assume that the chemical potential is
located above the $\nu^{\text{th}}$ bulk energy band.

Let us now consider making an adiabatic perturbation to the system that is
comprised of the change in the positions $\\{\mathbb{x}_{n}\\}_{n}$ and/or
vertex functions $\\{\tilde{V}^{(n)}_{0}\\}_{n}$ of the impurities. As the
span of the extended states’ energy bands is unaffected by such adiabatic
perturbations, the branch cut contribution to the winding number remains
invariant (up to the cases when the bound state merges with the band, as
discussed below). This remark is essentially true as such deformations of the
parameter space do not change the analytical structure of
$G^{(0)}(\mathbb{x}_{n},\ \mathbb{x}_{n^{\prime}})$, through the functionals
of which alone our topological invariant is expressed. We hence conclude that
such changes may only unleash themselves in the spectral flow of the bound
state energies.

As was anticipated in Section II.4, the bound state energies are energy-wise
located inside the energy gaps of the bulk system. This assertion also regards
the energy gap below the bottom of the lowest band $\omega\in(-\infty,\
\min_{\mathbb{k}}\epsilon_{1,\mathbb{k}}]$, which can accommodate the bound
states in the case of attractive impurities, for example.

The energies of the bound states $\epsilon_{\text{bs}}$ inside the energy gaps
$[\max_{\mathbb{k}}\epsilon_{\alpha,\mathbb{k}},\
\min_{\mathbb{k}}\epsilon_{\alpha+1,\mathbb{k}}]$ surrounded by a pair of
bands $\alpha,\ \alpha+1,\ (\alpha=1,\ ...,\ \nu-1)$, are solely characterized
by their location within the gap. The same holds true for the infinite gap
below the bottom of the lowest bulk energy band, with $\epsilon_{\text{bs}}$
now being energy-wise located in $(-\infty,\
\min_{\mathbb{k}}\epsilon_{1,\mathbb{k}}]$. This implies that the spectral
flow of these energies is constituted in the motion of $\epsilon_{\text{bs}}$
in between the top of $\epsilon_{\alpha,\mathbb{k}}$ and the bottom of
$\epsilon_{\alpha+1,\mathbb{k}}$, or between the negative infinity and
$\min_{\mathbb{k}}\epsilon_{1,\mathbb{k}}$ shall some states be also found in
there. When merging with one of the energy bands (either
$\epsilon_{\alpha,\mathbb{k}}$ or $\epsilon_{\alpha+1,\mathbb{k}}$, and
$\epsilon_{1,\mathbb{k}}$ solely when considering the gap preceding the entire
band structure), the value of the contour integral winding number (33)
relating to that band gets modified by unity footnote_2 . It follows that the
motion of the bound state poles, inside such energy gaps below the one hosting
the chemical potential, has absolutely no effect on the topological invariant
(33) (one may see this result as a form of charge conservation), as $B$, by
definition, resides below the lowest pole (effectively meaning that none of
the states are allowed to escape the occupied spectral region from below).

The flow of the energies of the impurity-localized bound states residing
inside the gap separating the conduction and the valence bands apart (the gap
where the chemical potential is located), on the other hand, affects the
winding number in Eq. (33). When a bound state crosses the chemical potential
from above or below, the number of poles encompassed by the integration
contour increases or decreases correspondingly. That means that the unit of
the electron charge gets either pumped in or out of the system, modifying the
topological invariant by $\pm 1$. This discussion is summarized in Fig. 1.

The elaboration above allows us to draw the following physical conclusion:

Localized adiabatic perturbations in insulators, may only result in the
localized charge redistributions, owing to the change in the occupancy of the
perturbation-localized bound states at the Fermi level.

This intuitive result is nothing but a direct consequence of the universal
nearsightedness principle of W. Kohn [Kohn_1996, , Prodan_2005, , Prodan_2006,
] stating that, at fixed chemical potential, the electronic charge density
depends on the external field (in our case being an assembly of localized
scattering centers) only at nearby points.

Another conclusion drawn by E. Prodan and W. Kohn in Ref. [Prodan_2005, ] (see
also Ref. [Prodan_2006, ] for the fine details in $d=1$) is that the adiabatic
perturbations to the external potential, no matter how strong, have a
negligible effect on the local charge density beyond a certain characteristic
length scale, which, in the insulating regime, is naturally provided by the
charge correlation length $\xi_{g}$. From the viewpoint of our topological
invariant (33), this means that in the case of well-separated impurities
$|\mathbb{x}_{n}-\mathbb{x}_{n^{\prime}}|/\xi_{g}\gg 1$, the topological
invariant is expected to approach a sum of the individual single-impurity
invariants, as distant impurities are not supposed to be able to “talk” with
one another on such scales. Indeed, in an insulating state, it is well-known,
that the two-point correlation functions $G^{(0)}(\mathbb{x}_{n},\
\mathbb{x}_{n^{\prime}})$ decay exponentially at large distances $\sim
e^{-|\mathbb{R}_{\mathbb{m}_{n}}-\mathbb{R}_{\mathbb{m}_{n^{\prime}}}|/\xi_{g}}$
(where $\mathbb{m}_{n}$ labels the unit cell accommodating the $n^{\text{th}}$
scattering center), meaning that we can approximate

$\displaystyle(\mathcal{G}^{(0)}(z))_{n,n^{\prime}}\simeq$
$\displaystyle\delta_{n,n^{\prime}}G^{(0)}(\mathbb{x}_{n},\mathbb{x}_{n}),$
(34)

implying that

$\displaystyle\mathcal{M}(z)\simeq$
$\displaystyle\bigoplus_{n=1}^{N}\left(1_{N_{c}}-G^{(0)}(\mathbb{x}_{n},\mathbb{x}_{n})V_{0}^{(n)}\right),$
(35)

and

$\displaystyle\delta{Q}\simeq$
$\displaystyle-\sum_{n=1}^{N}\oint_{C}\frac{dz}{2\pi
i}\frac{\partial}{\partial
z}\log\det\left\\{1_{N_{c}}-G^{(0)}(\mathbb{x}_{n},\mathbb{x}_{n})V_{0}^{(n)}\right\\}.$
(36)

This result may be seen as a form of the conventional Born approximation of
the linear transport theory, whereby, to the lowest order in the impurity
density, one considers impurities as independent.

###  III.2 An illustration: A pair of magnetic impurities in an illuminated
quantum wire

Figure 2: Panel a): A schematic illustration of a ballistic quantum wire
featuring Rashba-style spin-orbit coupling (defined by a spin-orbit vector
$\mathbb{k}_{R}$) and submersed into a spatially-periodic arrangement of
electric $\mathbb{E}$ and magnetic $\mathbb{B}$ fields of wavelength
$\lambda$. The two impurity atoms, separated by distance $R$ and carrying an
effective magnetic moment of $\mathbb{m}_{\text{eff}}^{(j)},\ j=1,2$, are
schematically shown by atomic symbols pierced with the magnetic moment-
symbolizing arrows. Panel b): The bulk energy spectrum of the two-impurity
problem. The energy bands are shown in dark blue, the chemical potential
located inside the second spectral gap (above the fourth energy band) is
depicted in orange, and the relevant spectral region is highlighted in light
blue.

To illustrate some of the points highlighted in the above discussion, we here
consider a simple model of a spin-orbit-interacting ballistic quantum wire,
submersed into the background of the spatially oscillating electromagnetic
field. The bulk Hamiltonian assumes the form of the Pauli Hamiltonian with an
extra Rashba-like term:

$\displaystyle H_{x}^{(0)}=$
$\displaystyle\frac{\left(p+\frac{e}{c}A_{x}(x)\right)^{2}}{2m}+\frac{\mathbb{k}_{R}\cdot{\bm{\sigma}}}{m}\left(p+\frac{e}{c}A_{x}(x)\right)$
$\displaystyle+\frac{\mu_{B}g_{e}}{2}{\bm{\sigma}}\cdot\mathbb{B}(x).$ (37)

Above, $\mathbb{k}_{R}=(k_{R,x},\ k_{R,y},\ k_{R,z})$ is the Rashba spin-orbit
vector, ${\bm{\sigma}}=(\sigma_{x},\ \sigma_{y},\ \sigma_{z})$ is the vector
of the Pauli spin matrices, $\mu_{B}=\frac{e}{2mc}$ is the Bohr magneton, $c$
is the speed of light in vacuum, $g_{e}$ is the electron’s Landé g-factor, and

$\displaystyle\mathbb{B}(x)=$
$\displaystyle\nabla\times\mathbb{A}(\mathbb{x})\Big{|}_{\mathbb{x}=x\hat{\mathbb{e}}_{x}},\quad
A_{x}(x)=\hat{\mathbb{e}}_{x}\cdot\mathbb{A}(\mathbb{x})\Big{|}_{\mathbb{x}=x\hat{\mathbb{e}}_{x}},$
(38)

with $\hat{\mathbb{e}}_{x}$ being the ort in the $x$-direction, and
$\mathbb{A}(\mathbb{x})$ being the electromagnetic vector potential of the
monochromatic plane-wave form

$\displaystyle\mathbb{A}(\mathbb{x})=\mathbb{A}_{0}\cos(\mathbb{q}\cdot\mathbb{x}+\varphi),$
(39)

in the Coulomb gauge

$\displaystyle\nabla\cdot\mathbb{A}(\mathbb{x})=0\iff\mathbb{q}\cdot\mathbb{A}_{0}=0.$
(40)

The wave vector of the background electromagnetic field defines the fictitious
lattice spacing

$\displaystyle L=\frac{2\pi}{\hat{\mathbb{e}}_{x}\cdot\mathbb{q}},$ (41)

where we have excluded the uninteresting case of the orthogonally propagating
wave $\hat{\mathbb{e}}_{x}\cdot\mathbb{q}=0$.

We note that the Hamiltonian in Eq. (37) falls into the class of systems
defined by the Hamiltonian (1), with $d=1$ and

$\displaystyle V(x)=$
$\displaystyle\frac{\mu_{B}g_{e}}{2}{\bm{\sigma}}\cdot\mathbb{B}(x)+\frac{e^{2}A_{x}^{2}(x)}{2mc^{2}}+\frac{e\mathbb{k}_{R}\cdot{\bm{\sigma}}A_{x}(x)}{mc},$
(42) $\displaystyle\tilde{A}_{x}(x)=$
$\displaystyle\frac{e}{c}A_{x}(x)+\mathbb{k}_{R}\cdot{\bm{\sigma}}.$ (43)

As this demonstration is assumed to be interpretative, it suffices to consider
the case of a pair of impurities, which we assume to be separated by distance
$R$:

$\displaystyle\tilde{V}(x)=\tilde{V}^{(1)}_{0}\delta(x)+\tilde{V}^{(2)}_{0}\delta(x-R).$
(44)

Note that we can place the first impurity at $x=0$ without loss of generality,
as its other positions inside the wire may be achieved by appropriate tuning
of the modulation’s phase $\varphi$. Furthermore, we assume the impurities to
exert both the electrostatic and the exchange “force” on the wire’s electrons,
which we encode in the following form of the impurities’ vertex functions

$\displaystyle\tilde{V}^{(j)}_{0}=U_{j}\sigma_{0}+\frac{\mu_{B}g_{e}}{2}{\bm{\sigma}}\cdot\mathbb{B}_{\text{eff}}^{(j)},$
(45)

where $\mathbb{B}_{\text{eff}}^{(j)}$ is the effective (also appropriately
screened to have a short-ranged effect only) magnetic field, produced by the
effective magnetic moment of the impurity atom
$\mathbb{m}_{\text{eff}}^{(j)}=\frac{q_{j}g_{j}}{2M_{j}c}\mathbb{S}^{(j)}$,
with $q_{j}$, $g_{j}$, and $M_{j}$ being the charge, g-factor, and mass of the
$j^{\text{th}}$ impurity. Furthermore, $U_{j}$ denotes the strength of the
electrostatic potential, defining the corresponding force exerted by the
impurity on the electrons. Not going into much of the microscopic details, in
the following, we treat $U_{j}$ and $\mathbb{B}_{\text{eff}}^{(j)}$ as some
constant parameters. The resulting setup is schematically illustrated in panel
a) of Fig. 2.

Further, to illustrate our point, we assume that the associated impurity
parameters $\\{R,\ \\{U_{j}\\}_{j},\
\\{\mathbb{B}_{\text{eff}}^{(j)}\\}_{j}\\}$ evolve with a fictitious
“adiabatic time” $\tau\in[0,T]$, in such a manner that their temporal
derivatives remain much smaller than the Fermi energy $\epsilon_{F}$ times
their value, for all $\tau\in[0,T]$.

Figure 3: The figure demonstrates the adiabatic flows of both the bound state
energy spectrum and the excess charge invariant in the toy model proposed in
Section III.2. Specifically, the spectral flow of the impurity-localized bound
state energies is shown in the upper row, while the second row is dedicated to
the invariant itself. As is explained in Appendix D, the position of the
second impurity is parametrized as $R(\tau)=(n_{R}-1)L+\bar{R}(\tau)$, where
$n_{R}$ denotes the number of the unit cell hosting the second imperfection,
and $\bar{R}(\tau)\in[0,L]$ describes its location within the unit cell. The
five distinct columns in the above figure correspond to five choices of
$n_{R}=1,\ 2,\ 3,\ 4,\ 5$. In all of the panels, red dashed lines correspond
to the actual solution, while black solid lines relate to the case of two
independent impurities (see the approximate formula (36)). As the separation
between the impurities becomes of the order of the charge localization length
$\xi_{g}=O(L)$ (see Appendix D), both adiabatic flows approach the limit of
two independent impurities.

The particular form of the pumping protocol used to produce the numerical data
and the concrete numerical values of the free model parameters are provided in
Appendix D. The resulting bulk energy spectrum is demonstrated in panel b) of
Fig. 2.

The numerical data for the excess charge-characterizing topological invariant,
as well as the spectral flow of bound state energies inside the chemical
potential-accommodating spectral gap, is shown in Fig. 3. In particular, using
the parametrization $R(\tau)=(n_{R}-1)L+\bar{R}(\tau)$, suggested in the
Appendix D, we present the data for five different values of
$n_{R}\in\\{1,\dots,5\\}$, as shown in five different columns of the
corresponding figure, with upper and lower rows corresponding to the spectral
flow and the topological invariant, respectively. Solid black and dashed
burgundy lines mark the cases of independent and “interacting” impurities,
correspondingly. By independent impurities, we here understand that the
separation between them is effectively infinite, so that the off-diagonal
blocks of the $\mathcal{M}(z)$ matrix (see Eq. (14) for the definition) may be
completely ignored. This means that the bound state spectrum of the
independent impurities is provided by the solution of
$\det\left(1_{2}-G^{(0)}(0,0)\tilde{V}^{(1)}_{0}\right)\det\left(1_{2}-G^{(0)}(R,R)\tilde{V}^{(2)}_{0}\right)\big{|}_{z\in\mathbb{R}}=0$,
while the topological invariant is given by the approximation (36). By the
“interacting” impurities, on the other hand, we understand that the exact
relations were used to produce the numerical data. The numerical technique for
evaluation of bulk position space Green’s functions, as well as the
topological indices of the form (33), is outlined in Ref. [Pias_2022, ]. In
our calculations, the values of the contour parameters were chosen as
$\eta=1,\ B=-30$ (such a choice of $B$ is motivated by the presence of the
bound states below the lowest band $\omega\in(-\infty,\
\min_{k}\epsilon_{k,1}]$ in our model (37)).

The central purpose of our demonstration is to show that upon the increase in
the impurity’s separation beyond the charge localization length $\xi_{g}=O(L)$
(see Appendix D), both the topological invariant and the bound state spectrum
approach that of a pair of independent impurities. This effect is a direct
consequence of the nearsightedness principle, telling us that a localized
cause leads to a localized effect. Furthermore, as one may anticipate, the
discontinuous jumps of the excess charge invariant occur precisely at the
points where bound states enter/leave the occupied part of the energy
spectrum, as is explained in Section III.1. Another interesting observation is
the non-zero value of the topological invariant at the beginning of the
adiabatic evolution in $\tau$, where the strengths of the electrostatic
repulsion are the smallest $0<U_{j}\ll 1$ (see Appendix D). This feature is a
consequence of the presence of impurity-localized bound states below the
bottom of the lowest energy band. Such an effect is well-known in the case of
attractive scalar impurities, whereas here, it is generated by the non-Abelian
structure of the model, and, to the best of our knowledge, was not reported
previously in the literature.

###  III.3 Topological invariants characterizing the boundary charge in
unidimensional crystals

In this section, we would like to comment on the topological invariants
characterizing boundary charges in unidimensional crystals, extensively
discussed in Refs. [Pletyukhov_2020, , Muller_2021, , Miles_2021, , Pias_2022,
]. In particular, let us consider a $d=1$ semi-infinite system described by
the Hamiltonian (1), with the boundary placed at $x=x_{b}$. An appropriate
restriction of $x$ defines the respective right and left subsystems:

$\displaystyle x$ $\displaystyle\in[x_{b},\infty),\quad\text{right sub-
system},$ (46) $\displaystyle x$
$\displaystyle\in(-\infty,x_{b}],\quad\text{left sub-system}.$ (47)

In our definition, the primitive unit cell is defined as the one starting at
the boundary of the right semi-infinite system $\text{UC}=[x_{b},\ x_{b}+L]$,
with $L$ being the lattice period. In this definition, the left half-system is
always obtained from the right one by a local inversion operation, which acts
by the inversion of local coordinates within each unit cell.

Now we define the boundary charge operators corresponding to right and left
semi-infinite systems as the envelope-weighted integrals of the expectation
values of the appropriate excess charge density operators:

$\displaystyle Q_{B}^{(R)}=$
$\displaystyle\int_{x_{b}}^{\infty}dxf(x)\braket{\delta\widehat{\rho}_{R}(x)},$
(48) $\displaystyle Q_{B}^{(L)}=$
$\displaystyle\int_{-\infty}^{x_{b}}dxf(x)\braket{\delta\widehat{\rho}_{L}(x)},$
(49)

where, in analogy with Eq. (17),
$\delta\widehat{\rho}_{S}(x)=\widehat{\rho}_{S}(x)-\bar{\rho}$, and
$\widehat{\rho}_{S}(x)$ is the density operator referring to the system $S=R,\
L$. Furthermore, the envelope function $f(x)$ is chosen in accordance with Eq.
(22), with $x_{p}=x_{b}$, and the range of $x$ being restricted according to
Eqs. (46) and (47).

Let us now consider measuring the total excess charge $\delta{Q}$ accumulated
around $x=x_{b}$ in a translationally invariant system $x\in(-\infty,\
\infty)$. By the polarization charge neutrality condition $Q_{P}=0$,
demonstrated in Appendix B, the total excess charge also vanishes
$\delta{Q}=0$. On the other hand, we may consider a translationally invariant
system as a sum of right and left semi-infinite systems with a coupling
corresponding to the bulk Hamiltonian switched in between them. This coupling
manifests itself as a local perturbation and, by the nearsightedness principle
of Kohn, is capable of affecting the total charge locally by at most
introducing or removing a number of additional bound states, resulting in an
integer contribution $Q_{I}$. In this connection, we conclude that
$\delta{Q}=Q_{B}^{(R)}+Q_{B}^{(L)}-Q_{I}=0$, where $Q_{I}$ is known as the
interface invariant. One of the central results of Ref. [Pias_2022, ], was to
demonstrate that

$\displaystyle Q_{I}=$ $\displaystyle
Q_{B}^{(R)}+Q_{B}^{(L)}=-\oint_{C}\frac{dz}{2\pi i}\frac{\partial}{\partial
z}\log\det\left\\{G^{(0)}(x_{b},x_{b})\right\\}.$ (50)

That is, the interface invariant, characterizing the boundary charge upon
local inversions, is a topological quantum number given by the winding of the
determinant of bulk position space Green’s function evaluated at the location
of the boundary.

Now let us proceed with the transformations of the boundary charge under
translations. First, we consider the right boundary charge of the so-called
reference system, starting at $x_{b}=0$:

$\displaystyle Q_{B}^{(R)}(0)=\int_{0}^{\infty}dxf(x)(\rho(x)-\bar{\rho}),$
(51)

and we would like to analyze the changes in this quantity upon the translation
of the boundary by $x_{\varphi}\in[0,\ L]$. Instead of shifting the boundary,
we consider adding the following potential
$\hat{V}(x)=\hat{V}_{0}\Theta(x)\Theta(x_{\varphi}-x),\
\hat{V}_{0}\rightarrow\infty$. By the Pauli principle, the charge density
becomes zero for $x\in[0,x_{\varphi}]$ as these states sit at infinite energy
above the chemical potential $\mu$. From the definition of the boundary
charge, we are left with the following contribution:

$\displaystyle\delta Q_{B}^{(R)}(x_{\varphi})$
$\displaystyle=\int_{0}^{x_{\varphi}}dxf(x)(0-\bar{\rho})\mod 1$
$\displaystyle=-\bar{\rho}x_{\varphi}\mod 1,$ (52)

where $\mod 1$ contribution again comes from the nearsightedness principle.
This analysis allows us to conclude that:

$\displaystyle
Q_{B}^{(R)}(x_{\varphi})-Q_{B}^{(R)}(0)=\bar{\rho}x_{\varphi}+I(x_{\varphi}),$
(53)

where $I(x_{\varphi})$ is known as the boundary invariant. Another important
result of Ref. [Pias_2022, ] was to show that

$\displaystyle I(x_{\varphi})=-\oint_{C}\frac{dz}{2\pi
i}\frac{\partial}{\partial z}\ln\det\mathcal{U}(x_{\varphi}),$ (54)

where $\mathcal{U}(x_{\varphi})$ is defined via the path-ordered exponential

$\displaystyle\mathcal{U}(x_{\varphi})=$
$\displaystyle\text{P}\\!\exp\left\\{\int_{0}^{x}dx^{\prime}\mathcal{L}(x^{\prime})\right\\},$
(55) $\displaystyle\mathcal{L}(x)=$
$\displaystyle[G^{(0)}(x,x)]^{-1}G_{2}^{(0)}(x,x^{+})-iA(x),$ (56)

and $G_{2}^{(0)}(x,x^{\prime})=\partial_{x^{\prime}}G^{(0)}(x,x^{\prime})$. In
other words, the boundary invariant is also a topological quantum number
expressed as a winding of the appropriate functional of bulk position space
Green’s functions.

In this way, we see that the quantization of the topological invariants
characterizing the boundary charge in one-dimensional insulators is a direct
consequence of the nearsightedness principle. As this intuitive physical
principle holds beyond the single spatial dimension, one expects the excess
charges accumulated on inhomogeneities of various spatial co-dimensions in
$d$-dimensional crystals to possess similar topological characterization
schemes. Indeed, linear scaling of the boundary charge, along with its
discontinuous jumps by a unit of the electron charge at the bound state
escape/entrance spectral points, was recently demonstrated in a two-
dimensional system [Hou_2022, ].

##  IV Conclusions and Outlook

In this paper, the quantization of the excess charges on localized scattering
centers in $d$-dimensional insulators was discussed. Our analysis reveals that
an assembly of such imperfections accumulates an integral excess charge, given
by a winding number expression. We find that an adiabatic perturbation (no
matter how strong) comprised of either relocation of the impurities or a
modification of their vertex functions (or both at the same time) results in
the change of the total charge by an integer, determined by the saldo of the
imperfection-localized bound states that entered or escaped the occupied
spectral region, inside the chemical potential-hosting bulk spectral gap. The
quantization of this topological invariant was shown to be a direct
consequence of the nearsightedness principle of the electronic matter,
limiting the range of the effect of a localized cause. Additionally, this
local behavior of the electronic matter in the insulating state was shown to
be responsible for the quantization of the topological invariants
characterizing the unidimensional boundary charge studied in [Pletyukhov_2020,
, Muller_2021, , Miles_2021, , Pias_2022, ]. Furthermore, our study confirms
the central paradigm of Ref. [Pletyukhov_2020_prr, ], namely that localized
perturbations in insulators specifically lead to the change in occupancy of
the corresponding perturbation-localized bound states, modifying the total
charge, defined as the macroscopic average on the scales significantly
exceeding both the unit cell size $L$ and the charge correlation length
$\xi_{g}$, by at most an integer.

As is now obvious, the present paper is of conceptual value only as the
evaluation of the suggested topological invariant (33) for a specific multi-
impurity ($N\gg 1$) system poses a challenge on its own. In particular, this
concerns questions regarding the regularization schemes for the higher-
dimensional equal-argument Green’s functions, as well as the basic questions
regarding the numerical feasibility of the problem. Furthermore, it would be
of future interest to study the expansion of the topological invariant in the
interaction between the individual impurities, as generated by the off-
diagonal blocks of the $\mathcal{M}(z)$ matrix, and analyze its ties with the
conventional Born series for the impurity-dressed $T$-matrix. As it is
suggested in the present study, in the insulating state, the impurity density
$\rho_{I}$ has to be always contrasted with the inverse charge localization
length $\xi_{g}$, in such a manner that the condition
$1\gg\rho_{I}\xi_{g}^{d}$ implies the validity of the Born approximation,
treating impurities as independent.

##  V Acknowledgments

The author gratefully acknowledges the durable exchange of ideas with M.
Pletyukhov and H. Schoeller. Further, the author generously thanks S. Miles
and M. Pletyukhov for their valuable comments.

Most of the present work was done at the Institut für Theorie der
Statistischen Physik of RWTH Aachen and was financially supported by the
Deutsche Forschungsgemeinschaft via RTG 1995.

##  Appendix A Contraction of two Green’s functions

Quite generically we may represent

$\displaystyle G=\sumint_{s}\frac{\ket{s}\bra{s}}{z-\epsilon_{s}},$ (57)

where the meta-index $s$ labels the eigenstates $\ket{s}$ and eigenenergies
$\epsilon_{s}$ of the Hamiltonian. Considering the product of the Green’s
function with itself

$\displaystyle GG$
$\displaystyle=\sumint_{s}\sumint_{s^{\prime}}\frac{\ket{s}\bra{s^{\prime}}}{(z-\epsilon_{s})(z-\epsilon_{s^{\prime}})}\underbrace{\braket{s}{s^{\prime}}}_{\delta(s,s^{\prime})}=\sumint_{s}\frac{\ket{s}\bra{s}}{(z-\epsilon_{s})^{2}}$
$\displaystyle=-\frac{\partial}{\partial\omega}\sumint_{s}\frac{\ket{s}\bra{s}}{z-\epsilon_{s}}=-\frac{\partial}{\partial\omega}G.$
(58)

Taking the position space matrix elements

$\displaystyle\braket{\mathbb{x}}{GG}{\mathbb{x}^{\prime\prime}}$
$\displaystyle=-\frac{\partial}{\partial\omega}G(\mathbb{x},\mathbb{x}^{\prime\prime}),$

and inserting

$\displaystyle
1=\int_{\mathbb{R}^{d}}d^{(d)}\mathbb{x}\ket{\mathbb{x}}\bra{\mathbb{x}},$
(59)

we obtain the desired identity

$\displaystyle\int_{\mathbb{R}^{d}}d^{(d)}\mathbb{x}^{\prime}G(\mathbb{x},\mathbb{x}^{\prime})G(\mathbb{x}^{\prime},\mathbb{x}^{\prime\prime})$
$\displaystyle=-\frac{\partial}{\partial\omega}G(\mathbb{x},\mathbb{x}^{\prime\prime}).$
(60)

##  Appendix B Polarization charge

We consider

$\displaystyle\int_{\mathbb{R}^{d}}d^{(d)}\mathbb{x}f(\mathbb{x})(\rho^{(0)}(\mathbb{x})-\bar{\rho})$
$\displaystyle=\sum_{\mathbb{m}}\int_{\text{UC}}d^{(d)}\bar{\mathbb{x}}f(\bar{\mathbb{x}}+\mathbb{R}_{\mathbb{m}})(\rho^{(0)}(\bar{\mathbb{x}})-\bar{\rho}).$
(61)

Above we parametrized the position space variable $\mathbb{x}$ as
$\mathbb{x}=\mathbb{R}_{\mathbb{m}}+\bar{\mathbb{x}}$, for some vector of
integers $\mathbb{m}$, and $\bar{\mathbb{x}}$ is the local coordinate within
the unit cell $\bar{\mathbb{x}}\in\text{UC}$. Furthermore, we used the
periodicity property of $\rho^{(0)}(\mathbb{x})$, implied by the periodicity
of the equal-argument Green’s function

$\displaystyle G^{(0)}(\mathbb{x},\
\mathbb{x})=G^{(0)}(\mathbb{x}+\mathbb{R}_{\mathbb{m}},\
\mathbb{x}+\mathbb{R}_{\mathbb{m}}),\quad\forall\mathbb{m}\in\mathbb{Z}^{d}.$
(62)

The envelope function varies significantly only in the crossover region
$|\mathbb{R}_{\mathbb{m}}|=O(L_{p})$, allowing us to approximate

$\displaystyle f(\bar{\mathbb{x}}+\mathbb{R}_{\mathbb{m}})\approx
f(\mathbb{R}_{\mathbb{m}})+\bar{\mathbb{x}}\cdot\nabla{f}(\mathbb{R}_{\mathbb{m}}),$
(63)

leading to

$\displaystyle\int_{\mathbb{R}^{d}}d^{(d)}\mathbb{x}f(\mathbb{x})(\rho^{(0)}(\mathbb{x})-\bar{\rho})$
$\displaystyle=\int_{\text{UC}}d^{(d)}\bar{\mathbb{x}}\sum_{\mathbb{m}}(\bar{\mathbb{x}}\cdot\nabla{f}(\mathbb{R}_{\mathbb{m}}))(\rho^{(0)}(\bar{\mathbb{x}})-\bar{\rho}).$
(64)

Now approximating

$\displaystyle\sum_{\mathbb{m}}(\bar{\mathbb{x}}\cdot\nabla{f}(\mathbb{R}_{\mathbb{m}}))\approx\frac{1}{V_{\text{UC}}}\int_{\mathbb{R}^{d}}d^{(d)}\mathbb{y}(\bar{\mathbb{x}}\cdot\nabla_{\mathbb{y}}{f}(\mathbb{y}))$
$\displaystyle=\frac{1}{V_{\text{UC}}}\int_{\mathbb{R}^{d}}d^{(d)}\mathbb{y}\nabla_{\mathbb{y}}\cdot(\bar{\mathbb{x}}{f}(\mathbb{y}))=0,$
(65)

where in the last step we used Gauss’ divergence theorem.

##  Appendix C Contour integral representation

First, we rewrite Eq. (32) as

$\displaystyle\delta{Q}=$ $\displaystyle-\frac{1}{2\pi
i}\int_{-\infty}^{\mu}d\omega\frac{\partial}{\partial\omega}\log\det\left\\{\mathcal{M}(z)\right\\}$
$\displaystyle+\frac{1}{2\pi
i}\int_{-\infty}^{\mu}d\omega\frac{\partial}{\partial\omega}\left(\log\det\left\\{\mathcal{M}(z)\right\\}\right)^{*}.$
(66)

Now we remind ourselves that

$\displaystyle\left(\log f(z)\right)^{*}=\log\left(f(z)\right)^{*}\equiv\log
f^{*}(z).$ (67)

Furthermore, one has

$\displaystyle\left(\det\left\\{\mathcal{M}(z)\right\\}\right)^{*}=\det\left\\{\mathcal{M}^{\dagger}(z^{*})\right\\},$
(68)

where, as before, the Hermitian conjugate does not affect the $z$-variable.
Now we have

$\displaystyle\det\left\\{\mathcal{M}^{\dagger}(z^{*})\right\\}=\det\left\\{\right(1-\mathcal{G}^{(0)}(z^{*})\tilde{\mathcal{V}}^{(0)}\left){}^{\dagger}\right\\}$
$\displaystyle=\det\left\\{1-\tilde{\mathcal{V}}^{(0)}\mathcal{G}^{(0)}(z^{*})\right\\}$
$\displaystyle=\det\left\\{\mathcal{M}(z^{*})\right\\},$ (69)

where to get from the pre-last to the last lines we employed the
Weinstein–Aronszajn identity.

It hence follows that

$\displaystyle\delta{Q}=$ $\displaystyle-\frac{1}{2\pi
i}\int_{-\infty}^{\mu}d\omega\frac{\partial}{\partial\omega}\log\det\left\\{\mathcal{M}(\omega+i\eta)\right\\}$
$\displaystyle-\frac{1}{2\pi
i}\int_{\mu}^{-\infty}d\omega\frac{\partial}{\partial\omega}\log\det\left\\{\mathcal{M}(\omega-i\eta)\right\\}$
$\displaystyle=$ $\displaystyle-\oint_{C}\frac{dz}{2\pi
i}\frac{\partial}{\partial z}\log\det\left\\{\mathcal{M}(z)\right\\}.$ (70)

Above, $C$ is the counterclockwise rectangular contour defined as a union of
four segments:

$\displaystyle C=$
$\displaystyle[B+i\eta,\mu+i\eta)\cup[\mu+i\eta,\mu-i\eta)\cup[\mu-i\eta,B-i\eta)$
$\displaystyle\cup[B-i\eta,B+i\eta),\quad
B\rightarrow-\infty,\quad\eta\rightarrow 0^{+}.$ (71)

We note that the integral in (70) remains unaffected under continuous contour
deformations, so long as the analytic structure of the integrand within the
patch of the complex plane enclosed by contour $C$ remains intact. In this
connection, we may replace $C$ with an arbitrary non-self-intersecting curve
crossing the real axis at two points only, at any energy below the lowest
eigenvalue of the full Hamiltonian, and at the chemical potential.

##  Appendix D Parameters and protocols

In the numerical example provided in Section III.2, the parameters of the
model were chosen according to

$\displaystyle\mathbb{q}$
$\displaystyle=\frac{2\pi(\mathbb{e}_{x}+\kappa\mathbb{e}_{y})}{\lambda\sqrt{1+\kappa^{2}}},\quad\mathbb{A}_{0}=\frac{A_{0}(\kappa\mathbb{e}_{x}-\mathbb{e}_{y})}{\sqrt{1+\kappa^{2}}},\quad
m=1,$ (72) $\displaystyle\kappa$
$\displaystyle=\frac{1+\sqrt{5}}{2},\quad\lambda=4,\quad\frac{e}{c}A_{0}=1.17,\quad\mathbb{k}_{R}=\begin{pmatrix}0.32\\\
1.39\\\ 1.24\end{pmatrix}.$ (73)

Note that as we have set the electron’s mass $m=1$ to unity (in addition to
the electric charge $e=1$ and reduced Plank’s constant $\hbar=1$), we work in
Hartree’s atomic units.

In this way, the electromagnetic wave is propagating in the $x-y$ plane, with
the corresponding magnetic field being

$\displaystyle\mathbb{B}(\mathbb{x})=\frac{2\pi
A_{0}}{\lambda}\mathbb{e}_{z}\sin(\mathbb{q}\cdot\mathbb{x}+\varphi).$ (74)

By definition, the corresponding lattice period is given by

$\displaystyle L=\lambda\sqrt{1+\kappa^{2}}=2\sqrt{2(5+\sqrt{5})}.$ (75)

To produce the data, we used the following pumping protocol for the
impurities’ separation

$\displaystyle
R(\tau)=(n_{R}-1)L+\bar{R}(\tau),\quad\bar{R}(\tau)=\frac{L}{T}\tau,\quad\frac{L}{T}\ll
v_{F},$ (76)

where $v_{F}$ is the Fermi velocity and $n_{R}$ is an integer specifying the
number of the unit cell hosting the second impurity. For the impurities’
vertex functions, we further make an assumption of the equivalent impurities:
$U^{(1)}(\tau)=U^{(2)}(\tau)=:U(\tau)$ and
$|\mathbb{B}_{\text{eff}}^{(1)}(\tau)|=|\mathbb{B}_{\text{eff}}^{(2)}(\tau)|=:B_{I}(\tau)$.
The direction of the magnetic moments, on the other hand, is allowed to be
different in two scattering centers and is parametrized in spherical polar
coordinates

$\displaystyle\frac{\mathbb{B}_{\text{eff}}^{(j)}(\tau)}{B_{I}(\tau)}=\begin{pmatrix}\cos\phi^{(j)}(\tau)\sin\theta^{(j)}(\tau)\\\
\sin\phi^{(j)}(\tau)\sin\theta^{(j)}(\tau)\\\
\cos\theta^{(j)}(\tau)\end{pmatrix}.$ (77)

In the following, we assume that, as is the case with the location of the
second impurity within the unit cell number $n_{R}$, the impurity strength
also grows linearly with $\tau$

$\displaystyle
U(\tau)=U_{0}+\delta{U}\frac{\tau}{T},\quad\frac{\delta{U}}{T}\ll\epsilon_{F}^{2}.$
(78)

On the other hand, we assume the effective magnetic field of the impurity to
oscillate as

$\displaystyle
B_{I}(\tau)=B_{0}+\delta{B}\sin\left(\frac{6\pi\tau}{T}\right).$ (79)

The direction of the spins is prescribed by

$\displaystyle\phi^{(1)}(\tau)$
$\displaystyle=\phi^{(2)}(\tau)=2\pi\sin\left(\frac{8\pi\tau}{T}\right),$ (80)
$\displaystyle\theta^{(j)}(\tau)$
$\displaystyle=\frac{\pi}{2}\left(1+(-1)^{j}\frac{\tau}{T}\right).$ (81)

The rest of the parameters are chosen as

$\displaystyle
U_{0}=0,\quad\delta{U}=10,\quad\frac{e}{c}B_{0}=3,\quad\frac{e}{c}\delta{B}=1.5.$
(82)

Now let us estimate the charge localization length $\xi_{g}$ for the second
bulk spectral gap, where the chemical potential $\mu$ is assumed to be placed.
According to Ref. [Pias_2022, ], the Fermi velocity may roughly be estimated
as $v_{F}\approx\frac{k_{F}}{m}\approx\frac{2\pi}{mL}\approx 0.825816$. The
energy gap at the Fermi level was numerically computed to be roughly
$E_{g}\approx 0.271394$, leading to the following estimate $\xi_{g}\approx
3=O(L)$.

## References

  * (1) K. von Klitzing, G. Dorda, and M. Pepper, “New method for high-accuracy determination of the fine-structure constant based on quantized Hall resistance”, Phys. Rev. Lett. 45, 494 (1980).
  * (2) D. J. Thouless, M. Kohmoto, M. P. Nightingale, and M. den Nijs, “Quantized Hall conductance in a two-dimensional periodic potential”, Phys. Rev. Lett. 49, 405 (1982).
  * (3) Y. Hatsugai, “Chern number and edge states in the integer quantum Hall effect”, Phys. Rev. Lett. 71, 3697 (1993); Y. Hatsugai, “Edge states in the integer quantum Hall effect and the Riemann surface of the Bloch function”, Phys. Rev. B 48, 11851, (1993).
  * (4) M. Z. Hasan and C. L. Kane, “Colloquium: topological insulators”, Rev. Mod. Phys. 82, 3045 (2010).
  * (5) X.-L. Qi and S.-C. Zhang, “Topological insulators and superconductors”, Rev. Mod. Phys. 83, 1057 (2011).
  * (6) A. Kitaev, “Periodic table for topological insulators and superconductors”, AIP Conf. Proc. 1134, 22 (2009).
  * (7) M. Pletyukhov, D. M. Kennes, J. Klinovaja, D. Loss, and H. Schoeller, “Surface charge theorem and topological constraints for edge states: An analytical study of one-dimensional nearest-neighbor tight-binding models” Phys. Rev. B 101, 165304 (2020); M. Pletyukhov, D. M. Kennes, J. Klinovaja, D. Loss, and H. Schoeller, “Topological invariants to characterize universality of boundary charge in one-dimensional insulators beyond symmetry constraints”, Phys. Rev. B 101, 161106(R) (2020).
  * (8) N. Müller, K. Piasotski, D. M. Kennes, H. Schoeller, and M. Pletyukhov, “Universal properties of boundary and interface charges in multichannel one-dimensional models without symmetry constraints”, Phys. Rev. B 104, 125447 (2021).
  * (9) S. Miles, D. M. Kennes, H. Schoeller, and M. Pletyukhov, “Universal properties of boundary and interface charges in continuum models of one-dimensional insulators”, Phys. Rev. B 104, 155409 (2021).
  * (10) K. Piasotski, N. Müller, D. M. Kennes, H. Schoeller, and M. Pletyukhov, “Universal properties of boundary and interface charges in multichannel one-dimensional continuum models”, Phys. Rev. B 106, 165405 (2022).
  * (11) W. Kohn and A. Yaniv, “Locality principle in wave mechanics”, PNAS 75(11), 5270 (1978).
  * (12) W. Kohn, “Density Functional and Density Matrix Method Scaling Linearly with the Number of Atoms”, Phys. Rev. Lett. 76, 3168 (1996).
  * (13) E. Prodan and W. Kohn, “Nearsightedness of electronic matter”, PNAS 102, 11635 (2005).
  * (14) E. Prodan, “Nearsightedness of electronic matter in one dimension”, Phys. Rev. B 73, 085108 (2006).
  * (15) H.-R. Trebin, “The topology of non-uniform media in condensed matter physics”, Adv. Phys. 31, 195 (1982).
  * (16) M. Pletyukhov, D. M. Kennes, K. Piasotski, J. Klinovaja, D. Loss, and H. Schoeller, “Rational boundary charge in one-dimensional systems with interaction and disorder”, Phys. Rev. Res. 2, 033345 (2020).
  * (17) J. Goldstone and F. Wilczek, “Fractional Quantum Numbers on Solitons”, Phys. Rev. Lett. 47, 986 (1981).
  * (18) C. S. Weber, K. Piasotski, M. Pletyukhov, J. Klinovaja, D. Loss, H. Schoeller, and D. M. Kennes, “Universality of Boundary Charge Fluctuations”, Phys. Rev. Lett. 126, 016803 (2021).
  * (19) K. Nomura and N. Nagaosa, “Electric Charging of Magnetic Textures on the Surface of a Topological Insulator”, Phys. Rev. B 82, 161401(R) (2010).
  * (20) This includes models featuring a momentum-cubic spin-orbit interaction, or even tight-binding models, describing crystals in terms of incomplete basis sets of localized Wannier orbitals.
  * (21) D. K. Park, “Green’s‐function approach to two‐ and three‐dimensional delta‐function potentials and application to the spin‐$1/2$ Aharonov–Bohm problem”, J. Math. Phys. 36, 5453 (1995).
  * (22) D. A. Atkinson, H. W. Crater, “An exact treatment of the Dirac delta function potential in the Schrödinger equation”, Amer. Jour. Phys. 43, 301 (1975).
  * (23) R. Jackiw, “Delta-function potentials in two- and three-dimensional quantum mechanics” MAB Bég memorial volume (1991).
  * (24) C. N. Friedman, “Perturbations of the Schrödinger equation by potentials with small support”, J. Funct. Anal. 10, 346 (1972).
  * (25) S. Kivelson and J. R. Schrieffer, “Fractional charge, a sharp quantum observable”, Phys. Rev. B 25, 6447 (1982).
  * (26) In the case of one-dimensional ($d=1$) single-channel ($N_{c}=1$) systems, it is not possible for the bound state poles to coexist with the continuum (or extended) states. As a result, the effect of the bound state pole merging with the band is to increase the order of the branching pole at the band edge. In multichannel ($N_{c}>1$) and/or multidimensional ($d>1$) systems, on the other hand, the poles may, in principle, coexist with the band continuum. In either case, as a result, the band contribution to the winding number gets modified by unity.
  * (27) Z. Hou, C. S. Weber, D. M. Kennes, D. Loss, H. Schoeller, J. Klinovaja, M. Pletyukhov, “Realization of a three-dimensional quantum Hall effect in a Zeeman-induced second order topological insulator on a torus”, arXiv preprint arXiv:2212.09053, (2022)