\documentclass[onecolumn,english,prl,showpacs,superscriptaddress]{revtex4}
\usepackage{epsfig}
\usepackage{graphicx}
\usepackage{amsfonts}
\usepackage[figuresright]{rotating}
\usepackage{amssymb}
\usepackage{amsmath}
%\usepackage{psfrag}
\usepackage{amsmath}
\usepackage{graphicx}

%%%%%%%%%% Start TeXmacs macros
\newcommand{\tmop}[1]{\ensuremath{\operatorname{#1}}}
%%%%%%%%%% End TeXmacs macros

\begin{document}

\title{Supplemental Material}
\author{Hsiang-Hsuan Hung}
\affiliation{Department of Physics, The University of Texas at
Austin, Austin, TX, 78712, USA }
\author{Lei Wang}
\affiliation{Theoretische Physik, ETH
Zurich, 8093 Zurich, Switzerland}
\author{Zheng-Cheng Gu}
\affiliation{Institute for Quantum Information, California Institute
of Technology, Pasadena, CA 91125, USA} \affiliation{Department of
Physics, California Institute of Technology, Pasadena, CA 91125,
USA}
\author{Gregory A. Fiete}
\affiliation{Department of Physics, The University of Texas at
Austin, Austin, TX, 78712, USA }

\date{\today}
\pacs{71.10.Fd,71.70.Ej}
%71.10.Fd   Lattice fermion models (Hubbard model, etc.)
%71.70.Ej       Spin-orbit coupling in in condensed matter,
\maketitle




%\bibliography{/Users/wanglei/Documents/Papers/papers}









\begin{appendix}
\begin{widetext}
\section{Sign-free determinant projector QMC}

%
%QMC stuff goes here
%

The determinant QMC has been shown to be an excellent and unbiased
approach to deal with strongly correlated system with Hubbard
interactions.\cite{sugiyama1986,sorella1989,white1989,scalettar1991,assaad1996,assaad1999,
meng2010,cai2012}
In the projector algorithm, the ground state wave function
$|\Psi_0\rangle$ can be obtained using standard projection
procedures on a trivial wave function $|\Psi_T\rangle$, as long as
one requires $\langle \Psi_T |\Psi_0\rangle \ne 0$. The expectation
value of an observable $A$ is obtained by
\begin{eqnarray}
\langle A\rangle =\lim_{\Theta \to \infty}\frac{\langle \Psi_T |e^{-
\frac{\Theta }{2}H} A e^{- \frac{\Theta
}{2}H}|\Psi_T\rangle}{\langle \Psi_T |e^{-\Theta H}| \Psi_T\rangle}.
\label{eqn:expection}
\end{eqnarray}
The projection operator $e^{-\Theta H}$ can be discretized into many
time slices $e^{-\Theta H}=[e^{-\Delta \tau H}]^M$ with
$\Theta=\Delta \tau M$ where $\Delta \tau \ll 1$ and $M$ is the
number of time slices with a large integer number; $e^{-\Delta \tau
H}=e^{-\Delta \tau (H_0+H_U)}$ is the imaginary time-evolution
propagator during $\Delta \tau$. The noninteracting ground state of
$H_0$ is a good candidate for the trial wave function
$|\Psi_T\rangle$. With this trial wave function, we have confirmed
that the determinant projector QMC is in a good agreement with our
exact diagonalization results on a $L\times L=3\times 3$ system. By
the first order Suzuki-Trotter decomposition, one can decompose
$e^{-\Delta \tau H}$ as
\begin{eqnarray}
e^{-\Delta \tau H} \simeq e^{-\Delta \tau H_0}e^{-\Delta \tau H_U},
\label{eqn:suzuki}
\end{eqnarray}
where $H_0$ is the tight-binding Hamiltonian of the generalized
Kane-Mele-Hubbard (KMH) model as shown in Eq. (1) of the main text.
$H_U=\frac{U}{2}\sum_i(n_i-1)^2$ involves 4 fermionic operators and
cannot be represented in terms of single-particle basis. However, by
the discrete $SU(2)$-invariant Hubbard-Stratonovich transformation,
\cite{assaad1998} the interacting imaginary time-evolution operator
$e^{-\Delta \tau H_U}$ (for $U>0$) can be decomposed as
\begin{eqnarray}
e^{-\Delta \tau
\frac{U}{2}(n_i-1)^2}=\frac{1}{4}\sum_{l=\pm1,\pm2}\gamma(l)e^{i\sqrt{\Delta
\tau \frac{U}{2}}\eta(l)(n_i-1)}+O(\Delta
\tau^4),\label{eqn:HStransformation}
\end{eqnarray}
where $\gamma(\pm1)=1+\sqrt{6}/3$, $\gamma(\pm2)=1-\sqrt{6}/3$;
$\eta(\pm1)=\pm\sqrt{2(3-\sqrt{6})}$ and
$\eta(\pm2)=\pm\sqrt{2(3+\sqrt{6})}$ are 4-component auxiliary
fields determined by Monte Carlo samplings. Ref.
[\onlinecite{assaad2002,assaad2003,hung2011}] provide pedagogical
introductions about the QMC method. In this work, we employ $\Delta
\tau t=0.05$ in all the QMC simulations.


In the determinant algorithm with the Suzuki-Trotter decomposition
Eq. (\ref{eqn:suzuki}) and the Hubbard-Stratonovich transformation
Eq. (\ref{eqn:HStransformation}), the denominator of Eq.
(\ref{eqn:expection}) reads as
\cite{sorella1989,assaad2002,hohenadler2012,zheng2011} (up to a
constant factor)
\begin{eqnarray}
\langle \Psi_T | e^{-\Theta H} |\Psi_T\rangle &=& \langle \Psi_T |
\prod^M_{\tau=1} e^{-\Delta \tau H_{\tau}} |\Psi_T\rangle=\langle
\Psi_T | \prod^M_{\tau=1} e^{-\Delta \tau H_0}e^{-\Delta \tau
H_{U,\tau}} |\Psi_T\rangle\\ \nonumber &=&\sum_{\lbrace l_{i,\tau}
\rbrace} \Big\{ \prod_{i,\tau} \gamma(l_{i,\tau}) \prod_{\sigma}
\textrm{Tr} \Big(\prod^M_{\tau=1} e^{-\Delta \tau
\sum_{i,j}c^{\dag}_{i,\sigma}[{\bf
H^{\sigma}_{0}}]_{ij}c_{j,\sigma}}e^{i\sqrt{\Delta \tau
\frac{U}{2}}\eta(l_{i,\tau})(n_{i,\sigma}-\frac{1}{2})}
\Big)\Big\}\\
&=&\sum_{\lbrace l_{i,\tau} \rbrace} \Big\{ \prod_{i,\tau}
\gamma(l_{i,\tau}) p[\lbrace \eta(l_i) \rbrace]\Big\}, \nonumber
\end{eqnarray}
where $\sum_{l_{i,\tau}}$ runs over possible auxiliary
configurations $\eta(l_{i,\tau})$, where $i=1-N$, $\tau=1-M$; ${\bf
H^{\sigma}_0}$ is the matrix kernel of $H_0$ with spin-$\sigma$. The
probability weight $p$ for a given auxiliary configuration $\lbrace
\eta(l_{i,\tau}) \rbrace$ is simply denoted as \cite{hirsch1985}
\begin{eqnarray}
p(\lbrace \eta \rbrace)=\det \Big(O_{\uparrow}[\eta(l_{i,\tau})]
\Big)\det \Big(O_{\downarrow}[\eta(l_{i,\tau})] \Big),
\end{eqnarray}
where $\det \Big( O_{\sigma}[\eta(l_{i,\tau})] \Big)=\textrm{Tr}
\Big(\prod^M_{\tau=1} e^{-\Delta \tau
\sum_{i,j}c^{\dag}_{i,\sigma}[{\bf
H^{\sigma}_{0}}]_{ij}c_{j,\sigma}}e^{i\sqrt{\Delta \tau
\frac{U}{2}}\eta(l_{i,\tau})(n_{i,\sigma}-\frac{1}{2})} \Big)$. When
$p<0$, QMC simulations meet notorious minus-sign problems.

It has been proven, that at half filling, there exists a
particle-hole symmetry in the Kane-Mele-Hubbard model ($H_0$ without
$t_{3N}$ terms), such that the probability is always
positive-definitive.\cite{zheng2011,hohenadler2012} To show the
positiveness of $p(\lbrace \eta \rbrace)$ in the generalized KMH
model, we employ the particle-hole transformation on the $H_0$ and
$H_U$ with $\downarrow$ but remain those with $\uparrow$ unchanged.
The particle-hole transformation acts as
\begin{eqnarray}
c_{i,\sigma} \to \xi_i d^{\dag}_{i,\sigma}, \ \ c^{\dag}_{i,\sigma}
\to \xi_i d_{i,\sigma},\nonumber
\end{eqnarray}
where $\xi_i=-1$ ($\xi_i=1$) if $i$ belongs to $A$ ($B$) sublattice.
In the Kane-Mele-Hubbard model, upon such a transformation, the
nearest-neighbor tight-binding term turns out to be
\begin{eqnarray} & &-tc^{\dag}_{i,\downarrow}c_{j,\downarrow}
-tc^{\dag}_{j,\downarrow}c_{i,\downarrow} \nonumber \\
&\to& -t\xi_i \xi_j d_{i,\downarrow}d^{\dag}_{j,\downarrow}-t\xi_i
\xi_j d_{j,\downarrow}d^{\dag}_{i,\downarrow}=-t(
d^{\dag}_{j,\downarrow}d_{i,\downarrow}+d^{\dag}_{i,\downarrow}d_{j,\downarrow}).\nonumber
\end{eqnarray}
Note that the $t$ hopping connects $A$ and $B$ sublattices, so we
have $(-1)\xi_i \xi_j=1$. The correspondence spin-up term on the
same $\langle i,j \rangle$ bond is
$-t(c^{\dag}_{i,\uparrow}c_{j,\uparrow}
+c^{\dag}_{j,\uparrow}c_{i,\uparrow})$. For the
second-nearest-neighbor hopping term with spin-down, the
particle-hole transformation acts it as
\begin{eqnarray} & & i\lambda_{SO}c^{\dag}_{i,\downarrow}c_{j,\downarrow}
-i\lambda_{SO}c^{\dag}_{j,\downarrow}c_{i,\downarrow} \nonumber \\
&\to& i\lambda_{SO}(-1)\xi_i \xi_j(
d^{\dag}_{j,\downarrow}d_{i,\downarrow}-d^{\dag}_{i,\downarrow}d_{j,\downarrow})=-i\lambda
_{SO}(
d^{\dag}_{j,\downarrow}d_{i,\downarrow}-d^{\dag}_{i,\downarrow}d_{j,\downarrow}).\nonumber
\end{eqnarray}
Here the $i\lambda_{SO}$ hopping connects the same sublattices, so
we have $(-1)\xi_i \xi_j=-1$. The correspondence spin-up term on the
same $\langle\langle i,j \rangle\rangle$ bond is
$-i\lambda_{SO}c^{\dag}_{i,\uparrow}c_{j,\uparrow}
+i\lambda_{SO}c^{\dag}_{j,\uparrow}c_{i,\uparrow}$. Thus, under the
particle-hole symmetry, the tight-binding matrix kernel with the
first and second nearest-neighbor hopping transforms as ${\bf
H^{\downarrow}_0} \to {\bf \tilde{H}^{\downarrow}_0}={\bf
H^{\uparrow}_0}^*$, and thus the Kane-Mele-Hubbard model is
particle-hole symmetric. For the third-nearest-neighbor hopping
$t_{3N}$, the particle-hole transformation provides
\begin{eqnarray} & &-t_{3N}c^{\dag}_{i,\downarrow}c_{j,\downarrow} \nonumber \\
&\to& -t_{3N}\xi_i \xi_j
d_{i,\downarrow}d^{\dag}_{j,\downarrow}=-t_{3N}
d^{\dag}_{j,\downarrow}d_{i,\downarrow}\nonumber
\end{eqnarray}
unchanged since the $t_{3N}$ hopping connects $A$ and $B$
sublattices. Therefore, the real-valued third-neighbor hopping
$t_{3N}$ in the generalized KMH model remains $H_0$ particle-hole
symmetric at half-filling.



The Hubbard interaction $H_U$ on $\downarrow$ transforms as
\begin{eqnarray}
& & i\sqrt{\Delta \tau
\frac{U}{2}}\eta(l_{i,\tau})(n_{i,\downarrow}-\frac{1}{2}) \nonumber \\
&\to& i\sqrt{\Delta \tau
\frac{U}{2}}\eta(l_{i,\tau}) \Big\{ (\xi_i)^2
d_{i,\downarrow}d^{\dag}_{i,\downarrow}-\frac{1}{2} \Big\} \nonumber \\
&=&-i\sqrt{\Delta \tau
\frac{U}{2}}\eta(l_{i,\tau})(d^{\dag}_{i,\downarrow}d_{i,\downarrow}-\frac{1}{2}),\nonumber
\end{eqnarray}
which is the complex conjugate of $H_U$ on $\uparrow$. Consequently,
upon the particle-hole symmetry, one can have $\det(O_{\downarrow})=
\det( O_{\uparrow} )^*$ and the probability weight
$p=\det({O_{\uparrow}})\det({O_{\downarrow}})=|\det(O_{\uparrow})|^2$
being real positive. The QMC simulation in the half-filled
generalized KMH model is sign-free and numerically exact.

\section{single particle Green's functions and $Z_2$ invariant}

Without sign problems, the QMC samplings provide highly accurate not
only in equal-time Green's functions but also in time-displaced
Green's functions\cite{assaad1996,feldbacher2001}
\begin{eqnarray}
G_{\sigma}(\vec{r},\tau)=\langle \Psi_0 |
c_{\sigma}(\vec{r},\tau)c^{\dag}_{\sigma}(0)|\Psi_0\rangle,\nonumber
\end{eqnarray}where $\tau >0$. By performing double
Fourier transformation we obtain the Green's functions in momentum
space and with Matsubara frequency, i. e.
$G_{\sigma}(\mathbf{k},i\omega_{n})$.


\begin{figure}[htb!]
\epsfig{file=BZ2.eps,clip=0.1,width=0.32\linewidth,angle=0}
\caption{(Color online)  The Brillouin zone of the honeycomb
lattice. The time-reversal invariant momentum (TRIM) points labeled
by the green dots are $\Gamma=(0,0)$,
$M_{1,2}=(\pm\frac{\pi}{\sqrt{3}a},\frac{\pi}{3a})$, and
$M_{3}=(0,\frac{2\pi}{3a})$. The open and solid circles denote
graphene Dirac points $K_{1,2}=(\pm\frac{4\pi}{3\sqrt{3}a},0)$}.
\label{fig:bandstructure}
\end{figure}
It has been shown that the zero frequency Green's functions are able
to evaluate the  $Z_2$ invariant index in the interacting
case.\cite{wang2012} The $Z_2$ invariant is determined by the parity
of the eigenvectors of the inverse Green's functions on
time-reversal invariant momentum (TRIM) points which are obtained
from
\begin{eqnarray}
[G({\mathbf k}_{i},0)]^{-1}|\mu_i \rangle=\mu_i
|\mu_i\rangle.\nonumber
\end{eqnarray}
Note that since there still exists an inversion symmetry in the
generalized KMH model, the inverse Green's functions and the parity
operator have simultaneous eigenvectors, \textit{i.e.}
$P|\mu_{i}\rangle=\eta_{\mu_{i}}|\mu_{i}\rangle$. In the
(generalized) KMH model,  the parity operator exchanges $A$, $B$
sublattices independent of spin index. Therefore, with the spinor
convention $\Psi^{\dag}=(c^{\dag}_{A,\uparrow} \
c^{\dag}_{B,\uparrow} \ c^{\dag}_{A,\downarrow} \
c^{\dag}_{B,\downarrow})$, the parity operator is defined as $P=I
\otimes \sigma^{x}$.\cite{fu2007} In the QMC simulations, the
particle-hole symmetry provides $G_{\uparrow}({\mathbf
k_i},0)=G_{\downarrow}({\mathbf k_i},0)$, while  ${\mathbf k}_{i}$
is at TRIM, i.e., ${\mathbf k}=-{\mathbf k}+{\mathbf G}$ for a
reciprocal vectors ${\mathbf G}$. Therefore, we can directly
diagonalize $G_{\sigma}({\mathbf k}_{i},0) =[-H_{\mathbf
k}-\Sigma({\mathbf k}_{i},0)]^{-1}$ instead of inverse Green's
functions for all ${\mathbf k}_{i}\in$ TRIM points
\begin{eqnarray}
G_{\sigma}({\mathbf k}_{i},0)|\tilde{\mu_i} \rangle=\tilde{\mu_i}
|\tilde{\mu_i}\rangle \nonumber,
\end{eqnarray} and choose the eigenvectors associated with positive
eigenvalues ($\tilde{\mu_i}>0$, denoting occupied bands and are
called right-zero \cite{Wang:prx12}). In the honeycomb lattice, the
TRIM points are $\Gamma$, $M_{1,2,3}$ as depicted in Fig.
\ref{fig:bandstructure}. Then we can employ the formalism proposed
by Fu and Kane\cite{fu2007,Wang:prx12} to evaluate the $Z_2$
invariant as
\begin{eqnarray}
(-1)^{\nu}=\prod_{{\mathbf k_i} \in TRIM} \tilde{\eta}_{\mu_i},
%\label{eqn:z2invariant}
\end{eqnarray}
where $\tilde{\eta}_{\mu_i}=\langle \tilde{\mu_i}
|P|\tilde{\mu_i}\rangle$. When $\nu=0$ for a trivial insulator,
whereas $\nu=1$ for a $Z_2$ topological insulator. In the case of
$U=0$, $\tilde{\eta}_{\mu_i}=\pm1$. In the cases of finite $U$, we
find that $\langle \tilde{\eta}_{\mu_i} \rangle=\pm1$ can be still
obtained by sufficient QMC simulations. As $t_{3N}$ approaches the
topological critical point, $(-1)^{\nu}$ will be smeared out and is
laid between $\pm1$. In this case, more QMC samplings are required
for more accurate values.


Note that since  $G_{\uparrow}({\mathbf
k}_{i},0)=G_{\downarrow}({\mathbf k}_{i},0)$, and $G({\mathbf
k}_{i},0)$ [$=G_{\uparrow}({\mathbf k}_{i},0) \oplus
G_{\downarrow}({\mathbf k}_{i},0)$] and $P$ ($=I \otimes \sigma^x$)
have the simultaneous eigenvector sets, one has a relation:
\begin{eqnarray}
G_{\uparrow,\downarrow}({\mathbf k}_{i},0)=\alpha_{\mathbf k_i}
\sigma^x. \label{eqn:gamma}
\end{eqnarray}
In the context we show that in addition to the $Z_2$ invariant, the
proportional coefficient $\alpha_{\mathbf k}$ also plays another
role to characterize the $Z_2$ topological insulator/trivial
insulator phase transition and even is more sensitive than $\nu$
numerically. Upon the topological phase transition, the bulk gap
will close at the TRIM points. Thus, the zero-frequency
single-particle Green's functions are divergent on the poles.
\cite{gurarie2011}




The relation Eq. (\ref{eqn:gamma}) should be expected both in the
noninteracting and interacting cases. However, as $U \ne 0$ Eq.
(\ref{eqn:gamma}) is not guaranteed in a single measurement in the
QMC simulations. The proportionality relation between the
zero-frequency Green's functions and the parity matrix $\sigma^x$
can be recovered only upon enough samplings. To interpret this, we
present the $6 \times 6$ benchmark results for the matrix elements
of the zero-frequency Green's functions at ${\mathbf k_i=M_1}$ as a
function of the number of measurements in Figs. \ref{fig:benchmark}.
\begin{figure}[!htb]
\epsfig{file=benchmark1.eps,clip=0.1,width=0.3\linewidth,angle=0}
\epsfig{file=benchmark2.eps,clip=0.1,width=0.24\linewidth,angle=0}
\caption{(Color online) The matrix elements of the zero-frequency
Green functions $G({\mathbf M_1},0)$ vs the number of samplings $m$
at (a) $t_{3N}=0.32t$ and (b) $t_{3N}=0.37t$. $\lambda_{SO}=0.4t$
and $U=4t$. $\textrm{Re}[g_{ij}]$ and $\textrm{Im}[g_{ij}]$ denote
the real part and imaginary part of $[G({\mathbf M_1},0)]_{ij}$,
respectively; $||g_{ii}||$ denotes the diagonal component of
$G({\mathbf M_1},0)$ in magnitudes.} \label{fig:benchmark}
\end{figure}
$g_{ij}=[G({\mathbf M_1},0)]_{ij}$ and $m$ denotes the number of
measurements. $\lambda_{SO}=0.4t$ and $U=4t$ are used. In this case,
the topological phase boundary is identified at $t_{3N}=0.348t$. We
choose the value of $t_{3N}$ close to the critical point. Fig.
\ref{fig:benchmark} (a) shows $t_{3N}=0.32t$ in the $Z_2$
topological insulator phase and (b) for $t_{3N}=0.37t$ in the
trivial insulator. From the panels, it is evident that the structure
of the Green's function does not satisfy Eq. (\ref{eqn:gamma})
without sufficient samplings. At small $m$, the real parts of
$g_{12}$ and $g_{21}$ are not equal; furthermore, $g_{12}$ and
$g_{21}$ have imaginary parts, and both of $g_{11}$ and $g_{22}$ are
finite. However, one can see that, upon sampling sufficient times,
$\textrm{Re}[g_{12}] \simeq \textrm{Re}[g_{21}]$, and meanwhile
$\textrm{Im}[g_{12}]$, $\textrm{Im}[g_{21}]$, $||g_{11(22)}||$ go to
zero. Thus, in the $m \to \infty$ limit, Eq. (\ref{eqn:gamma}) is
recovered. Also note that, $\alpha_{\mathbf
M_1}=\textrm{Re}[g_{12}]$ in both cases indicates opposite sign as
observed by the signature of the topological phase transition.
Moreover, by such $m$ scaling, we also confirm that, the value of
the $Z_2$ invariant also monotonically approaches to $\pm 1$. In our
paper, we choose the value of $m$ large enough to determine the
$\sigma^x$ structure and extract the coefficients, i.e.,
$\alpha_{{\mathbf k_i}}$ in Eq. (\ref{eqn:gamma}).




\section{Critical Hubbard interactions for antiferromagnetism}

In the generalized KMH model,  a strong Hubbard interaction can also
derive the antiferromagnetic (AF) ordering, due to the bipartite
lattice structure. Similarly to the KMH model (with $t_{3N}=0$)
\cite{rachel2010,zheng2011,hohenadler2011},  in the generalized KMH
model, finite values of $\lambda_{SO}$ also break the $SU(2)$
symmetry down to the  $U(1)$ symmetry and the dominant magnetism
behavior lies on x-y plane. The planar spin structure factor can be
defined as\cite{hohenadler2011,zheng2011}
\begin{eqnarray}
S_{AF}=\sum_{\vec{r},\vec{r}_j}(-1)^{\vec{r}_i+\vec{r}_j}\langle
S^{+}_iS^-_j +S^{-}_iS^+_j\rangle.\nonumber
\end{eqnarray}
$(-1)^{\vec{r}_i}=1(-1)$ for $i \in A(B)$ sublattice. This is
similar to determining the N\'eel type ordering using the
antiferromagnetic spin structure factor at ${\mathbf k}=(\pi,\pi)$
in a square lattice.

To identify whether there exists the antiferromagnetism in the
thermodynamic limit, we study the finite-size scaling behavior of
$S_{AF}$ at $L \to \infty$. Generally speaking, the spin-orbital
coupling will suppress AF ordering, and larger $\lambda_{SO}$'s are
associated with larger $U_c$'s to induce the AF ordering. Note that,
due to the presence of third nearest neighboring hopping $t_{3N}$
which favors the N\'eel pattern in the second order perturbation,
the threshold interaction $U_c$ in the generalized KMH model is
smaller than that in the KMH model.
\begin{figure}[!tb]
\epsfig{file=AF_L0.4U4.eps,clip=0.1,width=0.235\linewidth,angle=0}
\epsfig{file=AF_L0.4U5.eps,clip=0.1,width=0.23\linewidth,angle=0}
\epsfig{file=AF_L0.4U6.eps,clip=0.1,width=0.23\linewidth,angle=0}
\epsfig{file=AF_L1U6.eps,clip=0.1,width=0.235\linewidth,angle=0}
\caption{(Color online) (a)-(c) The finite size scaling of the
antiferromagnetic spin structure factor $S_{AF}/N$ vs $1/L$ at
$\lambda_{SO}=0.4t$ and different $U=4t,5t,6t$. (d) $S_{AF}/N$ vs
$1/L$ at $\lambda_{SO}=t$ and $U=6t$. Here $N=2\times L^2$.}
\label{fig:criticalU}
\end{figure}

The QMC results on $S_{AF}/N$ vs $1/L$ are shown in Figs.
\ref{fig:criticalU}. In (a), we can see that, for
$\lambda_{SO}=0.4t$,
 $U=4t$ is not sufficiently large to induce the AF
ordering. At $U=5t$, $S_{AF}$ is enhanced and the $U$ value is close
to the critical value to drive the AF ordering. In (c) under the
interaction $U=6t$, $S_{AF}$ saturates to a finite value at $1/L \to
0$, suggesting that the AF ordering exists in the thermodynamic
limit. Fig. \ref{fig:criticalU} (d) depicts the case of $U=6t$ but
at $\lambda_{SO}=t$. Compared to (c), where an AF ordering is
induced, the structure factor in (d) still goes to zero in the $L
\to \infty$ limit. Thus, stronger spin-orbital couplings obviously
suppress the existence of AF ordering and raise values of critical
interactions $U_c$ in the generalized KMH model.


\section{Single-particle excitation}

In this subsection, we present the approach to evaluate  the
single-particle excitation (charge gap) $\Delta_c$ in the QMC
simulations. The charge gap is defined as the energy cost to add a
particle into (or remove a particle from) the system composed of
$N_e$ fermions. Supposed that,  $\hat{H} |
\Psi^{N_e+1}_n\rangle=E^{N_e+1}_n | \Psi^{N_e+1}_n\rangle$ and
$\hat{H} | \Psi^{N_e}_n\rangle=E^{N_e}_n | \Psi^{N_e}_n\rangle$,
then the charge gap reads $\Delta_c \equiv E^{N_e+1}_0-E^{N_e}_0$.
It can be obtained via calculating the on-site time-displaced
Green's functions which are written as
\begin{eqnarray}
G(\vec{r}=0,\tau)&=&\frac{1}{N}\sum_{i,\sigma}G_{\sigma}(i,i;\tau)
\nonumber \\ &=&\frac{1}{N}\sum_{i,\sigma}\langle \Psi^{N_e}_0 |
c_{\sigma}(i,\tau)c^{\dag}_{\sigma}(i)|\Psi^{N_e}_0 \rangle \nonumber \\
&=& \frac{1}{N}\sum_{i,\sigma}\langle \Psi^{N_e}_0 | e^{\tau
\hat{H}}c_{\sigma}(i)e^{-\tau \hat{H}}c^{\dag}_{\sigma}(i)|\Psi^{N_e}_0 \rangle. \nonumber
\\
&=& \frac{1}{N}\sum_{n,i,\sigma} e^{-\tau (E^{N_e+1}_n-E^{N_e}_0)}
|\langle \Psi^{N_e}_0 | c_{\sigma}(i) |\Psi^{N_e+1}_n
\rangle|^2.\nonumber
\end{eqnarray}
Therefore, at large $\tau$, we have $G(\vec{r}=0,\tau)
\sim e^{-\tau \Delta_c}$
 and
then one can find the slope of  $\ln{G(\vec{r}=0,\tau)}$ at large
$\tau$ to determine the value of $\Delta_c$. Refs.
[\onlinecite{assaad1996,hung2011,feldbacher2001,assaad1996b}]
provide the detailed descriptions. The evaluation of the excitation
by the on-site single-particle Green's function can determine the
value of the single-particle excitation without concerning about
specific momentum points, e.g. $\Delta_c({\mathbf k})$. (Note that,
in the noninteracting limit, the gap of the KMH model with
$\lambda=0$ closes at the Dirac points $K_{1,2}$, whereas the gap of
the generalized KMH model with $t_c$ closes at $M_{1,2,3}$.)

\end{widetext}
\end{appendix}


%\bibliography{KM}

\begin{thebibliography}{23}
\expandafter\ifx\csname
natexlab\endcsname\relax\def\natexlab#1{#1}\fi
\expandafter\ifx\csname bibnamefont\endcsname\relax
  \def\bibnamefont#1{#1}\fi
\expandafter\ifx\csname bibfnamefont\endcsname\relax
  \def\bibfnamefont#1{#1}\fi
\expandafter\ifx\csname citenamefont\endcsname\relax
  \def\citenamefont#1{#1}\fi
\expandafter\ifx\csname url\endcsname\relax
  \def\url#1{\texttt{#1}}\fi
\expandafter\ifx\csname urlprefix\endcsname\relax\def\urlprefix{URL
}\fi \providecommand{\bibinfo}[2]{#2}
\providecommand{\eprint}[2][]{\url{#2}}

\bibitem[{\citenamefont{Sugiyama and Koonin}(1986)}]{sugiyama1986}
\bibinfo{author}{\bibfnamefont{G.}~\bibnamefont{Sugiyama}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{S.~E.} \bibnamefont{Koonin}},
  \bibinfo{journal}{Ann. Phys.} \textbf{\bibinfo{volume}{168}},
  \bibinfo{pages}{1} (\bibinfo{year}{1986}).

\bibitem[{\citenamefont{Sorella et~al.}(1989)\citenamefont{Sorella, Baroni,
  Car, and Parrinello}}]{sorella1989}
\bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{Sorella}},
  \bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{Baroni}},
  \bibinfo{author}{\bibfnamefont{R.}~\bibnamefont{Car}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Parrinello}},
  \bibinfo{journal}{Europhys. Lett.} \textbf{\bibinfo{volume}{8}},
  \bibinfo{pages}{663} (\bibinfo{year}{1989}).

\bibitem[{\citenamefont{White et~al.}(1989)\citenamefont{White, Scalapino,
  Sugar, Loh, Gubernatis, and Scalettar}}]{white1989}
\bibinfo{author}{\bibfnamefont{S.~R.} \bibnamefont{White}},
  \bibinfo{author}{\bibfnamefont{D.~J.} \bibnamefont{Scalapino}},
  \bibinfo{author}{\bibfnamefont{R.~L.} \bibnamefont{Sugar}},
  \bibinfo{author}{\bibfnamefont{E.~Y.} \bibnamefont{Loh}},
  \bibinfo{author}{\bibfnamefont{J.~E.} \bibnamefont{Gubernatis}},
  \bibnamefont{and} \bibinfo{author}{\bibfnamefont{R.~T.}
  \bibnamefont{Scalettar}}, \bibinfo{journal}{Phys. Rev. B}
  \textbf{\bibinfo{volume}{40}}, \bibinfo{pages}{506} (\bibinfo{year}{1989}).

\bibitem[{\citenamefont{Scalettar et~al.}(1991)\citenamefont{Scalettar,
  Scalapino, Sugar, and White}}]{scalettar1991}
\bibinfo{author}{\bibfnamefont{R.~T.} \bibnamefont{Scalettar}},
  \bibinfo{author}{\bibfnamefont{D.~J.} \bibnamefont{Scalapino}},
  \bibinfo{author}{\bibfnamefont{R.~L.} \bibnamefont{Sugar}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{S.~R.} \bibnamefont{White}},
  \bibinfo{journal}{Phys. Rev. B} \textbf{\bibinfo{volume}{44}},
  \bibinfo{pages}{770} (\bibinfo{year}{1991}).

\bibitem[{\citenamefont{Assaad and Imada}(1996{\natexlab{a}})}]{assaad1996}
\bibinfo{author}{\bibfnamefont{F.~F.} \bibnamefont{Assaad}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Imada}}, \bibinfo{journal}{J.
  Phys. Soc. Jpn.} \textbf{\bibinfo{volume}{65}}, \bibinfo{pages}{189}
  (\bibinfo{year}{1996}{\natexlab{a}}).

\bibitem[{\citenamefont{Assaad}(1999)}]{assaad1999}
\bibinfo{author}{\bibfnamefont{F.~F.} \bibnamefont{Assaad}},
  \bibinfo{journal}{Phys. Rev. Lett.} \textbf{\bibinfo{volume}{83}},
  \bibinfo{pages}{796} (\bibinfo{year}{1999}).

\bibitem[{\citenamefont{Meng et~al.}(2010)\citenamefont{Meng, Lang, Wessel,
  Assaad, and Muramatsu}}]{meng2010}
\bibinfo{author}{\bibfnamefont{Z.~Y.} \bibnamefont{Meng}},
  \bibinfo{author}{\bibfnamefont{T.~C.} \bibnamefont{Lang}},
  \bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{Wessel}},
  \bibinfo{author}{\bibfnamefont{F.~F.} \bibnamefont{Assaad}},
  \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Muramatsu}},
  \bibinfo{journal}{Nature} \textbf{\bibinfo{volume}{464}},
  \bibinfo{pages}{847} (\bibinfo{year}{2010}).

\bibitem[{\citenamefont{Cai et~al.}(2012)\citenamefont{Cai, Hung, Wang, Li, and
  Wu}}]{cai2012}
\bibinfo{author}{\bibfnamefont{Z.}~\bibnamefont{Cai}},
  \bibinfo{author}{\bibfnamefont{H.-H.} \bibnamefont{Hung}},
  \bibinfo{author}{\bibfnamefont{L.}~\bibnamefont{Wang}},
  \bibinfo{author}{\bibfnamefont{Y.}~\bibnamefont{Li}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{C.}~\bibnamefont{Wu}}, \bibinfo{journal}{ArXiv
  e-prints:1207.6843}  (\bibinfo{year}{2012}).

\bibitem[{\citenamefont{Assaad}(1998)}]{assaad1998}
\bibinfo{author}{\bibfnamefont{F.~F.} \bibnamefont{Assaad}},
  \bibinfo{journal}{ArXiv e-prints: cond-mat/9806307}  (\bibinfo{year}{1998}).

\bibitem[{\citenamefont{Assaad}(2002)}]{assaad2002}
\bibinfo{author}{\bibfnamefont{F.~F.} \bibnamefont{Assaad}},
  \emph{\bibinfo{title}{Quantum {Monte Carlo} methods on lattices: The
  determinantal approach in Quantum Simulations of Complex Many-Body Systems:
  >From Theory to Algorithms, Lecture Notes}} (\bibinfo{publisher}{NIC Series
  Vol. {\bf 10}}, \bibinfo{year}{2002}).

\bibitem[{\citenamefont{Assaad}(2003)}]{assaad2003}
\bibinfo{author}{\bibfnamefont{F.~F.} \bibnamefont{Assaad}},
  \bibinfo{journal}{AIP Conf. Proc.} \textbf{\bibinfo{volume}{678}},
  \bibinfo{pages}{117} (\bibinfo{year}{2003}).

\bibitem[{\citenamefont{Hung}(2011)}]{hung2011}
\bibinfo{author}{\bibfnamefont{H.-H.} \bibnamefont{Hung}},
  \emph{\bibinfo{title}{Exotic quantum magnetism and superfluidity in optical
  lattices}} (\bibinfo{publisher}{PhD thesis, University of California, San
  Diego}, \bibinfo{year}{2011}).

\bibitem[{\citenamefont{Zheng et~al.}(2011)\citenamefont{Zheng, Zhang, and
  Wu}}]{zheng2011}
\bibinfo{author}{\bibfnamefont{D.}~\bibnamefont{Zheng}},
  \bibinfo{author}{\bibfnamefont{G.-M.} \bibnamefont{Zhang}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{C.}~\bibnamefont{Wu}}, \bibinfo{journal}{Phys.
  Rev. B} \textbf{\bibinfo{volume}{84}}, \bibinfo{pages}{205121}
  (\bibinfo{year}{2011}).

\bibitem[{\citenamefont{Hohenadler et~al.}(2012)\citenamefont{Hohenadler, Meng,
  Lang, Wessel, Muramatsu, and Assaad}}]{hohenadler2012}
\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Hohenadler}},
  \bibinfo{author}{\bibfnamefont{Z.~Y.} \bibnamefont{Meng}},
  \bibinfo{author}{\bibfnamefont{T.~C.} \bibnamefont{Lang}},
  \bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{Wessel}},
  \bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Muramatsu}},
  \bibnamefont{and} \bibinfo{author}{\bibfnamefont{F.~F.}
  \bibnamefont{Assaad}}, \bibinfo{journal}{Phys. Rev. B}
  \textbf{\bibinfo{volume}{85}}, \bibinfo{pages}{115132}
  (\bibinfo{year}{2012}).

\bibitem[{\citenamefont{Hirsch}(1985)}]{hirsch1985}
\bibinfo{author}{\bibfnamefont{J.~E.} \bibnamefont{Hirsch}},
  \bibinfo{journal}{Phys. Rev. B} \textbf{\bibinfo{volume}{31}},
  \bibinfo{pages}{4403} (\bibinfo{year}{1985}).

\bibitem[{\citenamefont{Feldbacher and Assaad}(2001)}]{feldbacher2001}
\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Feldbacher}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{F.~F.} \bibnamefont{Assaad}},
  \bibinfo{journal}{Phys. Rev. B} \textbf{\bibinfo{volume}{63}},
  \bibinfo{pages}{073105} (\bibinfo{year}{2001}).

\bibitem[{\citenamefont{Wang and Zhang}(2012{\natexlab{a}})}]{wang2012}
\bibinfo{author}{\bibfnamefont{Z.}~\bibnamefont{Wang}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{S.-C.} \bibnamefont{Zhang}},
  \bibinfo{journal}{Phys. Rev. X} \textbf{\bibinfo{volume}{2}},
  \bibinfo{pages}{031008} (\bibinfo{year}{2012}{\natexlab{a}}).

\bibitem[{\citenamefont{Fu and Kane}(2007)}]{fu2007}
\bibinfo{author}{\bibfnamefont{L.}~\bibnamefont{Fu}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{C.~L.} \bibnamefont{Kane}},
  \bibinfo{journal}{Phys. Rev. B} \textbf{\bibinfo{volume}{76}},
  \bibinfo{pages}{045302} (\bibinfo{year}{2007}).

\bibitem[{\citenamefont{Wang and Zhang}(2012{\natexlab{b}})}]{Wang:prx12}
\bibinfo{author}{\bibfnamefont{Z.}~\bibnamefont{Wang}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{S.-C.} \bibnamefont{Zhang}},
  \bibinfo{journal}{Phys. Rev. X} \textbf{\bibinfo{volume}{2}},
  \bibinfo{pages}{031008} (\bibinfo{year}{2012}{\natexlab{b}}).

\bibitem[{\citenamefont{Gurarie}(2011)}]{gurarie2011}
\bibinfo{author}{\bibfnamefont{V.}~\bibnamefont{Gurarie}},
  \bibinfo{journal}{Phys. Rev. B} \textbf{\bibinfo{volume}{83}},
  \bibinfo{pages}{085426} (\bibinfo{year}{2011}).

\bibitem[{\citenamefont{Rachel and Le~Hur}(2010)}]{rachel2010}
\bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{Rachel}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{K.}~\bibnamefont{Le~Hur}},
  \bibinfo{journal}{Phys. Rev. B} \textbf{\bibinfo{volume}{82}},
  \bibinfo{pages}{075106} (\bibinfo{year}{2010}).

\bibitem[{\citenamefont{Hohenadler et~al.}(2011)\citenamefont{Hohenadler, Lang,
  and Assaad}}]{hohenadler2011}
\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Hohenadler}},
  \bibinfo{author}{\bibfnamefont{T.~C.} \bibnamefont{Lang}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{F.~F.} \bibnamefont{Assaad}},
  \bibinfo{journal}{Phys. Rev. Lett.} \textbf{\bibinfo{volume}{106}},
  \bibinfo{pages}{100403} (\bibinfo{year}{2011}).

\bibitem[{\citenamefont{Assaad and Imada}(1996{\natexlab{b}})}]{assaad1996b}
\bibinfo{author}{\bibfnamefont{F.~F.} \bibnamefont{Assaad}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Imada}},
  \bibinfo{journal}{Phys. Rev. Lett.} \textbf{\bibinfo{volume}{76}},
  \bibinfo{pages}{3176} (\bibinfo{year}{1996}{\natexlab{b}}).

\end{thebibliography}




\end{document}
