Designer Monte Carlo Simulation for Gross-Neveu Transition
DDesigner Monte Carlo Simulation for Gross-Neveu-Yukawa Transition
Yuzhi Liu,
1, 2
Wei Wang,
1, 2
Kai Sun, and Zi Yang Meng
4, 1, 51
Beijing National Laboratory for Condensed Matter Physics and Institute of Physics,Chinese Academy of Sciences, Beijing 100190, China School of Physical Sciences, University of Chinese Academy of Sciences, Beijing 100190, China Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA Department of Physics and HKU-UCAS Joint Institute of Theoretical and Computational Physics,The University of Hong Kong, Pokfulam Road, Hong Kong, China Songshan Lake Materials Laboratory, Dongguan, Guangdong 523808, China (Dated: February 28, 2020)In this manuscript, we study quantum criticality of Dirac fermions via large-scale numerical simulations,focusing on the Gross-Neveu-Yukawa (GNY) chiral-Ising quantum critical point with critical bosonic modescoupled with Dirac fermions. We show that finite-size effects at this quantum critical point can be efficientlyminimized via model design, which maximizes the ultraviolet cutoff and at the same time places the bare controlparameters closer to the nontrivial fixed point to better expose the critical region. Combined with the efficientself-learning quantum Monte Carlo algorithm, which enables non-local update of the bosonic field, we find thatmoderately-large system size (up to 16 ×
16) is already sufficient to produce robust scaling behavior and criticalexponents. The conductance of free Dirac fermions is also calculated and its frequency dependence is found tobe consistent with the scaling behavior predicted by the conformal field theory. The methods and model-designprinciples developed for this study can be generalized to other fermionic QCPs, and thus provide a promisingdirection for controlled studies of strongly-correlated itinerant systems.
I. INTRODUCTION
Quantum critical points (QCP) in the presence of fermionicfluctuations are among the most intriguing topics in the studyof strongly correlated systems, which pave the way towardsnew paradigms of quantum matter beyond the conventionalFermi-liquid theory of metal and the Landau-Ginzburg-Wilsontheory of phase transitions. Due to the strong coupling nature,to understand these fermionic QCPs remains a highly chal-lenging task, and many fundamental questions remain open.For example, after four decades of intensive studies, quan-tum criticality in the presence of Fermi surfaces still remainsa highly active research subject, from the early effort of theHertz-Mills-Moriya framework [1–3] based on the leading-order approximation to the more recent studies that reveal thecrucial and nontrivial contributions of higher order terms [4–7]. Even today, the fate of such fermionic QCPs still remainsan important open question and a source for exciting new ideasand insights [7–10].The situation at the numeric front is equally, if not more,challenging, as the divergent length scale and the fermionic na-ture of this problem prevents us from obtaining unbiased andreliable numerical solutions at the thermodynamic limit. Forexample, in contrast to bosonic QCPs, many of which can be ef-ficiently simulated via quantum Monte Carlo (QMC) method-ologies [11], the notorious sign problem makes it highly chal-lenging to access the same level of knowledge in the presenceof fermionic fluctuations. Only till very recent, new pathwaystowards large-scale unbiased simulations of fermoinc QCPsbecome available, where utilizing designers models, variousfermoinc QCPs become accessible to sign-problem-free QMCsimulations, such as nematic [12], ferromagnetic [13] andantiferromagnetic QCPs [14–17] (See Ref. [18] for a recentreview). From these efforts, it becomes possible to obtainaccurate and reliable information about the scaling behaviors in the close vicinity of these QCPs, to test and improve ourtheoretical knowledge about these challenging problems.In this manuscript, we study fermionic QCPs in strongly-correlated Dirac systems utilizing Gross-Neveu-Yukawa typeof models, where interactions between Dirac fermions are me-diated by a bosonic field. These QCPs are often refereed to asthe Gross-Neveu (GN) or Gross-Neveu-Yukawa (GNY) tran-sitions, which, as will be shown below, can be simulated viasign-problem-free QMC methods. For theoretical treatment,due to the absence of Fermi surfaces and the emergent Lorentzand conformal symmetry, some of the GNY QCPs are be-lieved to be among the “simplest” fermionic QCPs, and thushold the highest expectation for achieving agreement betweencontrolled analytical calculations and numerical simulations.In recent years, such efforts have been attempted from both nu-merical and theoretical sides, e.g., the chiral-Ising GNY tran-sitions have been studied via high-order expansions [19, 20]and quantum Monte Carlo simulations [21–24].In the numerical studies of GNY transitions and its quantumcriticality, the most significant challenge lies in the finite-sizeeffect. A typical fermionic quantum Monte Carlo requires acomplexity of O ( β N ) where N = L d is the lattice size with L being the linear span of the d dimensional system and β theinverse temperature. In addition to the typical challenge of crit-ical slowing down, quantum critical points in the GNY modelshave two additional sources of finite-size effects. First, in a lat-tice model, Dirac fermions emerge in the close vicinity of theDirac points. As far as critical phenomena are concerned, thepart of the energy band near the Dirac point (with linear disper-sion) contributes to the correct infrared (IR) scaling behavior,while the rest part of the energy band, which deviates fromthe Dirac linear dispersion, gives short-distance non-universalphysics. In other words, the scaling behavior comes from onlya portion of the Brillouin zone (less than 15% in a typical latticemodel as shown below), which effectively reduces the system a r X i v : . [ c ond - m a t . s t a t - m ec h ] F e b size by ∼
85% in numerical simulations and thus significantlyamplifies finite-size effects. To enlarge the linear-dispersionregion, special models are utilized, such as the SLAC fermionconstruction in the recent work of Ref. 23. However, such spe-cial models often involve infinite-long-range hopping, whichmakes the system non-local. Secondly, in a typically latticemodel with bosons and fermions coupled together, these twotypes of particles in general have different velocities. This ve-locity difference explicitly breaks the Lorentz symmetry, andresults in a very slow (logarithmic) RG flow of the velocitiestowards the fixed point [25]. Such a slow RG flow requiresextremely large system sizes to reach the vicinity of the fixedpoint, making it very challenge to produce the correct expo-nents on finite size lattice.In the study of fermionic quantum criticality, to overcomethe critical slowing down, various efforts have been made to ac-cess larger system sizes, such as the self-learning Monte Carlo(SLMC) [26, 27], which provides more efficient updates, andits momentum space extension (elective momentum ultra-sizeMonte Carlo EQMC) [16]. These improvements addressesthe numerical aspect, and are implemented a posteriori tothe model design. In this manuscript, we utilize a differentapproach to reduce the finite-size effect through a priori opti-mizing the model design to better utilize the emergent Lorentzsymmetry at the chiral Ising GNY transition. In particular,we develop designer model Hamiltonians with two proper-ties: (1) the Dirac linear-dispersion region occupies a largerportion of the Brillouin zone without sacrificing locality and(2) the bosons and fermions have the same velocity to avoidthe slow RG flow [25]. These efforts minimize the finite-sizeeffects mentioned above. Combined with advanced numericscheme such as SLMC, we show explicitly that finite-size ef-fects are efficiently suppressed in the obtained results, givenrise to robust critical exponents. Furthermore, we find ourdesigner Hamiltonian provide robust results in the behavior offree Dirac fermion conductivity at finite Matsubara frequency,such consistency gives one the confidence of obtaining thesimilar level analysis in the GNY QCP in the future works.The rest of the paper is organized as follows: in Sec. IIthe designer model of the N f = II. DESIGNER MODEL
The designer model is realized by coupling Dirac fermionswith dynamical bosonic field [21, 22, 28]. The model is definedon a square lattice shown in Fig. 1, with two bosonic sites andtwo fermionic sites per unit cell. The bosons acquire an Ising symmetry and the Lagrangian of the system is given as L = L Boson + L Fermion + L Coupling (1)with L Boson = (cid:213) p (cid:20) (cid:0) ∂φ p ∂τ (cid:1) + m φ p + φ p (cid:21) + (cid:213) ( p , q ) J pq ( φ p − φ q ) , L Fermion = (cid:213) ( i , j ) ,σ ψ † i ,σ [( i ∂ τ − µ ) δ ij − t ij e i σθ ij ] ψ j ,σ + h . c ., L Coupling = (cid:213) (cid:104)(cid:104) i , j (cid:105)(cid:105) ,σ λ ij φ p ψ † i ,σ ψ j ,σ + h . c ., (2)where ψ i ,σ ( ψ † i ,σ ) is the fermionic annihilation (creation) op-erator at the fermionic site i with spin σ = ↑ or ↓ , and φ p represents the scalar bosonic field at the bosonic site p . µ ischemical potential and we set µ = θ ij . Forthe nearest-neighbor hopping between different sublattices t and the third-nearest-neighbor hopping between different sub-lattices t , the phase factor is set to θ = π /
4, which introduces a π -magnetic-flux for each plaquette of the fermion lattice. Thistight-binding model gives rise to two Dirac cones at ( π, ) and ( , π ) in the Brillouin zone (BZ). The conventional π -fluxmodel with only t is depicted in Fig. 1 (a), which is the tight-binding model employed in previous studies [21, 22]. Themodel utilized in this work is shown in Fig. 1 (b), where oneadditional hopping term t is introduced with t = t /
27. Asshown in the right panels of Fig. 1 (a) and (b), this additionalhopping term greatly increases the area of the linear dispersionregion, denoted by the green contour lines, from 13.7% in (a)to 39.8% in (b) of the entire BZ. Here, this green contour marksmomentum points, at which the energy band deviates from theideal linear dispersion by 5%. It serves as the ultraviolet cutoff( Λ ) in the renormalization group (RG) analysis of this Diracsystem.For a critical system, the RG flow is governed by the dimen-sionless quantity L Λ , where L is the system size. A larger Λ effectively increases the system size, and thus pushes thesystem closer to the thermodynamic limit. Similar ideas ofincreasing the linear dispersion area has been used in anotherrecent QMC study [23], where the linear dispersion region isexpanded to the entire BZ via introducing infinite-long-rangehopping terms, known as the SLAC fermion model. Theselong-range hopping terms render the model non-local. In con-trast, in the model that we adopted here, although the linear-dispersion region doesn’t cover the entire BZ, the Hamiltonianis local with only short-range couplings.Fig. 1 (c) depicts the full model with circles representingthe fermionic sites and squares the bosonic ones. Similar tothe fermions, the bosons here also have extended, but local,interactions, denoted as J , J , J and J from the first-neighborto the fourth-neighbor couplings respectively. The values ofthese coupling strengths are shown in the caption of Fig. 1. Atthe QCP, the critical boson mode exhibits a linear dispersion atsmall wavevector. These values of the J ’s maximize the area E(x,y)
E(x,y) (-π,0) (-π/2,0) (0,0) (π/2,0) (π,0) (-2π,0) (-3π/2,0) (-π,0) (-π/2,0) (0,0)
FIG. 1. The fermion tight-binding models are shown in (a) and (b), with (a) only contains the nearest-neighbor hopping t (Ref. [21]) and(b) contain both the nearest-neighbor hopping t and the third-nearest-neighbor hopping between different sublattices t (this work). Thegreen dashed lines in the right panels mark the linear dispersion region, out of which the band structure starts to deviate from the Diraclinear dispersion. The t hopping in Fig. (b) greatly increases the area of the linear dispersion region. The black arrows in the left panelrepresent the phase of the complex hopping strength θ = π , and we set t = t = t . The value of t is obtained from the optimizationprocess of trying to enlarge the linear dispersion area in BZ. The lattice model is depicted in (c), where the coupling constants for bosons are J = − J , J = J , J = − J . The values of boson couplings are also obtained from the parameter optimization process. The bosons andfermions are coupled together via the λ i , j = ± λ i , j alternates from one plaque to another, as marked by in (c). Overall,one unit cell of this model contains two fermion sites and two boson sites. (d) and (e) show dispersions of fermions and bosons respectively(see Appendix C for detail explanation of how to obtain these dispersions). Here, we presented both the bare dispersions (ignoring interactions)and the renormalized ones at the chiral Ising transition via QMC simulations (with L = of this linear-dispersion region, which helps to suppress finite-size effects same as what we did for the fermions above. Withthis parameter choice, the bare boson velocity at the criticalpoint of L Boson is v b = (cid:112) J /
8, and the Fermi velocity atthe Dirac point is v f = √ t . Here, we set t = J such that bosons and fermions shareidentical velocity v b = v f (See Appendix A for details). Thisidentical velocity results in an emergent Lorentz symmetry,where v b = v f serves as the speed for light. For a Lorentz-invariant system, it is known that the speed of light would notflow under RG, even in the presence of strong interactions [29].This a priori setup avoids the slow logarithmic RG flow of thevelocity mentioned above [25].The couplings between bosons and fermions are denoted as λ ij , which couples to a pair of next-nearest-neighbor fermionsites with the boson site between them. The coupling constant λ ij takes the value of ±
1, where the two bosonic sublatticestake opposite signs. In the ordered phase (i.e. the boson fieldhas a nonzero expectation value), this alternating sign givesrise to the required Berry phase to the fermions, and thus gapsout the Dirac points, transforming the Dirac semi-metal intoa dynamically-generated quantum-spin-Hall insulator. Moredetails about such a topological phase transition and numericalrealization can be found in Ref. [21].Fig. 1 (d) and (e) depict the dispersions near the fermionicDirac point at ( π, ) and bosonic dispersion at ( , ) at the QCP.The dispersions are obtained from the fitting of the imaginarytime fermionic and bosonic Green’s functions, and are ex-plained in Appendix C. The green and brown lines in Fig. 1(d) are the bare fermion dispersions for the models shown inFig. 1(a) and (b), respectively. One can see that the linear- dispersion region for the brown line is indeed larger than thatof the green line. The blue stars are the QMC results obtainedat the chiral Ising critical point. Fig. 1 (e) presents the bosonicdispersion at the bare Ising critical point (the brown line) andthe chiral Ising critical point (the blue stars). It is worthwhileto highlight that for both fermions and bosons, their veloci-ties receive little renormalization. This absence of velocityrenormalization is because we set the bosons and fermions tohave the same bare velocity as mentioned above (we think thedeviation of the boson dispersion at the QCP from the brownline still comes from the inevitable finite size effect of L = III. RESULTS
In the simulation, we vary the value of the mass term m in L Boson in Eq. (1). Without fermions, tuning this mass termtriggers a ( + ) d Ising transition for the φ field. In the presenceof the fermions, this transition induces a fundamental change inthe fermion band structure. In the disordered (ordered) phase, (cid:104) φ (cid:105) = (cid:104) φ (cid:105) (cid:44) -0.2 -0.17 -0.14 -0.11 -0.08 -0.05 -0.02 -0.2 -0.17 -0.14 -0.11 -0.08 -0.05 -0.02 FIG. 2. (a) The Binder ratio U ( L , m ) for different system sizes L across the transition of m . Curves for two consecutive system sizes L and L + m c ( L ) . (b)The correlation ratio R ( L , m ) for different system sizes L across thetransition of m . Curves for two consecutive sizes L and L + m c ( L ) . The finite-size analysis forthese crossing points are shown in Fig. 3. the Dirac fermions. In addition, it can be shown that thisDirac mass carry a nontrivial topological index, and resultsin a dynamically generated quantum spin Hall insulator [21].Furthermore, the interplay between bosons and fermions alsochanges the universality of the quantum critical point, fromIsing to the GNY chiral Ising, which is believed to be amongthe simplest fermionic QCPs where possibly analytical andnumerical approaches might be able to give consistent criticalexponents to build a concrete CFT description [19, 24, 32, 33].For the model described above, the Dirac fermions acquire N f = N f = -0.20-0.16-0.12-0.08 -15 -10 -5 FIG. 3. (a) Crossing point analysis for the Binder and correlationratios obtained in Fig. 2. The two m c ( L ) curves extrapolate to thesame critical point m ∗ c = − . ( ) at the thermodynamic limit (notethat we set the two curves to meet at a single value the m ∗ c in thefitting) and the two powers in the extrapolation according to Eq. (6)consistently reveal 1 / ν GN + ω = . ( ) , with different non-universalcoefficients a . (b) Crossing point analysis of the Binder ratio U c ( L ) and the correlation ratio R c ( L ) according to Eqs. (7) and (8). Theobtained correction exponent ω = . ( ) . Combining the exponentsfrom (a) and (b), we find 1 / ν GN = . ( ) . (c) Data collapsing for theorder parameter measured at different m and system sizes L . Here, weutilized the the exponent ν GN and the critical point m ∗ c obtained aboveand the system sizes are L = , , , , ,
16. A good collapse isachieved with dynamic critical exponent z = η φ = . ( ) . A. Crossing-point analysis
First we determine the precise location of the quantum crit-ical point via the Binder ratio ( U ) [34] and correlation ratio( R ) [35] of the bosonic order parameter m φ = N (cid:213) (cid:104) p , q (cid:105) φ p φ q (3)where N is the number of boson lattice sites, and the twodimensionless ratios are given as U = (cid:32) − (cid:104) m φ (cid:105)(cid:104) m φ (cid:105) (cid:33) (4) R = − S ( Q + δ q ) S ( Q ) (5)here S ( q ) = N (cid:205) (cid:104) p , q (cid:105) e − i q ·( r p − r q ) φ p φ q is the structure factorof the bosonic correlation function and r p and r q stand forthe position of the bosonic field on the lattice in Fig. 1 (c).The ordering wavevector in our case is Q = ( , ) , and δ q = ( , π / L ) or ( π / L , ) is the smallest momentum away from Q on a finite lattice with linear span L .The results are shown in Fig. 2 (a) and (b). The finite-size transition point m c ( L ) is the crossing points of the twoconsecutive system sizes, and it is clear from the two panels ofFig. 2 that the position of the m c ( L ) -s are gradually drifting as L increases. According to the finite size scaling analysis [36, 37],the drift of both the critical point m c ( L ) and the corresponding R c ( L ) and U c ( L ) obey the following scaling behavior at nearthe critical point, m c ( L ) = m ∗ c + aL −( / ν GN + ω ) , (6) U c ( L ) = U ∗ c + bL − / ω , (7) R c ( L ) = R ∗ c + cL − / ω , (8)where ν GN and ω are the correlation length exponent and thecorrection expoent of the chiral Ising GNY transition. m ∗ c , U ∗ c and R ∗ c are critical point, correlation ratio and Binder ratioin the thermodynamic limit and a , b and c are non-universalfitting coefficients.As shown in Fig. 3(a), the crossing points of different pairs ofsystem size ( L , L + ) drift towards m ∗ c as a power-law functionas shown in Eq. (6). It is important to emphasize here that forboth the Binder ratio U and the correlation ratio R , the scalinganalysis generates the same critical point m ∗ c = − . ( ) and the same critical exponents 1 / ν GN + ω = . ( ) . Thisconsistency indicates that finite-size effects are under goodcontrol in the current simulation, even though our system sizes(upto L =
16) are just moderately large.In Fig. 3(b), we fit the two curves according to Eqs. (7) and(8), to trace the universal thermodynamic values of the Binderand correlation ratios U ∗ c and R ∗ c , as well as the correctionexponent ω . It turns out that both the coefficient b and c arepositive numbers (in general, their signs depend on the modeland could be opposite in certain systems as shown in Ref. [37]),such that both curves decreases as L increases. Interestingly,the two independent fits (for U and R ) yields the same sameexponent value ω = . ( ) . With the value of ω obtained,we can now determine 1 / ν GN = . ( ) from the combinedexponents in Fig. 3 (a).With m ∗ c , ν G N , we can directly collapse the bosonic or-der parameter with the universal scaling form (cid:104) m φ (cid:105) L z + η = f ( L / ν GN ( m − m ∗ c )/ m ∗ c ) . The results are shown in Fig. 3 (c),with z = ( + ) D chiral Ising universality, the dataof L = , , · · · ,
16 give a very good collapse and from which FIG. 4. The imaginary time decay of the fermion Greens’s function atthe chiral Ising QCP. From the power-law decay of G f ( τ ) ∼ / τ + η ψ we obtain the anomalous dimension of the Dirac fermions η ψ = . ( ) . we read the bosonic anomalous dimension η φ = . ( ) at thechiral Ising GNY QCP.The fermionic anomalous dimension η ψ , however, cannotbe obtained from the above bosonic crossing analysis but haveto be extracted from the fermion Green’s function at the QCP.For a QCP with the dynamic critical exponent z =
1, thereare two ways to obtain this anomalous dimension, by fittingeither the time dependence or the spatial dependence of theGreen’s function, and they should produce the same value.But since in the lattice simulation the spatial and time axesare actually anisotropic, one has more distance along the time-axis than that of the spatial-axis. Moreover, since the Lorentzsymmetry in a lattice model is an emergent symmetry insteadof an explicit one, time and space are not fully identical in ourmodel. In our previous work [21], with only nearest-neighborfermion hopping, we could not obtain power-law decay in real-space of the fermion Green’s function at the finite size QCP. Inthis work, although the results have been greatly improved, westill see that for the fermionic anomalous dimension, the time-dependence data has more data points and exhibiting betterpower-law behavior than the space-dependence data, and wethus perform a more accurate fitting with the time-dependencedata of the fermion Green’s function.In Fig. 4, we present the imaginary time decay of the fermionGreen’s function G f ( τ ) = N (cid:205) Ni = (cid:104) ψ † i ,σ ( τ ) ψ i ,σ ( )(cid:105) at the QCP m = m ∗ c and fit the data with its scaling form G f ( τ ) ∼ / τ + η ψ , (9)where the power 2 is the bare-scaling dimension for free Diracfermions in ( + ) D, and η ψ accounts for the anomalous partof the exponent at the nontrivial GNY fixed point. In Fig. 4,we plot G f ( τ ) at the QCP ( m = m ∗ c ) for several system sizesin the log-log scale. It is clear that at large imaginary time,the power-law decay of the fermion Green’s function manifestswith an anomalous dimension η ψ = . ( ) . / ν GN η φ η ψ ω This work 1.0(1) 0.59(2) 0.05(2) 0.8(1)previous QMC [21] 1.2(1) 0.65(3)previous QMC [22] 1.20(1) 0.62(1) 0.38(1)perturbative RG 4- (cid:15) loop, Padé [19] 0.931 0.7079 0.0539 0.794perturbative RG 4-loop, interpolation betwen 2+ (cid:15) and 4- (cid:15) [20] 1.004 0.735 0.042bootstrap estimation [32] 0.88 0.742 0.044TABLE I. Comparison of obtained critical exponents for the N f = ν G N is the correlation length exponent, η φ is theanomalous dimension of the bosonic field, η ψ is the anomalous dimension of the fermionic field and ω is the correction exponent in the bosonicsector. This completes our analysis of the critical exponents of the N f = v f = v b a priori and as will be shown in the next section, we find that this setupavoided the slow RG flow of velocity, such that the finite-sizesimulations are closer to the nontrivial fixed point of the QCP.Second, the 1 / ν GN , η φ in Tab. I are consistent with previousQMC results, and because finite-size effects have been effi-ciently supressed in the present study, although our systemsizes are only moderately large ( L = / ν GN , η φ andthe correction exponent ω , where the later cannot be accessedin the previous QMC simulations. Our fermionic anoma-lous dimension η ψ is different from the the previous QMCresult [22], and is significantly closer to the theory predictionsfrom perturbative RG and conformal bootstrap. In our previ-ous QMC work [21], the η ψ couldn’t be obtained through thefitting formula shown in Eq. (9), precisely because the strongfinite-size effect in the simulation renders the slow (logarith-mic) RG flow and consequently the Green’s function data atfinite sizes cannot produce a robust value of η ψ . Therefore weare more confident with η ψ obtained in this work. In fact, theexponents reported here, { / ν GN , η φ , η ψ , ω } offer the mostcomplete set of numerically measured critical exponents tillnow for the N f = B. Velocities at transition point
As have been mentioned above, our designer Hamiltonianenjoys the advantage of maintaining the Lorentz symmetryeven at the bare level, i.e., the linear dispersion region of bothfermion and boson are enlarged in the BZ and their velocitiesare made identical by choosing suitable values of t and J .With such optimization, finite-size effects at the chiral IsingGNY QCP are reduced and thus the exponents converges tostable values even with just moderately large system sizesup to L =
16. To further verify whether the fermion andboson velocities indeed remain identical and maintain theirbare values, here we directly measure these two velocity at thechiral Ising GNY QCP.Here, we measure the velocities via the dynamical Green’sfunctions (the dynamical fermionic and bosonic propagators) G f ( k , τ ) = (cid:104) c † ( k , τ ) c ( k , )(cid:105) and G b ( q , τ ) = (cid:104) φ ( q , τ ) φ (− q , )(cid:105) .In the projector QMC, these two Green’s functions are readilyobtained from the imaginary time correlation functions [30,31]. Since both fermion and boson have two sublattices, wemeasure the trace of the 2 × G ( k , τ ) ∝ e − ∆ ( k ) τ , the many-body excitation gaps are obtained at eachmomenta k and q . This relation between the excitation gapand the momentum marks the “mass shell”, in analogy to therelativistic quantum field theory. Then we fit this relation nearthe Dirac point K = ( π, ) for fermions and near the Γ point Q = ( , ) for bosons as ∆ f ∼ v f ( k − K ) and ∆ b ∼ v b ( k − Q ) respectively. This is how the velocities at the strong-couplingfixed point are obtained(see Appendix C for detail calculation).The results are shown in Figs. 1 (d) and (e). With system sizeup to L =
14, we find that v f is very close to its bare value withvery good linear behavior and v b is also close to its bare valuewith slightly larger error-bars and deviations. This is due to thefact that our boson is a continuous field and the Monte Carlomoves, even with non-local SLMC update scheme, still suffersfrom the slow Monte Carlo dynamics and rendering slightlyworse data quality. But, nevertheless, As shown in Fig. 1 (d)and (e), v f does’t flow and remains a constant, in agreementwith theory expectation. The boson dispersion suffers morefrom finite-size effect and slow QMC updates as shown inFig. 1 (e). Within the uncertainty of finite-size effects, it isbelieved to consist with theory expectation. However, to verifythe theory prediction at the quantitative level, larger systemssize are needed to overcome the finite-size effect and achievethe thermodynamic limit. This result confirms the absence ofthe slow velocity RG flow, which helps control the finite-sizeeffect in our simulations. FIG. 5. The frequency dependence of the current-current correlationwith L =
50 and T = /
200 for free Dirac fermions, where thered (blue) data points are obtained from the designer model with t = /
27 (the conventional model with t = σ ∞ = at the limit of ω n (cid:29) T . The inset shows thefitted σ ∞ as a function of the system size L . As L increases, σ ∞ fromthe designer model converges to the exact value (the dash line) muchfaster than that of the conventional model. C. Dirac fermion conductivity
With the basic scaling properties of the N f = ( + ) D quantum rotormodels at the Wilson-Fisher fixed point [38–41] that the quan-tum critical conductivity at finite Matsubara frequency indeedobeys the scaling law predicted by the conformal field theory.However, similar results in fermionic QCPs are still absent.From previous studies of the ( + ) D O(2) transition, itis known that very large L and β are needed to observe the correct scaling behavior [38–41] σ ( i ω n ) σ Q = − ω n (cid:104) J x ( ω n ) J x (− ω n )(cid:105) + · · · = σ ∞ + b (cid:0) T ω n (cid:1) − / ν + b (cid:0) T ω n (cid:1) + · · · , (10)at ω n (cid:29) T where σ ( i ω n ) is the conductivity at the Matsubara frequency ω n = π nT , and σ Q = ( e ∗ ) (cid:126) is the quantum unit of conductancewith e ∗ the effective charge in the model such that σ / σ Q becomes dimensionless. To simplify the discussion, we willset σ Q = J x ( ω n ) in Eq. (10) is the Fourier transformation of J x ( r i , τ ) = − it (cid:104)( e − i σθ c † ( r i , τ ) c ( r i + ˆ x , τ ) − e i σθ c † ( r i + ˆ x , τ ) c ( r i , τ ))(cid:105) whereˆ x denotes the current along the x -direction of the lattice. Wenote that here the ˆ x and ˆ y directions are equivalent due to thefour-fold rotational symmetry, and thus can be symmetrized inthe simulation to increase the data quality. σ ∞ is the limitingvalue of conductivity at T → b and b are coefficient relatedwith the operator product expansion of the current operatorsin terms of other operators of the CFT and can in principle beestimated in the large N limit. Coefficients at bosonic O ( N ) Wilson-Fisher CFT have been estimated at large N [38] andalso recent conformal bootstrap results for N = L =
16, we will restrain ourselves from thecomputation of σ ( i ω n ) at the GNY chiral Ising QCP and leaveit for a separated project. Here we will only demonstrate thecomputation for free Dirac fermions in the absence of inter-actions, i.e., bare fermions at the tree level, such that large L and β can be accessed. Because this limit is well understoodtheoretically, we can compare the numerical results with thevalues predicted by the conformal field theory, which serves asa benchmark for future investigations. From this study, we findthat the designer model discussed above indeed dramaticallysuppresses the finite-size effect, and produces much more ac-curate scaling function of the conductivity, requiring smallersystem sizes in comparison with that in the conventional mod-els. In addition, this free-fermion simulation also providesan estimation about the low-bound of system sizes needed forstudying the scaling of conductivity at the GNY QCP.At the tree level, i.e., free Dirac fermions, the b term inEq. (10) vanishes and the scaling law can be formulated as (cid:104) J x ( ω n ) J x (− ω n )(cid:105) = − σ ∞ ω n − b T ω n + · · · . (11)at ω n (cid:29) T Further more, for free Dirac fermions in ( + ) D, it is knownthat the σ ∞ = /
16 per each fermion species [38]. Be-cause our model contains four Dirac cones (two valleys andtwo spins), we will expect σ ∞ = / (cid:104) J x ( ω n ) J x (− ω n )(cid:105) at thebare Dirac fermion level are shown in Fig. 5, where the resultsin the main panel are for L =
50 system and temperature T = / t = t = /
27) and theconventional model with ( t = t = y = − . x + .
097 is also shown in figure and thecorresponding σ ∞ = .
246 (the constant 1 .
097 which doesnot appeared in Eq. (11) comes from the contribution of back-ground of (cid:104) J x ( ω n ) J x (− ω n )(cid:105) at T → , ω n → t , and this is dueto the fact that our extended linear dispersion in the BZ suc-cessfully reduced the finite energy cut-off due to the deviationof dispersion away from linear.We carried out the same computation for a few differentsystem sizes L (while keeping T = / σ ∞ is shown in the inset of Fig. 5). The dashed line in theinset indicates the exact value σ ∞ = / L → ∞ and it isclear that as L is increasing, the σ ∞ obtained from the fittingare approaching 1 /
4. More importantly, the (red) data pointsfrom the designer model approaches the exact value muchfaster than that of the conventional model (blue), and the onlydifference between these two model is that a small amountof next-nearest-neighbor hopping t = /
27 is introduced inthe designer model. Without t , σ ∞ deviates from the exactvalue by about 10%, even for the largest system size L = L =
30 produces more accurate σ ∞ than the conventional model at L =
50, which implies that thedesigner model allows us to access quantum criticality with thesame or even better accuracy at only about 1 / t .The calculation of the conductivity at the N f = IV. CONCLUSION
Over the years, controlled study of fermionic QCPs is thequestion that have been haunting the minds of physicists in thefield of strongly correlated systems. Although extensive effortshave been devoted and great progress has been achieved, evengreater challenges are still lying in front of the community.Among those, controlled analytical calculation of the criticalproperties and numerical approaches such that the computa-tional complexity ( O ( β N ) for fermionic QMC for example)can be overcome and the thermodynamic limit can be safelyreached, are the important steps towards the final solution ofthis challenging problem. For this objective, our efforts herefocuses on improving the situation at the numeric front, utiliz-ing new techniques in model design guided by insights gainedfrom the field theory studies.The designer Hamiltonian that we have engineered realizesa lattice model of the ( + ) d N f = L =
16, avoiding the heavy computational burden forpushing to larger sizes. In fact, we are able to obtain thecomplete set of critical exponents { / ν GN , η φ , η ψ , ω } thanprevious works (including that of ourselves [21]) and thereforeprovide the much needed benchmark results for the furtherdevelopments at the analytical front such as controlled RG cal-cuation and conformal bootstrap analysis. We also calculatedthe conductivity of the Dirac fermions and find its finite fre-quency scaling behavior consistent with conformal field theoryprediction. Our approach provides the promising direction to-wards the eventual controlled study of fermionic QCPs, andthe same concept and principles can be easily generalized toother fermionic QCPs as well, such as the GNY chiral XY [43]and Heisenberg [19, 23] and the fermion surface coupled to thecritical bosons [17], topological orders [44, 45] and Yukawa-SYK system [46]. Efforts along this line of thinking may helpshed new light on key open questions in various fermionicQCPs and may eventually lead to their final solutions. ACKNOWLEDGEMENT
We thank Michael Scherer, Lukas Jannsen, Joseph Ma-ciejko, Andreas Läuchli, Thomas Lang, Stefan Wessel, FakherAssaad and Xi Dai for insightful discussions on both thefermionic GNY transitions and QMC methodologies, WilliamWitczak-Krempa for enligthening conversations on the finitefrequency conductivity at QCPs, Ning Su and Junchen Rongfor explanations on the conformal bootstrap estimation onthe fermion QCPs and Da-Chuan Lu and Yang Qi and Yi-Zhuang You for the discussion on the calculation of con-ductivity of Dirac fermion at finite frequency with the lat-tice cut-off considered. YZL, WW and ZYM acknowledgethe supports from the Ministry of Science and Technology ofChina through the National Key Research and DevelopmentProgram (Grant No. 2016YFA0300502), the Strategic Prior-ity Research Program of the Chinese Academy of Sciences(Grant No. XDB28000000), the National Science Foundationof China (Grant Nos. 11421092,11574359,11674370) and Re-search Grants Council of Hong Kong Special AdministrativeRegion of China through 17303019. K.S. acknowledges sup-port from the National Science Foundation under Grant No.EFRI-1741618. We thank the Center for Quantum SimulationSciences in the Institute of Physics, Chinese Academy of Sci-ences, the Computational Initiative at the Faculty of Science atthe University of Hong Kong and the National SupercomputerCenter in Tianjin and the National Supercomputer Center inGuangzhou for their technical support and generous allocationof CPU time.
Appendix A: Velocity estimations
In this section, we provide analytic calculations of the ve-locities of free fermion and boson of our model in Eq. (1).Through Legendre transformation of L Fermion , the fermionc H Fermion can be written as H Fermion = (cid:213) ( ij ) ,σ − t ij e − i σθ c † i ,σ c j ,σ + h . c . (A1)In the momentum space, through Fourier transformation, weobtain H Fermion = − t (cid:18) A + BA ∗ + B ∗ (cid:19) − t (cid:18) C + DC ∗ + D ∗ (cid:19) where A = cos [ k x + k y ] exp [ i π ] B = cos [ k x − k y ] exp [− i π ] C = cos [ x + y ] exp [ i π ] D = cos [ x − y ] exp [− i π ] (A2)From the Eq. (A2), one sees the Dirac points are located inthe ( π ,0) and (0, π ) in BZ which we then Taylor expand in theneighborhood of them, H Fermion = √ t (cid:18) − k x − i k y − k x + i k y (cid:19) (A3)So the free fermion velocity v F = √ t is deduced.For the free boson which does not include the term of φ ,considering the following form of bosonic action, S = ∫ dt (cid:213) i ( ∂ t φ i ) − (cid:213) i , j J ij ( φ i − φ j ) − (cid:213) i m φ (A4)where J , J , J and J are given in Fig. 1 (c).Again going to the momentum space via Fourier transfor-mation, we obtain S = (cid:213) ω (cid:213) q [ ω − J [( − cos [ q x ] − cos [ q y ])− ( − cos [ q x ] − cos [ q y ]) + ( − cos [ q x ] − cos [ q y ])− ( − cos [ q x ] − cos [ q y ])] − m ] φ , w (A5)The Eular-Lagrange’s equation require S = φ ω = J [( − cos [ q x ] − cos [ q y ]) − ( − cos [ q x ] − cos [ q y ] + ( − cos [ q x ] − cos [ q y ]))− ( − cos [ q x ] − cos [ q y ])] + m (A6)close to Q = ( , ) , the dispersion is ω = ± (cid:114) J ( q x + q y ) + m (A7)So the velocity of boson is v b = (cid:112) J / v F = v b require the relation of J = . t . Appendix B: projective QMC and SLMC
As mention in the main text, we employ the projectiveQMC [21] to solve the model in Eq. (1). In the QMC simula-tion, the partition function is written as Z = (cid:104) Φ T | exp (− Θ H )| Φ T (cid:105) = (cid:213) { φ } W Boson (cid:214) σ = ↑ , ↓ det ( P † σ B σ m · · · B σ P σ ) , (B1)where Θ plays the role of inverse temperature and φ τ, i is thebosonic field at imaginary time τ and site i . For the bosonicpart of the configurational weight, the weight W Boson is of theform W Boson = exp [− ∆ τ ( (cid:213) i ,τ ∆ τ ( φ i ,τ + − φ i ,τ ) + (cid:213) i , j (cid:213) τ V ij ,τ )] , (B2)where V ij ,τ = J ij ( φ i ,τ − φ j ,τ ) + m φ i ,τ For the fermionic determinant, | Φ T (cid:105) is the trial wave func-tion which is comprised of the eigenvectors of the free feremionHamiltonian for the occupied states, represented by a P σ ma-trix with size N × N where N is the number of electrons atthe half-filling. B στ matrix at time slice τ is then B στ = exp (− ∆ τ ( H Fermion + H Coupling )) = exp (− ∆ τ H Fermion ) exp (− ∆ τ H Coupling ) . (B3)Let’s now compute the ratio for the fermion determinantfor the QMC update. Introduce the notation B ( τ , τ ) = B στ · · · B στ + , B < = P † σ B ( Θ , τ ) and B > = B ( τ, ) P σ , onehas W Fermion = det [ B < B > ] , (B4)where the Green’s function is1 − G ( τ ) = B > ( B < B > ) − B < (B5)0 J eff k J eff1 J eff2 J eff3 -0.0143045485581 -0.205175906625 -0.0539032527282 -0.111523376999TABLE II. Fitted values of J eff k , J eff1 , J eff2 and J eff3 when m = − .
06 for L = Θ = L . The local update of determinant QMC suffers from the criti-cal slowing-down close to the QCP, to overcome this problem,we make use of the recently developed self-learning MonteCarlo (SLMC) [26, 27] scheme to perform more efficientMonte Carlo moves in the configurational space. The coreidea of SLMC is to obtain the low-energy effective bosonicHamiltonian of the problem at hand, and this can be achievedby training an effective model from the configurations savedin a small lattice obtained from the conventional determinantQMC simulation. In our case, the Lagrangian of effectivemodel is of the following form L = L Boson + L eff (B6)with L Boson = (cid:213) p (cid:20) (cid:0) ∂φ p ∂τ (cid:1) + m φ p + φ p (cid:21) + (cid:213) ( p , q ) J pq ( φ p − φ q ) , L eff = J eff k (cid:213) i λ i φ i τ + (cid:213) ( ij , n ) J eff n φ i τ φ j τ (B7)The spirit of the learning process is to replace the effect of L Fermion and L Coupling with L eff . λ i take value of ± J eff k and J eff n -s (the n th-nearest interaction between the bosons) are the trainedeffective couplings which is shown in Tab II from a L = m = − .
06 and Θ = L . With the trained L eff and the bare boson L Boson , one can then use the Lagrangian inEq. (B6) to guide the Monte Carlo move of the boson fieldsand only after substantially many such updates, evalue thefermion determinant in Eq. (B1) to respect the detailed balancecondition of the original model in Eq. (1). For more detaileddescription of SLMC, the assiduous readers are referred to theRefs. [18, 26, 27].One more obstacle particularly associated with the currentproblem is that for training the effective model, the calculationof the fermion determinant is necessary, but the det [ B < B > ] canbe very large (as large as e ) in our model. To overcome suchnumerical difficulty, we note that the real important quantity weneed to calculate is log ( det [ B < B > ]) in the SLMC scheme [27].If the det [ B < B > ] is divided into product of several parts (suchas T = A A A ), then the large value in the determinant canbe breaked into sum several small parts (eg. log det [ T ] = log det [ A ] + log det [ A ] + log det [ A ] ), where each part issmall enough to be handled numerically, i.e.,det [ P † B s ( Θ , ) P ] = det [ P † B n B n − · · · B B P ] , (B8)the UDV decomposed can be used in the model B > = U > DV (B9) where U > is a rectangular matrix 2 N × N and D contains the N eigenvalues of scale from large to small and V is an upperunit triangular matrix N × N . In this way, we can obtaindet [ P + B n B n − · · · B B P ] = det [ P + B n B n − · · · B B P ] = det [ P + B n B n − · · · B U > D V ] = det [ P + B n B n − · · · U > D V D V ]· · · = det [ P + U > n D n V n D n − V n − · · · D V D V ] = det [ P + U > n ] det [ D n D n − D n − · · · D D ] det [ V n V n − · · · V V ] , (B10)such that we can break the determine into three parts andcompute the weight accordingly.With the determiant weight obtained for small system sizes,we can train and obtain the parameters for the effective model.As the example shown in Tab. II, one can see that interactions ofseveral different distances J eff1 , , are all important as they havethe same scale. With the cumulative update [27] in SLMC,it turns out that the simulation with such effective model butapplied to L ≥ Appendix C: Calculation of bosonic and fermionic dispersion
Here we take a closer look at how to obtain the fermionicand bosonic dispersion shown in the Fig. 1 (d) and (e) of themain text. In the projector QMC, dynamical Green’s function G f ( k , τ ) and G b ( q , τ ) in the the imaginary time correlationfollow the Lehmann spectral representation such that at longtime they satisfy the relation G ( k , τ ) ∝ e − ∆ ( k ) τ . So for eachfinite size system with L momenta along the lattice axis, ateach moment point k , one can obtain the excitation gap, i.e.the dispersion, from the following fitting scheme,ln G ( k , τ ) = − ∆ ( k ) τ + C , (C1)with C a fitting parameter. As an example, here we showthe data in Fig. 6 (a) and (b) for fermion and boson Green’sfunction at the closest (I) and next-closest (II) momenta awayfrom their gapless points for L =
14 system. For the fermionicpart (shown in Fig. 6 (a)), the obtained gap ∆ ( k = ( π, π L )) = . ( ) and ∆ ( k = ( π, π L )) = . ( ) , are the values shown inthe Fig. 1 (d) which are close to the the free Dirac dispersionof our lattice model. For the bosonic part (shown in Fig. 6(b)), the finite size effect is a bit stronger here and at longerimaginary time, the data exhibit faster decay (at such τ the G b ( τ ) -s are very close to zero which can not be used forfitting). Nevertheless, one can still see the clear exponential1 -5 -4 -3 -2 -1
60 0.5
FIG. 6. Fermion and Boson Green’s functions and their fitting ac-cording to Eq. (C1) from which the finite size gap ∆ ( k ) at differentmomenta are obtained, these gaps, prepared on all the momenta alongthe high-symmetry-path, give rise to the dispersions ω ( k ) shown inFig. 1 (d) and (e) of the main text. part of the G b ( q , τ ) and the fitting here gives rise to ∆ ( q = ( , π L )) = . ( ) and ∆ ( q = ( , π L )) = . ( ) , are the valuesshown in the Fig. 1 (e). We believe the deviation of ∆ ( k = ( , π L )) = . ( ) from the linear bosonic critical dispersion isa finite size effect and this will be overcome by similar analysisbut with even larger L , which we leave for future works. [1] J. A. Hertz, Phys. Rev. B , 1165 (1976).[2] A. J. Millis, Phys. Rev. B , 7183 (1993).[3] T. Moriya, Spin Fluctuations in Itinerant Electron Magnetism (Springer-Verlag Berlin Heidelberg, 1985).[4] A. Abanov, A. V. Chubukov, and J. Schmalian, Advances inPhysics , 119 (2003).[5] A. V. Chubukov, C. Pépin, and J. Rech, Phys. Rev. Lett. ,147003 (2004).[6] A. V. Chubukov and D. L. Maslov, Phys. Rev. Lett. , 216401(2009).[7] S.-S. Lee, Phys. Rev. B , 165102 (2009).[8] M. A. Metlitski and S. Sachdev, Phys. Rev. B , 075127 (2010).[9] M. A. Metlitski and S. Sachdev, Phys. Rev. B , 075128 (2010).[10] A. Schlief, P. Lunts, and S.-S. Lee, Phys. Rev. X , 021010(2017).[11] A. W. Sandvik, AIP Conference Proceedings , 135 (2010),https://aip.scitation.org/doi/pdf/10.1063/1.3518900. [12] Y. Schattner, S. Lederer, S. A. Kivelson, and E. Berg, Phys.Rev. X , 031028 (2016).[13] X. Y. Xu, K. Sun, Y. Schattner, E. Berg, and Z. Y. Meng, Phys.Rev. X , 031058 (2017).[14] M. H. Gerlach, Y. Schattner, E. Berg, and S. Trebst, Phys. Rev.B , 035124 (2017).[15] Z. H. Liu, X. Y. Xu, Y. Qi, K. Sun, and Z. Y. Meng, Phys. Rev.B , 045116 (2018).[16] Z. H. Liu, X. Y. Xu, Y. Qi, K. Sun, and Z. Y. Meng, Phys. Rev.B , 085114 (2019).[17] Z. H. Liu, G. Pan, X. Y. Xu, K. Sun, and Z. Y. Meng, Proceedingsof the National Academy of Sciences , 16760 (2019).[18] X. Y. Xu, Z. H. Liu, G. Pan, Y. Qi, K. Sun, and Z. Y. Meng,Journal of Physics: Condensed Matter , 463001 (2019).[19] N. Zerf, L. N. Mihaila, P. Marquard, I. F. Herbut, and M. M.Scherer, Phys. Rev. D , 096010 (2017). [20] B. Ihrig, L. N. Mihaila, and M. M. Scherer, Phys. Rev. B ,125109 (2018).[21] Y.-Y. He, X. Y. Xu, K. Sun, F. F. Assaad, Z. Y. Meng, and Z.-Y.Lu, Phys. Rev. B , 081110 (2018).[22] S. Chandrasekharan and A. Li, Phys. Rev. D , 021701 (2013).[23] T. C. Lang and A. M. Läuchli, Phys. Rev. Lett. , 137602(2019).[24] M. Schuler, S. Hesselmann, S. Whitsitt, T. C. Lang, S. Wessel,and A. M. Läuchli, arXiv e-prints , arXiv:1907.05373 (2019),arXiv:1907.05373 [cond-mat.str-el].[25] B. Roy, V. Juričić, and I. F. Herbut, Journal of High EnergyPhysics , 18 (2016).[26] J. Liu, Y. Qi, Z. Y. Meng, and L. Fu, Phys. Rev. B , 041101(2017).[27] X. Y. Xu, Y. Qi, J. Liu, L. Fu, and Z. Y. Meng, Phys. Rev. B ,041119 (2017).[28] Y. Otsuka, S. Yunoki, and S. Sorella, Phys. Rev. X , 011029(2016).[29] S. Weinberg, The quantum theory of fields (Cambridge Univer-sity Press, 2005).[30] F. Assaad and H. Evertz, “World-line and determinantal quan-tum monte carlo methods for spins, phonons and electrons,”in
Computational Many-Particle Physics , edited by H. Fehske,R. Schneider, and A. Weiße (Springer Berlin Heidelberg, Berlin,Heidelberg, 2008) pp. 277–356.[31] Z. Y. Meng, T. C. Lang, S. Wessel, F. F. Assaad, and A. Mura-matsu, Nature (2010), 10.1038/nature08942. [32] L. Iliesiu, F. Kos, D. Poland, S. S. Pufu, and D. Simmons-Duffin,Journal of High Energy Physics , 36 (2018).[33] J. Rong and N. Su, arXiv e-prints , arXiv:1807.04434 (2018),arXiv:1807.04434 [hep-th].[34] K. Binder, Zeitschrift für Physik B Condensed Matter , 119(1981).[35] S. Pujari, T. C. Lang, G. Murthy, and R. K. Kaul, Phys. Rev.Lett. , 086404 (2016).[36] H. Shao, W. Guo, and A. W. Sandvik, Science , 213 (2016).[37] Y. Q. Qin, Y.-Y. He, Y.-Z. You, Z.-Y. Lu, A. Sen, A. W. Sandvik,C. Xu, and Z. Y. Meng, Phys. Rev. X , 031052 (2017).[38] E. Katz, S. Sachdev, E. S. Sørensen, and W. Witczak-Krempa,Phys. Rev. B , 245109 (2014).[39] W. Witczak-Krempa, E. S. Sørensen, and S. Sachdev, NaturePhysics , 361 (2014).[40] K. Chen, L. Liu, Y. Deng, L. Pollet, and N. Prokof’ev, Phys.Rev. Lett. , 030402 (2014).[41] A. Lucas, S. Gazit, D. Podolsky, and W. Witczak-Krempa, Phys.Rev. Lett. , 056601 (2017).[42] S. M. Chester, W. Landry, J. Liu, D. Poland , D. Simmons-Duffin, N. Su, and A. r. Vichi, arXiv e-prints , arXiv:1912.03324(2019), arXiv:1912.03324 [hep-th].[43] Y. Da Liao, Z. Y. Meng, and X. Y. Xu, Phys. Rev. Lett.123