Universal quantum criticality in the metal-insulator transition of two-dimensional interacting Dirac electrons

The metal-insulator transition has been a subject of intense research since Nevil Mott has first proposed that the metallic behavior of interacting electrons could turn to the insulating one as electron correlations increase. Here, we consider electrons with massless Dirac-like dispersion in two spatial dimensions, described by the Hubbard models on two geometrically different lattices, and perform numerically exact calculations on unprecedentedly large systems that, combined with a careful finite size scaling analysis, allow us to explore the quantum critical behavior in the vicinity of the interaction-driven metal-insulator transition. We find thereby that the transition is continuous and determine the quantum criticality for the corresponding universality class, which is described in the continuous limit by the Gross-Neveu model, a model extensively studied in quantum field theory. We furthermore discuss a fluctuation-driven scenario for the metal-insulator transition in the interacting Dirac electrons: the metal-insulator transition is triggered only by the vanishing of the quasiparticle weight but not the Dirac Fermi velocity, which instead remains finite near the transition. This important feature cannot be captured by a simple mean-field or Gutzwiller-type approximate picture, but is rather consistent with the low energy behavior of the Gross-Neveu model.


I. INTRODUCTION
The metal-insulator transition is one of the most fundamental and yet profound physical phenomena of quantum mechanics, and, in the absence of correlations, is described by the conventional band theory [1,2]. A metal is such if electrons do not fill an integer number of bands, otherwise insulating behavior settles because an energy gap is required to excite an electron from a fully occupied to an empty band. With this simple criterion, most insulating and metallic properties were successfully explained [3]. However, it was soon realized by Mott [4] in 1949 that the electron correlations could play a major role in several materials, as they could become insulators even when the band theory predicts instead metals: these are the so called Mott insulators. Since then, many theoretical and numerical works have tried to shed lights on this issue, but our understanding of interaction-driven metal-insulator transitions still remains rather controversial because strongly correlated systems are hard to solve using both analytical and numerical methods, at least, when the spatial dimensionality is larger than one [5] but smaller than infinity [6].
Gutzwiller has introduced in the middle 60's a correlated framework [7], that was later used to derive the properties of the metal-insulator transition as a function of the correlation strength U . This framework predicts * otsukay@riken.jp for generic lattice models that for U below the critical point U c , the quasiparticle weight Z, which should be exactly one in the non-interacting band theory, is strongly renormalized by the correlation and vanishes as Z ≃ (U c − U ). At the same time, the bandwidth W , renormalized by the electron correlations, reduces to zero at the transition in the same way as Z vanishes. The prediction of the Brinkman-Rice approximation [8] has been later confirmed and further extended by the dynamical mean-field theory (DMFT) [9,10], an approach that is exact only in the limit of large spatial dimensions.
Here, we focus on a specific realization of the metalinsulator transition in two-dimensional lattice models which can be treated with a numerically exact method, i.e., the Hubbard models defined on the honeycomb lattice and on the square lattice with π flux penetrating each plaquette. These models are equipped with a free electron energy dispersion with nodal gapless points in the Brillouin zone and with linear dispersion (see Figs. 1), a very peculiar character of the so called "massless Dirac electrons". We set one electron per site, where the noninteracting band is half-filled and the Fermi surface is constituted by the Dirac points. Due to these gapless Dirac points, it is possible to have a non trivial metalinsulator transition at a finite value of the correlation strength even in such bipartite lattices [11].
Quite recently, a numerical simulation of the Hubbard model on the honeycomb lattice has provided evidences for a possible unconventional phase, a spin liquid phase with no classical order, close to the metal-insulator tran-sition occurring at sufficiently large U [12]. Although there is still activity [13], the possibility of such an intermediate phase between the semi-metal (SM) and the antiferromagnetic (AF) Mott insulator seems now rather unlikely, in view of the large scale simulations that we have reported recently [14], clearly showing that the AF moment develops continuously from zero once the insulating phase is entered. Later studies have also confirmed the simplest scenario of a direct and continuous transition [15][16][17][18].
Similarly, a stable spin liquid in interacting Dirac electrons represented by a different model has been also proposed in Ref. [19]. Here, the Hubbard model on the square lattice is studied, in which a flux π is added to each plaquette in order to obtain a massless Dirac dispersion in the non-interacting limit [Figs. 1(b) and 1(d)] (referred to as the π-flux model hereafter). Based on an approximate numerical simulation with relatively small clusters, a spin liquid phase has been observed between the SM and the AF Mott insulator, which in this case only has a finite charge gap with a vanishingly small spin gap [19]. This finding is significant also in the context of high-T c cuprate superconductors because the π-flux model is considered as one of the relevant models to understand the mechanism of superconductivity [20,21]. However, as in the case of the honeycomb lattice model, this quantum disordered state has also been disputed [17,18,22].
Although it is obviously very important to search for spin liquid phases in "realistic" models, here we take a different perspective. After several years of efforts on these strongly correlated systems, we feel that the time is mature to examine the quantum criticality in the metalinsulator transition of interacting electrons in two spatial dimensions, and in particular the interacting Dirac electrons described by these two models where U c is finite and their ground state properties can be explored by using an unbiased and formally exact numerical method. This is precisely the main purpose of this paper.
Moreover, it has been recently shown that the interacting Dirac electrons on the honeycomb lattice can be mapped in the continuous limit onto a model well known in quantum field theory, i.e., the Gross-Neveu (GN) model [23] in the chiral Heisenberg universality class with N = 8 fermion components [24,25]. Since the π-flux model is expected to have a similar effective theory in the continuous limit, it is reasonable to conjecture that the metal-insulator transition in these two specific lattice models belongs to the same universality class. In order to address this issue, here we perform large-scale quantum Monte Carlo (QMC) calculations and evaluate the critical exponents with a high degree of accuracy. This is indeed made possible because, with the help of the auxiliary field technique [26][27][28], these fermionic models can be studied without the notorious "sign problem" [29,30]. The careful finite size scaling analysis finds that the critical exponents for these two models are the same within statistical errors and thus confirms the conjecture. Our results represents the first accurate determination of the  1), where the lattice constant between the nearest neighbor sites is set to be one. The unit cells for both models thus contain two sites. The nearest neighbor hopping parameters are indicated by t (solid lines) and t ′ (dotted lines). t ′ is set to be −t for the π-flux model in (b). The cluster with L = 3 (L = 4) for the honeycomb (π-flux) model is indicated by blue dashed line. The unit vector ex ( ey) along the x (y) direction is also indicated in (b). The non-interacting energy dispersions ε k are shown for (c) the honeycomb lattice model and (d) the π-flux model. The Fermi level is at ε k = 0 for half-filling, and the Dirac points are located exactly at the Fermi level, where the valence and the conduction bands touch with opposite chiralities. Notice that in both models there are two distinct Dirac points at (c) k = 2π 3 (± 1 √ 3 , 1) and (d) k = (± π 2 , 0), corresponding to two distinct valleys and thus there are eight components of Dirac fermions in total due to different chiral, spin, and valley degrees of freedom.
critical exponents for the GN model in the chiral Heisenberg universality class with N = 8 [18,31,32].
The other interesting issue to be addressed in this paper is to explore the quantum critical behavior in both metallic and insulating phases at the vicinity of the metal-insulator transition, in particular, the fate of the quasiparticle weight Z and the Fermi velocity v F when approaching the critical point U c from the metallic side. For electrons with the usual energy dispersions such as the one in the square lattice, the Gutzwiller-type approximate description [7,8] and the simple DMFT approach [9,10] predict that Z and v F are both renormalized by the interaction, and vanish at U c . This scenario is valid for any lattice model and in any dimensionality within the Gutzwiller approximation since, within this method, the free electron dispersion is simply renormalized by a Gutzwiller factor Z that vanishes at the transition. Analogously, the same scenario holds within the single-site DMFT [33] because, once the self-energy is assumed to be momentum-independent, the free electron dispersion can be renormalized only through the quasiparticle weight Z [9]. Instead, our unbiased and numerically exact calculations support the qualitative prediction based on the renormalization group (RG) analysis for the GN model [34] and the recent numerical results for the honeycomb lattice model obtained by advanced quantum cluster methods [35,36]: with increasing the correlation strength, Z vanishes at the transition, while the Fermi velocity v F remains finite.
Our large-scale QMC calculations also provide a firm numerical evidence for the absence of a spin liquid phase in between the SM and the AF insulator for the π-flux model, thus ruling out the possibility of a spin liquid phase reported previously in Ref. [19]. This is very similar to the case for the honeycomb lattice model, where the originally proposed spin liquid phase [12] is turned out to be rather implausible after our large-scale calculations [14]. The metal-insulator transitions in both models are rather direct and continuous, and can be characterized by the quantum critical behavior of the quasiparticle weight in the metallic phase and the antiferromagnetic order parameter in the insulating phase. These results therefore suggest that the electron correlation alone is not enough but other factors such as geometrical frustration are required for a magnetically disordered spin liquid state [37].
The rest of this paper is organized as follows. The definition of the two models and a brief description of the QMC method employed are given in Sec. II. The ground state phase diagrams are first obtained in Sec. III by a rather conventional way of extrapolating order parameters to the thermodynamics limit. Section IV is devoted to more detailed analysis to determine the critical exponents with high accuracy. The fate of the Fermi velocity is investigated in Sec. V. Finally, the results are discussed in the context of the GN model, followed by an outlook and conclusions, in Sec. VI. The energy resolved momentum distribution function is described in Appendix A and the leading correction to the scaling analysis is discussed in Appendix B.

II. MODELS AND METHOD
We consider two variants of the Hubbard models in two spatial dimensions, whose low-lying energy states are described by the interacting Dirac fermions with spin 1/2 degree of freedom at half-filling. The Hamiltonian in standard notations readŝ where c † is is the creation operator of electron at site i and spin s (=↑, ↓), n is = c † is c is , and the sum i, j runs over all pairs of nearest neighbor sites i and j. The first model is defined on the honeycomb lattice with the uniform hopping t ij = t [see Fig. 1(a)]. The second one is on the square lattice with a flux of π penetrating on each square plaquette, represented, with an appropriate gauge transformation, by t i,i+ ex = t and t i,i+ ey = (−1) ix+iy t, where the position of site i is given as i x e x + i y e y and e x ( e y ) denotes the unit vector along the x (y) direction [see Fig. 1 The clusters considered here consist of (L τ 1 , L τ 2 ) with N s = 2L 2 sites for the honeycomb lattice model and (L e x , L e y ) with N s = L 2 sites for the π-flux model, as indicated in Figs. 1(a) and 1(b), respectively, with periodic boundary conditions. The number of electrons are set to be equal to the number of sites in both models. In order to include the Dirac points among the allowed momenta in the non-interacting energy dispersions, L is chosen to be a multiple of three (four) for the honeycomb lattice (π-flux) model. The smallest clusters are indicated by dashed lines in Figs. 1(a) and 1(b). The largest clusters considered here are N s = 2, 592 sites for the honeycomb lattice model and N s = 1, 600 for the π-flux model.
Although the two models are quite different, they are both characterized by the non-interacting energy dispersions ε k with two gapless Dirac cones, as shown in Figs. 1(c) and 1(d), leading to a semi-metallic behavior at half-filling and for small coupling U/t [29,30]. The effective low-energy HamiltonianĤ 0 eff in the non-interacting limit at the vicinity of the Dirac points for spin s is described asĤ where δk = (δk x , δk y ) is the momentum measured from the Dirac point, v 0 F = 3t/2 (2t) is the Dirac Fermi velocity in the non-interacting limit for the honeycomb lattice (πflux) model, and σ = (σ x , σ y , σ z ) are the Pauli matrices acting on the two different sublattices. As shown below, both models display the metal-insulator transitions at finite critical values U c /t from the non-magnetic SM to the AF long-range ordered insulating phase, in good agreement with previous numerical studies [29,30]. It should also be noted that in the simplest mean-field picture the insulating phase emerges because the mass term proportional to σ z is introduced Eq. (2) when an AF order sets in. Therefore, there is no unit cell doubling in the AF insulating phase for both models [38].
We employ the auxiliary field QMC method [26][27][28] to investigate the ground state properties of these two models. The expectation value of a physical observableÔ over the ground state |Ψ 0 ofĤ is obtained by projecting out trial wave functions to the ground state, i.e., where and |ψ L (|ψ R ) is the left (right) trial wave function, chosen to have finite overlap with the exact ground state. We choose |ψ L as a mean-field wave function of Eq. (1) with an AF order parameter in the x direction, while a Slater determinant of the non-interacting Hamiltonian is used for |ψ R to which a tiny perturbation term is added to remove the degeneracy at the two Dirac points. These choices for the trial wave functions have been shown to yield a particularly fast convergence in the imaginary time projection onto the ground state [14].
The imaginary time evolution operator e −τĤ with the projection time τ is divided into N τ pieces, i.e.,

Nτ
, where τ = ∆τ N τ and N τ is the Trotter number (integer). By setting ∆τ t ≪ 1, we can use the Suzuki-Trotter decomposition [39,40], e −∆τĤ = e − 1 2 ∆τĤ0 e −∆τĤI e − 1 2 ∆τĤ0 + O ∆τ 3 , whereĤ 0 is the hopping term andĤ I is the interacting term of the Hubbard modelĤ in Eq. (1). Notice that the systematic error introduced in this decomposition is O ∆τ 3 . The discrete Hubbard-Stratonovich transformation is applied to −∆τĤI , which introduces an auxiliary Ising field at each site as well as at each imaginary time slice [41]. As shown in Fig. 2, we have confirmed that the systematic errors due to finite τ and ∆τ are sufficiently small, compared to the statistical errors in Monte Carlo importance sampling, when we choose τ = L + 4 and ∆τ t = 0.1. More technical details are found in our previous report [14].

III. GROUND STATE PHASE DIAGRAM
In this section, we shall focus on the continuous nature of the quantum phase transition between the SM and the AF insulator. For this purpose, we calculate two fundamental quantities, the staggered magnetization and the quasiparticle weight, which characterize two different aspects across the transition. The former quantity reveals the magnetic transition to the AF state and the latter one directly captures the metal-insulator transition. These transitions are expected to occur at the same critical U c , unless there is an intermediate phase such as the spin liquid phases [12][13][14][15][16]19]. Here in this section we take a conventional and straightforward way, i.e., by first calculating the staggered magnetization and the quasiparticle weight on different finite clusters and then extrapolating them in the thermodynamic limit, to obtain the ground state phase diagram as a function of U/t.

A. staggered magnetization
The staggered magnetization on each finite cluster with a linear dimension L, expressed as m s (L), is calculated from the spin structure factor S AF (L), i.e., where At a fixed value of U/t, the staggered magnetizations m s (L) are calculated on clusters with L = 6, 9, 12, 15, 18, 24, and 36 for the honeycomb lattice model, and with L = 8, 12, 15, 18, 24, 32, and 40 for the π-flux model, and are extrapolated to the thermodynamic limit using polynomial functions in 1/L. The typical results are shown in Fig. 3. It is observed that such a simple functional form represents the data rather well and thus we can estimate the AF order parameter m s reasonably accurately. The results for the extrapolated AF order parameter m s are summarized in Fig. 4. The extrapolated values to the thermodynamic limit are also indicated at 1/L = 0.
The critical points U c , above which the AF order parameter m s is finite, and the critical exponents β are estimated by assuming a form of the AF order parameter as a function of U as The estimated U c /t is 3.85 ± 0.02 for the honeycomb lattice model and 5.65 ± 0.05 for the π-flux model, as indicated in Fig. 4. Notice that U c for the π-flux model is larger than that for the honeycomb lattice model. This is easily understood because v 0 F for the former model is larger than that for the latter model [see Eq.
(2) and also Figs. 1(c) and 1(d)] and thus a larger U is required to induce the AF order. Although U c /t is different for these two models, our calculations find that the critical exponents β for the two models are the same within statistical errors, i.e., β = 0.75 ± 0.06 for the honeycomb lattice model and 0.80 ± 0.09 for the π-flux model, as indicated in Fig. 4, which will be confirmed in more details by a careful and accurate finite-size scaling analysis in Sec. IV.

B. quasiparticle weight
Next, we study the metal-insulator transition by considering the momentum distribution function, i.e., the ground state occupation of the one electron states labeled by momentum k, spin s, and non-interacting en- ergy ε k . We calculate this occupation number n(ε k ) [defined in Eq. (A3)] as a function of ε k , where we average over equivalent momenta with the same energy (see Appendix A for the details). Typical results of the "energy resolved" momentum distribution function n(ε k ) are shown in Fig. 5 for both models. For U < U c , a jump in the momentum distribution function n(ε k ) occurs for ε k → 0 when k approaches the Dirac points. The singularity in n(ε k ) implies a long-distance power-law behavior in the density matrix c † is c js in real space, which is the fingerprint of a metal. From general grounds [42], by applying the well known "Migdal theorem" [43], the quasiparticle weight Z can be related to the jump in the momentum distribution function. Therefore, we can have direct access to Z in the thermodynamic limit. However, the finite size effects are rather significant and need to be carefully controlled in order to reach definite conclusions based on the available finite size calculations.
It should also be noticed in Fig. 5 that the energy resolved momentum distribution function n(ε k ) becomes smooth without a visible discontinuity at the Fermi level for large U/t. This is interpreted as an exponential decay of the density matrix at large distance, because the density matrix is just the Fourier transform of the momentum distribution function. This clearly indicates the presence of a gap in the charge sector.
We find that the following procedure works for estimating Z in the thermodynamic limit. Since the "energy resolved" momentum distribution function n(ε k ) is smooth near ε k = 0, we analyze n(ε k ) calculated on different finite clusters by extrapolating it to ε k = 0 with a polynomial function. More specifically, the quasiparticle weight on each finite cluster, Z L , is defined as the jump in n(ε k ), which is evaluated by extrapolating to the Fermi level (ε k = 0) the closest three data points of n(ε k ) for Solid curves are least-square fits of the three data points closest to the Fermi level with ε k < 0 and ε k > 0 using quadratic polynomials of ε k . Notice that, due to the particle-hole symmetry, n(ε k > 0) = 1 − n(ε k < 0) and thus n(ε k = 0) = 1/2. ZL and ∆n(u, L) are indicated in (a) for U/t = 4 and in (b) for U/t = 5.8. ε k < 0 or ε k > 0 with a quadratic polynomial, as shown in Fig. 5. Notice that due to the particle-hole symmetry the two extrapolations using the data points for ε k < 0 and for ε k > 0 are related and in practice only one of the two is necessary. The extrapolated quantity Z L for finite L is certainly much closer to the value in the thermodynamic limit and quite generally converges smoothly to the quasiparticle weight with a quadratic polynomial function of 1/L, as shown in Fig. 6. However, this extrapolation is valid only in the metallic regime, where the assumed polynomial convergence in 1/L is justified. Indeed, we find in Fig. 6 that, for larger U/t, Z L is extrapolated to a negative value, instead of being positive. This is clearly an inconsistent extrapolation because Z L in the insulating phase is expected to converge exponentially. Nevertheless, these inconsistent extrapolations are very useful as they allow us to identify the extension of the metallic region and determine the critical U MI c for the metal-insulator transition. It is indeed estimated in Fig. 6 that U MI c /t ∼ 3.9 for the honeycomb lattice model and U MI c /t ∼ 5.6 for the π-flux model, which are in excellent agreement with the critical U c for the AF order (see Fig. 4).
The results for the quasiparticle weight Z in the thermodynamic limit as a function of U/t are summarized in Fig. 4 for both models. The obtained Z for U < U MI c appears well behaved and can provide the critical U MI c by assuming that  for the π-flux model [44]. The critical values for U MI c are thus consistent with those for U c determined from the AF order parameter m s within the statistical errors. These results clearly imply that the transition from the SM to the AF insulator is continuous in both models, where the insulating behavior shows up immediately when the AF order is developed for U ≥ U c , i.e., U MI c = U c . Therefore, our large-scale calculations exclude the intermediate phases previously reported in Refs. [12] and [19], and reveal a continuous transition between the SM and the AF insulator. As shown in Sec. IV, the careful finite-size scaling analysis finds that the data collapse fits for both m s and Z are convincing with setting U MI c = U c , also suggesting that the quantum critical points for these two quantities are located at the same U value.
Finally, we remark a semantic issue. The "Mott transition" is very widely used to describe a metal-insulator transition driven by the electron correlation, regardless of whether the symmetry may or may not be broken in the insulating phase. There are certainly more confusion and ambiguity on how to define properly a "Mott insulator", especially in our cases when the unit cell contains an even number of electrons and the AF order found in the insulating phase does not break the translation symmetry. However, this semantic issue is irrelevant to the main purpose of our study and does not change our conclusions.

IV. FINITE-SIZE SCALING ANALYSIS
Having established the continuous character of the transition, let us now evaluate the critical exponents which characterize the quantum phase transition. For this purpose, here we employ the careful finite-size scaling analysis for staggered magnetization and quasiparticle weight.

A. staggered magnetization
For the staggered magnetization, we make use of the standard finite-size scaling ansatz [45][46][47], where ν is the critical exponent of the correlation length ξ ≃ |U − U c | −ν , f m (uL 1/ν ) denotes a model dependent scaling function and u = (U − U c )/U c is the reduced coupling. Notice that the u dependence of m s (L) given in Eq. (5) is explicitly indicated here in Eq. (11) as m s (u, L).
In the above finite-size scaling ansatz, we also take into account the leading correction term, a term proportional to cL −ω , with c and ω being additional fitting parameters [48], which is however expected less important for sufficiently large cluster sizes. The finite-size scaling analysis is performed with a recently proposed method based on the Bayesian statistics [49]. The remarkable advantages of this method are i) the weak dependence on the initial fitting parameters and ii) the applicability to a wide range of the reduced coupling u.
In order to estimate reliable error bars for the fitting parameters in Eq. (11), we adopt a straightforward resampling technique. The fitting procedure is summarized as follows. First, we prepare a data set to be fitted, based on the raw QMC data of m s (u, L) for various u and L, by adding to m s (u, L) a Gaussian-distributed noise with the zero average and the standard deviation estimated by the QMC calculations for m s (u, L). Second, we pick up at random initial values of the fitting parameters U c /t, ν, β, and ω around the optimal ones, namely, U c /t = 3.8 (5.5) for the honeycomb lattice (π-flux) model, ν = 1.0, β = 0.8, and ω = 0.8. Third, with these initial fitting parameters, we perform the data collapse analysis for each resampled data for m s (u, L) based on the Bayesian statistics [49] and obtain the best converged values of U c /t, ν, β, and ω. We repeat this procedure typically a few thousands times to average the converged parameters and estimate the statistical errors. Typical examples of the resampling procedure are shown in Figs. 7 and 8. The clear advantage of the resampling technique is that the degree of uncertainty of the fitting parameters can be immediately verified and therefore a reliable estimate of the error bars can be safely obtained.
The results of the data collapse for the staggered magnetization m s (u, L) are shown in Fig. 9, confirming that our numerical calculations are quite accurate for this quantity since the data for different L collapses almost perfectly into a universal curve. We find that the critical exponents are quite stable and converged to ν ≃ 1 and β ≃ 0.75 for both models, as indicated in Fig. 9. It should be noticed that, for both models, the values of U c , obtained from the data collapse plots in Fig. 9, agree within two standard deviations with the ones estimated straightforwardly by extrapolating m s (L) to the thermodynamic limit for each U (see Figs. 3 and 4).
It should be noted here that very recently the similar QMC method has been applied to the honeycomb lattice Ref. [18]. In their report, the critical exponents for the honeycomb lattice model are estimated as ν ≃ 0.84 and β ≃ 0.71, and are claimed to be consistent with those obtained by the ǫ-expansion for the GN model [25] (see also  Table I), while the critical exponents for the π-flux model are unavailable. We argue that the disagreement of the critical exponents between their estimations and ours for the honeycomb model and the difficulty to determine the critical exponents for the π-flux model in Ref. [18] are due to the limited lattice sizes used in Ref. [18], i.e., up to 648 sites for the honeycomb lattice model and 784 sites for the π-flux model. As shown in Figs. 10(a) and 10(b), our numerical data do not provide perfect data collapse plots for both models when we fix U c , ν, and β which are reported in Ref. [18]. However, the scattering of the data is particularly evident only for the two largest sizes, not studied in Ref. [18]. Indeed, if these largest data sets are excluded from the plots, the data collapse seems acceptable even by using the values of U c , ν, and β reported in Ref. [18], especially for the honeycomb lattice model shown in Fig. 10(a). On the other  Ref. [18]. In Ref. [18], the critical exponents ν and β are estimated, based on their numerical data and renormalization group analysis, only for the honeycomb lattice model, and are assumed to be the same for the π-flux model. For comparison, data collapse fits without the leading correction term are also shown in the lower panels for (c) the honeycomb lattice model and (d) the π-flux model, where Uc, ν, and β are determined by the resampling technique using the data sets for L ≥ 9 in (c) and L ≥ 12 in (d) (see also Table. II in Appendix B). In all figures, the same data are used as those in Fig. 9, but the leading correction term is not considered, i.e., c = 0 in Eq. (11), for a fair comparison.

B. quasiparticle weight
As far as the critical behavior in the charge sector is concerned, the scaling ansatz is applied to the jump ∆n(u, L) of the momentum distribution function n(ε k ) across the Fermi level with a form where η ψ is the anomalous dimension of the fermion field Ψ [50] and f n (uL 1/ν ) is a scaling function. The jump ∆n(u, L) is simply obtained as a difference of n(ε k ) at the two closest points to the Fermi level, i.e., the ones above and below the Fermi level (see Fig.5). Here we do not take Z L as a scaling quantity in Eq. (12) to avoid possible artifacts caused by the extrapolation procedure for Z L where n(ε k ) is extrapolated to ε k = 0. The finite-size scaling for ∆n(u, L) is certainly more difficult because ∆n(u, L) is the direct finite size jump of n(ε k ) across the Fermi level and is much larger than the quasiparticle weight Z in the thermodynamic limit, as shown in Fig.5. Therefore, we fix the critical U c and the expo-nent ν in Eq. (12) to the values already determined by the finite-size scaling analysis on m s (u, L). Moreover, we do not consider the correction term in the finite size scaling analysis (see Appendix B for the case with the leading correction term). The remaining parameter η ψ is determined using the resampling technique described in Sec. IV A. Despite the above simplifications, we find in Fig. 11 that the collapse plots are excellent also in this case, suggesting that the lattice sizes considered here are large enough and that we can faithfully describe the critical behavior also in the charge sector. The obtained values for the critical exponent η ψ are found to be the same for both models (η ψ ≃ 0.21 − 0.22) within statistical errors. However, these values differ from those obtained in Eq. (10) by taking the simple power law fit of Z L to extrapolate in the thermodynamic limit (see Figs. 4 and 6) [44]. We believe that the finite-size scaling analysis based on Eq. (12) results in a more accurate estimation of the exponent, since all data, not only in the metallic region for U < U c but also in the insulating region for U > U c , are used. The obtained critical exponents are summarized in Table I for the honeycomb lattice model and the π-flux model. As clearly shown in Table I, these two different models lead to the same critical exponents for the spin and charge sectors with a considerable degree of accuracy. Therefore, our numerical results firmly verify the universal quantum criticality in the apparently different lattice models, which share only the massless Dirac energy dispersion in the non-interacting limit. To our knowledge, this represents one of the first unbiased studies on the critical properties of the metal-insulator transition in two spatial dimensions. The highly accurate estimation of the critical exponents as well as the critical quantities follows from the unprecedentedly large-scale simulations that we are now able to perform [14], combined with the careful and accurate finite size scaling.

C. scaling functions
In order to firmly establish the universal character of the metal-insulator transition, we follow the Privman and Fisher's argument [51]. It states that if two different models belong to the same universality class, the scaling functions of a physical quantity for these two models are related by non-universal metric factors as where f α (x) and f ′ α (x) are the scaling functions for the two models. Figures 12(a) and 12(b) show the collapse fits without the correction term for the staggered magnetization and for the jump in n(ε k ), respectively. We clearly find in Fig. 12 that the collapse fits for the honeycomb lattice model can superpose onto the appropriately scaled collapse fits for the π-flux model with the nonuniversal constants c 1 and c 2 . This implies that the scaling functions of the two models are essentially equal and confirms in a robust and unambiguous way the existence of the universal critical behavior of the metal-insulator transition studied here.

V. CHARGE STRUCTURE FACTOR
Let us now investigate the long wavelength limit or equivalently the small |q| behavior of the static charge structure factor N (q). Being the static structure factor N (q) the integral over all the frequencies of the dynamical structure factor N (q, ω), i.e. N (q) = dωN (q, ω), it depends, e.g., as in the spatial dimensionality D > 1 within Fermi liquid theory [52,53], on the charge and the Fermi velocities, as well as on other low energy parameters such as the Landau parameter F s 0 in the standard Fermi liquid theory. This is therefore an interesting quantity and gives us information on how the dynamics evolves as we approach the metal-insulator transition point at U c from the metallic side. The static charge structure factor N (q) at momentum q is defined as where N U is the number of unit cells, n j = α s c † jαs c jαs , and c † jαs is the creation operator of electron at unit cell r j and sublattice α (= A, B) with spin s (=↑, ↓). In the metallic region for U < U c , the charge structure factor should behave as the non-interacting one, i.e., for small |q|, where α is a suitable constant that can be renormalized by the interaction, as in Fermi liquid theory [52]. Figure 13 shows the ratio of the charge structure factor for finite U and the noninteracting limit, denoted as N U=0 (q), at the smallest non-zero momentum q * available for a given system size.
Since N (q) ≃ |q| 2 in the insulating phase [54][55][56], we expect that R ∼ 1/ ln L in the insulating phase, whereas R ∼ const. in the metallic phase. Therefore, there is a change of behavior in R across the critical value of U c . Indeed, we find in Fig. 13 that R decreases with increasing L for U larger than U c . It should be emphasized that R remains finite as we approach the transition point from the metallic side. As shown in Fig. 14, R as a function of U/t for different system sizes crosses around the critical point at a finite value of R, which indicates that the coefficient α in Eq. (15) is neither singular nor critical at the critical point for U → U c , but approaches a finite constant in the thermodynamic limit. On the other hand, this coefficient α vanishes at the metal-insulator transition within the Gutwiller [7] or Brinkman-Rice [8] approximation on any lattice (including the honeycomb lattice) as a result of the vanishing of the double occupancyd = n i↑ n i↓ and the sum rule valid at half-filling |q| =0 where g NN is the nearest neighbor density correlation, i.e., g NN = n jA n jB − n jA n jB , and n jα = s c † jαs c jαs . This is simply because g NN → 0 when d → 0 at the metal-insulator transition described by the Gutzwiller approximation [57], and is therefore in sharp contrast with our results. FIG. 14. Charge structure factors N (q) at the smallest nonzero momentum |q| = q * divided by the the charge structure factor NU=0(q * ) in the non-interacting limit for the honeycomb lattice model with different system sizes indicated in the figure. The dashed line denotes the critical Uc determined by the data collapse fit of the staggered magnetization in Fig. 9(a).
Although our results do not represent a direct evidence for the non-singular behavior of v F close to the metalinsulator transition, v F is certainly related to the charge velocity, as in Fermi liquid theory, which in turn affects the value of α. Therefore, our results imply that v F remains finite at the metal-insulator transition. In this respect, it should be noted that the evolution of α as a function of U found in Fig. 14 is compatible with the expected behavior of the Dirac Fermi velocity v F for the GN universality class of the metal-insulator transition [34]. As discussed previously [36], the vanishing of the quasiparticle weight Z without a renormalization of v F is understood as a consequence of an equal divergence in the momentum k and frequency ω derivatives of the electron self-energy Σ(k, ω) at the Dirac point. The quasiparticle weight Z at the Dirac point with momentum k F is given as Therefore, in order to compensate the divergence of ∂ ∂ω ReΣ(k F , ω) ω=0 , i.e., Z → 0, at the metal-insulator transition point, ∂ ∂k ReΣ(k, 0) k=kF must diverge at the same time. This implies that the strong momentum dependence of Σ(k, ω) around the Dirac point, not included in the simplest DMFT approach, is an essential ingredient to describe the metal-insulator transition. Similar arguments are also found in earlier studies of the t-J model with the large-N expansion [58]. It is also interesting to note that the nodal Fermi velocity remains finite in the carrier number controlled Mott metal-insulator transition described by a Gutzwiller projected d-wave BCS state, which exhibits the massless Dirac dispersion at the nodal point [59].
In any event, our results are certainly useful to characterize the metal-insulator transition in the charge sector. Indeed, the critical value U c in the charge sector can be estimated directly by this simple analysis shown in Fig. 14 without performing rather elaborated finite size scaling of quantities such as the charge gap [12,19], which is usually very time consuming and difficult to compute with high accuracy.

VI. DISCUSSION AND CONCLUSIONS
The GN models have been extensively studied in quantum field theory [23]. In the standard field theoretical treatment of the transition, the critical behavior is studied in space-time dimension d = D + 1, where D is the spatial dimensionality. The critical exponents of the GN models have been evaluated by several standard techniques, such as the large N expansion [60,61] and the ǫ-expansion around the lower d = 2 + ǫ [62] or upper d = 4 − ǫ [31] critical dimensions, and these are summarized in Table I. It is clearly noticed in Table I that there exist sizable discrepancies among the critical exponents calculated by these different analytical techniques. Therefore, an unbiased numerical study is highly desired to clarify the critical behavior of the GN models.
In this paper, we have established, based on robust and reliable numerical simulations on fairly large clusters, the universal properties of the metal-insulator transition for the two different lattice models of interacting Dirac electrons in two spatial dimensions. We have determined the critical exponents which characterize the universal quantum critical behavior in both metallic and insulating phases at the vicinity of the transition for these models. Since it is expected that the effective low-energy theory of these lattice models are described in the continuous limit by the chiral Heisenberg universality class of the GN model with N = 8, our study represents currently the most accurate and reliable results also for this fundamental GN model to reveal the universal critical behavior. Indeed, our results resolve some of the inconsistency among different approximate approaches for the GN models shown in Table I, especially evident for N = 8.
We have also clarified how the quasiparticle in the SM phase collapses as the AF insulating phase is approached with increasing U for these two models. We have shown that the quasiparticle weight Z diminishes and becomes zero at the metal-insulator transition and found that the exponent η Z = νη ψ , characterizing the renormalization of the quasiparticle weight Z [see Eq. (10)], is much smaller for both models (η Z ≃ 0.2) than the one predicted by simple mean-field and dynamical mean-field correlated theories of the metal-insulator transition, for which η Z = 1 [8,10]. More interestingly, we have also found that small q limit of the static charge structure factor is not singular at the metal-insulator transition, suggesting that the Dirac Fermi velocity v F is not critical at the transition. These critical behaviors, small η Z and finite v F , are qualitatively the same as the ones expected for the GN universality class of the metal-insulator transition [34]. Therefore, our results provide a clear numerical evidence that the critical behavior of the GN model applies also to lattice models and describes correctly the metal-insulator transition in interacting Dirac electrons. It should be noted here that the non-criticality of v F for the honeycomb lattice model is also found by quantum cluster methods [35,36], although it was first overlooked within the single-site DMFT [33].
It should also be remarked that, strictly speaking, the metal-insulator transition studied here does not describe a genuine "Mott transition" because the insulating phase breaks SU (2) symmetry for U > U c . Indeed, as discussed below, each type of possible symmetry breaking in interacting Dirac fermions determines a different universality class of the transition characterized with different critical exponents. Nevertheless, the metal-insulator transition studied here does not originate from a Slater-type nesting instability at weak coupling since the density of states is zero at the Fermi level. The transition instead occurs at an intermediate or strong coupling region as evidenced by the fact that U c ≃ 5.55t for the π-flux model is almost the same as the non-interacting band width 4 √ 2t. Very recently, using the continuous time QMC, Wang et al. [68] have studied similar models of spinless fermions with the nearest neighbor repulsion V , where, with increasing V , a transition from the SM to a staggered charge-density-wave (CDW) state occurs. With a finite size scaling analysis based on the results for clusters up to 450 sites, they have estimated the critical exponents, ν TABLE I. Critical exponents, ν, β, and η ψ , of the interaction-driven phase transition in interacting Dirac fermions in d = 2 + 1 for the lattice models (honeycomb lattice and π-flux models) and the effective continuous models (Gross-Neveu models) with the total number N of fermion components. Different classes correspond to different symmetries broken in the ordered phases. Numbers in parentheses for ν, β, and η ψ indicate statistical errors in the last digits. For comparison, the critical exponents for other related models with N = 4 and 8, belonging to different universality classes, are also listed. FRG stands for the functional renormalization group method. and β, for the CDW order parameter [68]. These critical exponents have been recently revisited by the Majorana QMC method on larger clusters up to 1152 sites [69]. As seen in Table I, these exponents are clearly different from those for the spinful models that we have studied here. This is simply understood because these spinless and spinful models belong to different universality classes. Indeed, it is known that the spinless lattice models with the nearest neighbor interaction at half-filling are described in the continuous limit by the GN model with N = 4 and the chiral Ising universality class. It is also interesting to notice in Table I that the critical exponents for the spinless lattice models estimated numerically are rather different from the analytical results. This also demonstrates that numerically exact studies are highly valuable to accurately determine the quantum criticality and to remove the ambiguities that might arise from inadequate approximations. It should be emphasized that the various universality classes depend on the physics, namely, the symmetry that is broken in the ordered phase and the total number N of Dirac fermion components that describe the corresponding critical theory [64]. The well explored universality classes for the GN models in the continuous limit include the following three classes: 1. Chiral Ising universality class. Z 2 symmetry, i.e., a discrete order parameter is broken, for instance, when a commensurate CDW order settles.
2. Chiral XY universality class. U (1) symmetry is broken and the order parameter is characterized by an angle, as in a superconducting state.
3. Chiral Heisenberg universality class. SU (2) symmetry is broken. This should occur in the transition studied here, as the order parameter-the staggered magnetization-is characterized by a vector with three components [SU (2) is equivalent to SO (3)].
Among these, the two classes have been studied so far based on the lattice models with unbiased numerical techniques: the chiral Ising universality class by Wang et al. [68] and Li et al. [69], and the chiral Heisenberg universality class studied here and in Ref. [18]. We expect that the quantum criticality belonging to the chiral XY universality class emerges in a negative U Hubbard model with the Dirac points of the non-interacting energy dispersion at the Fermi level, provided that the SU (2) symmetry which relates the CDW order to the superconducting one is not satisfied (otherwise the chiral Heisenberg universality class applies again). For example, by adding the next nearest-neighbor hopping t ′ in the same lattice models studied here but with a negative U (no sign problem with t ′ in the negative U Hubbard model), the chiral XY universality class with different critical exponents can be investigated in the same unbiased numerical approach.
In conclusions, we have investigated the critical behaviors of the metal-insulator transition in the interacting Dirac electrons in two spatial dimensions, described by the two different lattice models, the Hubbard model on the honeycomb lattice and the π-flux Hubbard model on the square lattice, at half-filling. We have performed the unprecedentedly large-scale quantum Monte Carlo simulations to systematically calculate the staggered magnetization and the momentum distribution function. The calculation of the momentum distribution function is particularly important because it allows us to examine the quasiparticle weight and thus explore the critical behavior also in the metallic side, which had never been successful previously in the unbiased numerically exact calculations. The ground state phase diagrams determined by extrapolating the staggered magnetization and the quasiparticle weight to the thermodynamics limit have revealed the continuous nature of the transition between the SM and the AF insulator with no intermediate phase. Therefore, our results firmly rule out the possibility of a spin liquid phase for these two models proposed in the earlier studies. We have obtained the critical exponents by the careful and accurate finite-size scaling analysis and found that the two lattice models belong to the same universality class. Since the low-energy effective model in the continuous limit for these two models is described by the GN model with N = 8 and the chiral Heisenberg universality class, our results represent currently the most accurate determination of the quantum criticality of this universality class. Finally, we have shown that the quasiparticle weight monotonically decreases with increasing U and becomes zero at the metal-insulator transition, while the Fermi velocity seems non-critical. This qualitatively important feature is indeed in good agreement with the one expected for the GN universality class of the metal-insulator transition and cannot be captured by the simple mean-field or Gutzwiller-type approximate argument.
at momentum k with the energy ε k = ±|h k |. Here for the honeycomb lattice model and for the π-flux model, where τ 1 and τ 2 are the primitive translational vectors defined in Figs. 1(a) and 1(b). In principle, in order to determine the quasiparticle weight Z, the occupation number should be calculated in terms of the natural orbitals, i.e. the eigenvectors of the density matrix. To this end, we should note that the density matrix at momentum k is a 2 × 2 matrix and is represented as Therefore, the "dressed" quasiparticle operators simply readψ with the occupation ψ † k,±,sψ k,±,s = Notice that f k = − 1 2 h * k |h k | in the non-interacting limit where the bonding (anti-bonding) states are all occupied (empty). We have verified that even in the interacting case studied here the natural orbitalsψ k,±,s are almost indistinguishable from the non-interacting bonding and anti-bonding states ψ k,±,s Indeed, we have found in Fig. 15 that the difference between the quasiparticle weight Z calculated with the natural orbitals and the one with the non-interacting bonding and anti-bonding states is negligible for the models studied here both in the metallic and insulating phases. Therefore, the treatment for computing the quasiparticle weight Z in Sec. III B and Sec. IV B, i.e, the jump of the energy resolved momentum distribution function n(ε k ) given in Eq. (A3) at the Fermi level, is not only asymptotically valid to determine the corresponding critical exponent but also represents a good quantitative estimate of Z.
Appendix B: Effects of the leading correction term in the finite-size scaling analyses Here, we examine the robustness of the fitting parameters in Eqs. (11) and (12) with and without the leading correction term in the finite-size scaling analyses. Since the leading correction term is expected less important for sufficiently large clusters, we examine the system size dependence of the fitting parameters. Tables II and III summarize the fitting parameters in the staggered magnetization, which are obtained for clusters including the smallest lattice size L min in the data collapse without and with the leading correction term in the finite-size scaling ansatz of Eq. (11), respectively. These results are also compared in Fig. 16 as a function of L −1 min . It is noticed that U c /t systematically increases as the data sets with smaller L are removed when we use the finite-size scaling ansatz without the leading correction term, i.e., c = 0 in Eq. (11). Accordingly, the critical exponent β tends to decrease, although the critical exponent ν is less affected by including or not including the leading correction term in the finite-size scaling ansatz, at least, within the statistical errors. These systematic differences with varying L min imply that the leading correction term in the finite-size scaling ansatz is not negligible.
On the other hand, when the leading correction term is included, the results are robust against the choice of L min as shown in Table III and Fig. 16. In addition, U c and β evaluated in the data collapse analysis are statistically consistent with those obtained by the direct fit of thermodynamically extrapolated staggered magnetization, shown in Fig. 4. In Sec. IV A, we report the results of the data collapses with the leading correction term, because they are clearly more accurate and stable. The quality of our data and extrapolations is further supported by the fact that both results obtained with and without the leading correction term tend to be identical when the system sizes are large enough, as shown in Tables II and III, and in Fig. 16.
The finite-size scaling ansatz for the jump of n(ε k ) including the leading correction term is given as ∆n(u, L) = L −η ψ 1 + dL −ω ′ f n (uL 1/ν ), where d and ω ′ are additional fitting parameters. The obtained critical exponents η ψ for various choices of L min are summarized in Tables IV and V, and also in Fig. 17. We find that the stable solutions with ω ′ > 0 can not be obtained for small L min . Moreover, the estimated ω ′ tends to increase for larger L min (see Table V) and it is not possible to reach a converged value of ω ′ within the given statistical accuracy and the system sizes available. Nevertheless, we confirm that the estimated values of η ψ with the correction term are instead converged and fully consistent with those obtained without the correction term for large enough L min , as shown in Fig. 17. Therefore, we only show in Sec. IV B the results of the data collapse analysis without the leading correction term. Results of the critical points Uc/t and the critical exponents, ν and β, obtained from collapsing data of the staggered magnetization, ms(u, L), without the leading correction term, i.e., c = 0 in Eq. (11), for the honeycomb lattice model (upper rows) and the π-flux model (lower rows). Lmin refers to the smallest L used in the data collapse. The maximum L is 36 for the honeycomb lattice model and 40 for the π-flux model.