Ab initio quasiparticle bandstructure of ABA and ABC-stacked graphene trilayers

Ab initio quasiparticle bandstructure of ABA and ABC-stacked graphene trilayers


W̱e obtain the quasiparticle band structure of ABA and ABC-stacked graphene trilayers through ab initio density functional theory (DFT) and many-body quasiparticle calculations within the GW approximation. To interpret our results, we fit the DFT and GW bands to a low energy tight-binding model, which is found to reproduce very well the observed features near the K point. The values of the extracted hopping parameters are reported and compared with available theoretical and experimental data. For both stackings, the self energy corrections lead to a renormalization of the Fermi velocity, an effect also observed in previous calculations on monolayer graphene. They also increase the separation between the higher energy bands, which is proportional to the nearest neighbor interlayer hopping parameter . Both features are brought to closer agreement with experiment through the self energy corrections. Finally, other effects, such as trigonal warping, electron-hole asymmetry and energy gaps are discussed in terms of the associated parameters.


Graphene, a 2D sheet of carbon atoms in a honeycomb lattice, has attracted a lot of attention of the scientific community in the last few years due to its unique electronic properties, which lead to several potential applications in nanoelectronics Novoselov et al. (2004); Neto et al. (2009). However, since graphene is a zero-gap semiconductor, much of the current effort is directed in finding ways to open a gap for use in electronic devices. In particular, one way to do that is to consider graphene stacks, where a number of layers are stacked on top of each other with a particular arrangement. Much work has been done on bilayer graphene, where it was found that a tunable gap can be opened through application of an external electrical field or through doping Ohta et al. (2006); Zhang et al. (2009, 2008); Castro et al. (2007). In light of recent experimental progress, graphene trilayers are also attracting increasing attention, revealing electronic properties that depend on the stacking order of the three layers. The two most important stackings are ABA (Bernal) and ABC (rhombohedral), which are shown in Fig. 1. For ABA stacking, the low energy bands are predicted to consist of a set of monolayer and bilayer-like bands, with linear and quadratic dispersions, respectively Neto et al. (2009); Guinea et al. (2006); Partoens and Peeters (2006). Therefore, this trilayer is expected to show mixed properties from these two systems, which were already observed experimentally Taychatanapat et al. (2011). In the presence of an external electrical field perpendicular to the layers, these bands hybridize and a tunable overlap between the linear and parabolic bands is introduced Koshino and McCann (2009); Wu (2011). In contrast, for ABC stacking, the low energy bands consist of a pair of bands with cubic dispersion, which are very flat near the Fermi level. The large density of states associated with this behavior indicates that many-body interactions might play a crucial role in this case. In fact, there are already a few works in the literature investigating the possibility of different competing phases in this system, such as ferromagnetic order Olsen et al. (2013), charge and spin-density waves and quantum spin Hall phases Scherer et al. (2012), and even superconductivity Kopnin et al. (2011); Muñoz et al. (2013). Moreover, an external electrical field breaks the inversion symmetry of this trilayer and induces a tunable band gap, in a similar fashion to bilayer graphene Guinea et al. (2006); Zhang et al. (2010); Wu (2011) and also observed experimentally Lui et al. (2011); Zhang et al. (2011). We also point out that a similar type of dispersion and an associated quantum critical behavior was observed on a very different system, namely, the Laves phase of Neal et al. (2011); Alam and Johnson (2011).

Figure 1: Arrangement of carbon atoms on the ABA (left) and ABC (right) graphene trilayers. Atoms belonging to the A and B sublattices of each layer are represented by yellow and black spheres, respectively. The red dashed lines indicate the tight-binding parameters included in the model used to fit our first-principles results (see text).

In this work, we report first-principle calculations for ABA and ABC-stacked graphene trilayers, employing density functional theory (DFT) Hohenberg and Kohn (1964); Kohn and Sham (1965) and many-body quasiparticle corrections within the GW approximation. Following the framework of Hybertsen and Louie Hybertsen and Louie (1986), our calculations are done in two steps: First, a mean-field step is performed (DFT in our case), where the wavefunctions and energy bands are evaluated and stored. Then, in the next step, the quasiparticle corrections to the DFT energies are calculated within the GW approximation to the electron self energy, by using the DFT energies and wavefunctions 2. In our calculations, the DFT step is carried on using the LDA (Perdew-Zunger) approximation for the exchange-correlation functional Perdew and Zunger (1981), Troullier-Martins pseudopotentials for the electron-ion interaction Troullier and Martins (1991) and a plane wave basis set to expand the wavefunctions, with an energy cutoff of Ry, as implemented in the Quantum Espresso package et al (2009). The theoretical lattice constant for this pseudopotential is Å and the interlayer distance is set to the experimental value of graphite, Å. We also use a vacuum region of Å in the direction perpendicular to the layers in order to avoid interaction between periodic images. In the second step, the many-body calculations are perfomed following the scheme described below, which is implemented in the BerkeleyGW package Deslippe et al. (2012). For the calculation of the dielectric matrix in reciprocal space, we use an energy cutoff of Ry and a coarse Monkhorst-Pack k-point sampling Monkhorst and Pack (1976) with unnocuppied bands for general points away from zero and a fine grid and unnoccupied bands for . This matrix is first evaluated in the static limit by using the RPA approximation and then extended to non-zero frequencies by using a generalized plasmon pole (GPP) model. The self-energy operator is evaluated in the same coarse grid. In this step, the Coulomb interaction is truncated in the middle of the vacuum region of the slab. Finally, the bandstructure plots for the LDA and GW bands are generated along high-symmetry lines by using Wannier interpolation, as implemented in the Wannier90 package Mostofi et al. (2008).

The results of our calculations are shown in Fig. 2, where we compare the LDA and GW bandstructures in the energy range of the bands. For both calculations, we see that the bandstructures share the same qualitative features. For ABA stacking (left), the bandstructure near the Fermi level consist of a superposition of a pair of nearly linear bands, resembling those of monolayer graphene, and two pairs of parabolic bands, resembling those of bilayer graphene. However, unlike single and bilayer graphene, both sets of bands have small band gaps, due to the lack of inversion symmetry (or equivalently, A-B sublattice symmetry) in this trilayer. There is also a small offset between the linear and parabolic bands. The values of the gaps at the K point and energy offsets for LDA and GW are reported in Table 1, where they are also compared with recent experimental data. The agreement is good, specially for the energy offset. The GW gaps are systematically larger than the corresponding LDA gaps, which is a common trend observed in GW calculations Hybertsen and Louie (1986), and they are also larger than the experimental values. Possible reasons for this discrepancy are the unavoidable residual doping and substrate effects present in experiments, which tend to enhance screening and decrease the quasiparticle gap. We also point out that, although the LDA gaps seem to agree well with the experimental values, such an agreement is only fortuitous. It is well known that DFT underestimates quasiparticle gaps, even if the exchange-correlation functional was known exactly Martin (2004). Another important effect of the quasiparticle corrections is the renormalization of the Fermi velocity, which is visible from the increase of slope of the linear bands in GW when compared to LDA, a feature also observed in previous GW calculations in monolayer graphene Park et al. (2009); Yang et al. (2009); Trevisanutto et al. (2008).

For ABC stacking (right), the bandstructure is very different from the previous case. Near the Fermi level, there is a pair of bands with cubic dispersion, which are very flat near the K point. Moreover, in contrast with ABA, the ABC trilayer does have inversion symmetry, so the valence and conduction bands touch at the Fermi level. Due to the presence of a gap at the K point, these bands touch at three equivalent points located along the lines, where the dispersion is roughly linear. Further away from the Fermi level, there are two pairs of parabolic energy bands.

Figure 2: Comparison between the LDA (black lines) and GW (red lines) bandstructures obtained from our calculations for ABA (left) and ABC (right) stacking. The Brillouin Zone path is along the M-K- directions from left to right and it is centered in K. The Fermi level is set to zero in all cases, both in LDA and GW.
LDA GW Exp. Taychatanapat et al. (2011)
Monolayer gap 12 35 7
Bilayer gap 17 32 14
Offset 22 21 25
Table 1: Energy gaps and energy offset (meV) of the monolayer and bilayer-like bands in the ABA-stacked trilayer. The offset is defined as the energy difference between the middle of the two gaps.

Next, we fit the calculated ab initio bands to a low energy tight-binding (TB) model in order to extract the parameters that best describe the quasiparticle bandstructures. The hopping parameters included in our model are shown in Fig. 1, where the red dashed lines indicate the atoms connected by them. Atoms belonging to the inequivalent A and B sublattices are represented by yellow and dark spheres, respectively. Following common notation, is the nearest neighbor intralayer hopping and is the nearest neighbor interlayer hopping, connecting atoms that are right on top of each other in adjacent layers. These two parameters are sufficient to describe the main differences observed in the bandstructures, such as the type of dispersion of the low energy bands and their separations. The other parameters, also shown in Fig. 1, describe finer details and follow from a generalization of the Slonczewiski-Weiss-McClure (SWM) model of bulk graphite Dresselhaus and Dresselhaus (1981) and we use the same definition and sign convention for the parameters that are common to this model, which is extensively discussed in the literature Guinea et al. (2006); Partoens and Peeters (2006); Taychatanapat et al. (2011); Koshino and McCann (2009); Zhang et al. (2010). Another parameter, , not shown in Fig. 1, corresponds to the onsite energy difference between non-equivalent carbon atoms: High-energy sites correspond to carbon atoms on top of each other (connected by the hopping) and low-energy sites correspond to carbon atoms on top of hexagon centers in adjacent layers. Finally, we point out that surface effects are neglected in our TB calculations, since they were found to be very small in previous calculations on the ABA trilayer Koshino and McCann (2009).

We employ the least-squares method to fit the TB bands to either LDA or GW bands from our calculations, thus obtaining the set of parameters that best describes them. Since we want to describe the details of the bands near the K point (such as the small energy gaps), only k-points within a radius of from the K point are included in our fits. The comparison between GW bands and TB bands is shown in Fig. 3. We can see that the TB model describes very well the features observed in our calculations, even outside the radius of the fit. The fitted parameters are reported in Table 2. Different sign conventions are used in the literature, so we explicitly indicate when a different convention is being used. For both stackings, the quasiparticle corrections increase the value of , renormalizing the Fermi velocity () by about in ABA and in ABC trilayer, with respect to the LDA values. The GW value is in agreement with previous GW calculations on monolayer graphene Park et al. (2009); Yang et al. (2009); Trevisanutto et al. (2008) and with the experimental value for the monolayer et al. (2005). The parameter , which is associated with the distance between the higher energy bands and the Fermi level, is also increased by the quasiparticle corrections. In the absence of electron-hole asymmetry, these bands have energies for ABA and for ABC stacking at the K point Koshino and McCann (2009).

Figure 3: Comparison between the quasiparticle bands (black lines) and the corresponding adjusted TB bands (red lines) for ABA (left) and ABC (right) stacking. The path is the same as in Fig. 2 and the Fermi level is set to zero in both cases. Fitted parameters are shown in Table 2.
\arraybackslash \arraybackslashABA \arraybackslash \arraybackslash \arraybackslashABC \arraybackslash \arraybackslashBilayer \arraybackslashGraphite
Parameter \arraybackslashLDA \arraybackslashGW \arraybackslashExp. Taychatanapat et al. (2011) \arraybackslashLDA \arraybackslashGW \arraybackslashTheor. Zhang et al. (2010) \arraybackslashExp. Zhang et al. (2008) \arraybackslashExp. Dresselhaus and Dresselhaus (1981)
\arraybackslash2.590 \arraybackslash3.306 \arraybackslash3.100 \arraybackslash2.577 \arraybackslash3.188 \arraybackslash3.160 \arraybackslash3.000 \arraybackslash3.160
\arraybackslash0.348 \arraybackslash0.414 \arraybackslash0.390 \arraybackslash0.348 \arraybackslash0.415 \arraybackslash0.502 \arraybackslash0.400 \arraybackslash0.390
\arraybackslash-0.043 \arraybackslash-0.060 \arraybackslash-0.028 \arraybackslash-0.024 \arraybackslash-0.041 \arraybackslash-0.017 \arraybackslash- \arraybackslash-0.020
\arraybackslash0.283 \arraybackslash0.242 \arraybackslash0.315 \arraybackslash0.290 \arraybackslash0.323 \arraybackslash0.377 \arraybackslash0.300 \arraybackslash0.315
\arraybackslash0.162 \arraybackslash0.152 \arraybackslash0.041 \arraybackslash0.196 \arraybackslash0.287 \arraybackslash0.099 \arraybackslash0.150 \arraybackslash0.044
\arraybackslash0.024 \arraybackslash0.052 \arraybackslash0.050 \arraybackslash0.019 \arraybackslash0.126 \arraybackslash- \arraybackslash- \arraybackslash0.038
\arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash0.004 \arraybackslash0.084 \arraybackslash- \arraybackslash- \arraybackslash-
\arraybackslash0.010 \arraybackslash0.012 \arraybackslash0.046 \arraybackslash0.001 \arraybackslash0.023 \arraybackslash0.001 \arraybackslash0.018 \arraybackslash0.050
\arraybackslash0.002 \arraybackslash0.006 \arraybackslash- \arraybackslash0.002 \arraybackslash0.006 \arraybackslash- \arraybackslash- \arraybackslash-
This parameter was not obtained in this reference, it was set to a value from the literature.
The sign of this parameter was changed from the original reference in order to match our convention.
should not be confused with eV, a similar parameter used in graphite.
Table 2: Tight-binding parameters (in eV) obtained from our calculations (LDA or GW columns) and comparison with values from the literature, both theoretical and experimental. Wherever no value is shown, the corresponding parameter is either not applicable to that case or is set to zero. (in eV) is the standard deviation (error bars) of the adjustment inside its range.

We now discuss in detail the effects of the remaining TB parameters to the bandstructure features. The parameter causes a trigonal distortion of the low energy bands. As can be seen in Fig. 3, in ABA stacking, the parabolic bands have four energy minima: one at the K point and three along the equivalent lines. In our model, this is consistent with a positive sign for 3, in agreement with bulk graphite and bilayer graphene. In ABC, has a similar effect and is also positive, but there is an additional distortion caused by , as we discuss below. The parameter is responsible for a small electron-hole asymmetry in the bands. The ABA values for this parameter agree better with the experimental value from Bernal bilayer graphene than with the experimental value for this trilayer Taychatanapat et al. (2011); Zhang et al. (2008) (Table 2). On the other hand, the ABC values appear to be overestimated, since the TB bands show a larger asymmetry than predicted by LDA and GW, specially outside the range of the fit. This could be related to an insensitivity of the fit to this parameter, that is, since changes the curvature of the bands and they are mostly flat in the range of the adjustment, the curvature change may be not being correctly reproduced outside that range.

The remaining parameters have quite different roles for each trilayer. For ABA stacking, the parameters , and describe the inequivalency between A and B sublattices and are responsible for the opening of small band gaps in the linear and parabolic bands, also introducing an offset between them. As can be seen in Fig. 3, the TB model reproduces very well these features and the adjusted parameters are in good agreement with values from the literature. For ABC stacking, the parameter also opens a gap between the cubic bands at the K point, but the presence of inversion symmetry prevents the full opening of a band gap in this case. Hence, also induces a trigonal distortion, shifting the touching points of these bands from the K point to three points at the equivalent lines. This is consistent with a negative sign for in our model, in agreement with ABA trilayer and graphite. On the other hand, the parameters and don’t have any visible effects on the ABC bandstructure in the range considered, so they can be safely discarded in the description of the low energy bands, as was done previously Zhang et al. (2010). Nevertheless, we include in Table 2 the values given by the adjustment for future reference, but we stress that they could be also being affected by insensitivity. Note that the ABC value for cannot be directly compared with the values for ABA and graphite, since they have different definitions and roles. The parameter also plays a small role in ABC trilayer, introducing only a small electron-hole asymmetry in a similar fashion to . This could explain why the LDA value is smaller than the corresponding value for ABA, in agreement with the value obtained in a previous DFT calculation Zhang et al. (2010). On the other hand, the GW value is larger than the corresponding ABA value, which could indicate either a larger asymmetry of the bands or again an insensitivity of the adjustment.

In conclusion, we have studied the quasiparticle band structure of ABA and ABC trilayer graphene within the GW approximation. The bands obtained from these calculations were fitted to a TB model and the corresponding parameters were extracted. This model is found to reproduce very well the observed features in all cases, such as the type of dispersions and energy gaps. The fit parameters show a good agreement with available data from graphene bilayers, trilayers and graphite, and finer details of the bands are also correctly reproduced by them in the range of the fitting. The main effects of the quasiparticle corrections are the renormalization of the Fermi velocity and the increase of the separation between the higher energy bands, where both bring theory to closer agreement with experiment. Therefore, we expect this work to provide future reference for studies on graphene trilayers and other stackings, which can reveal more interesting properties and applications.

We thank Felipe Jornada for valuable discussions. This work was supported by the Brazilian funding agencies: CAPES, CNPq, FAPERJ and INCT-Nanomateriais de Carbono. Steven G. Louie acknowledges support from a Simons Foundation Fellowship in Theoretical Physics, and support from National Science Foundation Grant No. DMR10-1006184 (DFT and tight-binding analysis) and the Director, Office of Science, Office Basic Energy Sciences, Materials Sciences and Engineering Division, U.S. Department of Energy under contract No. DE-AC02-05CH11231 (GW simulations). We also thank NERSC for the computational resources employed in our calculations.


  1. preprint:
  2. This one-step procedure is known as a “” or a “one-shot” GW calculation. For more details, see Ref. Hybertsen and Louie (1986)
  3. Note that it’s common to find an opposite sign convention for in the literature Zhang et al. (2010); Koshino and McCann (2009); McCann and Koshino (2012); Taychatanapat et al. (2011), so care must be taken when comparing results and using reference values. A discussion of the conventions and their relation to the SWM model can be found in Ref. Partoens and Peeters (2006).


  1. K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Science 305, 666 (2004).
  2. A. H. C. Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009).
  3. T. Ohta, A. Bostwick, T. Seyller, K. Horn, and E. Rotenberg, Science 313, 951 (2006).
  4. Y. Zhang, T.-T. Tang, C. Girit, Z. Hao, M. C. Martin, A. Zettl, M. F. Crommie, Y. R. Shen, and F. Wang, Nature Lett. 459, 820 (2009).
  5. L. M. Zhang, Z. Q. Li, D. N. Basov, M. M. Fogler, Z. Hao, and M. C. Martin, Phys. Rev. B 78, 235408 (2008).
  6. E. V. Castro, K. S. Novoselov, S. Morozov, N. M. R. Peres, J. M. B. L. dos Santos, J. Nilsson, F. Guinea, A. K. Geim, and A. H. C. Neto, Phys. Rev. Lett. 99, 216802 (2007).
  7. F. Guinea, A. H. C. Neto, and N. M. R. Peres, Phys. Rev. B 73, 245426 (2006).
  8. B. Partoens and F. M. Peeters, Phys. Rev. B 74, 075404 (2006).
  9. T. Taychatanapat, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Nature Phys. 7, 621 (2011).
  10. M. Koshino and E. McCann, Phys. Rev. B 79, 125443 (2009).
  11. B.-R. Wu, Appl. Phys. Lett. 98, 263107 (2011).
  12. R. Olsen, R. van Gelderen, and C. M. Smith, Phys. Rev. B 87, 115414 (2013).
  13. M. M. Scherer, S. Uebelacker, D. D. Scherer, and C. Honerkamp, Phys. Rev. B 86, 155415 (2012).
  14. N. B. Kopnin, . T. T. Heikkilä, and G. E. Volovik, Phys. Rev. B 83, 220503(R) (2011).
  15. W. A. Muñoz, L. Covaci, and F. M. Peeters, Phys. Rev. B 87, 134509 (2013).
  16. F. Zhang, B. Sahu, H. Min, and A. H. MacDonald, Phys. Rev. B 82, 035409 (2010).
  17. C. H. Lui, Z. Li, K. F. Mak, E. Cappelluti, and T. F. Heinz, Nature Phys. 7, 944 (2011).
  18. L. Zhang, Y. Zhang, J. Camacho, M. Khodas, and I. A. Zaliznyak, Nature Phys. 7, 953 (2011).
  19. B. P. Neal, E. R. Ylvisaker, and W. E. Pickett, Phys. Rev. B 84, 085133 (2011).
  20. A. Alam and D. D. Johnson, Phys. Rev. Lett. 107, 206401 (2011).
  21. P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
  22. W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
  23. M. S. Hybertsen and S. G. Louie, Phys. Rev. B 34, 5390 (1986).
  24. J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981).
  25. N. Troullier and J. L. Martins, Phys. Rev. B 43, 1993 (1991).
  26. P. G. et al, J.Phys.:Condens.Matter 21, 395502 (2009).
  27. J. Deslippe, G. Samsonidze, D. A. Strubbe, M. Jain, M. L. Cohen, and S. G. Louie, Comput. Phys. Commun. 183, 1269 (2012).
  28. H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976).
  29. A. A. Mostofi, J. R. Yates, Y.-S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, Comput. Phys. Commun. 178, 685 (2008).
  30. R. M. Martin, Electronic Structure: Basic Theory and Practical Methods (Cambridge University Press, 2004).
  31. C.-H. Park, F. Giustino, C. D. Spataru, M. L. Cohen, and S. G. Louie, Nano Lett. 9, 4234 (2009).
  32. L. Yang, J. Deslippe, C.-H. Park, M. L. Cohen, and S. G. Louie, Phys. Rev. Lett. 103, 186802 (2009).
  33. P. E. Trevisanutto, C. Giorgetti, L. Reining, M. Ladisa, and V. Olevano, Phys. Rev. Lett. 101, 226405 (2008).
  34. M. S. Dresselhaus and G. Dresselhaus, Advances in Physics 30, 139 (1981).
  35. Y. Z. et al., Nature (London) 438, 201 (2005).
  36. E. McCann and M. Koshino, cond-mat.mes-hall (2012).
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