Impact of the electron-electron correlation on phonon dispersions:failure of LDA and GGA functionals in graphene and graphite

Impact of the electron-electron correlation on phonon dispersions: failure of LDA and GGA functionals in graphene and graphite


We compute electron-phonon coupling (EPC) of selected phonon modes in graphene and graphite using various ab-initio methods. The inclusion of non-local exchange-correlation effects within the GW approach strongly renormalizes the square EPC of the A K mode by almost 80 with respect to density functional theory in the LDA and GGA approximations. Within GW, the phonon slope of the A K mode is almost two times larger than in GGA and LDA, in agreement with phonon dispersions from inelastic x-ray scattering and Raman spectroscopy. The hybrid B3LYP functional overestimates the EPC at K by about 30%. Within the Hartree-Fock approximation, the graphene structure displays an instability under a distortion following the A phonon at K.

71.15.Mb, 63.20.kd, 78.30.Na, 81.05.Uw

The electron-phonon coupling (EPC) is one of the fundamental quantities in condensed matter. It determines phonon-dispersions and Kohn anomalies, phonon-mediated superconductivity, electrical resistivity, Jahn-Teller distortions etc. Nowadays, density functional theory within local and semi-local approximations (DFT) is considered the ”standard model” to compute ab-initio the electron-phonon interaction and phonon dispersions (1). Thus, a failure of DFT would have major consequences in a broad context. In GGA and LDA approximations (2), the electron exchange-correlation energy is a local functional of the charge density and the long-range character of the electron-electron interaction is neglected. These effects are taken into account by Green-function approaches based on the screened electron-electron interaction W, such as the GW method (3). GW is considered the most precise ab-initio approach to determine electronic bands but, so far, it has never been used to compute EPCs nor phonon dispersions. The semi-empirical B3LYP functional (2) partially includes long-range Hartree-Fock exchange. B3LYP has been used to compute phonon frequencies but, so far, not the electron-phonon coupling.

The electron-phonon coupling is a key quantity for graphene, graphite and carbon nanotubes. It determines the Raman spectrum, which is the most common characterization technique for graphene and nanotubes (4); (5), and the high-bias electron transport in nanotubes (6). Graphene and graphite are quite unique systems in which the actual value of the EPC for some phonons can be obtained almost directly from measurements. In particular, the square of the EPC of the highest optical phonon branch (HOB) at the symmetry K-point is proportional to the HOB slope near K (7). The HOB K slope can be measured by inelastic x-ray scattering (IXS) (8); (9) or by the dispersion of the D and 2D lines as a function of the excitation energy in a Raman experiment  (10); (11); (12); (13); (5). A careful look at the most recent data suggests that the experimental phonon slopes (and thus the EPC) are underestimated by DFT (5). The ability of DFT (LDA and GGA) in describing the EPC of graphene was also questioned by a recent theoretical work (14).

Here, we show that: i) the GW approach, which provides the most accurate ab-initio treatment of electron-correlation, can be used to compute the electron-phonon interaction and the phonon dispersion; ii) in graphite and graphene, DFT (LDA and GGA) underestimates by a factor 2 the slope of the phonon dispersion of the highest optical branch at the zone-boundary and the square of its electron-phonon coupling by almost 80%; iii) GW reproduces both the experimental phonon dispersion near K, the value of the EPC and the electronic band dispersion; iv) the B3LYP hybrid functional (2) gives phonons close to GW but overestimates the EPC at K by about 30 %; v) within Hartree-Fock the graphite structure is unstable.

In Fig. 1, we show the phonon dispersion of graphite computed with DFT  (15). In spite of the general good agreement with IXS data, the situation is not clear for the HOB near K. In fact, despite the scattering among experimental data, the theoretical HOB is always higher in energy with respect to measurements and the theoretical phonon slope (for the HOB near K) is underestimating the measured one. It is also remarkable that while the DFT K frequency is  1300 cm, the highest measured is much lower at  1200 cm.

Figure 1: Upper panel: Phonon dispersion of graphite. Lines are DFT calculations, dots and triangles are IXS measurements from Refs. (8); (9), respectively. Lower panel: phonon dispersion of graphene from DFT calculations. Dashed lines are obtained by subtracting from the dynamical matrix the phonon self-energy between the bands ( in the text).

The dispersion of the HOB near K can also be obtained by Raman measurements of the graphene and graphite D-line (1350 cm(12). The D-line frequency depends on the energy of the exciting laser . According to the double-resonance model (13); (12), activates a phonon of the HOB with momentum q=K+q along the K-M line (5) and energy . q is determined by where is the energy of the electronic state with momentum k. Thus, by measuring vs. and considering the electronic bands dispersion from DFT one can obtain the phonon dispersion vs. q (12). The phonon dispersion thus obtained is very similar to the one from IXS data and its slope is clearly underestimated by DFT (Fig. 2, upper panel). The same conclusion is reached by comparing the D-line dispersion vs. (directly obtained from measurements) with calculations (Fig. 2, lower panel). Note that the dispersions of the Raman 2D-line (5) is consistent with the dispersion of the D and thus in disagreement with DFT (LDA and GGA) as well.

Figure 2: (Color online). Upper panel: dispersion of the highest optical phonon in graphite near K. Calculations are from DFT, or corrected to include GW renormalization of the electron-phonon coupling. Here, the DFT dispersion is vertically shifted by -40 cm to fit measurements. Dots and triangles are IXS data from Refs. (8); (9), respectively. Squares, plus and diamonds are obtained from Raman data of Refs.  (10); (11); (12), respectively, using the double resonance model (13); (12). Lower panel: dispersion of the Raman D-line.

The steep slope of the HOB near K is due to the presence of a Kohn anomaly for this phonon (7). In particular, in Ref. (7), it was shown that the HOB slope is entirely determined by the contribution of the phonon self-energy between -bands, , to the dynamical matrix, . is the phonon pulsation, where is the mass. For a given phonon with momentum q,


where the sum is performed on wavevectors all over the Brillouin zone, is the EPC, is the derivative of the Kohn-Sham potential with respect to the phonon mode, is the Bloch eigenstate with momentum k, band index and energy . () identifies the occupied (empty) -band. In Fig. 1 we show a fictitious phonon dispersion obtained subtracting from the dynamical matrix () for each phonon. The HOB is the branch which is mostly affected and, for the HOB, becomes almost flat near K. Thus, DFT (LDA or GGA) fails in describing the HOB slope near K, slope which is determined by . is given by the the square EPC divided by -band energies. Thus, the DFT failure can be attributed to a poor description of the EPC or of the -band dispersion.

In graphene and graphite, it is known that standard DFT provides an underestimation of the and -band slopes of (17); (18). A very precise description of the bands, in better agreement with measurements, is obtained using GW (17); (18). We thus computed the -bands with DFT (both LDA and GGA) (16) and GW (19) and compared with Hartree-Fock (HF) (20) and B3LYP (20). Details are in  (21). The different methods provide band dispersions whose overall behavior can be described by a scaling of the energies (17). The different scaling factors can be obtained by comparing : the energy difference between the and bands at the symmetry point M (L) for graphene (graphite). is larger in GW than in DFT (Tab. 1). Thus, inclusion of the GW correction to the electronic bands alone results in a larger denominator in Eq. 1, providing a smaller phonon slope and a worse agreement with experiments. The underestimation of the K phonon slope in DFT is, thus, due to the EPC.

Figure 3: a, b): patterns of the -E and K-A phonons of graphene. Dotted and dashed lines are the Wigner-Seitz cells of the unit-cell and of the super-cell. c): Hartree-Fock equilibrium structure.

The EPC can be computed with linear response as, e.g., in Ref. (7) but, at present, the use of this technique within GW is not feasible. Alternatively, the EPC associated to a phonon mode can be determined by the variation of the electronic band energies by displacing the atoms according to the considered mode. In graphene, at K, there are doubly degenerate electronic states at the Fermi level. The HOB corresponds to the E phonon at and to the A at K. As an example, we consider the EPC associated to the -E phonon and we displace the atoms according to its phonon pattern (see Fig. 3). Following symmetry arguments (22), one can show that, in an arbitrary base of the two-dimensional space of the bands at K, the Hamiltonian is the 22 matrix:


where each atom is displaced by , , and , where the sum is performed on the two degenerate bands. Diagonalizing Eq. 2, we see that an atomic displacement following the -E phonon induces the splitting and


In analogous way, we define for the A phonon at K. Let us consider a graphene supercell. Such a cell can be used to displace the atoms following the K-A phonon (Fig. 3), since the K point is refolded in . Let us call the splitting of the bands induced by this displacement (since K is refolded in , here denotes the energies of the band of the supercell corresponding to the band at K in the unit cell). Considering the atomic distortion of Fig. 3, displacing each atom by , one can show that


In practice, by calculating band energies in the distorted structures of Fig. 3 and using Eqs. 34 one obtains the EPCs of the -E and K-A phonons between states. Similar equations can be used for graphite (23). Results are in Tab. 1 together with the computed phonon frequencies. The EPCs from DFT are in agreement with those from linear response (7). We also remark that, within the present ”frozen-phonon” approach, the Coulomb vertex-corrections are implicitly included within GW.

To study the effect of the different computational methods on the the phonon slope (which is determined by ) we recall that is the ratio of the square EPC and band energies (Eq. 1). Thus, we have to compare . As an example, assuming that the change of from DFT to GW is constant for q near K,


and provides the change in the K phonon slope going from DFT to GW. To understand the results, we recall that in standard DFT the exchange-correlation depends only on the local electron-density. In contrast, the exchange-interaction in HF and GW is non-local. Furthermore, in GW, correlation effects are non-local since they are described through a dynamically screened Coulomb interaction. The hybrid functional B3LYP gives results intermediate between DFT and HF.

DFT 4.03 44.4 11.0 1568 89.9 22.3 1275
DFT 4.08 45.4 11.1 1583 92.0 22.5 1303
GW 4.89 62.8 12.8 193 39.5
B3LYP 6.14 82.3 13.4 1588 256 41.7 1172
HF 12.1 321 26.6 1705 6020 498 960
DFT 4.06 43.6 10.7 1568 88.9 21.8 1299
DFT 4.07 44.9 11.0 1581 91.5 22.5 1319
GW 4.57 58.6 12.8 164.2 35.9 (1192)
Table 1: Electron-phonon coupling of the -E and K-A phonons computed with various approximations. (eV), (eV) and (eV/Å) are defined in the text. () is the phonon frequency of the E (A) mode (cm). The GW for graphite (in parenthesis) is not computed directly (see the text). is the imaginary unit.

Both and are heavily overestimated by HF, the K-EPC being so huge that graphene is no more stable (the KA phonon frequency is not real). Indeed, the HF equilibrium geometry is a reconstruction with alternating double and single bonds of 1.40 and 1.43 Å lengths as in Fig. 3 (with a gain of 0.9 meV/atom). These results demonstrate the major effect of the long-range character of the exchange for the K-EPC (14) but also the importance of the proper inclusion of the screening (included in GW but neglected in HF). Notice also that of graphite is smaller with respect to graphene by 10%. This is explained by the larger screening of the exchange in graphite (due to the presence of adjacent layers) than in graphene. On the contrary, within GGA and LDA, the graphite phonon frequencies and EPCs are very similar to those of graphene, since these functionals do not take into account the electron-electron interaction screening.

Concerning the phonon slope, is 15% larger than . Indeed, DFT reproduces with this precision the phonon frequency and dispersion of the HOB at . On the contrary, is 60% larger than , for graphite. This large increase with respect to DFT could explain the disagreement between DFT and the measured A phonon dispersion near K. To test this, we need to determine the GW phonon dispersion that, using Eq. 5 becomes , where . Moreover, we can assume since the component of the dynamical matrix (Eq. 1) is not expected to have an important dependence on q (Fig. 1). The value of is obtained as a fit to the measurements of Fig. 2 (24). The resulting K A phonon frequency is 1192 cm which is our best estimation and is almost 100 cm smaller than in DFT. The phonon dispersion thus obtained and the corresponding D-line dispersion are both in better agreement with measurements (Fig. 2).

The partial inclusion of long-range exchange within the semiempirical B3LYP functional leads to a strong increase of the EPC at K as compared to the LDA and GGA functionals. However, comparing to the GW value, the EPC is overestimated by 30% and the corresponding frequency for the K-A mode at 1172 cm falls well below the degenerate K-mode which is around 1200 cm in the experiment (8); (9) (Fig. 1) and at 1228 cm in our phonon calculation with B3LYP. We have checked that tuning the percentage of HF-exchange in the hybrid functional allows to match the EPC value of the GW approach (in which case, the K-A mode remains the highest mode. This may be a good way to calculate the full phonon dispersion of graphite/graphene within DFT, yet with an accuracy close to the one of the GW approach.

Concluding, GW is a general approach to compute accurate electron-phonon coupling where DFT functionals fail. Such a failure in graphite/graphene is due to the interplay between the two-dimensional Dirac-like band structure and the long-range character of the Coulomb interaction (14). However, GW can be also used in cases (in which the EPC is badly described by DFT) where the electron-correlation is short ranged (25).

Calculations were done at IDRIS (081202, 081827). C.A. and L.W. acknowledge French ANR PJC05_6741. We thank D.M. Basko, A. Rubio, J. Schamps, and C. Brouder for discussions and A. Marini for the code Yambo.


  1. S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, Rev. Mod. Phys. 73, 515 (2001).
  2. LDA, GGA and B3LYP refer, respectively, to D.M. Ceperley and B.J. Alder, Phys. Rev. Lett. 45, 566 (1980); J.P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996); A.D. Becke, J. Chem. Phys. 98, 5648 (1993).
  3. F. Aryasetiawan and D. Gunnarsson, Rep. Progr. Phys. 61, 237 (1998); W.G. Aulbur, L. Jönsson, and J.W. Wilkins, Solid State Phys. 54, 1 (2000); G. Onida, L. Reining, and A. Rubio, Rev. Mod. Phys. 74, 601 (2002).
  4. S. Reich and C. Thomsen, Phil. Trans. R. Soc. London A 362, 2271 (2004).
  5. A.C. Ferrari, J.C. Meyer, V. Scardaci, C. Casiraghi, M. Lazzeri, F. Mauri, S. Piscanec, D. Jiang, K.S. Novoselov, and A.K. Geim, Phys. Rev. Lett. 97, 187401 (2006); D. Graf, F. Molitor, K. Ensslin, C. Stampfer, A. Jungen, C. Hierold, and L. Wirtz, Nano Lett. 7, 238 (2007).
  6. Z. Yao, C. L. Kane, and C. Dekker, Phys. Rev. Lett. 84, 2941 (2000).
  7. S. Piscanec, M. Lazzeri, F. Mauri, A.C. Ferrari, and J. Robertson, Phys. Rev. Lett. 93, 185503 (2004).
  8. J. Maultzsch, S. Reich, C. Thomsen, H. Requardt, and P. Ordejón, Phys. Rev. Lett. 92, 075501 (2004).
  9. M. Mohr, J. Maultzsch, E. Dobardžić, I. Milošević, M. Damnjanović, A. Bosak, M. Krisch, and C. Thomsen, Phys. Rev. B 76, 035439 (2007).
  10. I. Pócsik, M. Hundhausen, M. Koós, and L. Ley, J. Non-Cryst. Solids 227, 1083 (1998).
  11. P. Tan, L. An, L. Liu, Z. Guo, R. Czerw, D.L. Carrol, P.M. Ajayan, N. Zhang, and H. Guo, Phys. Rev. B 66, 245410 (2002).
  12. J. Maultzsch, S. Reich, and C. Thomsen, Phys. Rev. B 70, 155403 (2004).
  13. C. Thomsen and S. Reich, Phys. Rev. Lett. 85, 5214 (2000).
  14. D. M. Basko and I. L. Aleiner, Phys. Rev. B 77, 041409(R) (2008).
  15. Calculations were done as in Ref. (1). Technical details and thus phonon dispersions are the same as in Ref. (7).
  16. LDA and GGA calculations were done with the code PWSCF (S. Baroni et al.,, with pseudopotentials of the type N. Troullier and J. L. Martins, Phys. Rev. B 43, 1993 (1991).
  17. S.Y. Zhou, G.H. Gweon, C.D. Spataru, J. Graf, D.H. Lee, S.G. Louie, and A. Lanzara Phys. Rev. B 71, 161403(R) (2005); S.G. Louie, in Topics in Computational Materials Science, C.Y. Fong editor (World Scientific, Singapore, 1997), p.96.
  18. A. Grüneis, C. Attaccalite, T. Pichler, V. Zabolotnyy, H. Shiozawa, S. L. Molodtsov, D. Inosov, A. Koitzsch, M. Knupfer, J. Schiessling, R. Follath, R. Weber, P. Rudolf, L. Wirtz, and A. Rubio Phys. Rev. Lett. 100, 037601 (2008).
  19. GW calculations were done with the code Yambo (A. Marini et al., the Yambo project,, within the non-self consistent GW approximation, starting from DFT-LDA wave-functions and using a plasmon-pole model for the screening, following M.S. Hybertsen and S.G. Louie, Phys. Rev. B 34, 5390 (1986).
  20. HF and B3LYP calculations were done with the code CRYSTAL (V.R. Saunders et al. CRYSTAL03 User’s Manual, University of Torino, Torino, 2003), using the TZ basis by Dunning (without the diffuse P-function).
  21. For graphite we used the experimental lattice parameters (=2.46 Å, =6.708 Å). For graphene we used =2.46 Å  and a vacuum layer of 20 a.u.. EPCs were calculated on a structure distorted by 0.01 a.u. For graphene, the electronic integration on the 11 cell was done with a 18181 grid for LDA/GGA, 36361 for GW, 66661 for B3LYP/HF. For graphite it was 18186. For the cell we used the nearest equivalent k-grid. Plane-waves are expanded up to 60 Ry cut-off. We used a Fermi-Dirac smearing with 0.002 Ry width for B3LYP/HF/GW and a Gaussian smearing with 0.02 Ry width for LDA/GGA.
  22. J.C. Slonczewski and P.R. Weiss Phys. Rev. 109, 272 (1958). See also Suppl. information to S. Pisana, M. Lazzeri, C. Casiraghi, K.S. Novoselov, A.K. Geim, A.C. Ferrari, and F. Mauri, Nature Materials 6, 198 (2007).
  23. In graphite, at the high-symmetry H point the four bands are degenerate two-by-two, being the energy difference. By displacing the atoms according to the E phonon, these bands remain degenerate and the energy difference is increased by . In analogy to Eqs. 34 we define . By displacing the atoms according to the K A phonon, the four bands are no longer degenerate, being () the two bands which are up(down)-shifted and . We define , where indicates the average between the four possible couples.
  24. In principle, a direct calculation of (and thus of ) could be obtained, e.g., by finite differences from a prohibitively expensive GW total energy calculation.
  25. P. Zhang, S.G. Louie and M.L. Cohen, Phys. Rev. Lett. 98, 067005 (2007).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description