††thanks: Electronic Address<EMAIL_ADDRESS>Electronic
Address<EMAIL_ADDRESS>Electronic Address:
<EMAIL_ADDRESS>
# Many-Body Quadrupolar Sum Rule for Higher-Order Topological Insulator

Wonjun Lee  Department of Physics, Pohang University of Science and Technology
(POSTECH), Pohang 37673, Republic of Korea  Center for Artificial Low
Dimensional Electronic Systems, Institute for Basic Science (IBS), Pohang
37673, Republic of Korea  Gil Young Cho  Department of Physics, Pohang
University of Science and Technology (POSTECH), Pohang 37673, Republic of
Korea  Center for Artificial Low Dimensional Electronic Systems, Institute for
Basic Science (IBS), Pohang 37673, Republic of Korea  Asia Pacific Center for
Theoretical Physics, Pohang 37673, Republic of Korea  Byungmin Kang  School of
Physics, Korea Institute for Advanced Study, Seoul 02455, Korea

###### Abstract

The modern theory of polarization establishes the bulk-boundary correspondence
for the bulk polarization. In this paper, we attempt to extend it to a sum
rule of the bulk quadrupole moment by employing a many-body operator
introduced in [Phys. Rev. B 100, 245134 (2019)] and [Phys. Rev. B 100, 245135
(2019)]. The sum rule that we propose consists of the alternating sum of four
observables, which are the phase factors of the many-body operator in
different boundary conditions. We demonstrate its validity through extensive
numerical computations for various non-interacting tight-binding models. We
also observe that individual terms in the sum rule correspond to the bulk
quadrupole moment, the edge-localized polarizations, and the corner charge in
the thermodynamic limit on some models.

##  I Introduction

Recent developments in topological insulators Hasan and Kane (2010); Qi and
Zhang (2011) discovered a new class of topological material, called higher-
order topological insulators Benalcazar  _et al._ (2017a, b); Schindler  _et
al._ (2018). Higher-order topological insulators are characterized by non-
trivial boundary-of-boundary states despite that the boundaries are trivial.
Among various theoretical tools, the bulk polarization has proven useful in
understanding many aspects of topological insulators Fu and Kane (2006); Coh
and Vanderbilt (2009); Hughes  _et al._ (2011); Yu  _et al._ (2011);
Alexandradinata  _et al._ (2014a). It is therefore natural to ask if the bulk
quadrupole or higher multipole moments could play the same role for higher-
order topological insulators.

The modern theory of polarization King-Smith and Vanderbilt (1993); Resta
(1993); Vanderbilt and King-Smith (1993) identifies the bulk polarization with
the sum of Wannier centers. Later, the bulk polarization of a many-body state
is shown to be identified as the phase factor of the expectation value of a
many-body operator Resta (1998). Despite the success of the modern theory of
polarization, the modern theory of quadrupole and higher multipole moments has
not been fully developed so far. One attempt is to employ the Wannier centers
to define the bulk quadrupole moment, but this approach has been successful
when there exist crystalline symmetries so that the bulk multipole moment is
quantized Benalcazar  _et al._ (2017b); Watanabe and Ono (2020). On the other
line of development, Refs. Kang  _et al._ , 2019; Wheeler  _et al._ , 2019
introduced a many-body operator $\hat{U}_{2}$ for the quadrupole moment
generalizing the many-body operator presented in Ref. Resta, 1998. This
approach can in principle be applied to systems without any crystalline
symmetries and to interacting systems, however, the observables given by
$\hat{U}_{2}$ have difficulties in identifying their physical meanings with
respect to their coordinate dependence Kang  _et al._ (2019) and lack of
periodicity Ono  _et al._ (2019) so their physical meaning have not been fully
understood yet.

One of the success of the modern theory of polarization is that it identifies
the bulk polarization with the boundary charge, which is called the bulk-
boundary correspondence Vanderbilt and King-Smith (1993); Rhim  _et al._
(2017). While the bulk polarization is determined by a full many-body electron
state, the boundary charge is determined by a simple one-body observable near
the boundary. Thus, the identification of the two is rather surprising and
also provides a firm justification of the modern theory of polarization. A
natural question is then whether such a correspondence for the bulk quadrupole
can be formulated. In recent studies Trifunovic (2020); Ren  _et al._ (2021),
bulk-boundary correspondences for the quadrupole moment were presented using
Wannier functions, where the corner charge is expressed as a sum of quantities
that depend on the choice of the bulk Wannier functions, although the corner
charge itself is not. We also note another recent study Ying , where the
corner charge of band insulators is computed using the adiabatic current
flowing along the edges.

In this paper, we propose another bulk-boundary correspondence in terms of the
operator $\hat{U}_{2}$, which is given by an alternating sum of four phase
factors, where the sum vanishes in the thermodynamic limit:

$\phi_{pp}-\phi_{op}-\phi_{po}+\phi_{oo}=0\quad\text{(mod 1)}.$

Four phase factors are the phase factors of the expectation values of
$\hat{U}_{2}$ with respect to the ground state in different boundary
conditions,

$\phi_{ab}=\frac{1}{2\pi}\text{Im}\Big{[}\log\big(\langle\textrm{GS}(a,b)|\hat{U}_{2}|\textrm{GS}(a,b)\rangle\big{missing})\Big{]},$

where ‘$a$’ and ‘$b$’ refer to as the boundary conditions along $x$\- and
$y$-direction, respectively, with ‘$a$’ and ‘$b$’ being either ‘$p$’
(periodic) or ‘$o$’ (open). We observe that for band insulators having well-
localized edge-localized polarizations defined in Ref. Benalcazar  _et al._ ,
2017b based on the hybrid Wannier function (HWF), the phase factors
$\phi_{pp}$, $\phi_{po/op}$, and $\phi_{oo}$ can be identified with the bulk
quadrupole moment, the edge-localized polarizations, and the corner charge,
respectively, which consequently elevates our bulk-boundary correspondence to
that for the quadrupole moment. We further observe that for some insulators,
there exists HWFs with the corresponding hybrid Wannier value at $0.5$ so that
one cannot directly apply the definition of the edge-localized polarization
presented in Ref. Benalcazar  _et al._ , 2017b due to its branch cut
dependence, as elaborated in Appendix D. For those insulators, we find that at
least the phase factor $\phi_{oo}$ can be regarded as the corner charge in the
thermodynamic limit while the physical meaning of the other phase factors
becomes unclear.

The paper is organized as follows. In Sec. II, we introduce sum rules encoding
the bulk-boundary correspondence for multipole moments using many-body
operators. In Sec. III, we provide numerical results presenting our
observations in Sec. II. We conclude in Sec. IV where the summary of main
results and possible future directions are discussed. In Appendix A, we derive
the dipolar and quadrupolar sum rules for classical systems. In Appendix B, we
provide a proof that the many-body operator can measure the corner charge in
one-dimensional geometry for band insulators. In Appendix C, we provide
numerical results for a $C_{3}$-symmetric insulator, which suggests that the
sum rule may hold in other types of symmetric insulators besides the
$C_{4}$-symmetric ones discussed in the main text. Finally, in Appendix D, we
provide technical remarks on the HWF-based edge-localized polarization
including the branch cut dependence which becomes more explicit for the models
with the hybrid Wannier value at $0.5$.

##  II Multipolar Sum Rules

In the following, we present the multipolar sum rules for quantum systems in
one- and two-dimension. For concreteness, we note here some remarks on the sum
rules. First, we only focus on one-dimensional systems having non-zero bulk
polarization in Sec. II.1 and two-dimensional systems with vanishing bulk
polarization in Sec. II.2. Second, we focus on circle/line geometry for one-
dimensional systems and torus/cylinder/rectangle geometry for two-dimensional
systems. Each geometry is characterized by open or periodic boundary
conditions. In addition, we always use hard-cutoff boundary conditions for
open boundaries. Third, for simplicity we restrict the position of each
orbital at a Bravais lattice site. Fourth, for all the quantum systems
presented below, we assume translation symmetry and $C_{2}$ symmetry, which
includes the cases where the corner charge is not quantized. Finally, for
convenience, we set the electron charge to 1 in the following discussions.

###  II.1 Dipolar Sum Rule

In this subsection, we introduce the dipolar sum rule for many-body quantum
systems as a generalization of the classical dipolar sum rule. The dipolar sum
rule relates the bulk polarization with the boundary charge. To motivate the
quantum mechanical dipolar sum rule, we first state the classical dipolar sum
rule,

$\mathcal{Q}_{c}=\mathcal{P},$ (1)

where $\mathcal{P}$ and $\mathcal{Q}_{c}$ are the classical bulk polarization
and boundary charge, respectively. The derivation of Eq. (1) and the
definitions of the electric moments can be found in Appendix A.1.

We now generalize the classical dipolar sum rule Eq. (1) to a one-dimensional
quantum mechanical system. To this end, we first need to define the bulk
polarization and the boundary charge for a quantum system. The boundary charge
can be defined as a one-body observable:

$Q^{(1)}_{c}\equiv\sum_{x=1}^{L_{x}/2}\big{(}\rho(x)-n_{e}\big{)}\quad\text{(mod
1)},$ (2)

where $L_{x}$ is the total system size,
$\rho(x)=\langle\hat{n}_{x}\rangle=\langle\sum_{\alpha=1}^{N_{\textrm{orb}}}c_{x,\alpha}^{\dagger}c_{x,\alpha}\rangle$
is the charge density at site $x$ with the number of orbitals per site $N_{\rm
orb}$, and $n_{e}$ is the average number of electrons per site. We then employ
the Resta’s many-body operator Resta (1998)

$\hat{U}_{1}\equiv\exp(\frac{2\pi
i}{L_{x}}\sum_{x=1}^{L_{x}}x(\hat{n}_{x}-n_{e})),$ (3)

where we include the background charge $n_{e}$ from ions sitting at each
lattice site. Using $\hat{U}_{1}$, the bulk polarization is given by

$P=\frac{1}{2\pi}\imaginary[\log\big(\langle\textrm{GS}|\hat{U}_{1}|\textrm{GS}\rangle\big{missing})\Big{]},$
(4)

where $|\textrm{GS}\rangle$ is the many-body ground state subject to the
periodic boundary condition. With these, the classical dipolar sum rule Eq.
(1) generalizes as Vanderbilt and King-Smith (1993); Rhim  _et al._ (2017)

$Q^{(1)}_{c}=P\quad\text{(mod 1)}$ (5)

in quantum systems. Unlike the dipole sum rule in classical systems, Eq. (5)
holds only modulo 1, the unit of electron charge.

Having obtained the sum rule for the dipole moment, we consider the further
characterization solely in terms of the many-body operator $\hat{U}_{1}$. For
this, we define the phase factor $\phi_{a}$ of the expectation value of
$\hat{U}_{1}$ as

$\phi_{a}\equiv\frac{1}{2\pi}\text{Im}\Big{[}\log\big(\langle\textrm{GS}(a)|\hat{U}_{1}|\textrm{GS}(a)\rangle\big{missing})\Big{]},$
(6)

where $\ket{\textrm{GS}(a)}$ is the ground state under the boundary conditions
of type ‘$a$’ which can be either ‘$p$’ (periodic) or ‘$o$’ (open).

While it seems not widely known, the boundary charge can also be captured from
$\phi_{o}$ in the thermodynamic limit, i.e.,

$Q^{(1)}_{c}=\phi_{o}\quad\text{(mod 1)}$ (7)

as the system size goes to infinity when the system is gapped. Our proof of
Eq. (7) for band insulators can be found in Appendix B. Thus, by combining
Eqs. (4) and (7), the dipole sum rule can be succinctly recast in terms of the
phase factors:

$\phi_{p}=\phi_{o}\quad\text{(mod 1)}.$ (8)

###  II.2 Quadrupolar Sum Rule

In this subsection, we introduce the quadrupolar sum rule for many-body
quantum systems, which generalizes the quantum mechanical dipolar sum rule to
the quadrupole case, and discuss its difficulties for certain cases in
identifying the physical meaning of individual terms in the sum rule.

The quadrupole sum rule would relate the bulk quadrupole moment, the edge-
localized polarization, and the corner charge Benalcazar  _et al._ (2017b);
Kang  _et al._ (2019) if these quantities are well-defined. To illustrate the
quadrupolar sum rule, we first need to state the classical quadrupolar sum
rule,

$\mathcal{Q}_{c}=-\mathcal{Q}_{xy}+\mathcal{P}^{\text{edge}}_{x}+\mathcal{P}^{\text{edge}}_{y},$
(9)

where $\mathcal{Q}_{xy}$, $\mathcal{P}^{\text{edge}}_{x/y}$, and
$\mathcal{Q}_{c}$ are the classical bulk-quadrupole moment, the edge-localized
polarizations, and the corner charge, respectively. The derivation of Eq. (9)
and the definitions of the classical electric moments can be found in Appendix
A.2.

A natural generalization of the classical quadrupolar sum rule to the quantum
mechanical one in two dimensions would be

$Q^{(2)}_{c}=-Q_{xy}+P^{\text{edge}}_{x}+P^{\text{edge}}_{y}\quad\text{(mod
1)},$ (10)

where $Q^{(2)}_{c}$ is the corner charge of a two dimensional system of the
size $L_{x}\times L_{y}$,

$Q^{(2)}_{c}\equiv\sum_{x=1}^{L_{x}/2}\sum_{y=1}^{L_{y}/2}\big{(}\rho(x,y)-n_{e}\big{)}\quad\text{(mod
1)},$ (11)

with $\rho(x,y)=\langle\hat{n}_{x,y}\rangle=\langle\sum_{\alpha=1}^{N_{\rm
orb}}c_{x,y,\alpha}^{\dagger}c_{x,y,\alpha}\rangle$ being the local charge
density, $Q_{xy}$ and $P^{\text{edge}}_{x/y}$ are the quantum mechanical bulk
quadrupole moment and the edge-localized polarizations, respectively. As a
side note, our corner charge Eq. (11) corresponds to the bare corner charge in
Ref. Ren  _et al._ , 2021. We remark that the edge-localized polarizations
$P^{\text{edge}}_{x/y}$ include the contribution from a dressing of polarized
one-dimensional chains along the boundaries Trifunovic  _et al._ (2019), which
does not affect the quadrupolar sum rule Eq. (10), and are fixed for a given
state. The precise forms of $Q_{xy}$ and $P^{\text{edge}}_{x/y}$ will be
presented below.

Similar to the dipole case, let us introduce the following many-body operator
Kang  _et al._ (2019); Wheeler  _et al._ (2019):

$\hat{U}_{2}\equiv\exp(\frac{2\pi
i}{L_{x}L_{y}}\sum_{x,y=1}^{L_{x},L_{y}}xy(\hat{n}_{x,y}-n_{e})),$ (12)

where $\hat{n}_{x,y}$ is the electron number operator at $(x,y)$ and $L_{x}$
($L_{y}$) is the linear system size in $x$-direction ($y$-direction). For the
open boundary, we assign the coordinate $x=1$ ($y=1$) to sites on the left
(bottom) boundary and the coordinate $x=L_{x}$ ($y=L_{y}$) to sites on the
right (top) boundary. For the periodic boundary, the $x$-coordinate
($y$-coordinate) could be assigned arbitrary as long as it starts from $1$ and
ends with $L_{x}$ ($L_{y}$) since it does not change the expectation value of
$\hat{U}_{2}$ with respect to a translation invariant state.

In order to formulate the quadrupolar sum rule in terms of the phase factors
of the expectation values of $\hat{U}_{2}$, we first introduce

$\phi_{ab}=\frac{1}{2\pi}\text{Im}\Big{[}\log\big(\langle\textrm{GS}(a,b)|\hat{U}_{2}|\textrm{GS}(a,b)\rangle\big{missing})\Big{]},$
(13)

where $\ket{\textrm{GS}(a,b)}$ is the ground state under the ‘$a$’ and ‘$b$’
boundary conditions along $x$\- and $y$-direction, respectively, with ‘$a$’
and ‘$b$’ being either ‘$p$’ (periodic) or ‘$o$’ (open). Each boundary
condition corresponds to different geometry, for example, ‘$pp$’ and ‘$po$’
mean the torus and the cylinder geometry, respectively. These phase factors
are gauge independent by construction.

Using these phase factors, we propose the bulk-boundary correspondence in
terms of the phase factors as

$\phi_{pp}-\phi_{op}-\phi_{po}+\phi_{oo}=0\quad\text{(mod 1)}$ (14)

and numerically find that this holds in the thermodynamic limit on band
insulators. We will call this as the sum rule of the phase factors, or simply
the sum rule. In addition, we also find that $\phi_{oo}$ agrees with
$Q^{(2)}_{c}$,

$Q^{(2)}_{c}=\phi_{oo}\quad\text{(mod 1)}.$ (15)

As a demonstration, we provide numerical confirmations of these on band
insulators in Sec. III. The above two equations, i.e., Eqs. (14) and (15), are
the main results of our present paper.

Based on the identification of the corner charge Eq. (15) and the comparison
between two sum rules Eqs. (10) and (14) suggest that $\phi_{pp}$ can be
naturally identified with the bulk quadrupole moment

$Q_{xy}=\phi_{pp}\quad\text{(mod 1)},$ (16)

which has been checked for the cases with non trivial bulk quadrupole moment
Kang  _et al._ (2019); Wheeler  _et al._ (2019). In addition, $\phi_{po/op}$
are identified with the edge-localized polarizations

$\begin{split}P^{\text{edge}}_{x}&=\phi_{po}\quad\text{(mod 1)}\\\
P^{\text{edge}}_{y}&=\phi_{op}\quad\text{(mod 1)}.\end{split}$ (17)

Whenever these identifications can be made, the sum rule Eq. (14) indeed
becomes the quadrupole sum rule for quantum systems. In Sec. III.1, we
numerically demonstrate the validity of the identifications of $\phi_{po/op}$
by comparing these with the HWF-based edge-localized polarizations defined in
Ref. Benalcazar  _et al._ , 2017b. However, it is important to note that these
identifications are not always possible. In Sec. III.2, we provide an example
on which these identifications fail as $\phi_{pp}$, $\phi_{op}$, and
$\phi_{po}$ are not convergent in the thermodynamic limit. Nonetheless, even
for these cases, the sum rule Eq. (14) and the identification of the corner
charge with $\phi_{oo}$ Eq. (15) always hold.

In the remainder of this section, we shortly discuss the HWF-based edge-
localized polarizationBenalcazar  _et al._ (2017b),

$\tilde{P}^{\textrm{edge}}_{x}=\sum_{j}\sum_{y=1}^{L_{y}/2}\nu^{j}\rho^{j}(y)\quad\text{(mod
1)},$ (18)

where $\rho^{j}(y)$ is the density, and $e^{2\pi i\nu_{j}}$ is the eigenvalue
of the $j$-th HWF $\psi^{j}(y)$ of the hybrid Wilson loop
$\mathcal{W}_{k_{x}}$ along the $x$-direction. We find that for band
insulators having well-localized $\tilde{P}^{\textrm{edge}}_{x/y}$ along the
boundaries, $\tilde{P}^{\textrm{edge}}_{x/y}$ and $\phi_{po/op}$ seem to agree
each other in the thermodynamic limit and the same quantized value for
insulators having $C_{4}$ symmetry,

$\begin{split}\tilde{P}^{\textrm{edge}}_{x}&=\phi_{op}\quad\text{(mod 1)}\\\
\tilde{P}^{\textrm{edge}}_{y}&=\phi_{po}\quad\text{(mod 1)}.\end{split}$ (19)

However, $\tilde{P}^{\textrm{edge}}_{x/y}$ crucially depends on the choice of
the branch cut of the hybrid Wannier values. In particular, for the models
with the hybrid Wannier value at $0.5$, this depdendance becomes more
explicit, thereby leading to a difficulty in computing
$\tilde{P}_{x/y}^{\textrm{edge}}$. We discuss this branch cut dependence more
in details in Appendix D.

##  III Numerical Demonstration

In this section, we provide numerical demonstrations of the sum rule Eq. (14)
and the observations, Eqs. (15) and (19). Due to size limitation, our numerics
are based on non-interacting tight-binding models on the square lattice.
However, we believe that the same sum rule should also hold in interacting
cases as well, since our formalism is based on many-body operators. In
addition, we separately discuss models without and with the hybrid Wannier
value at $0.5$ as there exists a difficulty in computing
$\tilde{P}^{\textrm{edge}}_{x/y}$ for the latter case. We therefore
numerically confirm the validity of Eq. (19) only for the former case while
the numerical confirmation of the sum rule Eq. (14) and Eq. (15) are presented
in all cases.

###  III.1 Models without the hybrid Wannier value at $0.5$

We present numerical results on models without the hybrid Wannier value at
$0.5$. Tested models are the quadrupole insulator, the edge-localized
polarization insulator, the quadrupole insulator with $\pi/2$-flux per
plaquette, and the two band model introduced in Ref. Zhou  _et al._ , 2015,
where details of these models can be found below. In the case of full open
boundary conditions, $C_{4}$ symmetry breaking term is always introduced in
order to split the possible degeneracy at the Fermi level. Due to such term,
the number of filled states is always equal to the filling
$n_{e}/N_{\textrm{orb}}$ times the number of sites $N$. The same $C_{4}$
breaking term often splits the degeneracy of the hybrid Wannier values at
$0.5$ in the case of mixed open and periodic boundary conditions.

For each model, we numerically compute the phase factors $\phi_{ab}$ Eq. (13),
the edge-localized polarizations $\tilde{P}_{x/y}^{\textrm{edge}}$ Eq. (18),
and the corner charge $Q_{c}^{(2)}$ Eq. (11) for various parameters. The
results are summarized in Tables 1, 2, and 3. All observables in the tables
are the extrapolated values in the thermodynamic limit $L\to\infty$ which are
obtained via quadratic extrapolations as a function of $1/L$, as shown in Fig.
1. In all cases, we find that the phase factors in the sum rule and the sum
rule itself are convergent in the thermodynamic limit. We also see that our
edge-localized polarizations, $\phi_{op/po}$, have values similar to the other
ones, $\tilde{P}_{x/y}^{\textrm{edge}}$, with differences at worst
$O(10^{-3})$. In addition, when a $C_{4}$ symmetry breaking term is small,
$\phi_{pp}$ has the same quantized value as the one predicted by the nested
Wilson-loop approach Benalcazar  _et al._ (2017b). Furthermore, the
identification of $\phi_{oo}$ with $Q^{(2)}_{c}$ and the sum rule hold up to
errors of $O(10^{-4})$. The errors may be attributed to the finite-size
effect.

####  III.1.1 Quadrupole insulator $H_{\text{quad}}$

As our first model with gapped Wannier spectrum, we consider the following
tight-binding model Benalcazar  _et al._ (2017b) having a non-zero bulk
quadrupole moment:

$\displaystyle H_{\textrm{quad}}(\boldsymbol{k})=$
$\displaystyle\big{(}\gamma_{x}+\lambda_{x}\cos(k_{x})\big{)}\Gamma_{4}+\lambda_{x}\sin(k_{x})\Gamma_{3}$
(20)
$\displaystyle+\big{(}\gamma_{y}+\lambda_{y}\cos(k_{y})\big{)}\Gamma_{2}+\lambda_{y}\sin(k_{y})\Gamma_{1}+\delta\,\Gamma_{0},$

where $\Gamma_{0}=\sigma_{3}\otimes\sigma_{0}$,
$\Gamma_{k}=-\sigma_{2}\otimes\sigma_{k}$ for $k=1,2,3$, and
$\Gamma_{4}=\sigma_{1}\otimes\sigma_{0}$ with $\sigma_{k=1,2,3}$ being Pauli
matrices and $\sigma_{0}$ being the $2\times 2$ identity matrix.

When $\delta=0$, there exists two anti-commuting mirror symmetries which
quantize the bulk quadrupole moment $Q_{xy}=0$ or $1/2$. The half-filled
ground state of Eq. (20) realizes topologically trivial quadrupole insulator
when $|\gamma_{x}|>|\lambda_{x}|$ and $|\gamma_{y}|>|\lambda_{y}|$ and
topologically non-trivial quadrupole insulator when
$|\gamma_{x}|<|\lambda_{x}|$ and $|\gamma_{y}|<|\lambda_{y}|$. When
$\delta\neq 0$, $C_{4}$ symmetry and the mirror symmetries are broken and
hence the bulk quadrupole moment is no longer quantized.

The numerical results for the quadrupole insulator are summarized in Table 1.

Model parameters | Phase factors | Electric moments | Sum rule
---|---|---|---
$H_{\text{quad}}(\gamma_{x},\gamma_{y},\lambda_{x},\lambda_{y},\delta)$ | $\phi_{pp}$ | $\phi_{po}$ | $\phi_{op}$ | $\phi_{oo}$ | $Q^{(2)}_{c}$ | $\tilde{P}_{x}^{\textrm{edge}}$ | $\tilde{P}_{y}^{\textrm{edge}}$ | $\sum(-1)^{ab}\phi_{ab}$
$(0.1,0.1,1.0,1.0,10^{-7})$ | 0.50000(0) | 0.50000(0) | 0.50000(0) | 0.50000(0) | 0.50000(0) | 0.50000(0) | 0.50000(0) | 0.00000(0)
$(0.5,0.5,1.0,1.0,0.7)$ | -0.11389(0) | -0.11385(0) | -0.11385(0) | -0.11392(0) | -0.11387(0) | -0.11386(0) | -0.11386(0) | -0.00010(0)
$(0.2,0.3,1.0,1.0,0.1)$ | -0.43156(0) | -0.43155(0) | -0.43155(0) | -0.43155(0) | -0.43156(0) | -0.43156(0) | -0.43155(0) | -0.00001(0)
$(0.1,0.1,1.0,1.2,0.1)$ | -0.44018(0) | -0.44017(0) | -0.44018(0) | -0.44019(0) | -0.44018(0) | -0.44018(0) | -0.44018(0) | -0.00002(0)
Table 1:  The phase factors $\phi_{ab}$ in Eq. (13) and the electric moments
$(Q_{c}^{(2)},\tilde{P}_{x}^{\text{edge}},\tilde{P}_{y}^{\text{edge}})$ in
Eqs. (11) and (18) are computed for the quadrupole insulator Eq. (20) with
various parameters. The sum of the phase factors $\sum(-1)^{ab}\phi_{ab}$ in
the last column refers to the combination
$\phi_{pp}-\phi_{po}-\phi_{op}+\phi_{oo}$ as per Eq. (14). All values are the
ones from extrapolating the observables as a function of $1/L$, where we
consider isotropic systems $L_{x}=L_{y}=L$ with $L$ from 23 to 30. An explicit
extrapolation procedure is presented in Fig. 1. We use round brackets to
denote the errors of the least significant digit. We see that the edge-
localized polarizations $\tilde{P}_{x}^{\textrm{edge}}$ and
$\tilde{P}_{y}^{\textrm{edge}}$ agree with phase factors $\phi_{po}$ and
$\phi_{op}$, corner charge $Q^{(2)}_{c}$ agrees with $\phi_{oo}$, and the sum
rule is satisfied up to small errors.  Figure 1: Extrapolation of the
observables (a) $\phi_{pp}$, (b) $\phi_{op}$ and
$\tilde{P}^{\text{edge}}_{y}$, (c) $\phi_{oo}$ and $Q^{(2)}_{c}$, as a
quadratic function in $1/L$ for the quadrupole insulator Eq. (20) with
parameters
$(\gamma_{x},\gamma_{y},\lambda_{x},\lambda_{y},\delta)=(0.5,0.5,1.0,1.0,0.7)$
and isotropic system sizes. The quadratic extrapolation works well and the
error from the extrapolation is less than $10^{-5}$ in all cases. (b)
$\phi_{op}$ and $\tilde{P}^{\text{edge}}_{y}$ agree with each other in
$O(10^{-5})$. (c) $\phi_{oo}$ and $Q^{(2)}_{c}$ agree with each other in the
same order of accuracy.

####  III.1.2 Edge-localized polarization insulator $H_{\text{edge}}$

Our next model is the edge-localized polarization insulator Benalcazar  _et
al._ (2017b). The tight-binding Hamiltonian in the momentum space can be
written as

$\displaystyle H_{\text{edge}}(\boldsymbol{k})=$
$\displaystyle\big{(}\gamma_{x}+\lambda_{x}\cos(k_{x})\big{)}\Gamma_{4}+\lambda_{x}\sin(k_{x})\Gamma_{3}$
(21)
$\displaystyle+\big{(}\gamma_{y}+\lambda_{y}\cos(k_{y})\big{)}\tilde{\Gamma}_{2}+\lambda_{y}\sin(k_{y})\tilde{\Gamma}_{1}+\delta\,\Gamma_{0},$

where $\tilde{\Gamma}_{2}=\sigma_{1}\otimes\sigma_{1}$,
$\tilde{\Gamma}_{1}=-\sigma_{1}\otimes\sigma_{2}$, and all the other $\Gamma$
matrices are the same as the ones in Eq. (20).

The half-filled ground state of Eq. (21) has a vanishing bulk quadrupole
moment but having a non-zero corner charge. In particular, when $\delta=0$,
there exist two mirror symmetries $M_{x}$ and $M_{y}$ which quantize the
corner charge to be $0$ or $1/2$, where the corner charge is originated sorely
from the edge-localized polarizations Benalcazar  _et al._ (2017b). The half-
filled ground state of Eq. (21) with
$\absolutevalue{\gamma_{x/y}}\ll\absolutevalue{\lambda_{x/y}}$ has
$(\tilde{P}_{x}^{\textrm{edge}},\tilde{P}_{y}^{\textrm{edge}})=(0,0.5)$ when
$\lambda_{x}>\lambda_{y}$ and
$(\tilde{P}_{x}^{\textrm{edge}},\tilde{P}_{y}^{\textrm{edge}})=(0.5,0)$ when
$\lambda_{x}<\lambda_{y}$.

The numerical results for the edge-localized polarization insulator are
summarized in Table 2.

Model parameters | Phase factors | Electric moments | Sum rule
---|---|---|---
$H_{\text{edge}}(\gamma_{x},\gamma_{y},\lambda_{x},\lambda_{y},\delta)$ | $\phi_{pp}$ | $\phi_{po}$ | $\phi_{op}$ | $\phi_{oo}$ | $Q^{(2)}_{c}$ | $\tilde{P}_{x}^{\textrm{edge}}$ | $\tilde{P}_{y}^{\textrm{edge}}$ | $\sum(-1)^{ab}\phi_{ab}$
$(0.1,0.1,1,1.2,10^{-7})$ | 0.00000(0) | 0.50000(0) | 0.00000(0) | 0.50000(0) | 0.50000(0) | 0.50000(0) | 0.00000(0) | 0.00000(0)
$(0.5,0.5,1.0,1.0,0.7)$ | -0.01661(0) | -0.09914(0) | -0.09914(0) | -0.18179(0) | -0.18174(0) | -0.09836(0) | -0.09836(0) | -0.00012(0)
$(0.2,0.3,1.0,1.0,0.1)$ | -0.05940(1) | -0.27871(1) | -0.31926(1) | 0.46139(1) | 0.46122(1) | -0.27833(1) | -0.31897(1) | -0.00004(1)
$(0.1,0.1,1.0,1.2,0.1)$ | -0.01585(0) | -0.44086(0) | -0.12064(0) | 0.45436(0) | 0.45434(0) | -0.44172(0) | -0.11936(0) | 0.00002(0)
Table 2:  The phase factors $\phi_{ab}$ and the electric moments
$(Q_{c}^{(2)},\tilde{P}_{x}^{\text{edge}},\tilde{P}_{y}^{\text{edge}})$ with
are computed for the edge-localized polarization insulator Eq. (21) with
various parameters. Here we use the same extrapolation procedure as in Table.
1. In this case as well, the edge-localized polarizations
$\tilde{P}_{x}^{\textrm{edge}}$ and $\tilde{P}_{y}^{\textrm{edge}}$ agree with
phase factors $\phi_{po}$ and $\phi_{op}$, corner charge $Q^{(2)}_{c}$ agrees
with $\phi_{oo}$, and the sum rule is satisfied up to small errors.

####  III.1.3 Quadrupole insulator with $\pi/2$-flux per plaquette
$H_{\text{quad}}^{(\pi/2)}$

Our third model is the quadrupole insulator with $\pi/2$-flux per plaquette
Kang  _et al._ (2019). The tight-binding Hamlitonian in momentum space can be
written as

$\displaystyle H_{\text{quad}}^{(\pi/2)}(\boldsymbol{k})=$
$\displaystyle\gamma_{x}\Gamma_{4}+\gamma_{y}(i\nu_{2}\otimes\nu_{2})+\delta\Gamma_{0}$
(22)
$\displaystyle-\lambda_{x}(\cos(k_{x})(\sigma_{y}\otimes\sigma_{0})+\sin(k_{x})(\sigma_{x}\otimes\sigma_{z}))$
$\displaystyle-\lambda_{y}(\cos(k_{y})(\nu_{1}\otimes\nu_{2})+\sin(k_{y})(i\nu_{1}\otimes\nu_{1})),$

where $\nu_{1}$ and $\nu_{2}$ are defined as

$\nu_{1}=\begin{bmatrix}0&1\\\
i&0\end{bmatrix}\quad\text{and}\quad\nu_{2}=\begin{bmatrix}0&1\\\
-i&0\end{bmatrix},$ (23)

and all $\Gamma$ and $\sigma$ matrices are the same as in Eq. (20).

Model parameters | Phase factors | Electric moments | Sum rule
---|---|---|---
$H_{\text{quad}}^{(\pi/2)}(\gamma_{x},\gamma_{y},\lambda_{x},\lambda_{y},\delta)$ | $\phi_{pp}$ | $\phi_{po}$ | $\phi_{op}$ | $\phi_{oo}$ | $Q^{(2)}_{c}$ | $\tilde{P}_{x}^{\textrm{edge}}$ | $\tilde{P}_{y}^{\textrm{edge}}$ | $\sum(-1)^{ab}\phi_{ab}$
$(0.1,0.1,1.0,1.0,10^{-7})$ | 0.50000(0) | 0.50000(0) | 0.50000(0) | 0.50000(0) | 0.50000(0) | 0.50000(0) | 0.50000(0) | 0.00000(0)
$(0.5,0.5,1.0,1.0,0.7)$ | -0.06272(0) | -0.09622(0) | -0.09622(0) | -0.12980(0) | -0.12975(0) | -0.09675(0) | -0.09675(0) | -0.00009(0)
$(0.2,0.3,1.0,1.0,0.1)$ | -0.40237(0) | -0.42169(0) | -0.42104(0) | -0.44037(0) | -0.44039(0) | -0.42192(0) | -0.42155(0) | -0.00001(0)
$(0.1,0.1,1.0,1.2,0.1)$ | -0.42101(0) | -0.43853(0) | -0.43160(0) | -0.44913(0) | -0.44913(0) | -0.43855(0) | -0.43165(0) | -0.00001(0)
Table 3:  The phase factors $\phi_{ab}$ and the electric moments
$(Q_{c}^{(2)},\tilde{P}_{x}^{\text{edge}},\tilde{P}_{y}^{\text{edge}})$ are
computed for the quadrupole insulator with $\pi/2$-flux per plaquette Eq. (22)
with various parameters. Here we use the same extrapolation procedure as in
Table. 1. In this case as well, the edge-localized polarizations
$\tilde{P}_{x}^{\textrm{edge}}$ and $\tilde{P}_{y}^{\textrm{edge}}$ agree with
phase factors $\phi_{po}$ and $\phi_{op}$, corner charge $Q^{(2)}_{c}$ agrees
with $\phi_{oo}$, and the sum rule is satisfied up to small errors.

Similar to the quadrupole insulator, which has $\pi$-flux per plaquette, the
bulk quadrupole moment of the half-filled ground state of Eq. (22) is
quantized when $\delta=0$. However, when $\delta\neq 0$, the bulk quadrupole
moment $Q_{xy}$ is non-zero and
$(Q_{xy},\tilde{P}_{x/y}^{\textrm{edge}},Q^{(2)}_{c})$ are all distinct unlike
in the quadrupole insulator. These features provide non-trivial checks of the
sum rule Eq. (14) and our observations Eqs. (15) and (19), which are
summarized in Table 3.

####  III.1.4 Two band insulator $H_{two}$

Our final model is the two band insulator introduced in Ref. Zhou  _et al._ ,
2015. The tight-binding Hamlitonian in momentum space can be written as

$\displaystyle
H_{\text{two}}(\boldsymbol{k})=\frac{t_{1}}{2}\sigma_{1}+e^{ik_{x}}\begin{bmatrix}t_{4}&t_{2}\\\
0&t_{4}\end{bmatrix}+e^{ik_{y}}\begin{bmatrix}t_{5}&t_{3}\\\
0&t_{5}\end{bmatrix}+h.c.,$ (24)

with the five hopping parameters $(t_{1},t_{2},t_{3},t_{4},t_{5})$. The half
filled ground state of this Hamiltonian has non-zero corner charge and edge-
localized polarizations.

The numerical results for the two band insulator with the same parameters as
in Ref. Zhou  _et al._ , 2015 are summarized in Table 4.

Model parameters |  Phase factors ($\times 10^{-6}$) |  Electric moments ($\times 10^{-6}$) |  Sum rule ($\times 10^{-6}$)
---|---|---|---
$H_{\text{two}}(t_{1},t_{2},t_{3},t_{4},t_{5})$ | $\phi_{pp}$ | $\phi_{po}$ | $\phi_{op}$ | $\phi_{oo}$ | $Q^{(2)}_{c}$ | $\tilde{P}_{x}^{\textrm{edge}}$ | $\tilde{P}_{y}^{\textrm{edge}}$ | $\sum(-1)^{ab}\phi_{ab}$
$(-2.2,-0.15,-0.1,-0.09,-0.06)$ | 0.000(0) | 1.459(0) | 1.459(0) | 2.920(0) | 2.918(0) | 0.851(0) | 0.851(0) | 0.002(0)
Table 4:  The phase factors $\phi_{ab}$ and the electric moments
$(Q_{c}^{(2)},\tilde{P}_{x}^{\text{edge}},\tilde{P}_{y}^{\text{edge}})$ are
computed for the two band insulator $H_{\text{two}}$ Eq. (24). Here we use the
same extrapolation procedure as in Table. 1. In this case as well, the edge-
localized polarizations $\tilde{P}_{x}^{\textrm{edge}}$ and
$\tilde{P}_{y}^{\textrm{edge}}$ agree with phase factors $\phi_{po}$ and
$\phi_{op}$, corner charge $Q^{(2)}_{c}$ agrees with $\phi_{oo}$, and the sum
rule is satisfied up to small errors.

###  III.2 Models with the hybrid Wannier value at $0.5$

We now turn to models with the hybrid Wannier value at $0.5$, which are all
based on $C_{4}$-symmetric insulators Benalcazar  _et al._ (2019). Note that
for the models presented below, a HWF with the corresponding hybrid Wannier
value at $0.5$ exists even after introducing a $C_{4}$ symmetry breaking term.
In these models, corner charges are originated from the filling anomaly $\eta$
Benalcazar  _et al._ (2019):

$\eta=\textrm{\\# ions}-\textrm{\\# electrons}\quad\textrm{mod}\,\,4,$ (25)

where the number of ions equals the number of lattice sites times the filling.
In each model, we take into account the filling anomaly $\eta$ which
determines the corner charge and the phase factor $\phi_{oo}$, where we tune
the chemical potential to include or not include the zero modes in the full
open boundary conditions.

For tested models which are specified below, we numerically compute the phase
factors $\phi_{ab}$ Eq. (13) and the corner charge $Q_{c}^{(2)}$ Eq. (11). The
results are summarized in Tables 5 and 6. In the tables,
$\tilde{P}_{x/y}^{\textrm{edge}}$ are not presented since they depend
sensitively on the choice of the branch cut for the considered models, as we
discuss in Appendix D. However, we do compare the corner charge $Q^{(2)}_{c}$
with $\phi_{oo}$ and check whether the sum rule Eq. (14) is satisfied. All
observables in the tables are the extrapolated values in the thermodynamic
limit as described in Sec. III.1.

For each model, we separate the cases based on the quadrupole moment
$\phi_{0}$ of the background ions,

$\phi_{0}\equiv-n_{e}(L_{x}+1)(L_{y}+1)/4\quad\text{(mod 1)},$ (26)

with $n_{e}$ being the charge of the ions at each lattice site and
$(L_{x},L_{y})$ being the linear system sizes. We then obtain $\phi_{ab}$ and
$Q_{c}^{(2)}$ in the thermodynamic limit $L\rightarrow\infty$ for each value
of $\phi_{0}$ and summarize the numerical results in Tables 5 and 6. We find
that for the model Eq. (28) defined below, the convergent values of the three
phase factors $\phi_{pp}$, $\phi_{op}$, and $\phi_{po}$ depend on $\phi_{0}$
while the phase factor $\phi_{oo}$ is convergent on the same model. One the
contrary, all the phase factors are convergent on the other model Eq. (27).
Thus, at least $\phi_{oo}$ is interpretable among the four phase factors for
these models. In addition, we find that in all cases, $\phi_{oo}$ converges to
corner charge $Q^{(2)}_{c}$ up to an error $O(10^{-5})$, and the sum rule
holds up to the same order of an error.

####  III.2.1 $C_{4}$-symmetric insulator $H_{2b}^{(4)}$

We first consider a $C_{4}$-symmetric insulator $H_{2b}^{(4)}$ introduced in
Ref. Benalcazar  _et al._ , 2019 with a $C_{4}$ symmetry breaking term
parameterized by $\delta$:

$H^{(4)}_{2b}(\boldsymbol{k})=\begin{bmatrix}\delta&t&e^{i(k_{x}+k_{y})}&t\\\
t&-\delta&t&e^{-i(k_{x}-k_{y})}\\\ e^{-i(k_{x}+k_{y})}&t&\delta&t\\\
t&e^{i(k_{x}-k_{y})}&t&-\delta\end{bmatrix}.$ (27)

When $\delta=0$ and $t<1$, the half-filled ground state of Eq. (27) respects
$C_{4}$ symmetry and has filling anomaly $\eta=2$. Therefore on full open
boundary conditions, filling $2N-2$ electrons, where $N$ is the number of
sites, results in the corner charge of $-1/2$ (in unit of electron charge) at
each corner. When $\delta\neq 0$, the ground state no longer respects $C_{4}$
symmetry and the corner charge is not quantized.

The numerical results for the $C_{4}$-symmetric insulator $H_{2b}^{(4)}$ are
summarized in Table 5.

Model parameters | Background | Phase factors | Electric moment | Sum rule
---|---|---|---|---
$H^{(4)}_{2b}(t,\delta)$ | $\phi_{0}$ | $\phi_{pp}$ | $\phi_{po}$ | $\phi_{op}$ | $\phi_{oo}$ | $Q^{(2)}_{c}$ | $\sum(-1)^{ab}\phi_{ab}$
$(0.1,0.001)$ | 0.5 | 0.50000(0) | 0.49992(0) | 0.49992(0) | 0.49983(0) | 0.49983(0) | 0.00000(0)
0
$(0.1,0.1)$ | 0.5 | 0.49972(0) | 0.49433(0) | 0.49433(0) | 0.48896(0) | 0.48895(0) | 0.00002(0)
0 | 0.49976(0) | 0.49435(0) | 0.49435(0)
$(0.1,0.5)$ | 0.5 | 0.49732(2) | 0.49095(1) | 0.49095(1) | 0.48461(0) | 0.48460(0) | 0.00003(0)
0 | 0.49732(0) | 0.49095(0) | 0.49095(0) | 0.00002(0)
Table 5:  The phase factors $\phi_{ab}$ and the corner charge $Q_{c}^{(2)}$
are computed for $H^{(4)}_{2b}$ in Eq. (27) with various parameters. The
number of filled electrons is taken as $2N-2$ with the number of sites $N$
when the system is under the full open boundary condition due to the filling
anomaly $\eta=2$. Here we use the same extrapolation procedure as in Table. 1.
The corner charge agrees with $\phi_{oo}$ and the sum rule Eq. (14) is
satisfied up to small errors.

####  III.2.2 $C_{4}$-symmetric insulator $H^{(4)}_{1b}\oplus H^{(4)}_{2c}$

Finally, we consider the $C_{4}$-symmetric insulator Benalcazar  _et al._
(2019) $H^{(4)}_{1b}\oplus H^{(4)}_{2c}$. This model was used as an example on
which $\phi_{pp}$ depends on the system size Ono  _et al._ (2019), which
indicates it is ill-defined on this model. We also see that $\phi_{pp}$ in
fact does as discussed below. The tight-binding Hamiltonian of this model in
the momentum space is given by

$H^{(4)}_{1b}\oplus
H^{(4)}_{2c}(\boldsymbol{k})=\begin{bmatrix}H^{(4)}_{1b}&\gamma^{(4)}(t)\\\
\gamma^{(4)}(t)^{\dagger}&H^{(4)}_{2c}\end{bmatrix},$ (28)

where

$H^{(4)}_{1b}(\boldsymbol{k})=\begin{bmatrix}0&t_{x}+e^{ik_{x}}&0&t_{y}+e^{ik_{y}}\\\
t_{x}+e^{-ik_{x}}&0&t_{y}+e^{ik_{y}}&0\\\
0&t_{y}+e^{-ik_{y}}&0&t_{x}+e^{-ik_{x}}\\\
t_{y}+e^{-ik_{y}}&0&t_{x}+e^{k_{x}}&0\end{bmatrix},$ (29)
$H^{(4)}_{2c}(\boldsymbol{k})=\begin{bmatrix}0&0&t_{x}&0\\\ 0&0&0&t_{y}\\\
t_{x}&0&0&0\\\ 0&t_{y}&0&0\end{bmatrix}+1.5\begin{bmatrix}0&0&e^{ik_{x}}&0\\\
0&0&0&e^{ik_{y}}\\\ e^{-ik_{x}}&0&0&0\\\ 0&e^{-ik_{y}}&0&0\end{bmatrix},$ (30)

and

$\gamma^{(4)}(t)=\begin{bmatrix}t&t&0&0\\\ 0&t&t&0\\\ 0&0&t&t\\\
t&0&0&t\end{bmatrix}.$ (31)

Both the $1/4$-filled ground state of $H^{(4)}_{1b}$ and the $1/2$-filled
ground state of $H^{(4)}_{2c}$ have non-zero bulk polarizations, and we stack
two models and introduce additional hopping term given by $\gamma^{(4)}(t)$ so
that the $3/8$-filled ground state of the stacked model has vanishing bulk
polarization while respecting $C_{4}$ symmetry. The stacked model
$H^{(4)}_{1b}\oplus H^{(4)}_{2c}$ has the filling anomaly $\eta=3$ which
results in $-3/4$ corner charge (in the unit of electron charge) at each
corner when $4N-3$ electrons are filled under full open boundary conditions,
where $N$ is the number of sites.

The numerical results for the $C_{4}$-symmetric insulator $H^{(4)}_{1b}\oplus
H^{(4)}_{2c}$ are summarized in Table 6.

Model parameters | Background | Phase factors | Electric moment | Sum rule
---|---|---|---|---
$H^{(4)}_{1b}\oplus H^{(4)}_{2c}(t_{x},t_{y},t)$ | $\phi_{0}$ | $\phi_{pp}$ | $\phi_{po}$ | $\phi_{op}$ | $\phi_{oo}$ | $Q^{(2)}_{c}$ | $\sum(-1)^{ab}\phi_{ab}$
$(0,0,0.1)$ | -0.25 | 0.33354(0) | -0.33994(1) | -0.07653(1) | 0.24992(1) | 0.25000(0) | -0.00007(1)
0.5 | 0.14434(1) | -0.49995(1) | -0.10573(0) | 0.24997(0) | -0.00002(0)
0.25 | 0.43112(2) | -0.15944(1) | -0.15944(1) | 0.24999(0) | 0.00000(0)
0 | -0.24983(2) | -0.49992(1) | -0.49992(1) | 0.24999(0) | 0.00000(0)
$(0.1,0.1,0.1)$ | -0.25 | 0.374(2) | -0.307(1) | -0.0683(3) | 0.24992(1) | 0.25000(0) | -0.00007(1)
0.5 | 0.1592(5) | 0.49999(0) | -0.0908(5) | 0.24997(0) | -0.00002(0)
0.25 | 0.489(2) | -0.131(1) | -0.131(1) | 0.24999(0) | 0.00000(0)
0 | -0.24984(2) | -0.49992(1) | -0.49992(1) | 0.24999(0) | 0.00000(0)
$(0.01,0.1,0.1)$ | -0.25 | 0.33455(4) | -0.33916(4) | -0.07630(0) | 0.24992(1) | 0.25000(0) | -0.00007(1)
0.5 | 0.14470(1) | 0.49999(0) | -0.10530(1) | 0.24997(0) | -0.00002(0)
0.25 | 0.43274(8) | -0.15863(4) | -0.15863(4) | 0.24999(0) | 0.00000(0)
0 | -0.24983(2) | -0.49992(1) | -0.49992(1) | 0.24999(0) | 0.00000(0)
Table 6:  The phase factors $\phi_{ab}$ and the corner charge $Q_{c}^{(2)}$
are computed for $H^{(4)}_{1b}\oplus H^{(4)}_{2c}$ in Eq. (28) with various
parameters. On each $\phi_{0}$, all observables are obtained by extrapolating
them as a function of $1/\sqrt{L_{x}L_{y}}$. The phase factors $\phi_{ab}$
converge to different values on each $\phi_{0}$, so we separate the
convergence values of $\phi_{ab}$ by each value of $\phi_{0}$. Each
extrapolation uses a sequence of four system sizes with the interval $(\Delta
L_{x},\Delta L_{y})=(4,4)$, and each sequence starts from
$(L_{x},L_{y})=(18,16),(18,17),(18,18),$ and $(19,19)$, which correspond to
$\phi_{0}=-0.25,0.5,0.25$, and $0$, respectively. Here, we use the same
extrapolation procedure as in Table. 1. The corner charge agrees with
$\phi_{oo}$ and the sum rule Eq. (14) is satisfied up to small errors.

##  IV Conclusion

In this paper, we have presented the bulk-boundary correspondence for the bulk
quadrupole moment which is expressed in terms of the phase factors of the
expectation values of the many-body operators with respect to the ground
states under various boundary conditions. Our bulk-boundary correspondence is
given by the cancellation between the four gauge-invariant phase factors,
which can be computed in fully interacting systems. We also have proposed that
one of them is expected to be identified with the corner charge. On these, we
have numerically observed that when the band insulator without the hybrid
Wannier value at $0.5$, each phase factor corresponds to a physical observable
including the bulk multipole moment and edge-localized polarizations. Whether
the same correspondence holds with the hybrid Wannier value at $0.5$ has not
been fully understood yet. Through extensive numerical computations on band
insulators, we have found that our bulk-boundary correspondence and the
identification of one of the phase factors with the corner charge hold up to
small errors of $O(10^{-4})$, which might be originated from the finite-size
effect of the numerics. Furthermore, we have also found that another two of
the phase factors have values similar to the edge-localized polarizations
defined in Ref. Benalcazar  _et al._ , 2017b with differences at most
$O(10^{-3})$, which supports our observations.

Let us conclude by making remarks on possible future directions of our work.
It would be interesting to prove the sum rules at least in the case of band
insulators. This would also extend our knowledge of bulk multipole moments in
solids. Since our formulation works for disordered and interacting systems,
one could also test our sum rules for those cases. Finally, finding
connections between previous works Trifunovic (2020); Ren  _et al._ (2021);
Ying  and our sum rule would be interesting.

###### Acknowledgements.

WL and GYC acknowledge the support of the National Research Foundation of
Korea (NRF) funded by the Korean Government (Grant No. 2020R1C1C1006048 and
2020R1A4A3079707), as well as Grant No. IBS-R014-D1. This work is also
supported by the Air Force Office of Scientific Research under Award No.
FA2386-20-1-4029. BK is supported by KIAS individual Grant PG069402 at Korea
Institute for Advanced Study and the National Research Foundation of Korea
(NRF) grant funded by the Korea government (MSIT) (No. 2020R1F1A1075569). GYC
acknowledges financial support from Samsung Science and Technology Foundation
under Project Number SSTF-BA2002-05.

##  Appendix A Classical Multipolar Sum Rules

In this appendix, we derive classical dipolar and quadrupolar sum rules based
on classical electrostatics. The derivations closely follow the construction
presented in Ref. Benalcazar  _et al._ , 2017b. Here, we consider continuum
systems for classical systems.

###  A.1 Classical Dipolar Sum Rule

The macroscopic polarization $\vec{\mathcal{P}}(\vec{R})$ is defined in the
classical electrostatics as the averaged first moment of the charge density
over a region $v(\vec{R})$, which is small compared with the whole system,

$\vec{\mathcal{P}}(\vec{R})=\frac{1}{|v(\vec{R})|}\int_{v(\vec{R})}d^{3}r\rho(\vec{r}+\vec{R})\vec{r},$
(32)

where $\rho(\vec{r})$ is the volume charge density and $|v(\vec{R})|$ is the
volume of the region $v(\vec{R})$. It is well-known that a macroscopic system
with the bulk polarization $\vec{\mathcal{P}}$ induces the charge density
$\rho$ via

$\rho=-\gradient\cdot\vec{\mathcal{P}}.$ (33)

Suppose that the system has a boundary to the vacuum so that the polarization
drops on the boundary, then the boundary charge $\mathcal{Q}_{c}$ is
accumulated on the boundary:

$\mathcal{Q}_{c}=-\int_{v}\gradient\cdot\vec{\mathcal{P}}dv=-\oint_{\partial
v}\vec{\mathcal{P}}\cdot d\vec{s},$ (34)

where $\partial v$ is the boundary of a volume $v$ which encircles the
boundary of the system as shown in Fig. 2, and $\vec{s}$ is the surface vector
normal to $\partial v$. If the system is one-dimensional and has a uniform
bulk polarization density $\mathcal{P}$ along the $x$-direction, then Eq. (34)
becomes

$\mathcal{Q}_{c}=\mathcal{P},$ (35)

which is depicted in Fig. 2, where the circle containing the plus symbol being
the boundary charge.

Figure 2: A schematics of the polarization in a one-dimensional system. Here,
gray arrows represent the bulk polarization. The bulk polarization $P$
accumulates charges at each boundary with the same amount but different signs.

###  A.2 Classical Quadrupolar Sum Rule

The macroscopic quadrupole moment $\mathcal{Q}_{ij}(\vec{R})$ is defined in
the classical electrostatics as the averaged second moment of the charge
density over $v(\vec{R})$,

$\mathcal{Q}_{ij}(\vec{R})=\frac{1}{|v(\vec{R})|}\int_{v(\vec{R})}d^{3}r\rho(\vec{r}+\vec{R})r_{i}r_{j},$
(36)

where $\rho(\vec{r})$ is the volume charge density. We first consider a
classical system with a uniform bulk quadrupole moment $\mathcal{Q}_{ij}$. We
assume that the system has a vanishing bulk polarization, which follows when
the system has the inversion symmetry. Then the quadrupole moment
$\mathcal{Q}_{ij}$ induces the charge density $\rho$ via

$\rho=\frac{1}{2}\sum_{i,j}\partial_{i}\partial_{j}\mathcal{Q}_{ij}.$ (37)

If we consider a rectangular geometry as in Fig. 3, then the corner charge
$\mathcal{Q}_{c}$ accumulated on a corner of a macroscopic two-dimensional
system due to quadrupole moment $\mathcal{Q}_{ij}$ drop is given by

$\mathcal{Q}_{c}=\frac{1}{2}\sum_{i,j}\int_{v}\partial_{i}\partial_{j}\mathcal{Q}_{ij}dv=\frac{1}{2}\sum_{i,j}\oint_{\partial
v}\partial_{i}\mathcal{Q}_{ij}n_{j}ds,$ (38)

where $n_{i}$ is the $i$-th component of the unit vector normal to $\partial
v$, $v$ is a volume envelops the corner and $\partial v$ is the boundary of
$v$, which are illustrated in Fig. 3. For a rectangular system, Eq. (38)
becomes

$\mathcal{Q}_{c}=\mathcal{Q}_{xy},$ (39)

where we used the symmetry of the quadrupole moment,
$\mathcal{Q}_{xy}=\mathcal{Q}_{yx}$. In addition to inducing corner charge,
the quadrupole moment also induces the polarizations along the edges, which we
call the edge-localized polarizations. To this end, we first consider the
geometry given in Fig. 3. The line charge density $\sigma_{w}$ at the edge $w$
in Fig. 3 is given as

$\sigma_{w}=-\sum_{i,j}\partial_{j}(n_{i}\mathcal{Q}_{ij}),$ (40)

where $n_{i}$ is the normal vector of the edge in Fig. 3. Since the charge
density is given as the divergence of $\sum_{i}n_{i}\mathcal{Q}_{ij}$, this
can be considered as the polarization induced by the quadrupole moment
localized on the edge $w$. We will call this polarization as
$\mathcal{P}^{\text{quad}}_{n}$. For a rectangular system, the quadrupole
induced polarizations at the $y$-boundaries $\mathcal{P}^{\text{quad}}_{x}$
and the $x$-boundaries $\mathcal{P}^{\text{quad}}_{y}$ are given by

$\mathcal{P}^{\text{quad}}_{x}=\mathcal{P}^{\text{quad}}_{y}=\mathcal{Q}_{xy}.$
(41)

On the other hand, one can dress one dimensional systems having edge
polarizations $\mathcal{P}^{\text{free}}_{x}$ and
$\mathcal{P}^{\text{free}}_{y}$ along the boundaries while respecting the
inversion symmetry, where $\mathcal{P}^{\text{free}}_{x}$ and
$\mathcal{P}^{\text{free}}_{y}$ do not come from the quadrupole moment. Note
that $\mathcal{P}^{\text{free}}_{x/y}$ also accumulates the charge
$\mathcal{Q}^{\textrm{free}}_{c}$ at the corner due to the dipolar sum rule

$\mathcal{Q}^{\textrm{free}}_{c}=\mathcal{P}^{\text{free}}_{x}+\mathcal{P}^{\text{free}}_{y}.$
(42)

Thus, if we define the edge-localized polarization
$\mathcal{P}^{\text{edge}}_{x/y}$ in classical systems as

$\mathcal{P}^{\text{edge}}_{x/y}\equiv\mathcal{P}^{\text{quad}}_{x/y}+\mathcal{P}^{\text{free}}_{x/y},$
(43)

then the corner charge is given by

$\mathcal{Q}_{c}=-\mathcal{Q}_{xy}+\mathcal{P}^{\text{edge}}_{x}+\mathcal{P}^{\text{edge}}_{y}.$
(44)

We call Eq. (44) the quadrupolar sum rule in classical systems, which is
summarized in Fig. 3.

Figure 3: Schematic understanding of the bulk quadrupole moment and the free
edge polarization. Here, gray arrows and orange arrows represent the
quadrupole moment and the free edge polarization, respectively. Both of them
accumulate charges at each corner, and the total amount of charge at each
corner is fixed by the charge conservation law Eq. (38).

##  Appendix B $\hat{U}_{1}$ operator and the corner charge

In this appendix, we show that for band insulators in one dimension, the phase
factor of the $\hat{U}_{1}$ expectation value with respect to the ground state
in the open boundary condition converges to the boundary charge in the
thermodynamic limit.

We prove the claim by constructing bulk and boundary Wannier functions. The
bulk Wannier functions are identical to the Wannier functions in the periodic
boundary condition and the boundary Wannier functions are localized near the
left and right boundaries. The phase factor of $\hat{U}_{1}$ expectation value
is determined by the bulk Wannier functions while boundary Wannier functions
do not contribute, as in the original bulk-boundary correspondence for the
polarization Vanderbilt and King-Smith (1993); Rhim  _et al._ (2017). To avoid
any ambiguity, we assume that no mode exists at the chemical potential, which
can always be done by either introducing a (small) symmetry breaking term
splitting the zero modes or tuning the chemical potential slightly.

In the open boundary condition, the position operator $\hat{x}$ is well-
defined, which we choose to take values in $\\{1,2,\cdots,L\\}$, and thus we
can diagonalize the position operator in the subspace of occupied single-
particle orbitals $P_{\textrm{occ}}$. We would like to show that the
eigenvectors of $P_{\textrm{occ}}(\hat{x})P_{\textrm{occ}}$ are exactly the
bulk and boundary Wannier functions with desirable properties.

We first recall the results in Refs. Kivelson, 1982; Alexandradinata  _et al._
, 2014b, which state that the Wannier functions $\\{|w_{\alpha,R}\rangle\\}$,
where $\alpha=1,\cdots,n_{e}$ with $n_{e}$ being the electron filling and
$R\in\mathbb{Z}$ labels the position at which the Wannier function is
localized, are the eigenvectors of $\hat{x}$ operator in the subspace of the
occupied single-particle orbitals of the infinite open system with the
corresponding eigenvalues $R+\nu_{\alpha}$. At the same time,
$\\{|w_{\alpha,R}\rangle\\}$ are the eigenvectors of $e^{\frac{2\pi
i\hat{x}}{L}}$ in the subspace of occupied single-particle orbitals of a
periodic system with the corresponding eigenvalues
$e^{i2\pi(R+\nu_{\alpha})/L}$, when the system size $L$ is sufficiently large.
Moreover, the Wannier functions $\\{|w_{\alpha,R}\rangle\\}$ are exponentially
localized and satisfy the usual property $\langle
x|w_{\alpha,R+1}\rangle=\langle x-1|w_{\alpha,R}\rangle$.

Let us take a positive integer $\delta$ which is much larger than the
localization lengths of $\\{|w_{\alpha,R}\rangle\\}$. Using the exponentially
localized nature,
$\\{|w_{\alpha,R}\rangle\\}_{R=\delta+1,\delta_{2},\cdots,L-\delta}$ are the
claimed bulk Wannier functions, i.e., eigenvectors of
$P_{\textrm{occ}}(\hat{x})P_{\textrm{occ}}$ localized on the bulk, where
$P_{\textrm{occ}}$ projects to the subspace of occupied single-particle
orbitals of a finite open system. Also, for a sufficiently large $L$, the bulk
Wannier functions
$\\{|w_{\alpha,R}\rangle\\}_{R=\delta+1,\delta_{2},\cdots,L-\delta}$ are the
eigenvectors of $P_{\textrm{occ}}(e^{\frac{2\pi
i\hat{x}}{L}})P_{\textrm{occ}}$ with the corresponding eigenvalues
$e^{i2\pi(R+\nu_{\alpha})/L}$. The remaining eigenvectors of
$P_{\textrm{occ}}(\hat{x})P_{\textrm{occ}}$ are localized near the left and
right boundaries, hence we call them the boundary Wannier functions.

We now evaluate the expectation value of $\hat{U}_{1}$ using the bulk and
boundary Wannier functions. Since the ground state of a band insulator is
given by the slater determinant of occupied single-particle orbitals, which we
choose to be the bulk and boundary Wannier functions,
$\langle\hat{U}_{1}\rangle$ can be expressed as

$\langle\hat{U}_{1}\rangle=e^{-(L+1)n_{e}\pi
i}\det\left[\left(\begin{array}[]{ccc}\mathcal{L}&0&0\\\ 0&\mathcal{D}&0\\\
0&0&\mathcal{R}\end{array}\right)\right],$ (45)

where $e^{-(L+1)n_{e}\pi i}$ is the phase factor from background ions,
$\mathcal{L}$ and $\mathcal{R}$ are $n_{e}\delta$-by-$n_{e}\delta$ matrices,
and

$\mathcal{D}=\textrm{diag}(e^{2\pi i\frac{\delta+1+\nu_{1}}{L}},\cdots,e^{2\pi
i\frac{L-\delta+\nu_{N_{\textrm{oc}}}}{L}})$ (46)

is a diagnal matrix. By expanding the diagonal part in the determinant, we get

$\displaystyle\langle U_{1}\rangle$ $\displaystyle=e^{-(L+1)n_{e}\pi
i}e^{\frac{2\pi
i}{L}\big{(}\frac{(L+1)(L-2\delta)}{2}n_{e}+L\sum_{\alpha=1}^{n_{e}}\nu_{\alpha}\big{)}}$
$\displaystyle\quad\times\det(\mathcal{L})\det(\mathcal{R})$
$\displaystyle=e^{2\pi i\frac{L+1}{L}(-\delta)}e^{2\pi
i\sum_{\alpha=1}^{n_{e}}\nu_{\alpha}}\det(\mathcal{L})\det(\mathcal{R})$
$\displaystyle=e^{-2\pi i\frac{\delta}{L}}e^{2\pi
i\sum_{\alpha=1}^{n_{e}}\nu_{\alpha}}\det(\mathcal{L})\det(\mathcal{R}),$ (47)

where we used the fact that $\delta$ is an integer in the last equality. Note
that $\mathcal{L}$ ($\mathcal{R}$) is the matrix representation of
$e^{\frac{2\pi i\hat{x}}{L}}$ with respect to the left (right) Wannier
functions which are localized within the range $\delta$. Thus, in the
thermodynamic limit $\hat{U}_{1}=I+O(\frac{1}{L})$ on the left (right) Wannier
functions and hence $\det(\mathcal{L})=1+O(\frac{1}{L})$
($\det(\mathcal{R})=1+O(\frac{1}{L})$). Thefore,

$\lim_{L\to\infty}\frac{1}{2\pi}\textrm{Im}\Big{[}\log\big(\langle
U_{1}\rangle\big{missing})\Big{]}=\sum_{\alpha=1}^{N_{\textrm{occ}}}\nu_{\alpha}=Q_{c}^{(1)},$
(48)

where the last equality follows from the original bulk-boundary correspondence
Vanderbilt and King-Smith (1993); Rhim  _et al._ (2017). This completes the
proof that $\hat{U}_{1}$ expectation value reproduces the boundary charge for
band insulators in the open boundary condition.

##  Appendix C Sum rule in $C_{3}$-symmetric insulator

As a demonstration that the sum rule Eq. (14) works for insulators having
other types of symmetry besides $C_{4}$ symmetry, we consider a
$C_{3}$-symmetric insulator Benalcazar  _et al._ (2019) and numerically
confirm the sum rule on this model.

###  C.1 $C_{3}$-symmetric insulator: $H^{(3)}_{2b}\oplus H^{(3)}_{2c}$

The tight-binding hamiltonian in the momentum space of the $C_{3}$-symmetric
insulator $H^{(3)}_{2b}\oplus H^{(3)}_{2c}$ is given by Benalcazar  _et al._
(2019)

$H^{(3)}_{2b}\oplus
H^{(3)}_{2c}(\boldsymbol{k})=\begin{bmatrix}H^{(3)}_{2b}&\gamma^{(3)}(t)\\\
\gamma^{(3)}(t)^{\dagger}&H^{(3)}_{2c}\end{bmatrix},$ (49)

where

$H^{(3)}_{2b}(\boldsymbol{k})=\begin{bmatrix}0&t_{0}+e^{i\boldsymbol{k}\cdot\boldsymbol{a_{2}}}&t_{0}+e^{-i\boldsymbol{k}\cdot\boldsymbol{a_{3}}}\\\
t_{0}+e^{-i\boldsymbol{k}\cdot\boldsymbol{a_{2}}}&0&t_{0}+e^{-i\boldsymbol{k}\cdot\boldsymbol{a_{1}}}\\\
t_{0}+e^{i\boldsymbol{k}\cdot\boldsymbol{a_{3}}}&t_{0}+e^{i\boldsymbol{k}\cdot\boldsymbol{a_{1}}}&0\end{bmatrix},$
(50)
$H^{(3)}_{2c}(\boldsymbol{k})=\begin{bmatrix}0&t_{0}+e^{i\boldsymbol{k}\cdot\boldsymbol{a_{1}}}&t_{0}+e^{i\boldsymbol{k}\cdot\boldsymbol{a_{2}}}\\\
t_{0}+e^{-i\boldsymbol{k}\cdot\boldsymbol{a_{1}}}&0&t_{0}+e^{-i\boldsymbol{k}\cdot\boldsymbol{a_{3}}}\\\
t_{0}+e^{-i\boldsymbol{k}\cdot\boldsymbol{a_{2}}}&t_{0}+e^{i\boldsymbol{k}\cdot\boldsymbol{a_{3}}}&0\end{bmatrix},$
(51)

and

$\gamma^{(3)}(t)=\begin{bmatrix}t&t&0\\\ 0&t&t\\\ t&t&0\end{bmatrix}.$ (52)

The $2/3$-filled ground states of $H^{(3)}_{2b}$ and $H^{(3)}_{2c}$ have non-
vanishing polarizations $1/3$ and $2/3$, respectively, so the stacked model
$H^{(3)}_{2b}\oplus H^{(3)}_{2c}$ at $2/3$-filling with the inter-layer
hopping term $\gamma^{(3)}(t)$ has vanishing bulk polarization and respects
$C_{3}$-symmetry.

We will consider the stacked model Eq. (49) on the square lattice. To do so,
we use the coordinates system $(x,y)$ with two basis vectors
$\vec{a}_{1}=(1,0)$ and $\vec{a}_{2}=(1/2,\sqrt{3}/2)$, i.e., the coordinates
$(x,y)$ refer to as the vector of $x\vec{a}_{1}+y\vec{a}_{2}$. All the phase
factors $\phi_{ab}$ defined in Eq. (13) and appears in Table 7 are computed
using these coordinates. We find that the stacked model on the square lattice
at $2/3$-filling has corner charges $1/3$ and $-1/3$ well localized near the
corners.

The numerical results are summarized in Table 7. We see that the corner charge
$Q^{(2)}_{c}$ agrees with $\phi_{oo}$ and the sum rule Eq. (14) is satisfied
up to small errors of $O(10^{-5})$.

Model parameters | Phase factors | Electric moment | Sum rule
---|---|---|---
$H^{(3)}_{2b}\oplus H^{(3)}_{2c}(t_{0},t)$ | $\phi_{pp}$ | $\phi_{po}$ | $\phi_{op}$ | $\phi_{oo}$ | $Q^{(2)}_{c}$ | $\sum(-1)^{ab}\phi_{ab}$
$(0.001,0.001)$ | -0.33340(0) | -0.33333(0) | -0.33333(0) | -0.33331(0) | -0.33333(0) | -0.00005(0)
$(0.1,0.001)$ | -0.33340(0) | -0.33335(0) | -0.33335(0) | -0.33334(0) | -0.33337(0) | -0.00005(0)
$(0.1,0.1)$ | -0.33377(0) | -0.33441(0) | -0.33441(0) | -0.33510(0) | -0.33513(0) | -0.00004(0)
$(0.2,0.1)$ | -0.33380(0) | -0.33504(0) | -0.33504(0) | -0.33631(0) | -0.33634(0) | -0.00004(0)
Table 7:  The phase factors $\phi_{ab}$ and the corner charge $Q_{c}^{(2)}$
are computed for $H^{(3)}_{2b}\oplus H^{(3)}_{2c}(t_{0},t)$ in Eq. (49) on a
square lattice with various parameters. The details about the lattice and
filling are discussed in the paragraph below Eq. (52). Here we use the same
extrapolation procedure as in Table. 1. The corner charge agrees with
$\phi_{oo}$ and the sum rule Eq. (14) is satisfied up to small errors.

##  Appendix D Branch cut dependence of the HWF-based edge-localized
polarization

In this section, we discuss the branch cut dependence in the hybrid Wannier
values of the HWF-based edge-localized polarization
$\tilde{P}^{\textrm{edge}}_{x/y}$ defined in Ref. Benalcazar  _et al._ , 2017b
and Eq. (18). To be self-contained, we first rewrite the definition of the
edge-localized polarization,

$\tilde{P}^{\textrm{edge}}_{x}=\sum_{j}\sum_{y=1}^{L_{y}/2}\nu^{j}\rho^{j}(y)\quad\text{(mod
1)},$ (53)

where $\rho^{j}(y)$ and $e^{2\pi i\nu_{j}}$ are the density and the $j$-th
eigenvalue of the hybrid Wilson loop $\mathcal{W}_{k_{x}}$ along the
$x$-direction, respectively. We will call $\nu^{j}$ the hybrid Wannier value.

The edge-localized polarization $\tilde{P}^{\textrm{edge}}_{x}$ can be well
defined after one fixes the branch cut of the hybrid Wannier value since the
hybrid Wannier value is defined by a phase angle, which has the modulo $2\pi$
ambiguity. Unless fixing the branch cut, each $\nu^{j}$ can be shifted by an
integer $n^{j}=\pm 1$, which results in shifting
$\tilde{P}^{\textrm{edge}}_{x}$ as

$\Delta\tilde{P}^{\textrm{edge}}_{x}=\sum_{j}\sum_{y=1}^{L_{y}/2}n^{j}\rho^{j}(y).$
(54)

Since $\sum_{y=1}^{L_{y}/2}\rho^{j}(y)$ is not quantized in general,
$\Delta\tilde{P}^{\textrm{edge}}_{x}$ may also not be quantized if $n^{j}$ is
not the same for all $j$. Consequently, the edge-localized polarization is
ambiguous before fixing the branch cut even if we take the modulo one
equivalence on both sides of Eq. (53). Here we fix the branch cut by fixing
the range of the hybrid-Wannier value as $\nu\in(-0.5,0.5]$. This choice is
made for our numerical results to be consistent with that of Ref. Benalcazar
_et al._ , 2017b. Note that another equally valid choice of the range is
$\nu\in[-0.5,0.5)$, which includes $-0.5$ instead of $+0.5$.

However, these possibilities in ranges again cause a problem when a hybrid
Wannier value $\nu^{j}$ is at the branch cut value since $\nu^{j}$ can be
either $+0.5$ or $-0.5$ depending on the choice of the range. In this case,
the edge-localized polarization is invariant under the choice of the range
only when $\sum_{y=1}^{L_{y}/2}\rho^{j}(y)$ is quantized, i.e., the HWF is
localized one of the two half systems. The situation gets worse when more than
two hybrid Wannier values $\nu^{1},\nu^{2},...,\nu^{n}$ are degenerated at the
brach cut since for all possible unitary transformations
$\tilde{\psi}^{i}(y)=\sum_{j=1}^{n}u_{ij}\psi^{j}(y)$ of the Wannier functions
$\\{\psi^{1},\psi^{2},...,\psi^{n}\\}$,
$\sum_{y=1}^{L_{y}/2}\absolutevalue{\tilde{\psi}^{i}(y)}^{2}$ should be
quantized for the edge-localized polarization not to depend sensitively on the
choice of the branch cut. To avoid this difficulty, we computed
$\tilde{P}^{\textrm{edge}}_{x/y}$ only for the models without the hybrid
Wannier value at $0.5$ in the main text.

## References

  * Hasan and Kane (2010) M. Z. Hasan and C. L. Kane, “Colloquium: Topological insulators,” Rev. Mod. Phys. 82, 3045–3067 (2010).
  * Qi and Zhang (2011) Xiao-Liang Qi and Shou-Cheng Zhang, “Topological insulators and superconductors,” Rev. Mod. Phys. 83, 1057–1110 (2011).
  * Benalcazar  _et al._ (2017a) Wladimir A Benalcazar, B Andrei Bernevig, and Taylor L Hughes, “Quantized electric multipole insulators,” Science 357, 61–66 (2017a).
  * Benalcazar  _et al._ (2017b) Wladimir A. Benalcazar, B. Andrei Bernevig, and Taylor L. Hughes, “Electric multipole moments, topological multipole moment pumping, and chiral hinge states in crystalline insulators,” Phys. Rev. B 96, 245115 (2017b).
  * Schindler  _et al._ (2018) Frank Schindler, Ashley M Cook, Maia G Vergniory, Zhijun Wang, Stuart SP Parkin, B Andrei Bernevig, and Titus Neupert, “Higher-order topological insulators,” Science advances 4, eaat0346 (2018).
  * Fu and Kane (2006) Liang Fu and C. L. Kane, “Time reversal polarization and a ${Z}_{2}$ adiabatic spin pump,” Phys. Rev. B 74, 195312 (2006).
  * Coh and Vanderbilt (2009) Sinisa Coh and David Vanderbilt, “Electric polarization in a chern insulator,” Phys. Rev. Lett. 102, 107603 (2009).
  * Hughes  _et al._ (2011) Taylor L. Hughes, Emil Prodan, and B. Andrei Bernevig, “Inversion-symmetric topological insulators,” Phys. Rev. B 83, 245132 (2011).
  * Yu  _et al._ (2011) Rui Yu, Xiao Liang Qi, Andrei Bernevig, Zhong Fang, and Xi Dai, “Equivalent expression of ${Z}_{2}$ topological invariant for band insulators using the non-abelian berry connection,” Phys. Rev. B 84, 075119 (2011).
  * Alexandradinata  _et al._ (2014a) A. Alexandradinata, Xi Dai, and B. Andrei Bernevig, “Wilson-loop characterization of inversion-symmetric topological insulators,” Phys. Rev. B 89, 155114 (2014a).
  * King-Smith and Vanderbilt (1993) R. D. King-Smith and David Vanderbilt, “Theory of polarization of crystalline solids,” Phys. Rev. B 47, 1651–1654 (1993).
  * Resta (1993) R Resta, “Macroscopic electric polarization as a geometric quantum phase,” Europhysics Letters (EPL) 22, 133–138 (1993).
  * Vanderbilt and King-Smith (1993) David Vanderbilt and R. D. King-Smith, “Electric polarization as a bulk quantity and its relation to surface charge,” Phys. Rev. B 48, 4442–4455 (1993).
  * Resta (1998) Raffaele Resta, “Quantum-mechanical position operator in extended systems,” Phys. Rev. Lett. 80, 1800–1803 (1998).
  * Watanabe and Ono (2020) Haruki Watanabe and Seishiro Ono, “Corner charge and bulk multipole moment in periodic systems,” Phys. Rev. B 102, 165120 (2020).
  * Kang  _et al._ (2019) Byungmin Kang, Ken Shiozaki, and Gil Young Cho, “Many-body order parameters for multipoles in solids,” Phys. Rev. B 100, 245134 (2019).
  * Wheeler  _et al._ (2019) William A. Wheeler, Lucas K. Wagner, and Taylor L. Hughes, “Many-body electric multipole operators in extended systems,” Phys. Rev. B 100, 245135 (2019).
  * Ono  _et al._ (2019) Seishiro Ono, Luka Trifunovic, and Haruki Watanabe, “Difficulties in operator-based formulation of the bulk quadrupole moment,” Phys. Rev. B 100, 245133 (2019).
  * Rhim  _et al._ (2017) Jun-Won Rhim, Jan Behrends, and Jens H. Bardarson, “Bulk-boundary correspondence from the intercellular zak phase,” Phys. Rev. B 95, 035421 (2017).
  * Trifunovic (2020) Luka Trifunovic, “Bulk-and-edge to corner correspondence,” Phys. Rev. Research 2, 043012 (2020).
  * Ren  _et al._ (2021) Shang Ren, Ivo Souza, and David Vanderbilt, “Quadrupole moments, edge polarizations, and corner charges in the wannier representation,” Phys. Rev. B 103, 035147 (2021).
  * (22) Xuzhe Ying, “Change of corner charge and adiabatic current distribution in two dimensional insulators with inversion symmetry,” arXiv:2109.10395 .
  * Trifunovic  _et al._ (2019) Luka Trifunovic, Seishiro Ono, and Haruki Watanabe, “Geometric orbital magnetization in adiabatic processes,” Phys. Rev. B 100, 054408 (2019).
  * Zhou  _et al._ (2015) Yuanjun Zhou, Karin M. Rabe, and David Vanderbilt, “Surface polarization and edge charges,” Phys. Rev. B 92, 041102 (2015).
  * Benalcazar  _et al._ (2019) Wladimir A. Benalcazar, Tianhe Li, and Taylor L. Hughes, “Quantization of fractional corner charge in ${C}_{n}$-symmetric higher-order topological crystalline insulators,” Phys. Rev. B 99, 245151 (2019).
  * Kivelson (1982) S. Kivelson, “Wannier functions in one-dimensional disordered systems: Application to fractionally charged solitons,” Phys. Rev. B 26, 4269–4277 (1982).
  * Alexandradinata  _et al._ (2014b) A. Alexandradinata, Xi Dai, and B. Andrei Bernevig, “Wilson-loop characterization of inversion-symmetric topological insulators,” Phys. Rev. B 89, 155114 (2014b).