Twisted bilayer graphene from first-principles: structural and electronic properties

TL;DR

This work addresses obtaining atomistically accurate structural and electronic properties of twisted bilayer graphene across a wide range of twist angles using density functional theory with an optimized local basis. By validating against plane-wave DFT and exact $\mathbf{k}\cdot\mathbf{p}$ models, it provides fully relaxed commensurate structures down to $0.987^{\circ}$ and demonstrates that lattice relaxation agrees with continuum elasticity, while electronic bands follow $\mathbf{k}\cdot\mathbf{p}$ trends albeit with a twist-angle offset around $\Delta\theta \approx 0.05^{\circ}$. The study also reveals moiré-scale wavefunction character and symmetry properties, including same chirality at the Dirac nodes, establishing an ab initio reference for future many-body analyses in tBLG. Overall, the approach enables detailed, accurate simulations of large moiré systems and lays the groundwork for incorporating correlation effects in a first-principles framework.

Abstract

We present a comprehensive first-principles study of twisted bilayer graphene (tBLG) for a wide range of twist angles, with a focus on structural and electronic properties. By employing density functional theory (DFT) with an optimized local basis set, we simulate tBLG, obtaining fully relaxed commensurate structures for twist angles down to 0.987°. For all angles the lattice relaxation agrees well with continuum elastic models. For angles accessible to plane-wave DFT (VASP), we provide a detailed comparison with our local basis DFT (SIESTA) calculations, demonstrating excellent agreement in both the atomic and electronic structure. The dependence of the Fermi velocity and band width on the twist angle shows qualitative agreement with results from an `exact' $\mathbf{k \cdot p}$ continuum model, but reveals a small twist angle offset. Additionally, we provide details of the low-energy wavefunction character, band inversion and symmetries. Our results provide an ab initio reference point for the microscopic structure and electronic properties of tBLG which will serve as the foundation for future studies incorporating many-body effects.

Figures (9)

• Figure 1: Comparison of the atomic structure and forces from DFT calculations using siesta and vasp, for tBLG with a twist angle of $\theta = 2.45^{\circ}$ and a monolayer lattice constant of $a_0 = 2.46$ Å. (a)-(b) The difference in positions obtained from siesta and vasp. (a) the in-plane direction ($\Delta\mathbf{r}$). The in-plane displacements due to relaxation are on the order of $10^{-2}$ Å for both codes, and the difference $\Delta\mathbf{r}$ is less than $1\%$ of the relaxation. (b) the out-of-plane direction ($\Delta z$), projected onto the unit cell diagonal, where the top layer is shown in blue and the bottom layer in red. (c)-(d) The forces obtained from vasp using the relaxed coordinates determined by siesta (red/blue) and vasp (green/purple) are shown for: (c) the in-plane direction ($\mathbf{F}_{\parallel}$); (d) the out-of-plane direction ($F_{\perp}$), projected onto the unit cell diagonal. The horizontal black lines show the force tolerance used in the vasp relaxation (10 meV/Å). (e) Comparison of the band structure of tBLG obtained using siesta (black lines) with the SZP basis set, and vasp (red, dashed lines).
• Figure 2: (a)-(b) Lattice relaxation in tBLG from DFT, for $\theta = 0.987^{\circ}$: (a) the in-plane displacement $\Delta\mathbf{r}$ of the top layer, obtained as the difference in positions between the rigid and relaxed structures, and (b) layer separation $d$. (c)-(d) GSFE and interlayer distance calculated using several different basis sets in siesta: (c) the stacking energy, and (d) the interlayer separation were calculated for untwisted graphene bilayers by sliding one layer over the other along the unit cell diagonal, using the SZP basis (red) and the DZPF basis (blue). The dashed black line shows the GSFE obtained using the fitting parameters in Ref. carr2018relaxation (Carr et al.). The circles represent results from DFT calculations, and the lines show fits to the data. (e)-(f) Elastic properties of monolayer graphene obtained from derivatives of total-energy curves: (e) energy versus lattice constants changes used to calculate the bulk modulus; (f) energy versus shear angle used to calculate the shear modulus. For the bulk modulus, the lattice constant $a$ was perturbed about its equilibrium value $a_0$. For the shear modulus, the angle $\gamma$ between the in-plane lattice vectors was modified about its equilibrium value $\gamma_0 = 60^{\circ}$.
• Figure 3: The absolute value $\left | \mathbf{u}_1 \right|$ of the first Fourier component of the displacement field (Eq. \ref{['eq:u-displacement']}) as a function of twist angle. The values obtained from lattice relaxation calculations are shown in blue, which were obtained by numerically minimizing Eq. \ref{['eq:V-tot']} using the SZP parameters in Table \ref{['tab:GSFE']}. The red circles show the corresponding Fourier components of the relaxed DFT structure, obtained by fitting the in-plane displacement to Eq. \ref{['eq:u-displacement']} for each twist angle.
• Figure 4: Electronic band structure of tBLG for commensurate twist angles near the magic angle. The black lines show the bands from DFT calculations with siesta and the red lines show the bands obtained using the $\mathbf{k}\cdot\mathbf{p}$ model. A small shift of the twist angle used in the $\mathbf{k}\cdot\mathbf{p}$ calculation, $\theta_{\mathbf{k}\cdot\mathbf{p}} = \theta_{\rm DFT} + \Delta\theta$ where $\Delta\theta \approx 0.05^{\circ}$, improves the agreement with the DFT bands.
• Figure 5: (a) Fermi velocity as a function of twist angle from siesta (black) and the exact $\mathbf{k}\cdot\mathbf{p}$ model (red), as a percentage of the Fermi velocity of graphene, $v^0_{\rm F} = 1\times 10^{6}$ m/s. (b) Band width $E_{w}$ (meV) at $\Gamma$ as a function of twist angle from siesta (black) and the exact $\mathbf{k}\cdot\mathbf{p}$ model (red). (c) Band gaps $E_{g}$ (meV) between the flat and dispersive bands at $\Gamma$ as a function of twist angle from siesta (black) and the exact $\mathbf{k}\cdot\mathbf{p}$ model (red). The upper gap is represented by the solid circles and the lower gap is represented by the hollow circles. In this figure $\theta_{\mathbf{k}\cdot\mathbf{p}}$ has not been shifted.