Fermion-sign-free Majarana-quantum-Monte-Carlo studies of quantum critical phenomena of Dirac fermions in two dimensions

Quantum critical phenomena may be qualitatively different when massless Dirac fermions are present at criticality. Using our recently-discovered fermion-sign-free Majorana quantum Monte Carlo (MQMC) method introduced by us in Ref. [1], we investigate the quantum critical phenomena of {\it spinless} Dirac fermions at their charge-density-wave (CDW) phase transitions on the honeycomb lattice having $N_s=2L^2$ sites with largest $L=24$. By finite-size scaling, we accurately obtain critical exponents of this so-called Gross-Neveu chiral-Ising universality class of {\it two} (two-component) Dirac fermions in 2+1D: $\eta=0.45(2)$, $\nu=0.77(3)$, and $\beta=0.60(3)$, which are qualitatively different from the mean-field results but are reasonably close to the ones obtained from renormalization group calculations.


Introduction
Quantum phase transitions have attracted enduring interests in condensed matter for many decades [2]. For contiguous phase transitions, scaling invariance and universality classes were introduced as fruitful concepts to understand them [3]. Renormalization group (RG) provides us a modern and powerful framework to characterize critical behavior and university class [4,5]. When gapless fermions are present and coupled with order parameter at criticality, quantum critical phenomena may be dramatically changed; characterizing them are often challenging due to strong couplings between gapless fermions and order parameter. Such fermionic quantum critical points include quantum spin-density-wave/charge-density-wave (CDW)/nematic transitions in metals with Fermi surface or in semimetals with massless Dirac fermions.
The smallest number of two-component Dirac fermions realizable in 2D lattice models is N = 2, a prototype example of which is from spinless fermions on the honeycomb lattice (note that N = 1 massless Dirac fermions in 2+1D can only be realized on surfaces of 3D topological insulators). In the presence of strong repulsion between fermions on nearest-neighbor (NN) sites, the ground state develops a finite CDW. Close to transition, chiral Dirac fermions strongly couple with Ising order parameter, rendering this chiral-Ising transition in 2 + 1D qualitatively different from the usual Ising transition in 2 + 1D. Recently, a remarkable continuous-time (CT) QMC study [31] of the CDW transition in the honeycomb model of N L 2 s 2 = sites with L up to 15 obtained 0.30(2) η = which is in some discrepancy with previous RG results ( 0.50 0.64 η ≈ ∼ ) [18][19][20][21]. We recently discovered an auxiliary-field QMC method employing Majorana representation which can solve the fermion-sign problem in such spinless fermion models [1]. Using this sign-free Majorana quantum Monte Carlo (MQMC), we performed highly-accurate simulations of spinless fermions on the honeycomb lattice of N L 2 s 2 = sites with largest L = 24 and obtained critical exponents of the N = 2 chiral-Ising transition in 2 + 1D: . Note that we would obtain a much smaller 0.32(2) η = which is consistent from previous CT-QMC simulations [31], if we were using simulation results on lattices only up to L = 15. Using MQMC, we also investigated the CDW transition in the π-flux square lattice model with NN repulsions on lattices of 2L 2 sites with largest L = 24 and obtained similar critical exponents. Our results are in reasonable agreement with the existing RG calculations of the Gross-Neveu-Yukawa theory: 0.50 0.64 η ≈ ∼ and 0.74 0.93 ν ≈ ∼ [13,18,20,21] even though there are still slight discrepancies. Our numerically-exact results may serve as a possible benchmark for higher-order RG calculations or other theoretical analysis in the future.

Models
To investigate the universal quantum critical behavior of the N = 2 chiral-Ising transition in 2 + 1D, we study the following simple interacting spinless fermion model on the honeycomb lattice: where c i † creates a fermion on site i, t t ij = is the hopping of fermions between NN sites, and V is the density interaction between fermions on NN sites. Here the hopping and interaction parts are labeled by H 0 and H int , respectively. Hereafter we set t = 1 as the unit of energy.
Noninteracting fermions on the honeycomb lattice features massless Dirac dispersion around two inequivalent Dirac points K ± ⃗ . Because of vanishing density of states, such N = 2 two-component mass Dirac fermions are stable against any weak interactions which preserving symmetries of the model but can be gapped when interactions are sufficiently strong. In this paper, we shall focus on the CDW transition induced by relatively-strong NN repulsions. Note that quantum anomalous Hall (QAH) ordering will be induced by strong next-nearest-neighbor (NNN) repulsions [35] and pair-density-wave ordering will be induced by relativelystrong NN attractions [36], both of which are not the focus of present numerical studies.

The MQMC method
As an intrinsically-unbiased numerical method, QMC [22][23][24][25][26] is an important and valuable tool to study quantum phases of matter and critical behaviors at phase transitions when it is free from the fermion-sign problem. To access the quantum critical behavior of the CDW transition induced by repulsions, it is desired to simulate the interacting model using fermion-sign-free QMC on lattices with relatively-large size. Nonetheless, conventional auxiliary-field QMC simulations of interacting spinless fermion model generally suffer from the fermion-sign problem because the usual strategies [37][38][39][40][41][42][43][44][45][46] employed in auxiliary-field QMC for solving fermion-sign problem of even species of fermions do not directly work. Recently, we discovered that certain interacting spinless fermion systems can be studied by fermion-sign-free QMC employing Majorana representation [1]. Upon observing that each complex fermion can be represented as two Majorana fermions, we perform Hubbard-Stratonovich (HS) transformations to decouple interactions in terms of bilinear Majorana fermions. Under certain conditions such as particle-hole symmetry, a symmetric treatment of two species of Majorana fermions is possible such that the Boltzmann weight is a product of two identical real quantities and is then positive definite.
To the best of our knowledge, MQMC is the first QMC approach based on auxiliary fields to solve fermion sign problem in a class of interacting spinless fermion models [1]. We emphasize that the auxiliary-field MQMC approach proposed here is qualitatively different from the sign-free continuum-time QMC method [47][48][49][50][51]. One important advantage of the fermion-sign-free MQMC is its high-efficiency of directly simulating interacting models at zero temperature by using projector. To be self-contained, here we briefly review the fermion-sign-free MQMC algorithm introduced by us in [1].
where n labels the discrete imaginary time, N Δτ β = τ , and the approximation is good for small enough Δτ. We observe that it will be advantageous to write complex fermion operators in Majorana representation: , which enable us to rewrite equation (1) as follows: Explicitly, we obtain the positive-definite Boltzmann weight or 2 giving rise to the same result. Because we are interested in the quantum phase transition at zero temperature, we can use projector [52][53][54] MQMC to explore the ground state properties of the system. In the projector QMC, a trial wave function T ψ | 〉is introduced and after projection the ground state 0 ψ | 〉can be obtained. Consequently, the expectation values of observable Ô is: In our simulations, we use an Slater-determinant wave function as the trial wave function T ψ . Similarly, we can prove that projector MQMC is also free from fermion-sign-problem in the parameter region where V 0 > , as shown in details in [1].

The CDW transition on the honeycomb lattice
We have performed fermion-sign-free MQMC simulations to study the spinless fermion model on the honeycomb lattice with NN repulsive interaction V. When V V c > , the ground state develops a finite CDW ordering which breaks the inversion symmetry of the model. To find the critical interaction V c of the CDW phase transition, we compute the CDW structure factor on a finite lattice using MQMC: V V c = , negative a 0 is generically obtained, which indicates that this fitting could somewhat overestimate V c when L is not large enough.
To obtain V c more accurately, we further compute the quartic of CDW order parameter M 4 defined as M N n n n n 1 2 According to the following scaling functions for M 2 and M 4 : where we have implicitly assumed the dynamical critical exponent z = 1 for the CDW transition in the current model, the Binder ratios of sufficiently large L should cross at V V c = . The calculated Binder ratios for different V and different L are shown in figure 1(b), which clearly show that V 1.355 c ≈ , which is consistent with the one obtained in [31].
After obtaining V c for the CDW transition, we can further compute the independent critical exponents η and ν which are critical exponents regarding correlation functions. Other critical exponents such as β may be obtained from η and ν through hyper-scaling relations. We use two different methods to obtain η. First, we perform a finite-size scaling analysis of M 2 according to the scaling function: M L L ( ) 2 1 Second, we compute the density-density correlation C L n n ( ) where r L L ( 2, 2) max ⃗ = is the largest possible separation between two sites in the lattice with periodic boundary conditions. At criticality, this correlation decays in power-law as C L L ( ) 1

Conclusions and remarks
Using the fermion-sign-free MQMC method, we have simulated the spinless fermion model on the honeycomb lattice as well as on the π-flux square lattice with NN repulsions. We focused on the quantum critical behavior at the CDW transition in these models, which is in the universality class of the chiral-Ising transition of N = 2 twocomponent mass Dirac fermion in 2 + 1D. The low energy physics of this CDW transition can be described by the N = 2 Gross-Neveu or Gross-Neveu-Yukawa theory in 2 + 1D. We numerically showed that 0.45(2) η = , 0.77(3) ν = , and 0.60(3) β = at this N = 2 chiral-Ising transition in 2 + 1D. The critical exponents obtained by the fermion-sign-free MQMC simulations are reasonably consistent with the ones obtained by the two-loop RG calculations in D 4 ϵ = − expansion even though there is still some slight discrepancy. The discrepancy should not come from the differences in the models used in numerical simulations and in RG calculations. In the RG calculations, the two massless Dirac fermions have the same chirality, which means that the symmetry breaking phase in its corresponding lattice model is the QAH state. The two Dirac fermions on the honeycomb lattice have opposite chirality and the broken symmetry phase is a CDW state. Nonetheless, the RG equations of the coupling constants in the two corresponding low-energy field theories describing the QAH and CDW transitions are identical. Consequently, the two seemly-different transitions should be in the same universality class, at least, in the sense of the critical exponents even though the two broken symmetry phases are topologically distinct. We believe that the critical exponents obtained in the present MQMC simulations could serve as a benchmark for more accurate higher-loop RG calculations in the ϵexpansions which are deferred to future studies.