# Topological invariants for interacting topological insulators: I. Efficient
numerical evaluation scheme and implementations

Yuan-Yao He  Han-Qing Wu  Department of Physics, Renmin University of China,
Beijing 100872, China  Zi Yang Meng  Beijing National Laboratory for Condensed
Matter Physics, and Institute of Physics, Chinese Academy of Sciences, Beijing
100190, China  Zhong-Yi Lu  Department of Physics, Renmin University of China,
Beijing 100872, China

###### Abstract

The aim of this series of two papers is to discuss topological invariants for
interacting topological insulators (TIs). In the first paper (I), we provide a
paradigm of efficient numerical evaluation scheme for topological invariants,
in which we demystify the procedures and techniques employed in calculating
$Z_{2}$ invariant and spin Chern number via zero-frequency single-particle
Green’s function in quantum Monte Carlo (QMC) simulations. Here we introduce a
periodization process to overcome the ubiquitous finite-size effect, so that
the calculated spin Chern number shows ideally quantized values. We also show
that making use of symmetry properties of the underlying systems can greatly
reduce the computational effort. To demonstrate the effectiveness of our
numerical evaluation scheme, especially the periodization process, of
topological invariants, we apply it on two independent two-dimensional models
of interacting topological insulators. In the subsequent paper (II), we apply
the scheme developed here to wider classes of models of interacting
topological insulators, for which certain limitation of constructing
topological invariant via single-particle Green’s functions will be presented.

###### pacs:

71.10.-w, 71.10.Fd, 71.27.+a

##  I Introduction

Topological insulators (TIs) Hasan and Kane (2010); Qi and Zhang (2011) in
free fermion systems can be characterized by topological invariants, such as
$Z_{2}$ invariant and spin Chern number, and the calculation of these
topological invariants is straightforward via noninteracting Hamiltonian
matrix. However, the characterization of interacting topological insulators
via these topological invariants is still not well understood, both
conceptually and technically in numerical calculation. Of course, there are
already various achievements. In the quantum (anomalous) Hall insulators
Klitzing  _et al._ (1980); Haldane (1988) with broken time-reversal symmetry,
the characterizing topological invariant is the Thouless-Kohmoto-Nightingale-
den Nijs (TKNN) invariant Thouless  _et al._ (1982); Avron  _et al._ (1983),
or first Chern number, which takes integer values. For TIs protected by time-
reversal symmetry and charge $U(1)$ symmetry, the $Z_{2}$ topological index is
applied to describe the system Kane and Mele (2005a, b); Fu and Kane (2006).
With additional $U(1)$ spin rotational symmetry in TIs, the
$U(1)_{\text{spin}}\times U(1)_{\text{charge}}\rtimes Z_{2}^{T}$ symmetry
($Z_{2}^{T}$ stands for time-reversal symmetry) results in $\mathbb{Z}$
classification and the appropriate topological invariant is the spin Chern
number $C_{s}$ Sheng  _et al._ (2006); Prodan (2009), which is actually the
case for Kane-Mele model without Rashba spin-orbit coupling Kane and Mele
(2005a, b). In such systems, the spin Hall conductivity $\sigma_{xy}^{spin}$
is related to spin Chern number as $\sigma_{xy}^{spin}=C_{s}\frac{e}{2\pi}$.

For noninteracting TIs, both $Z_{2}$ invariant Kane and Mele (2005b); Fu and
Kane (2006) and spin Chern Number Sheng  _et al._ (2006); Prodan (2009) can be
simply evaluated from Hamiltonian matrix in band basis. For example, if a
system has spatial inversion symmetry, the $Z_{2}$ invariant can be calculated
as product of parity eigenvalues of all occupied energy bands at all time-
reversal invariant momentum (TRIM) points in Brillouin zone (BZ) Fu and Kane
(2007). Spin Chern number, on the other hand, can be calculated by simply
integrating Berry curvature over the Brillouin zone (BZ).

For interacting TIs, the evaluations of topological invariants become much
more difficult and subtle. The proposals include constructing topological
invariants from single-particle Green’s function Avron  _et al._ (1983); So
(1985); Ishikawa and Matsuyama (1986); Volovik (1988); Wang  _et al._ (2010);
Gurarie (2011); Volovik (2009) or imposing twisted boundary phases to the
ground state wavefunction Niu  _et al._ (1985); Wang and Zhang (2014).
Recently, the constructions of topological invariants from single-particle
Green’s function, especially the zero-frequency version Wang and Zhang
(2012a); Wang and Yan (2013), have been actively investigated. Implementations
of topological invariants constructed from zero-frequency single-particle
Green’s function Wang and Zhang (2012a) in many-body numerical techniques have
been carried out in various studies. In one-dimensional Su-Schrieffer-Heeger
model Yoshida  _et al._ (2014), the winding number based on zero-frequency
single-particle Green’s function is calculated by DMRG method in
distinguishing topologically non-trivial and trivial phases. By LDA+Gutzwiller
and LDA+DMFT methods, the $Z_{2}$ invariant has been applied in identifying
the correlated TIs SmB6 Lu  _et al._ (2013) and PuB6 Deng  _et al._ (2013).
The $Z_{2}$ invariant has also been calculated in QMC Hung  _et al._ (2013);
Lang  _et al._ (2013); Hung  _et al._ (2014); Meng  _et al._ (2014) and
cluster perturbation theory Grandi  _et al._ (2015a, b) for various
generalizations of the Kane-Mele-Hubbard model. Moreover, there are dynamical
mean-field theory calculations of Bernevig-Hughes-Zhang model with
interactions in which the $Z_{2}$ invariant is calculated Budich  _et al._
(2012, 2013); Amaricci  _et al._ (2015). As for spin Chern number constructed
from zero-frequency single-particle Green’s function, it has been applied to
verify the topological phase transitions in Kane-Mele-Hubbard model by quantum
Monte Carlo Hung  _et al._ (2014); Meng  _et al._ (2014) and cellular
dynamical mean-field theory (CDMFT) Chen  _et al._ (2015).

The issue of evaluating spin Chern number for interacting topological
insulators in CDMFT Chen  _et al._ (2015), or more generally quantum cluster
methods Maier  _et al._ (2005), is that the correlation effects can only be
captured inside the small cluster, even through the calculated spin Chern
number is quantized due to the mean-field bath at thermodynamic limit.
Nevertheless, such approaches cannot faithfully monitor the topological phase
transitions where a length scale larger than the cluster size is involved. On
the other hand, QMC is more accurate in capturing both short- and long-range
correlation effects by handling supercell with much larger size, but, in
obtaining topological invariants, it suffers from finite-size effect and the
topological invariants calculated from QMC for interacting TIs are not exactly
integer-quantized Hung  _et al._ (2013, 2014). Actually, due to finite size
effect, spin Chern number can be even far away from expected integer results
Hung  _et al._ (2014). Thus, it will be a great improvement if one can bring
the merit of CDMFT (its thermodynamic limit) into QMC to overcome the finite-
size effect in topological invariants, since the integer quantization of
topological invariants is essential for achieving well-defined topological
phases and identifying topological phase transitions. Here, we provide such a
scheme.

In this work, employing large-scale quantum Monte Carlo simulations for
interacting TIs, we eliminate the severe finite-size effect and obtain
quantized topological invariants by introducing a periodization scheme and
imposing symmetries of the studied systems. The numerical evaluation scheme of
both $Z_{2}$ invariant and spin Chern number constructed from zero-frequency
single-particle Green’s function proposed in Ref. Wang and Zhang (2012a) is
systematically presented with all important details. To demonstrate the
strength of our calculation scheme, especially the periodization process, we
test it on two independent 2D models of interacting TIs, in which the
topological phase transitions driven by one-body model parameters are detected
by the integer-quantized topological invariants.

The rest of the paper is organized as follows. In Sec. II, despite of the
already existing literature Wang  _et al._ (2010); Gurarie (2011); Wang  _et
al._ (2011, 2012); Wang and Zhang (2012b, a, 2014), the construction of the
$Z_{2}$ invariants and spin Chern number via zero-frequency single-particle
Green’s function is discussed for the sake of a self-contained narrative. A
brief introduction of projector Quantum Monte Carlo method is also presented.
Then our numerical calculation scheme and the periodization technique for
these topological invariants are introduced in detail in Sec. III. After that,
we show the applications of our scheme in Sec. IV, based on QMC simulations
for two independent 2D model systems of interacting TIs. Finally, Sec. V
summarizes our method and makes connection with the paper (II) in this series
on identifying interaction-driven topological phase transitions without
symmetry breaking by the topological invariants calculated via the scheme
presented here for interacting topological insulators, where the limitation of
topological invariants constructed from the single-particle Green’s function
is clearly manifested.

##  II Topological Invariants and Quantum Monte Carlo method

###  II.1 $Z_{2}$ invariant and spin Chern number

In interacting fermion systems, the single-particle Green’s function is given
as
$G(i\omega,\mathbf{k})=[i\omega+\mu-\mathcal{H}(\mathbf{k})-\Sigma(i\omega,\mathbf{k})]^{-1}$,
where $\Sigma(i\omega,\mathbf{k})$ is the self-energy originating from
interaction. The zero-frequency single-particle Green’s function is Hermitian
Wang and Zhang (2012a) and we can obtain its real eigenvalues by diagonalizing
the Green’s function matrix $G(i\omega=0,\mathbf{k})_{N_{O}\times N_{O}}$,
where $N_{O}$ is the number of orbitals or bands. In the language of Ref. Wang
and Yan (2013), the topological invariants can be defined as follows. We
simply define the so-called topological Hamiltonian
$h_{\textrm{t}}(\mathbf{k})=-G^{-1}(i\omega=0,\mathbf{k})_{N_{O}\times
N_{O}}$, and then calculate the topological invariant as if
$h_{\textrm{t}}(\mathbf{k})$ is a noninteracting Bloch Hamiltonian. The
advantage of $h_{\textrm{t}}(\mathbf{k})$ is that it reduces to the free Bloch
Hamiltonian in the non-interacting limit. Equivalently, we can work with
$G(i\omega=0,\mathbf{k})_{N_{O}\times N_{O}}$, since the eigenvectors of
$G(i\omega=0,\mathbf{k})_{N_{O}\times N_{O}}$ and
$-G^{-1}(i\omega=0,\mathbf{k})_{N_{O}\times N_{O}}$ are the same.

In this paper we study systems with both time-reversal symmetry and $U(1)$
spin-rotational symmetry, indicating $S_{z}$ conservation. Taking the $U(1)$
spin rotational symmetry into account, we can see that the Green’s function is
diagonal with respect to the spin index, and the two sub-matrices are denoted
as $G_{\sigma}$ ($\sigma=\uparrow,\downarrow$), each of which can be
diagonalized as
$G_{\sigma}(0,\mathbf{k})|\phi_{m}(0,\mathbf{k})\rangle=\mu_{m}(0,\mathbf{k})|\phi_{m}(0,\mathbf{k})\rangle$.
Both $Z_{2}$ invariant and spin Chern number can be simply constructed from
the eigenvectors $|\phi_{m}(0,\mathbf{k})\rangle$. For time-reversal invariant
and spatial-inversion-symmetric systems with interactions, the $Z_{2}$
invariant can be constructed from $|\phi_{m}(0,\mathbf{k})\rangle$ at TRIM
points. The $Z_{2}$ invariant can be expressed as Wang and Zhang (2012a)

$\displaystyle(-1)^{\nu}=\prod_{\bm{\kappa}\in\text{TRIM}}\prod_{\mu_{m}>0}\eta_{m}(\bm{\kappa}),$
(1)

with
$\eta_{m}(\bm{\kappa})=\langle\phi_{m}(0,\bm{\kappa})|\hat{P}|\phi_{m}(0,\bm{\kappa})\rangle$,
$\bm{\kappa}$ stands for TRIM points and $\hat{P}$ is the spatial inversion
symmetry operator. Here we have already taken the Kramer’s degeneracy at TRIM
points into consideration, and we only incorporates the parity in one spin
sector in Eq. (1). This expression of $Z_{2}$ invariant is a generalization of
that for free fermion system to interacting systems. Note that
$\\{|\phi_{m}\rangle\\}$ ($\mu_{m}>0$) reduces to the filled Bloch bands in
the noninteracting limit. Numerical evaluation of this $Z_{2}$ invariant in
correlated systems is quite simple, and it has been demonstrated Budich  _et
al._ (2013); Meng  _et al._ (2014); Hung  _et al._ (2013); Lang  _et al._
(2013); Hung  _et al._ (2014); Chen  _et al._ (2015) that this topological
invariant works well in detecting topological phase transition with change of
$Z_{2}$ invariant in weakly correlated systems.

So far we only discussed the $Z_{2}$ invariant. Due to the $U(1)$ spin-
rotational symmetry, there is actually a $Z$ invariant, which contains more
information. To introduce this $Z$ invariant, let us recall the TKNN or Chern
number of fermion systems with broken time-reversal symmetry, which has been
generalized to interacting fermion systems as Wang and Zhang (2012a)

$\displaystyle\centering\mathcal{C}=\frac{1}{2\pi}\iint_{\mathbf{k}\in
BZ}d^{2}\mathbf{k}\mathcal{F}_{xy}(\mathbf{k})\@add@centering$ (2)

with

$\displaystyle\mathcal{F}_{xy}(\mathbf{k})=\partial_{k_{x}}\mathcal{A}_{y}(\mathbf{k})-\partial_{k_{y}}\mathcal{A}_{x}(\mathbf{k})$
$\displaystyle\mathcal{A}_{i}(\mathbf{k})=-i\sum_{\mu_{m}>0}\langle\phi_{m}(0,\mathbf{k})|\partial_{k_{i}}|\phi_{m}(0,\mathbf{k})\rangle.$
(3)

Due to the $U(1)$ spin rotational symmetry, the Green’s function is diagonal
with respect to the spin index, thus we can calculate Chern numbers by Eq.
(II.1) for both spin-up and spin-down sectors denoted as $C_{\uparrow}$ and
$C_{\downarrow}$. Then the spin Chern number $C_{s}$ is simply defined as
$C_{s}=(C_{\uparrow}-C_{\downarrow})/2$. For 2D time-reversal invariant
systems, $C_{\uparrow}+C_{\downarrow}=0$, which results in the relation
$C_{s}=C_{\uparrow}$, and we only need to calculate Chern number
$C_{\uparrow}$ from the spin-up part of zero-frequency single-particle Green’s
function $G_{\uparrow}(0,\mathbf{k})$. This expression of spin Chern number
has also been applied in QMC Hung  _et al._ (2014) and CDMFT Chen  _et al._
(2015) simulations for interacting two-dimensional topological insulators (or
the “quantum spin Hall insulators (QSHI)”, in the older terminology).

Through a simple derivation from Eq. (2) and Eq. (II.1), we can arrive at an
expression that is numerically more convenient Avron  _et al._ (1983); Hung
_et al._ (2014):

$\displaystyle\mathcal{C}=$ $\displaystyle\frac{1}{2\pi i}\iint_{\mathbf{k}\in
BZ}dk_{x}dk_{y}\cdot$
$\displaystyle\textrm{Tr}\Big{\\{}P(\mathbf{k})\left[\partial_{k_{x}}P(\mathbf{k})\partial_{k_{y}}P(\mathbf{k})-\partial_{k_{y}}P(\mathbf{k})\partial_{k_{x}}P(\mathbf{k})\right]\Big{\\}},$

where $P(\mathbf{k})$ is a projection operator matrix constructed from
eigenvectors $|\phi_{m}(0,\mathbf{k})\rangle$ of $G_{\sigma}(0,\mathbf{k})$:

$\displaystyle
P(\mathbf{k})=\sum_{\mu_{m}>0}|\phi_{m}(0,\mathbf{k})\rangle\langle\phi_{m}(0,\mathbf{k})|.$
(5)

The systems we are dealing with in this paper are all based on honeycomb
lattice, the detailed implementation of Eq. (II.1) on the honeycomb lattice
geometry is presented in Appendix A.

Turning off interactions, we can observe that Eq. (1) reduces to the $Z_{2}$
invariant for free fermion systems defined by Fu and Kane Fu and Kane (2007),
and both Eq. (2) and Eq. (II.1) reduce to the TKNN invariant (or Chern number)
Thouless  _et al._ (1982); Avron  _et al._ (1983) for noninteracting systems.

In the next section, we introduce our numerical evaluation scheme of $Z_{2}$
invariant in Eq. (1) and spin Chern number in Eq. (II.1) in QMC simulations,
for model systems of interacting TIs with $U(1)_{spin}\times
U(1)_{charge}\rtimes Z_{2}^{T}$ symmetry and the spatial inversion symmetry.
Especially, both $Z_{2}$ invariant and spin Chern number for interacting TIs
are necessary to be quantized to achieve well-defined topological phases.
Thus, both periodization technique and implementation of symmetry properties
during numerical calculations to reach quantized topological invariants are
mainly introduced in next section.

###  II.2 Quantum Monte Carlo method

In this series of work, we apply the projector quantum Monte Carlo (PQMC)
simulation Meng  _et al._ (2014), which is the zero-temperature version of
determinantal QMC algorithm Assaad and Evertz (2008). PQMC method obtains the
ground-state observables by carrying out an imaginary time evolution starting
from trial wavefunction that has overlap to the true many-body ground state.
The ground-state expectation value of physical observable is calculated as
follows,

$\langle\hat{O}\rangle=\lim\limits_{\Theta\to+\infty}\frac{\langle\psi_{T}|e^{-\Theta\hat{H}/2}\hat{O}e^{-\Theta\hat{H}/2}|\psi_{T}\rangle}{\langle\psi_{T}|e^{-\Theta\hat{H}}|\psi_{T}\rangle},$
(6)

where $|\psi_{T}\rangle$ is the trial wave function and $\Theta$ is projection
parameter. In all the simulations, to ensure that the algorithm arrives at the
truly converged ground state of finite size systems, we choose $\Theta=40/t$
and $\Delta\tau=0.05/t$, in which $\Delta\tau$ is the finite imaginary time
step applied in the Trotter decomposition of partition function.

We can obtain both the static and dynamic observables. Static ones include the
expectation values of energy densities, double occupancy, and spin-spin
correlation function. As for the dynamic properties, we can measure the
dynamic single-particle Green’s function and spin-spin correlation function,
from which we can determine both the single-particle gap and spin gap for the
many-body systems. Especially, we concentrate on single-particle Green’s
function $G_{\sigma}(\mathbf{k},\tau)$ to calculate the topological
invariants. Generally, $G_{\sigma}(\mathbf{k},\tau)$ (in spin sector $\sigma$)
is defined as

$[G_{\sigma}(\tau,\mathbf{k})]_{pq}=\frac{1}{N}\sum_{i,j}e^{i\mathbf{k}\cdot(\mathbf{R}_{i}-\mathbf{R}_{j})}\langle
T_{\tau}[c_{ip,\sigma}(\tau)c^{\dagger}_{jq,\sigma}]\rangle,$ (7)

where $i,j\in[1,N]$ are the unit cell indices and $p,q\in[1,N_{o}]$ are the
orbital indices inside a unit cell. In this manner, for each $\mathbf{k}$
point, $G_{\sigma}(\tau,\mathbf{k})$ is a $N_{o}\times N_{o}$ hermitian matrix
(according to Eq. 12), and for the $L\times L$ system, there are $N=L^{2}$
momentum points. From $G_{\sigma}(\tau,\mathbf{k})$ data, we can directly
obtain the single-particle gap. Most importantly, we need to construct the
zero-frequency single-particle Green’s function
$G_{\sigma}(i\omega=0,\mathbf{k})$ from $G_{\sigma}(\tau,\mathbf{k})$ by
combining Fourier transformation and symmetry analysis presented in Sec. III,
after which the topological invariants can be calculated according to Eq. 1
and Eq. II.1.

##  III Numerical evaluation scheme of topological invariants

This section is divided into four successive parts. First, in Sec. III.1, we
comment on the condition for zero-frequency single-particle Green’s function
$G_{\sigma}(i\omega=0,\mathbf{k})$ to be well-behaved in both free and
interacting fermion systems. Second, in Sec. III.2, we explain how to obtain
correct $G_{\sigma}(i\omega=0,\mathbf{k})$ data from QMC simulation on a
finite-size system. Third, Sec. III.3 clarifies some numerical details in
evaluations of $Z_{2}$ invariant and spin Chern number. Finally, as the most
important part, the periodization scheme of $G_{\sigma}(i\omega=0,\mathbf{k})$
is introduced in Sec. III.4 to achieve the ideally quantized topological
invariants from QMC simulation in a finite-size system.

###  III.1 Condition for well-behaved $G(i\omega=0,\mathbf{k})$

In this work, we shall be concerned with quantum many-body systems at zero
temperature. Suppose that the system under consideration is gapped and has a
unique ground state under periodic boundary condition (no intrinsic
topological order), thus there is a many-body energy gap $\Delta$ between the
ground state energy level $E_{0}$ and the first excited energy level $E_{1}$,
namely, $\Delta=E_{1}-E_{0}$. In the Lehmann representation, we can see that
the retarded Green’s function $G(z,\mathbf{k})$ with complex frequency
variable $z=\omega_{R}+i\omega_{I}$ is an analytical function in the
$\omega_{R}\in(-\Delta,\Delta)$ region. In fact, the Lehmann representation
for $G(z,\mathbf{k})$ under zero temperature reads

$\displaystyle G_{\alpha\beta}(z,\mathbf{k})$ $\displaystyle=\sum_{m\neq
0}\left[\frac{\langle 0|c_{\mathbf{k}\alpha}|m\rangle\langle
m|c_{\mathbf{k}\beta}^{\dagger}|0\rangle}{z-(E_{m}-E_{0})}+\frac{\langle
m|c_{\mathbf{k}\alpha}|0\rangle\langle
0|c_{\mathbf{k}\beta}^{\dagger}|m\rangle}{z+(E_{m}-E_{0})}\right].\hskip
8.5359pt$ (8)

Since the condition $E_{m}-E_{0}\geq\Delta$ is satisfied, $G(z,\mathbf{k})$
has no poles in the $\omega_{R}\in(-\Delta,\Delta)$ region on the real axis.
We can further observe that $G(z,\mathbf{k})$ is actually an analytical
function with $\omega_{R}\in(-\Delta,\Delta)$ and arbitrary $\omega_{I}$.
Topological invariants are defined in terms of zero-frequency Green’s
function, namely, the Green’s function at $z=0$. More precisely, since the
chemical potential has been absorbed into the definition of $E_{m}$ (namely,
$E_{m}$ is the eigenvalue of $\hat{H}-\mu\hat{N}$, where $\hat{N}$ is the
particle number operator and $\mu$ is the chemical potential), the zero-
frequency refers to the energy exactly at the chemical potential Wang and
Zhang (2012a). Tuning chemical potential within the energy gap is harmless in
an insulator. Thus, the zero-frequency single-particle Green’s function is
well-defined for TIs with the above mentioned properties under zero-
temperature and one can construct topological invariants from it.

From both Eq. (1) and Eq. (II.1), we can determine that there are two
scenarios of the jumping of topological invariants determined in terms of
Green’s function. The first scenario is the pole of Green’s function, which is
the conventional case, and the topological transitions in noninteracting
fermion systems belong to this class. The second scenario is the zero of
Green’s function (namely, an eigenvalue of the Green’s function matrix becomes
zero). From the Lehmann representation, we can see that the first scenario
must describe a phase transition, since $E_{m}-E_{0}=0$ implies gap closing.
On the other hand, the second scenario can be a topological phase transition
Gurarie (2011); Yoshida  _et al._ (2014), but it is not necessarily so.

###  III.2 Calculations of $G(i\omega=0,\mathbf{k})$

In QMC, we measure the imaginary-time displaced Green’s function
$G(\tau,\mathbf{k})$ and obtain $G(i\omega_{n},\mathbf{k})$ by Fourier
transformation as follows

$\displaystyle G(i\omega_{n},\mathbf{k})=\int_{0}^{\beta}d\tau
e^{i\omega_{n}\tau}G(\tau,\mathbf{k}).$ (9)

Note that Eq. (9) has already incorporated the anti-periodic condition for
$G(\tau,\mathbf{k})$ as $G(\tau+\beta,\mathbf{k})=-G(\tau,\mathbf{k})$.
However, Eq. (9) is only valid for finite temperature. At zero-temperature,
$i\omega_{n}$ becomes continuous on the imaginary frequency axis and the anti-
periodic condition for $G(\tau,\mathbf{k})$ is not quite meaningful. In such
case, one needs to apply the following Fourier transformation

$\displaystyle G(i\omega,\mathbf{k})=\int_{-\infty}^{+\infty}d\tau
e^{i\omega\tau}G(\tau,\mathbf{k}).$ (10)

The validity of using Eq. (10) under zero-temperature is presented in Appendix
B. From Eq. (10), we can directly obtain the zero-frequency single-particle
Green’s function $G(i\omega,\mathbf{k})$ by substituting $i\omega=0$ as,

$\displaystyle G(i\omega=0,\mathbf{k})=\int_{-\infty}^{+\infty}d\tau
G(\tau,\mathbf{k}).$ (11)

One can furthermore make use of the symmetry properties of
$G(\tau,\mathbf{k})$ to simplify the calculation. For a general multi-band
systems with $N_{O}$ orbitals (in each spin sector within a unit cell), one
can prove

$\displaystyle[G_{\sigma}(\tau,\mathbf{k})]_{pq}=[G_{\sigma}(\tau,\mathbf{k})]_{qp}^{\star}$
(12)

with $p,q\in[1,N_{O}]$. Eq. (12) explicitly shows that
$G_{\sigma}(\tau,\mathbf{k})$ is a $N_{O}\times N_{O}$ Hermitian matrix. If
the system preserves spatial inversion symmetry and the corresponding spatial
inversion operation transforms $p$ sublattice to $p^{\prime}$ sublattice, then
we can prove

$\displaystyle[G_{\sigma}(\tau,\mathbf{k})]_{pq}=[G_{\sigma}(\tau,-\mathbf{k})]_{p^{\prime}q^{\prime}}.$
(13)

Eq. (13) explicitly connects the $G_{\sigma}(\tau,\mathbf{k})$ data with
opposite wavevector points. If the system preserves particle-hole symmetry, we
can prove

$\displaystyle[G_{\sigma}(\tau,\mathbf{k})]_{pq}=-\xi_{p}\xi_{q}[G_{\sigma}(-\tau,-\mathbf{k})]_{qp},$
(14)

where $\xi_{p}$ and $\xi_{q}$ are the sign change during on-site particle-hole
transformation as $c_{p}\to\xi_{p}c_{p}^{\dagger}$ and
$c_{q}\to\xi_{q}c_{q}^{\dagger}$. Eq. (14) shows the connections of
$G_{\sigma}(\tau,\mathbf{k})$ data with positive and negative $\tau$. The
detailed proof for Eq. (12), Eq. (13) and Eq. (14) is demonstrated in Appendix
C. Combining these three symmetry properties of the
$G_{\sigma}(\tau,\mathbf{k})$ matrix, we can determine that the number of non-
zero and independent matrix elements in $G_{\sigma}(\tau,\mathbf{k})$ is much
smaller than $N_{O}^{2}$.

###  III.3 Numerical details in evaluating topological invariants

Besides the discussions in Sec. III.1 and Sec. III.2, there are still some
important details in numerical application of Eq. (II.1) for calculating the
topological invariants.

Firstly, the infinite integral over imaginary time $\tau$ in Eq. (11) can be
approximated by a cutoff $\theta$ as

$\displaystyle G(i\omega=0,\mathbf{k})\approx\int_{-\theta}^{+\theta}d\tau
G(\tau,\mathbf{k}).$ (15)

For TIs, the systems have $G_{\sigma}(\tau,\mathbf{k})\propto
e^{-\Delta_{sp}(\mathbf{k})\tau}$ at large $\tau$ with
$\Delta_{sp}(\mathbf{k})$ as the single-particle gap at $\mathbf{k}$ point. If
the gap $\Delta_{sp}(\mathbf{k})$ is large, the exponential decay of
$G_{\sigma}(\tau,\mathbf{k})$ in imaginary time will be very fast, and a
finite $\theta$ is sufficient for the system to evolve below the energy scale
of $\Delta_{sp}$. However, this approximation can induce considerable error
around the topological phase transition with single-particle gap closing,
since the decaying of $G_{\sigma}(\tau,\mathbf{k})$ is very slow around the
transition point. Later on, one will observe the nonmonotonic behavior in
results of Chern number $C_{\uparrow}$ close to the topological phase
transition points in Sec. IV, which is originating from the above $\tau$
cutoff. We then calculate the integral in Eq. (15) numerically by a simple
trapezoidal method, and the step size for the this simple method is actually
the $\Delta\tau$ used in QMC simulations.

Secondly, the first-order derivatives over $k_{x},k_{y}$ in Eq. (II.1) are
replaced by first-order finite-difference Hung  _et al._ (2014). Since the
interval between two adjacent $\mathbf{k}$ points is proportional to $1/L$ for
a $L\times L$ system, replacing the derivatives by finite-difference will
bring in error proportional to $1/L$, which is the key origin for finite-size
effect in Chern number calculated by Eq. (II.1). Also, the integral over
$\mathbf{k}$ in BZ region in Eq. (II.1) can only be performed by summations
over $L^{2}$ discrete $\mathbf{k}$ points for a $L\times L$ system. We can
imagine this can also contribute to finite-size effect in Chern number
calculation by Eq. (II.1). Overall, due to the finite-size effect in QMC
simulations, the spin Chern number result can be very far from the expected
integer Hung  _et al._ (2014).

###  III.4 Periodization of $G(i\omega=0,\mathbf{k})$

As discussed above, the topological invariants calculated from QMC simulations
in finite-size system suffer severe finite-size effect, and the obtained
result is usually far away from the expected integer values. To remove the
finite-size effect, we combine the idea of periodization process as the last
step in CDMFT algorithm with QMC simulations. CDMFT simulates a correlated
system within a cluster (with some other noninteracting bath energy levels)
and open boundary condition is applied for the finite-size cluster. To recover
the translational symmetry broken by the open boundary condition, a
periodization process Parcollet  _et al._ (2004); Biroli  _et al._ (2004);
Sakai  _et al._ (2012); Li  _et al._ (2015) is applied to construct the
Green’s function and self-energy with arbitrary momentum resolution. The
periodization process is generally realized by

$\displaystyle
Q^{L}(i\omega_{n},\mathbf{k})=\frac{1}{N_{c}}\sum_{i,j=1}^{N_{c}}Q_{ij}^{C}(i\omega_{n})e^{i\mathbf{k}\cdot(\mathbf{r}_{i}-\mathbf{r}_{j})},$
(16)

where $N_{c}$ is the size of cluster, $Q^{C}(i\omega_{n})$ is the quantity
calculated on the cluster and $Q^{L}(\mathbf{k},i\omega_{n})$ is the
periodization result with arbitrary $\mathbf{k}$ point in BZ. In CDMFT,
generally one can choose $Q^{C}(i\omega_{n})$ to be $M^{C}(i\omega_{n})$ and
$\Sigma^{C}(i\omega_{n})$, where
$M^{C}(i\omega_{n})=[i\omega_{n}+\mu-\Sigma^{C}(i\omega_{n})]^{-1}$ is the
cluster cumulant and $\Sigma^{C}(i\omega_{n})$ is the cluster self-energy.
Naively thinking, directly substituting $G^{C}(i\omega_{n})$ into Eq. (16)
seems to be the simplest way, since matrix inverse operations are needed
during numerical calculations if $Q^{C}(i\omega_{n})$ is chosen to be
$M^{C}(i\omega_{n})$ and $\Sigma^{C}(i\omega_{n})$ instead. However, due to
the breaking of lattice translational symmetry, the simple $G$ periodization
can generate artificial results and another more complicated version of $G$
periodization instead of Eq. (16) should be applied to obtain reasonable
results Sakai  _et al._ (2012). For our purpose of calculating topological
invariants in QMC, nevertheless, we can apply the $G$ periodization method
simply via Eq. (16) by choosing $Q^{C}(i\omega_{n})$ to be
$G^{C}(i\omega_{n})$, since the QMC simulations respect periodic boundary
condition on finite-size system, which serves as the difference between the
periodization methods applied here and CDMFT Sakai  _et al._ (2012).

The periodization process is performed as follows. We first carry out QMC
simulation for $L\times L$ system with $N_{O}$ orbitals inside a unit cell
(for one spin sector). Then with Eq. (15) we obtain the
$G_{\sigma}(i\omega=0,\mathbf{k})$ data at $L^{2}$ discrete momentum points in
BZ. After that we carry out the Fourier transformation to get the real space
$[G_{\sigma,ij}(i\omega=0)]_{pq}$ data. Finally, we apply the periodization
process as Eq. (16) to obtain $G_{\sigma}(i\omega=0,\mathbf{k})$ data at
arbitrary $\mathbf{k}$ point in BZ as

$\displaystyle[G_{\sigma}(0,\mathbf{k})]_{pq}=\frac{1}{L^{2}}\sum_{i,j=1}^{L^{2}}[G_{\sigma,ij}(i\omega=0)]_{pq}e^{i\mathbf{k}\cdot(\mathbf{R}_{i}-\mathbf{R}_{j})},\hskip
14.22636pt$ (17)

where $i,j$ stands for unit cells and $p,q\in[1,N_{O}]$ are the indexes for
orbital or band. Throughout the periodization process, the
$G_{\sigma}(i\omega=0,\mathbf{k})$ data at the original $L^{2}$ wavevector
points keep unchanged. Based on this fact, we can take the periodization
process in Eq. (17) as Lagrange’s interpolation of
$G_{\sigma}(i\omega=0,\mathbf{k})$. The validity of the periodization process
is based on the observation that the topology of the system will not change
during a continuous deformation if there is no gap closing, which makes the
periodization process natural and appropriate. This means that if the
continuous deformation do not generate any singularity of the single-particle
Green’s function, the periodization process will not alter the topological
property of the system, but only improve the quality of the calculated
topological invariants and allow us to achieve quantized topological
invariants, especially the spin Chern number.

During numerical calculation, the lattice size for periodization
(interpolation) process will be large but still discrete, so we choose the
interpolation lattice size $IL\times IL$ and construct the
$G_{\sigma}(i\omega=0,\mathbf{k})$ data at $IL\times IL$ wavevector points in
BZ by Eq. (17). $IL$ can be much bigger than original lattice size $L$ in QMC
simulation. Then we evaluate the spin Chern number for the $IL\times IL$
system by Eq. (II.1). In next section, one can clearly observe the converging
of spin Chern number to the ideal quantized value with increasing $IL$.

##  IV Applications in QMC

In this section, we apply the above mentioned numerical scheme to two
independent 2D interacting TIs. The topological phase transitions in these
systems are driven by single-particle parameters. Results of topological
invariants across topological phase transitions in generalized Kane-Mele-
Hubbard model (GKMH) are shown in Sec. IV.1, while Sec. IV.2 concentrates on
topological phase transitions in cluster Kane-Mele-Hubbard model (CKMH).

We set the following parameters in the QMC simulations,
$\Theta=40/t,\Delta\tau=0.05/t$. For the imaginary time integration in Eq.
(15), a cutoff $\theta=20/t$ is applied.

###  IV.1 Generalized Kane-Mele-Hubbard model

Generalized Kane-Mele-Hubbard (GKMH) model Hung  _et al._ (2013, 2014); Meng
_et al._ (2014) is given by

$\displaystyle\hat{H}$ $\displaystyle=$ $\displaystyle-t\sum_{\langle
i,j\rangle\sigma}t_{ij}(c^{\dagger}_{i\sigma}c_{j\sigma}+h.c)-t_{3}\sum_{\langle\\!\langle\\!\langle
i,j\rangle\\!\rangle\\!\rangle\sigma}(c^{\dagger}_{i\sigma}c_{j\sigma}+h.c.)$
(18) $\displaystyle+i\lambda\sum_{\langle\\!\langle
i,j\rangle\\!\rangle\alpha\beta}v_{ij}(c^{\dagger}_{i\alpha}\sigma^{z}_{\alpha\beta}c_{j\beta}-c^{\dagger}_{j\beta}\sigma^{z}_{\beta\alpha}c_{i\alpha})$
$\displaystyle+\frac{U}{2}\sum_{i}(n_{i\uparrow}+n_{i\downarrow}-1)^{2},$

Figure 1: (Color online) (a) Illustration of the GKMH model as in Eq. (18).
The unit cell is presented as the yellow shaded rectangle, consisting of A and
B sublattices denoted by black and gray dots. The lattice is spanned by
primitive vectors $\mathbf{a}_{1}=(\sqrt{3},0)a$,
$\mathbf{a}_{2}=(1/2,\sqrt{3}/2)a$ with $a$ the lattice constant. The black
and blue lines denote nearest-neighbor hopping $t$ and $t_{d}$, while the
$\lambda$ term and third-nearest-neighbor hopping $t_{3}$ are represented by
red and green lines. The arrows in red lines shows $\nu_{ij}=+1$ for spin-up
part. (b) BZ of GKMH model. $\mathbf{K}_{1},\mathbf{K}_{2}$ are Dirac Points,
while
$\bm{\kappa}=\mathbf{\Gamma},\mathbf{M}_{1},\mathbf{M}_{2},\mathbf{M}_{3}$ are
the four TRIM points.

For nearest-neighbor (NN) hopping, we have $t_{ij}=t_{d}$ for NN bonds inside
unit cells and $t_{ij}=t$ for the others, as demonstrated in Fig. 1 (a). The
$t_{3}$ term is the third-nearest-neighbor hopping. The fourth term represents
spin-orbit coupling ($\lambda$) connecting next-nearest-neighbor sites with a
complex (time-reversal symmetric) hopping. The factor $\nu_{ij}=-\nu_{ji}=\pm
1$ depends on the orientation of the two nearest-neighbor bonds that the
electron moves in going from site $i$ to $j$. The last term describes the on-
site Coulomb repulsion. For GKMH model, we set $t$ as energy unit.

Figure 2: (Color online) (a) $(t_{d}/t)-(t_{3}/t)$ Phase diagram for GKMH
model at $U=0$ and arbitrary $\lambda/t>0$. (b) Calculation results of Chern
number $C_{\uparrow}$ with $t_{d}/t=0.5$ and varying $t_{3}/t$, the linear
system sizes are $L=12,24,48,72,120,400$.

Due to the $U(1)_{charge}\times U(1)_{spin}\rtimes Z_{2}^{T}$ symmetry of the
system, GKMH model acquires a $\mathbb{Z}$ classification. Without interaction
$(U=0)$, the $(t_{d}/t)-(t_{3}/t)$ phase diagram for $\lambda/t>0$, determined
from both $Z_{2}$ invariant $(-1)^{\nu}$ and Chern Number $C_{s}$, is shown in
Fig. 2 (a). There are four phases and their phase boundaries in Fig. 2 (a) is
independent of the size of $\lambda/t$, as long as $\lambda/t>0$. In Fig. 2
(b), we show the spin Chern number $C_{\uparrow}$ for fixed $t_{d}/t=0.5$ and
increasing $t_{3}/t$, calculated from Eq. (II.1). Since the system is
noninteracting, we have
$G_{\sigma}(i\omega=0,\mathbf{k})=-[\mathcal{H}_{\sigma}(\mathbf{k})]^{-1}$
with $\mathcal{H}_{\sigma}(\mathbf{k})$ the $2\times 2$ noninteracting
Hamiltonian matrix for GKMH model. In Fig. 2 (b), one can clearly observe that
the Chern number $C_{s}$ converges to its expected quantized value, with
increasing lattice size.

With finite interaction strength $U$, we need to measure the single-particle
Green’s function $G_{\sigma}(\tau,\mathbf{k})$ in QMC simulation and then
construct the topological invariants from it. For GKMH model,
$G_{\sigma}(\tau,\mathbf{k})$ is a $2\times 2$ Hermitian matrix, according to
Eq. (12). Due to the presence of spatial inversion symmetry and Eq. (13), we
can obtain
$[G_{\sigma}(\tau,\mathbf{k})]_{11}=[G_{\sigma}(\tau,-\mathbf{k})]_{22}$ and
$[G_{\sigma}(\tau,\mathbf{k})]_{12}=[G_{\sigma}(\tau,-\mathbf{k})]_{21}$. Then
at TRIM points $\bm{\kappa}$, we have
$[G_{\sigma}(\tau,\bm{\kappa})]_{11}=[G_{\sigma}(\tau,\bm{\kappa})]_{22}$ and
$[G_{\sigma}(\tau,\bm{\kappa})]_{12}=[G_{\sigma}(\tau,\bm{\kappa})]_{21}$.
From Eq. (14), we can get
$[G_{\sigma}(\tau,\mathbf{k})]_{pp}=-[G_{\sigma}(-\tau,-\mathbf{k})]_{pp},p=1,2$
and $[G_{\sigma}(\tau,\mathbf{k})]_{12}=[G_{\sigma}(-\tau,-\mathbf{k})]_{21}$.
Combining these relations, we obtain the equations applied in numerical
calculations as
$[G_{\sigma}(-\tau,\mathbf{k})]_{11}=-[G_{\sigma}(\tau,\mathbf{k})]_{22}$,
$[G_{\sigma}(-\tau,\mathbf{k})]_{12}=[G_{\sigma}(\tau,\mathbf{k})]_{12}$. At
TRIM points, $G_{\sigma}(\tau,\bm{\kappa})$ is a real symmetric matrix with
equal diagonal elements. Based on these considerations, we deduce that there
are only two independent matrix elements in $G_{\sigma}(\tau,\mathbf{k})$ as

$\displaystyle[G_{\sigma}(i\omega=0,\mathbf{k})]_{11}=-[G_{\sigma}(i\omega=0,\mathbf{k})]_{22}$
$\displaystyle\approx\int_{0}^{+\theta}\Big{\\{}[G_{\sigma}(\tau,\mathbf{k})]_{11}-[G_{\sigma}(\tau,\mathbf{k})]_{22}\Big{\\}}d\tau,$
(19)

and

$\displaystyle[G_{\sigma}(i\omega=0,\mathbf{k})]_{12}=[G_{\sigma}(i\omega=0,\mathbf{k})]_{21}^{\star}$
$\displaystyle\approx
2\int_{0}^{+\theta}[G_{\sigma}(\tau,\mathbf{k})]_{12}d\tau.$ (20)

As mentioned above, $\theta$ is the cutoff for $\tau$ in the integral of Eq.
(13). Eq. (IV.1) and Eq. (IV.1) explicitly show that we only need to measure
two elements of $G_{\sigma}(\tau,\mathbf{k})$ matrix with $\tau>0$ at all
discrete $\mathbf{k}$ points. At TRIM points, we can obtain the simplified
relations as
$[G_{\sigma}(i\omega=0,\bm{\kappa})]_{11}=[G_{\sigma}(i\omega=0,\bm{\kappa})]_{22}=0$
and
$[G_{\sigma}(i\omega=0,\bm{\kappa})]_{12}=[G_{\sigma}(i\omega=0,\bm{\kappa})]_{21}$.
Namely, at TRIM points we have
$G_{\sigma}(i\omega=0,\bm{\kappa})=Z_{\bm{\kappa}}\sigma_{x}$ while
$\sigma_{x}$ is the Pauli Matrix and $Z_{\bm{\kappa}}$ is some
$\bm{\kappa}$-dependent coefficient Lang  _et al._ (2013); Hung  _et al._
(2013, 2014); Meng  _et al._ (2014).

For $Z_{2}$ invariant, we only deal with $G_{\sigma}(i\omega=0,\bm{\kappa})$.
The matrix representation of inversion symmetry operator $\hat{P}$ for GKMH
model is $P=\sigma_{x}$ for each spin sector. The parity of all unoccupied
eigenstates at some TRIM point $\bm{\kappa}$ is simply
$\eta(\bm{\kappa})=\text{sgn}\\{[G_{\sigma}(0,\bm{\kappa})]_{12}\\}$. Then we
can get the $Z_{2}$ invariant simply as
$(-1)^{\nu}=\eta(\bm{\Gamma})\eta(\mathbf{M}_{1})\eta(\mathbf{M}_{2})\eta(\mathbf{M}_{3})$.
As a result, we observe that the $Z_{2}$ invariant for GKMH model is integer-
quantized, free from finite-size effect. On the other hand, to calculate spin
Chern number $C_{s}$, we only need to obtain
$[G_{\sigma}(\tau,\mathbf{k})]_{11}$ and $[G_{\sigma}(\tau,\mathbf{k})]_{12}$
in $G_{\sigma}(\tau,\mathbf{k})$ matrix due to its Hermiticity. After that,
spin Chern number can be numerically evaluated through applying Eq. (II.1).

Figure 3: (Color online) (a) $Z_{2}$ invariant $(-1)^{\nu}$ and (b), (c) Chern
number $C_{\uparrow}$ for the $t_{3}$-driven topological phase transition in
GKMH model with $t_{d}/t=1.0$, $\lambda/t=0.2$ and $U/t=2.0$ from finite-size
QMC simulation and the periodization results. (a) $Z_{2}$ invariant
$(-1)^{\nu}$ is quantized as mentioned in the main text. (b) and (c), the
Chern number $C_{\uparrow}$ from finite-size QMC calculation indicated by the
red open hexagon with error bar acquires a drop with finite value at the
transition point, which can be taken as signature of topological phase
transition, but the $C_{\uparrow}$ itself is not quantized before or after the
phase transition, due to the finite-size effect. After the periodization with
the QMC data in $L=12$ (b) and $L=18$ (c) systems, the $C_{\uparrow}$ converge
to the ideal quantized integers, where $IL$ stands for interpolation lattice
size used in the periodization process.

To apply our numerical calculation scheme for topological invariants, we
choose two paths in the phase diagram of Fig. 2 (a). First, starting from
$t_{d}/t=1.0$, $\lambda/t=0.2$ and $U/t=2.0$, we calculate the topological
invariants to monitor the $t_{3}$-driven topological phase transition. Second,
we choose $t_{3}/t=0$, $\lambda/t=0.2$, $U/t=2.0$ and calculate the
topological invariants to monitor the $t_{d}$-driven topological phase
transition. The interaction is chosen to be at a small value $U/t=2.0$ to
avoid the appearance of antiferromagnetic state Meng  _et al._ (2014); Wu  _et
al._ (2015). Thus, there is no spontaneous symmetry breaking across these two
topological phase transitions driven by hopping parameters in the interacting
GKMH model.

For the $t_{3}$-driven topological phase transition, the results of both
$Z_{2}$ invariant and Chern number $C_{\uparrow}$ (equal to spin Chern number
$C_{s}$) from QMC simulations with finite size $L=6,12,18$ are shown in Fig. 3
(a), (b) and (c), at $t_{d}/t=1.0,\lambda/t=0.2$ and $U/t=2.0$. At $U=0$, the
GKMH model experiences a topological phase transition at $t_{3}/t=1/3$ from
$C_{s}=+1$ to $C_{s}=-2$ as indicated in Fig. 2 (a). At $U/t=2.0$, the
quantized $Z_{2}$ invariant in Fig. 3 (a) demonstrates that the topological
phase transition point changes little. Detailed calculations of parities at
all TRIM points show that the parities at
$\mathbf{M}_{1},\mathbf{M}_{2},\mathbf{M}_{3}$ change across the phase
transition, which is related to the fact that the single-particle gap closes
at $\mathbf{M}_{1},\mathbf{M}_{2},\mathbf{M}_{3}$ points Meng  _et al._
(2014). As for the Chern number from finite-size QMC calculation denoted by by
the red open hexagon symbols with error bar in Fig. 3 (b) and (c), we can
observe drop with finite values across the phase transition, the position
coincide with that in $Z_{2}$ invariant. Combining these results, the phase
transition point is $t_{3}/t\approx 0.334\thicksim~{}0.335$ for $L=6$ system
and $t_{3}/t\approx 0.333\thicksim 0.334$ for $L=12,18$ system.

The problem about the results of Chern number in Fig. 3 (b) and (c) is that
they are not integer-quantized, though the convergence with increasing $L$ can
be observed. To solve this problem, the periodization method described in Sec.
III.4 is applied based on the results of Chern number from QMC simulations of
$L=12,18$. As shown in Fig. 3 (b) and (c), one indeed observes the gradual
convergence of Chern number $C_{\uparrow}$ when the interpolation lattice size
$IL$ increases. In fact, for $IL=120$, the Chern number $C_{\uparrow}$ is
almost ideally quantized, which demonstrates that the periodization method
works very well.

Figure 4: (Color online) (a) $Z_{2}$ invariant $(-1)^{\nu}$ and (b), (c) Chern
number $C_{\uparrow}$ for the $t_{d}$-driven topological phase transition in
GKMH model with $t_{3}/t=0$, $\lambda/t=0.2$ and $U/t=2.0$ from finite-size
QMC simulation and periodization results. (a) $Z_{2}$ invariant is integer-
quantized and there is a gradual shifting of the transition point from
$t_{d}/t\approx 1.955$ in $L=6$ system to $t_{d}/t\approx 1.995$ in $L=18$
system. (b) and (c), the finite drop in Chern number from finite-size QMC
calculation with $L=12$ and $L=18$ (denoted by red open hexagon with error
bar) can be observed clearly. The periodization results of Chern number
$C_{\uparrow}$, with the QMC data in $L=12$ (b) and $L=18$ (c) systems, are
ideally quantized.

For the $t_{d}$-driven topological phase transition, the results of both
$Z_{2}$ invariant and Chern number $C_{\uparrow}$ (equal to spin Chern number
$C_{s}$) from finite-size QMC simulations are shown in Fig. 4, at $t_{3}/t=0$,
$\lambda/t=0.2$ and $U/t=2.0$. The noninteracting GKMH model obtains a
topological phase transition at $t_{d}/t=2.0$ for $t_{3}/t=0$, as long as
$\lambda/t>0$. From the results in Fig. 4, one sees the weak interaction
$U/t=2.0$ only give a small shift of the topological phase transition point.
Across this $t_{d}$-driven phase transition, both the parity change and
single-particle gap close happen only at the $\mathbf{M}_{2}$ point, due to
the anisotropy Hung  _et al._ (2013, 2014); Meng  _et al._ (2014) introduced
by $t_{d}$. From the integer-quantized $Z_{2}$ invariant, the phase transition
point only possesses a small shift in $t_{d}/t$, from $L=6$ system to $L=18$
system, in Fig. 4 (a). By means of the periodization process, with the QMC
simulation results of $L=12$ and $18$, we obtain the perfect, quantized Chern
number $C_{\uparrow}$ results in Fig. 4 (b) and (c).

###  IV.2 Cluster Kane-Mele-Hubbard model

The cluster Kane-Mele-Hubbard model Grandi  _et al._ (2015a); Wu  _et al._
(2012) (CKMH) has six honeycomb lattice sites as one unit cell, the
Hamiltonian is given as follows

$\displaystyle\hat{H}=$ $\displaystyle-\sum_{\langle
ij\rangle\sigma}t_{ij}(c^{\dagger}_{i\sigma}c_{j\sigma}+c^{\dagger}_{j\sigma}c_{i\sigma})$
(21) $\displaystyle+i\lambda_{I}\sum_{\langle\\!\langle
ij\rangle\\!\rangle\alpha\beta}v_{ij}(c^{\dagger}_{i\alpha}\sigma^{z}_{\alpha\beta}c_{j\beta}-c^{\dagger}_{j\beta}\sigma^{z}_{\beta\alpha}c_{i\alpha})$
$\displaystyle+i\lambda_{O}\sum_{\langle\\!\langle
ij\rangle\\!\rangle\alpha\beta}v_{ij}(c^{\dagger}_{i\alpha}\sigma^{z}_{\alpha\beta}c_{j\beta}-c^{\dagger}_{j\beta}\sigma^{z}_{\beta\alpha}c_{i\alpha})$
$\displaystyle+\frac{U}{2}\sum_{i}(n_{i\uparrow}+n_{i\downarrow}-1)^{2}\;.$

For nearest-neighbor (NN) hopping, we have $t_{ij}=t$ for NN bonds inside unit
cells and $t_{ij}=t_{d}$ for those connecting the six-site unit cells, as
demonstrated in Fig. 5 (a). The amplitudes for SOC term inside a unit cell and
between different unit cells are $\lambda_{I}$ and $\lambda_{O}$,
respectively. $U$ is the on-site Coulomb repulsion. Similar to the GKMH model,
the $U(1)_{change}\times U(1)_{spin}\rtimes Z_{2}^{T}$ symmetry is also
preserved in CKMH model, which results in $\mathbb{Z}$ classification.
Besides, both spatial inversion symmetry and particle-hole symmetry are also
present in CKMH model. Notice that the CKMH model has 6-site unit cell, in
this section, the linear system size $L$ in finite size QMC simulation
actually corresponds to that of the a 6-site unit cell, i.e., the total
lattice sites are $6\times L\times L$. Such large unit cell greatly increases
the QMC simulation efforts of CKMH model comparing to that of the GKMH model
in previous session, where the total lattice site is only $2\times L\times L$
(the computation efforts of QMC scale to the third power of the total lattice
sites).

Figure 5: (Color online) Illustration of the CKMH model in Eq. (21). The
yellow shaded region shows the six-site unit cell with primitive lattice
vectors $\mathbf{a}_{1}=(\sqrt{3},0)a$, $\mathbf{a}_{2}=(1/2,\sqrt{3}/2)a$
with the nearest-neighbor bond length $a/\sqrt{3}$. The six sublattices
$1,2,3,4,5,6$ are shown in different color. The black solid and black dotted
lines indicate the nearest-neighbor hopping term inside ($t$) and between
($t_{d}$) unit cells. The red solid and red dotted lines represent SOC terms
inside ($\lambda_{I}$) and between ($\lambda_{O}$) unit cells. The sign choice
for SOC hopping is the same as that in Fig. 1 (a). The on-site Coulomb
repulsion $U$ is shown by the blue shaded circle. Figure 6: (Color online)
Phase diagrams for noninteracting CKMH model under two different parameter
sets. (a) $(t_{d}/t)-(\lambda_{I}/t)$ phase diagram with $\lambda_{O}=0$; (b)
$\alpha-(\lambda_{I}/t)$ phase diagram with
$\alpha=t_{d}/t=\lambda_{O}/\lambda_{I}$.

Since there are three independent parameters $t_{d},\lambda_{I},\lambda_{O}$,
even the phase diagram for noninteracting CKMH model is already interesting.
To simplify the presentation, we only demonstrate results on two special
cases. First, we set $\lambda_{O}=0$, the $(t_{d}/t)-(\lambda_{I}/t)$ phase
diagram of this case is shown in Fig. 6 (a). We can observe that the QSHI
phases exist in the middle region of the $(t_{d}/t)-(\lambda_{I}/t)$ phase
diagram, with different spin Chern numbers $C_{s}=+1$ and $C_{s}=-2$, as a
function of $\lambda_{I}/t$. The presence of the $C_{s}=-2$ phase is
unexpected and very interesting, since the Chern number is changed by 3 when
going from $C_{s}=1$ to $C_{s}=-2$, further increase $\lambda_{I}/t$, QSHI
phase is destroyed. Second, we keep all three $t_{d},\lambda_{I},\lambda_{O}$
parameters finite and introduce a ratio of hopping
$\alpha=t_{d}/t=\lambda_{O}/\lambda_{I}$. The $\alpha-(\lambda_{I}/t)$ phase
diagram at $U=0$ is presented in Fig. 6 (b). In the plotted region of $\alpha$
and $t_{d}/t$, three phases with one non-trivial in $C_{s}$ are found.

With interaction, we apply the periodization process in Sec. III.4 to
calculate the topological invariants. Before presenting the data, let’s
discuss the structure of $G_{\sigma}(\tau,\mathbf{k})$ and
$G_{\sigma}(i\omega=0,\mathbf{k})$ for CKMH model as it is quite complicated.
First, both $G_{\sigma}(\tau,\mathbf{k})$ and
$G_{\sigma}(i\omega=0,\mathbf{k})$ are $6\times 6$ Hermitian matrices for CKMH
model. Second, using spatial inversion and particle-hole symmetries, we can
obtain useful relations among the matrix elements of
$G_{\sigma}(\tau,\mathbf{k})$. Combining these two symmetry properties,
$G_{\sigma}(\tau,\mathbf{k})$ and $G_{\sigma}(-\tau,\mathbf{k})$ are
explicitly related and we only need to calculate the
$G_{\sigma}(\tau,\mathbf{k})$ data with $\tau>0$. The detailed analysis is
presented in Appendix D. The $G_{\sigma}(i\omega=0,\mathbf{k})$ matrix for the
CKMH model only has 12 independent matrix elements and can be expressed as

$\displaystyle
G_{\sigma}(0,\mathbf{k})=\left(\begin{array}[]{cccccc}A_{1}&A_{4}&A_{5}&A_{6}&A_{7}&A_{8}\\\
A_{4}^{\star}&A_{2}&A_{9}&-A_{7}&A_{10}&A_{11}\\\
A_{5}^{\star}&A_{9}^{\star}&A_{3}&A_{8}&-A_{11}&A_{12}\\\
A_{6}^{\star}&-A_{7}^{\star}&A_{8}^{\star}&-A_{1}&A_{4}^{\star}&-A_{5}^{\star}\\\
A_{7}^{\star}&A_{10}^{\star}&-A_{11}^{\star}&A_{4}&-A_{2}&A_{9}^{\star}\\\
A_{8}^{\star}&A_{11}^{\star}&A_{12}^{\star}&-A_{5}&A_{9}&-A_{3}\end{array}\right),\hskip
11.38092pt$ (28)

where $A_{1},A_{2},A_{3}$ are real numbers and $A_{i},i=4,\cdots,12$ are
complex numbers. At TRIM points, Eq. (28) can be further simplified to only 9
independent matrix elements and $G_{\sigma}(i\omega=0,\bm{\kappa})$ matrix
obtains the following matrix structure as

$\displaystyle G_{\sigma}(0,\bm{\kappa})=$
$\displaystyle\left(\begin{array}[]{cccccc}0&B_{4}&iB_{5}&B_{6}&iB_{7}&B_{8}\\\
B_{4}&0&B_{9}&-iB_{7}&B_{10}&iB_{11}\\\
-iB_{5}&B_{9}&0&B_{8}&-iB_{11}&B_{12}\\\ B_{6}&iB_{7}&B_{8}&0&B_{4}&iB_{5}\\\
-iB_{7}&B_{10}&iB_{11}&B_{4}&0&B_{9}\\\
B_{8}&-iB_{11}&B_{12}&-iB_{5}&B_{9}&0\end{array}\right).\hskip 14.22636pt$
(35)

In Eq. (IV.2), $B_{i},i=4,\cdots,12$ are real numbers. The explicit
expressions for both $A_{i}$ and $B_{i}$ are given in Appendix D.

For numerical evaluation of $Z_{2}$ invariant in CKMH model, we adopt the
$G_{\sigma}(i\omega=0,\bm{\kappa})$ matrix in Eq. (IV.2) for TRIM points. The
matrix representation of spatial inversion symmetry operator for each spin
sector is

$\displaystyle P=\left(\begin{array}[]{cccccc}0&0&0&1&0&0\\\ 0&0&0&0&1&0\\\
0&0&0&0&0&1\\\ 1&0&0&0&0&0\\\ 0&1&0&0&0&0\\\ 0&0&1&0&0&0\end{array}\right).$
(42)

To calculate the $Z_{2}$ invariant, we only need to diagonalize the
$G_{\sigma}(i\omega=0,\bm{\kappa})$ matrix with
$\bm{\kappa}=\bm{\Gamma},\mathbf{M}_{1},\mathbf{M}_{2},\mathbf{M}_{3}$, and
then calculate the parities at these TRIM points by Eq. (1). To calculate spin
Chern number, we only need to adopt the matrix structure of
$G_{\sigma}(i\omega=0,\mathbf{k})$ in Eq. (28) and apply Eq. (II.1), first
with the finite-size QMC data and then with periodization process.

Figure 7: (Color online) (a) Z2 invariant $(-1)^{\nu}$ and (b), (c) Chern
number $C_{\uparrow}$ for the $\lambda_{I}$-driven topological phase
transition with $\lambda_{O}=0$, $t_{d}/t=0.5$ and $U/t=1.5$ from finite-size
calculations by QMC simulation (denoted by red open hexagon with error bar)
and periodization. The drop of integer-valued $Z_{2}$ invariant and the jump
of Chern Number $C_{\uparrow}$ can be taken as signature of topological phase
transition. In (b) and (c), with large interpolation lattice size $IL$ from
finite QMC data with $L=6$ and $12$, the Chern number $C_{\uparrow}$ reaches
its quantized value.

We concentrate on two independent paths in the parameter space of CKMH model.
First, we choose $\lambda_{O}=0$, $t_{d}/t=0.5$ and $U/t=1.5$ and study the
$\lambda_{I}$-driven topological phase transition in the interacting CKMH
model. For $U=0$, the transition point for this $\lambda_{I}$-driven
topological phase transition from $C_{s}=0$ to $C_{s}=+1$ is at
$\lambda_{I}/t\approx 0.289$. Second, we set $\alpha=1.8,U/t=2.0$ in the CKMH
model and study the $\lambda_{I}$-driven topological phase transition,
similarly, there is also a $\lambda_{I}$-driven topological phase transition
from $C_{s}=0$ to $C_{s}=+1$ at $\lambda_{I}/t\approx 0.1$ for $U=0$. In the
following, we calculate $Z_{2}$ invariant $(-1)^{\nu}$ and spin Chern number
$C_{s}$ for $\lambda_{I}/t$ parameter across these two phase transitions, to
demonstrate that our periodization process works for CKMH model as well.

Fig. 7 shows the $Z_{2}$ invariant $(-1)^{\nu}$ and Chern number
$C_{\uparrow}$ (equal to spin Chern number $C_{s}$), calculated from finite-
size QMC simulations of $L=6,12$ systems (both denoted by red open hexagon
with error bar) for the $\lambda_{O}=0$, $t_{d}/t=0.5$ and $U/t=1.5$ case. One
can observe that the $Z_{2}$ invariant is exactly integer-quantized to $\pm 1$
(Fig. 7 (a)). Such integer quantization in finite-size system is due to the
spatial inversion and particle-hole symmetries of the CKMH model. These two
symmetries result in the special matrix structure of
$G_{\sigma}(i\omega=0,\bm{\kappa})$ in Eq. (IV.2), which is sufficient to
guarantee that the parities at TRIM points are exactly $\pm 1$. As shown in
Fig. 7 (a), (b) and (c), for both $L=6$ and $L=12$, the $\lambda_{I}$-driven
topological phase transition point is at $\lambda_{I}/t=0.288\thicksim 0.289$
from the drop in $Z_{2}$ invariant $(-1)^{\nu}$ and the jump in Chern number
$C_{\uparrow}$. Across the topological phase transition, both single-particle
gap close and the parity change all happen at $\bm{\Gamma}$ point. As for the
Chern number $C_{\uparrow}$ from finite size $L=6$ and $L=12$ QMC calculation
(Fig. 7 (b) and (c)), they are still far from the ideal quantized result due
to the finite-size effect, although the trend of convergence with increasing
system size is present.

Figure 8: (Color online) (a) Z2 invariant $(-1)^{\nu}$ and (b), (C) Chern
number $C_{\uparrow}$ for the $\lambda_{I}$-driven topological phase
transition with $\alpha=1.8$ and $U/t=2.0$ from finite-size calculations by
QMC (denoted by red open hexagon with error bar) and periodization. The drop
of integer-valued $Z_{2}$ invariant and the jump of Chern number
$C_{\uparrow}$ can be taken as signature of topological phase transition. In
(b) and (c), with large $IL$ in periodization, the Chern number $C_{\uparrow}$
reaches its quantized value.

The calculation results of Chern number $C_{\uparrow}$ by the periodization
method for the $\lambda_{I}$-driven topological phase transition is presented
in Fig. 7 (b) and (c), from the QMC data in $L=6$ and $12$ systems. We can
observe that the periodization method can give good results of quantized
integer for $C_{\uparrow}$. The nonmonotonic behavior in $IL=12,24,48$ close
to the transition in Fig. 7 (b) is due to the inappropriate $\tau$ cutoff in
calculating Eq. (15), as mentioned in Sec. III.3. This behavior is absent for
large enough $IL$s. Through this periodization method, a sharp topological
phase transition from $C_{s}=0$ to $C_{s}=+1$ can be clearly seen.

Fig. 8 shows the results of $Z_{2}$ invariant $(-1)^{\nu}$ (Fig. 8 (a)) and
Chern number $C_{\uparrow}$ (Fig. 8 (b), (c)) for the case of
$\lambda_{I}$-driven topological phase transition with $\alpha=1.8$, $U/t=2$,
from finite-size QMC simulations of $L=6$, $12$ systems and periodization. The
sharp drop of integer-valued $Z_{2}$ invariant $(-1)^{\nu}$ defines the
topological phase transition, at $\lambda_{I}/t=0.103\thicksim 0.104$ for
$L=6$ system and $\lambda_{I}/t=0.100\thicksim 0.101$ for $L=12$ system. The
positions for the finite value jump of Chern number $C_{\uparrow}$ in $L=6$
and $L=12$ systems are consistent with those of $Z_{2}$ invariant
$(-1)^{\nu}$. Across the topological phase transition, both the single-
particle gap close and the parity change happen at $\bm{\Gamma}$ point. Still,
$C_{\uparrow}$ from QMC are not quantized integer value according to the
results in Fig. 8 (b) and (c).

The calculation results of Chern number $C_{\uparrow}$ after the periodization
for the $\lambda_{I}$-driven topological phase transition is presented in Fig.
8 (b) and (c), from the QMC data in $L=6$ and $12$ systems. With large enough
$IL$, the integer-valued Chern number $C_{\uparrow}$ can be obtained. Again,
very close to the transition point, the nonmonotonic behavior of
$C_{\uparrow}$ from periodization method appears both in Fig. 8 (b) and (c),
when $IL$ is not large, but it disappears after we increase $IL$.

##  V Conclusion

To conclude, we provide a toolkit to calculate the topological invariants
constructed from zero-frequency single-particle Green’s function for
interacting TIs. All the important numerical details are carefully documented,
hence serves the purpose of demystifying the numerical evaluation of $Z_{2}$
invariant and spin Chern number for interacting TIs. Most importantly, we
introduce a periodization process to eliminate the finite-size effect on spin
Chern number and obtain quantized topological invariants from finite-size QMC
simulations, which renders the topological phases well-defined.

To demonstrate the power of our calculation scheme, especially the
periodization process, both the topological phases and topological phase
transitions in two interacting TI models, namely, the Generalized Kane-Mele-
Hubbard model and the cluster Kane-Mele-Hubbard model, are identified by
numerical evaluation of topological invariants. The results show that the
numerical scheme works well in capturing the topological phases and their
phase transitions driven by one-body model parameter. The $Z_{2}$ invariant
are already integer-quantized by applying the symmetry properties of the
studied system during data process. With periodization process, the integer-
quantized spin Chern number is also achieved with QMC simulations in systems
with very small size. Through these calculations of topological invariants, we
can also determine the topological phase transition points accurately, at
least more accurate than those from the gap extrapolations.

The present work demonstrates that the numerical evaluation scheme, especially
the periodization process, of topological invariants for interacting TIs works
well, in distinguishing topologically phases and identifying their phase
transitions driven by the one-body model parameters. In paper (II) of this
series, we shall apply the numerical evaluation scheme developed here to wider
classes of models of interacting TIs, in which certain limitation of
constructing topological invariants from single-particle Green’s function is
manifested in a very interesting manner, calling for more versatile technique
to diagnose the interaction-driven topological phase transitions in
interacting TIs.

###### Acknowledgements.

We thank Ning-Hua Tong, Lei Wang, Yi-Zhuang You, Cenke Xu, Liang Fu and Xi Dai
for inspiring discussions, and in particular we acknowledge Zhong Wang for
helpful comments on the Green’s function formalism and the manuscript. The
numerical calculations were carried out at the Physical Laboratory of High
Performance Computing in Renmin University of China as well as the National
Supercomputer Center in TianJin on the platform Tianhe-1A. YYH, HQW and ZYL
acknowledge support from National Natural Science Foundation of China (Grant
Nos. 11474356 and 11190024) and National Program for Basic Research of MOST of
China (Grant No. 2011CBA00112). ZYM is supported by the National Natural
Science Foundation of China (Grant Nos. 11421092 and 11574359) and the
National Thousand-Young-Talents Program of China, and acknowledges the
hospitality of the KITP at the University of California, Santa Barbara, where
part of this work is completed.

##  Appendix A Detailed implementation of Eq. (II.1) for honeycomb lattice

In this appendix, we present the numerical implementation of Eq. (II.1) on
honeycomb lattice, provided that we have already obtained the zero-frequency
single-particle Green’s function data $G_{\sigma}(i\omega=0,\mathbf{k})$. The
reason for having this discussion is that the primitive vectors for honeycomb
lattice in real $(\mathbf{a}_{1},\mathbf{a}_{2})$ and momentum space
$(\mathbf{b}_{1},\mathbf{b}_{2})$ (see Fig. 1 (a), (b)) are not along the
$(x,y)$ or $(k_{x},k_{y})$ direction of the Cartesian coordinate. Hence, both
the integral over honeycomb lattice BZ and the derivatives over
$(k_{x},k_{y})$ in Eq. (II.1) can be Jacobian transformed from those on the
honeycomb lattice BZ to that on a square lattice as

$\displaystyle\left\\{\begin{array}[]{ll}q_{u}=k_{x}\\\
q_{v}=\frac{k_{x}+\sqrt{3}k_{y}}{2}\end{array}\right.\hskip
22.76228pt\left\\{\begin{array}[]{ll}0\leq q_{u}\leq\frac{2\pi}{\sqrt{3}a}\\\
0\leq q_{v}\leq\frac{2\pi}{\sqrt{3}a}\end{array}\right.\hskip
19.91684pt\Longrightarrow\hskip 19.91684ptk_{x}=q_{u}\hskip
19.91684ptk_{y}=\frac{-q_{u}+2q_{v}}{\sqrt{3}}.$ (47)

The transformed BZ is indeed a square one, and we can rewrite Eq. (II.1) by
substituting Eq. (47) as

$\displaystyle\mathcal{C}$ $\displaystyle=$ $\displaystyle\frac{1}{2\pi
i}\int_{0}^{2\pi/\sqrt{3}a}dq_{u}\int_{0}^{2\pi/\sqrt{3}a}dq_{v}\cdot\text{
Tr}\Big{\\{}Q(q_{u},q_{v})\Big{[}\partial_{q_{u}}Q(q_{u},q_{v})\partial_{q_{v}}Q(q_{u},q_{v})-\partial_{q_{v}}Q(q_{u},q_{v})\partial_{q_{u}}Q(q_{u},q_{v})\Big{]}\Big{\\}},$
(48)

where $Q(q_{u},q_{v})=P(q_{u},(-q_{u}+2q_{v})/\sqrt{3})$ and $P(k_{x},k_{y})$
is the projection matrix defined in Eq. (5). Comparing Eq. (II.1) and Eq.
(48), one observes the Jacobian transformation does not change the form of the
formula.

QMC simulates finite size system, so the integration and derivation in Eq.
(48) are discretized, we set $Q_{i,j}=Q(q_{u,i},q_{v,j})$ and
$(q_{u,i},q_{v,j})=(2i\pi/\sqrt{3}L_{1}a,2j\pi/\sqrt{3}L_{2}a)$ with
$i\in[0,L_{1}],j\in[0,L_{2}]$, then we have the expressions for finite
difference as

$\displaystyle\partial_{q_{u}}Q(q_{u},q_{v})|_{\mathbf{q}=(q_{u,i},q_{v,j})}=\frac{Q_{i+1,j}-Q_{i-1,j}}{2\delta_{q_{u}}}\hskip
28.45274pt\partial_{q_{v}}Q(q_{u},q_{v})|_{\mathbf{q}=(q_{u,i},q_{v,j})}=\frac{Q_{i,j+1}-Q_{i,j+1}}{2\delta_{q_{v}}},$
(49)

where $\delta_{q_{u}}=2\pi/\sqrt{3}L_{1}a,\delta_{q_{v}}=2\pi/\sqrt{3}L_{2}a$.
Due to the periodic boundary condition, we have
$Q_{0,j}=Q_{L_{1},j},Q_{L_{1}+1,j}=Q_{1,j}$ and
$Q_{i,0}=Q_{i,L_{2}},Q_{i,L_{2}+1}=Q_{i,1}$. Based on Eq. (49), we arrive at
the expression for the integrand in Eq. (48) as

$\displaystyle\text{Tr}\Big{\\{}$ $\displaystyle
Q(q_{u},q_{v})\Big{(}\partial_{q_{u}}Q(q_{u},q_{v})\partial_{q_{v}}Q(q_{u},q_{v})-\partial_{q_{v}}Q(q_{u},q_{v})\partial_{q_{u}}Q(q_{u},q_{v})\Big{)}\Big{\\}}$
$\displaystyle=\frac{1}{4\delta_{q_{u}}\delta_{q_{v}}}\text{Tr}\Big{\\{}Q_{i,j}\Big{(}[Q_{i+1,j},Q_{i,j+1}]+[Q_{i,j+1},Q_{i-1,j}]+[Q_{i-1,j},Q_{i,j-1}]+[Q_{i,j-1},Q_{i+1,j}]\Big{)}\Big{\\}}.$

Simultaneously, the integral over the square BZ in Eq. (48) is changed into
summation over discrete wavevector points as

$\int_{\mathbf{k}\in
BZ}f(\vec{k})d^{2}\mathbf{k}=\frac{\Omega_{BZ}}{L_{1}L_{2}}\sum_{\mathbf{k}\in
BZ}f(\mathbf{k}),$ (51)

where $L_{1}L_{2}$ is number of unit cells for the finite size system, and
$\Omega_{BZ}=4\pi^{2}/3a$ is the volume of BZ. For a finite-size system, the
result of the summation will deviate from the expected quantized integer,
which is the finite-size effect we have seen in the main text, i.e., Fig. 3,
Fig. 4, Fig. 7 and Fig. 8, but we have also seen that the summation results
converge to the quantized integer with increasing system size. Combining Eq.
(A) and Eq. (51), we deduce the constant coefficient for Eq. (48) as

$\displaystyle\frac{1}{2\pi
i}\cdot\frac{1}{4\delta_{q_{u}}\delta_{q_{v}}}\cdot\frac{\Omega_{BZ}}{N}=\frac{1}{2\pi
i}\cdot\frac{3L_{1}L_{2}a^{2}}{16\pi^{2}}\cdot\frac{4\pi^{2}}{3L_{1}L_{2}a^{2}}=\frac{1}{8\pi
i}.$ (52)

So finally, we can get the simplified expression of Eq. (48) as

$\displaystyle\mathcal{C}=$ $\displaystyle\frac{1}{8\pi
i}\sum_{i=1}^{L_{1}}\sum_{j=1}^{L_{2}}S(q_{u,i},q_{v,j})$ $\displaystyle
S(q_{u,i},q_{v,j})=\text{Tr}\Big{\\{}Q_{i,j}\Big{(}[Q_{i+1,j},Q_{i,j+1}]$
$\displaystyle+[Q_{i,j+1},Q_{i-1,j}]+[Q_{i-1,j},Q_{i,j-1}]+[Q_{i,j-1},Q_{i+1,j}]\Big{)}\Big{\\}}.$
(53)

During the calculation, we only need to prepare the projection matrix
$P(q_{u},(-q_{u}+2q_{v})/\sqrt{3})$ through the zero-frequency Green’s
function matrix obtained from QMC simulation. As Chern number should be an
integer for a gapped system, we can only calculate the imaginary part of
$S(q_{u,i},q_{v,j})$ in Eq. (A).

##  Appendix B Validity of Eq. (10) at zero-temperature

In this appendix, we validate the usage of Eq. (10) in obtaining
$G(i\omega,\mathbf{k})$ data from $G(\tau,\mathbf{k})$ at zero temperature.

Let’s start with a review on some basic facts about Matsubara frequency
Green’s function at finite-temperature. The Matsubara frequencies,
$i\omega_{n}=i(2n+1)\pi/\beta,n\in Z$ for fermion systems and
$i\omega_{n}=i2n\pi/\beta,n\in Z$ for boson systems, are actually the poles of
the corresponding Fermi-Dirac and Bose-Einstein distribution function as

$\displaystyle
n_{FD}(\varepsilon)=\frac{1}{e^{\beta\varepsilon}+1}=\frac{1}{2}+\frac{1}{\beta}\sum_{n=-\infty}^{+\infty}\frac{1}{i\omega_{n}-\varepsilon}\hskip
28.45274pt\omega_{n}=\frac{(2n+1)\pi}{\beta}$ (54) $\displaystyle
n_{BE}(\varepsilon)=\frac{1}{e^{\beta\varepsilon}-1}=-\frac{1}{2}-\frac{1}{\beta}\sum_{n=-\infty}^{+\infty}\frac{1}{i\omega_{n}-\varepsilon}\hskip
28.45274pt\omega_{n}=\frac{2n\pi}{\beta}.$

We can explicitly observe that the exact zero Matsubara Frequency $i\omega=0$
can only be physically reachable at zero temperature. Otherwise, the
imaginary-time single-particle Green’s function at finite temperature is
defined as

$\displaystyle
G(\tau,\hat{A}\hat{B})=-\Big{\langle}T_{\tau}[\hat{A}(\tau)\hat{B}(0)]\Big{\rangle}=-\theta(\tau)\Big{\langle}e^{\tau\hat{H}}\hat{A}e^{-\tau\hat{H}}\hat{B}\Big{\rangle}\pm\theta(-\tau)\Big{\langle}\hat{B}e^{\tau\hat{H}}\hat{A}e^{-\tau\hat{H}}\Big{\rangle},$
(55)

where $+$ is for fermionic operators and $-$ is for bosonic operators, while
$\hat{A}$ and $\hat{B}$ stands for single-particle fermionic or bosonic
operators. According to this definition, one can simply prove the periodic and
anti-periodic properties of $G(\tau,\hat{A}\hat{B})$ for fermionic and bosonic
systems as

Fermion $\displaystyle:G(\tau,\hat{A}\hat{B})=-G(\beta+\tau,\hat{A}\hat{B})$
(56) Boson
$\displaystyle:G(\tau,\hat{A}\hat{B})=G(\beta+\tau,\hat{A}\hat{B}).$

We can also write down the Lehmann representation of $G(\tau,\hat{A}\hat{B})$
for both fermionic and bosonic systems via expanding the expectation by all
many-body eigenstates of the system as

$\displaystyle
G(\tau,\hat{A}\hat{B})=\frac{1}{Z}\sum_{mn}e^{(E_{m}-E_{n})\tau}\Big{[}-\theta(\tau)e^{-\beta
E_{m}}\pm\theta(-\tau)e^{-\beta E_{n}}\Big{]}\langle m|\hat{A}|n\rangle\langle
n|\hat{B}|m\rangle.$ (57)

For finite temperature case, the Fourier transformation between the imaginary-
time Green’s function $G(\tau,\hat{A}\hat{B})$ and the Matsubara frequency
Green’s function $G(i\omega_{n},\hat{A}\hat{B})$

$\displaystyle
G(i\omega_{n},\hat{A}\hat{B})=\int_{0}^{\beta}G(\tau,\hat{A}\hat{B})e^{i\omega_{n}\tau}d\tau.$
(58)

So the Lehmann representation expression of $G(i\omega_{n},\hat{A}\hat{B})$ is

$\displaystyle G(i\omega_{n},\hat{A}\hat{B})=\sum_{mn}D_{mn}\frac{\langle
n|\hat{A}|m\rangle\langle m|\hat{B}|n\rangle}{i\omega_{n}-(E_{m}-E_{n})}\hskip
34.14322ptD_{mn}=\frac{\pm e^{-\beta E_{m}}+e^{-\beta E_{n}}}{Z},$ (59)

where the sign $\pm$ originates from the term $e^{i\omega_{n}\beta}=+1$ for
bosons and $e^{i\omega_{n}\beta}=-1$ for fermions. The above formulae
summarize the basic properties of Matsubara frequency Green’s function.

At exact zero-temperature, we have the Lehmann representation for
$G(\tau,\mathbf{k})$ as

$\displaystyle
G(\tau,\hat{A}\hat{B})=\sum_{m}\Big{[}-\theta(\tau)e^{(E_{0}-E_{m})\tau}\langle
0|\hat{A}|m\rangle\langle
m|\hat{B}|0\rangle\pm\theta(-\tau)e^{-(E_{0}-E_{m})\tau}\langle
m|\hat{A}|0\rangle\langle 0|\hat{B}|m\rangle\Big{]}.$ (60)

Taking $\beta\to+\infty$ limit in Eq. (59), the Lehmann representation for
$G(i\omega,\hat{A}\hat{B})$ at zero-temperature can also be reached as

$\displaystyle G(i\omega,\hat{A}\hat{B})=\sum_{m}\Big{[}\frac{\langle
0|\hat{A}|m\rangle\langle
m|\hat{B}|0\rangle}{i\omega-(E_{m}-E_{0})}\pm\frac{\langle
m|\hat{A}|0\rangle\langle 0|\hat{B}|m\rangle}{i\omega+(E_{m}-E_{0})}\Big{]}.$
(61)

Now the problem is, at zero temperature, Fourier transformation in Eq. (58)
cannot transfer $G(\tau,\hat{A}\hat{B})$ in Eq. (60) to
$G(i\omega,\hat{A}\hat{B})$ in Eq. (61), since Eq. (61) cannot be obtained if
one naively takes $\beta\to+\infty$ limit in Eq. (58), which corresponds to
integration domain $[0,+\infty]$. To solve this problem, we present a simple
generalization of finite temperature Matsubara frequency Green’s function
formalism into the zero-temperature case, by altering the Fourier
transformation in Eq. (58) to the following one as

$\displaystyle
G(i\omega,\hat{A}\hat{B})=\int_{-\infty}^{+\infty}G(\tau,\hat{A}\hat{B})e^{i\omega\tau}d\tau.$
(62)

With this new transformation, we can now obtain Eq. (61) from Eq. (60), at
zero-temperature. In practical QMC simulations, we do not take the integral in
Eq. (62) to $\pm\infty$, and instead we carry out a cutoff for $\tau$ in the
integral as Eq. (15) in the main text. Nevertheless, Eq. (62) validates the
usage of Eq. (11) in the main text calculating $G(i\omega=0,\mathbf{k})$ from
$G(\tau,\mathbf{k})$ matrix obtained from QMC simulations.

##  Appendix C Proof of Eq. (12), Eq. (13) and Eq. (14)

The definition of $[G_{\sigma}(\tau,\mathbf{k})]_{pq}$ with $p,q=1,2,\cdots,m$
as orbitals for each spin sector in a unit cell, is

$\displaystyle[G_{\sigma}(\tau,\mathbf{k})]_{pq}=-\Big{\langle}T_{\tau}[c_{\mathbf{k}p{\sigma}}(\tau)c_{\mathbf{k}q{\sigma}}^{\dagger}(0)]\Big{\rangle}=-\frac{1}{N}\sum_{i,j=1}^{N}e^{-i\mathbf{k}\cdot(\mathbf{R}_{i}-\mathbf{R}_{j})}\Big{\langle}T_{\tau}[c_{ip{\sigma}}(\tau)c_{jq{\sigma}}^{\dagger}(0)]\Big{\rangle},$
(63)

where $N=L^{2}$ is the number of unit cells. Then we have

$\displaystyle[G_{\sigma}(\tau,\mathbf{k})]_{qp}^{\star}=-\frac{1}{N}\sum_{i,j=1}^{N}e^{i\mathbf{k}\cdot(\mathbf{R}_{i}-\mathbf{R}_{j})}\Big{\langle}T_{\tau}[c_{iq{\sigma}}(\tau)c_{jp{\sigma}}^{\dagger}(0)]\Big{\rangle}^{\star}.$
(64)

For any operator, we have
$\langle\phi|\hat{A}|\phi\rangle^{\star}=\langle\phi|\hat{A}^{\dagger}|\phi\rangle$,
which means $\langle\hat{A}\rangle^{\star}=\langle\hat{A}^{\dagger}\rangle$.
Based on this relation, we have

$\displaystyle\Big{\langle}T_{\tau}[c_{iq{\sigma}}(\tau)c_{jp{\sigma}}^{\dagger}(0)]\Big{\rangle}^{\star}=\Big{\langle}T_{\tau}[c_{jp{\sigma}}(\tau)c_{iq{\sigma}}^{\dagger}(0)]\Big{\rangle}.$
(65)

Then we can prove

$\displaystyle[G_{\sigma}(\tau,\mathbf{k})]_{qp}^{\star}=$
$\displaystyle-\frac{1}{N}\sum_{i,j=1}^{N}e^{i\mathbf{k}\cdot(\mathbf{R}_{i}-\mathbf{R}_{j})}\Big{\langle}T_{\tau}[c_{jp{\sigma}}(\tau)c_{iq{\sigma}}^{\dagger}(0)]\Big{\rangle}=-\frac{1}{N}\sum_{i,j=1}^{N}e^{-i\mathbf{k}\cdot(\mathbf{R}_{i}-\mathbf{R}_{j})}\Big{\langle}T_{\tau}[c_{ip{\sigma}}(\tau)c_{jq{\sigma}}^{\dagger}(0)]\Big{\rangle}$
$\displaystyle=$ $\displaystyle[G_{\sigma}(\tau,\mathbf{k})]_{pq}.\hskip
14.22636pt$ (66)

So the proof for Eq. (12) is complete.

For Eq. (13), the inversion operation $\hat{\mathcal{I}}$ (in each spin
sector) transforms $p$ sublattice to $p^{\prime}$. For GKMH model,
$\hat{\mathcal{I}}$ transfers A sublattice into B as $1\leftrightarrow 2$,
shown in Fig. 1 (a), while it transfers $1\leftrightarrow 4$,
$2\leftrightarrow 5$ and $3\leftrightarrow 6$ for CKMH model, shown in Fig. 5.
We assume the relation $p\leftrightarrow p^{\prime}$ and $q\leftrightarrow
q^{\prime}$ under the spatial inversion symmetry operation for generally
multi-band systems. As the position vector should be inverse under
$\hat{\mathcal{I}}$, i.e., $\mathbf{R}_{i}\to-\mathbf{R}_{i}$ we obtain the
transformation for simple operators
$c_{\mathbf{k}p\sigma},c_{\mathbf{k}q\sigma}$ in reciprocal space as

$\displaystyle
c_{\mathbf{k}p\sigma}=\frac{1}{\sqrt{N}}\sum_{i}e^{-i\mathbf{k}\cdot\mathbf{R}_{i}}c_{ip\sigma}\hskip
11.38092pt\to$
$\displaystyle\mathcal{I}c_{\mathbf{k}p\sigma}\mathcal{I}^{-1}=\frac{1}{\sqrt{N}}\sum_{i}\mathcal{I}e^{-i\mathbf{k}\cdot\mathbf{R}_{i}}\mathcal{I}^{-1}\cdot\mathcal{I}^{-1}c_{ip\sigma}\mathcal{I}^{-1}=\frac{1}{\sqrt{N}}\sum_{i}e^{i\mathbf{k}\cdot\mathbf{R}_{i}}c_{ip^{\prime}\sigma}=c_{-\mathbf{k}p^{\prime}\sigma}.$
(67)

Similarly, we have

$\displaystyle\mathcal{I}c_{\mathbf{k}p\sigma}^{\dagger}\mathcal{I}^{-1}=c_{-\mathbf{k}p^{\prime}\sigma}^{\dagger}$
$\displaystyle\mathcal{I}c_{\mathbf{k}q\sigma}(\tau)\mathcal{I}^{-1}=c_{-\mathbf{k}q^{\prime}\sigma}(\tau),\hskip
22.76228pt\mathcal{I}c_{\mathbf{k}q\sigma}^{\dagger}(\tau)\mathcal{I}^{-1}=c_{-\mathbf{k}q^{\prime}\sigma}^{\dagger}(\tau).$
(68)

So for the imaginary-time Green’s function matrix, we have

$\displaystyle\mathcal{I}[G_{\sigma}(\tau,\mathbf{k})]_{pq}\mathcal{I}^{-1}=$
$\displaystyle-\Big{\langle}T_{\tau}[\mathcal{I}c_{\mathbf{k}p{\sigma}}(\tau)c_{\mathbf{k}q{\sigma}}^{\dagger}\mathcal{I}^{-1}]\Big{\rangle}=-\Big{\langle}T_{\tau}[\mathcal{I}c_{\mathbf{k}p{\sigma}}(\tau)\mathcal{I}^{-1}\mathcal{I}c_{\mathbf{k}q{\sigma}}^{\dagger}\mathcal{I}^{-1}]\Big{\rangle}=-\Big{\langle}c_{-\mathbf{k}p^{\prime}{\sigma}}(\tau)c_{-\mathbf{k}q^{\prime}{\sigma}}^{\dagger}\Big{\rangle}$
(69) $\displaystyle=$
$\displaystyle[G_{\sigma}(\tau,-\mathbf{k})]_{p^{\prime}q^{\prime}}$

Due to the spatial inversion symmetry, the $G(\tau,\mathbf{k})$ matrix should
be invariant under the inversion symmetry operation, from which we can get

$\displaystyle[G_{\sigma}(\tau,\mathbf{k})]_{pq}=[G_{\sigma}(\tau,-\mathbf{k})]_{p^{\prime}q^{\prime}}.$
(70)

This is exactly the Eq. (13).

Finally, we prove the relation Eq. (14). The standard definition for
$[G_{\sigma}(\tau,\mathbf{k})]_{pq}$ is

$\displaystyle[G_{\sigma}(\tau,\mathbf{k})]_{pq}$ $\displaystyle=$
$\displaystyle-\Big{\langle}T_{\tau}[c_{\mathbf{k}p{\sigma}}(\tau)c_{\mathbf{k}q{\sigma}}^{\dagger}(0)]\Big{\rangle}=-\theta(\tau)\Big{\langle}c_{\mathbf{k}p{\sigma}}(\tau)c_{\mathbf{k}q{\sigma}}^{\dagger}\Big{\rangle}+\theta(-\tau)\Big{\langle}c_{\mathbf{k}q{\sigma}}^{\dagger}c_{\mathbf{k}p{\sigma}}(\tau)\Big{\rangle}$
$\displaystyle=$
$\displaystyle-\theta(\tau)\Big{\langle}c_{\mathbf{k}p{\sigma}}c_{\mathbf{k}q{\sigma}}^{\dagger}(-\tau)\Big{\rangle}+\theta(-\tau)\Big{\langle}c_{\mathbf{k}q{\sigma}}^{\dagger}(-\tau)c_{\mathbf{k}p{\sigma}}\Big{\rangle}$
$\displaystyle=$
$\displaystyle-\theta(\tau)\frac{1}{N}\sum_{i,j=1}^{N}e^{-i\mathbf{k}\cdot(\vec{R}_{i}-\vec{R}_{j})}\Big{\langle}c_{ip{\sigma}}c_{jq{\sigma}}^{\dagger}(-\tau)\Big{\rangle}+\theta(-\tau)\frac{1}{N}\sum_{i,j=1}^{N}e^{-i\mathbf{k}\cdot(\vec{R}_{i}-\vec{R}_{j})}\Big{\langle}c_{jq{\sigma}}^{\dagger}(-\tau)c_{ip{\sigma}}\Big{\rangle},$

where we have applied the the following relations

$\displaystyle\Big{\langle}c_{\mathbf{k}p{\sigma}}(\tau)c_{\mathbf{k}q{\sigma}}^{\dagger}\Big{\rangle}=\Big{\langle}c_{\mathbf{k}p{\sigma}}c_{\mathbf{k}q{\sigma}}^{\dagger}(-\tau)\Big{\rangle}\hskip
28.45274pt\Big{\langle}c_{\mathbf{k}q{\sigma}}^{\dagger}c_{\mathbf{k}p{\sigma}}(\tau)\Big{\rangle}=\Big{\langle}c_{\mathbf{k}q{\sigma}}^{\dagger}(-\tau)c_{\mathbf{k}p{\sigma}}\Big{\rangle}.$
(72)

Then we carry out the particle-hole transformation for
$[G_{\sigma}(\tau,\mathbf{k})]_{pq}$ as $c_{p}\to\xi_{p}d_{p}^{\dagger}$ and
$c_{q}\to\xi_{q}d_{q}^{\dagger}$, we have

$\displaystyle\frac{1}{N}\sum_{i,j=1}^{N}e^{-i\mathbf{k}\cdot(\mathbf{R}_{i}-\mathbf{R}_{j})}\Big{\langle}c_{ip{\sigma}}c_{jq{\sigma}}^{\dagger}(-\tau)\Big{\rangle}\hskip
8.5359pt\longrightarrow\hskip 8.5359pt$
$\displaystyle\frac{1}{N}\sum_{i,j=1}^{N}e^{-i\mathbf{k}\cdot(\mathbf{R}_{i}-\mathbf{R}_{j})}\xi_{p}\xi_{q}\Big{\langle}d_{ip{\sigma}}^{\dagger}d_{jq{\sigma}}(-\tau)\Big{\rangle}$
$\displaystyle=\xi_{p}\xi_{q}\Big{\langle}d_{-\mathbf{k}p{\sigma}}^{\dagger}d_{-\mathbf{k}q{\sigma}}(-\tau)\Big{\rangle},$

and the other term

$\displaystyle\frac{1}{N}\sum_{i,j=1}^{N}e^{-i\mathbf{k}\cdot(\mathbf{R}_{i}-\mathbf{R}_{j})}\Big{\langle}c_{jq{\sigma}}^{\dagger}(-\tau)c_{ip{\sigma}}\Big{\rangle}\hskip
8.5359pt\longrightarrow\hskip 8.5359pt$
$\displaystyle\frac{1}{N}\sum_{i,j=1}^{N}e^{-i\mathbf{k}\cdot(\mathbf{R}_{i}-\mathbf{R}_{j})}\xi_{p}\xi_{q}\Big{\langle}d_{jq{\sigma}}(-\tau)d_{ip{\sigma}}^{\dagger}\Big{\rangle}$
$\displaystyle=\xi_{p}\xi_{q}\Big{\langle}d_{-\mathbf{k}q{\sigma}}(-\tau)d_{-\mathbf{k}p{\sigma}}^{\dagger}\Big{\rangle}.$

Combining the results of these two terms, under the particle-hole
transformation, $[G_{\sigma}(\tau,\mathbf{k})]_{pq}$ changes to

$\displaystyle[G_{\sigma}(\tau,\mathbf{k})]_{pq}\hskip
8.5359pt\longrightarrow\hskip
8.5359pt-\theta(\tau)\xi_{p}\xi_{q}\Big{\langle}d_{-\mathbf{k}p{\sigma}}^{\dagger}d_{-\mathbf{k}q{\sigma}}(-\tau)\Big{\rangle}+\theta(-\tau)\xi_{p}\xi_{q}\Big{\langle}d_{-\mathbf{k}q{\sigma}}(-\tau)d_{-\mathbf{k}p{\sigma}}^{\dagger}\Big{\rangle}.$
(75)

On the other hand, we can write down $[G_{\sigma}(-\tau,-\mathbf{k})]_{qp}$
according to original definition in Eq. (C)

$\displaystyle[G_{\sigma}(-\tau,-\mathbf{k})]_{qp}=-\theta(-\tau)\Big{\langle}c_{-\mathbf{k}q{\sigma}}(-\tau)c_{-\mathbf{k}p{\sigma}}^{\dagger}\Big{\rangle}+\theta(\tau)\Big{\langle}c_{-\mathbf{k}p{\sigma}}^{\dagger}c_{-\mathbf{k}q{\sigma}}(-\tau)\Big{\rangle}$
(76)

Comparing Eq. (75) and Eq. (76), and considering that the particle-hole
symmetry is preserved, we arrive at

$\displaystyle[G_{\sigma}(\tau,\mathbf{k})]_{pq}=-\xi_{p}\xi_{q}[G_{\sigma}(-\tau,-\mathbf{k})]_{qp},$
(77)

which is Eq. (14) in the main text.

##  Appendix D Derivation of Eq. (28) and Eq. (IV.2)

For the CKMH model, the spatial inversion operators as $1\leftrightarrow 4$,
$2\leftrightarrow 5$ and $3\leftrightarrow 6$ among the six sublattices. One
can show that spatial inversion symmetry result in the following properties of
$G_{\sigma}(\tau,\mathbf{k})$ matrix, according to Eq. (13).

$\displaystyle[G_{\sigma}(\tau,\mathbf{k})]_{11}=[G_{\sigma}(\tau,-\mathbf{k})]_{44}\hskip
28.45274pt[G_{\sigma}(\tau,\mathbf{k})]_{12}=[G_{\sigma}(\tau,-\mathbf{k})]_{45}\hskip
28.45274pt[G_{\sigma}(\tau,\mathbf{k})]_{13}=[G_{\sigma}(\tau,-\mathbf{k})]_{46}$
$\displaystyle[G_{\sigma}(\tau,\mathbf{k})]_{14}=[G_{\sigma}(\tau,-\mathbf{k})]_{41}\hskip
28.45274pt[G_{\sigma}(\tau,\mathbf{k})]_{15}=[G_{\sigma}(\tau,-\mathbf{k})]_{42}\hskip
28.45274pt[G_{\sigma}(\tau,\mathbf{k})]_{16}=[G_{\sigma}(\tau,-\mathbf{k})]_{43}$
$\displaystyle[G_{\sigma}(\tau,\mathbf{k})]_{22}=[G_{\sigma}(\tau,-\mathbf{k})]_{55}\hskip
28.45274pt[G_{\sigma}(\tau,\mathbf{k})]_{23}=[G_{\sigma}(\tau,-\mathbf{k})]_{56}\hskip
28.45274pt[G_{\sigma}(\tau,\mathbf{k})]_{24}=[G_{\sigma}(\tau,-\mathbf{k})]_{51}$
$\displaystyle[G_{\sigma}(\tau,\mathbf{k})]_{25}=[G_{\sigma}(\tau,-\mathbf{k})]_{52}\hskip
28.45274pt[G_{\sigma}(\tau,\mathbf{k})]_{26}=[G_{\sigma}(\tau,-\mathbf{k})]_{53}$
$\displaystyle[G_{\sigma}(\tau,\mathbf{k})]_{33}=[G_{\sigma}(\tau,-\mathbf{k})]_{66}\hskip
28.45274pt[G_{\sigma}(\tau,\mathbf{k})]_{34}=[G_{\sigma}(\tau,-\mathbf{k})]_{61}\hskip
28.45274pt[G_{\sigma}(\tau,\mathbf{k})]_{35}=[G_{\sigma}(\tau,-\mathbf{k})]_{62}$
$\displaystyle[G_{\sigma}(\tau,\mathbf{k})]_{36}=[G_{\sigma}(\tau,-\mathbf{k})]_{63}$
$\displaystyle[G_{\sigma}(\tau,\mathbf{k})]_{44}=[G_{\sigma}(\tau,-\mathbf{k})]_{11}\hskip
28.45274pt[G_{\sigma}(\tau,\mathbf{k})]_{45}=[G_{\sigma}(\tau,-\mathbf{k})]_{12}\hskip
28.45274pt[G_{\sigma}(\tau,\mathbf{k})]_{46}=[G_{\sigma}(\tau,-\mathbf{k})]_{13}$
$\displaystyle[G_{\sigma}(\tau,\mathbf{k})]_{55}=[G_{\sigma}(\tau,-\mathbf{k})]_{22}\hskip
28.45274pt[G_{\sigma}(\tau,\mathbf{k})]_{56}=[G_{\sigma}(\tau,-\mathbf{k})]_{23}$
$\displaystyle[G_{\sigma}(\tau,\mathbf{k})]_{66}=[G_{\sigma}(\tau,-\mathbf{k})]_{33}$
(78)

The particle-hole operators with
$\xi_{1}=+1,\xi_{2}=-1,\xi_{3}=+1,\xi_{4}=-1,\xi_{5}=+1,\xi_{6}=-1$, then
according to Eq. (14), we have the following properties for
$[G_{\sigma}(\tau,\mathbf{k})]_{1q},q=1,2,3,4,5,6$ as

$\displaystyle[G_{\sigma}(\tau,-\mathbf{k})]_{44}=-[G_{\sigma}(-\tau,\mathbf{k})]_{44}\hskip
28.45274pt[G_{\sigma}(\tau,-\mathbf{k})]_{45}=+[G_{\sigma}(-\tau,\mathbf{k})]_{54}\hskip
28.45274pt[G_{\sigma}(\tau,-\mathbf{k})]_{46}=-[G_{\sigma}(-\tau,\mathbf{k})]_{64}\hskip
14.22636pt$
$\displaystyle[G_{\sigma}(\tau,-\mathbf{k})]_{41}=+[G_{\sigma}(-\tau,\mathbf{k})]_{14}\hskip
28.45274pt[G_{\sigma}(\tau,-\mathbf{k})]_{42}=-[G_{\sigma}(-\tau,\mathbf{k})]_{24}\hskip
28.45274pt[G_{\sigma}(\tau,-\mathbf{k})]_{43}=+[G_{\sigma}(-\tau,\mathbf{k})]_{34}.$
(79)

With the condition of Hermitian matrix of $G_{\sigma}(\tau,\mathbf{k})$ in Eq.
(12), we have
$[G_{\sigma}(-\tau,\mathbf{k})]_{54}=[G_{\sigma}(-\tau,\mathbf{k})]_{45}^{\star},[G_{\sigma}(-\tau,\mathbf{k})]_{64}=[G_{\sigma}(-\tau,\mathbf{k})]_{46}^{\star}$.
With Eq. (D) and Eq. (D), we have

$\displaystyle[G_{\sigma}(\tau,\mathbf{k})]_{11}=-[G_{\sigma}(-\tau,\mathbf{k})]_{44}\hskip
28.45274pt[G_{\sigma}(\tau,\mathbf{k})]_{12}=+[G_{\sigma}(-\tau,\mathbf{k})]_{45}^{\star}\hskip
28.45274pt[G_{\sigma}(\tau,\mathbf{k})]_{13}=-[G_{\sigma}(-\tau,\mathbf{k})]_{46}^{\star}$
$\displaystyle[G_{\sigma}(\tau,\mathbf{k})]_{14}=+[G_{\sigma}(-\tau,\mathbf{k})]_{14}\hskip
28.45274pt[G_{\sigma}(\tau,\mathbf{k})]_{15}=-[G_{\sigma}(-\tau,\mathbf{k})]_{24}\hskip
28.45274pt[G_{\sigma}(\tau,\mathbf{k})]_{16}=+[G_{\sigma}(-\tau,\mathbf{k})]_{34}.$
(80)

Similarly, for $[G_{\sigma}(\tau,\mathbf{k})]_{2q},q=2,3,4,5,6$ from the
particle-hole symmetry, we have

$\displaystyle[G_{\sigma}(\tau,-\mathbf{k})]_{55}=-[G_{\sigma}(-\tau,\mathbf{k})]_{55}\hskip
28.45274pt[G_{\sigma}(\tau,-\mathbf{k})]_{56}=+[G_{\sigma}(-\tau,\mathbf{k})]_{65}\hskip
28.45274pt[G_{\sigma}(\tau,-\mathbf{k})]_{51}=-[G_{\sigma}(-\tau,\mathbf{k})]_{15}\hskip
14.22636pt$ (81)
$\displaystyle[G_{\sigma}(\tau,-\mathbf{k})]_{52}=+[G_{\sigma}(-\tau,\mathbf{k})]_{25}\hskip
28.45274pt[G_{\sigma}(\tau,-\mathbf{k})]_{53}=-[G_{\sigma}(-\tau,\mathbf{k})]_{35}.$

With
$[G_{\sigma}(-\tau,\mathbf{k})]_{65}=[G_{\sigma}(-\tau,\mathbf{k})]_{56}^{\star}$
and combining Eq. (D) and Eq. (81), we also have

$\displaystyle[G_{\sigma}(\tau,\mathbf{k})]_{22}=-[G_{\sigma}(-\tau,\mathbf{k})]_{55}\hskip
28.45274pt[G_{\sigma}(\tau,\mathbf{k})]_{23}=+[G_{\sigma}(-\tau,\mathbf{k})]_{56}^{\star}\hskip
28.45274pt[G_{\sigma}(\tau,\mathbf{k})]_{24}=-[G_{\sigma}(-\tau,\mathbf{k})]_{15}$
$\displaystyle[G_{\sigma}(\tau,\mathbf{k})]_{25}=+[G_{\sigma}(-\tau,\mathbf{k})]_{25}\hskip
28.45274pt[G_{\sigma}(\tau,\mathbf{k})]_{26}=-[G_{\sigma}(-\tau,\mathbf{k})]_{35}.$
(82)

For $[G_{\sigma}(\tau,\mathbf{k})]_{3q},q=3,4,5,6$ from the particle-hole
symmetry, we have

$\displaystyle[G_{\sigma}(\tau,-\mathbf{k})]_{66}=-[G_{\sigma}(-\tau,\mathbf{k})]_{66}\hskip
28.45274pt[G_{\sigma}(\tau,-\mathbf{k})]_{61}=+[G_{\sigma}(-\tau,\mathbf{k})]_{16}\hskip
28.45274pt[G_{\sigma}(\tau,-\mathbf{k})]_{62}=-[G_{\sigma}(-\tau,\mathbf{k})]_{26}\hskip
14.22636pt$
$\displaystyle[G_{\sigma}(\tau,-\mathbf{k})]_{63}=+[G_{\sigma}(-\tau,\mathbf{k})]_{36}.$
(83)

Combining Eq. (D) and Eq. (D), we also have

$\displaystyle[G_{\sigma}(\tau,\mathbf{k})]_{33}=-[G_{\sigma}(-\tau,\mathbf{k})]_{66}\hskip
28.45274pt[G_{\sigma}(\tau,\mathbf{k})]_{34}=+[G_{\sigma}(-\tau,\mathbf{k})]_{16}\hskip
28.45274pt[G_{\sigma}(\tau,\mathbf{k})]_{35}=-[G_{\sigma}(-\tau,\mathbf{k})]_{26}\hskip
14.22636pt$
$\displaystyle[G_{\sigma}(\tau,\mathbf{k})]_{36}=+[G_{\sigma}(-\tau,\mathbf{k})]_{36}.$
(84)

For $[G_{\sigma}(\tau,\mathbf{k})]_{4q},q=4,5,6$ from the particle-hole
symmetry, we have

$\displaystyle[G_{\sigma}(\tau,-\mathbf{k})]_{11}=-[G_{\sigma}(-\tau,\mathbf{k})]_{11}\hskip
28.45274pt[G_{\sigma}(\tau,-\mathbf{k})]_{12}=+[G_{\sigma}(-\tau,\mathbf{k})]_{21}\hskip
28.45274pt[G_{\sigma}(\tau,-\mathbf{k})]_{13}=-[G_{\sigma}(-\tau,\mathbf{k})]_{31}\hskip
14.22636pt$ (85)

With
$[G_{\sigma}(-\tau,\mathbf{k})]_{21}=[G_{\sigma}(-\tau,\mathbf{k})]_{12}^{\star},[G_{\sigma}(-\tau,\mathbf{k})]_{31}=[G_{\sigma}(-\tau,\mathbf{k})]_{13}^{\star}$
and combining Eq. (D) and Eq. (85), we also have

$\displaystyle[G_{\sigma}(\tau,\mathbf{k})]_{44}=-[G_{\sigma}(-\tau,\mathbf{k})]_{11}\hskip
28.45274pt[G_{\sigma}(\tau,\mathbf{k})]_{45}=+[G_{\sigma}(-\tau,\mathbf{k})]_{12}^{\star}\hskip
28.45274pt[G_{\sigma}(\tau,\mathbf{k})]_{46}=-[G_{\sigma}(-\tau,\mathbf{k})]_{13}^{\star}$
(86)

For $[G_{\sigma}(\tau,\mathbf{k})]_{5q},q=5,6$ from the particle-hole
symmetry, we have

$\displaystyle[G_{\sigma}(\tau,-\mathbf{k})]_{22}=-[G_{\sigma}(-\tau,\mathbf{k})]_{22}\hskip
28.45274pt[G_{\sigma}(\tau,-\mathbf{k})]_{23}=+[G_{\sigma}(-\tau,\mathbf{k})]_{32}$
(87)

With
$[G_{\sigma}(-\tau,\mathbf{k})]_{32}=[G_{\sigma}(-\tau,\mathbf{k})]_{23}^{\star}$
and combining Eq. (D) and Eq. (87), we also have

$\displaystyle[G_{\sigma}(\tau,\mathbf{k})]_{55}=-[G_{\sigma}(-\tau,\mathbf{k})]_{22}\hskip
28.45274pt[G_{\sigma}(\tau,\mathbf{k})]_{56}=+[G_{\sigma}(-\tau,\mathbf{k})]_{23}^{\star}$
(88)

For $[G_{\sigma}(\tau,\mathbf{k})]_{66}$ from the particle-hole symmetry, we
have

$\displaystyle[G_{\sigma}(\tau,-\mathbf{k})]_{33}=-[G_{\sigma}(-\tau,\mathbf{k})]_{33}$
(89)

Combining Eq. (D) and Eq. (89), we also have

$\displaystyle[G_{\sigma}(\tau,\mathbf{k})]_{66}=-[G_{\sigma}(-\tau,\mathbf{k})]_{33}$
(90)

To calculate the $G_{\sigma}(i\omega=0,\mathbf{k})$ matrix, we need to apply
the Fourier transformation as

$\displaystyle
G_{\sigma}(i\omega=0,\mathbf{k})=\int_{-\infty}^{+\infty}G_{\sigma}(\tau,\mathbf{k})d\tau.$
(91)

Then by Eq. (91), we can show that there are actually only 12 independent
matrix elements in $G_{\sigma}(i\omega=0,\mathbf{k})$ matrix for CKMH model,
which we need to calculate, combining Eq. (D), Eq. (D), Eq. (D), Eq. (86), Eq.
(88) and Eq. (90). For the diagonal matrix elements of
$G_{\sigma}(i\omega=0,\mathbf{k})$ matrix, we have 3 independent diagonal
elements as

$\displaystyle
A_{1}=[G_{\sigma}(i\omega=0,\mathbf{k})]_{11}=\int_{-\infty}^{+\infty}[G_{\sigma}(\tau,\mathbf{k})]_{11}d\tau=\int_{0}^{+\infty}\Big{\\{}[G_{\sigma}(\tau,\mathbf{k})]_{11}-[G_{\sigma}(\tau,\mathbf{k})]_{44}\Big{\\}}d\tau=-[G_{\sigma}(i\omega=0,\mathbf{k})]_{44}\hskip
17.07182pt$ $\displaystyle
A_{2}=[G_{\sigma}(i\omega=0,\mathbf{k})]_{22}=\int_{-\infty}^{+\infty}[G_{\sigma}(\tau,\mathbf{k})]_{22}d\tau=\int_{0}^{+\infty}\Big{\\{}[G_{\sigma}(\tau,\mathbf{k})]_{22}-[G_{\sigma}(\tau,\mathbf{k})]_{55}\Big{\\}}d\tau=-[G_{\sigma}(i\omega=0,\mathbf{k})]_{55}$
$\displaystyle
A_{3}=[G_{\sigma}(i\omega=0,\mathbf{k})]_{33}=\int_{-\infty}^{+\infty}[G_{\sigma}(\tau,\mathbf{k})]_{33}d\tau=\int_{0}^{+\infty}\Big{\\{}[G_{\sigma}(\tau,\mathbf{k})]_{33}-[G_{\sigma}(\tau,\mathbf{k})]_{66}\Big{\\}}d\tau=-[G_{\sigma}(i\omega=0,\mathbf{k})]_{66}$
(92)

For the off-diagonal matrix elements, we can determine that there are only 9
independent matrix off-diagonal elements in $G_{\sigma}(i\omega=0,\mathbf{k})$
matrix as

$\displaystyle
A_{4}=[G_{\sigma}(i\omega=0,\mathbf{k})]_{12}=\int_{-\infty}^{+\infty}[G_{\sigma}(\tau,\mathbf{k})]_{12}d\tau=\int_{0}^{+\infty}\Big{\\{}[G_{\sigma}(\tau,\mathbf{k})]_{12}+[G_{\sigma}(\tau,\mathbf{k})]_{45}^{\star}\Big{\\}}d\tau=+[G_{\sigma}(i\omega=0,\mathbf{k})]_{45}^{\star}\hskip
14.22636pt$ $\displaystyle
A_{5}=[G_{\sigma}(i\omega=0,\mathbf{k})]_{13}=\int_{-\infty}^{+\infty}[G_{\sigma}(\tau,\mathbf{k})]_{13}d\tau=\int_{0}^{+\infty}\Big{\\{}[G_{\sigma}(\tau,\mathbf{k})]_{13}-[G_{\sigma}(\tau,\mathbf{k})]_{46}^{\star}\Big{\\}}d\tau=-[G_{\sigma}(i\omega=0,\mathbf{k})]_{46}^{\star}$
$\displaystyle
A_{6}=[G_{\sigma}(i\omega=0,\mathbf{k})]_{14}=\int_{-\infty}^{+\infty}[G_{\sigma}(\tau,\mathbf{k})]_{14}d\tau=2\int_{0}^{+\infty}[G_{\sigma}(\tau,\mathbf{k})]_{14}d\tau$
$\displaystyle
A_{7}=[G_{\sigma}(i\omega=0,\mathbf{k})]_{15}=\int_{-\infty}^{+\infty}[G_{\sigma}(\tau,\mathbf{k})]_{15}d\tau=\int_{0}^{+\infty}\Big{\\{}[G_{\sigma}(\tau,\mathbf{k})]_{15}-[G_{\sigma}(\tau,\mathbf{k})]_{24}\Big{\\}}d\tau=-[G_{\sigma}(i\omega=0,\mathbf{k})]_{24}$
$\displaystyle
A_{8}=[G_{\sigma}(i\omega=0,\mathbf{k})]_{16}=\int_{-\infty}^{+\infty}[G_{\sigma}(\tau,\mathbf{k})]_{16}d\tau=\int_{0}^{+\infty}\Big{\\{}[G_{\sigma}(\tau,\mathbf{k})]_{16}+[G_{\sigma}(\tau,\mathbf{k})]_{34}\Big{\\}}d\tau=+[G_{\sigma}(i\omega=0,\mathbf{k})]_{34}$
$\displaystyle
A_{9}=[G_{\sigma}(i\omega=0,\mathbf{k})]_{23}=\int_{-\infty}^{+\infty}[G_{\sigma}(\tau,\mathbf{k})]_{23}d\tau=\int_{0}^{+\infty}\Big{\\{}[G_{\sigma}(\tau,\mathbf{k})]_{23}+[G_{\sigma}(\tau,\mathbf{k})]_{56}^{\star}\Big{\\}}d\tau=+[G_{\sigma}(i\omega=0,\mathbf{k})]_{56}^{\star}$
$\displaystyle
A_{10}=[G_{\sigma}(i\omega=0,\mathbf{k})]_{25}=\int_{-\infty}^{+\infty}[G_{\sigma}(\tau,\mathbf{k})]_{25}d\tau=2\int_{0}^{+\infty}[G_{\sigma}(\tau,\mathbf{k})]_{25}d\tau$
$\displaystyle
A_{11}=[G_{\sigma}(i\omega=0,\mathbf{k})]_{26}=\int_{-\infty}^{+\infty}[G_{\sigma}(\tau,\mathbf{k})]_{26}d\tau=\int_{0}^{+\infty}\Big{\\{}[G_{\sigma}(\tau,\mathbf{k})]_{26}-[G_{\sigma}(\tau,\mathbf{k})]_{35}\Big{\\}}d\tau=-[G_{\sigma}(i\omega=0,\mathbf{k})]_{35}$
$\displaystyle
A_{12}=[G_{\sigma}(i\omega=0,\mathbf{k})]_{36}=\int_{-\infty}^{+\infty}[G_{\sigma}(\tau,\mathbf{k})]_{36}d\tau=2\int_{0}^{+\infty}[G_{\sigma}(\tau,\mathbf{k})]_{36}d\tau$
(93)

Then taking the $\tau$ cutoff $\theta$ and transforming the integral to
summation over discrete $\tau$, we can get the following matrix structure of
$G_{\sigma}(i\omega=0,\mathbf{k})$ matrix as

$\displaystyle
G_{\sigma}(i\omega=0,\mathbf{k})=\left(\begin{array}[]{cccccc}A_{1}&A_{4}&A_{5}&A_{6}&A_{7}&A_{8}\\\
A_{4}^{\star}&A_{2}&A_{9}&-A_{7}&A_{10}&A_{11}\\\
A_{5}^{\star}&A_{9}^{\star}&A_{3}&A_{8}&-A_{11}&A_{12}\\\
A_{6}^{\star}&-A_{7}^{\star}&A_{8}^{\star}&-A_{1}&A_{4}^{\star}&-A_{5}^{\star}\\\
A_{7}^{\star}&A_{10}^{\star}&-A_{11}^{\star}&A_{4}&-A_{2}&A_{9}^{\star}\\\
A_{8}^{\star}&A_{11}^{\star}&A_{12}^{\star}&-A_{5}&A_{9}&-A_{3}\end{array}\right).$
(100)

As for the $Z_{2}$ invariant for CKMH model, we only need to obtain the
$G_{\sigma}(i\omega=0,\bm{\kappa})$ data at four TRIM points as
$\bm{\kappa}=\bm{\Gamma},\mathbf{M}_{1},\mathbf{M}_{2},\mathbf{M}_{3}$. From
the symmetry properties in Eq. (D), we can obtain that for the TRIM points, we
have

$\displaystyle[G_{\sigma}(\tau,\bm{\kappa})]_{11}=[G_{\sigma}(\tau,\bm{\kappa})]_{44}\hskip
28.45274pt[G_{\sigma}(\tau,\bm{\kappa})]_{12}=[G_{\sigma}(\tau,\bm{\kappa})]_{45}\hskip
28.45274pt[G_{\sigma}(\tau,\bm{\kappa})]_{13}=[G_{\sigma}(\tau,\bm{\kappa})]_{46}$
$\displaystyle[G_{\sigma}(\tau,\bm{\kappa})]_{14}=[G_{\sigma}(\tau,\bm{\kappa})]_{41}\hskip
28.45274pt[G_{\sigma}(\tau,\bm{\kappa})]_{15}=[G_{\sigma}(\tau,\bm{\kappa})]_{42}\hskip
28.45274pt[G_{\sigma}(\tau,\bm{\kappa})]_{16}=[G_{\sigma}(\tau,\bm{\kappa})]_{43}$
$\displaystyle[G_{\sigma}(\tau,\bm{\kappa})]_{22}=[G_{\sigma}(\tau,\bm{\kappa})]_{55}\hskip
28.45274pt[G_{\sigma}(\tau,\bm{\kappa})]_{23}=[G_{\sigma}(\tau,\bm{\kappa})]_{56}\hskip
28.45274pt[G_{\sigma}(\tau,\bm{\kappa})]_{24}=[G_{\sigma}(\tau,\bm{\kappa})]_{51}$
$\displaystyle[G_{\sigma}(\tau,\bm{\kappa})]_{25}=[G_{\sigma}(\tau,\bm{\kappa})]_{52}\hskip
28.45274pt[G_{\sigma}(\tau,\bm{\kappa})]_{26}=[G_{\sigma}(\tau,\bm{\kappa})]_{53}$
$\displaystyle[G_{\sigma}(\tau,\bm{\kappa})]_{33}=[G_{\sigma}(\tau,\bm{\kappa})]_{66}\hskip
28.45274pt[G_{\sigma}(\tau,\bm{\kappa})]_{34}=[G_{\sigma}(\tau,\bm{\kappa})]_{61}\hskip
28.45274pt[G_{\sigma}(\tau,\bm{\kappa})]_{35}=[G_{\sigma}(\tau,\bm{\kappa})]_{62}$
$\displaystyle[G_{\sigma}(\tau,\bm{\kappa})]_{36}=[G_{\sigma}(\tau,\bm{\kappa})]_{63}$
$\displaystyle[G_{\sigma}(\tau,\bm{\kappa})]_{44}=[G_{\sigma}(\tau,\bm{\kappa})]_{11}\hskip
28.45274pt[G_{\sigma}(\tau,\bm{\kappa})]_{45}=[G_{\sigma}(\tau,\bm{\kappa})]_{12}\hskip
28.45274pt[G_{\sigma}(\tau,\bm{\kappa})]_{46}=[G_{\sigma}(\tau,\bm{\kappa})]_{13}$
$\displaystyle[G_{\sigma}(\tau,\bm{\kappa})]_{55}=[G_{\sigma}(\tau,\bm{\kappa})]_{22}\hskip
28.45274pt[G_{\sigma}(\tau,\bm{\kappa})]_{56}=[G_{\sigma}(\tau,\bm{\kappa})]_{23}$
$\displaystyle[G_{\sigma}(\tau,\bm{\kappa})]_{66}=[G_{\sigma}(\tau,\bm{\kappa})]_{33}$
(101)

Then we can determine that at TRIM points, we have $A_{1}=A_{2}=A_{3}=0$ and
$[G_{\sigma}(\tau,\bm{\kappa})]_{14}$, $[G_{\sigma}(\tau,\bm{\kappa})]_{25}$
and $[G_{\sigma}(\tau,\bm{\kappa})]_{36}$ are all real numbers. From these
relations, we have $B_{6}=A_{6}$, $B_{10}=A_{10}$ and $B_{12}=A_{12}$, and
they are all real numbers as well. For other elements, we have

$\displaystyle
A_{4}=[G_{\sigma}(i\omega=0,\bm{\kappa})]_{12}=\int_{-\infty}^{+\infty}[G_{\sigma}(\tau,\bm{\kappa})]_{12}d\tau=\int_{0}^{+\infty}\Big{\\{}[G_{\sigma}(\tau,\bm{\kappa})]_{12}+[G_{\sigma}(\tau,\bm{\kappa})]_{12}^{\star}\Big{\\}}d\tau=B_{4}$
$\displaystyle
A_{5}=[G_{\sigma}(i\omega=0,\bm{\kappa})]_{13}=\int_{-\infty}^{+\infty}[G_{\sigma}(\tau,\bm{\kappa})]_{13}d\tau=\int_{0}^{+\infty}\Big{\\{}[G_{\sigma}(\tau,\bm{\kappa})]_{13}-[G_{\sigma}(\tau,\bm{\kappa})]_{13}^{\star}\Big{\\}}d\tau=iB_{5}$
$\displaystyle
A_{6}=[G_{\sigma}(i\omega=0,\bm{\kappa})]_{14}=\int_{-\infty}^{+\infty}[G_{\sigma}(\tau,\bm{\kappa})]_{14}d\tau=2\int_{0}^{+\infty}[G_{\sigma}(\tau,\bm{\kappa})]_{14}d\tau=B_{6}$
$\displaystyle
A_{7}=[G_{\sigma}(i\omega=0,\bm{\kappa})]_{15}=\int_{-\infty}^{+\infty}[G_{\sigma}(\tau,\bm{\kappa})]_{15}d\tau=\int_{0}^{+\infty}\Big{\\{}[G_{\sigma}(\tau,\bm{\kappa})]_{15}-[G_{\sigma}(\tau,\bm{\kappa})]_{15}^{\star}\Big{\\}}d\tau=iB_{7}$
$\displaystyle
A_{8}=[G_{\sigma}(i\omega=0,\bm{\kappa})]_{16}=\int_{-\infty}^{+\infty}[G_{\sigma}(\tau,\bm{\kappa})]_{16}d\tau=\int_{0}^{+\infty}\Big{\\{}[G_{\sigma}(\tau,\bm{\kappa})]_{16}+[G_{\sigma}(\tau,\bm{\kappa})]_{16}^{\star}\Big{\\}}d\tau=B_{8}$
$\displaystyle
A_{9}=[G_{\sigma}(i\omega=0,\bm{\kappa})]_{23}=\int_{-\infty}^{+\infty}[G_{\sigma}(\tau,\bm{\kappa})]_{23}d\tau=\int_{0}^{+\infty}\Big{\\{}[G_{\sigma}(\tau,\bm{\kappa})]_{23}+[G_{\sigma}(\tau,\bm{\kappa})]_{23}^{\star}\Big{\\}}d\tau=B_{9}$
$\displaystyle
A_{10}=[G_{\sigma}(i\omega=0,\bm{\kappa})]_{25}=\int_{-\infty}^{+\infty}[G_{\sigma}(\tau,\bm{\kappa})]_{25}d\tau=2\int_{0}^{+\infty}[G_{\sigma}(\tau,\bm{\kappa})]_{25}d\tau=B_{10}$
$\displaystyle
A_{11}=[G_{\sigma}(i\omega=0,\bm{\kappa})]_{26}=\int_{-\infty}^{+\infty}[G_{\sigma}(\tau,\bm{\kappa})]_{26}d\tau=\int_{0}^{+\infty}\Big{\\{}[G_{\sigma}(\tau,\bm{\kappa})]_{26}-[G_{\sigma}(\tau,\bm{\kappa})]_{26}^{\star}\Big{\\}}d\tau=iB_{11}$
$\displaystyle
A_{12}=[G_{\sigma}(i\omega=0,\mathbf{k})]_{36}=\int_{-\infty}^{+\infty}[G_{\sigma}(\tau,\mathbf{k})]_{36}d\tau=2\int_{0}^{+\infty}[G_{\sigma}(\tau,\mathbf{k})]_{36}d\tau=B_{12}$
(102)

We can observe that all the matrix elements of $G_{\sigma}(\tau,\bm{\kappa})$
matrix must be either purely real or purely imaginary, and the diagonal matrix
elements are all zero. Hence, the matrix structure as
$G_{\sigma}(\tau,\bm{\kappa})$ is as following,

$\displaystyle
G_{\sigma}(i\omega=0,\bm{\kappa})=\left(\begin{array}[]{cccccc}0&B_{4}&iB_{5}&B_{6}&iB_{7}&B_{8}\\\
B_{4}&0&B_{9}&-iB_{7}&B_{10}&iB_{11}\\\
-iB_{5}&B_{9}&0&B_{8}&-iB_{11}&B_{12}\\\ B_{6}&iB_{7}&B_{8}&0&B_{4}&iB_{5}\\\
-iB_{7}&B_{10}&iB_{11}&B_{4}&0&B_{9}\\\
B_{8}&-iB_{11}&B_{12}&-iB_{5}&B_{9}&0\end{array}\right),$ (109)

in which $\\{B_{i}|i=4,5,\cdots,12\\}$ are purely real numbers.

## References

  * Hasan and Kane (2010) M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * Qi and Zhang (2011) X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * Klitzing  _et al._ (1980) K. v. Klitzing, G. Dorda, and M. Pepper, Phys. Rev. Lett. 45, 494 (1980).
  * Haldane (1988) F. D. M. Haldane, Phys. Rev. Lett. 61, 2015 (1988).
  * Thouless  _et al._ (1982) D. J. Thouless, M. Kohmoto, M. P. Nightingale, and M. den Nijs, Phys. Rev. Lett. 49, 405 (1982).
  * Avron  _et al._ (1983) J. E. Avron, R. Seiler, and B. Simon, Phys. Rev. Lett. 51, 51 (1983).
  * Kane and Mele (2005a) C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 226801 (2005a).
  * Kane and Mele (2005b) C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 146802 (2005b).
  * Fu and Kane (2006) L. Fu and C. L. Kane, Phys. Rev. B 74, 195312 (2006).
  * Sheng  _et al._ (2006) D. N. Sheng, Z. Y. Weng, L. Sheng, and F. D. M. Haldane, Phys. Rev. Lett. 97, 036808 (2006).
  * Prodan (2009) E. Prodan, Phys. Rev. B 80, 125327 (2009).
  * Fu and Kane (2007) L. Fu and C. L. Kane, Phys. Rev. B 76, 045302 (2007).
  * So (1985) H. So, Progress of Theoretical Physics 74, 585 (1985).
  * Ishikawa and Matsuyama (1986) K. Ishikawa and T. Matsuyama, Z. Phys. C 33, 41 (1986).
  * Volovik (1988) G. E. Volovik, JETP 67, 1804 (1988).
  * Wang  _et al._ (2010) Z. Wang, X.-L. Qi, and S.-C. Zhang, Phys. Rev. Lett. 105, 256803 (2010).
  * Gurarie (2011) V. Gurarie, Phys. Rev. B 83, 085426 (2011).
  * Volovik (2009) G. E. Volovik,  _The universe in a helium droplet_ (Oxford University Press New York, 2009).
  * Niu  _et al._ (1985) Q. Niu, D. J. Thouless, and Y.-S. Wu, Phys. Rev. B 31, 3372 (1985).
  * Wang and Zhang (2014) Z. Wang and S.-C. Zhang, Phys. Rev. X 4, 011006 (2014).
  * Wang and Zhang (2012a) Z. Wang and S.-C. Zhang, Phys. Rev. X 2, 031008 (2012a).
  * Wang and Yan (2013) Z. Wang and B. Yan, Journal of Physics: Condensed Matter 25, 155601 (2013).
  * Yoshida  _et al._ (2014) T. Yoshida, R. Peters, S. Fujimoto, and N. Kawakami, Phys. Rev. Lett. 112, 196404 (2014).
  * Lu  _et al._ (2013) F. Lu, J. Zhao, H. Weng, Z. Fang, and X. Dai, Phys. Rev. Lett. 110, 096401 (2013).
  * Deng  _et al._ (2013) X. Deng, K. Haule, and G. Kotliar, Phys. Rev. Lett. 111, 176404 (2013).
  * Hung  _et al._ (2013) H.-H. Hung, L. Wang, Z.-C. Gu, and G. A. Fiete, Phys. Rev. B 87, 121113 (2013).
  * Lang  _et al._ (2013) T. C. Lang, A. M. Essin, V. Gurarie, and S. Wessel, Phys. Rev. B 87, 205101 (2013).
  * Hung  _et al._ (2014) H.-H. Hung, V. Chua, L. Wang, and G. A. Fiete, Phys. Rev. B 89, 235104 (2014).
  * Meng  _et al._ (2014) Z. Y. Meng, H.-H. Hung, and T. C. Lang, Modern Physics Letters B 28, 1430001 (2014).
  * Grandi  _et al._ (2015a) F. Grandi, F. Manghi, O. Corradini, and C. M. Bertoni, Phys. Rev. B 91, 115112 (2015a).
  * Grandi  _et al._ (2015b) F. Grandi, F. Manghi, O. Corradini, C. M. Bertoni, and A. Bonini, New Journal of Physics 17, 023004 (2015b).
  * Budich  _et al._ (2012) J. C. Budich, R. Thomale, G. Li, M. Laubach, and S.-C. Zhang, Phys. Rev. B 86, 201407 (2012).
  * Budich  _et al._ (2013) J. C. Budich, B. Trauzettel, and G. Sangiovanni, Phys. Rev. B 87, 235104 (2013).
  * Amaricci  _et al._ (2015) A. Amaricci, J. C. Budich, M. Capone, B. Trauzettel, and G. Sangiovanni, Phys. Rev. Lett. 114, 185701 (2015).
  * Chen  _et al._ (2015) Y.-H. Chen, H.-H. Hung, G. Su, G. A. Fiete, and C. S. Ting, Phys. Rev. B 91, 045122 (2015).
  * Maier  _et al._ (2005) T. Maier, M. Jarrell, T. Pruschke, and M. H. Hettler, Rev. Mod. Phys. 77, 1027 (2005).
  * Wang  _et al._ (2011) Z. Wang, X.-L. Qi, and S.-C. Zhang, Phys. Rev. B 84, 014527 (2011).
  * Wang  _et al._ (2012) Z. Wang, X.-L. Qi, and S.-C. Zhang, Phys. Rev. B 85, 165126 (2012).
  * Wang and Zhang (2012b) Z. Wang and S.-C. Zhang, Phys. Rev. B 86, 165116 (2012b).
  * Assaad and Evertz (2008) F. Assaad and H. Evertz, in  _Computational Many-Particle Physics_, Lecture Notes in Physics, Vol. 739, edited by H. Fehske, R. Schneider, and A. Weiße (Springer Berlin Heidelberg, 2008) pp. 277–356.
  * Parcollet  _et al._ (2004) O. Parcollet, G. Biroli, and G. Kotliar, Phys. Rev. Lett. 92, 226402 (2004).
  * Biroli  _et al._ (2004) G. Biroli, O. Parcollet, and G. Kotliar, Phys. Rev. B 69, 205108 (2004).
  * Sakai  _et al._ (2012) S. Sakai, G. Sangiovanni, M. Civelli, Y. Motome, K. Held, and M. Imada, Phys. Rev. B 85, 035102 (2012).
  * Li  _et al._ (2015) Q.-X. Li, R.-Q. He, and Z.-Y. Lu, Phys. Rev. B 92, 155127 (2015).
  * Wu  _et al._ (2015) H.-Q. Wu, Y.-Y. He, Y.-Z. You, C. Xu, Z. Y. Meng, and Z.-Y. Lu, arXiv 1506, 00549 (2015).
  * Wu  _et al._ (2012) W. Wu, S. Rachel, W.-M. Liu, and K. Le Hur, Phys. Rev. B 85, 205102 (2012).