# Effects of Boundary on Orbital Magnetization for a Bilayer System with
Different Chern Numbers

Si-Si Wang  School of Physics and Materials Science, Guangzhou University,
510006 Guangzhou, China  School of Mathematics and Information Science,
Guangzhou University, 510006 Guangzhou, China  Huangpu Research and Graduate
School of Guangzhou University, 510700 Guangzhou, China  Yan Yu  SKLSM,
Institute of Semiconductors, Chinese Academy of Sciences, P.O. Box 912,
Beijing 100083, China  School of Physical Sciences, University of Chinese
Academy of Sciences, Beijing 100049, China  Ji-Huan Guan  Beijing Academy of
Quantum Information Sciences, Beijing 100193, China  SKLSM, Institute of
Semiconductors, Chinese Academy of Sciences, P.O. Box 912, Beijing 100083,
China  Yi-Ming Dai  School of Physics and Materials Science, Guangzhou
University, 510006 Guangzhou, China  Hui-Hui Wang  School of Physics and
Materials Science, Guangzhou University, 510006 Guangzhou, China  Huangpu
Research and Graduate School of Guangzhou University, 510700 Guangzhou, China
Yan-Yang Zhang<EMAIL_ADDRESS>School of Physics and Materials Science,
Guangzhou University, 510006 Guangzhou, China  Huangpu Research and Graduate
School of Guangzhou University, 510700 Guangzhou, China  School of Mathematics
and Information Science, Guangzhou University, 510006 Guangzhou, China

###### Abstract

The real space formalism of orbital magnetization (OM) is an average of the
local OM over some appropriate region of the system. Previous studies prefer a
bulk average (i.e., without including boundaries). Based on a bilayer model
with an adjustable Chern number at half filling, we numerically investigate
the effects from boundaries on the real space expressions of OM. The size
convergence processes of its three constituent terms $M_{\mathrm{LC}}$,
$M_{\mathrm{IC}}$, $M_{\mathrm{BC}}$ are analysed. The topological term
$M_{\mathrm{BC}}$ makes a nonnegligible contribution from boundaries as a
manifestation of edge states, especially in the case of nonzero Chern numbers.
However, we show that the influence of the boundary on $M_{\mathrm{LC}}$ and
$M_{\mathrm{IC}}$ exactly compensates that on $M_{\mathrm{BC}}$. This
compensation effect leads to the conclusion that the whole sample average is
also a correct algorithm in the thermodynamic limit, which gives the same
value as those from the bulk average and the $k$ space formula. This
clarification will be beneficial to further studies on orbitronics, as well as
the orbital magnetoelectric effects in higher dimensions.

##  I I. Introduction

The magnetization of solids comes from two parts: spin magnetization and
orbital magnetization (OM). For the OM, the relevant theory had been confined
to calculating the change of OM as a bulk quantityMauri1996 ; MauriF1996 ;
Pickard2002 ; Sebastiani2002 . Until 2005, the modern theory of OM in solids
was developedVanderbilt2018 . The OM is then found to be closely related to
topological properties of solidsCeresoli2006 ; Malashevich2010 ; Bianco2011 ;
Bianco2013 ; XiDai2019 ; JPHanke2016 ; Lux2018 ; Helin2020 . Recently it is
proposed as the central quantity of a topological memory
deviceTopologicalMemory . How to tune the OM in solids is an important
objective of orbitronicsGo2017 ; Bhowal2020 ; DGo2021 . Besides the induction
of OM from external fieldsDong2020 ; Xiao2021 ; Hamada2020 ; Sharpe2019 ;
Serlin2020 ; Su2020 ; He2020 ; Huang2021 , another method is to produce a
spontaneous OM, such as, by engineering material structuresCinal2022 . These
potential applications raise demands of predicting materials’ OM based on a
deeper understanding of calculation methods.

Similar to Chern number as the integration of the Berry curvature in the
momentum space ($k$ space), the modern theory also expresses OM as a $k$ space
integral, which can be derived from the semi-classical methodXiao2005 ;
Xiao2010 , the Wannier representationThonhauser2005 ; Resta2005 ;
Thonhauser2011 , or the perturbation theoryShi2007 . This expression consists
of two terms associated with the local circulation (LC) and the itinerant
circulation (IC) motions respectively. With one additional term including the
contribution from Berry curvature (BC), it addresses on the same ground the
normal insulators, Chern insulators, and metalsCeresoli2006 .

Afterwards, there is another formalism of expressing OM and Chern number for a
bounded sample, which define them as local quantities in the real space ($r$
space)Bianco2011 ; Bianco2013 ; Bianco2016 ; Drigo2020 ; Sykes2021 . In the
thermodynamic limit, an average of these local quantities over an appropriate
region gives correct values consistent with the $k$ space formalismBianco2013
; MarrazzoResta2017 . Furthermore, the $r$ space formalism of the orbital
magnetoelectric coupling, i.e., the linear response of OM respect to an
electric field, has also been developedMalashevich2010 ; Olsen2017 ;
Varnava2018 ; Varnava2020 , which paves route towards deep understanding of
axion physics in solidsAxionReview1 ; AxionReview2 . While the $k$ space
theory requires translation symmetry of solids, the $r$ space method is
suitable for a system with arbitrary boundary conditions, or with spatial
inhomogeneity.

It has been recently realized that these $r$ space descriptions are not only
mathematical structures but also with physical meanings. Around the
topological phase transition, numerical simulations show that the local Chern
number plays the role as that of the order parameter in thermodynamic phase
transitions, with familiar properties such as universal scaling and divergence
of correlation lengthCaio2019 ; Ulcakar2020 . In twisted bilayer graphene, a
large but non-uniformly distributed local OM were theoretical
proposedXiDai2019 , which were experimentally detectedOMBilayerGraphene .

Nevertheless in this $r$ space formalism, there is still some vagueness on the
effects from boundaries. For the calculation of Chern number, the average
should be over the bulk region _only_ , and the average over the whole sample
(bulk plus boundaries) is always zero even for a Chern insulatorBianco2011 ;
Bianco2013 . Therefore the contribution from boundaries appears to be
nonnegligible but “unphysical”. For the calculation of OM for a normal
insulator, however, the average can be over the bulk _or_ the whole sample,
with the former a faster convergenceBianco2013 . Another option is to redefine
a boundary-determined term to the OM in the $k$ space formalismBianco2016 .
Back to the $r$ space, the whole sample average of OM for a Chern insulator
seems to be questionable due to the BC term closely related to Chern number,
for which the boundary contribution seems to be unphysical as mentioned above.
As for the $r$ space formalism of the orbital magnetoelectric coupling, it is
suggested that the average should be over the whole sample with the same
footing, even with nontrivial topologyOlsen2017 . Therefore, a detailed
understanding of boundary influence on OM is not very satisfactory so far. Is
it really necessary to or not to include boundaries on OM when one is
calculating in the $r$ space?

In this paper, we revisit these problems of the $r$ space formalism
numerically, by scrutinizing individual contributions to OM from the bulk and
boundary regions, as well as individual contributions from the LC, IC and BC
terms, respectively. Our calculations are based on a bilayer system with an
adjustable Chern number. Compared to a monolayer one, it has better tunability
of properties (symmetry, band structure and Chern number), and it is also a
good starting point for future investigations of higher dimensions, e.g., the
axion orbital effects in three dimensionsOlsen2017 ; Varnava2018 ; Varnava2020
; Wang2019 . We confirm that, the result of a whole sample average is still
consistent with the $k$ space method in the thermodynamic limit even for Chern
insulators. We also find that, contributions from different origins, i.e.,
bulk versus boundary, LC, IC versus BC terms, often carry different signs and
cancel each other to result in a small value of total OM. Therefore in the
materials engineering, a possible means to enhance the OM is to suppress one
of them.

This paper is organized as follows. In Sec.II we present the Hamiltonian. In
Sec. III we lay out the $r$ space method and the $k$ space method of OM. In
Sec. IV and V we study the influence of sample boundary on OM. Some small size
effects are discussed in Section VI. Conclusions are drawn in Sec. VII.

##  II II. The Model

The lattice structure of our model is composed of two square-lattice layers
extending in the $x$-$y$ plane as illustrated in Fig. 1. The first layer (red)
and the second layer (black) are stacked along the $z$ direction. Each layer
is a quantum anomalous Hall (QAH) model with a tunable Chern number and there
exist inter-layer couplings. Experimentally, such topological bilayers can be
realized by using state-of-the-art technologies of micro hetero-
structuresFabrication1 ; Fabrication2 ; Fabrication3 , ultracold
atomicFabricationOL1 ; FabricationOL2 ; FabricationOL3  or photonic
systemsFabricationPhot .

Figure 1:  The left panel denotes the lattice structure of our constructed
bilayer model. The first layer(Red) and the second layer (black) have
different model parameters and the same lattice constant $a$. $t$ is the
inter-layer coupling strength. The right panel is the top view of our model.
The region enclosed by blue lines represents the bulk region. $L_{B}$ denotes
the thickness of the boundaries. Here $L_{B}=4$.

The Hamiltonian of this bilayer system can be generally expressed asWang2019

$\mathcal{H}_{\mathrm{bi}}=\mathcal{H}_{1}+\mathcal{H}_{2}+\mathcal{H}_{c},$
(1)

where $\mathcal{H}_{L}$ ($L=1,2$) is the Hamiltonian for the $L$-th layer, and
$\mathcal{H}_{c}$ represents the coupling between them. Within each layer, we
choose the spin up component of the Bernevig-Hughes-Zhang modelBernevig2006
to construct a spinless QAH model. With one $s$ orbital and one $p$ orbital on
each site, this is one of the minimum models of the QAH state. In the $k$
space, the $L$-th layer Hamiltonian reads

$\mathcal{H}_{L}=\sum_{\bm{k,\alpha\beta}}c^{\dagger}_{L;\bm{k}\alpha}H_{L;\alpha\beta}(\bm{k})c_{L;\bm{k}\beta},$
(2)

where $c^{\dagger}_{L;\bm{k}\alpha}$ ($c_{L;\bm{k}\alpha}$) creates
(annihilate) an electron with wavenumber $\bm{k}$ and orbital
$\alpha\in\\{s,p\\}$ in layer $L\in\\{1,2\\}$. Here,
$H_{L;\alpha\beta}(\bm{k})$ is a $2\times 2$ matrix asBernevig2006

$\displaystyle H_{L}(\bm{k})$ $\displaystyle=$
$\displaystyle\varepsilon_{L}(\bm{k})I_{2\times
2}+\sum_{i}d^{i}_{L}(\bm{k})\sigma_{i}$ (3)
$\displaystyle\varepsilon_{L}(\bm{k})$ $\displaystyle=$
$\displaystyle-2D_{L}\big{[}2-\cos(k_{x}-q_{L}^{x})-\cos(k_{y}-q_{L}^{y})\big{]}$
$\displaystyle d^{1}_{L}(\bm{k})$ $\displaystyle=$ $\displaystyle
A_{L}\sin(k_{x}-q_{L}^{x}),\quad d^{2}_{L}(\bm{k})=A_{L}\sin(k_{y}-q_{L}^{y})$
$\displaystyle d^{3}_{L}(\bm{k})$ $\displaystyle=$ $\displaystyle
M_{L}-2B_{L}\big{[}2-\cos(k_{x}-q_{L}^{x})-\cos(k_{y}-q_{L}^{y})\big{]},$

where $\sigma_{i}$ are the Pauli matrices acting on the orbital space, and
$\bm{q}_{L}=(q_{x},q_{y})$ represents the momentum shift of the band structure
in the Brillouin zone for the $L$-th layer which helps us to visually
distinguish two groups of edge states (originated from two layers) as Fig. (2)
shown. The real space version of the $L$-th layer Hamiltonian
$\mathcal{H}_{L}=\sum_{ij,\alpha\beta}c^{\dagger}_{L;i\alpha}H_{L;\alpha\beta}(i,j)c_{L;j\beta}$
can be obtained from Eqs. (2) and (3) by performing a straightforward inverse
Fourier transformation
$c_{L;\bm{k}\beta}=\frac{1}{\sqrt{V}}\sum_{i}c_{L;i\beta}e^{-i\bm{k}\cdot{\bm{r}_{i}}}$,
where $i$ is the site index.

In the absence of the inter-layer coupling $H_{c}$, the band structure and
Chern number of each layer can be tuned independently by varying model
parameters. For example, the band gap is $2|M_{L}|$, and the Chern number

$C_{L}=\left\\{\begin{array}[]{lll}+1,&\quad 0<M_{L}/2B_{L}<2\\\ 0,&\quad
M_{L}/2B_{L}<0\\\ -1,&\quad.\end{array}\right.$ (4)

The inter-layer coupling is considered to be the simple form in real space as

$\mathcal{H}_{c}=\sum_{i\alpha}\big{(}tc^{\dagger}_{1;i\alpha}c_{2;i\alpha}+\mathrm{H.c.}\big{)},$
(5)

where $t$ is the inter-layer coupling strength. If the adiabatic tuning of the
inter-layer coupling term (5) does not close the bulk gap, the Chern number of
the bilayer system is just the sum of those of each layer, $C=C_{1}+C_{2}$.
Hence, it is no surprise that the Chern number of the system changes with the
band gap $|M_{L}|$ and the inter-layer coupling strength $t$Wang2019 .

In Fig. 2, we present band structures of the quasi-one-dimensional (quasi-1D)
ribbon with several typical settings of model parameters. The number of edge
states (red curves) in the bulk band gap reflects the total Chern number at
half filling. For example, there is one group of edge states in Fig. 2 (b) so
that the Chern number of the system is $1$. Similarly, the Chern numbers in
panel (a), (c) and (d) are $C=0$, $C=2$ and $C=2$ respectively. From the
parameter $q_{L}^{x}$, it can be deduced that the edge state on the left
($k_{x}<0$) originates from the second layer, and the other pair originates
from the first layer. The model parameters $D_{L}$ determine the e-h symmetry.
Only when both of $D_{L}$($L\in\\{1,2\\}$) are zero, the band structure is
symmetric respect to energy $E=0$. This symmetry is broken in Fig. 2 (d) with
$D_{1}=0.1$.

Figure 2:  The band structure of the bilayer QAH model in the quasi-one-
dimensional ribbon geometry with width $Ny=80$. Red lines are edge states. The
four panels correspond to different parameter settings. (a) $D_{1}=0$,
$M_{1}=-0.3$, $M_{2}=-0.5$. (b) $D_{1}=0$, $M_{1}=0.3$, $M_{2}=-0.5$. (c)
$D_{1}=0$, $M_{1}=0.3$, $M_{2}=0.5$. (d) $D_{1}=0.1$, $M_{1}=0.3$,
$M_{2}=0.5$. The other model parameters are identical: $A_{1}=A_{2}=0.3$,
$B_{1}=B_{2}=0.2$, $D_{2}=0$, $t=0.1$. Particle-hole symmetry is preserved in
panel(a)-(c), and is broken in panel(d).

##  III III. Methods

In the modern theory, the macroscopic OM of a general two dimensional
crystalline solid can be expressed as an integration over the Brillouin
zoneThonhauser2005 ; Resta2005 ; Ceresoli2006 ,

$\displaystyle M$
$\displaystyle=M_{\mathrm{LC}}+M_{\mathrm{IC}}+M_{\mathrm{BC}},$ (6)
$\displaystyle M_{\mathrm{LC}}$ $\displaystyle=\frac{e}{2\hbar
c}\mathrm{Im}\sum_{n}\int_{\varepsilon_{nk}\leqslant\mu}\frac{\mathrm{d}\bm{k}}{(2\pi)^{2}}\langle\partial_{\bm{k}}u_{n\bm{k}}|\times
H_{\bm{k}}|\partial_{\bm{k}}u_{n\bm{k}}\rangle$ $\displaystyle
M_{\mathrm{IC}}$ $\displaystyle=\frac{e}{2\hbar
c}\mathrm{Im}\sum_{n}\int_{\varepsilon_{nk}\leqslant\mu}\frac{\mathrm{d}\bm{k}}{(2\pi)^{2}}\varepsilon_{n\bm{k}}\langle\partial_{\bm{k}}u_{n\bm{k}}|\times|\partial_{\bm{k}}u_{n\bm{k}}\rangle$
$\displaystyle M_{\mathrm{BC}}$ $\displaystyle=-\mu\frac{e}{2\pi\hbar
c}\mathcal{C}$ $\displaystyle\mathcal{C}$
$\displaystyle=\frac{1}{2\pi}\mathrm{Im}\sum_{n}\int_{\varepsilon_{nk}\leqslant\mu}\mathrm{d}\bm{k}\langle\partial_{\bm{k}}u_{n\bm{k}}|\times|\partial_{\bm{k}}u_{n\bm{k}}\rangle$

where $|u_{n\bm{k}}\rangle$ is the cell-periodic Bloch function, and
$\varepsilon_{nk}$ is the Bloch eigenvalue which satisfies
$H_{\bm{k}}|u_{n\bm{k}}\rangle=\varepsilon_{nk}|u_{n\bm{k}}\rangle$. All
summations in equation (6) are over occupied bands $n$ up to $\mu$. Three
constituent terms of $M$ correspond to the local circulation (LC), the
itinerant circulation (IC), and the Berry curvature (BC) contributions
respectively. The dimensionless number $\mathcal{C}$ is quantized as the Chern
number $C$ (the topological invariant) when $\mu$ is in the bulk gap, but may
not be quantized if $\mu$ is not in the gap. In numerical calculations,
partial derivatives in the above equations cannot be simply replaced by finite
differences. Instead, a fixed gauge must be guaranteed along neighboring mesh
points. Details of this algorithm is introduced in the Appendix.

Analogous to the above $k$ space integration, OM and its constituent terms can
also be expressed as a real space integration asBianco2013 ; MarrazzoResta2017

$\displaystyle M_{\alpha}$
$\displaystyle=\frac{1}{A}\int_{A}d{\bm{r}}\;m_{\alpha}\left(\bm{r}\right),\quad\alpha\in\\{\mathrm{LC,\;IC,\;BC}\\}$
(7) $\displaystyle m\left(\bm{r}\right)$
$\displaystyle={m_{\mathrm{LC}}\left(\bm{r}\right)}+{m_{\mathrm{IC}}\left(\bm{r}\right)}+{m_{\mathrm{BC}}\left(\bm{r}\right)}$
$\displaystyle{m_{\mathrm{LC}}\left(\bm{r}\right)}$
$\displaystyle=\frac{e}{\hbar
c}\mathrm{Im}\langle\bm{r}|PxQHQyP|\bm{r}\rangle$
$\displaystyle{m_{\mathrm{IC}}\left(\bm{r}\right)}$
$\displaystyle=-\frac{e}{\hbar
c}\mathrm{Im}\langle\bm{r}|QxPHPyQ|\bm{r}\rangle$
$\displaystyle{m_{\mathrm{BC}}\left(\bm{r}\right)}$
$\displaystyle=-\mu\frac{e}{2\pi\hbar\ c}\emph{c}\left(\bm{r}\right)$
$\displaystyle{\emph{c}\left(\bm{r}\right)}$
$\displaystyle=4\pi\mathrm{Im}\langle\bm{r}|QxPyQ|\bm{r}\rangle$

where the ground-state projector $P$ and its complement $Q$ satisfy $Q=1-P$
(see Appendix for details). Here $m_{\mathrm{LC}}\left(\bm{r}\right)$,
$m_{\mathrm{IC}}\left(\bm{r}\right)$ and $m_{\mathrm{BC}}\left(\bm{r}\right)$
correspond to the local markers of the corresponding terms in Eq.(6). The
local Chern number $\emph{c}\left(\bm{r}\right)$ is also called the
topological marker, whose real space distribution has attracted some
attentionsCaio2019 ; Ulcakar2020 .

Now the total OM $M$ and its three constituent terms $M_{\alpha}$ have the
form of a real space average over some region with an area $A$, i.e.,

$\frac{1}{A}\int_{A}d{\bm{r}}\langle\cdots\rangle.$

This region can be the whole sample, or just the bulk region around the sample
center (i.e., excluding the boundary region), as illustrated in Fig. 1 (b). As
for the Chern number, it is worth emphasizing that
$\emph{c}\left(\bm{r}\right)$ have opposite signs in the boundary and bulk
regions: $\emph{c}\left(\bm{r}\right)$ in the bulk region are equal to the
Chern number $C$, while they have opposite values in the boundary region, so
that its average over the whole sample always vanishBianco2011 ; Bianco2013 .
This is related to the mathematical fact of a projector over a finite-
dimensional manifoldBianco2011 .

Hence for calculating the OM by using Eq. (7), this feature has no effect on a
normal insulator about whether the trace (or integration) should be taken over
the whole sample or the bulk region, but it is different for Chern insulators
due to the existence of the BC term. Previous works indicate that the trace of
$\mathit{m}\left(\bm{r}\right)$ should also be taken over the bulk region for
Chern insulatorsBianco2013 . Nevertheless, a remaining question is whether
there is any underlying physical meaning for ignoring boundaries. To figure
this out, it is necessary to study the influence of the sample boundary on
$M$, especially how the boundary affect the other two terms of $M$, namely,
$M_{\mathrm{IC}}$ and $M_{\mathrm{LC}}$.

In the following, we will investigate the influence of boundaries on $M$ and
its three constituent terms $M_{\alpha}$, and check whether the average can be
taken over the whole sample or not. This problem is important for the
calculation of $M$ itself, especially for samples with inhomogeneity or with
small sizes. Furthermore, it will offer insights on the real space algorithm
of axion magnetoelectric couplingWang2019 .

##  IV IV. Normal Insulator

Figure 3:  Bilayer with $C=0$. Upper row: band structure of a quasi-1D ribbon.
Middle row: total OM as a function of sample size, calculated from three
methods. Lower row: $M_{\mathrm{LC}}$(black lines) and $M_{\mathrm{IC}}$(red
lines) as functions of the sample size $N$. The left (right) column
corresponds to the case with (without) e-h symmetry, by setting
$D_{1}=D_{2}=0$ ($D_{1}=0.1,D_{2}=0$). Square: OM from the whole sample
average. Circle: from bulk region. Up triangle: from $k$ space method. Other
model parameters are identical for all panels: $A_{1}=A_{2}=0.3$,
$B_{1}=B_{2}=0.2$, $M_{1}=-0.3$, $M_{2}=-0.5$, $t=0.1$. The Fermi energy
$\mu=0.05$.

Let us first focus on the topologically trivial insulator as shown in Fig. 3.
The left (right) column of Fig. 3 corresponds to the case of e-h symmetric
(asymmetric) band structure respect to $E=0$. To see the convergence of the
$r$ space method with increasing sample size $N$ (corresponding to a bilayer
with $2N^{2}$ lattice sites), we plot corresponding OMs calculated from an
average over the whole sample (black line with square dots) and the bulk
region (red line with circle dots) in Fig. 3 (b) and 3 (e). Here, we fix
$\mu=0.05$ in the band gap. As illustrated in Fig. 1, the bulk region refers
to the inner part of the sample enclosed by the blue dashed square, while the
rest part is the boundary region. OMs from the $k$ space method is also shown
as a comparison. In the following, we select the central quarter
($\frac{N}{2}\times\frac{N}{2}$) of the sample as the bulk region, unless
otherwise specified (e.g., in Section VI). The thickness of the boundaries is
$L_{B}=\frac{N}{4}$.

One can see that $M$ calculated from $r$ space methods (black and red line)
converge to the thermodynamic limit obtained by $k$ space method (blue line),
and the average of $\mathit{m}\left(\bm{r}\right)$ over bulk region converges
faster than over the whole sample, which is consistent with previous
studiesBianco2013 . In addition, the value of $|M|$ is almost zero when the
system satisfies e-h symmetry in Fig. 3(b), while the value of $|M|$ is
greatly increased by breaking e-h symmetry in Fig. 3 (e).

To figure out the reason for the increase of $|M|$ as the e-h symmetry is
broken, we further plot its constituent terms, namely $M_{\mathrm{LC}}$ and
$M_{\mathrm{IC}}$ as functions of sample size $N$ in Fig. 3 (c) and 3 (f).
Notice now the term $M_{\mathrm{BC}}$=0. In Fig. 3 (c), $M_{\mathrm{LC}}$ and
$M_{\mathrm{IC}}$ are perfectly opposite to each other,
$M_{\mathrm{LC}}=-M_{\mathrm{IC}}$, resulting in a vanishing $M$. This can be
understood from Eq. (A3) in Appendix. When the system satisfies e-h symmetry
at half filling, the ground state projector $P$ becomes completely equivalent
to $Q$, and then
$\mathrm{Tr}\left\\{PxQHQyP\right\\}=\mathrm{Tr}\left\\{QxPHPyQ\right\\}$. So,
the expressions of $M_{\mathrm{LC}}$ and $M_{\mathrm{IC}}$ only differ by a
minus sign. This perfect cancelation does not hold when the e-h symmetry is
broken, since $P$ and $Q$ are no longer equivalent, as can be seen in Fig.
3(f). Hence, breaking e-h is an effective method to enhance the OM at half
filling.

Figure 4:  The OM $M$ and its constituents $M_{\mathrm{LC}}$,
$M_{\mathrm{IC}}$ and $M_{\mathrm{BC}}$ as functions of the Fermi energy
$\mu$. Upper row: OMs from the whole sample average. Middle row: OMs from the
bulk region average. Lower row: comparison of total OM from two methods of
real space average. The left (right) column corresponds to case with (without)
e-h symmetry by setting $D_{1}=D_{2}=0$ ($D_{1}=0.1,D_{2}=0$). The other model
parameters are the same as Fig. 3. The sample size $N=40$.

If the Fermi energy (chemical potential) $\mu$ is moved into the the valence
or conductance band, the system will be metallic. Now the expression of $M$
has an extra term $M_{\mathrm{BC}}$ depending on $\mu$ in Eq. (6) or (7),
which was zero when $\mu$ is in the band gap. In Fig. 4, we present the Fermi
energy dependence of OMs, for the e-h symmetric (left column) and asymmetric
(right column) cases, respectively. We find that in both cases, $|M|$
increases rapidly once $\mu$ enters either band until an obvious peak, as
shown in Fig. 4 (c) and (f). To see the origins of the OM, in Fig. 4, we still
display three constituent terms $M_{\mathrm{LC}}$, $M_{\mathrm{IC}}$ and
$M_{\mathrm{IC}}$, averaged over the whole sample (first row) and over the
bulk region (second row). Most contributions are from $M_{\mathrm{LC}}$ and
$M_{\mathrm{IC}}$ (which will not be cancelled out due to the inequivalence
between projectors $P$ and $Q$), but their weights may be different in
different energy regions. For example, as can be seen from panels (a) and (d),
the dominating contribution to the conduction (valence) band peak is from
$M_{\mathrm{IC}}$ ($M_{\mathrm{LC}}$).

More details can be noticed by comparing Fig. 4 (a) between (b), and (d)
between (e). Concrete values of $M_{\mathrm{LC}}$, $M_{\mathrm{IC}}$ and
$M_{\mathrm{IC}}$ may be different, when the average is over the whole sample
(first row) or over the bulk region (second row). However, their sums, the
total OM from both algorithms, matches perfectly at all Fermi energies (last
row). This preliminarily indicates that, even with nonzero contribution from
$M_{\mathrm{BC}}$, the boundary can safely be included in the calculation of
$M$ when using the real space method. We further verify that this is the case
for Chern insulators.

##  V V. Chern Insulator

###  V.1 A. Chern number $C=1$

In this part, we focus on the topologically non-trivial phase of the bilayer
system. Let us start from the Chern number $C=1$ case with band structures
shown in Fig. 5 (a) (e-h symmetric) and 5 (e) (e-h asymmetric). There is only
one group of edge states (red lines) originating from the first layer of the
system.

Figure 5:  Bilayer with $C=1$. First row: band structures of a quasi-1D
ribbon. Second row: the total OM as a function of the sample size $N$
calculated from different methods: whole sample average (black), bulk region
average (red), and $k$ space method (blue). Third row: $M_{\mathrm{LC}}$
(black lines) ,$M_{\mathrm{IC}}$ (red lines) and $M_{\mathrm{BC}}$ (blue
lines) as functions of the sample size $N$. Fourth row: The comparison of
$M_{\mathrm{BC}}$ (black line) with $\Delta(M_{\mathrm{LC}}+M_{\mathrm{IC}})$
(red line) as functions of the sample size $N$. The left (right) column
corresponds to the case with (without) e-h symmetry. Square: whole sample
average. Circle: bulk region average. Up triangle: $k$ space method.
$M_{1}=0.3$, and other model parameters are identical to Fig. 3.

We set the Fermi energy $\mu=0.05$ well inside the bulk gap, and calculate $M$
and its local distribution by using Eqs. Eq. (6) and (7). First, let us still
check their convergence with increasing size $N$. The second row of Fig. 5
show comparisons of $M$ averaged over the whole sample, the bulk, and from the
$k$ space method. Again, all three results agree in the thermodynamic limit,
even though the system is topologically non-trivial. This is different from
previous conclusions which indicate that the real space average should be
taken over the bulk region only for Chern insulatorsBianco2013 .

Figure 6:  The total OM as a function of the Fermi energy $\mu$ for the
bilayer with $C=1$, calculated from whole sample average (red) and $k$ space
method (black). (a) With e-h symmetry. (b) Without e-h symmetry. Model
parameters are identical to Fig. 5.

This can be attributed to the real space distribution of local Chern number
$c\left(\bm{r}\right)$. Only when averaged over the bulk region of the sample,
does the dimensionless number $\mathcal{C}$ equal to the Chern number $C$.
However, when we take the average of $c\left(\bm{r}\right)$ over the whole
sample, $\mathcal{C}$ will always vanish, and therefore the contribution from
the corresponding term $M_{\mathrm{BC}}=-\mu\frac{e}{2\pi\hbar c}\mathcal{C}$
remains zero.

To further clarify this convergence, we turn to details of three constituent
parts, $M_{\mathrm{LC}}$ and $M_{\mathrm{IC}}$ and $M_{\mathrm{BC}}$,
presented in the third row of Fig. 5. In previous case of normal insulator, we
have shown that in the special case of e-h symmetry, $M_{\mathrm{LC}}$ and
$M_{\mathrm{IC}}$ at half filling are opposite to each other. We find that
this is also the case here with $C=1$, but only when the average is over the
bulk region. As displayed in Fig. 5(c), when averaged over bulk,
$M_{\mathrm{LC}}$ (black circle line) and $M_{\mathrm{IC}}$ (red circle line)
are symmetric with respect to zero. But when the average is over the whole
sample, they are no longer symmetric. Meanwhile, $M_{\mathrm{BC}}$ changes
from a finite value (blue circle line) to zero (blue square line) as the
average changes from bulk region to the whole sample. This indicates that the
existence of boundary states not only affects the geometric term
$M_{\mathrm{BC}}$, but also affects two non-geometric terms, $M_{\mathrm{LC}}$
and $M_{\mathrm{IC}}$.

Now let us scrutinize the effect of boundaries. For these two types of real
space averages of OM, i.e., over the whole sample and over the bulk region,
their difference is that of the boundary. For example, let us define the
boundary contribution of non-geometric terms,

$\Delta(M_{\mathrm{LC}}+M_{\mathrm{IC}})=(M_{\mathrm{LC}}+M_{\mathrm{IC}})_{t}-(M_{\mathrm{LC}}+M_{\mathrm{IC}})_{b}$
(8)

where the subscripts $t$ and $b$ denote “whole sample” and “bulk region”,
respectively. Their convergence processes are plotted in the last row of Fig.
5. For comparison, the term $M_{\mathrm{BC}}$ averaged over the whole sample
and that over the bulk region are also plotted. It is interesting to notice
that the magnitude of $\Delta(M_{\mathrm{LC}}+M_{\mathrm{IC}})$ converges to
the finite value of $M_{\mathrm{BC}}$ averaged over the bulk region in Fig.
5(d). Remember that
$M_{\mathrm{BC}}^{\mathrm{bulk}}=-M_{\mathrm{BC}}^{\mathrm{boundary}}$Bianco2013
. Therefore in other words, including the boundary allows $M_{\mathrm{LC}}$
and $M_{\mathrm{IC}}$ to complement the loss of $M_{\mathrm{BC}}$ averaged
over the boundary, which we call it the compensation effect. It similarly
makes sense for the system without e-h symmetry. In Fig. 5(g), owing to e-h
symmetry broken, $M_{\mathrm{LC}}$ and $M_{\mathrm{IC}}$ are no longer
opposite, even for a bulk average. Along with the influence of the boundary,
the asymmetry between them is further enlarged. Nevertheless, the magnitude of
$\Delta(M_{\mathrm{LC}}+M_{\mathrm{IC}})$ still converges to the nonzero value
of $M_{\mathrm{BC}}$ with increasing sample size $N$ in Fig. 5(h). Therefore,
it make sure that the total OM calculated by the whole sample average method
in $r$ space [Eq. (7)] is consistent with the $k$ space method Eq. (6).

We plot the total OM $M$ as a function of the Fermi energy $\mu$ in Fig. 6,
calculated from the $k$ space method (black line) and the whole sample average
method (red line). For all $\mu$, and independent of e-h symmetry, they match
perfectly. This effectively verifies the correctness of including the
boundary. As for the magnitude of OM, an important observation is that its
maximum values sharply appear at the top of the valence band and the bottom of
the conduction band, where the topological edge states end. By comparing with
the case of $C=0$ in Fig. 4 (c) and (f), we can see that now the maximum OM
has been increased remarkably. This suggests that nontrivial topology is still
an efficient means to enhance OM. Similar phenomenon has been observed in
concrete materialsXiDai2019 ; JPHanke2016 ; Lux2018 ; Helin2020 .

Figure 7:  The OM $M$ and its constituents: $M_{\mathrm{LC}}$,
$M_{\mathrm{IC}}$ ,$M_{\mathrm{BC}}$, as functions of the Fermi energy $\mu$
for the QAH bilayer with $C=1$. Upper row: average over the whole sample.
Middle row: average over the bulk. Lower row: $k$ space method. The left
(right) column corresponds to the case with (without) e-h symmetry. Model
parameters are identical to Fig. 5.

For details, three terms as functions of the Fermi energy $\mu$ are further
investigated, by using the $k$ space method, the bulk average and the whole
sample average in $r$ space. The results are shown in Fig. 7. We firstly focus
on the left row, associated with the e-h symmetric case. From the comparison
of Fig. 7(b) and (c), we can see that for each term, the result from the bulk
average is consistent with that from the $k$ space method. For the whole
sample average, although the total OM still matches (Fig. 6), such a term-to-
term correspondence with $k$ space method does not hold any more [Fig. 7(a)
and (c)]. We will check what happens to each constituent term.

First, let us focus on the the bulk gap region, $\mu\in(-0.3,0.3)$. For the
bulk average, $M_{\mathrm{LC}}$ and $M_{\mathrm{IC}}$ remain unchanged with
varying $\mu$, while $M_{\mathrm{BC}}$ changes linearly. The slope satisfies
$\frac{\partial M}{\partial\mu}=\frac{\partial
M_{\mathrm{BC}}}{\partial\mu}\propto C$ which can be deduced from Eqs. (6) or
(7). When we switch to the whole sample average as shown in Fig. 7(a), on the
other hand, $M_{\mathrm{BC}}$ is always $0$ when $\mu$ is within the bulk gap.
This is consistent with the previous workBianco2013 . Surprisingly, within
this energy region $M_{\mathrm{LC}}$ and $M_{\mathrm{IC}}$ change linearly
with $\mu$, instead of remaining constant in Fig. 7 (b). The joint effect of
them results in a linear relationship between $M$ and $\mu$, and
$\frac{\partial M}{\partial\mu}\propto C$ still holds. Such a story is
similarly applicable to the system with broken e-h symmetry as presented in
Fig. 7(d)-7(f). This further verifies the above conclusion that the
compensation effect of the boundary ensures $M_{\mathrm{LC}}$ and
$M_{\mathrm{IC}}$ to supplement the loss of $M_{\mathrm{BC}}$, when the real
space average is switched from bulk region to whole sample. Thus, it is also
reliable to average over the whole sample for $r$ space method even for the
case of $C=1$.

###  V.2 B. Chern number $C=2$

Figure 8:  Similar to Fig. 5 but for the bilayer with $C=2$. First row: band
structures of a ribbon. Second row: the total OM as a function of the sample
size $N$ calculated from three methods. Third row: $M_{\mathrm{LC}}$
,$M_{\mathrm{IC}}$ and $M_{\mathrm{BC}}$ as functions of the sample size $N$.
Fourth row: comparison of $M_{\mathrm{BC}}$ with
$\Delta(M_{\mathrm{LC}}+M_{\mathrm{IC}})$ as functions of the sample size $N$.
The left (right) column corresponds to the case with (without) e-h symmetry.
$M_{2}=0.5$, and other model parameters are the same as Fig. 5.

The above conclusions can be generalized to the case with $C=2$. Similar to
Fig. 5, in Fig. 8 we present the size convergence of total OM from different
methods (second row), three constituent terms (third row) and the contribution
from the boundary (fourth row). Similar to the case of $C=1$, $M$ calculated
from the whole sample method fluctuates more than that from other methods.
However, it still converges to the correct value for sufficiently large
sample. As a matter of fact, as presented in Fig. 9, the total OM from the
whole sample average (red line) is perfectly consistent with that from the $k$
space method (black line) for all Fermi energies $\mu$.

Figure 9:  Similar to Fig. 6 but for the bilayer with $C=2$, total OM as a
function of the Fermi energy $\mu$, calculated from two methods. (a) With e-h
symmetry. (b) Without e-h symmetry. Model parameters are identical to Fig. 8.

Also similar to the $C=1$ case, the inclusion of boundary affects all three
components: $M_{\mathrm{LC}}$, $M_{\mathrm{IC}}$ and $M_{\mathrm{BC}}$. If the
system has e-h symmetry, $M_{\mathrm{LC}}$ (black circle line) and
$M_{\mathrm{IC}}$ (red circle line) are symmetric with respect to zero in Fig.
8(c) when averaged over bulk, while they are no longer symmetric when averaged
over the whole sample. Independent of e-h symmetry, the compensation effect we
met in the case of $C=1$ does work here. For example,
$\Delta(M_{\mathrm{LC}}+M_{\mathrm{IC}})$ still converges to the finite value
of $M_{\mathrm{BC}}$ averaged over the bulk region as shown in the last row of
Fig. 8. Furthermore, such a compensation effect works for all $\mu$ values,
including the bulk gap with large orbital magnetic moment from edge states, as
presented in Fig. 10.

Figure 10:  Similar to Fig. 7 but for the bilayer with $C=2$. The OM and its
constituents as functions of the Fermi energy $\mu$. Upper row: average over
the whole sample. Middle row: average over the bulk. Lower row: $k$ space
method. The left (right) column corresponds to the case with (without) e-h
symmetry. Model parameters are identical to Fig. 8.

One more thing should also be noticed by comparing the vertical scales in Fig.
9 and in Fig. 6. The peak value of the total OM with $C=2$ is even larger than
that in the case with $C=1$. This again confirms that nontrivial topology,
i.e., Chern number, or at least remarkably nonzero integrated Berry curvature,
is one of the most effective factor to enhance OMJPHanke2016 ; Lux2018 .

##  VI VI. Size Effects

In this section, we discuss some size effects of the real space algorithm. In
previous studies, when discussing the roles of boundaries, we separated the
bulk and boundary regions in a manner of a fixed fraction: the central quarter
of the sample is adopted as the bulk, and the rest part is the boundary. Now
we discuss the effect of different manners of separation by varying the
thickness of the boundary, $L_{B}$. In Fig. 1, bulk averaged OMs with
different definitions of $L_{B}$ are plotted, for the cases of $C=1$ (left
column) and $C=2$ (right column) respectively. Notice the $L_{B}=0$ case
corresponds to the whole sample average. Here, the values of constituent terms
may fluctuate with different $L_{B}$ but the total OM is invariant and
consistent with that from the $k$ space method. In other words, for a
sufficiently large sample, the details of separating the bulk and boundary
regions are not important, and again, even the whole sample average is also
valid.

Figure 11:  The bulk averaged OM and its constituents as functions of the
Fermi energy $\mu$ with different adoptions of boundary thickness $L_{B}$. The
left (right) column denotes the bilayer with $C=1$ ($C=2$). First row: the
bulk average of $M$. Second row: the bulk average of $M_{\mathrm{LC}}$. Third
row: the bulk average of $M_{\mathrm{IC}}$. Fourth row: the bulk average of
$M_{\mathrm{BC}}$. Different colors represent different thickness of
boundaries. Orange lines signify the results of $k$ space method. The sample
size $N=40$, and other model parameters are identical to Fig. 5.

The focus so far is the thermodynamic limit, i.e., real space results for a
sufficiently large sample. With the progress of nanostructure fabrication
technologies in thin films, clusters and manmade structures, it is also
practical to examine the OM of very small samples. To check the validity of
current real space algorithm in this case, we calculate the finite-size OM
from the “first principle” expressionBianco2013 ; Thonhauser2005

$M=-\frac{e}{2cA}\sum_{i}\langle\varphi_{i}|\bm{r}\times\bm{v}|\varphi_{i}\rangle,$
(9)

and then compare it with those from $r$ space (whole-sample average) and $k$
space method. Since the size is very small (less than 100 sites), the
separation of bulk and boundaries makes no sense, and therefore the $r$ space
algorithm is just a whole-sample average. The data are shown in Fig. 12. It is
interesting to see the perfect consistency between the results from Eq. (9)
(black line) and the whole-sample average (red line). This example again
confirms the importance of this whole-sample average in real space
calculations.

Figure 12:  The total OM as a function of the Fermi energy $\mu$ for very
small samples. The left (right) column represent the sample size $N=5$
($N=10$). Upper row: the bilayer with $C=0$. Middle row: the bilayer with
$C=1$. Lower row: the bilayer with $C=2$. Black lines: the OM calculated from
the definition Eq. (9). Red lines: the OM calculated from $r$ space method Eq.
(7). Blue lines: the OM calculated from $k$ space method Eq. (6). In all
panels, black curves have been covered by red curves due to perfect
coincidence between them. Model parameters are identical to Fig. 5.

##  VII VII. Summary

In summary, based on a bilayer model with an adjustable Chern number, we
numerically investigated the role of boundaries in the real space algorithm of
calculating orbital magnetization $M$, by scrutinizing its constituent terms
$M_{\mathrm{LC}}$, $M_{\mathrm{IC}}$, $M_{\mathrm{BC}}$, as well as two
different averages: over the bulk or over the whole sample. The first key
finding is that for the total OM, both average methods converge to the same
value calculated from $k$ space method, for sufficiently large sample size.
This conclusion does not depend on the topological nature of the system, i.e.,
Chern number at half filling.

The boundary region has a significant contribution to the geometric term
$M_{\mathrm{BC}}$, especially in the case of nontrivial topology. However, it
also has a non-negligible contribution to non-geometric terms
$M_{\mathrm{LC}}$ and $M_{\mathrm{IC}}$, and these two types of contributions
compensate each other exactly, which results in the consistency of two average
methods in the real space. The compensation effect is further manifested by
the property that $M_{\mathrm{LC}}$ and $M_{\mathrm{IC}}$ inherit the linear
relationship between $M_{\mathrm{BC}}$ and the Fermi energy $\mu$ for Chern
insulators when $\mu$ is in the bulk band gap. This is different from the
previous conclusion which advocates that the boundaries of Chern insulators
should not be considered in $r$ space expression because of the questionable
contribution of $M_{\mathrm{BC}}$ from boundaries. From the comparison among
systems with different Chern numbers, the systems with high Chern number are
suitable candidate materials with strong orbital magnetism. This provides a
theoretical reference for finding ways to enhance orbital magnetism. This $r$
space method of a whole-sample average is also important for nanostructures,
where the separation between the bulk and boundaries seems technologically
problematic. These clarifications will also be instructive to further studies
on the orbital magnetoelectric effects in real space.

##  VIII Acknowledgements

This work was supported by National Natural Science Foundation of China under
Grant Nos. 12104108, 11774336, and 12147112, and the Starting Research Fund
from Guangzhou University under Grant Nos. RQ2020082 and 62104360.

##  IX Appendix: Calculation of Orbital Magnetization

For a two dimensional material, the OM $M$ is defined as the orbital moment
per unit areaCeresoli2006 ,

$M=-\frac{e}{2cA}\sum_{i}\langle\varphi_{i}|\bm{r}\times\bm{v}|\varphi_{i}\rangle$
(A1)

where $e$ is the magnitude of the electronic charge, $c$ is the vacuum speed
of light and $A$ is the sample area. $|\varphi_{i}\rangle$ are the occupied
eigenstates, $\bm{r}$ is the position operator and $\bm{v}$ is the velocity
operator which satisfies

$v=-\frac{i}{\hbar}\left[\bm{r},H\right]$ (A2)

Here, $H$ is the Hamiltonian operator. It is obvious that the above expression
can be directly applied to open boundary conditions(OBCs). Due to the non-
periodic nature of the position operator $r$, it requires further derivation
to obtain $k$ space expression for periodic boundary conditions(PBCs) which is
known as the modern theory of orbital magnetization.

For a two dimensional gapped insulator, electrons are in the $x$-$y$ plane and
the OM $M$ is along $z$ direction. Eq. (A1) can be rewritten in the $r$ space
as follows,

$\displaystyle M=M_{\mathrm{LC}}+M_{\mathrm{IC}}$ (A3) $\displaystyle
M_{\mathrm{LC}}=\frac{e}{\hbar cA}\mathrm{ImTr}\left\\{PxQHQyP\right\\}$
$\displaystyle M_{\mathrm{IC}}=-\frac{e}{\hbar
cA}\mathrm{ImTr}\left\\{QxPHPyQ\right\\}$

where $M_{\mathrm{LC}}$ and $M_{\mathrm{IC}}$ correspond to the local
circulation (LC), the itinerant circulation (IC) respectively. The ground-
state projector $P$ and its complement $Q$ satisfy

$\displaystyle
P=\sum_{\varepsilon_{i}\leqslant\mu}|\varphi_{i}\rangle\langle\varphi_{i}|$
$\displaystyle Q=1-P$ (A4)

where $\mu$ is the Fermi energy, and $\varepsilon_{i}$ is the eigenvalue
evaluated by
$H({\bm{r}})|\varphi_{i}\rangle=\varepsilon_{i}|\varphi_{i}\rangle$. The
position operator $r$ becomes well-defined when ”sandwiched” between a $P$ and
a $Q$ within PBCs. Generalizing equation (A3) to Chern insulators and metals,
there is one more component coming from the Berry curvature(BC) contributions.
We call it $M_{\mathrm{BC}}$, which is associated with the topological
properties of materialsWang2020 . Hence,

$\displaystyle M=M_{\mathrm{LC}}+M_{\mathrm{IC}}+M_{\mathrm{BC}}$ (A5)
$\displaystyle M_{\mathrm{LC}}=\frac{e}{\hbar
cA}\mathrm{ImTr}\left\\{PxQHQyP\right\\}$ $\displaystyle
M_{\mathrm{IC}}=-\frac{e}{\hbar cA}\mathrm{ImTr}\left\\{QxPHPyQ\right\\}$
$\displaystyle M_{\mathrm{BC}}=-\mu\frac{e}{2\pi\hbar cA}\mathcal{C}$
$\displaystyle\mathcal{C}=4\pi\mathrm{ImTr}\left\\{QxPyQ\right\\}$

In addition, the $r$ space expression also gives the local description of $M$

$\displaystyle M=\frac{1}{A}\int d{\bm{r}}\mathit{m}\left(\bm{r}\right)$ (A6)
$\displaystyle{\mathit{m}\left(\bm{r}\right)}=m_{\mathrm{LC}}\left(\bm{r}\right)+m_{\mathrm{IC}}\left(\bm{r}\right)+m_{\mathrm{BC}}\left(\bm{r}\right)$
$\displaystyle{m_{\mathrm{LC}}\left(\bm{r}\right)}=\frac{e}{\hbar
c}\mathrm{Im}\langle\bm{r}|PxQHQyP|\bm{r}\rangle$
$\displaystyle{m_{\mathrm{IC}}\left(\bm{r}\right)}=-\frac{e}{\hbar
c}\mathrm{Im}\langle\bm{r}|QxPHPyQ|\bm{r}\rangle$
$\displaystyle{m_{\mathrm{BC}}\left(\bm{r}\right)}=-\mu\frac{e}{2\pi\hbar\
c}\emph{c}\left(\bm{r}\right)$
$\displaystyle{\emph{c}\left(\bm{r}\right)}=4\pi\mathrm{Im}\langle\bm{r}|QxPyQ|\bm{r}\rangle$

where $m_{\mathrm{LC}}\left(\bm{r}\right)$,
$m_{\mathrm{IC}}\left(\bm{r}\right)$ and $m_{\mathrm{BC}}\left(\bm{r}\right)$
correspond to the local properties of the three terms in Eq.(A5).
$\emph{c}\left(\bm{r}\right)$ is called the topological marker.

Using the relation between Wannier representation and Bloch representation,
Eq.(A5) can be linked with $k$ space expression which is called the modern
theory of OM for a Chern insulatorCeresoli2006 ; Malashevich2010 . In the $k$
space, the $M$ is expressed as a Brillouin-zone integral of Bloch-orbital
matrix elementsResta2010 ; Malashevich2010 ; Ceresoli2006 . The macroscopic OM
of a general two dimensional crystalline solid is

$\displaystyle M=M_{\mathrm{LC}}+M_{\mathrm{IC}}+M_{\mathrm{BC}}$ (A7)
$\displaystyle M_{\mathrm{LC}}=\frac{e}{2\hbar
c}\mathrm{Im}\sum_{n}\int_{\varepsilon_{nk}\leqslant\mu}\frac{\mathrm{d}\bm{k}}{(2\pi)^{2}}\langle\partial_{\bm{k}}u_{n\bm{k}}|\times
H_{\bm{k}}|\partial_{\bm{k}}u_{n\bm{k}}\rangle$ $\displaystyle
M_{\mathrm{IC}}=\frac{e}{2\hbar
c}\mathrm{Im}\sum_{n}\int_{\varepsilon_{nk}\leqslant\mu}\frac{\mathrm{d}\bm{k}}{(2\pi)^{2}}\varepsilon_{n\bm{k}}\langle\partial_{\bm{k}}u_{n\bm{k}}|\times|\partial_{\bm{k}}u_{n\bm{k}}\rangle$
$\displaystyle M_{\mathrm{BC}}=-\mu\frac{e}{2\pi\hbar c}\mathcal{C}$
$\displaystyle\mathcal{C}=\frac{1}{2\pi}\mathrm{Im}\sum_{n}\int_{\varepsilon_{nk}\leqslant\mu}\mathrm{d}\bm{k}\langle\partial_{\bm{k}}u_{n\bm{k}}|\times|\partial_{\bm{k}}u_{n\bm{k}}\rangle$

where $|u_{n\bm{k}}\rangle$ is the cell-periodic Bloch function, and
$\varepsilon_{nk}$ is the Bloch eigenvalue which satisfies
$H_{\bm{k}}|u_{n\bm{k}}\rangle=\varepsilon_{nk}|u_{n\bm{k}}\rangle$. All
summations in equation (A7) are over occupied bands $n$ up to $\mu$. In
numerical calculations, the derivatives
$|{\partial}_{\bm{k\alpha}}u_{n\bm{k}}\rangle$ are evaluated on a mesh of
discretized Brillouin zone. Since a simple finite-deference of the wave
functions $|u_{n\bm{k}}\rangle$ cannot ensure a fixed gauge on neighboring
grid points, we adopt the discretized covariant derivative
$|\tilde{\partial}_{\bm{k\alpha}}u_{n\bm{k}}\rangle$ which involves linear
combination of occupied states and a local fixed gauge around a definite grid
point. The definition of $|\tilde{\partial}_{\bm{k\alpha}}u_{n\bm{k}}\rangle$
is

${|\tilde{\partial}_{\bm{k\alpha}}u_{n\bm{k}}\rangle}=Q_{\bm{k}}|\partial_{\bm{k\alpha}}u_{n\bm{k}}\rangle$
(A8)

where $Q_{\bm{k}}$ is the Fourier transformation of Eq.(A4). For more details
see appendix A of reference Ceresoli2006 . Let us defineCeresoli2006 ;
Souza2008

$\displaystyle\tilde{f}_{\bm{k},i}=\frac{1}{v}\in_{ijl}q_{i}\sum_{n}\mathrm{Im}\langle\tilde{\partial}_{\bm{kj}}u_{n\bm{k}}\
|\tilde{\partial}_{\bm{kl}}u_{n\bm{k}}\rangle$ (A9)
$\displaystyle\tilde{g}_{\bm{k},i}=\frac{1}{v}\in_{ijl}q_{i}\sum_{n}\mathrm{Im}\langle\tilde{\partial}_{\bm{kj}}u_{n\bm{k}}\
\mathrm{H}_{\bm{k}}|\tilde{\partial}_{\bm{kl}}u_{n\bm{k}}\rangle$
$\displaystyle\tilde{h}_{\bm{k},i}=\frac{1}{v}\in_{ijl}q_{i}\sum_{n}\varepsilon_{nm\bm{k}}\mathrm{Im}\langle\tilde{\partial}_{\bm{kj}}u_{m\bm{k}}\
|\tilde{\partial}_{\bm{kl}}u_{n\bm{k}}\rangle$

where $q_{i}$ denotes the primitive reciprocal vectors that define the $k$
mesh along the $i$th direction, and $v$ is the volume of the unit cell of the
$k$ mesh. Thus, we can further rewrite Eq.(A7) as follows

$\displaystyle
M=\tilde{M}_{\mathrm{LC}}+\tilde{M}_{\mathrm{IC}}+M_{\mathrm{BC}}$ (A10)
$\displaystyle\tilde{M}_{\mathrm{LC}}=\frac{-e}{2\hbar
c}\int_{\mathrm{BZ}}\frac{\mathrm{d}\bm{k}}{(2\pi)^{2}}\tilde{g}_{\bm{k}}$
$\displaystyle\tilde{M}_{\mathrm{IC}}=\frac{-e}{2\hbar
c}\int_{\mathrm{BZ}}\frac{\mathrm{d}\bm{k}}{(2\pi)^{2}}\tilde{h}_{\bm{k}}$
$\displaystyle M_{\mathrm{BC}}=-\mu\frac{e}{2\pi\hbar c}\mathcal{C}$
$\displaystyle\mathcal{C}=\frac{1}{2\pi}\int_{\mathrm{BZ}}\tilde{f}_{\bm{k}}$

In this manuscript, we have adopted both the $k$ space and $r$ space
expression to study OMs. We will focus on the influence of the boundaries on
the $M$ and its three components: $M_{\mathrm{LC}}$, $M_{\mathrm{IC}}$ and
$M_{\mathrm{BC}}$.

## References

  * (1) F. Mauri and S.G. Louie, Phys. Rev. Lett. 76, 4246 (1996).
  * (2) F. Mauri, B. G. Pfrommer, and S. G. Louie, Phys. Rev. Lett. 77, 5300 (1996).
  * (3) C. J. Pickard and F. Mauri, Phys. Rev. Lett. 88, 086403 (2002).
  * (4) D. Sebastiani, G. Goward, I. Schnell, and M. Parrinello, Comput. Phys. Commun. 147, 707 (2002).
  * (5) D. Vanderbilt, Berry Phase in Electronic Structure Theory (Cambridge University Press, Cambridge, UK, 2018)
  * (6) D. Ceresoli, T. Thonhauser, D. Vanderbilt, and R. Resta, Phys. Rev. B 74, 024408 (2006).
  * (7) R. Bianco and R. Resta, Phys. Rev. B 84, 241106(R) (2011).
  * (8) R. Bianco and R. Resta, Phys. Rev. Lett. 110, 087202 (2013).
  * (9) A. Malashevich, I. Souza, S. Coh, N. J. Phys. 12, 053032 (2010).
  * (10) J. Liu, Z. Ma, J. Gao, and X. Dai, Phys. Rev. X 9, 031021 (2019).
  * (11) J.-P. Hanke, F. Freimuth, A. K. Nandy, H. Zhang, S. Blügel, and Y. Mokrousov, Phys. Rev. B 94, 121114(R) (2016).
  * (12) F. R. Lux, F. Freimuth, S. Blügel, and Y. Mokrousov, Commun. Phys. 1, 60 (2018).
  * (13) S.-Y. Li, Y. Zhang, Y.-N. Ren, J. Liu, X. Dai, and L. He, Phys. Rev. B 102, 121406(R) (2020).
  * (14) S. S. Pershoguba and V. M. Yakovenko, Phys. Rev. B 105, 064423 (2022).
  * (15) D. Go, J.-P. Hanke, P. M. Buhl, F. Freimuth, G. Bihlmayer, H.-W. Lee, Y. Mokrousov, and S. Blugel, Sci. Rep. 7, 46742 (2017).
  * (16) S. Bhowal and S. Satpathy, Phys. Rev. B 101, 121112(R) (2020).
  * (17) D. Go , D. Jo, H.-W. Lee, M. Kläui, and Y. Mokrousov, EPL 135, 37001 (2021).
  * (18) L. Dong, C. Xiao, B. Xiong, and Q. Niu, Phys. Rev. Lett. 124, 066601 (2020).
  * (19) C. Xiao, Y.-F. Ren, and B.-G. Xiong, Phys. Rev. B 103, 115432 (2021).
  * (20) M. Hamada and S. Murakami, Phys. Rev. Research 2,023275 (2020).
  * (21) A.L. Sharpe, E.J. Fox, A.W. Barnard, J. Finney, K. Watanabe, T. Taniguchi, M. Kastner, and D. Goldhaber-Gordon, Science 365, 605 (2019).
  * (22) M. Serlin, C. Tschirhart, H. Polshyn, Y. Zhang, J. Zhu, K. Watanabe, T. Taniguchi, L. Balents, and A. Young, Science 367, 900 (2020).
  * (23) Y. Su and S.-Z. Lin, Phys. Rev. Lett. 125, 226401 (2020).
  * (24) W.-Y. He, D. Goldhaber-Gordon, and K. Law, Nat. Commun. 11, 1 (2020).
  * (25) C.-L. Huang, N.-M. Wei, and A.H. MacDonald, Phys. Rev. Lett. 126, 056801 (2021).
  * (26) M. Cinal, Phys. Rev. B 105, 104403 (2022).
  * (27) D. Xiao, J. Shi, and Q. Niu, Phys. Rev. Lett. 95, 137204 (2005).
  * (28) D. Xiao, M.-C. Chang, and Q. Niu, Rev. Mod. Phys. 82, 1959 (2010).
  * (29) T. Thonhauser, D. Ceresoli, D. Vanderbilt, and R. Resta, Phys. Rev. Lett. 95, 137205 (2005).
  * (30) R. Resta, D. Ceresoli, T. Thonhauser, and D. Vanderbilt, Chem. Phys. Chem 6, 1815 (2005).
  * (31) T. Thonhauser, J. Mod. Phys. B 25, 1429 (2011).
  * (32) J. Shi, G. Vignale, D. Xiao, and Q. Niu, Phys. Rev. Lett. 99, 197202 (2007).
  * (33) R. Bianco and R. Resta, Phys. Rev. B 93, 174417 (2016).
  * (34) E. Drigo and R. Resta, Phys. Rev. B 101, 165120 (2020).
  * (35) J. Sykes and R. Barnett, Phys. Rev. B 103, 155134 (2021).
  * (36) A. Marrazzo and R. Resta, Phys. Rev. B 95, 121114(R) (2017).
  * (37) T. Olsen, M. Taherinejad, D. Vanderbilt, and I. Souza, Phys. Rev. B 95, 075137 (2017).
  * (38) N. Varnava and D. Vanderbilt, Phys. Rev. B 98, 245117 (2018).
  * (39) N. Varnava, I. Souza, and David Vanderbilt, Phys. Rev. B 101, 155130 (2020).
  * (40) D. M. Nenno, C. A. C. Garcia, J. Gooth, C. Felser, and P. Narang, Nat. Rev. Phys. 2, 682 (2020).
  * (41) A. Sekine and K. Nomura, J. Appl. Phys. 129, 141101 (2021).
  * (42) M. D. Caio, G. Möller, N. R. Cooper, and M. J. Bhaseen, Nat. Phys. 15, 257 (2019).
  * (43) L. Ulčakar, J. Mravlje, and T. Rejec, Phys. Rev. Lett. 125, 216601 (2020).
  * (44) C. L. Tschirhart, M. Serlin, H. Polshyn, A. Shragai, Z. Xia, J. Zhu, Y. Zhang, K. Watanabe, T. Taniguchi, M. E. Huber, and A. F. Young, Science 372, 1323 (2021).
  * (45) V. Vargiamidis, P. Vasilopoulos, M. Tahir, and N. Neophytou, Phys. Rev. B 102, 235426 (2020).
  * (46) D. Hara, M.S. Bahramy, and S. Murakami, Phys. Rev. B 102, 184404 (2020).
  * (47) O. Antebi, A. Stern, and E. Berg, Phys. Rev. B 105, 104423 (2022).
  * (48) S-S. Wang, Y-Y. Zhang, J-H. Guan, Y. Yu, Y. Xia, and S-S. Li, Phys. Rev. B 99, 125414 (2019).
  * (49) J. Ge, Y. Liu, J. Li, H. Li, T. Luo, Y. Wu, Y. Xu, and J. Wang, Natl. Sci. Rev. 7, 1280-1287 (2020).
  * (50) Y.-F. Zhao, R. Zhang, R. Mei, L.-J. Zhou, H. Yi, Y.-Q. Zhang, J. Yu, R. Xiao, K. Wang, N. Samarth, M. H. W. Chan, C.-X. Liu, and C.-Z. Chang, Nature 588, 419 (2020).
  * (51) M. Algarni, C. Tan, G. Zheng, S. Albarakati, X. Zhu, J. Partridge, Y. Zhu, L. Farrar, M. Tian, J. Zhou, X. Wang, Z. Mao, and L. Wang, Phys. Rev. B 105, 155407 (2022).
  * (52) N. Goldman, J. C. Budich, and P. Zoller, Nat. Phys. 12, 639 (2016).
  * (53) T. Hartke, B. Oreg , N. Jia, and M. Zwierlein, Phys. Rev. Lett. 125, 113601 (2020).
  * (54) M. Gall, N. Wurz, J. Samland, C. F. Chan, M. Köhl, Nature 589, 40 (2021).
  * (55) A. B. Khanikaev and G. Shvets, Nat. Photonics 11, 763 (2017).
  * (56) A. Bernevig, T. Hughes and S. C. Zhang, Science 314, 1757 (2006).
  * (57) S-S. Wang, Y-Y. Zhang, J-H. Guan, Y. Yu, Y. Xia, and S-S. Li, J. Phys. Condens. Matter 32, 335302 (2020).
  * (58) R. Resta, J. Phys. Condens. Matter 22, 123201 (2010).
  * (59) I. Souza, and D. Vanderbilt, Phys. Rev. B 77, 054438 (2008).
  * (60) G. Zhu, S-A. Yang, C. Fang, W-M. Liu, and Y. Yao, Phys. Rev. B 86, 214415 (2012).