Research Article

The role of electron-electron interactions in two-dimensional Dirac fermions

See allHide authors and affiliations

Science  10 Aug 2018:
Vol. 361, Issue 6402, pp. 570-574
DOI: 10.1126/science.aao2934

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


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 = Uc 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 Embedded Image where h.c. is the Hermitian conjugate, t is the tight-binding hopping, the second-quantized operator Embedded Image (Embedded Image) creates (annihilates) an electron of spin Embedded Image at position ri, and Embedded Image gives the electron density at position ri. The interaction Embedded Image consists of a short-range part acting between electrons on the same site with different spins, Embedded Image, and a long-range part depending on the distance between the electrons Embedded Image as Embedded Image, where 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 Uc0) 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 Embedded Image 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.

Fig. 1 Phase diagram for fermions on the honeycomb lattice with competing short-range and long-range Coulomb interactions.

For any value of long-range interaction α0, there is a critical value of the short-range interaction Uc0) calculated by using quantum Monte Carlo (QMC) (data points), for which the system undergoes a quantum phase transition to the Mott insulator. In the presence of long-range interactions, a larger value of on-site interaction is required to reach the quantum phase transition. The phase diagram can be understood by solving the renormalization group (RG) flow equations (red curve), including both on-site and nearest-neighbor interactions, where the effective on-site interactions are reduced by the long-range Coulomb tail. The solid blue line is a quartic interpolation of the data points. The shaded portion shows the region inaccessible to our numerical method. Error bars indicate our numerical uncertainty. The Fermi liquid regime, the weakly interacting semimetal, and the strongly interacting Mott antiferromagnet are marked by illustrations.

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 (vFvnon)/v0 against U/Uc0), which is the short-range component of the interactions normalized by the critical value obtained from the phase diagram in Fig. 1. Here vF is the interaction-renormalized Fermi velocity obtained in quantum Monte Carlo at the simulation scale Embedded Image, where k is the momentum, vnon is the tight-binding Fermi velocity at the simulation scale, and v0 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 Gk(τ) for large τ, where the single exponential decay time, Embedded Image determines the first excitation energy of the system [defined here as the energy difference between the ground state (GS) of 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).

Fig. 2 Dirac fermion Fermi velocity renormalized by electron-electron interactions.

Projective quantum Monte Carlo results for different short-range (U) and long-range (α0) components of the Coulomb interaction. Plotted is the relative change of the Fermi velocity with respect to the noninteracting value at the Dirac point. Small U/Uc0) defines the weak-coupling regime, where Monte Carlo data for different ratios γ of the long-range and short-range components collapse as a function of α0; here electron-electron interactions enhance the Fermi velocity in agreement with the RPA (left inset). Perturbation theory (PT) results are also shown. A metal-to-Mott insulator phase transition of the Gross-Neveu universality class occurs at U = Uc0), where a suppression of Fermi velocity can be understood as the coupling between Dirac fermions and the bosonic excitations of the nascent antiferromagnetic state (the brown star is an estimate of this Fermi velocity suppression determined by using spin-wave theory). The right inset shows quantum Monte Carlo data at various values of γ for the change in Fermi velocity from the value at the Gross-Neveu point collapse onto one curve as one moves away from the phase transition. Our numerics span the full crossover between the weak-coupling fixed point and the Gross-Neveu critical point. Estimates place topological insulators close to the phase transition, whereas quantum corral-like honeycomb lattices are in the weak-coupling limit. Graphene Dirac fermions lie somewhere in between these two regimes [see table S2 in (16)].

Our numerical data show that the Fermi velocity renormalization has notably different behavior for Embedded Image (which we call the “weak-coupling regime”) and Embedded Image, 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 Bi2Se3 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 Embedded Image, where Λs is our numerical scale and Λ = 6.2 ± 0.2 is obtained from fitting the data for small α. The RPA functions F10) and F00) 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/Uc 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 vGN(γ). Defining the change in interaction strength from the critical value (ε) = 1 − U/Uc0), we find (16) that (vFvnon)/v0 = C0 + C1ε + C2ε2 + mγγ with C0 = −0.384 ± 0.002, C1 = 1.35 ± 0.05, and C2 = −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.

Fig. 3 Determining the renormalization group flow parameters from our quantum Monte Carlo data.

Our simulations also provide data for momenta larger than Λs. We exploit this scale dependence to determine scales smaller than what we can simulate. (A) Representative data close to the weak-coupling fixed point where the on-site Hubbard model (blue data) shows no observable change in Fermi velocity. With long-range interactions, the Fermi velocity increases with decreasing momenta (red data), understood either by using a continuum perturbation theory (PT) (red curve) or lattice perturbation theory with no adjustable parameter (black curve), which diverge logarithmically. (B) The Gross-Neveu (GN) critical point is very different. Here, neither the Hubbard model (blue data) nor the data including the long-range Coulomb interaction (red data) show a logarithmic divergence at small Λk. A phenomenological fit captures the on-site Hubbard model (blue dashed line). The weak increase in renormalized Fermi velocity as one goes from Λk = 2 to Λk = 1 is seen in both the black curve (lattice perturbation theory) and the quantum Monte Carlo data with finite long-range interactions. (C) Logarithmic increase at Λk ≲ 2. The red curve in (B) is obtained by fitting for this increase using a first-order perturbation theory about the Gross-Neveu critical point.

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 Λk 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.

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 Λk ~ 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 (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 Λk, 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 (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 Λk. Figure 3C shows a blowup of the region around Λk = 2, where a single adjustable parameter Λ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 Λk for the case of γ = 0.5 [other values of γ are shown in fig. S10 of (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 C1 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, C1k) = C1(0.48) − (0.92 ± 0.06)(Λk − 0.48), where the slope was determined also by looking at the full quantum Monte Carlo dataset. Therefore, combining this final observation for C1k) with the fits already shown in Fig. 2 [i.e., dependence of Fermi velocity on U/Uc0)] and Fig. 3 (i.e., dependence on Λk), we can now extrapolate [see eq. S37 of (16)] our quantum Monte Carlo findings to any value of α0, U, and Λk, thereby making predictions for realistic Dirac fermions at any experimental scale.

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 = e2/(κℏv0) 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 Embedded Image and use the root mean square (rms) of the carrier density nrms in Fig. 4 (instead of Λk) for ease of comparison with experimental results. For the most typical situation of graphene on a substrate, typical measured Fermi velocities are 1.1 × 106 to 1.3 × 106 m/s (12); however, ab initio calculations predict 0.87 × 106 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 Embedded Image 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 nrms = 1 × 1010 cm–2, we predict α = 0.65 ± 0.03, corresponding to Fermi velocity of (1.34 ± 0.07) × 106 m/s [see fig. S11 in (16)].

Fig. 4 Theoretical prediction for experimental realizations of graphene.

Solid lines in the left panel (red, close to the Gross-Neveu fixed point; blue, close to the weak-coupling fixed point) show our theory for the interaction-induced change in the coupling constant at the energy scale of our numerics for realistic graphene with different α0, as determined by the choice of substrate. Pale lines show our theory for the experimental energy scale that is numerically inaccessible but can be obtained from our numerical data by using the renormalization group (RG) analysis. Quantum Monte Carlo simulations (data points) have the same values in the left and right panels and can differentiate between these two theories. The dashed arrow shows this RG flow schematically from the numerical scale (small circle) to the experimental scale (larger circle). Right panels show the flow for α0 = 0.1 and 1.2. For the most common realization of graphene on a dielectric substrate (16), we predict weak suppression of the coupling constant that changes only slightly under renormalization (the solid and pale lines are not too different for α0 ≈ 0.6). Most surprisingly, we predict that, for small α0 (e.g., graphene on metal substrates), there is an enhancement of the coupling constant (i.e., solid red line is larger than unity) that then further increases with renormalization to the experimental scale (red line increases in the right panel).

We predict that for topological insulator Bi2Se30 ≈ 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

Materials and Methods

Figs. S1 to S11

Tables S1 and S2

References (3541)

References and Notes

  1. Materials and methods are available as supplementary materials.
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 (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 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.
View Abstract

Stay Connected to Science

Navigate This Article