Interacting Electrons in Graphene: Fermi Velocity Renormalization and Optical Response

TL;DR

This paper presents a parameter-free Hartree-Fock model for electrons in graphene that accurately predicts Fermi velocity renormalization and optical response, validated by quantum Monte Carlo simulations.

Contribution

The authors develop a novel Hartree-Fock theory for graphene that uses a topological invariant, enabling parameter-free predictions of electronic properties.

Findings

Agreement with experimental Fermi velocity data

Derived explicit optical conductivity expression

Quantum Monte Carlo results agree with mean-field predictions

Abstract

We have developed a Hartree-Fock theory for electrons on a honeycomb lattice aiming to solve a long-standing problem of the Fermi velocity renormalization in graphene. Our model employs no fitting parameters (like an unknown band cutoff) but relies on a topological invariant (crystal structure function) that makes the Hartree-Fock sublattice spinor independent of the electron-electron interaction. Agreement with the experimental data is obtained assuming static self-screening including local field effects. As an application of the model, we derive an explicit expression for the optical conductivity and discuss the renormalization of the Drude weight. The optical conductivity is also obtained via precise quantum Monte Carlo calculations which compares well to our mean-field approach.

Equations

$$\frac{v_F^{\ast }}{v_F}=1+\frac{\alpha }{4}\ln \frac{\Lambda }{k},$$

$$H_{\mathbf{k}}=-{\mathcal{E}}_{\mathbf{k}}\left[\cos ({\phi }_{\mathbf{k}}){\sigma }_x-\sin ({\phi }_{\mathbf{k}}){\sigma }_y\right],$$

$${\mathcal{E}}_{\mathbf{k}}={\mathcal{E}}_{\mathbf{k}}^0+\frac{1}{2A}\sum _{{\mathbf{k}}^{\prime }\in 1.BZ}U(\mathbf{k}-{\mathbf{k}}^{\prime })e^{i({\phi }_{{\mathbf{k}}^{\prime }}-{\phi }_{\mathbf{k}})}{F}_{{\mathbf{k}}^{\prime }},$$

$$U(\mathbf{q})$$

$$\times {\left[{\delta }_{\mathbf{G},{\mathbf{G}}^{\prime }}-v_{\mathbf{G}}(\mathbf{q}){\chi }_{\mathbf{G},{\mathbf{G}}^{\prime }}(\mathbf{q})\right]}^{-1}{v}_{{\mathbf{G}}^{\prime }}(\mathbf{q}),$$

$$\frac{{\sigma }^{\ast }}{{\sigma }_0}=1+C(\alpha )\alpha \frac{v_F}{v_F^{\ast }},$$

$$D_{ii}={\left(\frac{e}{\hslash }\right)}^2\frac{g_s}{A}\sum _{\mathbf{k}\in 1.BZ,s=\pm }s[{\partial }_{k_i}^2{\mathcal{E}}_{\mathbf{k}}]n_F(s{\mathcal{E}}_{\mathbf{k}}),$$

$$\frac{D^{\ast }}{D}=\frac{v_F^{\ast }}{v_F}.$$

$${\mathbf{a}}_1=\frac{a}{2}\left(3,\sqrt{3}\right),{\mathbf{a}}_2=\frac{a}{2}\left(3,-\sqrt{3}\right),$$

$${\mathbf{\delta }}_1=\frac{a}{2}\left(-1,\sqrt{3}\right),{\mathbf{\delta }}_2=\frac{a}{2}\left(-1,-\sqrt{3}\right),{\mathbf{\delta }}_3=a\left(1,0\right).$$

$${\mathbf{b}}_1=\frac{2\pi }{3a}\left(1,\sqrt{3}\right),{\mathbf{b}}_2=\frac{2\pi }{3a}\left(1,-\sqrt{3}\right).$$

$${\mathbf{\Psi }}_{\mathbf{k}\lambda }(\mathbf{r})=\frac{1}{\sqrt{N_c}}\sum _j{e}^{i\mathbf{k}\cdot \left({\mathbf{R}}_j+{\mathbf{\eta }}_{\nu }\right)}\zeta (\mathbf{r}-{\mathbf{R}}_j-{\mathbf{\eta }}_{\nu }){\xi }_{\sigma },$$

$$H_{\bf k}^{0}=-t|\phi_{\bf k}|\left(\begin{array}[]{cc}m_{\bf k}^{0}&e^{i\varphi_{\bf k}}\\ e^{-i\varphi_{\bf k}}&-m_{\bf k}^{0}\end{array}\right)+E_{\bf k}^{0}1_{2\times 2},$$

$${\mathcal{E}}_{\mathbf{k}}^{0,\pm }=E_{\mathbf{k}}^0\pm t\left|{\varphi }_{\mathbf{k}}\right|\sqrt{1+{m_{\mathbf{k}}^0}^2}.$$

$$\displaystyle|\psi_{{\bf k}}^{-}\rangle=\left(\begin{array}[]{c}\cos\frac{\vartheta_{{\bf k}}^{0}}{2}\\ \sin\frac{\vartheta_{{\bf k}}^{0}}{2}e^{-i\varphi_{\bf k}}\end{array}\right)\;,\;|\psi_{{\bf k}}^{+}\rangle=\left(\begin{array}[]{c}\sin\frac{\vartheta_{{\bf k}}^{0}}{2}\\ -\cos\frac{\vartheta_{{\bf k}}^{0}}{2}e^{-i\varphi_{\bf k}}\end{array}\right)\;,$$

$$V=\frac{1}{2}\sum _{i,j;\lambda ,{\lambda }^{\prime }}{c}_{i\lambda }^{†}{c}_{j{\lambda }^{\prime }}^{†}⟨i\lambda ,j{\lambda }^{\prime }|V|i\lambda ,j{\lambda }^{\prime }⟩c_{j{\lambda }^{\prime }}{c}_{i\lambda }.$$

$$c_{i\lambda }^{†}=\frac{1}{N_c}\sum _{\mathbf{k}\in 1.BZ}{e}^{-i\mathbf{k}\cdot ({\mathbf{R}}_i+{\eta }_{\nu })}{c}_{\mathbf{k}\lambda }^{†},c_{\mathbf{k}\lambda }^{†}=\frac{1}{N_c}\sum _i{e}^{i\mathbf{k}\cdot ({\mathbf{R}}_i+{\eta }_{\nu })}{c}_{i\lambda }^{†}.$$

$$V=\frac{1}{2A}\sum _{\mathbf{k},{\mathbf{k}}^{\prime },\mathbf{q}}\sum _{\lambda ,{\lambda }^{\prime }}{U}_{\mathbf{q}}^{\lambda ,{\lambda }^{\prime }}{c}_{\mathbf{k}+\mathbf{q}\lambda }^{†}{c}_{{\mathbf{k}}^{\prime }-\mathbf{q}{\lambda }^{\prime }}^{†}{c}_{{\mathbf{k}}^{\prime }{\lambda }^{\prime }}{c}_{\mathbf{k}\lambda }$$

$$H_{\mathbf{k}}^{ee}$$

$$-\frac{1}{A}\sum _{{\mathbf{k}}^{\prime }\lambda ,{\lambda }^{\prime }}{U}_{\mathbf{k}-{\mathbf{k}}^{\prime }}^{\lambda {\lambda }^{\prime }}⟨c_{{\mathbf{k}}^{\prime }{\lambda }^{\prime }}^{†}{c}_{{\mathbf{k}}^{\prime }\lambda }⟩c_{\mathbf{k}\lambda }^{†}{c}_{\mathbf{k}{\lambda }^{\prime }}.$$

$$H_{\bf k}=-t|\phi_{\bf k}^{ee}|\left(\begin{array}[]{cc}m_{\bf k}&e^{i\varphi_{\bf k}}\\ e^{-i\varphi_{\bf k}}&-m_{\bf k}\end{array}\right)+E_{\bf k}1_{2\times 2}\;.$$

$${\mathcal{E}}_{\mathbf{k}}^{\pm }=E_{\mathbf{k}}\pm t\left|{\varphi }_{\mathbf{k}}^{ee}\right|\sqrt{1+m_{\mathbf{k}}^2}$$

$$t\left|{\varphi }_{\mathbf{k}}^{ee}\right|$$

$$t\left|{\varphi }_{\mathbf{k}}^{ee}\right|m_{\mathbf{k}}$$

$$E_{\mathbf{k}}$$

$$⟨i\lambda ,j{\lambda }^{\prime }|V|i\lambda ,j{\lambda }^{\prime }⟩=\int d^2{r}_1\int d^2{r}_2|\zeta ({\mathbf{r}}_1-{\mathbf{R}}_i-{\mathbf{\eta }}_{\nu }){|}^{2}V({\mathbf{r}}_1,{\mathbf{r}}_2)|\zeta ({\mathbf{r}}_2-{\mathbf{R}}_j-{\mathbf{\eta }}_{{\nu }^{\prime }}){|}^2.$$

$$V({\mathbf{r}}_1,{\mathbf{r}}_2)=\frac{1}{A^2}\sum _{\mathbf{q}\in 1.BZ,{\mathbf{G}}_1,{\mathbf{G}}_2}{e}^{i(\mathbf{q}+{\mathbf{G}}_1){\mathbf{r}}_1}{e}^{-i(\mathbf{q}+{\mathbf{G}}_2){\mathbf{r}}_2}V(\mathbf{q}+{\mathbf{G}}_1,\mathbf{q}+{\mathbf{G}}_2),$$

$$U_{\mathbf{q}}^{\nu ,{\nu }^{\prime }}=\frac{1}{A}\sum _{\mathbf{G},{\mathbf{G}}^{\prime }}{e}^{i\mathbf{G}{\mathbf{\eta }}_{\nu }}{e}^{-i{\mathbf{G}}^{\prime }{\mathbf{\eta }}_{{\nu }^{\prime }}}{f}^{\ast }(\mathbf{q}+\mathbf{G})f(\mathbf{q}+{\mathbf{G}}^{\prime })V(\mathbf{q}+\mathbf{G},\mathbf{q}+{\mathbf{G}}^{\prime })$$

$$V(\mathbf{q}+\mathbf{G},\mathbf{q}+{\mathbf{G}}^{\prime })=A{\left[{\delta }_{\mathbf{G},{\mathbf{G}}^{\prime }}-v_{\mathbf{G}}(\mathbf{q}){\chi }_{\mathbf{G},{\mathbf{G}}^{\prime }}(\mathbf{q})\right]}^{-1}{v}_{{\mathbf{G}}^{\prime }}(\mathbf{q}),$$

$${\chi }_{\mathbf{G},{\mathbf{G}}^{\prime }}(\mathbf{q},\omega )$$

$$⟨\mathbf{r}|\mathbf{k},s⟩=\sum _{\nu =a,b}{\psi }_{\nu ,\mathbf{k},s}{\Psi }_{\mathbf{k},\lambda }(\mathbf{r})$$

$$⟨\mathbf{k},s|e^{-i(\mathbf{q}+\mathbf{G})\hat{\mathbf{r}}}|\mathbf{k}+\mathbf{q},s^{\prime }⟩=f(\mathbf{q}+\mathbf{G})\sum _{\nu =a,b}{\psi }_{\nu ,\mathbf{k},s}^{\ast }{\psi }_{\nu ,\mathbf{k}+\mathbf{q},s^{\prime }}{e}^{-i\mathbf{G}\cdot {\mathbf{\eta }}_{\nu }}.$$

Full text

Introduction. The role of Coulomb interactions in graphene is still an open and important question,Kotov et al. (2012) also in view of the regime of hydrodynamic electron liquids in which the electron-electron interaction represents the dominant scattering process.Bandurin et al. (2016); Crossno et al. (2016) The influence becomes especially crucial around the neutrality point, and it has been manifested through the measurement of the effective cyclotron mass,Elias et al. (2011) by scanning tunneling spectroscopy,Li et al. (2009); Chae et al. (2012) by direct ARPES of the Dirac cones,Siegel et al. (2011) by quantum capacitance measurements,Yu et al. (2013) and also by Landau level spectroscopy,Faugeras et al. (2015) that there is a Fermi velocity renormalization when lowering the electronic density close to half-filling.

A one-loop renormalization group (RG) and analogous Hartree-Fock (HF) analysis based on the continuous Dirac model predicts the following scaling behavior González et al. (1994); Elias et al. (2011):

[TABLE]

where $\alpha=\frac{1}{4\pi\epsilon_{0}\epsilon}\frac{e^{2}}{\hbar v_{F}}\approx 2.2/\epsilon$ is the fine-structure constant of graphene with the bare Fermi velocity $v_{F}$, $\epsilon$ characterizes the static dielectric environment, and $\Lambda$ is the momentum cut-off.

Eq. (1) has been extended by several authors,Stauber et al. (2005); Mishchenko (2007); Sheehy and Schmalian (2007); Herbut et al. (2008); Vafek and Case (2008) but a recent multi-loop expansion claims that perturbation theory may be inadequate particularly for suspended graphene.Barnes et al. (2014) Nonetheless, within a non-perturbative functional renormalization group analysis, the perturbative series can be summed up to again yield Eq. (1) with almost the same prefactor $\alpha/4$ as obtained from the HF approach.Bauer et al. (2015); Sharma and Kopietz (2016) This suggests that a self-consistent mean-field theory will contain all the necessary ingredients to address interaction effects even close to the neutrality point.

The experimental data for the velocity renormalization can be fitted to Eq. (1) by adjusting the band cut-off $\Lambda$ as well as the effective dielectric screening constant $\epsilon$.Elias et al. (2011) Nevertheless, some ambiguities inherent to the renormalization group approach can only be resolved by resorting to a realistic tight-binding Hamiltonian rather than working with an effective low-energy theory.Abedinpour et al. (2011) This especially holds for the optical conductivity which has been the subject of persistent discussion regarding the constant $C$ in the expansion $\sigma^{*}/\sigma_{0}=1+C\alpha^{*}+O({\alpha^{*}}^{2})$, with $\sigma_{0}=e^{2}/4\hbar$ the universal conductivity and $\alpha^{*}/\alpha=v_{F}/v_{F}^{*}$. Herbut et al. (2008); Mishchenko (2008); Sheehy and Schmalian (2009); Juricic et al. (2010); Sodemann and Fogler (2012); Rosenstein et al. (2013); Teber and Kotikov (2014); Link et al. (2016); Boyda et al. (2016) After almost 10 years of debate, one way to resolve this controversy could be an alternative, but well-defined numerical approach which still allows for analytical insight. We have therefore performed detailed HF calculations on the honeycomb lattice that hopefully will be able to shed some light on this issue from a different angle. We complement this with state-of-the-art quantum Monte Carlo calculations which now precisely determine the optical conductivity at energies of the order of the hopping parameter.

In contrast to earlier HF calculations preformed for a graphene quantum dot,Ozfidan et al. (2015) the Dirac model,Trushin and Schliemann (2011) multilayer graphene,Trushin and Schliemann (2012) and graphene on a latticeJung and MacDonald (2011) we now take into account self-screening and finite electronic densities which are shown to be crucial to explain the experimental data without the need for a fitting parameter. We further take advantage of the fact that the HF wave function is independent of the interaction strength even for a lattice model and relate this to a topological invariant which protects the chirality of the Dirac Fermions around the nodal points. This reduces the numerical cost and further results in HF equations that in some limits are identical to the ones obtained from RG equationsGiuliani et al. (2012) and Hubbard-Stratonovich transformationAstrakhantsev et al. (2015). The knowledge of the HF wave function further enables us to derive analytical expressions for the optical conductivity in the unscreened case with $C=1/4$, close to the value of Ref. Juricic et al., 2010. Including self-screening, we obtain $C\approx 0.05$ for suspended samples in agreement with our Monte Carlo calculations.

The Drude weight renormalization is another motivation of this work, because studies of the electromagnetic response of various classes of correlated electron materials are often based on the f-sum rule.Basov et al. (2011) Integrating the optical conductivity over the spectral range is then related to the Drude weight $D$ which is independent of the interaction in a Galilean invariant system. However, this is not the case for Dirac systems anymore, and electron-electron interactions modify the Drude weight in a non-trivial way which is larger than the Drude weight of the non-interacting system.Abedinpour et al. (2011); Levitov et al. (2013) Sum rule analysis in Dirac systems Wu et al. (2015); Post et al. (2015) have thus to be taken with care. Renormalization of the Drude weight is also of interest for plasmonics in Dirac systems as the plasmon energy scales as $\sqrt{D}$.Stauber (2014) Within our approach, we can analytically discuss the Drude weight for electronic densities close to half-filling.

The model. We will model interacting Dirac fermions in graphene within a nearest-neighbor tight-binding model using HF theory. An important insight of this work is that the crystallographic structure factor is a topological invariant which does not depend on the ground-state as long as it obeys the three-fold symmetry of the underlying honeycomb lattice, see Supplementary Information (SI).SI If we do not allow for chiralDrut and Lähde (2009) or time-reversalMarino et al. (2015) symmetry breaking, the HF wave function of the interacting system is then given by the non-interacting wave function. This is reminiscent to the absence of wave function renormalization in the RG approach.González et al. (1994)

In order to consistently include electron-electron interactions, the Coulomb potential needs to be periodic as discussed by Jung and MacDonald.Jung and MacDonald (2011) This introduces screening of the Coulomb interaction at small distances. If one further incorporates screening at large distances due to tight-binding electrons, the interaction potential is not translationally invariant anymore and additional non-diagonal local field effects need to be considered, see SI.SI To calculate the atomic orbital form factor $f(\mathbf{q})=\int d\mathbf{r}e^{-i\mathbf{q}\cdot\mathbf{r}}\left|\zeta({\mathbf{r}})\right|^{2}$ with $\zeta({\mathbf{r}})$ being the one-electron atomic wavefunction, we take the full angular dependence of the wave function into account in contrary to Refs. Jung and MacDonald, 2011 or Link et al., 2016.

Hartree-Fock theory. As mentioned above, the relative phase between the spinor components is a topological invariant and the single particle HF Hamiltonian $H_{\bf k}$ for any interaction strength can be written as

[TABLE]

where $e^{i\varphi_{\bf k}}=\phi_{\bf k}/|\phi_{\bf k}|$ with the crystallographic structure factor $\phi_{\bf k}=\sum_{i}e^{i\mathbf{k}\cdot\mathbf{\delta}_{i}}$, and the three nearest-neighbor lattice vectors $\mathbf{\delta}_{i}$, $i=1,2,3$. In the SI SI , we present the mean-field theory of a more general Hamiltonian which also includes a momentum dependent mass and energy shift for sake of generality.SI

The above Hamiltonian is characterized by the renormalized energy dispersion $\mathcal{E}_{\bf k}$ which is determined self-consistently by the following equation:

[TABLE]

where we introduced $\mathcal{E}_{\bf k}^{0}=t\left|\phi_{\bf k}\right|$ as the non-interacting dispersion relation, $t$ is the tunnel-matrix element between nearest carbon atoms and $A$ denotes the sample area. Further, we have $F_{\bf k}=n_{F}(-\mathcal{E}_{\bf k})-n_{F}(\mathcal{E}_{\bf k})$ with the Fermi distribution function $n_{F}(\epsilon)=(e^{\beta(\epsilon-\mu)}+1)^{-1}$ and the chemical potential $\mu$ at a finite temperature $\beta=1/(k_{B}T)$. The Coulomb potential conserving the lattice symmetry and including the local field effects reads

[TABLE]

where ${\bm{G}},{\bm{G}}^{\prime}$ are the reciprocal lattice vectors, $f({\bm{q}})$ is the form factor, $v_{\bm{G}}({\bm{q}})=\frac{e^{2}}{2\epsilon_{0}\epsilon}\frac{1}{|{\bm{q}}+{\bm{G}}|}$ is the Fourier-transformed screened Coulomb potential, and $\chi_{{\bm{G}},{\bm{G}}^{\prime}}({\bm{q}})$ is the static polarizability matrix with local field effects, see SI SI . Neglecting the self-consistency by replacing $\mathcal{E}_{\bf k}\rightarrow\mathcal{E}_{\bf k}^{0}$ on the right hand side of Eq. (3), we arrive at the same equation that was obtained from a Hubbard-Stratonovich transformation on the lattice.Astrakhantsev et al. (2015)

In Fig. 1, the renormalized band structure of neutral graphene is shown for $t$=3.1eV ($v_{F}=10^{6}$m/s) between the high symmetry points of the Brillouin zone for various coupling constant $\alpha$, i.e., for different dielectric environments $\epsilon$. By this, we can discuss suspended graphene ($\epsilon=1$, $\alpha=2.2$), graphene on top of silicon ($\epsilon=2.45$, $\alpha=0.9$) or hBN-encapsulated graphene ($\epsilon=4.9$, $\alpha=0.45$). The solid lines refer to self-screened interaction which are compared to the dispersion due to bare interaction for $\alpha=2.2$ (dashed line). The inset shows the region close to the Dirac point where only slight deviations from the linear behavior can be seen.

Fermi velocity renormalization. At half-filling $\mu=0$, $T=0$ and no self-screening, Eq. (3) is an explicit equation. This yields the analytical expression of Eq. (1) when assuming $U({\bm{q}})=v_{{\bm{G}}=0}({\bm{q}})$ and converting the summation over the Brillouin zone by the Dirac cone approximation. Solving Eq. (3) numerically, we obtain a fit for the cut-off parameter with $\Lambda\approx 1.75\mathrm{\mathring{A}}^{-1}$. This is in contrast to Ref. Jung and MacDonald, 2011, where $\Lambda\approx 20\mathrm{\mathring{A}}^{-1}$ was obtained, but agrees well with the usual argument of fixing $\Lambda$ by conserving the total number of states in the Brillouin zone when compared to the tight-binding model, yielding $\Lambda\approx 1.58\mathrm{\mathring{A}}^{-1}$. The precise value of $\Lambda$ depends only weakly on the non-universal short-ranged Coulomb interaction, see SI.SI

For unscreened interaction, the correction to the Fermi velocity $v_{F}^{*}/v_{F}-1$ is proportional to $\alpha$ and the result is scale-independent of the hopping parameter $t$. At zero doping, self-screening can be incorporated by $\alpha\rightarrow\alpha/\epsilon^{RPA}$ because $\epsilon^{RPA}=1+\frac{\pi}{2}\alpha$ is momentum independent within the Dirac-cone approximation Kan04 . This yields good agreement with the experimental data of Ref. Elias et al., 2011 for small densities $n\lesssim 20\times 10^{10}$cm*-2**without* the need for any fitting parameter, see magenta curve of Fig. 2.

For densities $n\gtrsim 20-60\times 10^{10}$cm*-2*, there is a decrease of the Fermi velocity which cannot be accounted for by the results for neutral graphene. Since it might be due to screening at finite densities, we incorporated the momentum dependent polarization function as outlined above. Even though the renormalized Fermi velocity now depends on $\alpha$ as well as on $t$ in a non-trivial way, in the asymptotic limit it becomes independent of $t$ for $\mu=0$.

On the left hand side of Fig. 2, we show the solution of Eq. (3) using the bare (black squares) and the self-consistent (blue stars) polarization function. The bare solution agrees well with the experimental data for suspended graphene up to $n\lesssim 40\times 10^{10}$cm*-2*, but at higher densities the experimental data drops whereas the theoretical value remains approximately constant (on a logarithmic scale). This has to be contrasted with the experimental data for hBN-encapsulated graphene,Yu et al. (2013) where good agreement is obtained over the whole density range up to $n\sim 5\times 10^{12}$cm*-2*, see Fig. 2 (right). For the particular choice of the bare hopping amplitude in the two cases, i.e., $t=3.1$eV and $t=2.6$eV, respectively, we were guided by the original Refs. Elias et al. (2011) and Yu et al. (2013) where similar values were used, see also Ref. Rei02 .

Optical response. Let us now turn to the interaction effects on the optical response, first discussed for Dirac electrons in Ref. Stauber and Peres, 2008 in the case of electron-phonon coupling. To do so, we will couple the gauge field via the Peierls substitution by replacing ${\bf k}\rightarrow{\bf k}+\frac{e}{\hbar}\mathbf{A}$ in the mean-field Hamiltonian $H_{\bf k}$. This procedure provides the correct vertex correction such that the optical f-sum rule is satisfied.

Since the HF wave function is known, the optical conductivity can be deduced from the non-interacting tight-binding modelPeres and Stauber (2008) by replacing the bare dispersion by the renormalized one. For momenta close to the Dirac point, the dispersion is isotropic with $v^{*}_{F}/v_{F}=1+C(\alpha)\alpha\ln\Lambda/k$. We then obtain for small frequencies $\omega\ll t/\hbar$ the following result:SI

[TABLE]

where $\sigma_{0}=\frac{e^{2}}{4\hbar}$ denotes the universal conductivity. For unscreened interaction, we obtain explicitely $C=1/4$ (see SISI ) as mentioned in the introduction. This compares well with $C\approx 0.26$ obtained in Refs. Juricic et al., 2010; Rosenstein et al., 2013.

For the self-screened interaction, $C\rightarrow C(\alpha)$ becomes a function with $C(\alpha\to 0)\to 0.25$ (unscreened limit) and $C(\alpha=2.5)\approx 0.05$ for suspended graphene with $t=2.7$eV. It is interesting to note that the universal factor of the scaling law $C(\alpha)$ is independent of all considered hopping matrix elements from $t=2.6$eV to $t=3.1$eV, and it compares well if self-screening is incorporated via RPA within the Dirac cone approximation, i.e., $C(\alpha)=[4(1+\frac{\pi}{2}\alpha)]^{-1}$. This is shown on the left of Fig. 3; assuming a small, but finite electronic density would further decrease the constant $C(\alpha)$.

On the right hand side, we plot the same for the conductivity of Eq. (5) with $v_{F}^{*}/v_{F}=1$. This is compared to the conductivity at $\hbar\omega=0.7t$ as obtained from HF and quantum Monte Carlo calculations, see SI SI . We obtain good agreement between the two approaches for suspended graphene which justifies our mean-field approach.

Drude weight. In a Galilean invariant system, the Drude weight $D$ is independent of the interaction. However, this is not the case for Dirac systems and electron-electron interactions modify the Drude weight in a non-trivial way which is larger than the Drude weight of the non-interacting system. Abedinpour et al. (2011); Levitov et al. (2013)

Making use of the fact that the HF wave function is given by the non-interacting wave function, one can obtain an analytical expression for the Drude weight from the optical $f$-sum rule (see SISI ):

[TABLE]

with $g_{s}=2$ the spin degeneracy, and $i=x,y$. This expression for the Drude weight in the presence of Coulomb interactions generalizes the known result for non-interacting systems: $D=e^{2}n/m$ for Schrödinger particles with the density $n$ and mass $m$ and $D=g_{s}(e/h)^{2}\pi\mu$ for for Dirac particles in graphene.Stauber (2014) As it might have been expected, within our mean-field theory the interacting electrons behave as independent quasi-particles with renormalized dispersion $\lambda\mathcal{E}_{{\bf k}}$ and we obtainSI

[TABLE]

A similar relation was obtained in Ref. Abedinpour et al., 2011 in the case of unscreened interactions. Changes due to trigonal warping and finite temperature can also be discussed by our approach.

Summary. We presented a realistic tight-binding approach to the electron band structure of graphene renormalized by the Coulomb interaction. We identified a topological invariant which leaves the HF wave function unchanged even in the presence of the Coulomb interaction and found analytical expressions for the optical conductivity as well as for the Drude weight. By this, we were able to link our findings to the measured optical conductivity which shows only little renormalization due to self-screened interactions. Precise Monte Carlo calculations yield good agreement for suspended samples and support our mean-field approach. We also show that the Fermi velocity and Drude weight renormalization are the same according to the expectations of a mean-field theory.

Our results compare well with experiments for suspended as well as for hBN-encapsulated graphene without invoking fitting parameters. But in the case of suspended graphene, we were not able to account for the velocity renormalization in the case of larger densities $n\gtrsim 40\times 10^{10}$cm*-2*. This effect cannot be explained by our tight-binding model and we expect the influence of ripples and corrugations, partially due to the applied gate, to be responsible for this effective screening of the long-ranged Coulomb interaction. This would also imply the absence of interaction renormalization of (high-density) plasmons in suspended graphene.

Acknowledgements. We acknowledge interesting discussions with Igor Herbut, Jeil Jung, Peter Kopietz, Vieri Mastropietro, and Anand Sharma. We also thank Kostya Novoselov for providing us with the experimental data. This work has been supported by Spain’s MINECO under Grant No. FIS2014-57432-P, by the Comunidad de Madrid under Grant No. S2013/MIT-3007 MAD2D-CM, and by Germany’s Deutsche Forschungsgemeinschaft (DFG) via SFB 689. The work of M. Ulybyshev was supported by DFG grant BU 2626/2-1. D. Boyda acknowledges the support by RFBR Grant No. 16-32-00362-mol-a.

Supplementary Information

I The non-interacting tight-binding model

We will model interacting Dirac fermions in graphene within a nearest-neighbor tight-binding model. The honeycomb lattice of graphene is a non-Bravais bipartite lattice, which means that it consists of two interpenetrating sublattices, each of them forming a triangular Bravais lattice. The lattice vectors are chosen as follows:

[TABLE]

where a $\sim 1.42\,\mathrm{\mathring{A}}$ is the C-C covalent bond length. The three nearest-neighbor vectors in real space are given by

[TABLE]

The reciprocal lattice vectors ${\bm{b}}_{1}$ and ${\bm{b}}_{2}$ defined by the condition ${\bm{a}}_{i}\cdot{\bm{b}}_{j}=2\pi\delta_{ij}$ are then given by

[TABLE]

Within the nearest-neighbor tight-binding model, the general Bloch basis state is given by

[TABLE]

where $N_{c}$ is the number of unit cells, $\xi_{\sigma}$ denotes the spin part of the wave function, $\zeta$ the one-electron atomic wavefunction ($p_{z}$ orbital of carbon) at $\bm{R}_{j}+\bm{\eta}_{\nu}$, $\bm{\eta}_{\nu}$ is the position of sublattice $\nu$ in the crystallographic basis, and $\lambda=(\nu,\sigma)$ shall include sublattice and spin degrees of freedom. In the following, we will mostly choose $\bm{\eta}_{a}=(0,0)$, $\bm{\eta}_{b}=(a,0)$.

The free tight-binding Hamiltonian can then be written as

[TABLE]

where $t$ = 3.1eV, $\phi_{\bf k}=\sum_{i=1,2,3}e^{i\bm{k}\cdot\bm{\delta}_{i}}$ the structure factor, and $e^{i\varphi_{\bf k}}=\phi_{\bf k}/|\phi_{\bf k}|$. We also included a mass term and a constant energy term for sake of generality. The eigenenergies read

[TABLE]

The eigenvectors are given by

[TABLE]

with $\cos\vartheta_{\bf k}^{0}=m_{\bf k}^{0}/\sqrt{1+{m_{\bf k}^{0}}^{2}}$ and $\sin\vartheta_{\bf k}^{0}=1/\sqrt{1+{m_{\bf k}^{0}}^{2}}$. We note that by defining the Bloch state as usual, i.e., $\bm{\Psi}_{{{\bf k}}\lambda}(\bm{r}+\bm{R}_{i})=e^{i{\bf k}\cdot\bm{R}_{i}}\bm{\Psi}_{{{\bf k}}\lambda}(\bm{r})$, we arrive at the eigenvectors of the Hamiltonian with a relative phase $e^{-i\varphi_{\bf k}}$ between the first and second spinor component.

II Mean-Field theory

Let us introduce the Hamiltonian for the electron-electron interaction as

[TABLE]

We define the Fourier transformation as follows:

[TABLE]

This is consistent with Eq. (11), i.e., ${\Psi}_{{{\bf k}}\lambda}(\bm{r})=\langle{\bm{r}}|c_{{\bm{k}}\lambda}^{\dagger}|0\rangle$. With the Coulomb propagator $U_{\bm{q}}^{\lambda,\lambda^{\prime}}$ defined in the next section, we then have

[TABLE]

where $A=A_{c}N_{c}$ is the system area with $A_{c}=3\sqrt{3}a^{2}/2$.

The interaction term within the mean-field approximation reads

[TABLE]

The two terms are the Hartree and Fock (exchange) term, respectively. In the following, we will only consider the exchange interaction of the HF Hamiltonian as the Hartree term is neutralized by a positive background.

We numerically solved the total Hamiltonian $H=\sum_{\bf k}H_{\bf k}$, $H_{\bf k}=H_{\bf k}^{0}+H_{\bf k}^{ee}$, and found that the in-plane (azimuthal) angle of the pseudospin phase $\varphi$ is not changed by the interaction strength. This can also be shown analytically, see below. The relative phase between the spinor components can thus be seen as a topological invariant and with $\phi_{\bm{k}}^{ee}=\langle c_{{\bm{k}},b,\sigma}^{\dagger}c_{{\bm{k}},a,\sigma}\rangle/2\propto\phi_{\bm{k}}$, the single particle $H_{\bf k}$ for any interaction strength can be written as

[TABLE]

This parametrization will allow for a more efficient numerical analysis, crucial in order to address the physics close to the neutrality point. The eigenenergies of the interacting system now read

[TABLE]

and the eigenvectors are again given by Eq. (18) after replacing $\vartheta_{\bf k}^{0}\rightarrow\vartheta_{\bf k}$ where $\cos\vartheta_{\bf k}=m_{\bf k}/\sqrt{1+m_{\bf k}^{2}}$ and $\sin\vartheta_{\bf k}=1/\sqrt{1+m_{\bf k}^{2}}$.

In the above equations, we have introduced the renormalized energy dispersion $|\phi_{\bf k}^{ee}|$ if $m_{\bf k}=0$ and $E_{\bf k}=0$, the total dimensionless, ${\bf k}$-dependent mass $m_{\bf k}$, and the total ${\bf k}$-dependent energy shift $E_{\bf k}$. These quantities are determined self-consistently by the following relations:

[TABLE]

Above, we defined $F_{\bf k}^{\pm}=f({\mathcal{E}}_{{\bf k}}^{-})\pm f({\mathcal{E}}_{{\bf k}}^{+})$ with $f(\epsilon)=(e^{\beta(\epsilon-\mu)}+1)^{-1}$ the Fermi distribution function and $\mu$ the chemical potential at a finite temperature $\beta=1/(k_{B}T)$.

During the iteration process to obtain self-consistency, the chemical potential relative to the neutrality point will be kept constant. For $\mu=0$, Eq. (27) becomes trivial because $\sum_{{\bf k}^{\prime}}U^{11}_{{\bf k}-{\bf k}^{\prime}}=0$. If we also neglect the self-consistency, we arrive at the same set of equation if one-loop corrections to the electron propagator and to the interaction potential are considered.Astrakhantsev et al. (2015)

III Screening and local field effects

We will now specify the Coulomb propagator $U^{\nu,\nu^{\prime}}_{\bm{q}}$ defined in Eq. (21). If one incorporates screening due to tight-binding electrons, the interaction potential is not translationally invariant anymore and we have instead ($\lambda=(\nu,\sigma)$)

[TABLE]

With the condition $V(\bm{r}_{1},\bm{r}_{2})=V(\bm{r}_{1}+{\bm{R}}_{i},\bm{r}_{2}+{\bm{R}}_{i})$, we obtain the following representation:Giuliani and Vignale (2005)

[TABLE]

where the sum runs over all reciprocal lattice vectors ${\bm{G}}=n\bm{b}_{1}+m\bm{b}_{2}$ with integers $n,m\in\bm{Z}$. The interaction potential $H^{ee}$ of Eq. (21) is then defined with the Coulomb propagator given by

[TABLE]

with $f({\bm{q}})=\int d^{2}r|\zeta(\bm{r})|^{2}e^{-i{\bm{q}}\cdot\bm{r}}$ and where we used $c_{{\bm{k}}+{\bm{G}},\lambda}=e^{i{\bm{G}}\cdot{\bm{\eta}}_{\nu}}c_{{\bm{k}},\lambda}$. We note that the definition is independent of the choice of the crystallographic basis ${\bm{\eta}}_{\nu}$ and that $U^{\nu,\nu^{\prime}}_{\bm{q}}=U^{\nu^{\prime},\nu}_{-{\bm{q}}}$. It can also be shown that $U^{1,1}_{\bm{q}}=U^{2,2}_{{\bm{q}}}$ expected from the equivalence of the two sub-lattices.

The screened potential shall be calculated within the RPA-approximation:

[TABLE]

with $v_{\bm{G}}({\bm{q}})=v({\bm{q}}+{\bm{G}})=\frac{e^{2}}{2\epsilon_{0}\epsilon|{\bm{q}}+{\bm{G}}|}$. The (dynamical) response function is given by

[TABLE]

The eigenfunctions are given by

[TABLE]

with the general Bloch eigenstate $\Psi_{{\bm{k}},\lambda}(\bm{r})$ given in Eq. (11) and $\psi_{a,{\bm{k}},-}=\cos\frac{\vartheta_{{\bf k}}}{2}$, $\psi_{b,{\bm{k}},-}=\sin\frac{\vartheta_{{\bf k}}}{2}e^{-i\varphi_{\bf k}}$, $\psi_{a,{\bm{k}},+}=\sin\frac{\vartheta_{{\bf k}}}{2}$, and $\psi_{b,{\bm{k}},+}=-\cos\frac{\vartheta_{{\bf k}}}{2}e^{-i\varphi_{\bf k}}$. The matrix overlap function is given by

[TABLE]

Finally, we note that due to time-reversal symmetry we have $\chi_{{\bm{G}},{\bm{G}}^{\prime}}({\bm{q}},\omega)=\chi_{{\bm{G}}^{\prime},{\bm{G}}}(-{\bm{q}},\omega)$.

We note that the condition

[TABLE]

is fulfilled in all cases which is needed to fulfill the self-consistency equations as we will show in the next section.

III.1 Dirac cone approximation

For small in-plane momentum $q$, we can use the Dirac cone approximation for the evaluation of the sum over the Brillouin zone. This can be done analytically and we get for the doped case:

[TABLE]

with $\gamma_{1}=(e^{-i{\bm{G}}\cdot\eta}+e^{i{\bm{G}}^{\prime}\cdot\eta})/2$, $\gamma_{2}=(e^{-i{\bm{G}}\cdot\eta}+e^{i{\bm{G}}^{\prime}\cdot\eta})/2-(1-e^{-i{\bm{G}}\cdot\eta})(1-e^{i{\bm{G}}^{\prime}\cdot\eta})/4$, $\gamma_{3}=(1+e^{-i{\bm{G}}\cdot\eta})(1+e^{i{\bm{G}}^{\prime}\cdot\eta})/4$, $\gamma_{4}=(1-e^{-i{\bm{G}}\cdot\eta})(1-e^{i{\bm{G}}^{\prime}\cdot\eta})/8$. For the undoped case, this reduces to

[TABLE]

The band-cut off $\Lambda$ shall not be confused with the one used in the scaling law of the Fermi velocity and can be obtained from the exact numerical solution. We get $\Lambda=2.29\,\mathrm{\mathring{A}}^{-1}$.

IV Topological invariant leads to phase locking between the spinor components

Here, we show that the phase between the components of the spinor is invariant with respect to the Coulomb interaction. This property is due to the boundary condition $c_{{\bm{k}}+{\bm{G}},\lambda}=e^{-i{\bm{G}}\cdot{\bm{\eta}}_{\nu}}c_{{\bm{k}},\lambda}$ and the underlying three-fold and at the Dirac point rotational symmetry, i.e., we neglect possible symmetry breaking due to interaction.

Let us define $\phi_{\bm{k}}^{ee}=\langle c_{{\bm{k}},b,\sigma}^{\dagger}c_{{\bm{k}},a,\sigma}\rangle/2$. The boundary condition implies

[TABLE]

where we choose $\bm{\eta}_{a}=(0,0)$, $\bm{\eta}_{b}=(a,0)$.

Due to the three-fold symmetry, we have $\phi_{\bm{k}}^{ee}=\phi_{O{\bm{k}}}^{ee}=\phi_{O^{2}{\bm{k}}}^{ee}$ where $O$ is a rotation about the angel $\varphi=2\pi/3$. Thus we can formally set $f_{1}({\bm{k}})=\phi_{\bm{k}}^{ee}/3$, $f_{2}({\bm{k}})=\phi_{O{\bm{k}}}^{ee}/3$, $f_{3}=\phi_{O^{2}{\bm{k}}}^{ee}/3$ and write $\phi_{\bm{k}}^{ee}=\sum_{i=1..3}f_{i}({\bm{k}})$. At the Dirac point, the system is further rotationally invariant and we can set $f_{i}({\bm{k}})=f({\bm{v}}_{i}\cdot{\bm{k}})$, with an arbitrary function $f$ and ${\bm{v}}_{i}$ being the three vectors related by a rotation of $\varphi=2\pi/3$. In order to obtain the boundary condition Eq. (38), the arbitrary function must be related to the exponential and we thus have $f(x)=Ae^{ix}$ and ${\bm{v}}_{i}={\bm{\delta}}_{i}$ the nearest-neighbor vectors with an arbitrary constant $A$.

We thus have $\phi_{\bm{k}}^{ee}\propto\sum_{i}e^{i{\bm{\delta}}_{i}\cdot{\bm{k}}}=\phi_{\bm{k}}={}_{0}\langle c_{{\bm{k}},b,\sigma}^{\dagger}c_{{\bm{k}},a,\sigma}\rangle_{0}/2$. The crystallographic structure factor in the case of the non-interacting system can thus be obtained from general symmetry arguments with the same boundary conditions as in Eq. (38) and it is independent of the interaction or other perturbations which preserve the basic symmetries.

Finally, let us note that Eq. (4.21) of Ref. Giuliani et al. (2012) comes to the same conclusion based on the Feynman graphs with two external fermionic lines with the same wave vector ${\bf k}$ and a more formal proof of this relation can be found in App. C of the same reference.

V Derivation of the form factor

Including lattice scale effects, the explicit form of $f(q)$ can be calculated considering $\zeta(r)$ to be a hydrogenic 2$p_{z}$ orbital type, $\zeta(r,\vartheta)=\frac{1}{4\sqrt{2\pi}}(\frac{Z}{a_{0}})^{3/2}\frac{Zr}{a_{0}}e^{\frac{-Zr}{2a_{0}}}\cos(\vartheta)$ with $Z$ the effective (screened) atomic charge. By taking the Fourier transformation of the charge distribution, the form factor is defined by

[TABLE]

We choose q to be in the $xy$-plane with q = q(1,0,0) and define $\tilde{a}_{0}=\frac{a_{0}}{Z}$. This yields

[TABLE]

We thus obtain the final result:

[TABLE]

With $\langle r^{2}\rangle=30\tilde{a}_{0}^{2}$, we can reproduce the covalent bond radius of carbon by choosing $\tilde{a}_{0}=a/2/\sqrt{30}$. A larger effective radius could also consider screening effects from the sp2-orbitals and we will choose $\tilde{a}_{0}=3a/2/\sqrt{30}$.Jung and MacDonald (2011) The short-distance cut-off has only little influence on the scaling behaviour of the Fermi velocity, i.e., the cut-off parameter only slightly increases with increasing short-distance screening corresponding to the absolute band width with $\Lambda=1.75\,\mathrm{\mathring{A}}^{-1}$ and $\Lambda=1.82\,\mathrm{\mathring{A}}^{-1}$ for the two parameters mentioned above.

VI Optical response

Let us now discuss the interaction effects on the optical response. For this, we will couple the gauge field employing the Peierls substitution by replacing ${\bf k}\rightarrow{\bf k}+\frac{e}{\hbar}\bm{A}$ in the mean-field Hamiltonian $H_{\bf k}$. This procedure provides the correct vertex correction such that the optical f-sum rule is satisfied.

Assuming a linearly polarized light along the $i$-direction with $i=x,y$, the current operator reads

[TABLE]

where $j_{i}^{P}$ and $j_{i}^{D}$ are the paramagnetic and diamagnetic contribution, respectively.

The real part of the conductivity $\sigma_{ii}(\omega)$ represents the optical absorption of the light. Using the Kubo formula within the linear response formalism, $\Re{\sigma}_{ii}$ can be split into two terms containing a regular part and a delta singularity

[TABLE]

where $D_{ii}$ is the Drude weight corresponding to the charge stiffness. The regular part of the conductivity reads for $\omega>0$

[TABLE]

where $P_{\bf k}^{i}=\langle\psi_{{\bf k}}^{+}|\hbar v_{\bf k}^{i}|\psi_{{\bf k}}^{-}\rangle$ is the interband momentum matrix element with the velocity operator along $i$-direction defined as $\hbar v_{\bf k}^{i}=\partial_{k_{i}}H_{\bf k}$, see Eq. (42). For the mean-field Hamiltonian $H_{\bf k}$, Eq. (23), we finally obtain

[TABLE]

If $m_{\bf k}=0$,

[TABLE]

where $|P_{\bf k}^{i,0}|^{2}=\frac{t^{2}a^{2}}{16}g_{\bf k}$ corresponds to the value in non-interacting case and

[TABLE]

with $\tilde{\phi}_{\bf k}=e^{-i\bm{k}\cdot\bm{\delta}_{3}}\phi_{\bf k}$.Peres and Stauber (2008)

VI.1 Optical conductivity

There has been considerable work on the renormalization of the conductivity due to electron-electron interaction. To leading order in $\alpha^{*}=\alpha\frac{v_{F}}{v_{F}^{*}}$, one gets via a renormalization group analysis

[TABLE]

with $C_{1}=25/12-\pi/2\approx 0.51$ (hard cut-off),Herbut et al. (2008) $C_{2}=19/12-\pi/2\approx 0.01$ (soft cut-off),Mishchenko (2008); Sheehy and Schmalian (2009) and $C_{3}=11/6-\pi/2\approx 0.26$ (dimensional regularization).Juricic et al. (2010) Chiral anomalies have been claimed to be responsible for these discrepancies and a perturbative analysis based on the tight-binding model yields $C=C_{3}$.Rosenstein et al. (2013)

Since for particle-hole symmetric Hamiltonian there is no sign-problem, the optical conductivity of graphene can also be studied using quantum Monte Carlo calculations and values $C\approx 0.05$ were obtained for different interaction strengths.Boyda et al. (2016) Large error bares were attached to these calculations, but we have implemented an improved scheme for the analytical continuation, obtaining significantly more precise results (see below).

From our Hartree-Fock calculations, we obtain the analytical result $C=0.25$ for unscreened interaction, but numerical results suggest lower values for self-screened interaction with $C\approx 0.05$ for suspended graphene ($\alpha=2.5$ for $t=2.7$eV).

For $m_{\bm{k}}=0$ and momenta close to the Dirac point, we have $|\phi_{\bm{k}}|=\frac{3}{2}ak$ and an isotropic dispersion. The regular part of the real optical conductivity then reads:

[TABLE]

As shown by Sharma and Kopietz,Sharma and Kopietz (2016) close to the neutrality point we have the following functional behavior valid to all orders in $\alpha$:

[TABLE]

From Eq. (25), we get $A=1$ and integrating the above scaling law yields $\mathcal{E}_{\bm{k}}^{+}=\hbar v_{F}k(1+B\ln(\Lambda/k)+B)$. For small $\omega$, we then get the following result:

[TABLE]

with $\sigma_{0}=\frac{g_{s}g_{v}e^{2}}{16\hbar}$ the universal conductivity. For unscreened interaction, we have the analytical result $B=\alpha/4$ and thus $C=0.25$. For screened interaction, a good approximation is given by $B=\alpha/4\epsilon^{RPA}$ and thus $C(\alpha)=0.25/(1+\frac{\pi}{2}\alpha)$. These results also holds in the Dirac cone approximation.

VI.2 Drude weight

The optical $f$-sum rule reads

[TABLE]

where $j_{D}^{i}\equiv-\frac{e^{2}}{\hbar^{2}}\sum_{{\bf k}}\partial_{k_{i}}^{2}H_{\bf k}$ is the diamagnetic current operator. From the above $f$-sum rule we obtain the expression for the Drude weight:

[TABLE]

Because of the knowledge of the HF-wave function, we can derive an analytical formula for the Drude weight in the presence of electron-electron interaction. With the spin-degeneracy $g_{s}=2$, we obtain the general result, valid also in the presence of an arbitrary mass term,

[TABLE]

Note that the expectation value of the diamagnetic current operator, $\langle j_{D}^{i}\rangle$, is not equal to the Drude weight since the Hellmann-Feynman theorem does not apply for the second derivative. For particles with mass $m$ and parabolic dispersion, the well-known Drude weight $D=e^{2}n/m$ with $n$ the particle density is obtained.Stauber (2014) The above formula generalizes this result to two bands and we expect it to hold also for general multi-band systems (index $s=1...n$).

Using partial integration, one gets at $T=0$ in the thermodynamic limit $A_{s}\to\infty$

[TABLE]

where the last equation holds for isotropic Fermi surfaces.We thus obtain a direct link between the Drude weight and the Fermi energy $\mu$. Comparing systems with the same density, i.e., with the same Fermi wave number $k_{F}$, we obtain the same renormalization of the Drude weight as for the Fermi velocity:

[TABLE]

The effect of temperature, large doping and disorder might lead to important changes which can be discussed in the presented formalism.

VII Numerical implementation and results

The Hartree-Fock equations with unscreened interaction are numerically solved using a grid with up to $N_{c}=15000^{2}$ lattice sites. This allows us to accurately discuss the scaling behavior for the undoped system as well as the Fermi-velocity renormalization for small finite electronic densities down to $n=10^{10}$cm*-2*.

In order to match the experimental results, it is crucial to incorporate self-screening effects. For this, the static polarization function $\chi_{{\bm{G}},{\bm{G}}^{\prime}}({\bm{q}})$ has been calculated for momenta on a small grid $N_{c}=300^{2}$ and later we used bilinear interpolation to obtain the polarization on a larger grid $N_{c}=3000^{2}$. The momentum summation for the polarization function has also been performed on a larger grid $N_{c}=600^{2}$. Finally, we approximated $\chi$ for small momenta by the analytical formula of the Dirac cone approximation, see Eq. (36). In both cases, i.e., for the unscreened as well as for the screened interaction, the sum and matrix over the reciprocal lattice vectors ${\bm{G}}=n\bm{b}_{1}+m\bm{b}_{2}$ is truncated by $|n|,|m|\leq n_{max}$ with $n_{max}=4,6$.

For the self-consistent solution, we approximated the static polarization by the bare polarization divided by the renormalized Fermi velocity. This is justified following the work by Sodemann and FoglerSodemann and Fogler (2012) who found that $\epsilon(q)=1+\frac{\pi}{2}\alpha^{*}+O({\alpha^{*}}^{2})$ where $\alpha^{*}=\frac{1}{4\pi\epsilon_{0}\epsilon}\frac{e^{2}}{\hbar v_{F}^{*}}$ is the fine-structure constant with respect to the renormalized velocity. The self-consistent dispersion is scaled by a factor 1.4 compared to the dispersion coming from the bare polarization and is usually larger than the experimental data.

VII.1 Electronic dispersion

In Fig. 4, we show the universal scaling behavior $(\frac{v_{F}^{*}}{v_{F}}-1)/\alpha$ for different grid sizes $N_{c}=3000^{2},15000^{2}$. A fit to

[TABLE]

with $B=\alpha/4$ yields $\Lambda=1.75\,\mathrm{\mathring{A}}^{-1}$. Screening can be incorporated by $\alpha\rightarrow\alpha/\epsilon^{RPA}$ which agrees well with the full calculation.

For finite densities, we observe a kink in $v_{F}^{*}$ at $k_{F}$ which is smeared out at finite temperatures. This kink is clearly an artifact of the unscreened interaction and screening cannot be incorporated afterwards. It is thus crucial to calculate the Fermi velocity renormalization including the self-screened interaction. The results are shown in Fig. 5 for the two orthogonal directions of the Bloch momentum. As expected, finite doping acts as a cutoff of the scaling law, but sizable effects are still observed if we trace the Fermi velocity at the Fermi level $k_{F}$. We also obtain a larger prefactor of the logarithm (black dashed line) compared to the undoped case (magenta dotted-dashed line) along the $KM$-direction.

Finally, in Fig. 6, we present the numerical curves for the renormalized Fermi velocity at the neutrality point for various coupling strengths $\alpha$. We used $t=3.1$eV, but the asymptotic results do not depend on the hopping parameter. The fits to Eq. (57) are shown as blue lines. They are slightly different for the two orthogonal directions, and we will consider the one along the $K\Gamma$-direction since it is more close related to Dirac cone physics.Stauber et al. (2010); Stauber (2010)

VII.2 Optical conductivity

In Fig. 7, we show the optical conductivity for suspended graphene with $t=2.7$eV in the vicinity of the Dirac point, based on a system with $N_{c}=3000^{2}$ lattice sites. Convergence is reached for energies up to $\hbar\omega\approx 1.2t$ for $N_{\Delta K}=513^{2}$ grid points of the Brillouin zone around one Dirac cone (due to time-reversal symmetry it is sufficient to only consider one Dirac cone). We also show the optical conductivity for graphene on a substrate with $\epsilon=2$ or $\alpha=1.25$ (blue) and $\epsilon=5$ or $\alpha=0.5$ (green) for $N_{\Delta K}=413^{2}$.

In order to compare the mean-field results with the quantum Monte Carlo calculations, we choose an energy of $\hbar\omega=0.7t$ where the plateau of the conductivity has developed (see below). We obtain for the interacting conductivity $\sigma^{*}=1.117$ for $\epsilon=1$, $\sigma^{*}=1.109$ for $\epsilon=2$, and $\sigma^{*}=1.095$ for $\epsilon=5$. These values show excellent agreement with the improved quantum Monte Carlo calculations, see below.

VIII Precise Quantum Monte Carlo study of graphene conductivity

We perform Quantum Monte Carlo (QMC) calculations for the following interacting tight-binding model for electrons at $\pi$-orbitals:

[TABLE]

Here $\sum_{\langle x,y\rangle}$ means the summation over all pairs of nearest-neighbor sites. $\hat{a}^{{\dagger}}_{x,\sigma}$, $\hat{a}_{x,\sigma}$ are the creation/annihilation operators for electrons with spin-up and spin-down, respectively. $t$=2.7 eV is the hopping amplitude. $\hat{q}_{x}=\sum_{\sigma}\hat{a}^{{\dagger}}_{x,\sigma}\hat{a}_{x,\sigma}-1$ is the charge operator at site $x$. $V_{xy}$ is the electron-electron interaction potential. We use the potentials calculated with the constrained RPA method Wehling et al. (2011) for suspended graphene (see Ulybyshev et al. (2013) for details). These potentials correspond to the case $\varepsilon=1.0$ in our calculations. For $\varepsilon\neq 1$ we uniformly rescale potentials at all distances by a factor of $1/\varepsilon$. Since we can treat only finite sample, we impose periodic spatial boundary conditions as in Refs. Buividovich and Polikarpov (2012); Ulybyshev et al. (2013); Smith and von Smekal (2014).

The calculation of graphene conductivity is performed with the same numerical procedure as was already described in Boyda et al. (2016). We calculate Euclidean current-current correlator $G(\tau)$ using Hybrid Monte Carlo algorithm, which is especially advantageous for systems of large spatial extent in the presence of long-range interaction. Since this particular Hamiltonian 58 doesn’t lead to sign problem in QMC calculations, we obtain numerically exact results without any further physical assumptions. Optical conductivity $\sigma(\omega)$ is extracted from the Green-Kubo relation:

[TABLE]

This integral equation is unstable thus we need a special numerical algorithm with some kind of regularization to find the stable solution. It’s especially important in our case when the input data for correlator $G(\tau)$ is defined up to some statistical errors because it was obtained in statistical Monte Carlo procedure. In general we rely on the Backus-Gilbert (BG) method proposed in Brandt et al. (2015), which we already employed in Boyda et al. (2016).

In order to improve accuracy and reduce statistical and systematic errors, we switched to substantially lower temperatures in comparison with our previous study (0.0625 and 0.125 eV in comparison with 0.5 eV in Boyda et al. (2016)). According to the general idea of Euclidean time formalism, it leads to increased number of time slices. Unfortunately, it’s impossible to rely on the old version of the BG method in this case, since the regularization scheme used in Brandt et al. (2015) and Boyda et al. (2016) doesn’t work well for lattices with large ($\sim 100$) number of Euclidean time slices. Significant modifications of the algorithm were introduced in order to increase it’s stability and to improve the resolution in frequency. These modifications are discussed below alongside with the careful study of possible systematic errors.

We start from the brief description of the BG method. The estimator of the spectral function $\sigma(\omega)$ is defined as the convolution of the exact spectral function $\tilde{\sigma}(\tilde{\omega})$ with the resolution function $\delta(\omega,\tilde{\omega})$.

[TABLE]

The resolution function is defined as a linear combination of the kernel profiles:

[TABLE]

The coefficients $q_{j}(\omega)$ are determined by minimizing the width of resolution function around $\omega$: $\partial_{q_{j}}D=0$, where $D$ is defined as

[TABLE]

In practice, we minimize the width imposing the normalization condition:

[TABLE]

and sometimes the additional condition $\delta(\omega,0)=0$. The second condition is introduced in order to eliminate the influence of Drude peak on the estimator for optical conductivity. The result of this minimization yields

[TABLE]

where

[TABLE]

The matrix $W$ is extremely ill-conditioned, with $C(W)\equiv\frac{\lambda_{\text{max}}}{\lambda_{\text{min}}}\approx O(10^{20})$. Thus, one needs to regularize the method in order to obtain stable results for given set of data $G(\tau_{i})$.

Previous studies employing the BG algorithm Brandt et al. (2015); Boyda et al. (2016) have used the regularization scheme based on the addition of the covariance matrix $C_{j,k}$ of the Euclidean correlator $G$ to the kernel (66):

[TABLE]

where $\lambda$ is a small regularization parameter. The method worked well in Boyda et al. (2016) for lattices with 20 Euclidean time steps, but we failed with this approach in the calculations at small temperatures where the typical number of time steps is of the order of 100. Moreover, there is an important type of calculations where we apply the analytical continuation to formally exact correlator computed for the free fermions. It should be done in order to check the validity of the method and to study the systematic errors. Despite the fact that the numerical values of exact data contain only round-off errors, it’s still impossible to use them in the Green-Kubo relation without regularization. Unfortunately, the covariance matrix is not defined for this kind of data, thus we again need some alternative regularization methods.

In this paper we used combination of the so-called Tikhonov regularization and averaging of the correlator over intervals in Euclidean time. The Tikhonov regularization is widely used for the ill-posed problems of the form $Ax=b$. In this method, one seeks a solution to the modified least-squares function

[TABLE]

where $\Gamma$ is an appropriately chosen matrix. In the standard Tikhonov regularization we choose $\Gamma=\lambda{\bm{1}}$, where $\lambda$ is again a small regularization parameter. As usually, for small $\lambda$ the solutions fit the data well but are oscillatory, while at large $\lambda$ the solutions are smooth but do not fit the data.

In practice the Tikhonov regularization is introduced during the inversion of the kernel matrix $W$, where we use the singular value decomposition (SVD):

[TABLE]

where $\Sigma=\text{diag}(\sigma_{1},\sigma_{2},\dots,\sigma_{N}),~{}\sigma_{1}\geq\sigma_{2}\geq\cdots\geq\sigma_{N}$. The inverse is thus easily expressed as

[TABLE]

Applying standard Tikhonov regularization simply modifies the matrix $\Sigma$ in the following way

[TABLE]

Thus, one can see that the small eigenvalues which satisfy $\lambda\gg\sigma_{i}$ are smoothly cut off.

Now let’s turn to the second regularization technique, namely, the averaging of the correlator over some intervals in time. In this procedure, we take the correlator data $\{G(\tau_{i});~{}i=0,1,\dots,N_{\tau}-1\}$ and map this to a new set $\{\tilde{G}(\tilde{\tau}_{j});~{}j=1,\dots,N_{\text{int}}\}$ where

[TABLE]

Due to the linearity of the Green-Kubo relation, one can construct $\{\bar{K}(\bar{\tau}_{j});~{}j=1,\dots,N_{\text{int}}\}$ in an analogous manner and use the new kernel and the new correlator in the BG method without any further modifications of the algorithm. The averaging increases the signal-to noise ratio in the input data. From the point of view of the ill-defined inversion problem for the kernel matrix $W$, we reduce the matrix size thus reducing the number of extremely small numbers $\sigma_{i}$ in $\Sigma$ (see 69 and 70) and making the SVD decomposition more stable.

We use lattices with spatial sizes $24\times 24$, $36\times 36$, $48\times 48$, $72\times 72$ and $96\times 96$ (the two largest sizes only in the test calculations for the free fermions). In the most of our QMC runs the temperature is equal to 0.125 eV with 80 steps in Euclidean time. As it was shown in Buividovich et al. (2016) this discretization is good enough if we are far from antiferromagnetic phase transition. In analytical continuation we use Tikhonov regularization with constant $\lambda=10^{-12}...10^{-11}$ and additional averaging over Euclidean time starting from the 10th step in time. The length of averaging intervals is equal to 10 time slices.

In the beginning we should study the systematic errors appeared due to finite lattice size, non-zero temperature and regularization during analytical continuation. We start from the check that the method applied to the free current-current correlator really reproduces the analytic profile of conductivity for the free tight-binding model. Fig. 8 shows the results of this study. First of all, one can see that the profiles become stable in the limit of large lattices. Special attention should be paid to the plateau with center around $\omega=0.7t$ (the point is defined as a position corresponding to the minimal value of the derivative $d\sigma/d\omega$). In the limit of large lattices this plateau reproduces correct value of non-interacting conductivity at $\omega=0.7t$ ($\sigma_{NI}=1.058\sigma_{0}$) with systematic error less than one percent: see fig. 11.

Another argument why we should look at this plateau is that it evolves into the plateau corresponding to the Dirac fermions in the low temperature limit. This feature is demonstrated in the figure 9. Indeed, once the temperature is decreased, the plateau shifts towards lower frequencies. The value of $\sigma(\omega)$ at plateau simultaneously goes towards the standard limit of 2D Dirac fermions $\sigma_{0}$. Thus we can conclude that the correct value of conductivity for Dirac fermions can be reproduced in our approach once we simultaneously increase the lattice size and decrease the temperature.

However, in practice it’s impossible to make QMC simulation at very large lattices and very small temperature (the latter needs increased number of steps in Euclidean time). Thus we’ll stop at lattice size equal to $48\times 48$, where correct value of $\sigma|_{\omega=0.7t}$ for the free fermions is already reproduced with good precision once we perform calculations at the temperature T=0.125 eV. This frequency is used in all further calculations in the paper where we compare QMC results with analytic calculations.

The last potential source of systematic errors is the regularization. We checked the dependence of conductivity on the constant $\lambda$ in the whole interval $\lambda=10^{-12}...10^{-11}$ and the relative variation of conductivity $\sigma|_{\omega=0.7t}$ is smaller than 1%. Thus we can conclude that systematic errors are under control since all of them are actually smaller than statistical uncertainties.

Conductivity profiles obtained in real QMC calculations are shown in the figure 10. Once can see that the general feature of $\sigma(\omega)$ profiles: the plateau around $\omega=0.7t$, is preserved. Actually, for suspended graphene the plateau is even substantially wider. Results for conductivity at $\omega=0.7t$ are summarized in the figure 11. Interestingly, the finite size effects for interacting systems become milder with increased interaction strength. This feature provides additional support for the statement that the finite size effects are under control, because the case of free fermions shows the upper limit for the finite size effects. Results for the largest accessible in QMC lattice ($48\times 48$) are used in the main text for comparison with other approaches.