Dynamical Generation of Topological Masses in Dirac Fermions

TL;DR

This paper reports the first unbiased large-scale numerical evidence of a topological Mott insulator in strongly-correlated Dirac semimetals, showing interaction-driven topological mass generation and a quantum phase transition in the (2+1)D N=8 chiral Ising class.

Contribution

It provides the first unbiased numerical demonstration of a topological Mott insulator in Dirac semimetals with interactions mediated by Ising spins, revealing a dynamical topological mass generation.

Findings

Topological mass term is dynamically generated.

Quantum phase transition belongs to (2+1)D N=8 chiral Ising universality class.

Large scale sign-free quantum Monte Carlo simulations confirm the results.

Abstract

We report discovery of a topological Mott insulator in strongly-correlated Dirac semimetals. Such an interaction-driven topological state has been theoretically proposed but not yet observed with unbiased large scale numerical simulations. In our model, interactions between electrons are mediated by Ising spins in a transverse field. The results indicate that the topological mass term is dynamically generated and the resulting quantum phase transition belongs to the (2+1)D $N=8$ chiral Ising universality class. These conclusions stem from large scale sign free quantum Monte Carlo simulations.

Tables

Table 1. Table 1: Examples of the acceptance (1) or rejectance (0) for cluster updates (first column), number of spin to be flipped (second column), size of Wolff cluster (third column) and size of the geometric cluster (fourth column).

Accept/Reject	$n_{\text{flip}}$	$n_{\text{Wolff}}$	$n_{\text{geometric}}$
0	2400	2398	2
0	174	170	4
1	993	11	982
1	119	7	112
1	74	46	28
1	3436	2	3434
1	330	4	326
0	3421	3411	10
1	1203	23	1180
1	38	24	14

Equations

$$H$$

$$H_{\text{Fermion}}$$

$$H_{\text{Ising}}$$

$$H_{\text{Coupling}}$$

$$H$$

$$H_{\text{Fermion}}$$

$$H_{\text{Ising}}$$

$$H_{\text{Coupling}}$$

$$\displaystyle\mathcal{Z}=\langle\psi_{T}|Tr_{\{\mathbf{s}\}}\big{\{}e^{-\Theta H}\big{\}}|\psi_{T}\rangle$$

$$\mathcal{Z}$$

$$\mathcal{W}_{\text{Ising}}=\lambda^{MN_{s}}\exp\Big{[}-\Big{(}-\Delta\tau J\sum_{\ell=1}^{M}\sum_{\langle pq\rangle}s_{p,\ell}^{z}s_{q,\ell}^{z}-\gamma\sum_{p}\sum_{\ell=1}^{M}s_{p,\ell+1}^{z}s_{p,\ell}^{z}\Big{)}\Big{]}$$

$$\displaystyle\left\{\begin{array}[]{ll}\lambda=\sqrt{\sinh(\Delta\tau h)\cosh(\Delta\tau h)}\\ \gamma=-\frac{1}{2}\ln[\tanh(\Delta\tau h)]\end{array}.\right.$$

$$B_{\ell }^{\sigma }=e^{{\mathcal{C}}_{\ell ,4}}{e}^{{\mathcal{C}}_{\ell ,3}}{e}^{{\mathcal{C}}_{\ell ,2}}{e}^{{\mathcal{C}}_{\ell ,1}}{e}^{{\mathcal{K}}_2}{e}^{{\mathcal{K}}_1},$$

$${\beta }^{\prime }{\hat{H}}^{\prime }=-\Delta \tau J\sum _{\ell =1}^M\sum _{⟨pq⟩}{s}_{p,\ell }^z{s}_{q,\ell }^z-\gamma \sum _p\sum _{\ell =1}^M{s}_{p,\ell +1}^z{s}_{p,\ell }^z,$$

$$\pi (a)P(a\to b)\mathcal{A}(a\to b)=\pi (b)P(b\to a)\mathcal{A}(b\to a),$$

$$\frac{P(b\to a)}{P(a\to b)}=\frac{{\pi }_s(a)}{{\pi }_s(b)}.$$

$$\displaystyle\mathcal{A}(a\to b)=\min\Big{\{}\frac{P(b\to a)}{P(a\to b)}\frac{\pi(b)}{\pi(a)},1\Big{\}}=\min\Big{\{}\frac{\pi_{f}(b)}{\pi_{f}(a)},1\Big{\}}.$$

$$\displaystyle\mathcal{A}(a_{1}\to a_{n})=\mathcal{A}(a_{1}\to a_{2})\mathcal{A}(a_{2}\to a_{3})\cdots\mathcal{A}(a_{n-1}\to a_{n})=\min\Big{\{}\frac{\pi_{f}(a_{n})}{\pi_{f}(a_{1})},1\Big{\}}.$$

$$U_2$$

$$R_{Corr}$$

$$S^{\text{Ising}}(\mathbf{k})=\frac{1}{2N}\sum _{a=A,B}\sum _{mn}{e}^{-i\mathbf{k}\cdot ({\mathbf{R}}_m-{\mathbf{R}}_n)}⟨s_{ma}^z{s}_{na}^z⟩,$$

$$\displaystyle S^{\text{Ising}}(\boldsymbol{\Gamma})/L^{2}=\frac{1}{2N^{2}}\sum_{a=A,B}\sum_{mn}\langle s_{ma}^{z}s_{na}^{z}\rangle=\frac{1}{2}\Big{\langle}\Big{(}\frac{1}{N}\sum_{p\in A}s_{p}^{z}\Big{)}^{2}+(\frac{1}{N}\sum_{p\in B}s_{p}^{z}\Big{)}^{2}\Big{\rangle}=\frac{1}{2}\Big{(}\langle m_{A}^{2}\rangle+\langle m_{B}^{2}\rangle\Big{)},$$

$$\displaystyle\langle m^{2}\rangle=\frac{1}{4}\Big{(}\langle m_{A}^{2}\rangle+\langle m_{B}^{2}\rangle+2\langle m_{A}m_{B}\rangle\Big{)}.$$

$${\hat{H}}_{\text{Fermion}}=-t\sum _{⟨ij⟩\sigma }(e^{+i\sigma \varphi }{c}_{i\sigma }^{†}{c}_{j\sigma }+e^{-i\sigma \varphi }{c}_{j\sigma }^{†}{c}_{i\sigma })+\sum _{⟨⟨ij⟩⟩\sigma }{t}_{ij}^{\prime }\xi (c_{i\sigma }^{†}{c}_{j\sigma }+c_{j\sigma }^{†}{c}_{i\sigma }).$$

$$[{\mathbf{G}}_{\sigma }(\tau ,\mathbf{k}){]}_{pq}=-⟨T_{\tau }[c_{\mathbf{k}p\sigma }(\tau )c_{\mathbf{k}q\sigma }^{†}(0)]⟩=-\frac{1}{N}\sum _{i,j=1}^N{e}^{-i\mathbf{k}\cdot ({\mathbf{R}}_i-{\mathbf{R}}_j)}⟨T_{\tau }[c_{ip\sigma }(\tau )c_{jq\sigma }^{†}(0)]⟩$$

$$\tau >0\to [{\mathbf{G}}_{\sigma }^{>}(\tau ,\mathbf{k}){]}_{pq}=-\frac{1}{N}\sum _{i,j=1}^N{e}^{-i\mathbf{k}\cdot ({\mathbf{R}}_i-{\mathbf{R}}_j)}⟨c_{ip\sigma }(\tau )c_{jq\sigma }^{†}⟩$$

$$\tau <0\to [{\mathbf{G}}_{\sigma }^{<}(\tau ,\mathbf{k}){]}_{pq}=+\frac{1}{N}\sum _{i,j=1}^N{e}^{-i\mathbf{k}\cdot ({\mathbf{R}}_i-{\mathbf{R}}_j)}⟨c_{jq\sigma }^{†}(-\tau )c_{ip\sigma }⟩$$

$$\displaystyle\mathbf{G}_{\sigma}(i\omega=0,\mathbf{k})=\int_{-\infty}^{+\infty}\mathbf{G}_{\sigma}(\tau,\mathbf{k})d\tau=\int_{0}^{+\infty}\Big{[}\mathbf{G}_{\sigma}^{>}(\tau,\mathbf{k})+\mathbf{G}_{\sigma}^{<}(\tau,\mathbf{k})\Big{]}d\tau\simeq\int_{0}^{+\theta}\Big{[}\mathbf{G}_{\sigma}^{>}(\tau,\mathbf{k})+\mathbf{G}_{\sigma}^{<}(\tau,\mathbf{k})\Big{]}d\tau.$$

$$\mathcal{C}=$$

$$P(\mathbf{k})=\sum _{{\mu }_m>0}|{\varphi }_m(0,\mathbf{k})⟩⟨{\varphi }_m(0,\mathbf{k})|,$$

Full text

Introduction. Combination of the richness of quantum many-body effects and the elegance of topological physics Klitzing et al. (1980); Thouless et al. (1982); Haldane (1988); Hasan and Kane (2010); Qi and Zhang (2011) has revealed remarkable phenomena and new principles of physics, such as the fractional quantum Hall effect Tsui et al. (1982); Laughlin (1983) and topological order WEN (1990). Among these discoveries, one intriguing example is interaction-driven topological states, where strong correlations among particles convert a conventional state of matter into a topological one. One pathway towards such states is to utilize the phenomenon of spontaneous symmetry breaking Raghu et al. (2008); Sun et al. (2009); Herbut and Janssen (2014); Zhang et al. (2009); Yu et al. (2010), i.e. in a system where nontrivial topological structures are prohibited by symmetry, strong interactions can spontaneously break symmetry and thus stabilize a topologically nontrivial ground state. As proposed in Ref. Raghu et al. (2008), such a phenomenon can arise in a 2D Dirac semimetal (DSM) through a quantum phase transition that breaks spontaneously the time-reversal or the spin rotational symmetry, resulting in an interaction-driven, quantum-Hall or quantum-spin-Hall (QSH), topological insulator, dubbed topological Mott insulators (TMI).

Although the general principle about TMI has been well understood, finding such a state via unbiased theoretical/numerical methods turns out to be challenging due to the strong coupling nature of the problem and the presence of competing orders. Extensive numerical efforts on interacting DSMs García-Martínez et al. (2013); Daghofer and Hohenadler (2014); Guo and Jia (2014); Motruk et al. (2015); Capponi and Läuchli (2015); Scherer et al. (2015) report negative results, suggesting that in all explored parameter regimes, topologically-trivial competing states always have lower energy and thus the proposed TMI states cannot be stabilized. The successful alternative came lately, by substituting the DSM by a semimetal with a quadratic band crossing Sun et al. (2009); Herbut and Janssen (2014); Sun et al. (2012), an interaction-driven quantum Hall state is observed numerically Wu et al. (2016). Furthermore, experimental realization of such scenario has very recently been proposed in functionalized $\alpha$-Fe2O3 nanosheet Liang et al. (2017). However, whether a TMI can emerge from a DSM without the assistance of a quadratic band crossing point, as in the original proposal Raghu et al. (2008), still remains an open question. It is also worthwhile to highlight that between the two possible types of TMIs, quantum Hall and quantum spin Hall Raghu et al. (2008), only the former has been observed in numerical studies Wu et al. (2016). Hence, to find a time-reversal invariant TMI is one key objective of this study.

On the other hand, in a seemingly unrelated research area, recent developments in sign-problem-free quantum Monte Carlo (QMC) approaches for itinerant fermions coupled to fluctuating bosonic fields open the door to investigate many intriguing strongly-correlated systems, such as antiferromagnetic fluctuations mediated superconductivity in metals Berg et al. (2012); Schattner et al. (2016a), nematic quantum critical points in itinerant systems Lederer et al. (2015); Schattner et al. (2016b), as well as non-fermi liquid in itinerant quantum critical regions Wang et al. (2016a); Lederer et al. (2016); Xu et al. (2016). The strong-coupling nature of these systems makes analytical approach challenging Lee (2009); Metlitski and Sachdev (2010a, b); Dalidovich and Lee (2013); Schlief et al. (2016), and hence sign-problem-free QMC solutions pave a new avenue towards quantitative understanding about these systems. These QMC approaches also offer a new platform for studying strongly-correlated topological states, and have recently been utilized to study topological phase transitions in DSM Xu et al. (2017) and exotic states with topological order Assaad and Grover (2016); Sachdev et al. (2016); Gazit et al. (2017).

In this Letter, we study interaction-driven topological Mott insulators in Dirac semimetals with the aforementioned QMC approach. Instead of bare interactions, our model utilizes fluctuating bosonic fields to mediate interactions between fermions. At the level of the effective field theory, the model is equivalent to the originally proposed TMI model in Ref. Raghu et al. (2008), except for a minor difference in symmetry irrelevant to topology. For the study of TMI, our modified model shows two advantages: (1) other competing orders are strongly suppressed, allowing a clear TMI phase; (2) the sign-problem is avoided and thus the model can be solved via QMC techniques. Comparing to previous exact diagonalization studies García-Martínez et al. (2013); Daghofer and Hohenadler (2014); Guo and Jia (2014); Capponi and Läuchli (2015); Wu et al. (2016), the QMC approach can access much larger system size and reveals detailed information about the critical properties associated with the interaction-driven topological transition. Our QMC results show a continuous quantum phase transition from a DSM state to a QSH-type TMI phase, with the critical scaling at the quantum critical point agreeing nicely with the $N=8$ chiral Ising universality Chandrasekharan and Li (2013); Otsuka et al. (2016).

Model and Method. Our model describes Dirac fermions coupled to a transverse field Ising model. As illustrated in Fig. 1(a), fermions in this model reside on the lattice sites (disks), while Ising spins are placed on each dual lattice site (squares) at the plaquette centers. The Hamiltonian consists of three parts,

[TABLE]

where indices $i,j$ represent fermion sites and $p,q$ label the dual lattice sites for Ising spins $s^{z}$. Fermion spins are labeled by subindex $\sigma$. $H_{\text{Fermion}}$ describes the nearest-neighbor (NN) hopping for fermions, which contains a staggered flux $\pm 4\phi$ for each plaquette. Here, we request spin-up and spin-down fermions to carrier opposite flux patterns to preserve the time-reversal symmetry. The Ising spins are governed by $H_{\text{Ising}}$, which describes a ferromagnetic ($J>0$) transverse-field Ising model. The last term $H_{\text{Coupling}}$ couples the Ising spins with the next-nearest-neighbor (NNN) fermion hoppings, where the coupling constant $\xi_{ij}=\pm\xi t$ has a staggered sign structure alternating between neighboring plaquette, i.e., $+$ ($-$) for solid (dashed) NNN bonds as illustrated in Fig. 1(a). Up to a basis change, the low-energy physics in this model can be described by the following effective field theory $S=\sum_{\sigma}\int\mathbf{dr}dt\bar{\Psi}_{\sigma}(i\gamma^{\mu}\partial_{\mu}+g\sigma\varphi\gamma^{3}\gamma^{5})\Psi_{\sigma}+S_{\varphi}$, where $\gamma^{\mu}$ are gamma matrices and $\varphi$ is a bosonic field governed by the $\varphi^{4}$-theory $S_{\varphi}$. Here, $\sigma=\pm 1$ (up or down) is the fermion spin index, and $g$ is the coupling constant for the boson-fermion interactions. This effective field theory is in strong analogy to the model proposed early on in Ref. Raghu et al. (2008), provided that we decouple the fermion-fermion interactions with a Hubbard-Stratonovich auxiliary field, as appropriate in the limit $h/J\to\infty$ Assaad and Grover (2016). It is also worthwhile to emphasize that in our model, the fermion spins only preserve U(1) symmetry, while the model in Ref. Raghu et al. (2008) has a SU(2) spin symmetry. This difference has little effect on topological properties, but as discussed below it changes the critical scaling as well as the finite temperature phase diagram.

As in the original model of TMI, our Hamiltonian also contains a symmetry which prohibits nontrivial topology. It is easy to verify that our Hamiltonian is invariant under the following $Z_{2}$ transformation, $\hat{P}=\hat{R}_{x}(\pi)\times\hat{T}_{A\to B}$, where $\hat{R}_{x}(\pi)$ stands for $\pi$-rotation along $x$-axis for both Ising and fermion spins, and $\hat{T}_{A\to B}$ represents space translation from sublattice $A$ to $B$ inside a unit cell. Because the topological index (the spin Chern number) flips sign under this transformation, this symmetry requires the index to vanish and thus any (quantum spin Hall) topological insulator is prohibited, unless this $Z_{2}$ symmetry is broken spontaneously.

To explore the ground-state phase diagram of this model, we employ the projector quantum Monte Carlo (PQMC) method Assaad and Evertz (2008), with details presented in Sec.I.A of the supplemental material (SM) Sup . In addition to the usual local updates of Ising spins, both Wolff Wolff (1989) and geometric cluster updates Heringa and Blöte (1998) are applied in our simulations, as shown in Sec.I.B of SM Sup . Our QMC simulations are free of the sign problem at and away from half filling Wu and Zhang (2005). In this Letter, we focus on the coupling strength $0\leq\xi\leq 1$ with $J=t=1$ and the system sizes simulated in this work are $L=4,6,8,10,12,14$ with $N=L^{2}$ unit cells and $N_{s}=2L^{2}$ lattice sites.

Ground state phase diagram. The ground state phase diagram in the $\xi-h$ plane is shown in Fig. 1(c). Several regimes in the phase diagram can be solved exactly. At $\xi=0$, the fermions and Ising spins decouple: the fermions form a non-interacting Dirac semimetal, and the Ising spins undergo a paramagnetic to ferromagnetic (PM-FM) quantum phase transition at $h_{c}=3.046(3)$ in the 3D Ising universality class Pfeuty and Elliott (1971); Xu et al. (2017). At $h=0$, quantum fluctuations of Ising spins vanish and Ising spins form a fully-polarized FM state. As a result, the fermions turn into a non-interacting quantum-spin-Hall topological insulator, whose Hamiltonian is $H_{\text{Fermion}}+H_{\text{Coupling}}$ with fully polarized Ising spins $s^{z}=+1(-1)$  Sun et al. (2009); Hou (2013); Wu et al. (2016) (See Sec. V.A in the SM Sup for details). At $h\to\infty$, the Ising spins are aligned along the $x$-axis. Second order perturbation theory around this point, gives rise to an interaction of order $\xi^{2}/h$ between the fermions. Since the Dirac semimetal is a stable state of matter, we expect that it will be realized in the limit $h\rightarrow\infty$.

At $\xi>0$ and intermediate $h$, we find a direct second-order quantum transition between the PM and FM phases. This transition is also the topological phase transition for the fermions, in which the Dirac semimetal acquires a topological mass gap corresponding to the quantum spin Hall topological insulator. This conclusion is consistent with the symmetry analysis above, where the PM (FM) phase preserve (spontaneously breaks) the $Z_{2}$ symmetry and thus a quantum spin Hall insulator is prohibited (allowed). At $\xi>0$, the scaling exponents at the transition deviates from the 3D Ising universality class. Due to the coupling between fermions and bosons, the $\xi>0$ phase transition flows to a different universality class, namely the $N=8$ component chiral Ising universality class Chandrasekharan and Li (2013); Otsuka et al. (2016).

FM-PM phase transition for Ising spins. We determine the location of QCP via the Binder cumulant Binder (1981): $U_{2}=\frac{1}{2}(3-\frac{\langle m^{4}\rangle}{\langle m^{2}\rangle^{2}})$ and correlation ratio Kaul (2015); Pujari et al. (2016): $R_{\text{Corr}}=1-\frac{S^{\text{Ising}}(\mathbf{Q+\mathbf{q}})}{S^{\text{Ising}}(\mathbf{Q})}$, where $m=\frac{1}{N_{s}}\sum_{p}s_{p}^{z}$ and $S^{\text{Ising}}(\mathbf{k})$ is the trace of the structure factor matrix ($2\times 2$) of Ising magnetic order at $\mathbf{k}$ point. Here, $\mathbf{Q}=\boldsymbol{\Gamma}=(0,0)$ is the ordering vector for Ising spin, and $\mathbf{q}$ is the smallest momentum on the lattice, i.e., $(0,\frac{2\pi}{L})$ or $(\frac{2\pi}{L},0)$. Both $U_{2}$ and $R_{\text{Corr}}$ converge to 0 (1) in the PM (FM) phase at the thermodynamic limit. The crossing points for finite-size results of $U_{2}$ and $R_{\text{Corr}}$, respectively, provide the location of QCP. In this way, we first determine the position of QCP and then perform finite-size scaling analysis of $\langle m^{2}\rangle$ close to it to extract the critical exponents.

The results of $U_{2}$ and $R_{\text{Corr}}$, as well as the data collapse of $\langle m^{2}\rangle$ for $\xi=0.5$ and $\phi=\pi/4$ ($\pi$-flux in each plaquette) are presented in Fig. 2. Up to system size $L=12$, we can obtain the finite size crossing points $h=4.06$ for $U_{2}$ and $h=4.10$ for $R_{\text{Corr}}$ as the approximate location of QCP. In Fig. 2(c), we collapse the data as $\langle m^{2}\rangle L^{z+\eta}=f(L^{1/\nu}(h-h_{c})/h_{c})$ for $L=6,8,10,12,14$ and $L=10,12,14$, respectively. The critical exponents extracted from these two collapses are slightly different especially in $\eta$, indicating some finite-size effect. As will be discussed below, this shifting of exponents is due to a crossover phenomenon. Combining both collapses, we take the exponents as $\nu=0.85(2),\eta=0.61(7)$ (taking $z=1$) with $h_{c}=4.11(1)$, which are well consistent with the results presented in Ref. Chandrasekharan and Li (2013) as $\nu=0.83(1),\eta=0.62(1)$ for $N=8$ components chiral Ising universality class.

We employed two additional measurements to further corroborate the critical exponents. First, we performed finite-size scaling analysis for $S^{\text{Ising}}(\mathbf{k})$ at $\xi=0.50$ and $\phi=\pi/4$, which is shown in Sec. II.B of the SM Sup , with the extracted critical exponents $\nu=0.84(4),\eta=0.62(6)$. Second, we also simulated the model with $\xi=0.50$ and $\phi=\pi/8$ (half-$\pi$ flux) and obtained the critical exponents from the finite-size scaling of $\langle m^{2}\rangle$, and the results are presented in Sec. III.A of SM Sup . The obtained critical exponents are $\nu=0.85(3),\eta=0.63(7)$ with $h_{c}=4.242(3)$. These exponents are well consistent with those in Fig. 2(c), rendering the $N=8$ components chiral Ising universality class.

The properties of QCPs for the PM-FM phase transitions of Ising spins for $\xi=0.25,0.75,1.00$ as presented in the phase diagram of Fig. 1(c), are also determined with $U_{2}$ and $R_{\text{Corr}}$, as well as the finite-size scaling of $\langle m^{2}\rangle$ and excitation gaps of fermions.

Topological phase transition for fermions. Our numerical results further show a single phase transition from Dirac semimetal to topological Mott insulator with decreasing transverse-field $h$, which comes hand in hand with the PM-FM phase transition of Ising spins. As shown in Fig. 3, we find that the fermions remains gapless in the PM phase with vanishing gap at the Dirac point. Here, both $\Delta_{sp}(\mathbf{X})$, the average single particle gap at two Dirac points $\mathbf{X}_{1}$ and $\mathbf{X}_{2}$, and $\Delta_{s}(\mathbf{M})$, the two-particle spin gap at $M$, vanish at the thermodynamic limit, consistent with the Dirac semimetal spectrum. In the FM phase, both gaps start to merge at $h/J=4.10\sim 4.15$, consistent with the location of PM-FM phase transition point for Ising spins. It is worthwhile to highlight that the gaps remain finite in the whole FM phase with $h<h_{c}$, indicating the absence of topological phase transition. Since the fermions form a quantum spin Hall insulator in the exactly solvable limit at $h=0$, this finite gap implies that the whole FM phase shares the same nontrivial topology. To further verify this conclusion, we compute directly the topological invariant, the spin Chern number $C_{s}=(C_{\uparrow}-C_{\downarrow})/2$. As shown in Sec. V.B. of SM Sup , we obtain $C_{s}=+1$ for whole $h<h_{c}$ region, indicating that FM phase is a quantum-spin-Hall topological insulator.

In Sec. IV of SM Sup , we present raw data of dynamic quantities $G(\mathbf{k},\tau)$ and $S^{xy}(\mathbf{k},\tau)$, from which $\Delta_{sp}(\mathbf{X})$ and $\Delta_{s}(\mathbf{M})$ are extrapolated. The comparisons between $2\Delta_{sp}(\mathbf{X})$ and $\Delta_{s}(\mathbf{M})$ are also shown to reveal the effect of electron-electron interactions. Furthermore, the gap opening of $\Delta_{sp}(\mathbf{X})$ and $\Delta_{s}(\mathbf{M})$ at $\xi=0.25,0.75,1.00$ match the QCPs of PM-FM phase transition for Ising spins, thus supporting the picture of a semimetal-TMI topological phase transition.

Finite-size scaling crossover. As discussed above at $\xi=0$ and $\xi>0$, the PM-FM transition belongs to two different universality classes, 3D Ising and $N=8$ chiral Ising. As a result, in the thermodynamic limit, the scaling exponents will change discontinuously as we change the value of $\xi$ away from [math]. In numerical studies, because of the finite size, such a discontinuous change will not show up. Instead, a crossover behavior is expected, i.e. at small $\xi$, a crossover length scale $L_{c}(\xi)$ shall arise. For $L<L_{c}$ ($L>L_{c}$), the scaling behavior merges towards the 3D Ising ($N=8$ chiral Ising) universality class. As $\xi$ approaches zero (increases), $L_{c}$ diverge to infinity (decreases to microscopic values) and thus the 3D Ising ($N=8$ chiral Ising) universality class is fully recovered. Such an effect is indeed observed in our data. In Sec. VI in SM Sup , we present the finite-size scalings of $\langle m^{2}\rangle$ from $L=6,8,10,12$ and $L=8,10,12$, respectively, for $\xi=0.25,0.50,0.75$. At $\xi=0.25$, the data collapse suffers strongly from the finite-size effect, and chiral Ising exponents only arise in very large system sizes, especially for $\eta$. However, as $\xi$ increases, the chiral Ising exponents emerge even if the smallest size $L=6$ is included in the fitting.

Discussions. Because the fermion spin in our model only preserves a $U(1)$ symmetry, instead of $SU(2)$, our topological Mott insulator breaks a $Z_{2}$ symmetry in contrast to the $SU(2)$ symmetry breaking in Ref. Raghu et al. (2008). This difference in symmetry breaking patterns is irrelevant for topology. However, this leads to different scaling exponents at the transition Chandrasekharan and Li (2013). Furthermore, at finite temperature, the symmetry breaking phase in our model survives, while the $SU(2)$ symmetry breaking arises only at $T=0$.

To the best of our knowledge, our study demonstrates the first interaction-driven quantum-spin-Hall topological Mott insulator from unbiased numerical method, and for the first time, this novel topological phenomenon becomes accessible to large-scale lattice QMC simulations. Our work points out a new route to realize interaction-driven topological phases and phase transitions. It has experimental relevance since the interaction-driven quantum anomalous Hall effect has recently being suggested in functionalized $\alpha$-Fe2O3 nanosheet Liang et al. (2017).

We (YYH, XYX, ZYM and ZYL) acknowledge fundings from the Ministry of Science and Technology of China through National Key Research and Development Program under Grant No. 2016YFA0300502 and from the National Science Foundation of China under Grant Nos. 91421304, 11421092, 11474356, 11574359, 11674370 as well as the National Thousand-Young Talents Program of China. Y.Y.H is also supported by the Outstanding Innovative Talents Cultivation Funded Programs 2016 of Renmin University of China. K.S. acknowledges support from the National Science Foundation under Grant No. PHY1402971 and the Alfred P. Sloan Foundation. F.F.A thanks the German Research Foundation (DFG) for financial support through the SFB 1170 ToCoTronics. We thank the Physical Laboratory of High Performance Computing in Renmin University of China, the Center for Quantum Simulation Sciences in the Institute of Physics, Chinese Academy of Sciences and the Tianhe-1A platform at the National Supercomputer Center in Tianjin for their technical support and generous allocation of CPU time.

I I. Projector QMC Algorithm for the model

I.1 A. PQMC Formalism

The fermion-Ising model in the main text is expressed as

[TABLE]

In the projector QMC setup, the partition function of this model can be expressed as

[TABLE]

where the trace is taken over all Ising spin configurations, and $|\Psi_{T}\rangle$ is the Slater-determinant trial wave function.

In projector QMC framework, we first break the projection length $\Theta$ into $M$ slices ($\Theta=M\Delta\tau$). Since the $s_{p}^{z}$ operator in $H_{\text{Coupling}}$ term is diagonal under Ising spin configurations, the Boltzmann weight from $H_{\text{Ising}}$ is separable. After tracing out fermion degrees of freedom, we can obtain the partition function of the fermion-Ising coupling model as Xu et al. (2017)

[TABLE]

where

[TABLE]

is from $H_{\text{Ising}}$, when the 2D transverse-field Ising model is mapped to a 3D classical Ising model Isakov and Moessner (2003); Blöte and Deng (2002); Wang et al. (2016b); Xu et al. (2017), with

[TABLE]

Due to the spin-staggered phase $e^{i\sigma\phi}$ in $H_{\text{Fermion}}$ term in Eq. (S1), the spin-up determinant $\det(P_{\uparrow}^{\dagger}B_{M}^{\uparrow}\cdots B_{2}^{\uparrow}B_{1}^{\uparrow}P_{\uparrow})$ in Eq. (I.1) is complex-conjugate to the spin-down one $\det(P_{\downarrow}^{\dagger}B_{M}^{\downarrow}\cdots B_{2}^{\downarrow}B_{1}^{\downarrow}P_{\downarrow})$. Since the configuration weight for the classical Ising spin part is always positive, the configuration weight containing both the fermion and classical spin parts is non-negative thus there is no sign problem in the QMC simulations of the model Hamiltonian in Eq. (S1). To calculate the $B_{\ell}^{\sigma}$, we have applied Trotter decomposition for both the $e^{-\Delta\tau H_{\text{Fermion}}}$ and $e^{-\Delta\tau H_{\text{Coupling}}}$ terms. For the $H_{\text{Fermion}}$ term, we have applied the checkerboard decomposition, which divides all the NN hopping terms into two parts such that in each part all the hopping terms commute with one another. For the $H_{\text{Coupling}}$ term as the fermion-Ising spin coupling term, it’s separated into four parts such that in each part all the hopping term commute with one another. Including both of these decompositions, we have the $B_{\ell}^{\sigma}$ matrix as

[TABLE]

where $\mathcal{C}_{\ell,4},\mathcal{C}_{\ell,3},\mathcal{C}_{\ell,2},\mathcal{C}_{\ell,1}$ correspond to the $H_{\text{Coupling}}$ term while $\mathcal{K}_{2},\mathcal{K}_{1}$ represent the $H_{\text{Fermion}}$ term. Here, all of $\mathcal{C}_{\ell,4},\mathcal{C}_{\ell,3},\mathcal{C}_{\ell,2},\mathcal{C}_{\ell,1}$ and $\mathcal{K}_{2},\mathcal{K}_{1}$ are effectively $4\times 4$ matrices.

For the results presented in main text, we set $\Delta\tau=0.04/t$ and purposely increase the projection parameter $\Theta$ for increasing system size, such as $\Theta=20$ for $L=4$ and $\Theta=50$ for $L=14$, to ensure convergence to ground state.

I.2 B. Cluster Updates for Ising spins

As for the updates of the classical Ising spins during the PQMC simulations, in one sweep, we first apply the local updates to scan through the space-time of configuration and then take one cluster update. For the case studied here, the 3D classical Ising model (2D transverse field Ising) has ferromagnetic interactions in both space and time directions as

[TABLE]

with $J>0$ and $\gamma=-\frac{1}{2}\ln[\tanh(\Delta\tau h)]>0$. This is a highly anisotropic model and the ratio between the coupling strengths in time and space is $r=\gamma/(\Delta\tau J)$, which can reach $10^{1}\sim 10^{3}$. To guarantee the efficiency and correctness, we employ local updates, Wolff cluster updates Wolff (1989) and geometric cluster update Heringa and Blöte (1998) in the simulation. The local update is, of course, the simplest one to be implemented. The two cluster updates, greatly suppress the critical slowing down, are complementary. As close to the quantum critical point, the size of Wolff cluster can sometimes be as large as the entire space-time lattice, rendering the update actually with low efficiency around the quantum criticality. The geometric cluster updates implicitly impose a restriction on the size of the cluster, making it not to grow too large. On the other hand, the geometric cluster updates can not change the magnetization of the system, but Wolff cluster update can explicitly change the value of $\sum_{p}s_{p,\ell}^{z}$.

In PQMC, a sweep contains the propagation of the configurations in the imaginary time from $\tau=\Theta$ to $\tau=0$ and successive propagation from $\tau=0$ to $\tau=\Theta$. In the propagations of both $\tau=\Theta\to 0$ and $\tau=0\to\Theta$, the local updates of classical Ising spin are carried out sequentially according to the space-time lattice. Then at $\tau=0$, the cluster updates are performed to flip spins in a cluster constructed by Wolff and geometric algorithms.

The detailed balance condition of Wolff or geometric updates is given by

[TABLE]

where $a,b$ are configurations of Ising spins, and $P(a\to b)$, $\mathcal{A}(a\to b)$ are the priori and acceptance probabilities. With $\pi(a)=\pi_{s}(a)\pi_{f}(a)$ (where $\pi_{s}(a)$ is the weight of spin part and $\pi_{f}(a)$ is the weight of fermion part) as the total configuration weight, both Wolff and geometric cluster algorithms carefully choose the cluster such that the priori probability ratio has the following property

[TABLE]

Then we can obtain the acceptance probability as

[TABLE]

One can see that the acceptance probability is simply the ratio of fermion part of weights for the classical Ising spin configurations. At $\tau=0$, arbitrary number of steps of Wolff and geometric cluster updates can be performed and the overall acceptance probability is

[TABLE]

Hence the implementation of cluster updates in PQMC simulations goes as follows, first construct the cluster of Ising spins according to Wolff and geometric cluster algorithms, and then calculate the ratio of fermion part of weights to obtain the acceptance probability, such that the cluster update can be decided to be accepted or rejected.

From our simulations, we find that, in one sweep, performing one Wolff cluster update and one geometric cluster update reaches the highest acceptance ratio for the overall cluster updates. For example, the acceptance ratio for cluster updates is typically four times larger than that for the local updates. Table 1 lists the number of spins flipped, size of the Wolff and geometric clusters, for linear system size $L=12$ with $\xi/t=0.50,h/J=3.35$ close to FM-PM phase transition. One can see that in one step of cluster update, thousands of Ising spins are flipped.

II II. Quantum Critical Points of Ising spins

In this section, we demonstrate some additional results for the FM-PM phase transition of the Ising spins. These include Sec. II. A., the determination of the QCPs for $\xi=0.25,0.75,1.00$, via the Binder cumulant $U_{2}$ and correlation ratio $R_{\text{Corr}}$. Sec. II. B, data collapse of structure factor $S^{\text{Ising}}(\boldsymbol{\Gamma})/L^{2}$, to extract the critical exponents.

II.1 A. QCPs for different $\xi$

As shown in the main text, we first determine the FM-PM phase transition points from the crossing of the Binder cumulant $U_{2}$ and correlation ratio $R_{\text{Corr}}$, defined as

[TABLE]

with the Ising spin magnetization $m=\frac{1}{N_{s}}\sum_{p}s_{p}^{z}$ (here we have $N_{s}=2L^{2}$). And the structure factor of classical Ising spin is defined on checkerboard lattice (two sublattices) as

[TABLE]

with $m,n$ as unit cell indexes and $a=A,B$ indicating the sublattice ($N=L^{2}$ as number of unit cells). In Eq. (S14), $\mathbf{Q}$ is the ordered wave vector, which is $(0,0)$ for the ferromagnetic order, and $\mathbf{q}$ is the smallest distance vector away from $\mathbf{Q}$, equal to $(0,2\pi/L)$ or $(2\pi/L,0)$. It is expected that both $U_{2}$ and $R_{Corr}$ converge to 1 inside the ordered phase in the thermodynamic limit, while they are both 0 in disordered phase. So there is a jump from 1 to 0 when going from ordered phase to the disordered phase in thermodynamic limit, and for finite-size system they both show crossing between data with different system sizes, and the location of the crossing point defines the quantum critical point. Also, both $\langle m^{2}\rangle$ and $S^{\text{Ising}}(\boldsymbol{\Gamma})/L^{2}$ can tell whether the FM-PM phase transition is continuous or of first-order, depending on whether they are smooth or have drops close to the phase transition point.

In Fig. S1,  S2 and  S3, we present the results of both $U_{2}$ and $R_{\text{Corr}}$ for $\xi=0.25,0.75,1.00$ across the QCPs with $L=4,6,8,10,12$. Similar to the results of $\xi=0.50$, shown in the main text, both $U_{2}$ and $R_{\text{Corr}}$ converges to $1$ deeply inside the ordered phase and [math] inside the disordered phase. The crossing points, for example $h/J=3.368,4.103,5.210,6.78$ for $L=10,12$ of $R_{\text{Corr}}$ for $\xi=0.25,0.50,0.75,1.00$, are good estimation of the QCPs, from which, we can further determine the more precise position of the thermodynamic QCPs by combining the data collapse of $\langle m^{2}\rangle$ and extrapolation of excitation gaps.

II.2 B. Data Collapse for Structure factor $S^{\text{Ising}}(\boldsymbol{\Gamma})/L^{2}$

In the main text, we have obtained the critical exponents of the FM-PM phase transition of Ising spins from the data collapse of $\langle m^{2}\rangle$ data with the relation $\langle m^{2}\rangle L^{z+\eta}=f(\frac{h-h_{c}}{h_{c}}L^{1/\nu})$. In this part, we want to show that we can obtain the same results for $\nu,\eta$ exponents from the data collapse of $S^{\text{Ising}}(\boldsymbol{\Gamma})/L^{2}$ with similar relation as $L^{z+\eta}\cdot[S^{\text{Ising}}(\boldsymbol{\Gamma})/L^{2}]=f(\frac{h-h_{c}}{h_{c}}L^{1/\nu})$.

According to the definition of $S^{\text{Ising}}(\mathbf{k})$ in Eq. (S15), we have

[TABLE]

with $m_{\alpha}=\frac{1}{N}\sum_{p\in\alpha}s_{p}^{z}$ ($\alpha=A,B$) sublattice magnetization for Ising spins. Similarly, we can express the $\langle m^{2}\rangle$ by $m_{A}$ and $m_{B}$ as

[TABLE]

Then we can observe that at $h\to 0$, the system is almost classically ordered and we have $\langle m_{A}^{2}\rangle=\langle m_{B}^{2}\rangle=\langle m_{A}m_{B}\rangle$ due to the classical decoupling of the correlation as $\langle m_{A}m_{B}\rangle=\langle m_{A}\rangle\langle m_{B}\rangle$. Thus, in the $h\to 0$ limit, we have $\langle m^{2}\rangle=S^{\text{Ising}}(\boldsymbol{\Gamma})/L^{2}$. However, for finite $h$, the equality doesn’t hold anymore and $\langle m_{A}^{2}\rangle$ can be different from $\langle m_{A}m_{B}\rangle$ across the QCPs. However, they represent the same physical meaning of magnetization. Thus, we can surely perform the data collapse for them independently and the critical exponents extracted from them are expected to be the same.

Below we perform the data collapse for $S^{\text{Ising}}(\boldsymbol{\Gamma})/L^{2}$ for $\xi=0.50$ at $\pi$-flux case ($\phi=\pi/4$). The results are presented in Fig. S4. The collapse of data from $L=6,8,10,12,14$ yields $\nu=0.81(1),\eta=0.57(2),h_{c}/J=4.084(5)$, while the collapse of data from $L=10,12,14$ gives $\nu=0.85(2),\eta=0.67(1),h_{c}/J=4.108(3)$ . Combining them, we conclude $\nu=0.84(4),\eta=0.62(6)$ for the FM-PM phase transition, which is well consistent with the results extracted from data collapse of $\langle m^{2}\rangle$, as shown in Fig.2 in main text.

III III. Quantum Phase Transitions for half-$\pi$ flux case ($\phi=\pi/8$)

To further confirm the general properties of both Ising spins and fermions across the QCPs, we have also simulated the half-$\pi$-flux case as $\phi=\pi/8$. In Sec. III. A., we show that the FM-PM phase transition has almost the same critical exponents as that for the $\phi=\pi/4$ case, indicating the same universality class. What’s more, the excitation gaps of the fermions also show gap opening behavior with decreasing $h$, suggesting the same DSM-TMI phase transition as that for $\phi=\pi/4$ case, which are presented in Sec. III. B. For this half-$\pi$ flux case ($\phi=\pi/8$), we only concentrate on $\xi=0.50$.

III.1 A. FM-PM phase transition for Ising spins

Again, the location of the QCP is determined from the crossings of $U_{2}$ and $R_{\text{Corr}}$ for different system sizes, which are shown in Fig. S5. Considering the convergence to thermodynamic limit, the QCP of this $\xi=0.50$ with $\phi=\pi/8$ should be $h_{c}/J>4.262$ (see the insets of Fig. S5), which is larger than that for $\phi=\pi/4$ case with the same $\xi$ parameter.

Then we take a closer look at the QCP and perform the data collapse of $\langle m^{2}\rangle$ around the QCP. The results are shown in Fig. S6. Again, we have performed the collapse of the data for $L=6,8,10,12$ and $L=8,10,12$ to obtain the reasonable values of $\nu$ and $\eta$. The collapse of data from $L=6,8,10,12$ yields $\nu=0.84(2),\eta=0.58(2),h_{c}/J=4.242(3)$, while the collapse of data from $L=8,10,12$ gives $\nu=0.86(2),\eta=0.67(2),h_{c}/J=4.272(6)$. Combining them, we conclude $\nu=0.85(3),\eta=0.63(7)$ for the FM-PM phase transition, which is well consistent with the results extracted from data collapse of $\phi=\pi/4$ case, as shown in Fig.2 in main text as well as Fig. S4 in Sec. II. B. These results imply that both the FM-PM phase transitions of Ising spins for $\phi=\pi/4$ and $\phi=\pi/8$ cases belong to the $N=8$ Chiral Ising universality class, with critical exponents consistent with previous work in Ref. Chandrasekharan and Li, 2013; Otsuka et al., 2016.

III.2 B. DSM-TMI phase transition for fermions

Besides the PM-FM phase transition for the Ising spins, we have also confirmed it’s accompanied by the DSM-TMI phase transition, via monitoring the opening of excitation gaps. We have calculated the single-particle gap $\Delta_{sp}(\mathbf{X})$ (average over $\mathbf{X}_{1}$ and $\mathbf{X}_{2}$ points) and spin gap $\Delta_{s}(\mathbf{M})$.

The extrapolations of $\Delta_{sp}(\mathbf{X})/t$ and $\Delta_{s}(\mathbf{M})/t$ over $1/L$ are shown in Fig. S7. We can observe that with decreasing $h/J$, both $\Delta_{sp}(\mathbf{X})/t$ and $\Delta_{s}(\mathbf{M})/t$ have gap opening at $h/J\in[4.25,4.30]$, suggesting the DSM-TMI phase transition at $h_{c}/J\in[4.25,4.30]$. This location of the QCP is consistent with the results of Binder cumulant $U_{2}$ and Correlation Ratio $R_{\text{Corr}}$ shown in Fig. S5, and also the QCP from the data collapse presented in Fig. S6. These consistency all suggest that the PM-FM phase transition for Ising spins and DSM-TMI phase transition for fermions happens simultaneously, and in the $N=8$ Chiral Ising universality class.

IV IV. Raw Data of dynamic properties and Excitation Gaps

In this section, we present raw data on dynamic properties, including the dynamic single-particle Green’s function $G(\mathbf{X},\tau)$ and dynamic spin-spin correlation function $S^{xy}(\mathbf{M},\tau)$ in Sec. IV. A. and the comparisons of single-particle gap $\Delta_{sp}(\mathbf{X})/t$ and spin gap $\Delta_{s}(\mathbf{M})/t$ in Sec. IV. B.

IV.1 A. Dynamic single-particle Green’s function and spin-spin correlation function

Since we extract the excitation gaps $\Delta_{sp}(\mathbf{X})/t$ and $\Delta_{s}(\mathbf{M})/t$ from $G(\mathbf{X},\tau)\propto e^{-\Delta_{sp}(\mathbf{X})\tau}$ and $S^{xy}(\mathbf{M},\tau)\propto e^{-\Delta_{s}(\mathbf{M})\tau}$ at large $\tau$ limit, we need to ensure that the data of $G(\mathbf{X},\tau)$ and $S^{xy}(\mathbf{M},\tau)$ have high quality. Here, we want to demonstrate that the dynamic data obtained from our simulations indeed have very good quality.

The raw data of $G(\mathbf{X},\tau)$ and $S^{xy}(\mathbf{M},\tau)$ for $\xi=0.50,\phi=\pi/4$ and $L=10$ are shown in Fig. S8 with both linear and semi-log coordinates. The perfect straight lines of $\ln[G(\mathbf{X},\tau)]$ and $\ln[S^{xy}(\mathbf{M},\tau)]$ over $\tau$ allow us to extract the gaps $\Delta_{sp}(\mathbf{X})/t$ and $\Delta_{s}(\mathbf{M})/t$ with high precision. This good quality of the dynamic data also holds for different system size $L$, different $\xi$ parameter and also different $\phi$ such as $\phi=\pi/8$.

Notice that here we have calculated the single-particle gap at $\mathbf{X}_{1}$ and $\mathbf{X}_{2}$ points in BZ and the spin gap at $\mathbf{M}$ point. Since in the large $h$ limit, the fermions has DSM ground state and the Dirac cones locate at $\mathbf{X}_{1}$ and $\mathbf{X}_{2}$ points. With decreasing $h$, the fermions open gaps at $\mathbf{X}_{1}$ and $\mathbf{X}_{2}$ points and enter into the TMI phase. As for the spin gap, it has smallest value at $\mathbf{M}$ point, which can be understood as follows: since the coupling term $\hat{H}_{\text{Coupling}}$ in the model Hamiltonian of Eq. (S1) couples the Ising spin with next-nearest-neighbor (NNN) hopping of fermions, it effectively induces a NNN density-density interaction of fermions when integrate out the Ising spins, and that NNN density-density interaction favors collinear order on square lattice, which has ordered wave vector $\mathbf{Q}=\mathbf{M}$ on checkerboard lattice. We note, that this $\mathbf{Q}=\mathbf{M}$ spin order of the fermions should not be confused with the usual antiferromagnetic long-range order on the square lattice.

IV.2 B. Excitation Gaps

In Fig. 3 in the main text, the extrapolations of $\Delta_{sp}(\mathbf{X})/t$ and $\Delta_{s}(\mathbf{M})/t$ over $1/L$ are presented, and one can observe that the values of $2\times\Delta_{sp}(\mathbf{X})/t$ seem to be very close to those of $\Delta_{s}(\mathbf{M})/t$ within the same model parameters and system size. It’s well known that the two-particle excitation gaps are twice of single-particle gap in noninteracting fermion systems. What’s more, we have indeed observed $\Delta_{s}(\mathbf{M})/t\approx 2\times\Delta_{sp}(\mathbf{X})/t$ at small $h$, since it’s very close to noninteracting fermion system. However, close to the QCP in the phase diagram, this relation does not hold any longer.

In Fig. S9, the comparison of twice of single-particle gap $2\times\Delta_{sp}(\mathbf{X})/t$ and spin gap $\Delta_{s}(\mathbf{M})/t$ for $\xi=0.50$ under $\pi$-flux case ($\phi=\pi/4$), for systems $L=4,6,8,10,12$ across the DSM-TMI phase transition are shown. One can see that $2\times\Delta_{sp}(\mathbf{X})/t$ is larger than $\Delta_{s}(\mathbf{M})/t$. This indicates the presence of effective electron-electron interactions, mediated by the Ising spin flucutations in the model Hamiltonian. The observation that $2\times\Delta_{sp}(\mathbf{X})/t$ is larger than $\Delta_{s}(\mathbf{M})/t$ also holds for different $\xi$ parameter and $\phi=\pi/8$ cases close to the DSM-TMI transition.

In Fig. 2 in the main text, we have shown the extrapolations of excitation gaps $\Delta_{sp}(\mathbf{X})/t$ and $\Delta_{s}(\mathbf{M})/t$ over $1/L$ as the proof of DSM-TMI phase transition for fermions only for $\xi=0.50$ with $\phi=\pi/4$. Here, we also present the extrapolations of $\Delta_{sp}(\mathbf{X})/t$ for $\xi=0.25,0.75,1.00$ with $\phi=\pi/4$ and show the gap opening during the DSM-TMI phase transition. The results are shown in Fig. S10. We only present the data of $\Delta_{sp}(\mathbf{X})/t$ while the spin gap $\Delta_{s}(\mathbf{M})/t$ has almost the same behaviour of gap opening with decreasing $h$. For $\xi=0.25,0.75,1.00$ cases with $\phi=\pi/4$, the excitation gaps opens at $h_{c}/J\in[3.35,3.40]$, $h_{c}/J\in[5.20,5.25]$ and $h_{c}/J\in[6.75,6.80]$, respectively, which are well consistent with data crossing points of Binder cumulant and correlation ratio for Ising spins in Fig. S1, Fig. S2 and Fig. S3.

V V. Topological nature of TMI

In addition to the opening of excitation gaps, in this section we demonstrate the fermions in TMI indeed have QSHI ground state at $h<h_{c}$. At $h=0$ limit, Ising spins are classically ordered and the fermions become non-interacting, here one can analytically show the system has QSHI ground state, this part is presented in Sec. V. A. On the other hand, at finite $h<h_{c}$, we can also provide theoretical arguments and more importantly numerical evidence that the system is in QSHI ground state as well, demonstrated in Sec. V. B.

V.1 A. $h=0$ Limit

At $h=0$, the Ising spin in Eq. (S1) ordered in a classical way – ferromagnetic order without quantum fluctuation. Thus, Hamiltonian becomes non-interacting and remaining fermion part can be written as

[TABLE]

With the choices of $t_{ij}^{\prime}=\pm t_{2}$, this simple tight-binding model has QSHI ground state Sun et al. (2009); Hou (2013); Wu et al. (2016). This model has $U(1)_{\text{spin}}\times U(1)_{\text{charge}}\rtimes Z_{2}^{T}$ symmetry (here $Z_{2}^{T}$ stands for time-reversal symmetry), which renders the $Z$ classification of the topological index. In Fig. S11, we present the spectral functions for spin-up and spin-down part along one of the edges for the model in Eq. (S18) on a ribbon geometry (periodic boundary condition in $x$-direction, open boundary condition in $y$-direction) shown in Fig. S11(a). We can observe that it has helical edge states with a large band gap in the bulk. We note here the topological invariant, spin Chern Number $C_{s}=(C_{\uparrow}-C_{\downarrow})/2=+1$, can also be calculated via the standard zero-frequency Green’s function formalism He et al. (2016). Hence, combining these information, one can see that at $h=0$ the system is indeed in a QSHI ground state for any finite $\xi$.

V.2 B. Finite $h$ with $h<h_{c}$

As shown both in the main text as well as in Sec. III. B and Sec. VI. A, the fermions in the coupling model Eq. (S1) has DSM ground state for $h>h_{c}$, namely the Ising spins in paramagnetic state. Decreasing $h$ to $h<h_{c}$, we have observed the opening of single-particle and two-particle excitation gaps, the system enters TMI phase with QSHI character. But since here fermions are interacting, mediated by quantum fluctuations of the Ising spins, we cannot obtain the band structure analytically as in Sec. V. A., hence will relie on the QMC results.

To numerically verify the QSHI ground state of fermions in $h<h_{c}$ with finite $\xi$, we calculate the topological invariant: the spin Chern number $C_{s}=(C_{\uparrow}-C_{\downarrow})/2$ with $C_{\sigma}$ as the Chern number for spin-$\sigma$ channel. Due to the time-reversal symmetry, $C_{\uparrow}=-C_{\downarrow}$, thus $C_{s}=C_{\uparrow}$. We have applied the method with zero-frequency single-particle Green’s function to calculate $C_{s}$, which was successfully demonstrated by some of us in Ref. He et al., 2016. As mentioned in the main text, since the coupling model Eq. 1 has $Z_{2}$ symmetry, we need to add a pinning field $\hat{H_{z}}=B_{z}\sum_{p}s_{p}^{z}$ to Ising spins to break this $Z_{2}$ symmetry, incorporating the effect of spontaneously $Z_{2}$ symmetry breaking in thermodynamic limit of TMI phase. In practical simulations, we choose $B_{z}=0.001J$. As for the calculation of $C_{s}$, we first obtain the $\mathbf{G}_{\sigma}(\tau,\mathbf{k})$ data with both $\tau\geq 0$ and $\tau<0$ from the QMC simulations as

[TABLE]

where $p,q=A,B$ standing for sublattice indexes. Then we can obtain the zero-frequency single-particle Green’s function $\mathbf{G}_{\sigma}(i\omega=0,\mathbf{k})$ as

[TABLE]

In the last step of Eq. (S20), a cut-off is applied to the integral, due to exponential decaying of $\mathbf{G}_{\sigma}(\tau,\mathbf{k})$ at large $\tau$. Then we can use the following formula with $\mathbf{G}_{\sigma}(i\omega=0,\mathbf{k})$ to calculate $C_{\sigma}$ for finite-size system as

[TABLE]

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

[TABLE]

and $\mu_{m}$ is the corresponding eigenvalue of $\mathbf{G}_{\sigma}(0,\mathbf{k})$ with eigenvector $|\phi_{m}(0,\mathbf{k})\rangle$.

Due to the finite momentum mesh in $L\times L$ system, the spin Chern number $C_{s}$ calculated from Eq. (S21) suffers finite-size effect and can be away from quantized integer for small $L$. Thus, we apply the interpolation scheme of $\mathbf{G}_{\sigma}(i\omega=0,\mathbf{k})$ to achieve denser momentum resolution and approach quantized spin Chern number, as demonstrated in Ref. He et al., 2016.

The results of spin Chern number $C_{s}$ for $L=4,8,12$ systems are shown in Fig. S12(a) for $\xi=0.50$ with $\phi=\pi/4$. Since the spin Chern number in Eq. (S21) is not well defined for gapless DSM phase, we only measure it in the QSHI phase with $h/J\in[0,4]$, i.e., $h<h_{c}$. As one can see, with increasing linear system size, $L=8$ to $L=12$, $C_{\uparrow}(=C_{s})$ increases gradually though it’s not quantized. And since $h/J=4$ is close to the QCP at $h_{c}=4.11(1)$ and the single-particle gap $\Delta_{sp}(\mathbf{X})/t$ is small, the cut-off $\theta$ causes the dip in $C_{\uparrow}$ though the system is still in the QSHI phase. After the interpolation with $IL=512$ from $L=8,12$ systems, $C_{\uparrow}=1$, reaches quantized integer perfectly, meaning that the system is indeed inside TMI phase with QSHI character.

We have also measured the excitation gaps for $\xi=0.50,\phi=\pi/4$ with decreasing $h$ from $h/J=4$ to $h/J=0$, as shown in Fig. S12(b), in which we have observed that both single-particle gap $\Delta_{sp}(\mathbf{X})/t$ and spin gap $\Delta_{s}(\mathbf{M})/t$ increase monotonously to the analytical values $\Delta_{sp}(\mathbf{X})/t=2$ and $\Delta_{s}(\mathbf{M})/t=4$ at $h=0$ point. This demonstrates the fact that there is no further topological phase transition in $h<h_{c}$ region, the QSHI at $h=0$ smoothly cross over into TMI at finite $h$. Furthermore, we have scanned a path inside the gapped region at $h/J=3$ with $\xi\in[0,1]$ and the results for energy and structure factors for several fermion bilinears (not shown) have no signature of phase transition, suggesting that the fermions in the whole $h<h_{c}$ region are in TMI ground state.

VI VI. Finite-size scaling crossover of FM-PM phase transition for Ising spins

Without coupling to fermions, 2D transverse-field Ising model has a continuous $h$-tuned quantum phase transition at 3D Ising universality class. Thus, at $\xi=0$ point of model in Eq. (1), the $h$-tuned quantum phase transition belongs to the 3D Ising universality class with $\nu=0.629971(4),\eta=0.026298(2)$ Hasenbusch et al. (1999). With coupling to fermion, however, the universality class is altered to the $N=8$ Chiral Ising. Therefore, there are two different kinds of universality classes along the phase transition line (the red line) in the ground-state phase diagram in Fig. 1 (c) in main text. A natural question here is, how the 3D Ising universality class evolves into the $N=8$ Chiral Ising universality class. In this section, we present some results on this problem.

We think that there exists a finite-size scaling crossover behavior of the universality class for the QCPs of Ising spins. In the thermodynamic limit, the QCP should belong to 3D Ising universality at $\xi=0$ point and $N=8$ Chiral Ising universality class for infinitesimally small $\xi$. However, in the finite-size systems simulated, we are expected to see the following features. First, at small $\xi$, the critical exponents obtained from finite-size scaling applied to data of small system sizes are closer to 3D Ising universality class, while those obtained from data of larger system sizes are closer to $N=8$ Chiral Ising universality class. Second, the finite-size scalings of $\langle m^{2}\rangle$ data near QCPs from the same system size for increasing $\xi$ parameters should arrive at critical exponents closer to those of $N=8$ Chiral Ising universality class. These features imply that there exist a length scale $L_{c}$ for each $\xi$, below which the critical exponents obtained from scaling should be close to 3D Ising universality class and otherwise close to $N=8$ Chiral Ising universality class. Clearly, the $L_{c}$ should be larger for smaller $\xi$ and becomes smaller for larger $\xi$.

The finite-size scalings of $\langle m^{2}\rangle$ data from $\xi=0.00,0.25,0.50,0.75$ are shown in Fig. S13. First of all, we have applied the critical exponents $\nu=0.629971(4),\eta=0.026298(2)$ to collapse the $\langle m^{2}\rangle$ data for $\xi=0$ and we have obtained $h_{c}/J=3.046(3)$, which is well-consistent with previous results Xu et al. (2017, 2016) and the data collapse also have high quality. Then for $\xi=0.25,0.50,0.75$, we have performed the data collapses with free $\nu,\eta,h_{c}$ and data from $L=6,8,10,12$ and $L=8,10,12$, respectively. For $\xi=0.25$, we can observe a dramatic change of $\eta$ exponent from $\eta=0.30(2)$ to $\eta=0.61(5)$, by simply adding the $\langle m^{2}\rangle$ data of $L=12$ system. This fact explicitly shows that for smaller system sizes like $L=6,8$, the critical behavior is closer to 3D Ising universality, while it’s more likely to be $N=8$ Chiral Ising universality class for $L=12$ system. We can conclude $L_{c}\approx 10$ for $\xi=0.25$, signifying the finite-size scaling crossover behavior. The critical exponents from data collapses of $L=6,8,10,12$ and $L=8,10,12$ are both converging to the numbers in $N=8$ Chiral Ising universality class, for increasing $\xi=0.25,0.50,0.75$ parameters. This indicates the decreasing $L_{c}$ length scale in the finite-size scaling for increasing $\xi$. Third, we can also observe that the $\eta$ exponent suffers much stronger finite-size effect than that of $\nu$ exponent, especially for small $\xi$. This is simply due to the fact that the two universality classes have similar $\nu$ exponents while $\eta$ exponent differs a lot. As a result, the crossover between these two universality classes gives much larger deviation for $\eta$ than $\nu$. All of these numerical results support the finite-size scaling crossover behavior in the ground-state phase diagram.