Hybrid-functional electronic structure of multilayer graphene

Multilayer graphene with rhombohedral and Bernal stacking is supposed to be metallic, as predicted by density functional theory calculations using semilocal functionals. However, recent angular resolved photoemission and transport data have questioned this point of view. In particular, rhombohedral flakes are suggested to be magnetic insulators, a view supported also by hybrid-functional calculations. Bernal flakes composed of an even number of layers are insulating (for N <= 6), while those composed of an odd number of layers are pseudogapped (for N <= 7). Here, by systematically benchmarking with plane-waves codes, we develop very accurate all-electron Gaussian basis sets for graphene multilayers, allowing a precise description of the electronic structure in the 100 meV energy range from the Fermi energy at the hybrid-functional level. We find, in agreement with our previous calculations, that rhombohedral stacked multilayers are gapped and magnetic. However, the valence band curvature and the details of the electronic structure at the similar to 10 meV scale show a dependence on the basis set. A substantially extended basis set is needed to describe the long-range interlayer interactions and, consequently, to correctly reproduce the effective mass of the valence band top at the K point. In the case of Bernal stacking, we show that exact exchange gaps the flakes composed by four layers and opens pseudogaps for N = 3, 6, 7, 8. However, the gap or pseudogap size and its behavior as a function of thickness are not compatible with experimental data. Moreover, hybrid functionals lead to a metallic solution for five layers and a magnetic ground state for five, six, and eight layers. Magnetism is very weak with practically no effect on the electronic structure and the magnetic moments are mostly concentrated in the central layers. Our hybrid-functional calculations on trilayer Bernal graphene are in excellent agreement with GW results. For thicker multilayers, our calculations are a benchmark for many-body theoretical modeling of the low energy electronic structure.