# Monte-Carlo study of quasiparticle dispersion relation in monolayer graphene

###### Abstract

The density of electronic one-particle states in monolayer graphene is studied by performing the Hybrid Monte-Carlo simulations of the tight-binding model for electrons on the orbitals of carbon atoms which make up the graphene lattice. Density of states is approximated as a derivative of the number of particles over the chemical potential at sufficiently small temperature. Simulations are performed in the partially quenched approximation, in which virtual particles and holes have zero chemical potential. It is found that the Van Hove singularity becomes much sharper than in the free tight-binding model. Simulation results also suggest that the Fermi velocity increases with interaction strength up to the transition to the phase with spontaneously broken chiral symmetry.

Monte-Carlo study of quasiparticle dispersion relation in monolayer graphene

P. V. Buividovich^{†}^{†}thanks: Speaker. ^{†}^{†}thanks: This work was supported by the S. Kowalewskaja award from the Alexander von Humboldt Foundation.

Institute for Theoretical Physics, Regensburg University, Universitätsstrasse 31, D-93053 Regensburg, Germany

\abstract@cs

## 1 Introduction

Electron transport properties of graphene are of great importance both for the industrial applications of this novel material as well as for our understanding of the physics of strongly correlated electrons. It turns out that at low energies electronic excitations in graphene effectively behave as free massless Dirac fermions [1]. However, they propagate with the Fermi velocity , which is much less than the speed of light . There is also a natural way to make these Dirac fermions massive: one should break the symmetry between the two simple rhombic sublattices of hexagonal graphene lattice by introducing the so-called “staggered potential”, so that the potential is lifted by some value on the sites of one sublattice and lowered by on the sites of other sublattice. This “staggered potential” then plays the role of the Dirac mass. At higher energies of order of the hopping energy the quasiparticle dispersion relation becomes nonlinear. At , it has a saddle point, which results in the logarithmic divergence in the density of states. This Van Hove singularity leads to the so-called Lifshitz transition as the Fermi energy of the system approaches the saddle point. The total width of the valence band is . The dispersion relation and the density of states at different values of the staggered potential are illustrated on Fig. 1.

Charge carriers in graphene are also subject to electromagnetic interactions. The treatment of these interactions can be significantly simplified due to the smallness of the Fermi velocity. First, retardation effects can be neglected and, second, the Lorentz force between the electrons also is suppressed by a factor as compared to the electric force. Thus it is sufficient to consider only the instantaneous Coulomb interaction between Dirac quasiparticles.

An important question is how the parameters of the free tight-binding model such as the widths of the valence band, the Fermi velocity and the mass gap change when one takes into account the interactions between quasiparticles. In order to pose such a question, one has to adapt the Fermi-Landau liquid model and assume that despite the interactions, in some sense quasiparticles can be still described as free fermions with some modified dispersion relation. Strictly speaking, graphene is an example of the so-called marginal Fermi liquid [2], for which the Fermi-Landau liquid approximation becomes inapplicable in the vicinity of the Fermi point (Fermi point is a degenerate Fermi surface with zero radius) . However, as discussed in [3], even a small doping brings monolayer graphene back into the Fermi-Landau liquid regime. One can expect that the introduction of the “staggered potential”, which eliminates the Fermi point, should have the same effect.

In this paper the density of one-particle states in monolayer graphene is studied numerically. It is found that the Lifshitz phase transition, associated with the Van Hove singularity around the saddle point in the dispersion relation of the tight-binding model, becomes much more pronounced in the interacting theory. Such behavior is in agreement with the analytical calculations based on renormalization-group arguments [4]. There are also indications that the Fermi velocity increases with the interaction strength up to the transition to the phase with spontaneously broken chiral symmetry, again in agreement with renormalization-group arguments [5]. On the other hand, it seems that the width of the artificially induced (with the help of the staggered potential) mass gap remains practically constant. However, measurements with a better energy resolution are required to clarify the situation completely.

In analytical calculations, it is most convenient to study the quasiparticle dispersion relation and the corresponding density of states by analyzing the poles of the fermionic Green functions. In numerical Monte-Carlo simulations, however, fermionic Green functions can only be calculated for a finite number of discrete values of imaginary (Euclidean) time. Extraction of the real-time Green functions from such Euclidean field correlators is an ill-defined numerical problem, which has no unique solution. The most commonly used method for extracting real-time quantities from Euclidean correlators is the Maximal Entropy Method (MEM) [6]. However, this method introduces significant systematic uncertainties, in particular, due to ambiguity of the choice of the “model function” which incorporates our prior knowledge about the expected form of the corresponding spectral function. Typically, MEM tends also to smear the singularities (such as thresholds or Van Hove singularities) of the spectral functions at the scales of order of temperature. An attempt to extract the AC conductivity of graphene from numerically calculated current-current correlators was already reported by the author in [7]. It was found that MEM indeed was not able to reproduce the AC conductivity with precision which would be sufficient for quantitative analysis.

It should be noticed, however, that in a strict sense the density of states which can be formally extracted from the fermionic Green function does not correspond to the density of any eigenstates of the interacting Hamiltonian. The reason is that the notion of quasiparticle is not strictly defined in this case. The interpretation of the spectral function of interacting fermionic gas in terms of the density of quasiparticle states only makes sense in the framework of the Landau-Fermi liquid model. However, in this framework one can also think of many other possible definitions of the density of states, which are only constrained by the requirement to correctly reproduce the density of states for a free Fermi gas with an arbitrary dispersion relation for fermions.

In this work, the density of states is defined as the derivative of the number of particles over chemical potential. For a gas of free fermions with the density of one-particle states , the particle density and its derivative can be written as:

(1.0) | |||

(1.0) |

In the limit of zero temperature, the Fermi factor becomes the Heaviside step function , and its derivative with respect to the chemical potential becomes the -function . Thus in the limit of zero temperature the derivative is exactly the density of one-particle states. At finite temperature the -function is smeared over a finite energy range with the width of order of temperature and becomes an exponentially decaying function. A direct estimate of the distribution width yields:

(1.0) |

Thus by performing simulations at sufficiently small temperatures and by measuring the derivative , one can obtain a reasonably accurate estimate of the density of one-particle states. The resulting function is smeared as compared to the zero-temperature limit, similarly to the results which can be obtained by using MEM. On the other hand, such an estimate is not based on any “model function” and is thus free from possible bias in the measurements caused by the particular choice of this function.

In practice, a direct numerical implementation of such measurements is a formidably difficult task due to the so-called “sign problem” of the Monte-Carlo simulations at finite chemical potential [8]. At zero chemical potential (that is, for graphene at half-filling) the sign problem is absent due to the symmetry between particles and holes [9, 7, 10]. For this reason, in this work the quantity is calculated in the partially quenched approximation, in which virtual particles are not influenced by finite chemical potential. Such approximation can be introduced as follows: consider the tight-binding model with independent fermion flavors and assume that chemical potential is nonzero only for fermions. For the time being it is convenient to assume that the staggered potential is equal to zero. Transforming the partition function of such a system into the path integral representation along the lines of [7, 9, 10], it is straightforward to obtain the following result:

(1.0) |

where is the Hubbard-Stratonovich [10] or the electrostatic potential [7, 9] field with the action , is the Euclidean time, variable labels the sites of the graphene hexagonal lattice, is the fermionic hopping operator and is the one-particle Hamiltonian of the tight-binding model. Consider now the linear response of the system to adding more flavours of fermions at finite chemical potential, that is, the derivative . This yields the partially quenched partition function

(1.0) |

This expression was derived by taking into account the invariance of the path integral under the reflection and the anti-unitary symmetry of the operator at zero chemical potential: with , where the plus and the minus signs are taken for even and odd sites of the graphene hexagonal lattice, respectively. Clearly, the square of is just the identity operator. At low energies, when quasiparticles in graphene can be described as Dirac fermions, plays the role of the Dirac matrix.

Finally, taking the derivatives of the quenched partition function over the chemical potential, one finds the partially quenched approximation for and :

(1.0) | |||

(1.0) |

where is the number of lattice sites in the system. These expressions differ from the expression for the full, unquenched functions and by the absence of chemical potential in the factor . However, since the density of one-particle states is not a rigorously defined quantity in the interacting theory, there is no reason to believe that is a better estimate of then . Indeed, for a free fermion gas both definitions are completely equivalent.

The observables (1) and (1) are now well-suited for numerical simulations and can be calculated by averaging the quantities and over the equilibrium ensemble of the fields with the weight [9, 10, 7]. To this end Hybrid Monte-Carlo simulations of the tight-binding model of graphene on the lattice with toric topology which consisted of hexagonal cells were performed. Temporal size of the lattice is with lattice spacing . Coulomb interactions are modelled by coupling the tight-binding model to the noncompact gauge field in dimensions, as in [7]. The corresponding coupling constant is controlled by the dielectric permittivity of substrate on which the graphene monolayer is placed: , where [7]. Lattice size in the direction perpendicular to the graphene plane is . Simulation setup and the discretization of the fermionic operator are the same as in [7]. Despite the fact that the estimate of the density of states in (1) can be expressed in terms of the partially quenched partition function (1) only at zero staggered potential , the simulations have to be performed at some nonzero (in this work and ) in order to ensure the invertibility of the fermionic operator . The derivative is found by first calculating according to (1) for equidistant values of separated by and then numerically differentiating this function using the symmetric difference

(1.0) |

The resulting dependence of on is illustrated on Fig. 2 for and . For comparison, the function for the free theory and the density of states for free Dirac fermions (which corresponds to in the limit of zero temperature) are also plotted:

(1.0) |

This expression takes into account that on the hexagonal lattice with periodic boundary conditions the number of states in the element of the momentum space is given by [7].

First of all one can note the sharp rise of the peaks at , which are associated with the Van Hove singularities in the free theory. This is an indication that the Lifshitz phase transition (associated with the crossing of the Van Hove singularity by the Fermi level) might become much sharper in the interacting theory. Indeed, such phase transitions are known to be unstable with respect to even small interaction between particles. Application of renormalization-group techniques to quasiparticles in the vicinity of the Van-Hove singularities also show that as the interactions are switched on, the density of states tends to increase [4]. In other words, the quasiparticle dispersion relation becomes more flat in the vicinity of the point. This sharpening of the Van Hove singularity can be clearly seen for both values of starting from the smallest nonzero coupling constant (which corresponds to substrate dielectric permittivity ). It is also interesting to note a small asymmetry between the Van Hove singularities at and at and a slight shift of the position of the minimum of the density of states from to positive values of .

Another important prediction of the renormalization-group analysis, which has been recently confirmed experimentally [11], is the logarithmic divergence of the Fermi velocity in the interacting tight-binding model [5]. As follows from (1), the increase of Fermi velocity should result in a depletion of the density of states near . Superficially, since the slope of near clearly increases towards smaller , it seems that our results point to the decrease of the Fermi velocity with interaction strength. However it might well be that due to the smearing (1) this increase of the density of states is just caused by the strong growth of the density of states at the Van Hove singularity. In order to demonstrate that the density of states near indeed decreases, on Fig. 3 the minimal value of for in the range is plotted as a function of . This minimum is always situated close to . One can see that indeed this minimal value decreases as the Coulomb interaction becomes stronger for . This is an indication of the decrease of the density of states in the vicinity of . As demonstrated in [9, 7], at smaller the symmetry between simple sublattices of graphene hexagonal lattice breaks down spontaneously. In this phase there is no reason to believe that the Fermi-Landau liquid model is a good approximation. A glance at Fig. 2 suggest also that the width of the energy gap which is artificially induced by the staggered potential is practically independent of the interaction strength.

To conclude, the presented results are consistent both with the sharpening of the Van Hove singularity and with the increase of the Fermi velocity predicted by renormalization-group calculations [5, 4]. It seems, however, that the former modification the dispersion relation is much stronger than the latter, and that much better energy resolution is required in order to study quantitatively the evolution of the Fermi velocity and the mass gap in the tight-binding model with interactions. At present larger-scale simulations with significantly smaller temperatures (so that the smearing factor in (1) is much narrower) and with larger lattices are in progress, which would hopefully help to study the quasiparticle dispersion relation in more details.

## Acknowledgements

The author is grateful to Dr. M. I. Polikarpov, Dr. C. Popovici and Dr. M. V. Ulybyshev for many interesting discussions and also to Dr. L. von Smekal for his enlightening comments on the Van Hove singularity and the Lifshitz phase transition. Numerical calculations were performed at the ITEP computer systems “Graphyn” and “Stakan” (the author is grateful to A. V. Barylov, A. A. Golubev, V. A. Kolosov, I. E. Korolko, M. M. Sokolov for the help). The author was supported by the S. Kowalewskaja award from the Alexander von Humboldt foundation.

## References

- [1] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, A. K. Geim, Rev. Mod. Phys. 81 (2009), 109 - 162, [ArXiv:0709.1163]; G. W. Semenoff, Phys. Scr. 2012 (2012), 014016, [ArXiv:1108.2945]; G. W. Semenoff, Phys. Rev. Lett. 53 (1984), 2449 - 2452.
- [2] C. M. Varma, P. B. Littlewood, S. Schmitt-Rink, Phys. Rev. Lett. 63 (1989), 1996 - 1999; J. Gonzalez, F. Guinea, M. A. H. Vozmediano, Nucl. Phys. B 424 (1994), 595 - 618, [ArXiv:hep-th/9311105].
- [3] S. Das Sarma, E. H. Hwang, W. Tse, Phys. Rev. B 75 (2007), 121406, [ArXiv:cond-mat/0610581].
- [4] J. Gonzalez, F. Guinea, M. A. H. Vozmediano, Nucl. Phys. B 485 (1997), 694 - 726, [ArXiv:cond-mat/9602033].
- [5] R. Shankar, Rev. Mod. Phys. 66 (1994), 129 - 192; J. Gonzalez, F. Guinea, M. A. H. Vozmediano, Phys. Rev. B 59 (1999), 2474, [ArXiv:cond-mat/9807130]; M. S. Foster, I. L. Aleiner, Phys. Rev. B 77 (2008), 195413, [ArXiv:0802.0283].
- [6] M. Jarrell, J. E. Gubernatis, Phys. Rep. 269 (1996), 133 - 195; M. Asakawa, T. Hatsuda, Y. Nakahara, Prog. Part. Nucl. Phys. 46 (2001), 459, [ArXiv:hep-lat/0011040].
- [7] P. V. Buividovich, M. I. Polikarpov, Monte-Carlo study of the electron transport properties of monolayer graphene within the tight-binding model (2012), [ArXiv:1206.0619].
- [8] M. Troyer, U. Wiese, Phys. Rev. Lett. 94 (2005), 170201, [ArXiv:cond-mat/0408370]; P. de Forcrand, O. Philipsen, J. Phys. G: Nucl. Part. Phys. 35 (2008), 104098, [ArXiv:0807.0860].
- [9] J. E. Drut, T. A. Lähde, Phys. Rev. B 79 (2009), 165425, [ArXiv:0901.0584].
- [10] R. Brower, C. Rebbi, D. Schaich, PoS(Lattice 2011)056 (2012), [ArXiv:1204.5424].
- [11] D. C. Elias et al., Nature Phys. 7 (2011), 701 - 704, [ArXiv:1104.1396].