## Abstract

Tang *et al*. (Research Articles, 10 August 2018, p. 570) report on the properties of Dirac fermions with both on-site and Coulomb interactions. The substantial decrease, up to ~40%, of the Fermi velocity of Dirac fermions with on-site interaction is inconsistent with the numerical data near the Gross-Neveu quantum critical point. This results from an inappropriate finite-size extrapolation.

The low-energy excitations of many condensed matter systems, such as electrons on the honeycomb lattice of graphene, can be described by massless Dirac fermions with a Dirac cone-like dispersion relation and a corresponding Fermi velocity. The inclusion of interactions among the fermions eventually leads to a breakdown of this description, once the system undergoes a quantum phase transition to an insulating phase beyond a critical interaction strength. Below this interaction-induced quantum critical point (QCP), the system is characterized by massless Dirac fermions with a renormalized Fermi velocity. The quantification of this velocity renormalization constitutes a challenge in numerical simulations: Crossover effects strongly alter finite-size system estimates close to critical points, and a careful analysis of the actual excitation energies is required to extract reliable results.

Tang *et al*. (*1*) extract the momentum-resolved one-particle excitation energies from imaginary-time correlation functions obtained by projective quantum Monte Carlo (QMC) simulations. Upon approaching the Dirac points, the lattice dispersion of the noninteracting (tight-binding) fermion system takes on a linear, relativistic form that defines the tight-binding Fermi velocity *v*_{0} at the Dirac point. The inclusion of either on-site (Hubbard) interactions or extended Coulomb interactions leads to changes of these excitation energies. Below the interaction-induced Gross-Neveu QCP, the dispersion remains gapless at the Dirac point in the thermodynamic limit (TDL) at infinite lattice size, defining the semimetallic (SM) regime. For the case of the Hubbard model, the Gross-Neveu QCP is known to be located at an on-site repulsion of *U*_{c}(γ = 0) = 3.85(2)*t*, beyond which the model exhibits antiferromagnetic order (*2*). Here, *t* denotes the nearest-neighbor hopping strength on the honeycomb lattice, and γ = 3α_{0}/*U* in terms of the Coulomb interaction strength α_{0}. Throughout this comment, we follow the notation used in (*1*).

In order to extract the interaction-induced renormalization of the Fermi velocity within the SM phase, the excitation gaps obtained from the QMC data for finite-size systems need to be extrapolated to the TDL. Finite-size effects are observed in all excitation energies, but in particular close to the QCP. This is seen in Fig. 1, which shows the bare finite-size excitation gaps, extracted from the imaginary-time QMC data as detailed in the supplementary materials of (*1*), based on the datasets made available online by the authors of (*1*). We observe that the finite-size effects are most pronounced at the Dirac points themselves (Fig. 1), where the gap vanishes in the TDL within the SM regime for *U < U*_{c}(0) and at the Gross-Neveu QCP *U* = *U*_{c}(0). On the other hand, for momenta in the immediate vicinity of the Dirac points, the finite-size effects are seen to be much weaker (Fig. 1), and one may estimate the TDL values of the excitation energies at these momenta from the values on the largest system sizes accessed in (*1*).

In Fig. 1 we also include data provided by Tang *et al*., showing their finite-size extrapolated gaps. This processed data (based on the interpolation scheme used in their figure S2) are seen to be incompatible with the behavior of the excitation energies for small values of *a*Δ*k* extracted with our scheme. Moreover, as shown in Fig. 2A, the excitation energies close to, but excluding, the Dirac point exhibit only a weak dependence on *U*. Thus, for γ = 0, the low-energy Dirac dispersion, and hence the Fermi velocity, is in fact only weakly modified by the on-site interactions. In particular, the low-energy dispersion traced by our data in Fig. 1 for *U* = 3.75*t* is clearly inconsistent with the ~40% decrease of the Fermi velocity from *v*_{0} reported in (*1*), which is indicated by the lower red line in Fig. 1.

A reliable estimate for the Fermi velocity at the Dirac point for values of *U* inside the SM regime can be obtained from a finite-size analysis of the rescaled lowest particle-excitation energy *E*/(*a*Δ*k*) at the closest momentum to the Dirac point on each finite lattice. The corresponding finite-size values are compared to *v*_{0} in Fig. 2B, and they demonstrate a remarkably weak renormalization of the Fermi velocity throughout the SM phase. A reduction by ~40% from the value *v*_{0} is not compatible with the observed steady approach of *E*/(*a*Δ*k*) toward *v*_{0} with increasing system size for all considered values of *U* within the SM regime.

The substantial overestimation of the Fermi velocity suppression by the on-site interaction reported in (*1*) (see also Fig. 2C) is in fact due to an inappropriate finite-size extrapolation procedure, which is documented in figure S2 of (*1*): The authors of (*1*) use the slope between the finite-size excitation energies at the Dirac point and the closest point to the Dirac point [with a linear interpolation to the simulation scale] as estimator. The finite-size energies at the Dirac point suffer from particularly large finite-size effects near the Gross-Neveu QCP, and the strong suppression of the Fermi velocity that is reported in (*1*) near the Gross-Neveu QCP merely reflects the enhanced finite-size effects of the excitation energy at the Dirac point, but not the renormalization of the actual low-energy dispersion. The extraction of velocities based on the softest excitations is also reported to be subtle for related quantum phase transitions [see, e.g., (*3*–*5*)].

Their means of data analysis therefore did not allow the authors of (*1*) to faithfully reproduce the Fermi velocity renormalization beyond the weak-coupling regime. The Fermi velocity renormalization shown in figure 2 of (*1*) is affected strongly by their finite-size analysis scheme, in particular in the vicinity of the Gross-Neveu QCP at *U*_{c}(γ), which calls for a revised analysis and interpretation of the numerical data along the lines outlined in this comment.

**Acknowledgments:**We thank H.-K. Tang and colleagues for making their data openly available.

**Funding:**Supported by FWF projects I-2868-N27 and F4018 and by DFG projects RTG 1995 and FOR 1807.

**Author contributions:**S.H., T.C.L., and M.S. performed the data analyses and prepared the figures; S.W. and A.M.L. directed the investigation; the manuscript reflects the contributions of all authors.

**Competing interests:**The authors declare no competing interests.

**Data and materials availability:**Data and computer scripts are available at Harvard Dataverse (

*6*).