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.
pacs: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.
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.
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.
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. 3, 4 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.
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.
- S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, Rev. Mod. Phys. 73, 515 (2001).
- 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).
- 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).
- S. Reich and C. Thomsen, Phil. Trans. R. Soc. London A 362, 2271 (2004).
- 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).
- Z. Yao, C. L. Kane, and C. Dekker, Phys. Rev. Lett. 84, 2941 (2000).
- S. Piscanec, M. Lazzeri, F. Mauri, A.C. Ferrari, and J. Robertson, Phys. Rev. Lett. 93, 185503 (2004).
- J. Maultzsch, S. Reich, C. Thomsen, H. Requardt, and P. Ordejón, Phys. Rev. Lett. 92, 075501 (2004).
- 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).
- I. Pócsik, M. Hundhausen, M. Koós, and L. Ley, J. Non-Cryst. Solids 227, 1083 (1998).
- 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).
- J. Maultzsch, S. Reich, and C. Thomsen, Phys. Rev. B 70, 155403 (2004).
- C. Thomsen and S. Reich, Phys. Rev. Lett. 85, 5214 (2000).
- D. M. Basko and I. L. Aleiner, Phys. Rev. B 77, 041409(R) (2008).
- Calculations were done as in Ref. (1). Technical details and thus phonon dispersions are the same as in Ref. (7).
- LDA and GGA calculations were done with the code PWSCF (S. Baroni et al., http://www.quantum-espresso.org), with pseudopotentials of the type N. Troullier and J. L. Martins, Phys. Rev. B 43, 1993 (1991).
- 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.
- 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).
- GW calculations were done with the code Yambo (A. Marini et al., the Yambo project, http://www.yambo-code.org/), 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).
- 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).
- 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.
- 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).
- 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. 3- 4 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.
- 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.
- P. Zhang, S.G. Louie and M.L. Cohen, Phys. Rev. Lett. 98, 067005 (2007).