# Universal properties of boundary and interface charges in multichannel one-
dimensional models without symmetry constraints

Niclas Müller  Institut für Theorie der Statistischen Physik, RWTH Aachen,
52056 Aachen, Germany and JARA - Fundamentals of Future Information Technology
Kiryl Piasotski  Institut für Theorie der Statistischen Physik, RWTH Aachen,
52056 Aachen, Germany and JARA - Fundamentals of Future Information Technology
Dante M. Kennes  Institut für Theorie der Statistischen Physik, RWTH Aachen,
52056 Aachen, Germany and JARA - Fundamentals of Future Information Technology
Max Planck Institute for the Structure and Dynamics of Matter, Center for Free
Electron Laser Science, 22761 Hamburg, Germany  Herbert Schoeller  Institut
für Theorie der Statistischen Physik, RWTH Aachen, 52056 Aachen, Germany and
JARA - Fundamentals of Future Information Technology  Mikhail Pletyukhov
<EMAIL_ADDRESS>Institut für Theorie der Statistischen Physik,
RWTH Aachen, 52056 Aachen, Germany and JARA - Fundamentals of Future
Information Technology

###### Abstract

The boundary charge that accumulates at the edge of a one-dimensional single-
channel insulator is known to possess the universal property, that its change
under a lattice shift towards the edge by one site is given by the sum of the
average bulk electronic density and a topologically invariant contribution,
restricted to the values $0$ and $-1$ [Phys. Rev. B 101, 165304 (2020)]. This
quantized contribution is associated with particle-hole duality, ensures
charge conservation and fixes the mod(1) ambiguity appearing in the Modern
Theory of Polarization. In the present work we generalize the above-mentioned
single-channel results to the multichannel case by employing the technique of
boundary Green’s functions. We show that the topological invariant associated
with the change in boundary charge under a lattice shift in multichannel
models can be expressed as a winding number of a certain combination of
components of bulk Green’s functions as function of the complex frequency, as
it encircles the section of the energy axis that corresponds to the occupied
part of the spectrum. We observe that this winding number is restricted to
values ranging from $-N_{c}$ to $0$, where $N_{c}$ is the number of channels
(orbitals) per site. Furthermore, we consider translationally invariant one-
dimensional multichannel models with an impurity and introduce topological
indices which correspond to the quantized charge that accumulates around said
impurity. These invariants are again given in terms of winding numbers of
combinations of components of bulk Green’s functions. Through this
construction we provide a rigorous mathematical proof of the so called
nearsightedness principle formulated by W. Kohn [Phys. Rev. Lett. 76, 3168
(1996)] for noninteracting multichannel lattice models.

##  I Introduction

Over the last decades, the study of topological insulators, sparked by the
discoveries of the quantum Hall effect [Klitzing_1980, , Thouless_1982, ] and
later of the quantum spin Hall effect [Konig_2007, ], has attracted much
interest. Topological insulators (TI) are states of matter with a gapped bulk
and symmetry-protected gapless edge states [He_2019, ]. These dissipationless
edge states are at the center of interest in the field of TIs, due to their
numerous promising applications in quantum computing [Bomantara_2018, ;
Mong_2014, ; Stern_2013, ; Miyake_2010, ] and spintronics [He_2019, ,
Pesin_2012, ], and their exponential localization poses questions concerning
the nature of the charge distribution in their spatial vicinity
[pletyukhov_etal_prb_20, ]. As particularly interesting in this respect appear
both the boundary charge $Q_{B}$, i.e. the charge that accumulates at the edge
of a finite insulator, and the interface charge $Q_{I}$, i.e. the charge that
accumulates around the interface shared by a pair of insulators
[Pletyukhov_etal3_2020, ]. The study of both of these quantities has a long
history [jackiw_rebbi_prd_76, ; su_schrieffer_heeger, ;
jackiw_schrieffer_nuclphys_81, ; rice_mele_prl_82, ; kivelson_prb_83, ;
su_schrieffer_prl_81, ; heeger_etal_review_RMP_88, ; witten_physlett_79, ;
takayama_etal_prb_80, ; goldstone_wilczek_prl_81, ; jackiw_semenoff_prl_83, ],
connecting them to the field of TIs however is a rather young endeavor.

With the emergence of the so-called modern theory of polarization (MTP),
interest in the study of the boundary charge reawakened [King_1993, ;
Vanderbilt_1993, ; Resta_1994, ; Kudin_2007, ; Marzari_2012, ; Spaldin-2012, ;
Rhim_2017, ; Miert, ]. One of the major achievements of the MTP lies in the
so-called surface charge theorem, relating the surface (boundary) charge to
the bulk polarization, which in turn is related [King_1993, , Vanderbilt_1993,
] to the Zak-Berry phase [zak_89, ]. However, the MTP is restricted by the
fact that the Zak-Berry phase is defined modulo an integer, since upon an
Abelian gauge transformation it changes by the winding number of the
corresponding phase. This fact complicates bridging the field of TIs with MTP,
since MTP is unable to predict the number of edge states, the central quantity
in the field of TIs.

Topological indices related to the boundary charge in one-dimensional, single
channel models without symmetry constraints were introduced recently
[pletyukhov_etal_prb_20, , Pletyukhov_etal2_2020, ]. In particular it was
shown, that knowledge of the exact eigenstates of a given model allows to find
the gauge in which there is a unique relation between the boundary charge
$Q_{B,\alpha}$ of a given band $\alpha$ and the corresponding Zak-Berry phase,
fixing the unknown integer of the surface charge theorem. A new topological
invariant underpinning the universal behavior of the boundary charge upon
continuous shifting of the lattice towards the boundary was introduced
[pletyukhov_etal_prb_20, , Pletyukhov_etal2_2020, ]. It was shown, that the
boundary charge $Q_{B}$ is connected to universal long-wavelength properties
of topological insulators, making the study of the boundary charge an
invaluable tool for the characterization of TIs beyond symmetry constraints
[pletyukhov_etal_prb_20, ; Park_etal_2016, ; Thakurathi_etal_2018, ;
Pletyukhov_etal2_2020, ; Lin_etal_2020, ; Weber_etal_2021, ; Lin_etal_2021, ;
Laubscher_etal_2021, ].

In another recent work [Pletyukhov_etal3_2020, ] the universal properties of
boundary and interface charges were embedded in a general framework, relating
them generically to the nearsightedness principle [Kohn_1996, ,
Prodan_Kohn_2005, ]. This principle states a fundamental property of
insulators that local perturbations by external fields lead only to local
charge redistributions resulting in an addition/removal of an integer number
of electrons to/from a perturbed region. Based on this theorem two invariants
were established, with quantized values in generic systems, including
interactions, random disorder, and multi-channel systems. The first invariant
is related to the change of the boundary charge by discrete translations of
the lattice, similar to the one described above for single-channel systems.
The second invariant is related to the sum of the two boundary charges left
and right to a barrier separating a translational invariant lattice in two
subparts. The quantization of this invariant is related to a local inversion
of the lattice and, together with the first invariant, provides the basis for
the quantization of interface charges at domain walls [jackiw_rebbi_prd_76, ;
su_schrieffer_heeger, ; jackiw_schrieffer_nuclphys_81, ; rice_mele_prl_82, ;
kivelson_prb_83, ; su_schrieffer_prl_81, ; heeger_etal_review_RMP_88, ] and
the generalization of the Goldstone-Wilczek formula [witten_physlett_79, ;
takayama_etal_prb_80, ; goldstone_wilczek_prl_81, ; jackiw_semenoff_prl_83, ]
to arbitrary tight-binding lattices.

The topic of the present paper is to analyse the two invariants introduced in
Ref. [Pletyukhov_etal3_2020, ] in all detail for noninteracting multi-channel
tight-binding models. Via an explicit representation in terms of a winding
number defined purely in terms of bulk quantities of the infinite system, we
rigorously prove the quantization and provide a mean to explicitly calculate
their integer values from bulk properties. This generalizes the results of
Refs. [pletyukhov_etal_prb_20, , Pletyukhov_etal2_2020, ] to systems with more
than a single orbital per lattice site, establishes a bulk-boundary
correspondence between universal properties of the boundary charge and bulk
properties, and provides a rigorous proof of the nearsightedness principle.

To calculate the invariants we use the method of boundary Green’s functions
(BGF) [Arrachea_etal_2009, ; Zazunov_etal_2016, ; Peng_etal_2017, ;
Zazunov_etal_2017, ; Komnik_etal_2017, ; Zazunov_etal_2018, ; Pinon_etal_2020,
; Alvarado_etal_2020, ]. Complementary to the method of
[pletyukhov_etal_prb_20, ] this technique allows for the direct construction
of the open boundary lattice Green’s function, avoiding the computation of
eigenstates [Komnik_etal_2017, ]. The BGF method proved to be useful in
analyzing electronic transport in superconducting systems [Zazunov_etal_2017,
, Zazunov_etal_2018, ], thermal transport in spin heterostructures
[Arrachea_etal_2009, ], as well as transport in topological superconductors
hosting Majorana bound states [Zazunov_etal_2016, , Alvarado_etal_2020, ].
Lately, the boundary Green’s function technique was generalized to higher
dimensional systems [Pinon_etal_2020, ] where the boundaries manifest as lines
and planes in two and three dimensions, respectively. Other noteworthy
applications of the BGF formalism include the efficient construction of
topological phase diagrams [Peng_etal_2017, ], as well as bulk-boundary
correspondence related investigations in both non-interacting
[Essin_Gurarie_2011, ] and interacting topological insulators [Gurarie_2011,
]. Similar methods were also extensively used in the 60s and 70s with respect
to determining the electronic spectrum as well as conductivity in one-
dimensional metallic systems with strong disorder [bychkov_1966, ;
Bychkov_1967, ; Kazamanyan_19721, ].

In the present work the method of boundary Green’s functions is used to study
boundary and interface charges and their associated topological invariants in
multichannel extensions of generalized Aubry-André-Harper models [aa_model, ,
harper_model, ]. Employing the BGF technique, we express the boundary charge
of a one-dimensional multichannel insulator in terms of an energy integral of
the local spectral density which, in turn, is expressed via the bulk Green’s
functions. Such a representation is beneficial since it avoids the direct
diagonalization of the semi-infinite system’s Hamiltonian. We introduce the
topological invariant $I$ associated with the change of boundary charge
$\Delta{Q}_{B}$ (relative to the average charge per site in the bulk
$\bar{\rho}$) under lattice translations. This invariant is expressed as a
winding number of a particular combination of components of bulk Green’s
functions, clearly demonstrating the bulk-boundary correspondence discussed
within the realm of topological insulators [Fidkowski_etal_2011, ;
Mong_Shivamoggi_2011, ; Fukui_etal_2012, ; Yu_etal_2017, ]. When the whole
lattice is shifted by a single site some number of edge states (each
contributing a unit of electron charge) may either cross the chemical
potential from above or below respectively, hence maintaining the integer
valuedness of this invariant. Similar arguments are frequently used in
discussions of adiabatic charge pumping [Thouless_83, , Niu_Thouless_1984, ].
As opposed to the single channel case, where this invariant is limited to the
values $0$ and $-1$, in a system with $N_{c}$ orbitals per site $I$ can take
on integers in the range $-N_{c}\leq I\leq 0$, naturally generalizing the
$N_{c}=1$ case (see Refs. [pletyukhov_etal_prb_20, , Pletyukhov_etal2_2020,
]). It is worth noting that the form of bulk-boundary correspondence
encompassed by this invariant is not limited to any particular symmetry
classes and is purely based on such physical principles as charge conservation
and Pauli exclusion principle. To the best of our knowledge this invariant was
not established rigorously in the literature and thus provides a novel
contribution to the fields of topological insulators, boundary Green’s
functions and the modern theory of polarization.

We also adopt the BGF formalism to study the quantization of interface charge
and establish rigorously another invariant related to the sum of the boundary
charges left and right to a barrier. In particular, we study two types of
interface models: The first is obtained from a translationally invariant model
by adding a finite potential barrier on a given site, whereas the second type
is obtained from a translationally invariant model by weakening the link
between two adjacent sites. We derive expressions for the interface charge in
both impurity and weak link models in terms of winding numbers of combinations
of components of bulk Green’s functions, thus demonstrating the quantization
of the interface charge. Since the unit cells are perfectly matched, there is
no dipole moment in the interface charge distribution, which is thus generated
solely by removing a number of valence band electrons to the emerging
interface-localized bound states. This result is of fundamental importance in
the theory of topological insulators, since it provides the direct analytical
proof of the nearsightedness principle [Kohn_1996, , Prodan_Kohn_2005, ]. This
principle assigns the robustness property to excess charges localized near a
boundary or at an interface of two insulators (cf. [Pletyukhov_etal3_2020, ]
for the analytic proof of the interface charge quantization in a special
case).

Finally, we substantiate our findings with a number of numerical examples,
comparing the winding numbers with their physical counterparts (i.e. boundary
and interface charges) as computed from exact diagonalization of finite
systems. We use randomly generated Hamiltonians for this comparison,
demonstrating the validity of these novel invariants.

This paper is organized as follows. In Sec. II we introduce the class of
models under consideration and define their boundary Green’s functions. In
Sec. III we define the boundary charge for multichannel one-dimensional models
and express it in terms of boundary Green’s functions. We then use results
from Sec. II to express the change in boundary charge upon lattice shifts via
the winding number of a specific combination of components of bulk Green’s
functions. In Sec. IV we define the interface charge for two different
impurity models: 1) finite potential on a single site and 2) a weakened link
between two adjacent unit cells. We express the interface charge in both cases
in terms of boundary Green’s functions and use this representation to cast the
expression for the interface charge in the form of a winding number. In Sec. V
we explain how to efficiently evaluate the various winding numbers and then
show examples demonstrating the validity of all three invariants as applied to
randomly generated models. Finally, in Sec. VI we state our summary.

##  II Boundary Green’s function

Let us consider a class of translationally invariant one-dimensional (1D)
lattice models with $Z$ sites per unit cell (labelled by $1\leq j\leq Z$ in
the following) and $N_{c}$ states per site (called channels and labelled by
$\sigma$). We use the global coordinate index $m=Z(n-1)+j\sim(n,j)$, which as
well contains the unit cell index $n$.

In particular, we focus on the nearest-neighbour hopping models

$\displaystyle\hat{H}_{0}=$
$\displaystyle-\sum_{m=-\infty}^{\infty}\sum_{\sigma,\sigma^{\prime}=1}^{N_{c}}(t_{m,\sigma\sigma^{\prime}}\,\,|m+1,\sigma\rangle\langle
m,\sigma^{\prime}|+h.c.)$ (1)
$\displaystyle+\sum_{m=-\infty}^{\infty}\sum_{\sigma,\sigma^{\prime}=1}^{N_{c}}v_{m,\sigma\sigma^{\prime}}\,\,|m,\sigma\rangle\langle
m,\sigma^{\prime}|,$ (2) $\displaystyle t_{m}$ $\displaystyle=t_{m+Z},\quad
v_{m}=v_{m}^{\dagger}=v_{m+Z},$ (3)

where $t_{m}$ and $v_{m}$ are hopping and potential matrices, respectively,
both of the size $N_{c}\times N_{c}$.

Solutions $|\psi_{k\alpha}^{(0)}\rangle$ of the eigenvalue problem
$\hat{H}_{0}|\psi_{k\alpha}^{(0)}\rangle=\epsilon_{k\alpha}^{(0)}|\psi_{k\alpha}^{(0)}\rangle$
are labelled by the Bloch momentum $k$, $-\pi\leq k<\pi$, and the band index
$\alpha$, $1\leq\alpha\leq N_{c}Z$. The Bloch Hamiltonian

$\displaystyle
h_{k}=\left(\begin{array}[]{cccccc}v_{1}&-t_{1}^{\dagger}&0&\ldots&0&-t_{Z}e^{-ik}\\\
-t_{1}&v_{2}&-t_{2}^{\dagger}&\ldots&0&0\\\ 0&-t_{2}&&&\vdots&\vdots\\\
\vdots&\vdots&&&-t_{Z-2}^{\dagger}&0\\\
0&0&\ldots&-t_{Z-2}&v_{Z-1}&-t_{Z-1}^{\dagger}\\\
-t_{Z}^{\dagger}e^{ik}&0&\ldots&0&-t_{Z-1}&v_{Z}\end{array}\right)$ (10)

is related to $\hat{H}_{0}$ via

$\displaystyle\hat{H}_{0}=\int_{-\pi}^{\pi}dk\,\,|k\rangle\langle k|\otimes
h_{k},\quad\langle n|k\rangle=\frac{e^{ikn}}{\sqrt{2\pi}},$ (11)

with
$h_{k}=\sum_{jj^{\prime}\sigma\sigma^{\prime}}(h_{k})_{j\sigma,j^{\prime}\sigma^{\prime}}|j,\sigma\rangle\langle
j^{\prime},\sigma^{\prime}|$. The normalized eigenstates obeying
$h_{k}|\chi_{k\alpha}\rangle=\epsilon_{k\alpha}^{(0)}|\chi_{k\alpha}\rangle$,
$\langle\chi_{k\alpha}|\chi_{k\alpha^{\prime}}\rangle=\delta_{\alpha\alpha^{\prime}}$,
help us express

$\displaystyle\langle
m,\sigma|\psi_{k\alpha}^{(0)}\rangle=\frac{e^{ikn}}{\sqrt{2\pi}}\langle
j,\sigma|\chi_{k\alpha}\rangle.$ (12)

We note the completeness relation and the identity resolution

$\displaystyle\hat{1}=\sum_{\alpha}\int_{-\pi}^{\pi}dk\,\,|\psi_{k\alpha}^{(0)}\rangle\langle\psi_{k\alpha}^{(0)}|=\sum_{m,\sigma}|m,\sigma\rangle\langle
m,\sigma|.$ (13)

Using the notation $\psi_{k\alpha}^{(0)}(m,\sigma)=\langle
m,\sigma|\psi_{k\alpha}^{(0)}\rangle$, we also quote the component-wise form
of (13)

$\displaystyle\delta_{mm^{\prime}}\delta_{\sigma\sigma^{\prime}}=\sum_{\alpha}\int_{-\pi}^{\pi}dk\,\,\psi_{k\alpha}^{(0)}(m,\sigma)\,\,\psi_{k\alpha}^{(0)\,*}(m^{\prime},\sigma^{\prime}).$
(14)

The bulk retarded Green’s function
$\hat{G}^{(0)}(\omega)=\frac{1}{\omega-\hat{H}_{0}+i\eta}$ has the matrix
expression

$\displaystyle G^{(0)}_{m\sigma,m^{\prime}\sigma^{\prime}}$
$\displaystyle=\langle
m,\sigma|\hat{G}^{(0)}(\omega)|m^{\prime},\sigma^{\prime}\rangle$
$\displaystyle=\int_{-\pi}^{\pi}\frac{dk}{2\pi}\,\,\langle
j,\sigma|\frac{e^{ik(n-n^{\prime})}}{\omega+i\eta-
h_{k}}|j^{\prime},\sigma^{\prime}\rangle.$ (15)

In the following we make use of the reduced notations

$\displaystyle
G^{(0)}_{m,m^{\prime}}=\int_{-\pi}^{\pi}\frac{dk}{2\pi}\,\,\langle
j|\frac{e^{ik(n-n^{\prime})}}{\omega+i\eta-h_{k}}|j^{\prime}\rangle\equiv
G^{(0)}_{j,j^{\prime}}(n-n^{\prime}),$ (16)

implying that each element $[(\omega+i\eta-h_{k})^{-1}]_{jj^{\prime}}$ is a
$N_{c}\times N_{c}$ matrix block with internal indices
$\sigma,\sigma^{\prime}$.

Adding an arbitrary potential $\hat{V}$ to $\hat{H}_{0}$, we break the
translational invariance of the 1D lattice model. The retarded Green’s
function $\hat{G}(\omega)=\frac{1}{\omega+i\eta-\hat{H}}$ of the perturbed
model $\hat{H}=\hat{H}_{0}+\hat{V}$ satisfies the Dyson equation

$\displaystyle\hat{G}(\omega)=\hat{G}^{(0)}(\omega)+\hat{G}^{(0)}(\omega)\hat{V}\hat{G}(\omega).$
(17)

To mimic a boundary model which is defined in the right half-space ($m\geq
1$), we can choose an infinitely high potential for $m\leq 0$, which would
block an occupation of sites in the left half-space. In models which only
allow for nearest neighbour hopping it is however sufficient to put a high
potential just on the single site $m=0$: It will play the role of the
impenetrable barrier between the right and left half-spaces. In the following
we restrict ourselves to this class of models and choose $\langle
m,\sigma|\hat{V}|m^{\prime},\sigma^{\prime}\rangle=V_{0}\,\,\delta_{\sigma\sigma^{\prime}}\,\delta_{m,0}\,\delta_{m^{\prime},0}$,
aiming to perform the limit $V_{0}\to\infty$ afterwards. A case of longer
ranged hoppings will be discussed elsewhere [Pias_new, ].

Choosing the ultra-local potential of the above stated form has the advantage
that this choice allows us to study the left boundary model with $m\leq-1$ in
the same setting. The Dyson equation (17) acquires then the special form

$\displaystyle
G_{m\sigma,m^{\prime}\sigma^{\prime}}=G_{m\sigma,m^{\prime}\sigma^{\prime}}^{(0)}+G_{m\sigma,0\sigma_{1}}^{(0)}V_{0}\,\,G_{0\sigma_{1},m^{\prime}\sigma^{\prime}},$
(18)

where we omit for brevity the $\omega$-dependence and implicitly assume a
summation over the repeated index $\sigma_{1}$ (this convention is also used
in the following). In the matrix notation, Eq. (18) is equivalent to

$\displaystyle
G_{m,m^{\prime}}(\omega)=G_{m,m^{\prime}}^{(0)}(\omega)+V_{0}\,\,G_{m,0}^{(0)}(\omega)\,G_{0,m^{\prime}}(\omega).$
(19)

To solve this equation for $G_{m,m^{\prime}}(\omega)$, we first set $m=0$ and
find

$\displaystyle
G_{0,m^{\prime}}(\omega)=[1-V_{0}\,G_{0,0}^{(0)}(\omega)]^{-1}\,G_{0,m^{\prime}}^{(0)}(\omega).$
(20)

Inserting this result back into (19) yields the expression for
$\hat{G}(\omega)$ in terms of $\hat{G}^{(0)}(\omega)$:

$\displaystyle G_{m,m^{\prime}}(\omega)$
$\displaystyle=G_{m,m^{\prime}}^{(0)}(\omega)$
$\displaystyle+V_{0}\,G_{m,0}^{(0)}(\omega)[1-V_{0}\,G_{0,0}^{(0)}(\omega)]^{-1}\,G_{0,m^{\prime}}^{(0)}(\omega).$
(21)

Now it is appropriate to take the limit of the infinite barrier height
$V_{0}\to\infty$. It leads to the so called boundary Green’s function

$\displaystyle G_{m,m^{\prime}}(\omega)$
$\displaystyle=G_{m,m^{\prime}}^{(0)}(\omega)$ $\displaystyle-
G_{m,0}^{(0)}(\omega)[G_{0,0}^{(0)}(\omega)]^{-1}G_{0,m^{\prime}}^{(0)}(\omega).$
(22)

##  III Boundary charge and the associated topological invariant

In this section we study an excess charge which is accumulated near the hard-
wall boundary of a semi-infinite lattice. In particular, we derive an
expression for a change $\Delta Q_{B}$ of the boundary charge $Q_{B}$ under
the lattice shift by one site towards the wall in terms of a topological
invariant $I$. The latter is given by a winding number and therefore takes
integer values. Thereby we achieve a multichannel generalization of our
earlier single-channel result [pletyukhov_etal_prb_20, ,
Pletyukhov_etal2_2020, ]. However, instead of explicitly constructing the
boundary problem eigenstates, which was feasible in the single-channel
consideration, we resort now to the representation in terms of the boundary
Green’s functions, which was introduced in the previous section.

###  III.1 Boundary charge definition

The solution of the boundary model reads

$\displaystyle\hat{H}|\psi_{s}\rangle=\epsilon_{s}|\psi_{s}\rangle.$ (23)

Here, $s$ is a general index to label all eigenstates, and it can have both
continuous and discrete domains with extended and localized eigenstates,
respectively. We note expressions for the completeness and the boundary
Green’s function in the basis $\psi_{s}(m,\sigma)=\langle
m,\sigma|\psi_{s}\rangle$

$\displaystyle\delta_{mm^{\prime}}\delta_{\sigma\sigma^{\prime}}$
$\displaystyle=\sumint_{s}\,\psi_{s}(m,\sigma)\,\psi_{s}^{*}(m^{\prime},\sigma^{\prime}),$
(24) $\displaystyle G_{m\sigma,m^{\prime}\sigma^{\prime}}(\omega)$
$\displaystyle=\sumint_{s}\,\frac{\psi_{s}(m,\sigma)\,\psi_{s}^{*}(m^{\prime},\sigma^{\prime})}{\omega+i\eta-\epsilon_{s}}.$
(25)

At zero temperature, the charge density equals

$\displaystyle\rho_{m}$
$\displaystyle=\sumint_{s}\sum_{\sigma}\Theta(\mu-\epsilon_{s})\,\psi_{s}^{*}(m,\sigma)\,\psi_{s}(m,\sigma)$
(26) $\displaystyle=\sum_{\sigma}\int
d\omega\,\Theta(\mu-\omega)\sumint_{s}\psi_{s}^{*}(m,\sigma)\,\psi_{s}(m,\sigma)\,\delta(\omega-\epsilon_{s})$
$\displaystyle=-\frac{1}{\pi}\sum_{\sigma}\int
d\omega\,\Theta(\mu-\omega)\,\text{Im}\,G_{m\sigma,m\sigma}(\omega),$ (27)

where $\mu$ is the chemical potential. Inserting the solution (22) into the
above expression we find

$\displaystyle\rho_{m}-\rho_{m}^{(0)}$
$\displaystyle=\frac{1}{\pi}\text{Im}\int d\omega\Theta(\mu-\omega)$
$\displaystyle\times\text{tr}\left\\{G_{m,0}^{(0)}(\omega)[G_{0,0}^{(0)}(\omega)]^{-1}G_{0,m}^{(0)}(\omega)\right\\},$
(28)

where the trace operation is performed in the channel space, and

$\displaystyle\rho_{m}^{(0)}=-\frac{1}{\pi}\text{Im}\int
d\omega\Theta(\mu-\omega)\text{tr}\left\\{G_{m,m}^{(0)}(\omega)\right\\}\equiv\rho_{j}^{(0)}$
(29)

is the density in the translationally invariant model (and thus in the bulk),
which only depends on the site index $j$ within the unit cell.

The boundary charge (in units of electron’s charge) is defined as

$\displaystyle Q_{B}=\sum_{m=1}^{\infty}(\rho_{m}-\bar{\rho}^{(0)})f_{m},$
(30)

where $\bar{\rho}^{(0)}=\frac{1}{Z}\sum_{j=1}^{Z}\rho_{j}^{(0)}$ is the unit-
cell averaged density in the bulk, and $f_{m}$ is an envelope function
mimicking a charge probe (see [pletyukhov_etal_prb_20, ,
Pletyukhov_etal3_2020, ] for details). In particular, $f_{m}\approx 1$ for
$m\lesssim M$, and gradually falls off to zero value on the interval
$M\lesssim m\lesssim M+N$, with $M,N\gg 1$. Splitting

$\displaystyle Q_{B}$
$\displaystyle=\sum_{m=1}^{\infty}(\rho_{m}-\rho_{m}^{(0)})f_{m}$ (31)
$\displaystyle+\sum_{m=1}^{\infty}(\rho_{m}^{(0)}-\bar{\rho}^{(0)})f_{m},$
(32)

one can show that $Q_{B}=Q^{\prime}_{B}+Q_{P}$, with

$\displaystyle Q^{\prime}_{B}$ $\displaystyle=\frac{1}{\pi}\text{Im}\int
d\omega\,\,\Theta(\mu-\omega)$
$\displaystyle\times\text{tr}\left\\{[G_{0,0}^{(0)}(\omega)]^{-1}\sum_{m=1}^{\infty}G_{0,m}^{(0)}(\omega)\,G_{m,0}^{(0)}(\omega)\right\\},$
(33) $\displaystyle Q_{P}$
$\displaystyle=-\sum_{j=1}^{Z}\frac{j}{Z}\left(\rho_{j}^{(0)}-\bar{\rho}^{(0)}\right).$
(34)

The contribution (33) is obtained from (31) by approximating $f_{m}\approx 1$
(which is justified, since (28) decays on the scale of a typical localization
length $\xi\ll M$). This contribution arises from the density modulation close
to the boundary. It contains contributions from the Friedel density
oscillations of extended states as well as integer-valued contributions from
edge states whose energies might reside in spectral gaps. In Appendix A.1 we
derive the expression (111) for the corresponding integrand in which the sum
over $m$ is performed.

In turn, the contribution (34) represents the dipole moment of the unit cell
(also called the polarization charge). It is induced by the spatial variation
of $f_{m}$, which takes place far away from the boundary (see Ref.
[pletyukhov_etal_prb_20, ] for details of deriving (34) from (32)). For the
following analysis it is convenient to express

$\displaystyle\rho_{j}^{(0)}=$ $\displaystyle-\frac{1}{\pi}\int
d\omega\,\,\Theta(\mu-\omega)$ $\displaystyle\times\text{Im
tr}\left\\{\int_{-\pi}^{\pi}\frac{dk}{2\pi}\,\,\langle
j|\frac{1}{\omega+i\eta-h_{k}}|j\rangle\right\\}$ (35) $\displaystyle=$
$\displaystyle\sum_{\sigma}\sum_{\alpha=1}^{\nu}\int_{-\pi}^{\pi}\frac{dk}{2\pi}\,\,|\chi_{k\alpha}(j,\sigma)|^{2},$
(36)

where $\nu$ is a number of the fully occupied bands, which depends on the
level of the chemical potential $\mu$. Using the normalization of the states
$|\chi_{k\alpha}\rangle$ it is also straightforward to show that

$\displaystyle\bar{\rho}^{(0)}$
$\displaystyle=\frac{1}{Z}\sum_{j=1}^{Z}\sum_{\sigma}\sum_{\alpha=1}^{\nu}\int_{-\pi}^{\pi}\frac{dk}{2\pi}\,\,|\chi_{k\alpha}(j,\sigma)|^{2}=\frac{\nu}{Z}.$
(37)

###  III.2 Topological invariant for boundary charge change under the lattice
shift

We shift the lattice by one site towards the boundary and study $\Delta
Q_{B}=\tilde{Q}_{B}-Q_{B}$, where $\tilde{Q}_{B}$ is the boundary charge in
the shifted system. For its expression it is sufficient to replace
$t_{j}\to\tilde{t}_{j}=t_{j+1}$, $v_{j}\to\tilde{v}_{j}=v_{j+1}$ for $1\leq
j\leq Z-1$, as well as $t_{Z}\to\tilde{t}_{Z}=t_{1}$,
$v_{Z}\to\tilde{v}_{Z}=v_{1}$.

First we evaluate the change in the polarization charge (34)

$\displaystyle\Delta Q_{P}$
$\displaystyle=\tilde{Q}_{P}-Q_{P}=-\sum_{j=1}^{Z}\frac{j}{Z}\left(\tilde{\rho}_{j}^{(0)}-\rho_{j}^{(0)}\right)$
(38)
$\displaystyle=-\sum_{j=1}^{Z-1}\frac{j}{Z}\left(\rho_{j+1}^{(0)}-\rho_{j}^{(0)}\right)-\rho_{1}^{(0)}+\rho_{Z}^{(0)}$
(39) $\displaystyle=\frac{\nu}{Z}-\rho_{1}^{(0)}.$ (40)

Defining the quantity

$\displaystyle I=\Delta Q_{B}-\frac{\nu}{Z}=\Delta
Q^{\prime}_{B}-\rho_{1}^{(0)},$ (41)

we show that it can take only integer values and thereby appears to be a
topological invariant. In particular, in Appendix A.3 we derive the winding
number expression

$\displaystyle I=$ $\displaystyle-\frac{1}{\pi}\int
d\omega\,\,\Theta(\mu-\omega)$
$\displaystyle\times\,\text{Im}\,\frac{\partial}{\partial\omega}\,\ln\det\left([G^{(0)}_{0,0}]^{-1}G_{0,1}^{(0)}t_{Z}\right).$
(42)

As is explained in Sec. V, this integral can be interpreted as a winding
number of the function
$\det\left([G^{(0)}_{0,0}]^{-1}G_{0,1}^{(0)}t_{Z}\right)$, as the complex
frequency $\omega$ encircles the occupied part of the spectrum, as indicated
in Fig. 1. Contributions to the winding number stem from branch cuts (bands)
and from poles (edge states).

The result (41) is a multichannel generalization of the analogous invariant
derived in [pletyukhov_etal_prb_20, , Pletyukhov_etal2_2020, ] for single-
channel models. In Appendix C.1 we demonstrate the equivalence of the
representation (42) with those quoted for $I$ in the above cited papers.

Figure 1: Schematic illustration of the winding of the quantity
$\det([G_{0,0}^{(0)}]^{-1}G_{0,1}^{(0)}t_{Z})$ in the complex $\omega$ plane.
Here, the complex contour traversed while changing the frequency $\omega$ is
indicated by the blue lines, with arrows indicating the direction. In turn,
the branch cuts and bound-state poles, which give the spectral ranges of the
boundary problem, are indicted as bold red lines and dots, respectively. It is
only these non-analytic features (lying below the chemical potential $\mu$)
which contribute to the winding number integral.

##  IV Interface charge

In this section we study an excess charge, which is accumulated at the
interface of two semi-infinite models. It is generally defined by

$\displaystyle Q_{I}$
$\displaystyle=\sum_{m=-\infty}^{\infty}(\rho_{m}-\bar{\rho}^{(0)})f_{m}$ (43)
$\displaystyle=\sum_{m=-\infty}^{\infty}(\rho_{m}-\rho_{m}^{(0)})f_{m}$ (44)
$\displaystyle+\sum_{m=-\infty}^{\infty}(\rho_{m}^{(0)}-\bar{\rho}^{(0)})f_{m}.$
(45)

Here for negative $m$ the envelope function is defined as a mirror image of
its positive part, and in addition we introduce a splitting into the two sums
which is analogous to the splitting made in Eqs. (31), (32).

Starting from the translationally invariant model, we consider two types of
interfaces: 1) adding a finite potential barrier $V_{0}$ on site $m=0$ (Sec.
IV.1) and 2) weakening the link between sites $m=0$ and $m=1$ by
$0\leq\lambda\leq 1$ (Sec. IV.2). Similar models have been recently studied in
Ref. [rhim_2018, ] in the search of a unified bulk-boundary correspondence for
band insulators.

Since in the underlying model unit cells of the right and left semi-infinite
lattices perfectly match, there is no dipole moment in the overall charge
distribution, and therefore the contribution (45) identically vanishes. The
net interface charge is thus created by removing electrons from the valence
band by emerging interface localized states. It is entirely given by the
contribution (44), which is alternatively represented by

$\displaystyle Q_{I}=$ $\displaystyle-\frac{1}{\pi}\,\text{Im}\int
d\omega\Theta(\mu-\omega)$
$\displaystyle\times\sum_{m=-\infty}^{\infty}\text{tr}\left\\{G_{m,m}-G_{m,m}^{(0)}\right\\},$
(46)

where $G$ is the Green’s function of the interface model, and $G^{(0)}$
corresponds to the bulk model.

For the above physical reason $Q_{I}$ is expected to be quantized in integer
units of the electron’s charge. In the following subsections we derive winding
number expressions for $Q_{I}$, which provide a quantitative confirmation of
our expectations. In other words, we analytically prove — now in the
multichannel setting — the nearsightedness principle postulated in Ref.
[Kohn_1996, ] and used in Ref. [Pletyukhov_etal3_2020, ] to prove the
quantization of the invariant. Details of our intermediate evaluation are
provided in Appendix B.

###  IV.1 Potential impurity

In this subsection we consider a model characterized by the Green’s function
(21), that is featuring the additional impurity potential $V_{0}$ at the site
$m=0$. On the basis of (46) it holds

$\displaystyle Q_{I}=$ $\displaystyle-\frac{1}{\pi}\,\text{Im}\int
d\omega\,\,\Theta(\mu-\omega)$ $\displaystyle\times
V_{0}\,\,\text{tr}\left\\{[1-V_{0}\,G_{0,0}^{(0)}]^{-1}\sum_{m=-\infty}^{\infty}G_{0,m}^{(0)}\,G_{m,0}^{(0)}\right\\}.$
(47)

Making use of (224), we express

$\displaystyle Q_{I}=$ $\displaystyle-\frac{1}{\pi}\,\text{Im}\int
d\omega\,\,\Theta(\mu-\omega)$
$\displaystyle\times\text{tr}\left\\{(1-V_{0}F_{0})^{-1}\,\frac{\partial(1-V_{0}F_{0})}{\partial\omega}\right\\}$
(48) $\displaystyle=$ $\displaystyle-\frac{1}{\pi}\int
d\omega\,\,\Theta(\mu-\omega)$
$\displaystyle\times\,\text{Im}\,\frac{\partial}{\partial\omega}\ln\det(1-V_{0}\,G_{0,0}^{(0)}),$
(49)

where $F_{0}=G_{0,0}^{(0)}$ [cf. Eq. (96)]. Like $I$ in (42), $Q_{I}$ acquires
integer values. They are equal to winding numbers of the function
$\det(1-V_{0}\,G_{0,0}^{(0)})$ in the complex $\omega$-plane [cf. Fig. 3].

At $V_{0}=0$ we restore the translational invariance, and therefore
$Q_{I}(V_{0}=0)=0$.

In the limit $V_{0}\to\infty$ (two isolated subsystems excluding site $m=0$)
we obtain

$\displaystyle Q_{I}(V_{0}\to\infty)=$ $\displaystyle-\frac{1}{\pi}\int
d\omega\,\,\Theta(\mu-\omega)\,\text{Im}\,\frac{\partial}{\partial\omega}\ln\det
G_{0,0}^{(0)}.$ (50)

###  IV.2 Link weakening

In this subsection we generalize the result of Ref. [Pletyukhov_etal3_2020, ]
(see Appendix C therein) to the multichannel case.

We add the perturbation
$V_{\lambda}=(1-\lambda)\,\,[t_{Z}^{\dagger}\,|0\rangle\langle
1|+t_{Z}\,|1\rangle\langle 0|]$, with $0\leq\lambda\leq 1$, to the
translationally invariant model. The corresponding Dyson equation reads

$\displaystyle
G_{m,m^{\prime}}=G_{m,m^{\prime}}^{(0)}+G_{m,1}^{(0)}\,t_{\lambda}\,G_{0,m^{\prime}}+G_{m,0}^{(0)}\,t_{\lambda}^{\dagger}\,G_{1,m^{\prime}}.$
(51)

where $t_{\lambda}=(1-\lambda)t_{Z}$. It follows

$\displaystyle G_{0,m^{\prime}}$
$\displaystyle=G_{0,m^{\prime}}^{(0)}+G_{0,1}^{(0)}\,t_{\lambda}\,G_{0,m^{\prime}}+G_{0,0}^{(0)}\,t_{\lambda}^{\dagger}\,G_{1,m^{\prime}},$
(52) $\displaystyle G_{1,m^{\prime}}$
$\displaystyle=G_{1,m^{\prime}}^{(0)}+G_{1,1}^{(0)}\,t_{\lambda}\,G_{0,m^{\prime}}+G_{1,0}^{(0)}\,t_{\lambda}^{\dagger}\,G_{1,m^{\prime}}.$
(53)

Solving these equations for $G_{0,m^{\prime}}$ and $G_{1,m^{\prime}}$ (see
Appendix B) and inserting the obtained solutions into (51), we evaluate (46).
Performing the sum over all sites, we derive the winding number expression

$\displaystyle Q_{I}=$ $\displaystyle-\frac{1}{\pi}\int
d\omega\,\,\Theta(\mu-\omega)$
$\displaystyle\times\text{Im}\,\frac{\partial}{\partial\omega}\ln\det\left[(1-\lambda^{2})\left(\Lambda+\frac{1}{2}\right)+\lambda^{2}\right],$
(54)

where $\Lambda=\frac{1}{2}-F_{0}^{-1}\,G_{0,1}^{(0)}\,t_{Z}\,F_{0}$ [cf. Eq.
(184)].

At $\lambda=1$ we restore the translational invariance, and therefore
$Q_{I}(\lambda=1)=0$.

At $\lambda=0$ (two isolated subsystems including site $m=0$) we obtain

$\displaystyle Q_{I}(\lambda=0)$ $\displaystyle=-\frac{1}{\pi}\int
d\omega\,\,\Theta(\mu-\omega)$
$\displaystyle\times\text{Im}\,\frac{\partial}{\partial\omega}\,\ln\det\left(\Lambda+\frac{1}{2}\right).$
(55)

Using Eqs. (187) and (194), we deduce

$\det\left[(1-\lambda^{2})\left(\Lambda+\frac{1}{2}\right)+\lambda^{2}\right]=\det\left[1-\left(1-\lambda^{2}\right)\,G_{0,1}^{(0)}t_{Z}\right],$
(56)

which is a convenient representation for evaluating the invariant (54) at
finite $\lambda$.

In Appendix C.2 we demonstrate how to reproduce from (54) the single-channel
expression for the interface charge previously derived in Ref.
[Pletyukhov_etal3_2020, ].

##  V Numerical Results

In this section we demonstrate the validity of the above defined novel
topological invariants, particularly Eqs. (42), (49), and (54). To this end we
evaluate these winding numbers for some randomly generated multichannel models
and compare their values to the corresponding boundary and interface charges,
evaluated from their original definitions. As it turns out, even for the
minimal nontrivial (that is, allowing for the ”spin-orbit” coupling)
multichannel models ($N_{c}=Z=2$), the analytical evaluation of these winding
numbers is unfeasible, so we have to resort to a numerical calculation.

This section is organzied as follows: in Sec. V.1 we describe our strategy of
the numerical evaluation of the winding numbers. In particular, we describe an
algorithm which allows for an efficient computation of these invariants. In
Sections V.2 and V.3 we show representative data for the winding numbers
corresponding to the boundary and interface charges, respectively.

###  V.1 Numerical implementation of the winding numbers

All three winding numbers discussed in this work are defined in terms of bulk
Green’s functions, which in turn are defined via a quasimomentum integral, see
Eq. (15). In the numerical evaluation, we approximate this integral by a sum
over $N_{k}$ evenly spaced (by $\delta_{k}=\frac{2\pi}{N_{k}}$) momenta
$k_{l}=-\pi+l\delta_{k}$, with $l=0,\dots,N_{k}-1$, i.e.

$G^{(0)}_{m\sigma,m^{\prime}\sigma^{\prime}}=\lim_{N_{k}\to\infty}\frac{\delta_{k}}{2\pi}\sum_{k_{l}}\left<j,\sigma\left|\frac{e^{ik_{l}(n-n^{\prime})}}{\omega+i\eta-
h_{k_{l}}}\right|j^{\prime},\sigma^{\prime}\right>.$ (57)

The various winding numbers computed from these Green’s functions are then
compared to exact diagonalization data computed from some finite size
Hamiltonian, suitably defined according to the boundary configuration in
question, i.e. hard wall, potential impurity or link weakening [see Fig. 2].
These finite size Hamiltonians have some dimension $D=N_{\text{sites}}N_{c}$,
where $N_{\text{sites}}$ denotes the number of sites (chosen in the following
to accommodate an integer number $N_{Z}$ of unit cells, that is
$N_{\text{sites}}=ZN_{Z}$), and are characterized by some spectrum
$\epsilon_{s}$ (where $s=1,\dots,D$) and bandwidth
$B=\max_{s}\epsilon_{s}-\min_{s}\epsilon_{s}$. We define the average level
spacing of a finite size Hamiltonian as

$\delta_{\epsilon}=\frac{B}{D}.$ (58)

The Green’s function evaluated via the sum over momenta converges to that of
an infinite system, if

$\sup_{k}||\delta_{k}\partial_{k}h_{k}||\ll\delta_{\epsilon}.$ (59)

Here, $||\dots||$ denotes any matrix norm. In the following we will thus
assume that all Green’s functions are computed according to Eq. (57) where the
momentum spacing $\delta_{k}$ is not sent to zero, but kept finite and chosen
according to Eq. (59). In practice we choose $\delta_{k}=\frac{1}{D}$ and find
that with this choice Eq. (59) is always fulfilled.

Figure 2: Diagram depicting the various boundary configurations considered
here: hard wall (a), potential impurity (b) and weakened link (c). In all
three cases we consider a closed ring of $N_{Z}$ unit cells, as indicated by
the bent dotted lines underneath the diagrams which contain $N_{Z}-2$
additional unit cells respectively. The labels above the diagrams correspond
to the numbering of the unit cells as discussed in Sec. II. The feature
distinguishing a given boundary configuration from the translationally
invariant ring is highlighted in red, respectively. a) shows the boundary
configuration of a hard wall, which is used in the computation of the shift in
boundary charge. The link between unit cells $n=0$ and $n=1$ is cut completely
and the boundary charge is measured to the right, i.e. starting with unit cell
$n=1$. The envelope function falls off to zero to the right of $n=1$, around
$n\sim N_{Z}/4$. b) shows the potential impurity boundary configuration, where
the onsite potential $V_{0}$ is added onto site $m=0$ ($n=0,j=Z$) of the
translationally invariant ring. Note that in the limit $V_{0}\to\infty$ we
effectively obtain a hard wall again (as in a)), but with the site $m=0$
missing. This boundary configuration is used to evaluate the interface charge,
that is the envelope function extends across the impurity to the left and
right, over the range between $-N_{Z}/4$ and $N_{Z}/4$. c) shows the boundary
configuration of a weakened link, where the hopping matrices between unit
cells $n=0$ and $n=1$ are reduced by a factor $0\leq\lambda\leq 1$. For
$\lambda=1$ the translationally invariant ring is recovered, the limit
$\lambda=0$ corresponds to the hard wall configuration of a). The interface
charge is again computed using an envelope function extending across the
impurity to the left and right.

Having settled the question of how to compute the bulk Green’s functions we
now turn to the question of how to evaluate the winding numbers themselves.
All three winding numbers defined in this work are given in terms of energy
integrals of the following kind (up to overall minus signs)

$\text{wn}\left[K\right]=\int\frac{d\omega}{\pi}\Theta(\mu-\omega)\text{Im}\,\partial_{\omega}\ln
K(\omega),$ (60)

where

$K\in\begin{cases}\det\left([G^{(0)}_{0,0}]^{-1}G_{0,1}^{(0)}t_{Z}\right)\\\
\det\left(1-V_{0}G_{0,0}^{(0)}\right)\\\
\det\left[1-\left(1-\lambda^{2}\right)G_{0,1}^{(0)}t_{Z}\right]\end{cases},$
(61)

depending on the quantity in question. Note that in all three cases
$K^{\ast}=K\rvert_{\omega+i\eta\to\omega-i\eta}$. The winding number
expression can thus be recast as follows:

$\displaystyle\text{wn}\left[K\right]$ $\displaystyle=\frac{1}{2\pi
i}\left[\int\displaylimits_{-\infty}^{\mu}\partial_{\omega}\ln
K\rvert_{\omega+i\eta}+\int\displaylimits_{\mu}^{-\infty}\partial_{\omega}\ln
K\rvert_{\omega-i\eta}\right]$ (62)
$\displaystyle=\oint_{\mathcal{C}}\frac{d\omega}{2\pi i}\partial_{\omega}\ln
K(\omega)=\lim_{n_{c}\to\infty}\sum_{n=1}^{n_{c}}\Delta_{n},$ (63)

where $\mathcal{C}$ denotes a rectangular contour with infinitesimal width
$2\eta$ in the imaginary direction, ranging (in the real direction) from the
lowest band edge $\min_{s}\epsilon_{s}$ (instead of $-\infty$) to the chemical
potential $\mu$, and where

$\Delta_{n}=\begin{cases}\frac{\arg K(\omega_{n+1})-\arg
K(\omega_{n})}{2\pi}\quad\text{if}\quad n+1\leq n_{c}\\\ \frac{\arg
K(\omega_{1})-\arg K(\omega_{n_{c}})}{2\pi}\quad\text{if}\quad
n=n_{c}\end{cases},$ (64)

with the $n_{c}$ ordered samples
$\left\\{\omega_{n}\right\\}_{n=1,\dots,n_{c}}$ from the contour $\mathcal{C}$
(see panel a) of Fig. 3).

Note that the value of the integral does not change upon a continuous
deformation of the rectangular contour $\mathcal{C}$, as long as one does not
cross any non-analytic features of $K$ (e.g. poles or branch-cuts) in the
process. Due to the symmetry under complex conjugation mentioned above it is
however numerically advantageous to resort to a rectangular contour, since any
evaluation of $K$ automatically yields two samples needed for the computation
of the winding number (at $\omega$ and at $\omega^{\ast}$). In the following
we will thus assume that we always exploit this symmetry.

Since the evaluation of $K$ for a single sample $\omega_{n}$ involves the
computationally expensive sum over $N_{k}$ momenta mentioned above, it is
desirable to limit the number of samples $n_{c}$ to a minimum. We achieve this
by not using a fixed grid of evenly spaced samples along the energy contour,
but by starting with a small number of evenly spaced samples (at least three),
and then iteratively adding samples in between those points, for which the
absolute value of the phase difference $|\Delta_{n}|$ is largest. The
algorithm consists of the following steps:

  * •

Choose $N_{\text{start}}$ evenly spaced samples $\omega_{n}$ on the energy
contour $\mathcal{C}$, compute the corresponding values of $K(\omega_{n})$,
and with those the phase differences $\Delta_{n}$. At each stage of the
algorithm, the sum $\sum_{n=1}^{n_{c}}\Delta_{n}$ provides an estimate of the
value of the winding number and $\max_{n}|\Delta_{n}|$ provides a measure for
the error in the winding number, which is expected to behave like
$\sim\mathcal{O}(n_{c}^{-1})$.

  * •

As long as there are neighbouring points on the $K$-contour, for which the
absolute value of the phase difference $|\Delta_{n}|$ is larger than some
threshold $\delta$: Determine $n_{0}$, such that
$|\Delta_{n_{0}}|=\max_{n}|\Delta_{n}|$, evaluate $K$ at the value
$\omega_{\text{new}}$ on the energy contour which lies in between (along the
energy contour) $\omega_{n_{0}}$ and $\omega_{n_{0}+1}$. Compute the phase
differences for the two new pairs of neighbours
($\left\\{\omega_{n_{0}},\omega_{\text{new}}\right\\}$ and
$\left\\{\omega_{\text{new}},\omega_{n_{0}+1}\right\\}$) which emerge due to
the addition of this new sample, recompute the sum over all phase differences.

Once the absolute values of all phase differences are below a certain
threshold, the (rounded) final estimate $\sum_{n=1}^{n_{c}}\Delta_{n}$ is
taken as the result of the winding number computation. Note that whenever we
evaluate $K$ at a given $\omega$, we also get the value at $\omega^{\ast}$ for
free, by exploiting the above-mentioned symmetry.

Fig. 3 shows an example of the computation of a winding number for
$K=\det\left([G^{(0)}_{0,0}]^{-1}G_{0,1}^{(0)}t_{Z}\right)$, $N_{c}=Z=3$,
three occupied bands and a randomly generated Hamiltonian whose matrix
elements are listed in the Supplemental Material [supp_mat, ].

Panel a) of Fig. 3 shows the energy contour with samples $\omega_{n}$ (in
orange) and the extent of the bands as indicated by the colored regions (white
areas correspond to gaps, and the red star denotes the pole). In practice, the
extent $\eta$ of the contour in the imaginary direction is not infinitesimally
small, but rather finite and small compared to the bandwidth. In this plot,
and in what follows, we always use $\eta=10\delta_{\epsilon}$. One can clearly
see that the density of samples along the contour varies greatly, with densely
populated, but also completely blank sections. Together with the quick
convergence of the winding number discussed below, this demonstrates the
usefulness of the above described algorithm, as compared to a fixed, uniform
grid, which would lead to many unnecessary evaluations of $K$.

Panels b) and c) show the estimate of the winding number
$\sum_{n=1}^{n_{c}}\Delta_{n}$ and the scaled error
$n_{c}\max_{n}|\Delta_{n}|$ as functions of the number of samples $n_{c}$,
respectively. It is apparent that convergence is reached quickly and that the
error falls off as $\sim\mathcal{O}(n_{c}^{-1})$, as expected. Occasional
peaks in the scaled error correspond to the discovery of new parts of the
$K$-contour by the algorithm. Once such a new section has been discovered it
is quickly smoothed out by considering more samples in the corresponding
region of the energy contour, reducing the scaled error again.

Panel d) shows the contour described by $K(\omega_{n})$ in the limit where the
value of the winding number has converged. The colors of the various sections
of the contour correspond to the equally colored bands in panel a), while
sections corresponding to gaps are colored in gray.

In this example (and in all following calculations), we use
$N_{\text{start}}=10$ and $\delta=\frac{1}{360}$, i.e. in the converged limit
no pair of neighbouring points has an absolute phase difference of more than
one degree as measured from the origin. Although we have discussed the special
case of $K=\det\left([G^{(0)}_{0,0}]^{-1}G_{0,1}^{(0)}t_{Z}\right)$ here, we
report that other combinations of Green’s functions, and in particular the
relevant functions defined in Eq. (61), show qualitatively similar behavior as
the example shown in Fig. 3.

Figure 3: Representative data of the evaluation of a
$K=\det\left([G^{(0)}_{0,0}]^{-1}G_{0,1}^{(0)}t_{Z}\right)$ winding number for
a randomly generated $N_{c}=Z=3$ model with three occupied bands (components
of the Hamiltonian are listed in the Supplemental Material [supp_mat, ]).
Panel a) shows samples comprising the energy contour (in orange) and the
distribution of bands (colored) and gaps (white) and an edge state pole marked
by the red star. The solid black rectangle underneath the samples denotes the
complete contour, the curved arrow in the upper right indicates its
orientation. Panels b) and c) show the estimate of the winding number
$\sum_{n=1}^{n_{c}}\Delta_{n}$ and the scaled error
$n_{c}\max_{n}|\Delta_{n}|$ as functions of the number of samples $n_{c}$,
respectively. Panel d) shows the $K-$contour with arrows indicating the
orientation. It is characterized by a winding number equal to two, which
accords with the value obtained in panel b). The colors of the various
sections of the $K-$contour correspond to the bands, while gaps are denoted by
gray color.

###  V.2 Numerical validation of the boundary charge invariant

The boundary charge invariant [Eq. (42)] is related to the change in boundary
charge under a shift of the lattice by one site towards the boundary. In order
to compute the boundary charge we define a finite size Hamiltonian with
$N_{Z}$ unit cells and cut the system, as depicted in Fig. 2(a).

In order to evaluate the change in the boundary charge upon continuous shift
of the system towards the left boundary we consider the following form for the
onsite potentials $v_{n}$ and hopping matrices $t_{n}$:

$\displaystyle v_{n}$
$\displaystyle=v_{0}+\cos\left(\varphi+\frac{2\pi(n-1)}{Z}\right)\delta_{v},$
(65) $\displaystyle t_{n}$
$\displaystyle=t_{0}+\cos\left(\varphi+\frac{2\pi(n-1)}{Z}\right)\delta_{t},$
(66)

where $v_{0}=v_{0}^{\dagger},\delta_{v}=\delta_{v}^{\dagger},t_{0}$ and
$\delta_{t}$ are random $N_{c}\times N_{c}$-dimensional matrices. Here, the
phase variable $\varphi$ controls the shift of the right subsystem toward the
boundary, in particular a shift of $\varphi\to\varphi+\frac{2\pi}{Z}$ shifts
it by one site.

Note that in the previous work [pletyukhov_etal_prb_20, ] a more general class
of phase dependencies for the components of the Hamiltonian was considered,
with higher Fourier coefficients and additional random phase shifts.
Restricting ourselves to the form stated above however suffices to demonstrate
the validity of the novel winding number invariant defined in this work, since
no conceptual differences arise in the spectral properties of the Hamiltonian
upon including such higher Fourier coefficients.

We diagonalize the finite size Hamiltonian with open boundary conditions and
$N_{Z}=100$ unit cells in order to compute the boundary charge
$Q_{B}(\varphi)$ in the right subsystem (i.e. between the cut link on the left
$n=1$ and the envelope’s function fall-off region $n\sim N_{Z}/4$ on the
right) for a given chemical potential $\mu$ and phase $\varphi$. The change in
boundary charge upon shift by one site is then simply given by $\Delta
Q_{B}(\varphi)=\tilde{Q}_{B}(\varphi)-Q_{B}(\varphi)$ where
$\tilde{Q}_{B}(\varphi)=Q_{B}\left(\varphi+\frac{2\pi}{Z}\right)$.

Fig. 4 shows a comparison of the boundary charge $Q_{B}$ (computed in the way
stated above), the change in boundary charge $\Delta Q_{B}$ and the winding
number invariant $I$, all as functions of the phase $\varphi$ for a randomly
generated $N_{c}=3,Z=4$ model with two occupied bands. The matrices
$v_{0},\delta_{v},t_{0},\delta_{t}$ used in this example are listed in the
Supplemental Material [supp_mat, ].

Panel a) of Fig. 4 shows a band structure of the translationally invariant
ring as a function of the phase $\varphi$. Some of the $N_{c}Z=12$ bands touch
at special phases, but the gap above the second band remains open throughout
the whole phase cycle. The dashed line at the top of the second band
corresponds to the chemical potential and the rectangular box denotes the
relevant region of the band structure.

Panel b) shows a zoom on the relevant bands and gaps (especially the second
gap), the chemical potential as a dashed black line and the energies of edge
states localized on the right from the cut link in both the initial (orange)
and shifted (green) finite system with $N_{Z}=100$ unit cells, all as
functions of the phase $\varphi$.

Panel c) shows the boundary charge $Q_{B}$ and the boundary charge of the
shifted system $\tilde{Q}_{B}$ as a function of the phase $\varphi$. One can
see that the boundary charge behaves almost linearly $\sim\frac{\varphi}{\pi}$
interrupted by two jumps by $-1$ at those phases, where edge states leave the
topmost occupied band.

Panel d) shows the difference in boundary charge $\Delta Q_{B}$ and the
winding number $I$ [Eq. (42)], both as functions of the phase $\varphi$. One
can see that the two lineshapes are identical, up to an offset of
$-\frac{\nu}{Z}=-\frac{1}{2}$, demonstrating the validity of the boundary
charge invariant for this randomly generated model.

We report that in preparation of this work we have simulated numerous of these
random models with varying $N_{c},Z$ and number of occupied bands and have
never seen failure of Eq. (42). We furthermore report that in no case we have
seen $I<-N_{c}$ or $I>0$.

Figure 4: Representative data of the comparison of boundary charge $Q_{B}$ and
boundary charge invariant $I$ for a randomly generated $N_{c}=3,Z=4$ model
with two occupied bands (components of the Hamiltonian are listed in the
Supplemental Material [supp_mat, ]). Panel a) shows the extent of the bands
and gaps of the model as a function of the phase $\varphi$, the chemical
potential (dashed line) and a rectangular box indicating the relevant energy
window. Panel b) shows a zoom of this region with additional orange (green)
lines indicating edge states localized on the right from the cut link in the
system (the shifted system), as computed from the diagonalization of a finite
size Hamiltonian with $N_{Z}=100$ unit cells respectively. The chemical
potential is indicated as a dashed line. Panel c) shows both the boundary
charge $Q_{B}$ and shifted boundary charge
$\tilde{Q}_{B}=Q_{B}\left(\varphi+\frac{2\pi}{Z}\right)$ as functions of
$\varphi$. Panel d) shows the change in boundary charge $\Delta Q_{B}$ and the
winding number invariant $I$ [Eq. (42)]. The offset between these two curves
is exactly $-\frac{\nu}{Z}=-\frac{1}{2}$ for all values of $\varphi$.

###  V.3 Numerical validation of the interface charge invariants

In addition to the winding number relating to the change in boundary charge
discussed above, this work also introduces two winding numbers relating to the
interface charge that accumulates around an impurity in an otherwise
translationally invariant system. We consider two types of impurities: a local
potential on site $m=0$ and a weakened link between sites $m=0$ and $m=1$.

In the following we discuss one example for each of these two cases
respectively, comparing the value of the winding numbers to the interface
charges both as functions of the strength of the impurity. To this end we
model the translationally invariant system as a ring with $N_{Z}$ unit cells
and add the respective impurity [see Fig. 2(b,c)]. We then diagonalize these
finite size Hamiltonians and compute the sum of the boundary charges to the
left and to the right of the impurity. We denote these two types of boundary
charges by $Q_{B}^{(R)}$ (to the right of the impurity) and $Q_{B}^{(L)}$ (to
the left of the impurity, also including the impurity site $m=0$), such that
the interface charge is given by $Q_{I}^{\text{def}}=Q_{B}^{(R)}+Q_{B}^{(L)}$,
where we use the label ”def” in order to distinguish the interface charge
calculated on the basis of the definition (43) from the winding number
invariants, which we label by $Q_{I}$ [see Eqs. (49), (54)].

####  V.3.1 Potential impurity

We consider a finite size Hamiltonain with $N_{Z}=100$ unit cells in a ring
configuration and add the potential impurity $V_{0}$ on site $m=0$, as
depicted in Fig. 2(b). We diagonalize this Hamiltonian and use the resulting
wavefunctions to compute the interface charge.

Fig. 5 shows a comparison of the such computed interface charge (which we
denote by $Q_{I}^{\text{def}}$) and the corresponding topological invariant
$Q_{I}$ (as defined in Eq. (49)) for a randomly generated $N_{c}=4,Z=3$ model
with four occupied bands. The Hamiltonian blocks are listed in the
Supplemental Material [supp_mat, ].

Panel Fig. 5(a) shows the band structure of the translationally invariant
model (i.e. for $V_{0}=0$) versus $k$, with a rectangle indicating the
relevant region around the fourth gap.

Panel b) shows a zoom onto this region with additional lines indicating the
energies of bound states as functions of the impurity strength $V_{0}$ in
units of the bandwidth $B$. One can see that as the impurity strength is
increased, in total three bound states emerge from the fourth band, two of
which remain in the gap for all values of $V_{0}$, while the third bound state
joins the overlying band.

Panel c) shows the interface charge $Q_{I}^{\text{def}}$ and the invariant
$Q_{I}$ [Eq. (49)], both as functions of the impurity strength. The lineshapes
of $Q_{I}$ and $Q_{I}^{\text{def}}$ perfectly overlap, demonstrating the
validity of this topological invariant.

Again, we report that in preparation of this work we have tested many random
models, always finding agreement of (43) with Eq. (49).

Figure 5: Representative data for the comparison of the interface charge
computed on the basis of Eq. (43) and the interface charge invariant (49) for
a randomly generated $N_{c}=4,Z=3$ model with a potential impurity of strength
$V_{0}$ on site $m=0$ and four occupied bands. The components of the
Hamiltonian are listed in the Supplemental Material [supp_mat, ]. Panel a)
shows the band structure of the translationally invariant system ($V_{0}=0$)
as a function of quasimomentum $k$ and a rectangle indicating the relevant
energy window above the fourth band. Panel b) shows the gap above the fourth
band and the energies of bound states in a $N_{Z}=100$ model as a function of
the impurity strength $V_{0}$ measured in units of the bandwidth $B$. Panel c)
shows the lineshapes of $Q_{I}^{\text{def}}$ and $Q_{I}$, both as functions of
$V_{0}$, which perfectly overlap.

####  V.3.2 Link weakening

Next, we consider a finite size Hamiltonain with $N_{Z}=100$ unit cells in a
ring configuration and weaken the link between sites $m=0$ and $m=1$ by
$0\leq\lambda\leq 1$, as depicted in Fig. 2(c). We diagonalize this
Hamiltonian and use the resulting wavefunctions to compute the interface
charge.

Fig. 6 shows a comparison of the interface charge $Q_{I}^{\text{def}}$
(computed in terms of (43) as described above) with the corresponding winding
number invariant $Q_{I}$ given in Eq. (54) for a randomly generated
$N_{c}=Z=5$ model with six occupied bands (the Hamiltonian parameters are
listed in the Supplemental Material [supp_mat, ]).

Panel a) of Fig. (6) shows the band structure of the translationally invariant
system (i.e. for $\lambda=1$) as a function of quasimomentum $k$ and the
relevant gap as denoted by the black rectangle.

Panel b) shows a zoom onto the relevant gap with the chemical potential
(dashed black line) and the energy of eigenstates of the finite system with
$N_{Z}=100$ unit cells as a function of the link strength $\lambda$. One can
see that in total, as the strength of the link is weakened, two eigenmodes
leave the sixth band of the system turning into bound states.

Panel c) shows the boundary charges $Q_{B}^{(L)}$ and $Q_{B}^{(R)}$ on the
left ($m\leq 0$) and right ($m\geq 1$) sides of the impurity, the interface
charge $Q_{I}^{\text{def}}$ computed on the basis of (43), and the winding
number invariant $Q_{I}$, all as functions of the link strength $\lambda$. One
can see that $Q_{I}^{\text{def}}$ and $Q_{I}$ agree throughout the range of
$\lambda$, confirming the validity of Eq. (54) for this random model.

Figure 6: Representative data for the comparison of the interface charge
computed on the basis of Eq. (43) and the interface charge invariant (54) for
a randomly generated $N_{c}=Z=5$ model with a weakened link between sites
$m=0$ and $m=1$ with six occupied bands. The components of the Hamiltonian are
listed in the Supplemental Material [supp_mat, ]. Panel a) shows the band
structure of the translationally invariant system (i.e. for $\lambda=1$) as a
function of quasimomentum $k$ with the relevant gap marked by the black
rectangle. Panel b) shows a zoom into this gap with dots depicting the energy
of localized eigenstates of the finite Hamiltonian with $N_{Z}=100$ unit cells
as function of the link strength $\lambda$. The dashed line indicates the
value of the chemical potential. Panel c) shows the lineshapes of the boundary
charges on the left and right side of the impurity ($Q_{B}^{(L/R)}$), the
interface charge $Q_{I}^{\text{def}}$ and the invariant $Q_{I}$, all as
functions of the link strength $\lambda$. Note that $Q_{I}$ and
$Q_{I}^{\text{def}}$ perfectly overlap.

Again, we report that we tested numerous random models and always found
agreement of (43) with Eq. (54).

##  VI Summary

For generic one-dimensional insulators with $Z$ sites per unit cell and
$N_{c}$ channels (orbitals) per site we have developed a general theoretical
framework for boundary and interface charge investigation, employing the
method of boundary Green’s functions. Using this approach, we represent the
boundary and interface charges in terms of bulk Green’s functions of the
system, which proves particularly advantageous when dealing with multi-channel
systems, in which the construction of exact eigenstates of the system with
broken translational invariance is usually a rarely achievable goal due to the
possibility of overlapping bands.

We use this representation of the boundary charge in order to establish the
topological invariant associated with the changes of boundary charge under
lattice shifts. It was shown that upon a shift of the lattice as a whole by a
single lattice site $Q_{B}$ (relative to the average density in the bulk
$\bar{\rho}$) changes by an integer $-N_{c}\leq I\leq 0$, which can be
expressed as a winding number of a particular combination of bulk Green’s
function components. This result is a natural generalization of our previous
findings regarding single-channel systems, where $I$ was shown to take on
values in $0$ or $-1$. The present finding proves the conjecture in Sec. VI A
of Ref. [pletyukhov_etal_prb_20, ]. The character of the bulk-boundary
correspondence described by this invariant is not limited to any particular
symmetry class and holds for arbitrary systems. This invariant quantifies the
spectral flow of the boundary eigenvalue problem, the change by $-1$ in it
indicating every time that an electron is taken away by an emerging edge state
(which can happen at most $N_{c}$ times during one pumping cycle of the phase
$\varphi$).

In addition to the quantization of $\Delta Q_{B}-\bar{\rho}$ we established
the quantization of charge accumulated on local scattering centers, also known
as interface charge. Specifically, we demonstrate that the charge accumulating
either around a single potential impurity or a weak link is given by a
corresponding winding number and is thus a topological invariant. This
topological invariant also quantifies the spectral flow — now of the interface
eigenvalue problems. This observation generalizes the findings of Ref.
[Pletyukhov_etal3_2020, ] (see Appendix C therein) and provides an analytical
proof of the nearsightedness principle for generic non-interacting one-
dimensional tight-binding models with translational invariance to which a
single impurity is added.

In future work [Pias_new, ] it would be also interesting to extend the
obtained results to a class of models with longer than nearest-neighbor range
hoppings. Although longer ranged hoppings can always be rewritten effectively
in terms of nearest-neighbor hopping by increasing the number of channels per
site, it is of interest to avoid this artificial procedure and to find winding
number expressions for the invariants including longer ranged hoppings.

##  VII Acknowledgments

We acknowledge fruitful discussions with J. Klinovaja and D. Loss on the
subject of the boundary charge. The work was supported by the Deutsche
Forschungsgemeinschaft via RTG 1995.

##  Appendix A Properties of boundary Green’s functions and boundary charge

###  A.1 Boundary charge expression

Let us introduce the following short-hand notation for the Bloch Hamiltonian
defined in Eq. (10)

$\displaystyle h_{k}$ $\displaystyle=\left(\begin{array}[]{cc}A&b_{k}\\\
b_{k}^{\dagger}&v_{Z}\end{array}\right),$ (69) $\displaystyle
h_{k}-(\omega+i\eta)$ $\displaystyle=\left(\begin{array}[]{cc}\bar{A}&b_{k}\\\
b_{k}^{\dagger}&\bar{v}_{Z}\end{array}\right),$ (72)

where $\bar{A}=A-(\omega+i\eta)$, $\bar{v}_{Z}=v_{Z}-(\omega+i\eta)$,

$\displaystyle b_{k}=-\left(\begin{array}[]{c}t_{Z}e^{-ik}\\\ 0\\\ \vdots\\\
0\\\ t_{Z-1}^{\dagger}\end{array}\right),\quad
b_{k}^{\dagger}=-\left(e^{ik}t_{Z}^{\dagger},0,\dots,0,t_{Z-1}\right)$ (78)

and

$A=\begin{pmatrix}v_{1}&-t_{1}^{\dagger}\\\ -t_{1}&v_{2}&\ddots\\\
&\ddots&\ddots&-t_{Z-2}^{\dagger}\\\ &&-t_{Z-2}&v_{Z-1}\end{pmatrix}.$ (79)

Using the matrix identity

$\displaystyle\left(\begin{array}[]{cc}M_{11}&M_{12}\\\
M_{21}&M_{22}\end{array}\right)^{-1}=\left(\begin{array}[]{cc}(M_{11}-M_{12}M_{22}^{-1}M_{21})^{-1}&-M_{11}^{-1}M_{12}(M_{22}-M_{21}M_{11}^{-1}M_{12})^{-1}\\\
-(M_{22}-M_{21}M_{11}^{-1}M_{12})^{-1}M_{21}M_{11}^{-1}&(M_{22}-M_{21}M_{11}^{-1}M_{12})^{-1}\end{array}\right),$
(84)

where $M_{11}$ and $M_{22}$ are square invertible matrices, which may
eventually have different sizes, we evaluate the inverse of (72) and find

$\displaystyle\langle Z|\frac{1}{\omega+i\eta-h_{k}}|Z\rangle$
$\displaystyle=-(\bar{v}_{Z}-b_{k}^{\dagger}\bar{A}^{-1}b_{k})^{-1},$ (85)
$\displaystyle\langle Z|\frac{1}{\omega+i\eta-h_{k}}|j\rangle$
$\displaystyle=(\bar{v}_{Z}-b_{k}^{\dagger}\bar{A}^{-1}b_{k})^{-1}(b_{k}^{\dagger}\bar{A}^{-1})_{j},$
(86) $\displaystyle\langle j|\frac{1}{\omega+i\eta-h_{k}}|Z\rangle$
$\displaystyle=(\bar{A}^{-1}b_{k})_{j}(\bar{v}_{Z}-b_{k}^{\dagger}\bar{A}^{-1}b_{k})^{-1},$
(87)

where the last two relations are valid for $1\leq j\leq Z-1$. Using (78), we
express

$\displaystyle b_{k}^{\dagger}\bar{A}^{-1}b_{k}$
$\displaystyle=t_{Z}^{\dagger}\bar{A}^{-1}_{1,1}t_{Z}+t_{Z-1}\bar{A}^{-1}_{Z-1,Z-1}t_{Z-1}^{\dagger}$
$\displaystyle+e^{ik}t_{Z}^{\dagger}\bar{A}^{-1}_{1,Z-1}t_{Z-1}^{\dagger}+t_{Z-1}\bar{A}^{-1}_{Z-1,1}t_{Z}e^{-ik},$
(88) $\displaystyle(b_{k}^{\dagger}\bar{A}^{-1})_{j}$
$\displaystyle=-e^{ik}t_{Z}^{\dagger}\bar{A}^{-1}_{1,j}-t_{Z-1}\bar{A}^{-1}_{Z-1,j},$
(89) $\displaystyle(\bar{A}^{-1}b_{k})_{j}$
$\displaystyle=-\bar{A}^{-1}_{j,1}t_{Z}e^{-ik}-\bar{A}^{-1}_{j,Z-1}t_{Z-1}^{\dagger}.$
(90)

Introducing the notations

$\displaystyle c$
$\displaystyle=t_{Z}^{\dagger}\bar{A}^{-1}_{1,Z-1}t_{Z-1}^{\dagger},$ (91)
$\displaystyle c^{\dagger}$ $\displaystyle=t_{Z-1}\bar{A}^{-1}_{Z-1,1}t_{Z},$
(92) $\displaystyle r$
$\displaystyle=\bar{v}_{Z}-t_{Z}^{\dagger}\bar{A}^{-1}_{1,1}t_{Z}-t_{Z-1}\bar{A}^{-1}_{Z-1,Z-1}t_{Z-1}^{\dagger},$
(93) $\displaystyle m(k)$
$\displaystyle=\left(ce^{ik}+c^{\dagger}e^{-ik}-r\right)^{-1},$ (94)
$\displaystyle F_{n}$
$\displaystyle=\int_{-\pi}^{\pi}\frac{dk}{2\pi}m(k)e^{ikn},$ (95)

we establish the following relations for $j=Z$

$\displaystyle G_{0,0}^{(0)}(\omega)$ $\displaystyle=F_{0},$ (96)
$\displaystyle G_{0,m}^{(0)}(\omega)$ $\displaystyle=F_{-n},\quad
G_{m,0}^{(0)}(\omega)=F_{n},$ (97)

as well as for $j\neq Z$

$\displaystyle G_{0,m}^{(0)}(\omega)$
$\displaystyle=F_{-n+1}t_{Z}^{\dagger}\bar{A}^{-1}_{1,j}+F_{-n}t_{Z-1}\bar{A}^{-1}_{Z-1,j},$
(98) $\displaystyle G_{m,0}^{(0)}(\omega)$
$\displaystyle=\bar{A}^{-1}_{j,1}t_{Z}F_{n-1}+\bar{A}^{-1}_{j,Z-1}t_{Z-1}^{\dagger}F_{n}.$
(99)

Note that defining the hermitian conjugation in (92) we assume that it does
not conjugate the spectral parameter $\omega+i\eta$. We also remark that from
(98) it follows for $m=1$, i.e. for $n=1$ and $j=1$,

$\displaystyle
G_{0,1}^{(0)}(\omega)=(F_{0}t_{Z}^{\dagger}\bar{A}^{-1}_{1,1}t_{Z}+F_{-1}c^{\dagger})t_{Z}^{-1}.$
(100)

Using this representation we perform the sum

$\displaystyle\sum_{m=1}^{\infty}G_{0,m}^{(0)}(\omega)G_{m,0}^{(0)}(\omega)=\sum_{n=1}^{\infty}\sum_{j=1}^{Z}G_{0,m}^{(0)}(\omega)G_{m,0}^{(0)}(\omega)$
$\displaystyle=\sum_{n=1}^{\infty}(F_{-n+1}t_{Z}^{\dagger}\bar{A}^{-2}_{1,1}t_{Z}F_{n-1}+F_{-n}t_{Z-1}\bar{A}_{Z-1,Z-1}^{-2}t_{Z-1}^{\dagger}F_{n})$
$\displaystyle+\sum_{n=1}^{\infty}(F_{-n+1}t_{Z}^{\dagger}\bar{A}^{-2}_{1,Z-1}t_{Z-1}^{\dagger}F_{n}+F_{-n}t_{Z-1}\bar{A}^{-2}_{Z-1,1}t_{Z}F_{n-1})$
$\displaystyle+\sum_{n=1}^{\infty}F_{-n}F_{n}.$ (101)

With help of the differential identity

$\displaystyle\frac{\partial\bar{A}^{-1}}{\partial\omega}=-\bar{A}^{-1}\frac{\partial\bar{A}}{\partial\omega}\bar{A}^{-1}=\bar{A}^{-2}$
(102)

we simplify Eq. (101) to

$\displaystyle\sum_{m=1}^{\infty}G_{0,m}^{(0)}(\omega)G_{m,0}^{(0)}(\omega)$
$\displaystyle=$ $\displaystyle\sum_{n=1}^{\infty}(F_{-n+1}\frac{\partial
c}{\partial\omega}F_{n}+F_{-n}\frac{\partial
c^{\dagger}}{\partial\omega}F_{n-1}-F_{-n}\frac{\partial
r}{\partial\omega}F_{n})$
$\displaystyle+F_{0}t_{Z}^{\dagger}\frac{\partial\bar{A}^{-1}_{1,1}}{\partial\omega}t_{Z}F_{0}.$
(103)

Using the identities

$\displaystyle F_{-n+1}c+F_{-n-1}c^{\dagger}-F_{-n}r$
$\displaystyle=\delta_{n,0},$ (104) $\displaystyle
cF_{n+1}+c^{\dagger}F_{n-1}-rF_{n}$ $\displaystyle=\delta_{n,0},$ (105)

we find two equivalent representations

$\displaystyle\sum_{m=1}^{\infty}G_{0,m}^{(0)}(\omega)G_{m,0}^{(0)}(\omega)$
$\displaystyle=$ $\displaystyle-\sum_{n=1}^{\infty}(\frac{\partial
F_{-n+1}}{\partial\omega}c+\frac{\partial
F_{-n-1}}{\partial\omega}c^{\dagger}-\frac{\partial
F_{-n}}{\partial\omega}r)F_{n}$
$\displaystyle+F_{0}t_{Z}^{\dagger}\frac{\partial\bar{A}^{-1}_{1,1}}{\partial\omega}t_{Z}F_{0}+F_{-1}\frac{\partial
c^{\dagger}}{\partial\omega}F_{0}$ $\displaystyle=$ $\displaystyle
F_{0}t_{Z}^{\dagger}\frac{\partial\bar{A}^{-1}_{1,1}}{\partial\omega}t_{Z}F_{0}+\frac{\partial(F_{-1}c^{\dagger})}{\partial\omega}F_{0}-\frac{\partial
F_{0}}{\partial\omega}cF_{1}$ (106)

and

$\displaystyle\sum_{m=1}^{\infty}G_{0,m}^{(0)}(\omega)G_{m,0}^{(0)}(\omega)$
$\displaystyle=$ $\displaystyle-\sum_{n=1}^{\infty}F_{-n}(c\frac{\partial
F_{n+1}}{\partial\omega}+c^{\dagger}\frac{\partial
F_{n-1}}{\partial\omega}-r\frac{\partial F_{n}}{\partial\omega})$
$\displaystyle+F_{0}t_{Z}^{\dagger}\frac{\partial\bar{A}^{-1}_{1,1}}{\partial\omega}t_{Z}F_{0}+F_{0}\frac{\partial
c}{\partial\omega}F_{1}$ $\displaystyle=$ $\displaystyle
F_{0}t_{Z}^{\dagger}\frac{\partial\bar{A}^{-1}_{1,1}}{\partial\omega}t_{Z}F_{0}+F_{0}\frac{\partial(cF_{1})}{\partial\omega}-F_{-1}c^{\dagger}\frac{\partial
F_{0}}{\partial\omega}.$ (107)

Introducing

$\displaystyle
p=\bar{v}_{Z}+t_{Z}^{\dagger}\bar{A}^{-1}_{1,1}t_{Z}-t_{Z-1}\bar{A}^{-1}_{Z-1,Z-1}t_{Z-1}^{\dagger}=r+2t_{Z}^{\dagger}\bar{A}^{-1}_{1,1}t_{Z}$
(108)

we show that

$\displaystyle\sum_{m=1}^{\infty}G_{0,m}^{(0)}(\omega)G_{m,0}^{(0)}(\omega)=\frac{1}{2}F_{0}\frac{\partial
p}{\partial\omega}F_{0}-\frac{1}{2}\frac{\partial F_{0}}{\partial\omega}$
$\displaystyle+\frac{1}{2}\frac{\partial(F_{-1}c^{\dagger}-F_{1}c)}{\partial\omega}F_{0}-\frac{1}{2}\frac{\partial
F_{0}}{\partial\omega}(cF_{1}-c^{\dagger}F_{-1}),$ (109)
$\displaystyle\sum_{m=1}^{\infty}G_{0,m}^{(0)}(\omega)G_{m,0}^{(0)}(\omega)=\frac{1}{2}F_{0}\frac{\partial
p}{\partial\omega}F_{0}-\frac{1}{2}\frac{\partial F_{0}}{\partial\omega}$
$\displaystyle+\frac{1}{2}F_{0}\frac{\partial(cF_{1}-c^{\dagger}F_{-1})}{\partial\omega}-\frac{1}{2}(F_{-1}c^{\dagger}-F_{1}c)\frac{\partial
F_{0}}{\partial\omega}.$ (110)

Summing halves of each expression we obtain

$\displaystyle\text{tr}\left\\{[G_{0,0}^{(0)}]^{-1}\sum_{m=1}^{\infty}G_{0,m}^{(0)}(\omega)G_{m,0}^{(0)}(\omega)\right\\}$
$\displaystyle=\frac{1}{2}\text{tr}\left\\{\frac{\partial
p}{\partial\omega}F_{0}\right\\}-\frac{1}{2}\text{tr}\left\\{F_{0}^{-1}\frac{\partial
F_{0}}{\partial\omega}\right\\}$
$\displaystyle-\frac{1}{4}\text{tr}\left\\{F_{0}^{-1}\frac{\partial
F_{0}}{\partial\omega}(cF_{1}-c^{\dagger}F_{-1})\right\\}$
$\displaystyle-\frac{1}{4}\text{tr}\left\\{(F_{-1}c^{\dagger}-F_{1}c)\frac{\partial
F_{0}}{\partial\omega}F_{0}^{-1}\right\\}.$ (111)

###  A.2 Identities relating components of bulk Green’s functions

In this subsection we review properties of the matrices $F_{n}$ defined in Eq.
(95) and and prove some useful identities associated with them.

Rewriting the identity in Eq. (105) in the matrix form we find

$\displaystyle
F^{-1}_{n,n^{\prime}}F_{n^{\prime},n^{\prime\prime}}=\delta_{n,n^{\prime\prime}},$
(112)

with $F_{n^{\prime},n^{\prime\prime}}\equiv F_{n^{\prime}-n^{\prime\prime}}$
and

$\displaystyle
F^{-1}_{n,n^{\prime}}=c\delta_{n+1,n^{\prime}}+c^{\dagger}\delta_{n-1,n^{\prime}}-r\delta_{n,n^{\prime}}.$
(113)

We view $F^{-1}_{n,n^{\prime}}$ as a Hamiltonian of the effective tight-
binding model, $F_{n^{\prime},n^{\prime\prime}}$ being the corresponding
translationally invariant Green’s function. Perturbing this model by cutting
the link between the zeroth and first sites,

$\displaystyle(F_{S})^{-1}_{n,n^{\prime}}$
$\displaystyle=F^{-1}_{n,n^{\prime}}+(V_{S})_{n,n^{\prime}}$ (114)
$\displaystyle=F^{-1}_{n,n^{\prime}}-c\delta_{n,0}\delta_{n^{\prime},1}-c^{\dagger}\delta_{n,1}\delta_{n^{\prime},0},$
(115)

we split this system into two disconnected parts. Since the two parts are
completely independent we find $(F_{S})_{n,n^{\prime}}=0$ for $n\leq 0$,
$n^{\prime}\geq 1$ and for $n\geq 1$, $n^{\prime}\leq 0$. The function $F_{S}$
is called the surface Green’s function and obeys the following Dyson equations

$\displaystyle F_{S}$ $\displaystyle=F-FV_{S}F_{S}=F-F_{S}V_{S}F,$ (116)
$\displaystyle(F_{S})_{n,n^{\prime}}$
$\displaystyle=F_{n,n^{\prime}}+F_{n,0}c(F_{S})_{1,n^{\prime}}+F_{n,1}c^{\dagger}(F_{S})_{0,n^{\prime}}$
(117)
$\displaystyle=F_{n,n^{\prime}}+(F_{S})_{n,0}cF_{1,n^{\prime}}+(F_{S})_{n,1}c^{\dagger}F_{0,n^{\prime}}.$
(118)

Choosing $n=0$, $n^{\prime}=1$ we find

$\displaystyle 0=F_{-1}+F_{0}c(F_{S})_{1,1}=F_{-1}+(F_{S})_{0,0}cF_{0}$ (119)

and choosing $n=1$, $n^{\prime}=0$ we find

$\displaystyle
0=F_{1}+F_{0}c^{\dagger}(F_{S})_{0,0}=F_{1}+(F_{S})_{1,1}c^{\dagger}F_{0}.$
(120)

Comparing the two expressions, we find the following identities

$\displaystyle F_{0}c^{\dagger}F_{-1}$ $\displaystyle=F_{1}cF_{0},$ (121)
$\displaystyle F_{0}cF_{1}$ $\displaystyle=F_{-1}c^{\dagger}F_{0}.$ (122)

Next, focusing on the two semi-infinite Green’s functions
$(F_{R})_{n,n^{\prime}}=(F_{S})_{n\geq 1,n^{\prime}\geq 1}$,
$(F_{L})_{n,n^{\prime}}=(F_{S})_{n\leq 0,n^{\prime}\leq 0}$, and using the
condition that they remain invariant under adding one more site to the
corresponding semi-infinite models, i.e.

$\displaystyle F_{R}=\left(\begin{array}[]{cc}F_{R}^{-1}&c^{\dagger}\\\
c&-r\end{array}\right)^{-1},\quad
F_{L}=\left(\begin{array}[]{cc}-r&c^{\dagger}\\\
c&F_{L}^{-1}\end{array}\right)^{-1},$ (127)

we compare the $(1,1)$ matrix elements on both sides of the first relation (in
the lower right corner), and the $(0,0)$ matrix elements on both sides of the
second relation (in the upper left corner). To this end we use Eq. (84). We
find

$\displaystyle(F_{S})_{1,1}$
$\displaystyle=(-r-c(F_{S})_{1,1}c^{\dagger})^{-1},$ (128)
$\displaystyle(F_{S})_{0,0}$
$\displaystyle=(-r-c^{\dagger}(F_{S})_{0,0}c)^{-1}.$ (129)

Eliminating $(F_{S})_{1,1}$, $(F_{S})_{0,0}$, and $r$ with the help of Eqs.
(119), (120), (105) we obtain the following quadratic matrix equations

$\displaystyle(c^{\dagger}F_{-1})^{2}-c^{\dagger}F_{-1}-c^{\dagger}F_{0}cF_{0}$
$\displaystyle=0,$ (130)
$\displaystyle(F_{1}c)^{2}-F_{1}c-F_{0}c^{\dagger}F_{0}c$ $\displaystyle=0,$
(131) $\displaystyle(cF_{1})^{2}-cF_{1}-cF_{0}c^{\dagger}F_{0}$
$\displaystyle=0,$ (132)
$\displaystyle(F_{-1}c^{\dagger})^{2}-F_{-1}c^{\dagger}-F_{0}cF_{0}c^{\dagger}$
$\displaystyle=0.$ (133)

Their solutions yields

$\displaystyle c^{\dagger}F_{-1}$
$\displaystyle=\frac{1-\sqrt{1+4c^{\dagger}F_{0}cF_{0}}}{2},$ (134)
$\displaystyle F_{1}c$
$\displaystyle=\frac{1-\sqrt{1+4F_{0}c^{\dagger}F_{0}c}}{2},$ (135)
$\displaystyle cF_{1}$
$\displaystyle=\frac{1-\sqrt{1+4cF_{0}c^{\dagger}F_{0}}}{2},$ (136)
$\displaystyle F_{-1}c^{\dagger}$
$\displaystyle=\frac{1-\sqrt{1+4F_{0}cF_{0}c^{\dagger}}}{2},$ (137)

where the sign choice for the square roots is substantiated by the
perturbative expansion of $F_{1}$ and $F_{-1}$ in $c$, $c^{\dagger}$.

Combining Eq. (131) with Eq. (121), we show that

$\displaystyle F_{0}=(1-F_{1}c)(F_{0}-F_{1}F_{0}^{-1}F_{-1}).$ (138)

Analogous identities hold for the shifted system, that is for the tilded
functions.

###  A.3 Derivation of (42)

Performing a shift of the lattice leftwards by one site, as suggested in Sec.
III.2, we effectively redefine the unit cell (i.e. making it begin from site
$j=2$ and terminating it after site $j=1$). The corresponding Bloch
Hamiltonian defined in Eq. (72) changes to

$\displaystyle\tilde{h}_{k}-(\omega+i\eta)=\left(\begin{array}[]{cc}\tilde{A}&\tilde{b}_{k}\\\
\tilde{b}_{k}^{\dagger}&\bar{v}_{1}\end{array}\right),$ (141)

where $\bar{v}_{1}=v_{1}-(\omega+i\eta)$ and

$\displaystyle\tilde{b}_{k}=-\left(\begin{array}[]{c}t_{1}e^{-ik}\\\ 0\\\
\vdots\\\ 0\\\
t_{Z}^{\dagger}\end{array}\right),\quad\tilde{b}_{k}^{\dagger}=-\left(e^{ik}t_{1}^{\dagger},0,\ldots,0,t_{Z}\right).$
(147)

The matrix $\tilde{A}$ has the structure similar to that of $\bar{A}$ — they
differ from each other only by the labelling of sites. Moreover, they share a
common block matrix $a$ of the size $N_{c}(Z-2)\times N_{c}(Z-2)$. This
becomes obvious in the following representation

$\displaystyle\bar{A}=\left(\begin{array}[]{cc}\bar{v}_{1}&-T_{1}^{\dagger}\\\
-T_{1}&a\end{array}\right),\quad\tilde{A}=\left(\begin{array}[]{cc}a&-T_{Z-1}\\\
-T_{Z-1}^{\dagger}&\bar{v}_{Z}\end{array}\right),$ (152)

with

$\displaystyle T_{1}$ $\displaystyle=\left(\begin{array}[]{c}t_{1}\\\
\vdots\\\ 0\end{array}\right),\quad
T_{1}^{\dagger}=\left(t_{1}^{\dagger},\ldots,0\right),$ (156) $\displaystyle
T_{Z-1}$ $\displaystyle=\left(\begin{array}[]{c}0\\\ \vdots\\\
t_{Z-1}^{\dagger}\end{array}\right),\quad
T_{Z-1}^{\dagger}=\left(0,\ldots,t_{Z-1}\right).$ (160)

Using the identity in Eq. (84) and the additional relation

$\displaystyle-M_{11}^{-1}M_{12}(M_{22}-M_{21}M_{11}^{-1}M_{12})^{-1}$
$\displaystyle=$
$\displaystyle-(M_{11}-M_{12}M_{22}^{-1}M_{21})^{-1}M_{12}M_{22}^{-1},$ (161)

we evaluate

$\displaystyle\bar{A}_{1,1}^{-1}$
$\displaystyle=(\bar{v}_{1}-t_{1}^{\dagger}a_{2,2}^{-1}t_{1})^{-1}\equiv
V_{1}^{-1},$ (162) $\displaystyle\bar{A}_{Z-1,Z-1}^{-1}$
$\displaystyle=a_{Z-1,Z-1}^{-1}+a_{Z-1,2}^{-1}t_{1}V_{1}^{-1}t_{1}^{\dagger}a_{2,Z-1}^{-1},$
(163) $\displaystyle\bar{A}_{1,Z-1}^{-1}$
$\displaystyle=V_{1}^{-1}t_{1}^{\dagger}a_{2,Z-1}^{-1},$ (164)
$\displaystyle\bar{A}_{Z-1,1}^{-1}$
$\displaystyle=a_{Z-1,2}^{-1}t_{1}V_{1}^{-1},$ (165)

and

$\displaystyle\tilde{A}_{1,1}^{-1}$
$\displaystyle=a_{2,2}^{-1}+a_{2,Z-1}^{-1}t_{Z-1}^{\dagger}V_{Z}^{-1}t_{Z-1}a^{-1}_{Z-1,2},$
(166) $\displaystyle\tilde{A}_{Z-1,Z-1}^{-1}$
$\displaystyle=(\bar{v}_{Z}-t_{Z-1}a_{Z-1,Z-1}^{-1}t_{Z-1}^{\dagger})^{-1}\equiv
V_{Z}^{-1},$ (167) $\displaystyle\tilde{A}_{1,Z-1}^{-1}$
$\displaystyle=a_{2,Z-1}^{-1}t_{Z-1}^{\dagger}V_{Z}^{-1},$ (168)
$\displaystyle\tilde{A}_{Z-1,1}^{-1}$
$\displaystyle=V_{Z}^{-1}t_{Z-1}a_{Z-1,2}^{-1}.$ (169)

Note that we label the blocks of the matrix $a$ beginning with $2$ and ending
with $Z-1$.

Using these representations we express Eqs. (91), (92), (93), and (108) as
follows

$\displaystyle c$
$\displaystyle=t_{Z}^{\dagger}V_{1}^{-1}\bar{t}_{1}^{\dagger},$ (170)
$\displaystyle c^{\dagger}$ $\displaystyle=\bar{t}_{1}V_{1}^{-1}t_{Z},$ (171)
$\displaystyle r$
$\displaystyle=V_{Z}-t_{Z}^{\dagger}V^{-1}_{1}t_{Z}-\bar{t}_{1}V_{1}^{-1}\bar{t}_{1}^{\dagger},$
(172) $\displaystyle p$
$\displaystyle=V_{Z}+t_{Z}^{\dagger}V^{-1}_{1}t_{Z}-\bar{t}_{1}V_{1}^{-1}\bar{t}_{1}^{\dagger},$
(173)

where we have introduced the notation
$\bar{t}_{1}=t_{Z-1}a_{Z-1,2}^{-1}t_{1}$. In addition, we find analogs of
these quantities for the shifted system

$\displaystyle\tilde{c}$
$\displaystyle=t_{1}^{\dagger}\tilde{A}^{-1}_{1,Z-1}t_{Z}^{\dagger}=\bar{t}_{1}^{\dagger}V_{Z}^{-1}t_{Z}^{\dagger},$
(174) $\displaystyle\tilde{c}^{\dagger}$
$\displaystyle=t_{Z}\tilde{A}^{-1}_{Z-1,1}t_{1}=t_{Z}V_{Z}^{-1}\bar{t}_{1},$
(175) $\displaystyle\tilde{r}$
$\displaystyle=\bar{v}_{1}-t_{1}^{\dagger}\tilde{A}^{-1}_{1,1}t_{1}-t_{Z}\tilde{A}^{-1}_{Z-1,Z-1}t_{Z}^{\dagger}$
$\displaystyle=V_{1}-\bar{t}_{1}^{\dagger}V_{Z}^{-1}\bar{t}_{1}-t_{Z}V^{-1}_{Z}t_{Z}^{\dagger},$
(176) $\displaystyle\tilde{p}$
$\displaystyle=\bar{v}_{1}+t_{1}^{\dagger}\tilde{A}^{-1}_{1,1}t_{1}-t_{Z}\tilde{A}^{-1}_{Z-1,Z-1}t_{Z}^{\dagger}$
$\displaystyle=2\bar{v}_{1}-(V_{1}-\bar{t}_{1}^{\dagger}V_{Z}^{-1}\bar{t}_{1}+t_{Z}V^{-1}_{Z}t_{Z}^{\dagger})\equiv
2\bar{v}_{1}-\tilde{p}_{*}.$ (177)

As one can see, they all are expressed just in terms of the four matrices
$V_{1}$, $V_{Z}$, $\bar{t}_{1}$, and $t_{Z}$ (besides $\tilde{p}$, which has
an extra contribution $2\bar{v}_{1}$).

Defining

$\displaystyle\tilde{m}(k)$
$\displaystyle=\left(\tilde{c}e^{ik}+\tilde{c}^{\dagger}e^{-ik}-\tilde{r}\right)^{-1},$
(178) $\displaystyle\tilde{F}_{n}$
$\displaystyle=\int_{-\pi}^{\pi}\frac{dk}{2\pi}\tilde{m}(k)e^{ikn},$ (179)

and using Eqs. (33) and (111) for both the initial and the shifted systems we
evaluate

$\displaystyle\Delta Q^{\prime}_{B}$
$\displaystyle=\tilde{Q}^{\prime}_{B}-Q^{\prime}_{B}$
$\displaystyle=\frac{1}{\pi}\int d\omega\,\,\Theta(\mu-\omega)\,\text{Im}\,B$
(180) $\displaystyle-\frac{1}{\pi}\int d\omega\,\,\Theta(\mu-\omega)\,\text{Im
tr}\,\tilde{F}_{0},$ (181)

where

$\displaystyle B=$
$\displaystyle-\frac{1}{2}\text{tr}\left\\{\tilde{F}_{0}^{-1}\frac{\partial\tilde{F}_{0}}{\partial\omega}\right\\}+\frac{1}{2}\text{tr}\left\\{\tilde{F}_{0}\frac{\partial}{\partial\omega}\left(\tilde{\Lambda}^{\dagger}\tilde{F}_{0}^{-1}+\tilde{F}_{0}^{-1}\tilde{\Lambda}\right)\right\\}$
$\displaystyle+\frac{1}{2}\text{tr}\left\\{F_{0}^{-1}\frac{\partial
F_{0}}{\partial\omega}\right\\}+\frac{1}{2}\text{tr}\left\\{F_{0}\frac{\partial}{\partial\omega}\left(\Lambda
F_{0}^{-1}+F_{0}^{-1}\Lambda^{\dagger}\right)\right\\},$ (182)

and

$\displaystyle\tilde{\Lambda}$
$\displaystyle=\frac{\tilde{F}_{-1}\tilde{c}^{\dagger}-\tilde{F}_{1}\tilde{c}}{2}-\tilde{F}_{0}\frac{\tilde{p}_{*}}{2},$
(183) $\displaystyle\Lambda$
$\displaystyle=\frac{c^{\dagger}F_{-1}-cF_{1}}{2}-\frac{p}{2}F_{0}.$ (184)

Noticing that the contribution in Eq. (181) equals $\rho_{1}^{(0)}$ (since
$\tilde{F}_{0}=\tilde{G}_{0,0}^{(0)}=G_{1,1}^{(0)}$) we achieve exact
cancellation of the corresponding term in (41). Thereby we get

$\displaystyle I=\frac{1}{\pi}\int
d\omega\,\,\Theta(\mu-\omega)\,\text{Im}\,B.$ (185)

On the basis of Eqs. (121), (122) we conclude that

$\displaystyle\tilde{F}_{0}\tilde{\Lambda}^{\dagger}$
$\displaystyle=\tilde{\Lambda}\tilde{F}_{0},$ (186) $\displaystyle
F_{0}\Lambda$ $\displaystyle=\Lambda^{\dagger}F_{0}.$ (187)

This observation allows us to write down $B$ in a more compact form

$\displaystyle B=$
$\displaystyle-\frac{1}{2}\text{tr}\left\\{\tilde{F}_{0}^{-1}\frac{\partial\tilde{F}_{0}}{\partial\omega}\right\\}+\text{tr}\left\\{\tilde{F}_{0}\frac{\partial}{\partial\omega}\left(\tilde{F}_{0}^{-1}\tilde{\Lambda}\right)\right\\}$
$\displaystyle+\frac{1}{2}\text{tr}\left\\{F_{0}^{-1}\frac{\partial
F_{0}}{\partial\omega}\right\\}+\text{tr}\left\\{F_{0}\frac{\partial}{\partial\omega}\left(\Lambda
F_{0}^{-1}\right)\right\\}.$ (188)

With help of the identities (see the next section for the proof)

$\displaystyle\Lambda$
$\displaystyle=t_{Z}^{\dagger}\tilde{\Lambda}(t_{Z}^{\dagger})^{-1},$ (189)
$\displaystyle t_{Z}^{\dagger}\tilde{F}_{0}t_{Z}$
$\displaystyle=\left(\Lambda^{2}-\frac{1}{4}\right)F_{0}^{-1},$ (190)

we express

$\displaystyle B$ $\displaystyle=\text{tr}\left\\{F_{0}^{-1}\frac{\partial
F_{0}}{\partial\omega}\right\\}-\text{tr}\left\\{(\Lambda-\frac{1}{2})^{-1}\frac{\partial}{\partial\omega}(\Lambda-\frac{1}{2})\right\\}$
(191)
$\displaystyle=\frac{\partial}{\partial\omega}\ln\det\left[F_{0}(\Lambda-\frac{1}{2})^{-1}\right]$
(192)
$\displaystyle=-\frac{\partial}{\partial\omega}\ln\det\left[F_{0}^{-1}(\Lambda^{\dagger}-\frac{1}{2})\right],$
(193)

where we have applied Jacobi’s formula in the second step. Taking into account
that $F_{0}=G_{0,0}^{(0)}$ as well as the relation

$\displaystyle G_{0,1}^{(0)}=(\frac{1}{2}-\Lambda^{\dagger})t_{Z}^{-1},$ (194)

following from Eq. (100), we finally obtain Eq. (42).

###  A.4 Identities relating components of boundary Green’s functions

Our goal is to prove the relations in Eqs. (189), (190).

Let us treat the shifted system as the union of the first $Z-1$ sites and the
rest of the semi-infinite lattice (which coincides with the initial system up
to relabeling of sites $m\to m-Z+1$). They are coupled with each other by the
hoppings $-t_{Z}$ and $-t_{Z}^{\dagger}$, which occur between sites $Z-1$ and
$Z$. Formally, this separation is expressed as

$\displaystyle\tilde{G}_{m,m^{\prime}}^{-1}$
$\displaystyle=(\tilde{G}^{\prime})_{m,m^{\prime}}^{-1}$
$\displaystyle+t_{Z}\delta_{m,Z}\delta_{m^{\prime},Z-1}+t_{Z}^{\dagger}\delta_{m,Z-1}\delta_{m^{\prime},Z},$
(195) $\displaystyle(\tilde{G}^{\prime})_{m,m^{\prime}}^{-1}$
$\displaystyle=-\tilde{A}_{m,m^{\prime}}\delta_{1\leq m,m^{\prime}\leq Z-1}$
$\displaystyle+G_{m-Z+1,m^{\prime}-Z+1}^{-1}\delta_{m,m^{\prime}\geq Z},$
(196)

where $\tilde{G}^{\prime}$ is the Green’s function of the system in the
absence of coupling between its two subsystems. It is convenient to rewrite
Eq. (195) in the form of the Dyson equation

$\displaystyle(\tilde{G}^{\prime})_{m,m^{\prime}}=\tilde{G}_{m,m^{\prime}}$
$\displaystyle+\tilde{G}_{m,Z}t_{Z}(\tilde{G}^{\prime})_{Z-1,m^{\prime}}$
$\displaystyle+\tilde{G}_{m,Z-1}t_{Z}^{\dagger}(\tilde{G}^{\prime})_{Z,m^{\prime}}.$
(197)

Choosing $m=Z$ and $m^{\prime}=Z-1,Z$, and using that
$(\tilde{G}^{\prime})_{Z,Z-1}=(\tilde{G}^{\prime})_{Z-1,Z}=0$ (the case of two
isolated subsystems), we obtain the following identities

$\displaystyle 0$
$\displaystyle=\tilde{G}_{Z,Z-1}-\tilde{G}_{Z,Z}t_{Z}V^{-1}_{Z},$ (198)
$\displaystyle G_{1,1}$
$\displaystyle=\tilde{G}_{Z,Z}+\tilde{G}_{Z,Z-1}t_{Z}^{\dagger}G_{1,1}.$ (199)

Using the general expression Eq. (22) for the boundary Green’s functions we
express

$\displaystyle\tilde{G}_{Z,Z}$
$\displaystyle=\tilde{G}_{Z,Z}^{(0)}-\tilde{G}_{Z,0}^{(0)}[\tilde{G}_{0,0}^{(0)}]^{-1}\tilde{G}_{0,Z}^{(0)}$
(200)
$\displaystyle=\tilde{F}_{0}-\tilde{F}_{1}\tilde{F}_{0}^{-1}\tilde{F}_{-1}.$
(201)

By virtue of the tilded analog of the identity in Eq. (138) we establish

$\displaystyle\tilde{G}_{Z,Z}^{-1}$
$\displaystyle=\tilde{F}_{0}^{-1}(1-\tilde{F}_{1}\tilde{c}).$ (202)

To prove Eq. (189) we represent $\tilde{G}_{Z,Z-1}$ by means of the identity
in Eq. (22)

$\displaystyle\tilde{G}_{Z,Z-1}$
$\displaystyle=\tilde{G}_{Z,Z-1}^{(0)}-\tilde{G}_{Z,0}^{(0)}[\tilde{G}_{0,0}^{(0)}]^{-1}\tilde{G}_{0,Z-1}^{(0)}$
$\displaystyle=G_{1,0}^{(0)}-\tilde{F}_{1}\tilde{F}_{0}^{-1}\tilde{G}_{0,Z-1}^{(0)},$
(203)

where we have used
$\tilde{G}_{m-1,m^{\prime}-1}^{(0)}=G_{m,m^{\prime}}^{(0)}$, and hence
$\tilde{G}_{Z,Z-1}^{(0)}=G_{Z+1,Z}^{(0)}=G_{1,0}^{(0)}$. In turn, representing
$\tilde{G}_{0,Z-1}^{(0)}$ with help of the tilded analog of Eq. (98), we
obtain

$\displaystyle\tilde{G}_{0,Z-1}^{(0)}$
$\displaystyle=\tilde{F}_{0}\tilde{t}_{Z}^{\dagger}\tilde{A}^{-1}_{1,Z-1}+\tilde{F}_{-1}\tilde{t}_{Z-1}\tilde{A}^{-1}_{Z-1,Z-1}$
(204)
$\displaystyle=(\tilde{F}_{0}\tilde{c}+\tilde{F}_{-1}t_{Z}V_{Z}^{-1}t_{Z}^{\dagger})(t_{Z}^{\dagger})^{-1},$
(205)

where we used have $\tilde{t}_{Z}=t_{1}$ and $\tilde{t}_{Z-1}=t_{Z}$.
Combining Eqs. (203), (205) with Eqs. (198), (201), and (194) we obtain the
relation

$\displaystyle(t_{Z}^{\dagger})^{-1}(\frac{1}{2}-\Lambda)t_{Z}^{\dagger}$
$\displaystyle=\tilde{F}_{0}t_{Z}V^{-1}_{Z}t_{Z}^{\dagger}+\tilde{F}_{1}\tilde{c}$
(206) $\displaystyle=\frac{1}{2}-\tilde{\Lambda},$ (207)

which proves Eq. (189). To prove Eq. (190) we eliminate $\tilde{G}_{Z,Z-1}$
from Eqs. (198) and (199). We obtain

$\displaystyle
G_{1,1}=\left(\tilde{G}_{Z,Z}^{-1}-t_{Z}V^{-1}_{Z}t_{Z}^{\dagger}\right)^{-1}.$
(208)

Then we substitute Eq. (202) into this relation and obtain

$\displaystyle G_{1,1}$
$\displaystyle=\left(1-\tilde{F}_{1}\tilde{c}-\tilde{F}_{0}t_{Z}V^{-1}_{Z}t_{Z}^{\dagger}\right)^{-1}\tilde{F}_{0}$
(209)
$\displaystyle=\left(\frac{1}{2}+\tilde{\Lambda}\right)^{-1}\tilde{F}_{0}=(t_{Z}^{\dagger})^{-1}\left(\frac{1}{2}+\Lambda\right)^{-1}t_{Z}^{\dagger}\tilde{F}_{0}.$
(210)

On the other hand, on the basis of Eq. (22), it holds

$\displaystyle G_{1,1}$
$\displaystyle=G_{1,1}^{(0)}-G_{1,0}^{(0)}[G_{0,0}^{(0)}]^{-1}G_{0,1}^{(0)}$
(211)
$\displaystyle=\tilde{F}_{0}-(t_{Z}^{\dagger})^{-1}(\Lambda-\frac{1}{2})F_{0}^{-1}(\Lambda^{\dagger}-\frac{1}{2})t_{Z}^{-1}$
(212)
$\displaystyle=(t_{Z}^{\dagger})^{-1}\left[t_{Z}^{\dagger}\tilde{F}_{0}t_{Z}-(\Lambda-\frac{1}{2})^{2}F_{0}^{-1}\right]t_{Z}^{-1}.$
(213)

Comparing Eq. (210) with Eq. (213) we deduce Eq. (190).

##  Appendix B Green’s functions of the interface model with a weakened link

Below we establish the Green’s function of the model introduced in Sec. IV.2.

Eliminating $G_{0,m^{\prime}}$ from Eq. (52)

$\displaystyle G_{0,m^{\prime}}$
$\displaystyle=[1-G_{0,1}^{(0)}t_{\lambda}]^{-1}[G_{0,m^{\prime}}^{(0)}+G_{0,0}^{(0)}t_{\lambda}^{\dagger}G_{1,m^{\prime}}],$
(214)

and inserting this result into Eq. (53), we obtain the following linear
equation for $G_{1,m^{\prime}}$:

$\displaystyle G_{1,m^{\prime}}$
$\displaystyle=G_{1,m^{\prime}}^{(0)}+G_{1,1}^{(0)}t_{\lambda}[1-G_{0,1}^{(0)}t_{\lambda}]^{-1}G_{0,m^{\prime}}^{(0)}$
$\displaystyle+\\{G_{1,1}^{(0)}t_{\lambda}[1-G_{0,1}^{(0)}t_{\lambda}^{\dagger}]^{-1}G_{0,0}^{(0)}t_{\lambda}^{\dagger}+G_{1,0}^{(0)}t_{\lambda}^{\dagger}\\}G_{1,m^{\prime}}.$
(215)

Its solution reads

$\displaystyle G_{1,m^{\prime}}$
$\displaystyle=S^{-1}\\{G_{1,m^{\prime}}^{(0)}+G_{1,1}^{(0)}t_{\lambda}[1-G_{0,1}^{(0)}t_{\lambda}]^{-1}G_{0,m^{\prime}}^{(0)}\\},$
(216)

where

$\displaystyle S$
$\displaystyle=1-G_{1,0}^{(0)}t_{\lambda}^{\dagger}-G_{1,1}^{(0)}t_{\lambda}[1-G_{0,1}^{(0)}t_{\lambda}]^{-1}G_{0,0}^{(0)}t_{\lambda}^{\dagger}.$
(217)

Analogously we find

$\displaystyle G_{0,m^{\prime}}=$
$\displaystyle[1-G_{0,1}^{(0)}t_{\lambda}]^{-1}G_{0,0}^{(0)}t_{\lambda}^{\dagger}S^{-1}$
$\displaystyle\times\\{[1-G_{1,0}^{(0)}t_{\lambda}^{\dagger}](t_{\lambda}^{\dagger})^{-1}G_{0,0}^{(0)-1}G_{0,m^{\prime}}^{(0)}+G_{1,m^{\prime}}^{(0)}\\}.$
(218)

We note the useful matrix identity

$\displaystyle
1+G_{0,0}^{(0)}t_{\lambda}^{\dagger}S^{-1}G_{1,1}^{(0)}t_{\lambda}[1-G_{0,1}^{(0)}t_{\lambda}]^{-1}$
$\displaystyle=G_{0,0}^{(0)}t_{\lambda}^{\dagger}S^{-1}[1-G_{1,0}^{(0)}t_{\lambda}^{\dagger}](t_{\lambda}^{\dagger})^{-1}G_{0,0}^{(0)-1},$
(219)

which helps us achieve various representations for the above solutions. We
also recall that $G_{0,0}^{(0)}=F_{0}$, $G_{1,1}^{(0)}=\tilde{F}_{0}$, and
$G_{0,1}^{(0)}$ can be expressed via Eq. (194).

Using the identities in Eqs. (187), (190) we simplify the expression in Eq.
(217) as follows

$\displaystyle
S=(t_{Z}^{\dagger})^{-1}\frac{1-(1-\lambda^{2})\left(\frac{1}{2}-\Lambda\right)}{1-(1-\lambda)(\frac{1}{2}-\Lambda)}t_{Z}^{\dagger}=\frac{1-(1-\lambda^{2})\left(\frac{1}{2}-\tilde{\Lambda}\right)}{1-(1-\lambda)(\frac{1}{2}-\tilde{\Lambda})}.$
(220)

Consequently, we find

$\displaystyle G_{1,m^{\prime}}$
$\displaystyle=\frac{1-(1-\lambda)(\frac{1}{2}-\tilde{\Lambda})}{1-(1-\lambda^{2})\left(\frac{1}{2}-\tilde{\Lambda}\right)}G_{1,m^{\prime}}^{(0)}$
$\displaystyle+\frac{1-\lambda}{1-(1-\lambda^{2})(\frac{1}{2}-\tilde{\Lambda})}\tilde{F}_{0}t_{Z}G_{0,m^{\prime}}^{(0)},$
(221) $\displaystyle G_{0,m^{\prime}}=$
$\displaystyle\frac{1-(1-\lambda)(\frac{1}{2}-\Lambda^{\dagger})}{1-(1-\lambda^{2})(\frac{1}{2}-\Lambda^{\dagger})}G_{0,m^{\prime}}^{(0)}$
$\displaystyle+\frac{1-\lambda}{1-(1-\lambda^{2})(\frac{1}{2}-\Lambda^{\dagger})}F_{0}t_{Z}^{\dagger}G_{1,m^{\prime}}^{(0)}.$
(222)

Next, we are interested in

$\displaystyle(1-\lambda)^{-1}\sum_{m=-\infty}^{\infty}\text{tr}\\{G_{m,m}-G_{m,m}^{(0)}\\}$
$\displaystyle=\sum_{m=-\infty}^{\infty}\text{tr}\\{G_{m,1}^{(0)}t_{Z}G_{0,m}+G_{m,0}^{(0)}t_{Z}^{\dagger}G_{1,m}\\}$
$\displaystyle=\sum_{m=-\infty}^{\infty}\text{tr}\\{G_{m,1}^{(0)}t_{Z}\frac{1-(1-\lambda)(\frac{1}{2}-\Lambda^{\dagger})}{1-(1-\lambda^{2})(\frac{1}{2}-\Lambda^{\dagger})}G_{0,m}^{(0)}\\}$
$\displaystyle+\sum_{m=-\infty}^{\infty}\text{tr}\\{G_{m,1}^{(0)}t_{Z}\frac{1-\lambda}{1-(1-\lambda^{2})(\frac{1}{2}-\Lambda^{\dagger})}F_{0}t_{Z}^{\dagger}G_{1,m}^{(0)}\\}$
$\displaystyle+\sum_{m=-\infty}^{\infty}\text{tr}\\{G_{m,0}^{(0)}t_{Z}^{\dagger}\frac{1-(1-\lambda)(\frac{1}{2}-\tilde{\Lambda})}{1-(1-\lambda^{2})\left(\frac{1}{2}-\tilde{\Lambda}\right)}G_{1,m}^{(0)}\\}$
$\displaystyle+\sum_{m=-\infty}^{\infty}\text{tr}\\{G_{m,0}^{(0)}t_{Z}^{\dagger}\frac{1-\lambda}{1-(1-\lambda^{2})(\frac{1}{2}-\tilde{\Lambda})}\tilde{F}_{0}t_{Z}G_{0,m}^{(0)}\\}.$
(223)

Considering

$\displaystyle\sum_{m=-\infty}^{\infty}G_{0,m}^{(0)}G_{m,0}^{(0)}$
$\displaystyle=\int_{-\pi}^{\pi}\frac{dk}{2\pi}\sum_{j=1}^{Z}\langle
Z|\frac{1}{\omega+i\eta-h_{k}}|j\rangle$ $\displaystyle\times\langle
j|\frac{1}{\omega+i\eta-h_{k}}|Z\rangle$ $\displaystyle=-\frac{\partial
F_{0}}{\partial\omega},$ (224)
$\displaystyle\sum_{m=-\infty}^{\infty}G_{1,m}^{(0)}G_{m,1}^{(0)}$
$\displaystyle=\sum_{m=-\infty}^{\infty}\tilde{G}_{0,m-1}^{(0)}\tilde{G}_{m-1,0}^{(0)}=-\frac{\partial\tilde{F}_{0}}{\partial\omega},$
(225) $\displaystyle\sum_{m=-\infty}^{\infty}G_{0,m}^{(0)}G_{m,1}^{(0)}$
$\displaystyle=\int_{-\pi}^{\pi}\frac{dk}{2\pi}e^{-ik}\sum_{j=1}^{Z}\langle
Z|\frac{1}{\omega+i\eta-h_{k}}|j\rangle$ $\displaystyle\times\langle
j|\frac{1}{\omega+i\eta-h_{k}}|1\rangle$ $\displaystyle=-\frac{\partial
G_{0,1}^{(0)}}{\partial\omega}=\frac{\partial\Lambda^{\dagger}}{\partial\omega}t_{Z}^{-1},$
(226)

we obtain

$\displaystyle\sum_{m=-\infty}^{\infty}\text{tr}\\{G_{m,m}-G_{m,m}^{(0)}\\}$
$\displaystyle=(1-\lambda)\text{tr}\left\\{\frac{1-(1-\lambda)(\frac{1}{2}-\Lambda^{\dagger})}{1-(1-\lambda^{2})(\frac{1}{2}-\Lambda^{\dagger})}\frac{\partial\Lambda^{\dagger}}{\partial\omega}\right\\}$
$\displaystyle-(1-\lambda)\text{tr}\left\\{t_{Z}\frac{1-\lambda}{1-(1-\lambda^{2})(\frac{1}{2}-\Lambda^{\dagger})}F_{0}t_{Z}^{\dagger}\frac{\partial\tilde{F}_{0}}{\partial\omega}\right\\}$
$\displaystyle+(1-\lambda)\text{tr}\left\\{\frac{1-(1-\lambda)(\frac{1}{2}-\Lambda)}{1-(1-\lambda^{2})\left(\frac{1}{2}-\Lambda\right)}\frac{\partial\Lambda}{\partial\omega}\right\\}$
$\displaystyle-(1-\lambda)\text{tr}\left\\{\frac{1-\lambda}{1-(1-\lambda^{2})(\frac{1}{2}-\Lambda)}t_{Z}^{\dagger}\tilde{F}_{0}t_{Z}\frac{\partial
F_{0}}{\partial\omega}\right\\}.$ (227)

Using again the identities in Eqs. (187), (190) we derive the expression

$\displaystyle\sum_{m=-\infty}^{\infty}\text{tr}\\{G_{m,m}-G_{m,m}^{(0)}\\}$
$\displaystyle=\frac{1}{2}\text{tr}\left\\{\frac{1}{\frac{1+\lambda^{2}}{2(1-\lambda^{2})}+\Lambda}\frac{\partial\Lambda}{\partial\omega}\right\\}+\frac{1}{2}\text{tr}\left\\{\frac{1}{\frac{1+\lambda^{2}}{2(1-\lambda^{2})}+\Lambda^{\dagger}}\frac{\partial\Lambda^{\dagger}}{\partial\omega}\right\\}$
$\displaystyle=\frac{\partial}{\partial\omega}\ln\det\left[(1-\lambda^{2})\Lambda+\frac{1+\lambda^{2}}{2}\right],$
(228)

from which follows Eq. (54).

##  Appendix C Topological invariants in the single-channel case

In this section we show how to recover the earlier obtained results for the
single-channel case from the presently analysed multichannel expressions.

###  C.1 Boundary charge invariant

In the single-channel case, $t_{j}$ and $v_{j}$ are scalar parameters. As
explained in Ref. [pletyukhov_etal_prb_20, ], by an appropriate gauge choice,
one can make all hoppings real-valued, i.e. $t_{j}=t_{j}^{*}$. Therefore,
$c=c^{*}=\frac{1}{\det\bar{A}}\prod_{j=1}^{Z}t_{j}$, $F_{1}=F_{-1}$, and

$\displaystyle\Lambda^{*}$ $\displaystyle=-\frac{p}{2}F_{0},$ (229)
$\displaystyle F_{0}$
$\displaystyle=\int_{-\pi}^{\pi}\frac{dk}{2\pi}\frac{1}{2c\cos
k-r}=\frac{f_{0}}{2c},$ (230) $\displaystyle f_{0}$
$\displaystyle=\int_{-\pi}^{\pi}\frac{dk}{2\pi}\frac{1}{\cos
k-D}=-\frac{1}{D\sqrt{1-\frac{1}{D^{2}}}},$ (231) $\displaystyle D$
$\displaystyle=\frac{r}{2c}=\frac{1}{2\bar{t}^{Z}}\left(d_{1,Z}-t_{Z}^{2}d_{2,Z-1}\right),$
(232) $\displaystyle\frac{p}{2c}$
$\displaystyle=D+\frac{t_{Z}^{2}}{\bar{t}^{Z}}d_{2,Z-1},$ (233)

where $\bar{t}^{Z}=\prod_{j=1}^{Z}t_{j}$, and $d_{j,j^{\prime}}$ are the
determinants of tridiagonal matrices which start from $\bar{v}_{j}$ and end up
with $\bar{v}_{j^{\prime}}$, such that $j\leq j^{\prime}$. The function
$D(\omega)$ determines the dispersion of the bulk problem by virtue of the
equation $D(\omega=\epsilon_{k\alpha})=\cos k$. The edge state equation in the
boundary problem reads $s\equiv d_{1,Z-1}=0$.

Then we express

$\displaystyle\partial_{\omega}\ln\det\left([G_{0,0}^{(0)}]^{-1}G_{0,1}^{(0)}t_{Z}\right)=\partial_{\omega}\ln\left(p+\frac{1}{F_{0}}\right)$
$\displaystyle=$
$\displaystyle\partial_{\omega}\left[\ln\left(\frac{t_{Z}^{2}}{\bar{t}^{Z}}d_{2,Z-1}+D+\frac{1}{f_{0}}\right)-\ln
s\right].$ (234)

The argument in the first term is complex valued only above/below branch cuts,
which correspond to bands. For each band $\epsilon_{k\alpha}$, lying in the
energy range
$\epsilon_{\alpha}^{(-)}<\epsilon_{k\alpha}<\epsilon_{\alpha}^{(+)}$, we can
alternatively express this term as an integral over the Bloch momentum $k$ in
the Brillouin zone $-\pi<k<\pi$. We notice that

$\displaystyle\frac{1}{f_{0}(\omega+i\eta)}=$
$\displaystyle-D(\omega)\sqrt{1-\frac{1}{D^{2}(\omega+i\eta)}}$ (235)
$\displaystyle=$
$\displaystyle-i\,\text{sign}[D(\omega)D^{\prime}(\omega)]D(\omega)\sqrt{\frac{1}{D^{2}(\omega)}-1}$
$\displaystyle=$ $\displaystyle\,i\,\sin
k\,\text{sign}\left(\frac{d\epsilon_{k\alpha}}{dk}\right).$ (236)

For bands with odd $\alpha$, the inequality $\frac{d\epsilon_{k\alpha}}{dk}>0$
holds for $0<k<\pi$, while for even $\alpha$ it holds for $-\pi<k<0$.
Therefore, the integration along the upper branch (see Fig. 1) gives

$\displaystyle-\frac{1}{\pi}\int_{\epsilon_{\alpha}^{(-)}}^{\epsilon_{\alpha}^{(+)}}d\omega\frac{\partial}{\partial\omega}\ln\left(\frac{t_{Z}^{2}}{\bar{t}^{Z}}d_{2,Z-1}(\omega)+D(\omega)+\frac{1}{f_{0}(\omega+i\eta)}\right)$
$\displaystyle=\begin{cases}-\frac{1}{\pi}\int_{0}^{\pi}dk\frac{\partial}{\partial
k}\ln\left(\frac{t_{Z}^{2}}{\bar{t}^{Z}}d_{2,Z-1}(\epsilon_{k\alpha})+e^{ik}\right),&\text{odd}\,\,\alpha,\\\
-\frac{1}{\pi}\int_{-\pi}^{0}dk\frac{\partial}{\partial
k}\ln\left(\frac{t_{Z}^{2}}{\bar{t}^{Z}}d_{2,Z-1}(\epsilon_{k\alpha})+e^{ik}\right),&\text{even}\,\,\alpha.\end{cases}$
(237)

Taking the imaginary part of these expressions and identifying (up to the
normalization)
$\frac{t_{Z}^{2}}{\bar{t}^{Z}}d_{2,Z-1}(\epsilon_{k\alpha})+e^{ik}\sim\chi_{k\alpha}(1)e^{ik}$
(see Ref. [pletyukhov_etal_prb_20, ] for details), we obtain the $\alpha$th
band’s contribution to the invariant in Eq. (42)

$\displaystyle I_{\alpha}$ $\displaystyle=-\frac{1}{2\pi
i}\int_{-\pi}^{\pi}dk\frac{\partial}{\partial
k}\ln\left(\chi_{k\alpha}(1)e^{ik}\right)$ (238)
$\displaystyle=-1-\text{wn}[\chi_{k\alpha}(1)],$ (239)

that is $I_{\alpha}$ is expressed in terms of the winding number of the first
Bloch state vector component.

In band gaps, the term

$\displaystyle-\frac{1}{\pi}\int
d\omega\frac{\partial}{\partial\omega}\ln\left(\frac{t_{Z}^{2}}{\bar{t}^{Z}}d_{2,Z-1}+D+\frac{1}{f_{0}}\right)$
(240)

picks up the contribution $+1$ from the so called left edge state pole of the
initial system (i.e. with $\kappa<0$) and from the right edge state pole of
the shifted system (i.e. with $\tilde{\kappa}>0$). In turn, the term

$\displaystyle\frac{1}{\pi}\int d\omega\frac{\partial}{\partial\omega}\ln s$
(241)

picks up the contribution $+1$ from both the left ($\kappa<0$) and right
($\kappa>0$) edge state poles of the initial system.

Thus, the invariant defined in Eq. (42) in the single channel case amounts to

$\displaystyle I=\sum_{\alpha=1}^{\nu}I_{\alpha}+\Delta Q_{E},$ (242)

where $\Delta Q_{E}$ denotes the difference between numbers of right edge
states of the shifted and initial systems summed up over all gaps below the
chemical potential $\mu$. Thereby we recover our earlier result expressed in
Eqs. (257), (258) of Ref. [pletyukhov_etal_prb_20, ].

Remarkably, on the basis of (190) we establish the identity

$\displaystyle\left(\frac{t_{Z}^{2}}{\bar{t}^{Z}}d_{2,Z-1}+D+\frac{1}{f_{0}}\right)\left(\frac{t_{Z}^{2}}{\bar{t}^{Z}}d_{2,Z-1}+D-\frac{1}{f_{0}}\right)=\frac{t_{Z}^{2}}{\bar{t}^{2Z}}s\tilde{s},$
(243)

where $\tilde{s}=d_{2,Z}$. It allows us to relate

$\displaystyle\ln\left(\frac{t_{Z}^{2}}{\bar{t}^{Z}}d_{2,Z-1}+D+\frac{1}{f_{0}}\right)-\ln
s$ $\displaystyle=$
$\displaystyle\ln\tilde{s}-\ln\left(\frac{t_{Z}^{2}}{\bar{t}^{Z}}d_{2,Z-1}+D-\frac{1}{f_{0}}\right).$
(244)

The equation $\tilde{s}=0$ is known to determine both right
($\tilde{\kappa}>0$) and left ($\tilde{\kappa}<0$) edge states in the shifted
system. In turn, $\frac{t_{Z}^{2}}{\bar{t}^{Z}}d_{2,Z-1}+D-\frac{1}{f_{0}}=0$
holds at the left edge state of the shifted system ($\tilde{\kappa}<0$) and at
the right edge of the initial system ($\kappa>0$). Therefore, this observation
additionally confirms the result in Eq. (242).

###  C.2 Interface charge invariant

In the single-channel case we find

$\displaystyle\Lambda+\frac{1}{2}$
$\displaystyle=-\frac{1}{2}\left(\frac{t_{Z}^{2}}{\bar{t}^{Z}}d_{2,Z-1}+D-\frac{1}{f_{0}}\right)f_{0}.$
(245)

Therefore, the contribution to Eq. (54) from band $\alpha$ may be written as

$\displaystyle\tilde{Q}_{I,\alpha}=$
$\displaystyle-\frac{1}{\pi}\int_{\epsilon_{\alpha}^{(-)}}^{\epsilon_{\alpha}^{(+)}}d\omega\,\,\text{Im}\,\frac{\partial}{\partial\omega}\ln
f_{0}$ (246)
$\displaystyle-\frac{1}{\pi}\int_{\epsilon_{\alpha}^{(-)}}^{\epsilon_{\alpha}^{(+)}}d\omega\,\,\text{Im}\,\frac{\partial}{\partial\omega}\ln\left[-\left(\frac{t_{Z}^{2}}{\bar{t}^{Z}}d_{2,Z-1}+D-\frac{1}{f_{0}}\right)\right.$
$\displaystyle\left.\qquad+\lambda^{2}\left(\frac{t_{Z}^{2}}{\bar{t}^{Z}}d_{2,Z-1}+D+\frac{1}{f_{0}}\right)\right].$
(247)

The term in Eq. (246) receives contributions $-\frac{1}{2}$ from each band
edge. In the term Eq. (247) we make the same change of integration variable
$\omega\to k$ as discussed in Sec. C.1. This results in the expression

$\displaystyle Q_{I,\alpha}+1$ $\displaystyle=$ $\displaystyle-\frac{1}{2\pi
i}\int_{-\pi}^{\pi}dk\,\frac{\partial}{\partial
k}\ln\left[-e^{-i\varphi_{k\alpha}(1)-ik}+\lambda^{2}e^{i\varphi_{k\alpha}(1)+ik}\right]$
$\displaystyle=$ $\displaystyle 1+\text{wn}[\chi_{k\alpha}(1)]$
$\displaystyle-\frac{1}{2\pi i}\int_{-\pi}^{\pi}dk\,\frac{\partial}{\partial
k}\ln\left[-1+\lambda^{2}e^{2i\varphi_{k\alpha}(1)+2ik}\right],$ (248)

where $e^{i\varphi_{k\alpha}(1)}$ is the phase of $\chi_{k\alpha}(1)$. For
$\lambda^{2}<1$ the argument in the last line does not encircle the origin,
therefore the corresponding contribution vanishes. We finally get

$\displaystyle Q_{I,\alpha}=\text{wn}\,[\chi_{k\alpha}(1)],$ (249)

thereby reproducing our earlier result in Eq. (C35) of Ref.
[Pletyukhov_etal3_2020, ] (assuming the gauge fixed by the condition
$\text{Im}\,[\chi_{k\alpha}(Z)]=0$).

## References

  * (1) K. von Klitzing, G. Dorda, and M. Pepper, Phys. Rev. Lett. 45, 494 (1980).
  * (2) D. J. Thouless, M. Kohmoto, M. P. Nightingale, and M. den Nijs, Phys. Rev. Lett. 49, 405 (1982).
  * (3) M. Konig, S. Wiedmann, C. Brune, A. Roth, H. Buhmann, L. W. Molenkamp, X.-L. Qi, and S.-C. Zhang, Science 318, 766 (2007).
  * (4) M. He, H. Sun, and Q. L. He, Frontiers of Physics 14, 43401 (2019).
  * (5) R. W. Bomantara and J. Gong, Phys. Rev. B 98, 165421 (2018).
  * (6) R. S. K. Mong, D. J. Clarke, J. Alicea, N. H. Lindner, P. Fendley, C. Nayak, Y. Oreg, A. Stern, E Berg, K. Shtengel, and M. P. A. Fisher, Phys. Rev. X 4, 011036 (2014).
  * (7) A. Stern and N. H. Lindner, Science 339, 1179 (2013).
  * (8) A. Miyake, Phys. Rev. Lett. 105, 040501 (2010).
  * (9) D. Pesin and A. H. MacDonald, Nat. Mater. 11, 409 (2012).
  * (10) M. Pletyukhov D. M. Kennes, J. Klinovaja, D. Loss, and H. Schoeller, Phys. Rev. B 101, 165304 (2020).
  * (11) M. Pletyukhov, D. M. Kennes, K. Piasotski, J. Klinovaja, D. Loss, and H. Schoeller, Phys. Rev. Research 2, 033345 (2020).
  * (12) R. Jackiw and C. Rebbi, Phys. Rev. D 13, 3398 (1976).
  * (13) W.P. Su, J.R. Schrieffer, and A.J. Heeger, Phys. Rev. Lett. 42, 1698 (1979); Phys. Rev. B 22, 2099 (1980).
  * (14) R. Jackiw and J.R. Schrieffer, Nucl. Phys. B 190, 253 (1981).
  * (15) M.J. Rice and E.J. Mele, Phys. Rev. Lett. 49, 1455 (1982).
  * (16) S. Kivelson, Phys. Rev. B 28, 2653 (1983).
  * (17) W.P. Su and J.R. Schrieffer, Phys. Rev. Lett. 46, 738 (1981).
  * (18) A.J. Heeger, S. Kivelson, J.R. Schrieffer, and W.-P. Su, Rev. Mod. Phys. 60, 781 (1988).
  * (19) E. Witten, Phys. Lett. B 86, 283 (1979).
  * (20) H. Takayama, Y.R. Lin-Liu, and K. Maki, Phys. Rev. B 21, 2388 (1980).
  * (21) J. Goldstone and F. Wilczek, Phys. Rev. Lett. 47, 986 (1981).
  * (22) R. Jackiw and G. Semenoff, Phys. Rev. Lett. 50, 439 (1983).
  * (23) R. D. King-Smith and D. Vanderbilt, Phys. Rev. B 47, 1651(R) (1993).
  * (24) D. Vanderbilt and R. D. King-Smith, Phys. Rev. B 48, 4442 (1993).
  * (25) R. Resta, Rev. Mod. Phys. 66, 899 (1994).
  * (26) K. N. Kudin and R. Car, J. Chem. Phys. 126, 234101 (2007).
  * (27) N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and D. Vanderbilt, Rev. Mod. Phys. 84, 1419 (2012).
  * (28) N. A. Spaldin, J. Solid State Chem. 195, 2 (2012).
  * (29) J.-W. Rhim, J. Behrends, and J. H. Bardarson, Phys. Rev. B 95, 035421 (2017).
  * (30) G. van Miert and C. Ortix, Phys. Rev. B 96, 235130 (2017).
  * (31) J. Zak, Phys. Rev. Lett. 62, 2747 (1989).
  * (32) M. Pletyukhov, D. M. Kennes, J. Klinovaja, D. Loss, and H. Schoeller, Phys. Rev. B 101, 161106(R) (2020).
  * (33) J.-H. Park, G. Yang, J. Klinovaja, P. Stano, and D. Loss, Phys. Rev. B 94, 075416 (2016).
  * (34) M. Thakurathi, J. Klinovaja, and D. Loss, Phys. Rev. B 98, 245404 (2018).
  * (35) Y.-T. Lin, D. M. Kennes, M. Pletyukhov, C. S. Weber, H. Schoeller, and V. Meden, Phys. Rev. B 102, 085122 (2020).
  * (36) C. S. Weber, K. Piasotski, M. Pletyukhov, J. Klinovaja, D. Loss, H. Schoeller, and D. M. Kennes, Phys. Rev. Lett. 126, 016803 (2021).
  * (37) Y.-T. Lin, C. S. Weber, D. M. Kennes, M. Pletyukhov, H. Schoeller, and V. Meden, Phys. Rev. B 103, 195119 (2021).
  * (38) K. Laubscher, C. S. Weber, D. M. Kennes, M. Pletyukhov, H. Schoeller, D. Loss, and J. Klinovaja, arXiv:2101.10301 (2021), accepted in Phys. Rev. B (2021).
  * (39) W. Kohn, Phys. Rev. Lett. 76, 3168 (1996).
  * (40) E. Prodan and W. Kohn, PNAS 102, 11635 (2005).
  * (41) L. Arrachea, G. S. Lozano, and A. A. Aligia, Phys. Rev. B 80, 014425 (2009).
  * (42) A. Zazunov, R. Egger, and A. Levy Yeyati, Phys. Rev. B 94, 014502 (2016).
  * (43) Y. Peng, Y. Bao, and F. von Oppen, Phys. Rev. B 95, 235143 (2017).
  * (44) A. Zazunov, R. Egger, M. Alvarado, and A. Levy Yeyati, Phys. Rev. B 96, 024516 (2017).
  * (45) A. Komnik and S. Heinze, Phys. Rev. B 96, 155103 (2017).
  * (46) A. Zazunov, A. Iks, M. Alvarado, A. Levy Yeyati, and R. Egger, Beilstein J. Nanotechnol. 9, 1659 (2018).
  * (47) S. Pinon, V. Kaladzhyan, and C. Bena, Phys. Rev. B 101, 115405 (2020).
  * (48) M. Alvarado, A. Iks, A. Zazunov, R. Egger, and A. Levy Yeyati, Phys. Rev. B 101, 094511 (2020).
  * (49) A. M. Essin and V. Gurarie, Phys. Rev. B 84, 125132 (2011).
  * (50) V. Gurarie, Phys. Rev. B 83, 085426 (2011).
  * (51) Yu. A. Bychkov and A. M. Dikhne, JETP Lett. 3, 202 (1966).
  * (52) Yu. A. Bychkov and A. M. Dykhne, Sov. Phys. JETP 24, 1285 (1967).
  * (53) Z. A. Kazamanyan, Sov. Phys. JETP 61, 1215 (1971).
  * (54) S. Aubry and G. André, Ann. Israel Phys. Soc. 3, 133 (1980).
  * (55) P. G. Harper, Proc. Phys. Soc. A 68, 874 (1955).
  * (56) L. Fidkowski, T. S. Jackson, and I. Klich, Phys. Rev. Lett. 107, 036601 (2011).
  * (57) R. S. K. Mong and V. Shivamoggi, Phys. Rev. B 83, 125109 (2011).
  * (58) T. Fukui, K. Shiozaki, T. Fujiwara, and S. Fujimoto, J. Phys. Soc. Jpn. 81, 114602 (2012).
  * (59) Y. Yu, Y.-S. Wu, and X. Xie, Nucl. Phys. B 916, 550 (2017).
  * (60) D. J. Thouless, Phys.Rev. B 27, 6083 (1983).
  * (61) Q. Niu and D. J. Thouless, J. Phys. A: Math. Gen. 17, 2453 (1984).
  * (62) K. Piasotski, N. Müller, D. Kennes, H. Schoeller, and M. Pletyukhov, unpublished.
  * (63) J.-W. Rhim, J. H. Bardarson, and R.-J. Slager, Phys. Rev. B 97, 115143 (2018).
  * (64) See Supplemental Material for additional details of the models’ parameters.