Edge stability and edge quantum criticality in 2D interacting topological insulators

TL;DR

This paper uses quantum Monte Carlo simulations to study the stability of helical edge states in 2D topological insulators with interactions, revealing conditions for robustness and edge magnetic order, and uncovering universal critical properties.

Contribution

First sign-free QMC simulation of an interacting 2D TI model showing edge state robustness and magnetic quantum criticality with fractionalized charges at domain walls.

Findings

Helical edge states are robust against two-particle backscattering with moderate Coulomb repulsion.

Strong Coulomb interactions can gap edge states and induce magnetic order.

Fractionalized e/2 charge appears at magnetic domain walls on edges.

Abstract

Robustness of helical edge states in 2D topological insulators (TI) against strong interactions remains an intriguing issue. Here, by performing the first sign-free quantum Monte Carlo (QMC) simulation of the Kane-Mele-Hubbard-Rashba model which describes an interacting 2D TI with two-particle backscattering on edges, we verify that the gapless helical edge states are robust against a finite range of two-particle backscattering when the Coulomb repulsion is not strong. However, when the Coulomb repulsion is strong enough, the helical edge states can be gapped by infinitesimal two-particle backscattering, resulting in edge magnetic order. We further reveal universal properties of the magnetic edge quantum critical point (EQCP). At magnetic domain walls on edges, we find that a fractionalized charge of e/2 emerges. Implications of our results to recent transport experiments in the InAs/GaSb quantum well, which is a 2D TI with strong interactions, will also be discussed.

Equations

$$H=-t\sum _{⟨ij⟩\sigma }{c}_{i\sigma }^{†}{c}_{j\sigma }+\sum _{⟨⟨ij⟩⟩}(i\lambda {\nu }_{ij}{\sigma }_{\alpha \beta }^z{c}_{i\alpha }^{†}{c}_{j\beta }+H.c.)$$

$$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ +U\sum_{i}n_{i\uparrow}n_{i\downarrow}+V\sum_{\left\langle ij\right\rangle}\big{(}c^{\dagger}_{i\uparrow}c^{\dagger}_{j\uparrow}c_{i\downarrow}c_{j\downarrow}+H.c.\big{)},$$

$$\displaystyle H=\int dx\Big{[}v_{F}(\psi^{\dagger}_{R\uparrow}i\partial_{x}\psi_{R\uparrow}-\psi^{\dagger}_{L\downarrow}i\partial_{x}\psi_{L\downarrow})$$

$$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ +g_{2}\psi^{\dagger}_{R\uparrow}\psi_{R\uparrow}\psi^{\dagger}_{L\downarrow}\psi_{L\downarrow}+g(\psi^{\dagger}_{R\uparrow}\psi^{\dagger}_{R\uparrow}\psi_{L\downarrow}\psi_{L\downarrow}+H.c.)\Big{]},\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\$$

$$\displaystyle{\cal L}=\frac{1}{2\pi K}\left[\frac{1}{v}(\partial_{\tau}\phi)^{2}+v(\partial_{x}\phi)^{2}\right]-\frac{g}{(2\pi\delta)^{2}}\cos(4\phi),\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\$$

$$C(r)=\frac{1}{L_x}\sum _{x_i=1}^{L_x}⟨S_{i,A}^x{S}_{i+r,A}^x⟩,$$

$$H$$

$$\hat{h}=\sum _{⟨ij⟩}-\frac{t{\Delta }_{\tau }}{2}{\gamma }_i^T{\sigma }^0{\tau }^t{\gamma }_j+\sum _{⟨⟨ij⟩⟩}\frac{i\lambda {\nu }_{ij}{\Delta }_{\tau }}{2}{\gamma }_i^{T}i{\sigma }^z{\tau }^0{\gamma }_j+\sum _i{\lambda }_U{\phi }_i{\gamma }_i^T{\sigma }^0{\tau }^y{\gamma }_i+\sum _{⟨ij⟩,\alpha }{\lambda }_V{\varphi }_{ij}^{\alpha }{\gamma }_i^T{B}_{\alpha }{\gamma }_j$$

$$\xi =\frac{L_x}{2\pi }\sqrt{\frac{S(\vec{0})}{S(\delta \vec{k})}-1},$$

$$G_k^f(\tau )=\sum _{\sigma =\uparrow ,\downarrow }⟨c_{k\sigma }^{†}(\tau )c_{k\sigma }(0)⟩,$$

Full text

Topology has played an increasingly important role in condensed matter physics in the past few decades XGWenbook ; Fradkinbook , especially after discoveries of quantum Hall effect TKNN ; Laughlin1983 , high-temperature superconductors Anderson1987 ; Kivelson1987 ; Xiaogang-2006RMP , and recently topological insulators Hasan-Kane-10 ; Qi-Zhang-11 . Especially, tremendous progress has been made in both theoretical and experimental understandings of topological insulators (TI) with negligible or weak interactions Kane2005a ; Kane2005b ; Bernevig2006prl ; Bernevig2006science ; Fu2007a ; Moore2007 ; Roy2009 ; Qi2008 ; Schnyder2008 ; Kitaev2009 ; Molenkamp2007 ; Hasan2008 ; Chen2009 . For instance, 2D TIs are topological state of matter supporting the gapless helical edge states which are protected by time reversal (TR) symmetry when interactions are relatively weak. As strong interactions may destabilize the gapless excitations on boundaries Fidkowski2010 ; Qi2013 ; Yao2013 ; Ryu2012 ; Ashvin2013prx ; Wang2014prb ; Wang2014science ; DHLee2011 ; Xu2015prb ; Ziani2015 ; Mudry2015 ; Gefen2017 ; Sangiovanni2017 of symmetry-protected topological phases Hasan-Kane-10 ; Qi-Zhang-11 ; ZCGu2009prb ; Pollmann2010prb ; Chen2011prb ; Chen2012science ; YMLu2012 , understanding the interplay between topology and interactions, especially in experimentally accessible topological materials, remains an intriguing problem.

2D TIs were first discovered in HgTe quantum wells (QWs) Bernevig2006science ; Molenkamp2007 and later in InAs/GaSb QWs Liu2008 ; Du2011 ; Moler2014 . Although the correlation effect in the helical edge states of the HgTe QWs is relatively weak (a Luttinger liquid Haldane1981 ; Kane1992 ; Nagaosa1994 ; Giamarchibook ; Wu2006 ; Xu2006 with edge Luttinger parameter $K\approx 0.81$ Joseph2009 ; Kane2009prb , close to the noninteracting limit of $K=1$), recent transport experiments Du2015a ; Du2015b reported evidences of strong interactions in the helical edge states of the InAs/GaSb QWs (with edge Luttinger parameter $K\approx 0.22$ Du2015b ). Consequently, it is desired to thoroughly investigate the robustness of gapless helical edges states of 2D TIs with strong interactions due to its importance in both fundamental physics and future applications.

It was predicted that helical edge states are perturbatively stable against weak two-particle backscattering if the edge Luttinger parameter (which is mainly determined by local Coulomb interaction) K$$>$$\frac{1}{2} Wu2006 ; Xu2006 . When the local Coulomb interaction is strong enough such that K$$<$$\frac{1}{2}, the two-particle backscattering interaction becomes relevant and an infinitesimal two-particle backscattering can lead to edge magnetic order, which breaks time-reversal (TR) symmetry spontaneously and opens up a gap in helical edge states Wu2006 ; Xu2006 . However, such perturbative analysis of weak two-particle backscattering may not be directly applied for the case of strong interactions. Consequently, non-perturbative methods are desired to explore the effect of strong interactions in topological insulators such as the InAs/GaSb QWs.

To fill in this gap, by performing quantum Monte Carlo (QMC), we study the stability of helical edge states of 2D interacting TIs in the presence of both strong Coulomb (Hubbard) repulsion and Rashba interactions that cause two-particle backscattering on edges. Here, it is the Rashba interactions that break the $U$(1) spin-rotational symmetry and render a finite two-particle backscattering in helical edge states. Usual QMC simulations of fermionic quantum models often encounter the notorious fermion-sign-problem Loh1990 ; Troyer2005 ; Wu2005 , which prevents the application of conventional QMC BSS1981 ; Hirsch1985 ; Assaad2008 ; Berg2012 from accurately studying systems with large size and at low temperature. However, using the novel Majorana algorithm introduced by us LJY2015 , sign-free Majorana QMC simulations LJY2015 ; LJY2016 ; Xiang2016 can be performed on the 2D interacting TI model on the honeycomb lattice with both Hubbard repulsion and Rashba interaction [see Eq. (Edge stability and edge quantum criticality in 2D interacting topological insulators) below], which we call the Kane-Mele-Hubbard-Rashba (KMHR) model. Note that previous sign-free QMC can only study 2D interacting TIs without edge two-particle backscattering Wu2011 ; Assaad2011 ; Assaad2013 which cannot accurately address the issue of quantum critical behaviors of the edge magnetic ordering and that the effect of single-particle Rashba hopping term was studied by QMC in the presence of the sign problem Assaad2014prb or by other methods which involve approximations Rachel2014 .

From our sign-free QMC simulations of the KMHR model, we have verified that the helical edge states of an interacting 2D TI are robust against a finite range of two-particle backscattering when the Hubbard interaction is relatively weak. We show that the edge states spontaneously break time-reversal symmetry when the Rashba interaction exceeds a critical value and obtain the critical exponents of the edge quantum critical point (EQCP) in 1+1 dimensions. Moreover, we find that a fractionalized charge of $e/2$ emerges at the magnetic domain wall on the edges when magnetic ordering occurs in the edges of the interacting TI. To the best of our knowledge, this is the first time that the effect of two-particle backscattering can be simulated numerically-exactly in 2D interacting TI by the unbiased and nonperturbative approach of sign-free QMC.

The 2D interacting TI model: We first introduce a minimal model on the honeycomb lattice describing a 2D interacting topological insulator allowing two-particle backscattering on edges:

[TABLE]

where $c^{\dagger}_{i\sigma}$ creates an electron on site $i$ with spin polarization $\sigma=\uparrow,\downarrow$, $n_{i\sigma}=c^{\dagger}_{i\sigma}c_{i\sigma}$ is the particle number operator. The first term in Eq. (Edge stability and edge quantum criticality in 2D interacting topological insulators) represents the usual nearest-neighbor hopping while the second one with $\nu_{ij}$=$\pm 1$ is the Kane-Mele spin-orbital coupling (SOC) Kane2005a . The Kane-Mele SOC reduces the spin rotational symmetry from $SU(2)$ to $U(1)$. The Hubbard $U$ term mimics the effects of Coulomb repulsion. The Hamiltonian in Eq. (Edge stability and edge quantum criticality in 2D interacting topological insulators) without the last term is called Kane-Mele-Hubbard (KMH) model, for which a sufficiently strong $U$ would favor AF ordering with magnetic moments in the $xy$ planeRachel2010 ; JXLi2011 ; Wu2011 ; Assaad2011 ; Assaad2013 ; Sangiovanni2015 . For the $V$ term in Eq. (Edge stability and edge quantum criticality in 2D interacting topological insulators), we call it “Rashba interaction” because this two-particle spin-flip term can be induced by the Rashba SOC. This is why the Hamiltonian in Eq. (Edge stability and edge quantum criticality in 2D interacting topological insulators) is denoted as the “Kane-Mele-Hubbard-Rashba” model. The Rashba interaction only respects a residual $Z_{2}$ spin symmetry; namely $(-1)^{N_{\uparrow}}$ is conserved where $N_{\uparrow}$ is the number of spin-up electrons. In other words, a finite Rashba interaction $V$ further reduces the spin symmetry from $U(1)$ to $Z_{2}$, which introduces two-particle spin-flip backscattering in the helical edge states. Note that the single-particle backscattering in helical edge states is forbidden even in the presence of finite Rashba interaction due to the TR symmetry.

Among all the two-particle scattering processes respecting the TR symmetry, only the two-particle spin-flip backscattering has the potential to open up a gap in the helical edge states. Consequently, this Rashba interaction is the leading term which can destabilize the gapless helical edge states of the 2D TI. It was thought for many years that this model with both Hubbard repulsion and spin-flip interaction cannot be simulated by sign-problem-free QMC Wu2011 ; Assaad2011 ; Assaad2013 . However, employing the Majorana algorithm recently introduced in Ref. LJY2015 , we show that the KMHR model in Eq. (Edge stability and edge quantum criticality in 2D interacting topological insulators) can be simulated by QMC without encountering the notorious fermion sign problem such that accurate large-scale QMC simulations of the model can be performed to investigate the edge stability of the 2D interacting TI.

Here we briefly discuss how Majorana representation can help solve the sign problem of the KMHR model. First, we introduce the Majorana representation of spin-1/2 electrons: $c_{i\sigma}=\frac{1}{2}(\gamma^{1}_{i\sigma}+i\gamma^{2}_{i\sigma})$, $c^{\dagger}_{i\sigma}=\frac{1}{2}(\gamma^{1}_{i\sigma}-i\gamma^{2}_{i\sigma})$, where $\gamma^{\tau}_{i\sigma}$ are Majorana fermions operators with $\tau=1,2$ representing Majorana index and $\sigma=\uparrow,\downarrow$ spin index. In the Majorana representation, we perform Hubbard-Stratonovich transformations such that the decoupled Hamiltonian at imaginary time $\tau$ respects two anti-commuting anti-unitary symmetries: $T^{-}=i\sigma^{y}\tau^{x}K$ and $T^{+}=\sigma^{x}\tau^{x}K$ (see the SM for details). According to the Majorana TR principle LJY2016 ; Xiang2016 , the QMC simulation using the Majorana algorithm is sign-free!

The quantum phase diagram: We perform large-scale sign-free projector QMC Sugiyama1986 ; Sorella1989 ; White1989 to study the ground state properties of the 2D interacting topological insulator with two-particle spin-flip scattering. In our simulations, we set $\lambda=0.1t$ unless noted otherwise. To investigate the stability of helical edge states with respect to two-particle spin-flip scattering, we consider the system with periodic boundary condition (PBC) along the $x$-direction but open boundary condition (OBC) along the $y$-direction, as shown in Fig. 1(a). Since the bulk is fully gapped, we focus on the helical edge states, which can be effectively described by the LL theory Haldane1981 ; Kane1992 ; Giamarchibook with the following Hamiltonian in the continuum:

[TABLE]

where $\psi_{R\uparrow}$ and $\psi_{L\downarrow}$ are annihilation operators of spin-$\uparrow$ right-moving and spin-$\downarrow$ left-moving electrons on the edge, respectively. In Eq. (Edge stability and edge quantum criticality in 2D interacting topological insulators), the last term represents $g(\psi^{\dagger}_{R\uparrow}(x)\psi^{\dagger}_{R\uparrow}(x+\delta)\psi_{L\downarrow}(x)\psi_{L\downarrow}(x+\delta)+H.c.)$, where $\delta$ is the short-distance cutoff. Here $v_{F}$ is the Fermi velocity of non-interacting electrons on edges and is approximately proportional to the Kane-Mele SOC $\lambda$. $g_{2}$ is the two-particle forward-scattering and is proportional to $U$; $g$ is the two-particle backscattering induced by the Rashba interaction $V$. Note that single-particle backscattering is forbidden by the TR symmetry even in the presence of Rashba SOC and that in Eq. (Edge stability and edge quantum criticality in 2D interacting topological insulators) we neglect the interactions $g_{4}[\psi^{\dagger}_{R\uparrow}(x)\psi_{R\uparrow}(x)\psi^{\dagger}_{R\uparrow}(x+\delta)\psi_{R\uparrow}(x+\delta)+R\!\to\!L]$ which only renormanizes the Fermi velocity.

The fermionic edge theory in Eq. (Edge stability and edge quantum criticality in 2D interacting topological insulators) can be bosonized with the following sine-Gordon Lagrangian density:

[TABLE]

where $\phi$ is the boson field obtained from bosonizing the fermions on the edge Wu2006 ; Xu2006 , $v\!=\!v_{F}\sqrt{1-(\frac{g_{2}}{v_{F}})^{2}}$ is the renormalized velocity of bosons, $K=\sqrt{\frac{1-g_{2}/v_{F}}{1+g_{2}/v_{F}}}$ is the Luttinger parameter which decreases as $U$ increases, and $g$ characterizes the two-particle backscattering whose scaling dimension is $\Delta_{g}=4K$. Consequently, two-particle spin-flip backscattering is irrelevant when the Luttinger parameter K$$>$$\frac{1}{2} but relevant when K$$<$$\frac{1}{2} Kane1992 ; Wu2006 ; Xu2006 . We expect that the helical edge states are stable against infinitesimal two-particle backscattering when the Hubbard interaction $U$ is not sufficiently strong because of the irrelevance of weak two-particle backscattering for K$$>$$\frac{1}{2}. When the two-particle backscattering is relevant for $K<\frac{1}{2}$, the mass gap opens even for weak Rashba interaction $V$ and the gap is given by $\Delta\approx\frac{1}{\delta}g^{\frac{1}{2-4K}}$ for small $g$. This conclusion comes from the perturbative analysis of $g$. But, for $K>\frac{1}{2}$, strong enough $g$ should also drive a mass gap. This can be obtained from the QMC calculations below.

We first compute the values of Luttinger parameter of the helical edge states for different Hubbard interaction $U$ while setting the Rashba interaction $V=0$. The edge Luttinger parameter $K$ can be extracted by measuring the correlation function of magnetic order parameter:

[TABLE]

where $S^{x}_{i,A}=\frac{1}{2}(c^{\dagger}_{i\uparrow}c_{i\downarrow}+c^{\dagger}_{i\downarrow}c_{i\uparrow})$ are the spin operators in the $A$-sublattice of the unit cell $i=(x_{i},y_{i})$ on the edge $y_{i}=L_{y}$ ($L_{x}$ and $L_{y}$ are the number of unit cells along the $x$ and $y$ directions, respectively). Note that due to the U(1) spin-rotational symmetry at $V=0$, using $S^{x}$, $S^{y}$, or $\vec{S}\cdot\hat{n}$ with any $\hat{n}$ in Eq. (4) would give rise to the same value of $K$. In other words, we compute the correlation functions between edge sites of the $A$ sublattice, which are the tips on the top zigzag edge of the honeycomb lattice, as shown in Fig. 1(a). In our QMC simulations, we choose $L_{x}=2L_{y}=L$. According to the LL theory, the magnetic correlation function in Eq. (4) scales as $C(r)\sim r^{-2K}$ for large $r$ . In order to reduce the finite-size effect, we extract the Luttinger parameter from fitting C(r_{\textrm{max}}$$=$$\frac{L}{2})\sim L^{-2K}, which is the correlation function at the largest possible separation between two sites along the edge of the lattice. The obtained values of $K$ as a function of $U$ is shown in Fig. 2(a). It is clear that $K$ decreases as $U$ is increased, as expected. At $U=U^{\textrm{edge}}_{c}\approx 1.85t$, the Luttinger parameter $K=K_{c}=\frac{1}{2}$. The bulk AF magnetic ordering can occur for $U>U^{\text{bulk}}_{c}$ where $U^{\textrm{bulk}}_{c}\approx 5.15t$.

From bosonization analysis, the helical edge states should be stable against a finite range of two-particle backscattering for U$$<$$U^{\textrm{edge}}_{c}. Thus, for U$$<$$U^{\textrm{edge}}_{c}, we expect that a sufficiently large Rashba interaction can destabilize the helical edge states by inducing magnetic order. To illustrate this, we set $U=1.0t$ (the corresponding $K\sim 0.81$ in the limit of $V=0$) and vary the Rashba interaction $V$ to study the edge instability against the two-particle backscattering. The Rashba interaction with $V>0$ would favor magnetic ordering in the $x$ direction more than the $y$ direction. We use the RG-invariant quantity $R^{\textrm{edge}}_{c}$ in terms of second-moment correlation length Assaad2015 of edge correlation function $C(r)$ to identify the EQCP ($R^{\textrm{edge}}_{c}$ is defined in the SM). At the putative EQCP, $R^{\textrm{edge}}_{c}$ should cross at the same point for different system sizes. From the results of RG-invariant quantity $R_{c}^{\textrm{edge}}$, as shown in Fig. 2(b), we obtain the edge critical value $V^{\textrm{edge}}_{c}\approx 0.515t$. For V$$>$$V^{\textrm{edge}}_{c}, the edge ferromagnetic (FM) ordering spontaneously occurs (the magnetic moments on sites in the $A$ sublattice point to the same direction Wu2015 ). When $V$ is further increased such that $V>V^{\textrm{bulk}}_{c}\approx 0.73t$, the bulk AF ordering occurs. The magnetic order parameter M$$=$$\left\langle S^{x}_{i}\right\rangle for $U=1.0t$ as a function of $V$ is shown in Fig. 2(c).

For U^{\textrm{edge}}_{c}$$<$$U$$<$$U^{\textrm{bulk}}_{c} or equivalently K$$<$$\frac{1}{2}, an infinitesimal two-particle backscattering on the edges can gap out helical edge states by inducing magnetic order Wu2006 ; Xu2006 ; Wu2011 . For instance, we set $U=2.25t$ and vary the Rashba interaction $V$ to study the edge instability. As shown in Fig. 2(d), the edge magnetic ordering already occurs for a very weak $V$, indicating that the two-particle backscattering is relevant and that edge magnetic order is induced as long as V>V^{\textrm{edge}}_{c}$$=[math]. When $V$ is further increased such that $V>V^{\textrm{bulk}}_{c}\approx 0.45t$, the bulk AF will be induced and the whole system breaks TR symmetry. The global quantum phase diagram is summarized in Fig. 1(b). There is a line of edge QCP and a separated line of bulk QCP. In the intermediate range, the TR symmetry breaking occurs only on the edge while the bulk is still TR invariant.

Edge quantum criticality: It would be interesting to analyze the edge quantum critical behaviors Subirbook in topological phases to study exotic phenomena such as the emergence of supersymmetry Ashvin2014 ; Ponte2014 ; Li2016arxiv ; Jian2017prl . The EQCP can be driven by either $U$ or $V$. Here, we consider a weak but finite spin-flip interaction $V=0.3t$ and vary $U$ in the QMC simulations to study the EQCP. Similarly, we use RG-invariant quantity $R^{\textrm{edge}}_{c}$ to identify the EQCP. The results of our sign-free MQMC simulation on $R^{\textrm{edge}}_{c}$ (see the SM for details) clearly show that the edge quantum phase transition occurs at $U^{\textrm{edge}}_{c}\approx 1.72t$. Moreover, we computed the evolution of the single-particle gap as the function of $U$ (see the SM for details). When the TR symmetry is spontaneously broken at $U>U^{\textrm{edge}}_{c}$, the helical edge state becomes gapped. Close to the $U^{\textrm{edge}}_{c}$, the gap (equivalently the inverse correlation length) scales as \xi^{-1}$$\sim$$\Delta$$\sim$$\exp[-\frac{A}{\sqrt{U-U^{\textrm{edge}}_{c}}}], where $A$ is some constant.

We then investigate the quantum critical behaviours of this EQCP. When $U>U^{\textrm{edge}}_{c}$, the edge develops a finite magnetic order breaking the TR symmetry. The transition should belong to the universality class of Kosterlitz-Thouless transitions in 1+1 dimensions Fradkinbook . By fitting the correlation function of magnetic order at the EQCP: $C(\frac{L}{2})\sim L^{-\eta}$, we obtain the anomalous dimension of order parameter bosons: $\eta=0.94\pm 0.03$, as shown in Fig. 3(a). According to the bosonization analysis, the correlation length critical exponent $\nu$ of this EQCP should be infinite. In order to verify it, we perform data collapse of the edge magnetic structure factor $M_{2}=\frac{1}{L}\sum_{r}C(r)$ by the scaling function $M_{2}L^{\eta}=F(L^{1/\nu}(U-U^{\textrm{edge}}_{c}))$, $F$ is an unknown function. When $1/\nu=0$ and $\eta=0.94$, various points $(M_{2}L^{\eta},L^{1/\nu}(U-U^{\textrm{edge}}_{c}))$ of different $U$ around $U^{\textrm{edge}}_{c}$ and different $L$ can collapse to a single curve, as shown in Fig. 3(b). Consequently, we verified that the EQCP belongs to the Kosterlitz-Thouless universality in 1+1 dimensions with $\eta=0.94\pm 0.03$ and $\nu=\infty$.

Fractionalized charge: The edge magnetic ordering can render some exotic physics because the helical edge states of 2D TI has only half degrees of freedom of a one-dimensional system of spin-1/2 electrons respecting the TR symmetry. For instance, it was predicted theoretically that a fractionalized charge of $e/2$ can be induced at the magnetic domain wall in the helical edge states of a 2D TI system Qi20082 , which is a realization of the Jackiw-Rebbi mass soliton of the 1+1 Dirac theory Jackiw1976 . To verify the domain-wall fractionalized charge, an infinitesimal external magnetic field polarized in opposite directions in magnetic ordered phase to create a magnetic domain wall, as shown in Fig. 4(a).

We set $V=0.3t$ and $U=4.0t>U^{\textrm{edge}}_{c}$ such that the edge is magnetically ordered while the bulk still respects the TR symmetry. The lattice sizes we choose are $L_{x}=60$ and $L_{y}=6$. Two magnetic edge domain walls are generated at $x$=15 and $x$=45, as shown in Fig. 4(a), by imposing a weak external magnetic field $B_{z}=\pm 0.05$. We compute distribution of electron density $\rho(x)$ on the edge while the average electron density per site is $\bar{\rho}=1$. The excess charge localized at the magnetic domain wall is given by $n_{\textrm{dw}}=\sum_{x}\rho(x)-\bar{\rho}$, where $x$ is summed only around the domain wall. Our QMC results show that $n_{\textrm{dw}}\approx e/2$, as shown in Fig. 4(b). This provides robust numerical evidences of the emergent fractionalized charge at magnetic domain wall in the interacting helical edge states with spontaneous TR symmetry breaking. Such magnetic domain walls carrying fractionalized $e/2$ charge can be potentially measured in future experiments.

Discussions and conclusions: We now discuss the implication of our results to recent experiments done in the 2D TI InAs/GaSb QW Du2015a ; Du2015b . Transport experiments in Ref. Du2015b reported the temperature dependence of its edge conductance: $G_{xx}\sim T^{0.32}$ for sufficiently low temperature, which indicates that its edge Luttinger parameter is $K\approx 0.22$ assuming that the conductance is dominated by instantons at an impurity with fractionalized charges $e/2$ Joseph2009 . For $K<1/2$, the helical edge states can also be gapped by spontaneously breaking the TR symmetry on edges when the two-particle backscattering is allowed (namely when the Fermi level is tuned to near the Dirac point of the helical edge states). At finite temperature, charge transport could have contributions from magnetic domain wall with fractionalized charge of $e/2$, as shown above. It would be interesting to study transport properties of helical Luttinger liquids Oreg2012 ; Crepin2012 ; Schmidt2012prl ; Mirlin2014 ; Glazman2016prb ; Ziani2017np with two-particle backscattering in the region of $K<\frac{1}{2}$ by sign-free QMC simulations, which will be deferred to future works. As the helical edge states coupled with superconductors and with magnetism can support Majorana zero modes Hasan-Kane-10 ; Qi-Zhang-11 , it would also be interesting to study the effect of strong interactions on Majorana zero modes Kane2014prl ; Schmidt2015prb .

In conclusion, we have proposed the KMHR model to describe 2D interacting TIs allowing two-particle spin-flipping backscattering. It is the first two dimensional model with two-particle backscattering which can be simulated by sign-free QMC using the Majorana representation. Our large-scale QMC simulations of this model have shown that the helical edge states are robust under weak Rashba spin-flipping interaction when the Hubbard (or Coulomb) repulsion is not too strong. When Hubbard repulsion is strong enough, the gapless helical edge states are unstable again even infinitesimal two-particle backscattering by spontaneously breaking the TR symmetry and forming magnetic order at edges. The critical behaviors of this EQCP has also been obtained. Our work may provide a promising new direction to study the boundary stability and quantum criticality in topological phases of matter by non-perturbative approaches.

Acknowledgements: We sincerely thank Yi-Fan Jiang, Dung-Hai Lee and Shou-Cheng Zhang for helpful discussions. This work is supported in part by the NSFC under Grant No. 11474175 (ZXL and HY) and by the MOST of China under Grant No. 2016YFA0301001 (HY).

I Supplementary Materials

I.1 I. Proof of Sign-problem-free of interacting topological insulator model in Majorana representation

To prove the absence of sign-problem in QMC simulations of the KMHR model in Eq. (1) of the main text, we introduce the Majorana representation of spin-1/2 electrons: $c_{i\sigma}=\frac{1}{2}(\gamma^{1}_{i\sigma}+i\gamma^{2}_{i\sigma})$, $c^{\dagger}_{i\sigma}=\frac{1}{2}(\gamma^{1}_{i\sigma}-i\gamma^{2}_{i\sigma})$, where $\gamma^{\tau}_{i\sigma}$ are Majorana fermions operators with $\tau=1,2$ representing Majorana index and $\sigma=\uparrow,\downarrow$ representing spin index. First, we perform the particle-hole transformation on $\sigma=\downarrow$ electrons: $c_{i\downarrow}\rightarrow(-1)^{i}c_{i\downarrow}^{\dagger}$. Under this transformation, the NN hopping term and NNN Kane-Mele SOC term are invariant. The Hubbard interaction term changes sign and Rashba interaction becomes $V(c^{\dagger}_{i\uparrow}c^{\dagger}_{i\downarrow}c^{\dagger}_{j\downarrow}c^{\dagger}_{j\uparrow}+h.c)$. After the particle-hole transformation, the Hamiltonian in Eq. (1) can be rewritten in Majorana representation as:

[TABLE]

where $\gamma^{T}_{i}\equiv(\gamma^{1}_{i\uparrow},\gamma^{2}_{i\uparrow},\gamma^{1}_{i\downarrow},\gamma^{2}_{i\downarrow})$, $B_{1}=\sigma^{z}\tau^{z},B_{2}=i\sigma^{0}\tau^{z},B_{3}=\sigma^{0}\tau^{x},B_{4}=i\sigma^{z}\tau^{x}$. Upon Trotter decomposition and usual Hubbard-Stratonovich transformations, the decoupled Hamiltonian at imaginary time $\tau$ can be written in Majorana representation as:

[TABLE]

where $\Delta_{\tau}$ is imaginary time slice of the Trotter decomposition, $\varphi_{i}$ are imaginary-time dependent auxiliary fields on site $i$ and $\phi^{\alpha}_{ij}$ are imaginary-time dependent auxiliary fields on bond $\left\langle ij\right\rangle$. The decoupled Hamiltonian in Eq. (S2) possesses two anti-commuting Majorana-time-reversal symmetries: $T^{-}=i\sigma^{y}\tau^{x}K$ and $T^{+}=\sigma^{x}\tau^{x}K$. According to the Majorana TRS principle for sign-problem-free QMC, it belongs to Majorana class and is then sign-problem-free.

I.2 II. Details of the projector Majorana QMC simulations

We use projector QMC in the Majorana representation to investigate the ground state properties of the 2D interacting topological insulator described by the Hamiltonian in Eq. (1). In the projector QMC, the expectation value of an observable $O$ in the ground state can be evaluated as: $\frac{\left\langle\psi_{0}\right|O\left|\psi_{0}\right\rangle}{\left\langle\psi_{0}\mid\psi_{0}\right\rangle}=\lim_{\Theta\rightarrow\infty}\frac{\left\langle\psi_{T}\right|e^{-\Theta H}Oe^{-\Theta H}\left|\psi_{T}\right\rangle}{\left\langle\psi_{T}\right|e^{-2\Theta H}\left|\psi_{T}\right\rangle}$, where $\psi_{0}$ is the true ground state wave function and $\left|\psi_{T}\right\rangle$ is a trial wave function which should have a finite overlap with the true ground state wave function. Note that $\Theta$ is not the inverse temperature but the projection parameter. Although $\Theta\to\infty$ is needed to reach the exact ground state, in numerically calculations a sufficient large $\Theta$ works for practical purposes of obtaining physical quantities with required accuracy. Because of the absence of sign-problem, we can perform large-scale QMC simulations with large system sizes and sufficiently large $\Theta$. In the study of the bulk quantum phase transition, we use periodic boundary condition. In the study of edge quantum phase transition and edge Luttinger parameters, we use the periodic boundary condition in the $x$ direction but open boundary condition in the $y$ direction. In our QMC simulations, the imaginary-time projection parameter is $\Theta=60/t$ for the systems with torus boundary conditions. In the cases of cylinder boundary conditions, most systems are computed using $\Theta=75/t$ and some systems with large systems size or near critical points are computed using $\Theta=100/t$. We have checked that all the results stay nearly the same when larger $\Theta$ are used, which ensures desired convergence to the limit of $\Theta\to\infty$. We set $\Delta_{\tau}=0.05/t$ and the results do not change if we use smaller $\Delta_{\tau}$.

I.3 III. Details of the RG-invariant quantity $R^{\textrm{edge}}_{c}$ and single-particle gap

We use the RG-invariant ratio $R^{\textrm{edge}}_{c}$ in terms of second-momentum correlation length of edge AF magnetic order to identify the EQCP. The second-momentum correlation length of edge AF magnetic order is defined as

[TABLE]

where $\delta\vec{k}$ is the minimum lattice momentum internal for systems of size $L_{x}$. Here $S(\vec{k})$ is the structure factor of magnetic order at edge: $S(\vec{k})=\frac{1}{L_{x}}\sum_{r}C(r)e^{i\vec{k}\cdot\vec{r}}$. The ratio $R^{\textrm{edge}}_{c}$ is defined as $R^{\textrm{edge}}_{c}=\frac{\xi}{L}$. At the putative EQCP, $R^{\textrm{edge}}_{c}$ is RG-invariant, such that it should cross at the same point for different system sizes. In the disordered phase, the ratio $R^{\textrm{edge}}_{c}=\frac{\xi}{L}$ should decrease as the system size is increased. The trend is opposite in the ordered phase. Thus, this RG-invariant ratio is a powerful tool to accurate identify the QCP. Here, we fix the strength of spin-flip Rashba interaction $V=0.3t$ and vary $U$ to study the EQCP. The RG-invariant ratio $R^{\textrm{edge}}_{c}$ is evaluated by sign-free QMC. The results (shown in Fig. S1) explicitly show that the EQCP occurs at $U\approx 1.72t$.

When the TR symmetry is spontaneously broken at edge by form magnetic ordering, the helical edge state becomes gapped. The single-particle gap can be obtained in QMC through measuring the tails of the imaginary-time displaced Green’s function:

[TABLE]

where $\tau$ represents imaginary time, $c^{\dagger}_{k\sigma}(\tau)=e^{\tau H}c_{k\sigma}^{\dagger}e^{-\tau H}$, and $k$ is the edge momentum. Here $k$ is the momentum in $x$ direction and $c^{\dagger}_{k}=\frac{1}{L_{x}}\sum_{x=1}^{L_{x}}c^{\dagger}_{(x,y=L_{y})}e^{ikx}$. The single-particle gap $\Delta_{sp}$ corresponds to the single-particle excitation energy at $k=0$, which can be obtained from $G^{f}_{k=0}(\tau)\propto e^{-\tau\Delta_{sp}}$ when $\tau$ is large enough. We evaluate the single-particle gap as the function of $U$. The result clearly shows that when TR symmetry is spontaneously broken at $U>U_{c}^{\textrm{edge}}$, the single-particle gap in the helical edge states is opened.

I.4 IV. Numerical results of bulk AF quantum phase transition

When the Hubbard interaction $U$ or spin-flip Rashba interaction $V$ is strong enough, the spontaneous TR symmetry breaking should also occurs in the bulk. We use similar techniques including the RG-invariant ratio and data collapse to study the bulk quantum phase transition and analyze the quantum critical point. We fix $V=0.3t$ and vary the value of $U$. Similar to the case of EQCP, we measure RG-invariant ratio $R^{\textrm{bulk}}_{c}$ of the bulk AF magnetic order to determine the bulk quantum phase transition point. From the crossing point of $R^{\textrm{bulk}}_{c}$ in different systems sizes, we identify the bulk phase transition point $U_{c}^{\textrm{bulk}}=5.15t$, which is much larger than the edge one $U_{c}^{\textrm{edge}}=1.72t$. It clearly indicates that the edge phase transition occurs when bulk is still disordered. We also employed the data collapse technique to study the quantum critical behaviour of this bulk AF phase transition. The structure factor of bulk AF order $M_{2}$ at different values of $U$ and for different systems should collapse to a single curve $M_{2}L^{1+\eta}=F(L^{1/\nu}(U-U^{\textrm{bulk}}_{c}))$. Consequently, we verified that the bulk AF quantum phase transition belongs to the Ising universality in 2+1 dimensions with $1+\eta=1.01\pm 0.02$ and $\nu=0.63\pm 0.03$.