## Computers tease out interaction effects

Although graphene is often thought of as a material in which electron-electron interactions are negligible, some of its properties cannot be explained by such a simple picture. Tang *et al.* undertook comprehensive quantum Monte Carlo numerical calculations that consider both long-range and contact interactions in systems that, like graphene, have two-dimensional (2D) Dirac electrons. Different 2D Dirac materials systems, such as topological insulators and graphene on various substrates, reside in different parts of the resulting phase diagram.

*Science*, this issue p. 570

## Abstract

The role of electron-electron interactions in two-dimensional Dirac fermion systems remains enigmatic. Using a combination of nonperturbative numerical and analytical techniques that incorporate both the contact and long-range parts of the Coulomb interaction, we identify the two previously discussed regimes: a Gross-Neveu transition to a strongly correlated Mott insulator and a semimetallic state with a logarithmically diverging Fermi velocity accurately described by the random phase approximation. We predict that experimental realizations of Dirac fermions span this crossover and that this determines whether the Fermi velocity is increased or decreased by interactions. We explain several long-standing mysteries, including why the observed Fermi velocity in graphene is consistently about 20% larger than values obtained from ab initio calculations and why graphene on different substrates shows different behaviors.

In 1952, Freeman Dyson made the argument that all theoretical solutions to problems in quantum electrodynamics have the form of a perturbative asymptotic series expansion in the fine-structure constant α, where the Coulomb potential is of the form α/*r*, with *r* being the distance between two electrons. Furthermore, he showed that such solutions are uncontrolled beyond the order of perturbation theory given by the inverse fine-structure constant (*1*). For condensed matter realizations of two-dimensional (2D) Dirac fermions, the equivalent of the fine-structure constant is the long-range Coulomb coupling constant α ~ 1, which implies that any perturbative theory is potentially uncontrolled (*2*). Indeed, first-order perturbation theory for such systems gives rather peculiar results. To first order in α, the inverse coupling constant (proportional to the Fermi velocity) itself diverges both in the infrared (long distances) and ultraviolet (small distances). Introducing a lattice scale fixes the ultraviolet divergence but not the infrared one (*3*, *4*).

This divergence is what lead Ye and Sachdev (*5*) to describe the effects of long-range electron-electron interactions as “dangerously irrelevant”: Although the system flows under the renormalization group to a noninteracting theory, physical observables are strongly renormalized by the Coulomb interaction. Since then, dozens of theoretical works (*6*) have confirmed this basic picture: In the absence of disorder, and precisely at half filling, the role of long-range Coulomb interactions is to renormalize the electron Fermi velocity to infinity (limited only by the speed of light, if one includes a dynamical interaction).

In a parallel development, the Hubbard model on a honeycomb lattice provided a low-energy realization of Dirac fermions interacting through a short-range contact interaction *U*. Increasing the short-range interaction results in a quantum phase transition at a critical value *U* = *U _{c}* from a semimetal to an antiferromagnetic Mott insulator. This phase transition was predicted to be of the Gross-Neveu universality class (

*7*), which was recently confirmed numerically (

*8*,

*9*). In diagrammatic perturbation theory, typically both the Fermi velocity and quasiparticle residue vanish at a quantum phase transition. However, for the Gross-Neveu critical point, the Fermi velocity remains finite (despite the vanishing of the quasiparticle residue), and this suppression of the Fermi velocity to a finite value has been observed numerically (

*10*). Renormalization group studies (

*11*) have shown that the Fermi velocity is not modified by weak short-range interactions, and, in our numerics below, we verify all these features for the Hubbard model on a honeycomb lattice. However, because 2D Dirac fermions are unable to screen the long-range Coulomb potential, it is widely believed (

*6*,

*12*) that this pure on-site Hubbard model has limited applicability to experiments done in real materials.

Experimentally, 2D Dirac fermions can be realized in a variety of condensed matter systems, including on the surfaces of 3D topological insulators (*13*, *14*) and in artificial graphene made from quantum corrals of carbon monoxide arranged in a honeycomb lattice on a copper substrate (*15*), as well as in other systems (*16*). For concreteness, we focus our attention on graphene, the most studied and versatile realization of 2D Dirac fermions. Experiments have been unable to realize the precise configuration necessary to probe this strange interacting metallic state that features electron quasiparticles moving at the speed of light, despite their quasiparticle character smearing away; however, several probes of ultraclean graphene—including magnetotransport (*17*), infrared spectroscopy (*18*), capacitance (*19*), angle-resolved photoemission spectroscopy (*20*), tunneling spectra (*21*), and Raman scattering (*22*)—all reveal a clear breakdown of the noninteracting theory.

In this work, we address the competing effects of short-range and long-range parts of any realistic model of the Coulomb interaction. We use a nonperturbative, numerically exact, projective quantum Monte Carlo method to study the evolution of physical observables in a controllable manner. We find that, in the regime dominated by long-range interactions, there is an enhancement of Fermi velocity consistent with perturbation theory. Conversely, close to the phase transition dominated by short-range interactions, we find a suppression of Fermi velocity and a collapse of the numerical data for different values of the ratio between long-range and short-range interactions onto one curve. Our numerical results interpolate between these two limits and are valid for all interaction strengths, but are constrained by the finite system sizes in our simulations. Therefore, we use a renormalization group scheme to extrapolate the quantum Monte Carlo results to experimentally relevant energy scales, where we predict that observables will depend on both the short-range and long-range components of the Coulomb interaction as well as the energy scale of the observation (as we explain below, all these parameters can be tuned in current experiments). Moreover, we find that the lattice scale not only regularizes the ultraviolet divergence of the Fermi velocity but can also make the infrared divergence unnoticeable in the experimental window.

## Theoretical model

To accomplish this, we study interacting fermions on a honeycomb lattice with competing short-range and long-range interactions. The model is described by the Hamiltonian where h.c. is the Hermitian conjugate, *t* is the tight-binding hopping, the second-quantized operator () creates (annihilates) an electron of spin at position **r*** _{i}*, and gives the electron density at position

**r**

*. The interaction consists of a short-range part acting between electrons on the same site with different spins, , and a long-range part depending on the distance between the electrons as , where*

_{i}*a*is the lattice constant. We define γ = 3α

_{0}/

*U*as the ratio between the long-range and short-range components. The inclusion of the long-range component was made possible by recent developments in lattice quantum chromodynamics (

*23*), which we adapt for our purposes here (

*24*).

Including electron-electron interactions can do one of two things: The Dirac fermions can remain metallic but with a modified Fermi velocity or the interactions can gap the system, giving a Mott insulator. We use our quantum Monte Carlo method to map out the phase transition between the semimetallic phase and the Mott insulating phase by calculating the antiferromagnetic structure factor [see eq. S7 in (*16*)], which, in the thermodynamic limit, is finite for the Mott insulator and vanishing in the metallic phase. For any given strength of the long-range interaction, the structure factor shows a unique crossing point *U _{c}*(α

_{0}) for various system sizes [see fig. S3 in (

*16*) for a representative example with γ = 0.5 and system size (

*L*) = 6, 9, 12, and 15]. By tracing this crossing point for 10 different choices of γ, we can map out the phase diagram of this model shown in Fig. 1, highlighting the competing effects of the short-range and long-range parts of the Coulomb potential. We confirm an earlier finding (

*25*) that the critical point remains in the Gross-Neveu universality class, even with the presence of the long-range Coulomb interaction. It is commonly believed (

*6*,

*12*) that increasing the strength of the Coulomb interaction will favor the Mott insulating phase. By contrast, we find that reducing the relative role of Coulomb interaction versus on-site interaction, for example, through dielectric screening (

*26*) or biaxial strain (

*24*), provides the most favorable route for realizing antiferromagnetic order in graphene experiments. Renormalization group calculations in dimensions are the preferred method to understand such phase transitions (

*27*). In section 2.2 of (

*16*), we extend this technique to calculate the spin-full Gross-Neveu model, from which the red curve in Fig. 1 is obtained. This gives us two intuitive ways of thinking about why the phase transition shifts to the right with increasing long-range interactions. First, as one increases the long-range Coulomb tail, the effective on-site potential decreases. (One can think of this, qualitatively, as the difference between the on-site potential and the nearest-neighbor potential.) Therefore, with the inclusion of the long-range piece, one needs a larger on-site potential to get the same effective critical Hubbard potential (

*28*). Second, although the Hubbard potential favors an antiferromagnetic ground state, the nearest-neighbor potential favors instead a charge density wave ground state. By including the long-range piece, one needs a larger on-site potential to stabilize the antiferromagnetic phase.

Figure 2 shows our main numerical results on the Fermi velocity, which is the defining property of the massless Dirac spectrum. We plot the renormalization of Fermi velocity (*v _{F}* −

*v*

_{non})/

*v*

_{0}against

*U*/

*U*(α

_{c}_{0}), which is the short-range component of the interactions normalized by the critical value obtained from the phase diagram in Fig. 1. Here

*v*is the interaction-renormalized Fermi velocity obtained in quantum Monte Carlo at the simulation scale , where

_{F}*k*is the momentum,

*v*

_{non}is the tight-binding Fermi velocity at the simulation scale, and

*v*

_{0}is the tight-binding Fermi velocity at the Dirac point. Each data point in Fig. 2 is an extrapolation to the thermodynamic limit from four lattice sizes, each with an average of ~100,000 quantum Monte Carlo sweeps. The interacting Fermi velocity at the simulation scale is obtained by first determining the convergent ground state value of the unequal time Green function

*G*(τ) for large τ, where the single exponential decay time, determines the first excitation energy of the system [defined here as the energy difference between the ground state (GS) of

_{k}*N*+ 1 fermions with a total momentum

*k*and the ground state of

*N*fermions with zero total momentum]. Full details of our quantum Monte Carlo scheme (figs. S1 to S3), datasets, the analysis, and the discussion on the use of twisted boundary conditions (figs. S4 to S7) are provided in (

*16*).

Our numerical data show that the Fermi velocity renormalization has notably different behavior for (which we call the “weak-coupling regime”) and , which is in the vicinity of the Gross-Neveu critical point. In the weak-coupling regime, we observe an increase in the Fermi velocity, and all the quantum Monte Carlo data for different ratios of short-range and long-range components γ collapse when plotted as a function of the long-range interaction α_{0} (see Fig. 2, left inset). By contrast, in the vicinity of the Gross-Neveu critical point, the Hubbard model (α_{0} = 0) shows a 40% decrease in Fermi velocity. Even after including the long-range component of the Coulomb interaction, by subtracting the intercept of the Fermi velocity at the Gross-Neveu critical point, all the numerical data collapse to the Hubbard model function form (see Fig. 2, right inset). This shows that interacting fermions on a honeycomb lattice are governed by two very different fixed points: one controlled by the long-range interaction, giving an enhancement in Fermi velocity, and the other governed by the short-range interaction, giving a suppression of the Fermi velocity. As discussed in table S2, estimates for the realistic Coulomb potential in graphene place it in the crossover between these two regimes, whereas topological insulator Bi_{2}Se_{3} is close to the Mott transition and artificial graphene using quantum corrals is in the weak-coupling regime (*16*). Our numerical results span the full crossover between these two regimes.

The emergence of a stable weak-coupling fixed point and an unstable Gross-Neveu fixed point is anticipated by renormalization group studies. In the weak-coupling regime, all the quantum Monte Carlo data (Fig. 2, lines plotted on the left-hand side of the main panel) can be reproduced by using a one-parameter random phase approximation (RPA) theory , where Λ_{s} is our numerical scale and Λ = 6.2 ± 0.2 is obtained from fitting the data for small α. The RPA functions *F*_{1}(α_{0}) and *F*_{0}(α_{0}) have been rederived several times in the literature (*4*, *5*, *16*, *29*). Our nonperturbative numerics verify this functional dependence on α_{0}. This agreement with the nontrivial functional dependence on α_{0} gives us confidence to trust the logarithmic dependence on scale also predicted by the RPA calculation. Our full numerical dataset is consistent with this logarithmic increase in the weak-coupling limit [*P* < 10^{−3}; see figs. S8 and S9 in (*16*) for the full dataset and additional evidence].

Although there is currently no analytical theory for the dependence on *U*/*U _{c}* close to the Gross-Neveu fixed point, we can adequately describe the quantum Monte Carlo data by using a phenomenological model that has three parameters for the pure Hubbard model data and one additional parameter for linear dependence on γ = 3α

_{0}/

*U*of the Fermi velocity at the Gross-Neveu critical point

*v*

_{GN}(γ). Defining the change in interaction strength from the critical value (ε) = 1 −

*U*/

*U*(α

_{c}_{0}), we find (

*16*) that (

*v*−

_{F}*v*

_{non})/

*v*

_{0}=

*C*

_{0}+

*C*

_{1}ε +

*C*

_{2}ε

^{2}+

*m*

_{γ}γ with

*C*

_{0}= −0.384 ± 0.002,

*C*

_{1}= 1.35 ± 0.05, and

*C*

_{2}= −1.2 ± 0.1 determined just from the Hubbard; the single additional parameter

*m*

_{γ}= 0.333 ± 0.005 captures the effects of the long-range Coulomb potential. This fit is shown as the lines plotted on the right-hand side of Fig. 2. In section 5 of (

*16*), we describe the suppression of the Fermi velocity at the Gross-Neveu fixed point using an approximate spin-wave theory with short-range interactions, giving the brown star in Fig. 2. Approaching the Gross-Neveu transition, antiferromagnetic fluctuations become progressively slower and couple to the Dirac fermions, thereby reducing the Fermi velocity (

*30*).

## Flow beyond numerical scales

The analysis so far has been at the simulation scale Λ_{s} = 0.48 that is limited by the largest system size that we can simulate. However, the experimental scale is set mostly by the degree of disorder in the system (currently, energy scales as small as Λ* _{k}* = 10

^{−3}can be measured). How, then, do we extrapolate our numerical findings to the experimental regime? Although we cannot numerically probe scales smaller than Λ

_{s}, we can probe larger energy scales, thereby providing the inputs for a renormalization group flow from the numerical scale to the experimental scale. This procedure is illustrated in Fig. 3, which we now explain.

Figure 3A shows results for *U* = 0.1 in the weak-coupling regime. For the pure Hubbard model (α_{0} = γ = 0), there is no renormalization of the Fermi velocity, even at larger energy scales (blue circles). For small α_{0} (shown in Fig. 3A as red data points for quantum Monte Carlo results, with γ = 3α_{0}/*U* = 0.1 as a representative example), the values of the numerical data increase weakly with decreasing Λ* _{k}*. To show that this increase is consistent with a logarithm, we show two theoretical analyses. The red curve shows the first-order perturbation theory of a Dirac spectrum with the same single global parameter discussed already in Fig 2. We attribute the disagreement with the quantum Monte Carlo data to the fact that, at such large energy scales, the lattice model has nonlinear terms. To check this, we also solve numerically the first-order perturbation theory in α

_{0}on a finite lattice [see section 6 of (

*16*)]. Because the lattice is specified, there is no adjustable parameter in this calculation. This is shown as the black curve in Fig. 3A, and it agrees with the quantum Monte Carlo at large energy scales (reduced χ

^{2}= 1.96), giving us confidence in our analysis. In this method, we can go to lattice sizes as large as

*L*= 1500, and, in the small Λ

*window, the lattice perturbation theory results look similar to the logarithmic divergence in the continuum perturbation theory. Additional numerical results would be needed to confirm that this observed logarithmic dependence in our quantum Monte Carlo data persists all the way to the Dirac point.*

_{k}We now turn to the Gross-Neveu critical point. Figure 3B shows quantum Monte Carlo data for the Hubbard model (blue circles) and a representative example γ = 0.5 (red circles). For only contact interactions γ = 0, the renormalized Fermi velocity decreases with decreasing Λ* _{k}*, with a kink at Λ

*~ 1. There is no available theory for this behavior. The blue dashed line shows a phenomenological three-parameter fitting function used to capture this numerical finding (*

_{k}*16*). Including the long-range component, the dominant effect is an overall upward shift in the curve (as already discussed in Fig. 2). However, at large Λ

*, we notice another subtle difference. The Hubbard model always shows the renormalized velocity to be a monotonically decreasing function with decreasing Λ*

_{k}*, but with long-range interactions, the dependence is nonmonotonic. This observation suggests that the logarithmic increase in Fermi velocity persists into the strong-coupling regime, but the weak increase is masked by the strong decrease induced by the short-range Hubbard interaction. To check this, we obtain the black curve in Fig. 3B by adding the Hubbard results (blue dashed line in Fig. 3B) to the parameter-free lattice perturbation theory results (black line in Fig. 3A) shifted linearly with γ just like we did for the black lines on the right-hand side of Fig. 2, but here, instead of using the lowest-energy quantum Monte Carlo results to do this shift, we use the highest-energy data. Details of the fitting procedure are explained in section 8 of (*

_{k}*16*). The key insight here is that, although the effect of the logarithm is not observable at Λ

_{s}, it can be extracted at larger values of Λ

*.*

_{k}In Fig. 3B, the increase in the quantum Monte Carlo data (red circles) with decreasing Λ* _{k}*—associated with the infrared divergence—is weaker than that predicted by the lattice perturbation theory (black curve). This weaker divergence is also anticipated by looking at the structure of the terms in a perturbation theory expansion about the Gross-Neveu fixed point [see eq. S32 of (

*16*)]. These insights suggest fitting for the logarithm at large Λ

*. Figure 3C shows a blowup of the region around Λ*

_{k}*= 2, where a single adjustable parameter Λ*

_{k}_{RHS}= 70 ± 20 captures this increase. The subscript RHS (right-hand side) indicates that the scale dependence of the Fermi velocity close to the Gross-Neveu critical point (or right-hand side of Fig. 2) is different from that at weak coupling, that is, Fig. 3A or the left-hand side of Fig. 2. As a representative example, the red data points and curve in Fig. 3B show our quantum Monte Carlo data and an analytical first-order perturbation theory at the Gross-Neveu critical point with the fit parameter Λ

*for the case of γ = 0.5 [other values of γ are shown in fig. S10 of (*

_{k}*16*)]. The agreement between the theory and quantum Monte Carlo data is good. Moving away from the Gross-Neveu critical line, we observe that for the Hubbard model data (γ = 0), the parameter

*C*

_{1}discussed in association with Fig. 2 also gets a linear shift with energy scale. This introduces one additional parameter in the phenomenological model, that is,

*C*

_{1}(Λ

*) =*

_{k}*C*

_{1}(0.48) − (0.92 ± 0.06)(Λ

*− 0.48), where the slope was determined also by looking at the full quantum Monte Carlo dataset. Therefore, combining this final observation for*

_{k}*C*

_{1}(Λ

*) with the fits already shown in Fig. 2 [i.e., dependence of Fermi velocity on*

_{k}*U*/

*U*(α

_{c}_{0})] and Fig. 3 (i.e., dependence on Λ

*), we can now extrapolate [see eq. S37 of (*

_{k}*16*)] our quantum Monte Carlo findings to any value of α

_{0},

*U*, and Λ

*, thereby making predictions for realistic Dirac fermions at any experimental scale.*

_{k}## Implications for experiments

Our results on the role of electron-electron interactions apply to any Dirac fermion system for which one can define the short-range component of the Coulomb interaction *U*, its long-range tail α_{0}, and the experimental probe energy scale Λ* _{k}*. However, in Fig. 4, we use our results to make predictions for graphene, because the coupling constant for graphene Dirac fermions α

_{0}=

*e*

^{2}/(κ

*ℏv*

_{0}) can be tuned by using substrates of different dielectric constants κ, with

*e*and

*ℏ*being the elementary charge and Planck’s constant

*h*divided by 2π, respectively. The on-site potential for realistic graphene is about

*U*= 3.0 (

*16*). Our results apply only at half-filling. However, in real materials, carrier density inhomogeneity makes the Dirac point inaccessible experimentally (

*31*). We use this degree of spatial inhomogeneity to define the experimental probe scale and use the root mean square (rms) of the carrier density

*n*

_{rms}in Fig. 4 (instead of Λ

*) for ease of comparison with experimental results. For the most typical situation of graphene on a substrate, typical measured Fermi velocities are 1.1 × 10*

_{k}^{6}to 1.3 × 10

^{6}m/s (

*12*); however, ab initio calculations predict 0.87 × 10

^{6}m/s [see, for example, (

*32*)]. This notable discrepancy between theory and experiment has been largely unresolved in the literature—a notable exception is (

*32*), where an ab initio density functional theory–GW calculation was used that agrees precisely with our result at our simulation scale, but their spectrum became unphysical for cm

^{–2}. In our calculations, this interaction-caused enhancement of the Fermi velocity (or, equivalently, suppression in coupling constant) can be seen directly in the left panel of Fig. 4. For example, at α

_{0}= 1 and typical experimental energy scale

*n*

_{rms}= 1 × 10

^{10}cm

^{–2}, we predict α = 0.65 ± 0.03, corresponding to Fermi velocity of (1.34 ± 0.07) × 10

^{6}m/s [see fig. S11 in (

*16*)].

We predict that for topological insulator Bi_{2}Se_{3} (α_{0} ≈ 0.05) or for graphene on metallic substrates, interactions enhance, rather than suppress, the coupling constant in the experimental window. This enhancement originates from the suppression of Fermi velocities in the Hubbard model at the Gross-Neveu fixed point. The phenomenological theory predicts that the coupling constant changes by less than 5% over the entire experimental window. Most surprisingly, for sufficiently small α_{0}, we predict that flowing Λ* _{k}* closer to the Dirac point will result in an increase in α rather than the expected decrease. These predictions challenge the conventional wisdom on the role of electron-electron interactions in 2D Dirac fermions. On the other hand, for suspended graphene, or quantum corral-like graphene, we expect the experiments to be close to the weak-coupling fixed point, and, therefore, the RPA should work well in this regime. The mechanism for this nonuniversal behavior is curious: The lattice-scale physics that depends on experimentally tunable parameters like strain and choice of substrate regularizes the Dirac theory not only in the ultraviolet, as expected, but also at the infrared, which was unexpected. All these predictions can readily be tested with current experimental capabilities.

## Supplementary Materials

www.sciencemag.org/content/361/6402/570/suppl/DC1

Materials and Methods

Figs. S1 to S11

Tables S1 and S2

This is an article distributed under the terms of the Science Journals Default License.

## References and Notes

**Acknowledgments:**We acknowledge allocation of computational resources at the CA2DM (Singapore) and the Gauss Centre for Supercomputing (SuperMUC at the Leibniz Supercomputing Center).

**Funding:**This project was supported by the Singapore National Research Foundation (NRF-NRFF2012-01), Deutsche Forschungsgemeinschaft (SFB 1170 ToCoTronics, project C01), NSERC of Canada, and Singapore Ministry of Education (MOE2014-T2-1-112 and MOE2017-T2-1-130).

**Author contributions:**S.A. and F.F.A. conceived this project. H.-K.T. and F.F.A. developed and optimized the quantum Monte Carlo codes. H.-K.T. and J.N.B.R. analyzed the data under the guidance of P.S. and F.F.A. J.N.L. developed the lattice perturbation theory and analytical results under the guidance of J.N.B.R., S.A., and I.F.H. All authors worked together to interpret and understand the results and to write the paper.

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

**Data and materials availability:**The full collection of data developed in this work is available at https://figshare.com/articles/Source_data_of_Green_s_function_and_figure_files_/5131840 (

*33*). For these calculations, we used a projective version of the auxiliary field quantum Monte Carlo approach. All the data presented in this article can be reproduced using the Algorithms for Lattice Fermions (ALF) open-source general implementation of the finite temperature auxiliary field quantum Monte Carlo available at https://alf.physik.uni-wuerzburg.de and documented in (

*34*). There, detailed documentation on input-output and error analysis can be found. The ALF implementation allows one to simulate very general lattice models, including the long-range Coulomb repulsion. All the data presented in this paper can be reproduced with ALF by using a nearest-neighbor honeycomb lattice in the low-temperature limit.