Combined spontaneous symmetry-breaking and symmetry-protected topological order from cluster charge interaction

The study of symmetry-protected topological states in presence of electron correlations has recently aroused great interest as rich and exotic phenomena can emerge. Here, we report a concrete example by employing large-scale unbiased quantum Monte Carlo study of the Kane-Mele model with cluster charge interactions. The ground-state phase diagram for the model at half filling is established. Our simulation identifies the coexistence of a symmetry-protected topological order with a symmetry-breaking Kekul$\acute{e}$ valence bond order and shows that the spontaneous symmetry-breaking is accompanied by an interaction-driven topological phase transition (TPT). This TPT features appearance of zeros of single-particle Green's function and gap closing in spin channel rather than single-particle excitation spectrum, and thus has no mean-field correspondence.

On the other hand, it's well known that electronelectron interactions can also be the driving force of various GL phase transitions and corresponding long-range ordered phases, which are characterized by spontaneous symmetry breaking and appearance of long-range orders.Representative examples are charge-density-wave [30][31][32] , antiferromagetic [33][34][35] , and valence-bond-solid 36 phases in fermion Hubbard models, spontaneously breaking inversion, spin rotation, and translation symmetries, respectively.These disordered-ordered quantum phase transitions can be either of first order or continuous, which, for the latter case, belong to certain universality class, depending on the dimension and symmetry.
With these in mind, it will be of great interest to demystify whether topological symmetry protection and spontaneous symmetry breaking can coexist in a single quantum phase in fermion systems, namely an SPT phase with symmetry breaking and long-range order.Furthermore, the corresponding TPT and disordered-ordered phase transition in such a system should be quite different from the conventional ones mentioned above.Such exotic phases have already been studied in bosonic systems 37 , termed as symmetry-breaking topological insulators, while the fermionic counterpart is missing.
In this paper, we present a positive answer to above question via a concrete example of interacting fermion model studied by numerically exact, large-scale quantum Monte Carlo (QMC) simulations.Our QMC results have revealed a topologically nontrivial Kekulé valence bond solid (KVBS) phase, which spontaneously breaks Z 3 symmetry (three ways to form a KVBS long-range order) and has spin Chern number C s = −1.We have also found coinciding KVBS phase transition and TPT from C s = +1 to C s = −1 inside quantum spin-Hall insulator (QSHI) phases.Moreover, this TPT is accompanied by appearance of zeros of single-particle Green's function.The spin channel becomes critical while the single-particle gap remains finite across this exotic phase transition, which makes this TPT very different from the standard TPT in free fermion systems.
The rest of this paper is organized as follows.In Sec.II, we first present the interacting fermion model we studied and the QMC algorithm we employed.The physical quantities calculated in this work is also briefly reviewed.Then, the QMC simulation results, as the main part of this work, are presented and discussed in details in Sec.III.Finally, Sec.IV summarizes this work, and discusses the possible extensions in future works.

A. Kane-Mele Model with cluster charge interaction
Our model describes a topological insulator of fermions with cluster charge interaction on a honeycomb lattice.The Hamiltonian contains tight-binding and interaction parts Ĥ = Ĥ0 + ĤU as where i, j represent the lattice sites, α, β =↑, ↓ label fermion spins, and λ is the strength of spin-orbit coupling.Ĥ0 is the Kane-Mele model, consisting of the nearest-neighbor (NN) hopping and the intrinsic spinorbit coupling (SOC), and the factor ν ij = −ν ji = ±1 depends on the orientation of the next-nearest-neighbor bonds as demontrated in Fig. 1(a).ĤU stands for the cluster charge interaction in which the summation runs over all the hexagons on honeycomb lattice.Q = i∈ ni is the total charge operator in a hexagon with ni = α c + iα c iα .Throughout this work, we set t as the energy unit for simplicity.
The competition between the nontrivial band topology and the electron correlation in our model in Eq. ( 1) is quite explicit.First of all, the Kane-Mele model has a QSHI ground state with counterpropagating edge states 14,38 .This QSHI phase is an SPT phase protected by U (1) spin × U (1) charge Z T 2 (with Z T 2 as the timereversal) symmetry, which results in Z classification.The appropriate topological invariant for describing the QSHI phase is the spin Chen number C s = (C ↑ −C ↓ )/2 with C ↑ and C ↓ as the Chern numbers in spin-up and down channels 24,25 , respectively.Furthermore, this QSHI phase is stable against weak and local interactions, and thus it survives in small U region in our model.As for the cluster charge interaction, it has been shown 39 that it favors a KVBS and antiferromagnetic insulating phases at intermediate and strong interactions, respectively, when the SOC term in Eq. ( 1) is turned off.In the thermodynamic limit, the KVBS phase breaks the Z 3 symmetry, and forms the pattern of alternating strong and weak bonds like the one (of three) shown in the inset of Fig. 1(c).In this work, we only concentrate on the KVBS phase induced by intermediate cluster charge interaction.It's explicit that the broken symmetry in the KVBS phase and the symmetry protecting the QSHI phase in Kane-Mele model are in different symmetry sectors.This means that these two phases may coexist in the phase diagram, which is the key point that we engaged to demystify in this paper by employing numerically unbiased large-scale QMC simulations.

B. Projector quantum Monte Carlo method
We apply the projector QMC (PQMC) method, the zero-temperature version of the determinantal QMC algorithm, to study the ground-state properties of the model in Eq. ( 1) at half-filling, i.e., one electron per site on average.The PQMC algorithm calculates the groundstate expectation results of static (equal-time) observables within the projected wavefunction as where |Ψ T is a trial wavefunction nonorthogonal to the true ground state of the many-body system.With a large enough but finite projection parameter Θ, the many-body ground state as |Ψ g = e −Θ Ĥ/2 |Ψ T can be achieved for a finite-size system.Three steps need to be done for carrying out the imaginary-time projection e −Θ Ĥ/2 |Ψ T in Eq. ( 2) by the PQMC method.First, we divide the projection parameter into M slices as Θ = M ∆ τ and e −Θ Ĥ = (e −∆τ Ĥ ) M , where ∆ τ needs to be small enough.Second, we apply the Trotter decomposition, e −∆τ Ĥ = e −∆τ Ĥ0 e −∆τ ĤU + O(∆ 2 τ ), to separate the free fermion and interaction parts in the interacting model.The Trotter error O(∆ 2 τ ) in this step is fully controlled.Third, we decouple the interaction term into free fermions coupled to auxiliary fields by Hubbard-Stratonovich (HS) transformation.After that, standard importance sampling like the Metropolis algorithm can be performed to the auxiliary-field configurations and physical observables can be computed through singleparticle Green's function.Similar formula like Eq. ( 2) exists for dynamic quantities, including the imaginarytime correlation function in both fermionic and bosonic channels.
For the model we studied in Eq. (1), we adopt the HS transformation with four-component Ising fields to decouple the cluster charge interaction 40 .Fixed at halffilling, the model in Eq. ( 1) is free of minus sign problem as we can prove that the weight of every single auxiliaryfield configuration are nonnegative due to the particlehole symmetry.In this work, we choose ∆ τ t = 0.05 and Θt ≥ 80 for the linear system sizes we simulated as L = 9, 12, 15, 18, 21 (with number of lattice sites N s = 2L 2 ).These choices of ∆ τ and Θ parameters have been tested to fully get rid of the Trotter error and converge to the true many-body ground state, respectively.
We have measured various physical observables to determine the ground-state phase diagram of model (1).We first detect the formation of long-range KVBS order by measuring the correlation functions and its structure factor in reciprocal space as where Bm = α (c + mα c m+δα + h.c.) with m, n as indices for unit cells and δ standing for one of the three NN directions.The corresponding ordering vector for the KVBS order is K and K in the Brillouin zone (shown in Fig. 1(b)).Thus, the order parameter of the KVBS can be defined as ∆ K = m e iK•Rm Bm , where the phase factor e iK•Rm signifies the Z 3 symmetry breaking.Then the location of the quantum phase transition between the disordered and KVBS ordered phases can be determined via the correlation ratio R corr = 1 − CVBS(K+q) CVBS(K) where q represents the smallest momentum in the Brillouin zone for the corresponding finite-size lattice.The histogram of the KVBS order around the phase transition is also computed to better characterize the transition.
To obtain the information of excitations of the system, we measure the dynamic single-particle Green's function and spin-spin correlation functions as where γ = A, B for the sublattices.Then the singleparticle gap ∆ sp (k) and spin gap ∆ s (k) can be extrapolated from the large-τ behavior of G(k, τ ) and S(k, τ ) as G(k, τ ) ∝ e −∆sp(k)τ and S(k, τ ) ∝ e −∆s(k)τ , respectively.For model (1), the global minimum of singleparticle gap is either at k = K(K ) or k = M depending on the SOC strength λ, while the minimal spin gap is at The TPTs driven by interactions can be dramatically different from those in the free fermion systems.To characterize the topological nature of the phases and quantum phase transition for model (1), we employ the technique of computing spin Chern number via zerofrequency single-particle Green's function [19][20][21] .All the details of this calculation are carefully demonstrated and tested in Ref. 24.For model ( 1), the time-reversal symmetry guarantees C ↑ = −C ↓ , and thus we have C s = C ↑ .

III. NUMERICAL RESULTS AND DISCUSSIONS A. Ground-state phase diagram
We first briefly summarize the λ-U phase diagram obtained by our PQMC simulations for model (1), as shown in Fig. 1(c).With finite λ, the disordered QSHI phase with C s = +1 at U = 0 extends to weak but finite interaction regime as expected.Strikingly, a QSHI phase with C s = −1 coexisting with the long-range KVBS order is discovered at intermediate interactions.Furthermore, the TPT from C s = +1 to C s = −1, which is accompanied by the presence of zeros of single-particle Green's function, coincides with the disordered-ordered phase transition characterized by the appearance of the long-range KVBS order breaking Z 3 symmetry.Our simulation results also suggest that the quantum phase transition between these two phases is of first order (at least for λ/t >= 0.2).Across the phase transition point, the excitation gap in fermionic channel (single-particle gap) remains finite while the spin gap closes and reopens.These interesting behaviors render the interaction-driven TPT in our model rather exotic and dramatically different from that in free fermion system, where pole of singleparticle Green's function appears and single-particle gap vanishes.As a comparison, a quantum phase transition between the Dirac semimetal and the KVBS phase with emergent U (1) symmetry at the transition point has been (a) The total energy density derivative ∂ Ĥ /∂U/2L 2 , and (b) the structure factor CVBS(K)/L 2 of KVBS for λ/t = 0.2 with U/t crossing the quantum phase transition.The results for linear system sizes with L = 9, 12, 15, 18, and 21 are shown.
established 39 for the λ = 0 case, which is also shown in Fig. 1(c).Thus, the presence of SOC alters the ground states and the corresponding phase transitions to a great extent.Besides, we have also confirmed the antiferromagnetic Mott insulator phase is absent in the interaction range U/t = 0 ∼ 6, for arbitrary λ/t.
In this work, we have carried out PQMC simulations for λ/t = 0.05, 0.2, 0.4 as presented in Fig. 1(c).In the following sections, we mainly demonstrate the PQMC results for λ/t = 0.05 and 0.2 cases as the only difference for λ/t = 0.4 case is the net shift in phase boundary.

B. The KVBS phase transitions
The disordered-KVBS ordered phase transition, characterized by the formation of long-range KVBS order and breaking Z 3 symmetry, is first determined from the energies and the bond-bond correlation functions defined in Eq. (3).
The results for the total energy density derivative ∂ Ĥ /∂U/2L 2 and the structure factor λ/t = 0.2 with increasing U/t are presented in Fig. 2. The explicit kinks in the results of the total energy density derivative around U/t = 3.30 indicate the location of the quantum phase transition.Similar kinks can also be observed in the results of the structure factor.An important observation is that with increasing system size, the kinks in all these results tend to evolve into jumps.This can be taken as the indication of first-order disorderedordered phase transition, as the magnitude of the KVBS order can be evaluated as m VBS = C VBS (K)/L 2 .Besides, the finite-size effect turns out to be quite strong for U/t < 3.30 region, and all the results shown in Fig. 2 converge within L = 21 for U/t > 3.32.
To explicitly determine the location of this phase transition, Fig. 3 shows the results of correlation ratio R corr and finite-size scaling for KVBS order parameter around the phase transition for λ/t = 0.2.As we know, R corr saturates to identity and vanishes in the ordered and dis- ordered phases, respectively, in the thermodynamic limit.
The crossings of R corr results of systems with different sizes can be taken as the phase transition point 13 .As shown in Fig. 3(a), though the finite-size effect is quite significant, we can observe the convergence of the crossings from U/t 3.27 for L = 9 and L = 12 to U/t 3.29 for L = 18 and L = 21.Similar to the results shown in Fig. 2, the smooth curve of R corr also tends to change into jump with increasing system size, which is also an indication of first-order phase transition.A careful finitesize scaling of KVBS order parameters within U/t = 3.28, 3.29, 3.30 are shown in Fig. 3(b).These results explicitly show that in the thermodynamic limit, there is no long-range KVBS order for U/t = 3.28, 3.29, while m VBS = 0.12(2) can be extrapolated for U/t = 3.30.This dramatic change in m VBS within such small variation of U/t parameter also suggests a first-order phase transition.We have also tried the data collapse for results of structure factors for KVBS order and no reliable results can be extracted, which contradicts with the hypothesis of continuous phase transition.From the results in Fig. 3, we can extract the location of the phase transition as U c /t = 3.295(5) for λ/t = 0.2 case.In Fig. 4, the histogram of the KVBS order parameter ∆ K for λ/t = 0.2 across the phase transition is shown with L = 12.Before the transition (U/t = 3.275), ∆ K completely distributes around zero indicating absence long-range order.As a comparison, ∆ K mainly resides in three different patches after the transition, which represents three types of ordering connected by Z 3 symmetry (and the system falls into one of them in the thermodynamic limit as the symmetry breaking).The most important observation is the coexistence of the peak around zero and the other three peaks centered around finite values with U/t = 3.32, which is actually the fingerprint of first-order phase transition.Due to the finite-size effect, the value of U/t, where this coexistence is observed, is slightly different from the transition point in the thermodynamic limit.
All the above results are well consistent and suggest that the disordered-KVBS ordered phase transition is of first order.We obtained similar results for λ/t = 0.4 with U c /t = 4.652 (5).
Here we also present the understanding of the above quantum phase transition from theoretical aspect.The Landau cubic criterion shows that the phase transitions cannot be continuous if Ginzburg-Landau (GL) free energy contains cubic terms of order parameters.This is indeed the case of the above KVBS phase with Z 3 symmetry breaking, where the cubic terms is generally allowed in GL free energy.In Ref. 41, it was shown that such kind of disordered-ordered quantum phase transition can actually be continuous due to the coupling between the KVBS order parameter and gapless Dirac fermions, namely the fermion-induced quantum critical points.This is the case for λ = 0 as shown in the phase diagram in Fig. 1(c).With finite λ, the fermions are gapped out (as shown in Sec.III D) and Landau cubic criterion applies.Thus, this can explain the first-order quantum phase transition we observed above.

C. The topological phase transitions
To characterize the topological nature of the quantum states in the phase diagram in Fig. 1(c), we calculate spin Chern number C s from G α (iω = 0, k) (with α =↑, ↓), and we also examine the eigenvalues of G α (0, k) to get more information of the structure of poles and zeros of the single-particle Green's function across the TPT.
Figures.5(a) and (b) show the results of C s and the positive eigenvalue of G ↑ (0, K) matrix for λ/t = 0.2 across the TPT.As expected, the spin Chern number obtained from direct, finite-size PQMC calculations is smooth across the transition and it's far from integer quantized values.After the interpolation calculations 24 , we indeed observe integer jumps from C s = +1 to C s = −1 across the TPT.Again, the finite-size effect is quite significant for the first several system sizes, and the best estimate of the location of TPT is U c /t = 3.29 (1), which is well consistent with the critical value obtained for disordered-KVBS ordered phase transition in Sec.III B. Thus, these two transitions coincide.More importantly, we can conclude that the KVBS phase in the phase diagram in Fig. 1(c) is also a QSHI phase with C s = −1, combining the results in Sec.III B and Fig. 5(a).On the other hand, the spin Chern number changes only because of appearance of poles or zeros in the single-particle Green's function.We have indeed observed the zeros of G ↑ (0, k) at both K and −K points at the transition point from the positive eigenvalue of G ↑ (0, K), as shown in Fig. 5 ros in G ↑ (0, k) at K and −K.In addition, we have also confirmed the absence of poles in single-particle Green's function, which is consistent with the fact that the singleparticle gap remains open across the TPT as shown in Sec.III D. Similar results are obtained for λ/t = 0.05, as presented and discussed in Appendix IV.
The topologically nontrivial KVBS phase as well as the exotic TPT discussed above is quite different from those in free fermion (or mean-field) systems in several aspects.We have also studied the corresponding meanfield system by replacing the cluster charge interaction ĤU in Eq. ( 1) by the term m VBS ij α (c + iα c jα + h.c.) with finite and zero m VBS for the strong and weak bonds, respectively, as shown in the set of Fig. 1(c).The results are presented in Appendix B. With increasing m VBS , the QSHI phase with C s = +1 evolves into a topologically trivial phase with C s = 0 at a finite value of m VBS .And the TPT in this mean-field system is accompanied by single-particle gap closing and reopening.Thus, the appearance of the KVBS phase with C s = −1 in model (1) originates from electron correlations.Theoretically, the nontrivial KVBS phase should change into a trivial KVBS phase if m VBS in the thermodynamic limit is large enough.However, the antiferromagnetic Mott insulator, which breaks the spin U (1) symmetry and time-reversal symmetry, takes over for strong interactions before the KVBS order completely breaks the nontrivial topology.This is consistent with our numerical results that the trivial KVBS phase is absent in the phase diagram.

D. Excitation gaps
Both the single-particle gap and the spin gap for λ/t = 0.05 around the quantum phase transition are presented in Fig. 6.In Appendix C, the raw data of imaginarytime single-particle Green's function and spin-spin correlation function defined in Eq. ( 4) is also presented and discussed.There is a general dip in the curves of both gap when varying parameter U/t for all the system sizes.The location of the dips also changes with increasing system size because of finite-size effect.The spin gap is significantly smaller than the single-particle gap around the ( ) Finite-size scaling of the single-particle gap ∆sp(K)/t and the spin gap L∆s(Γ)/t with the dip values of different system sizes in Fig. 6 for λ/t = 0.05.A quadratic and linear fitting are applied for ∆sp(K)/t and L∆s(Γ)/t, respectively.
transition point, highlighting the fact that the low-energy excitations of this interacting system are indeed bosonic as mentioned in Sec.III B. In Fig. 6(b), the value of spin gap is even larger in the system with larger system size at some specific parameter U/t before the transition.This effect has also been observed in the Kane-Mele-Hubbard model with small SOC parameter 9 .
To reliably extrapolate the excitation gaps at the transition point and in the thermodynamic limit, we carry out the finite-size scaling of ∆ sp /t and L∆ s /t with the dip value of every system size shown in Fig. 6, for which the results are presented in Fig. 7.The extrapolated results show a finite single-particle gap 0.055(5)t and a vanishing spin gap at the transition point in the thermodynamic limit.Similar results are also obtained for λ/t = 0.2 (as presented and discussed in Appendix IV) with larger single-particle gaps at the corresponding transition points.Thus, we can conclude that across the TPT the critical mode is bosonic (spin) instead of fermionic, similar to that of free fermion systems.The finite singleparticle excitation gap is also consistent with the fact that poles of single-particle Green's function at the transition point are absent.

IV. SUMMARY AND DISCUSSION
In this work, we have studied the ground-state properties of the Kane-Mele model with cluster charge interaction in a numerically exact way via large-scale PQMC simulations.We have discovered a topologically nontrivial KVBS phase, as the combination of a QSHI phase with spin Chern number C s = −1 and a long-range KVBS order with spontaneous Z 3 symmetry breaking.We have also identified a quantum phase transition from a QSHI phase with C s = +1 to this KVBS phase.Across this transition, the spin excitation gap shows closing and reopening while the single-particle gap remains finite.We have also observed appearance of two zeros of the singleparticle Green's function at the transition point, which is consistent with the change of the spin Chern number.Our PQMC results of correlation functions for KVBS order suggest the quantum phase transition to be of first order.
Our work in this paper explicitly shows a coexistence of topological symmetry protection and spontaneous symmetry breaking in interacting fermion systems.We also note that the coexistence of intrinsic topological orders and charge orders has recently be studied and confirmed by exact diagonalization 42,43 in correlated fermion systems.Combined into a complete piece, these exotic quantum phases and corresponding phase transitions have greatly extended our knowledge of topological phases of matter as well as their quantum phase transitions.As a future study, it will be interesting to investigate whether a spontaneous symmetry breaking quantum phase transition, which is accompanied by a topological phase transition, can be continuous.Appendix A: Additional QMC results for λ/t = 0.05 and λ/t = 0.2 In the main test, we have shown joint PQMC results for λ/t = 0.05 and λ/t = 0.2 with the ground-state λ-U phase diagram in Fig. 1(c).In this appendix, we present the additional QMC results for λ/t = 0.05 and λ/t = 0.2 which are not shown in the main text.
For λ/t = 0.05 case, the results of KVBS structure factor, the histogram of KVBS order and the spin Chern number are shown in Figs. 8, 9 and 10, respectively.First, from the finite-size scaling of C VBS (K)/L 2 , we can determine the phase transition point U c /t = 2.09(1), as shown in Fig. 8. From the histogram of ∆ K presented in Fig. 9, we can clearly observe the evolution of KVBS order, which is absent for U/t = 2.05 and appears for U/t = 2.15, in consistence with Fig. 8. On the other side, the spin Chern number C s as shown in Fig. 10(a) shows finite-size effect, and it's converging to the value of U c /t obtained from Fig. 8.After the interpolation calculations, we obtain integer quantized values of Cs changing from C s = +1 to C s = −1, suggesting the same TPT as the one for λ/t = 0.2 presented in Fig. 5. Besides, the zeros of single-particle Green's function also appears at the transition point, which is also consistent with λ/t = 0.2 case.
As discussed in the main text, our QMC data suggests first-order quantum phase transition from the disordered QSHI phase to the QSHI phase with long-range KVBS order for λ/t = 0.2, especially the results of histogram of ∆ K .However, we are unfortunately not able to draw such a reliable conclusion of the quantum phase transition for λ/t = 0.05 case.We should note that for λ = 0 case, the quantum phase transition from Dirac semimetal to KVBS phase is continuous 39 .This means that there might exists finite-size scaling crossover behavior 13 when λ evolves from zero to 0.2 (first-order phase transition according to our results).Furthermore, it means that even if the phase transition for λ/t = 0.05 is of first order, significantly larger system sizes (than the ones shown in Fig. 8) are needed to resolve that.However, the results of KVBS structure factor shown in Fig. 8 fail to give a good data collapse, which actually contradicts the hypothesis of continuous phase transition.Theoretically, the results of excitation gaps shown in Fig. 6 and Fig. 7 explicitly shows the finite single-particle gap, which fits the physical picture discussed in Sec.III B and thus also supports a first order phase transition.Finite-size scaling of the single-particle gap ∆sp(K)/t and the spin gap ∆s(Γ)/t with the dip values of different system sizes in Fig. 11 for λ/t = 0.2.Quadratic fittings are applied.The lattice structure and the Brillouin zone of the mean-field Hamiltonian.(a) Each site in an unit cell is marked with an integer and the lattice vectors are also given as â1 and â2.The bold lines stand for the electron hoppings that couple to the order parameter mKVBS.We also draw the three nearest unit cells t R i ,R j which contribute the hopping matrixes T R i ,R j to the Bloch Hmailtonian H k .(b) The reciprocal lattice vectors b1, b2 are given corresponding to the lattice vectors â1 and â2, respectively.The first Brillouin zone is marked by the solid green line.By comparison, We also draw the reciprocal lattice vector and the Brillouin zone of Hamiltonian (1) with B1, B2 and the black dashed line, respectively.
For the λ/t = 0.2 case, the results of excitation gaps in fermion and spin channels as well as their finite-size scalings are shown in Figs.11 and 12. Similar to λ/t = 0.05 case shown in Fig. 6, we can observe that the spin gap is significantly smaller than the single-particle gap.After the extrapolation of the dip values of both gaps to the thermodynamic limit as shown in Fig. 12, the singleparticle gap reaches a value close to 0.2t while the spin gap vanishes at the transition point.Thus, the spin degrees of freedom becomes critical while the fermion channel is gapped across the phase transition for λ/t = 0.2.

FIG. 1 .
FIG. 1.The lattice structure and ground-state λ-U phase diagram.(a) Illustration of the honeycomb lattice and Hamiltonian (1).A unit cell of the lattice, indicated by â1, â2, contains A (red dot) and B (green dot) sublattices.The black lines represents NN hopping, and the red lines shows the SOC with arrows indicating ν ij = +1.(b) Brillouin zone of the model and the reciprocal lattice vector b1 and b2 are given corresponding to â1 and â2, respectively.(c) The λ-U phase diagram from PQMC simulations.Three cases of λ/t = 0.05, 0.2, and 0.4 are studied.The disordered-KVBS ordered phase transition is verified to coincide with the TPT from Cs = +1 to Cs = −1, as shown by the solid, blue line.Thus, a topologically nontrivial KVBS phase is established.As a comparison, the quantum phase transition with Uc/t = 1.666(8) for λ = 0 from Dirac semimetal to a trivial KVBS phase from Ref. 39 is also indicated by the solid, green dot.

FIG. 3 .
FIG. 3. (a)The correlation ratio Rcorr for the KVBS order across the phase transition, and (b) finite-size scaling of structure factors of KVBS order with U/t = 3.28, 3.29, and 3.30 for λ/t = 0.2.

FIG. 5 .
FIG. 5. (a) The spin Chern number Cs, and (b) positive eigenvalue η K of Gσ(iω = 0, K) matrix across the TPT for λ/t = 0.2.In (a), both results of direct QMC calculations in finite-size systems and application of the interpolation technique (see Ref. 24) are shown.
ACKNOWLEDGMENTS Y.Y.H thanks Zi-Xiang Li for valuable discussions.This work was supported by the National Science Foundation of China (Grants No. 11874421 and No. 11774422).R. Q. H. was supported by the Fundamental Research Funds for the Central Universities, and the Research Funds of Renmin University of China (Grant No. 18XNLG11).Computational resources were provided by the Physical Laboratory of High Performance Computing at Renmin University of China and National Supercomputer Center in Guangzhou with Tianhe-2 Supercomputer.The Flatiron Institute is a division of the Simons Foundation.
FIG. 11.(a) Singel-particle and (b) spin gap versus U/t near the phase transition point with λ/t = 0.2.
FIG. 13.The lattice structure and the Brillouin zone of the mean-field Hamiltonian.(a) Each site in an unit cell is marked with an integer and the lattice vectors are also given as â1 and â2.The bold lines stand for the electron hoppings that couple to the order parameter mKVBS.We also draw the three nearest unit cells t R i ,R j which contribute the hopping matrixes T R i ,R j to the Bloch Hmailtonian H k .(b) The reciprocal lattice vectors b1, b2 are given corresponding to the lattice vectors â1 and â2, respectively.The first Brillouin zone is marked by the solid green line.By comparison, We also draw the reciprocal lattice vector and the Brillouin zone of Hamiltonian (1) with B1, B2 and the black dashed line, respectively.