Band Unfolding Made Simple

TL;DR

This paper introduces a straightforward method for band unfolding in supercell calculations, utilizing the local density of states in reciprocal space, applicable to both periodic and non-periodic systems, and demonstrates its implementation in the SIESTA package with applications to Si and graphene defects.

Contribution

It presents a simple, intuitive approach to band unfolding based on qLDOS, valid for systems with or without translational symmetry, and provides an implementation in SIESTA.

Findings

Effective band unfolding for defective systems demonstrated

Method applicable to both periodic and non-periodic structures

Implementation successfully applied to Si and graphene defects

Abstract

We present a simple view on band unfolding of the energy bands obtained from supercell calculations. It relies on the relationship between the local density of states in reciprocal space (qLDOS) and the fully unfolded band structure. This provides an intuitive and valid approach not only for periodic, but also for systems with no translational symmetry. By refolding into the primitive Brillouin zone of the pristine crystal we recover the conventional unfolded bands. We implement our algorithm in the SIESTA package and apply it to defects on Si and graphene.

Equations

$$n_{SBZ}(\mathbf{K},ϵ)=\sum _i\delta (ϵ-{ϵ}_{\mathbf{K},i})$$

$$n(\mathbf{q},ϵ)$$

$${\psi }_{\mathbf{K},i}(\mathbf{r})=\frac{1}{(2\pi {)}^{3/2}}{\int }_{\infty }d\mathbf{q}{\tilde{\psi }}_{\mathbf{K},i}(\mathbf{q})e^{i\mathbf{qr}},$$

$${\tilde{\psi }}_{\mathbf{K},i}(\mathbf{q})=(2\pi {)}^{3/2}\sum _{\mathbf{G}}\delta (\mathbf{K}+\mathbf{G}-\mathbf{q}){\tilde{u}}_{\mathbf{K},i,\mathbf{G}},$$

$${\tilde{u}}_{\mathbf{K},i,\mathbf{G}}=\frac{1}{V_{SC}}{\int }_{SC}d\mathbf{r}{u}_{\mathbf{K},i}(\mathbf{r})e^{-i\mathbf{Gr}},$$

$$n(\mathbf{r},ϵ)=\frac{V_{SC}}{(2\pi {)}^3}\sum _i{\int }_{SBZ}d\mathbf{K}|{\psi }_{\mathbf{K},i}(\mathbf{r}){|}^2\delta (ϵ-{ϵ}_{\mathbf{K},i}).$$

$$n(ϵ)=\frac{1}{V_{SC}}{\int }_{SC}d\mathbf{r}n(\mathbf{r},ϵ)=\frac{1}{(2\pi {)}^3}{\int }_{\infty }d\mathbf{q}n(\mathbf{q},ϵ).$$

$$n_{RBZ}(\mathbf{k},ϵ)=\sum _{\mathbf{g}}n(\mathbf{k}+\mathbf{g},ϵ),$$

$${\psi }_{\mathbf{K},i}(\mathbf{r})=\sum _{\mathbf{R}}\sum _{\mu }{c}_{\mathbf{K},i,\mu }{\varphi }_{\mu }(\mathbf{r}-\mathbf{R}-{\mathbf{r}}_{\mu })e^{i\mathbf{K}(\mathbf{R}+{\mathbf{r}}_{\mu })},$$

$${\tilde{u}}_{\mathbf{K},i,\mathbf{G}}=\frac{(2\pi {)}^{3/2}}{V_{SC}}\sum _{\mu }{c}_{\mathbf{K},i,\mu }{\tilde{\varphi }}_{\mu }(\mathbf{K}+\mathbf{G})e^{-i\mathbf{G}{\mathbf{r}}_{\mu }},$$

$${\varphi }_{\mu }(\mathbf{r})={\varphi }_{\mu }(r)Y_{l_{\mu },m_{\mu }}(\hat{\mathbf{r}})$$

$${\tilde{\varphi }}_{\mu }(\mathbf{q})={\tilde{\varphi }}_{\mu }(q)Y_{l_{\mu },m_{\mu }}(\hat{\mathbf{q}})$$

$${\tilde{\varphi }}_{\mu }(q)=\sqrt{\frac{2}{\pi }}(-1{)}^{l_{\mu }}{\int }_0^{\infty }{r}^{2}dr{j}_{l_{\mu }}(qr){\varphi }_{\mu }(r),$$

Full text

I Introduction

Plots of so called energy bands are the most basic and used tool in interpreting the calculated electronic structure of simple crystals. Such plots represent the energy of the Bloch orbitals as a function of their crystal momentum in the primitive Brillouin zone (PBZ), that is, $E\left(\mathbf{k}\right)$. These theoretical band structures can be obtained within the tight-binding approximation or the density functional theory (DFT) and they have a direct connection with the results of angle resolved photoemission spectroscopy (ARPES) experiments. However, the simplicity of this approach disappears when the calculations involve large supercells with many atoms. As the size of the cell in real space increases, the first Brillouin zone in reciprocal space shrinks and more lines populate the band structure, hindering the extraction of useful information and the comparison with experiments.

Several authors have already developed techniques to unfold the supercell Brillouin zone into the primitive Brillouin zone. Some of the existing works focus on algorithms within the tight-binding approximation Boykin2005 ; Boykin2007 ; Boykin2018 ; Dargam or first principles calculations, employing as basis sets linear combinations of atomic orbitals (LCAO) Lee , plane waves Chen2018 or Wannier functions Ku . Another methodology studies the electronic structures of alloy systems Popescu10 ; Popescu12 . Other authors have delved into the theory of the problem and they have developed general formulations of it Allen2013 ; Huang ; Medeiros ; Rubel ; Kosugi .

Most of these methods search a direct relationship between Brillouin zones of the primitive crystal cell and of the simulation supercell, by expanding the supercell states in a basis set with the periodicity of the crystal. In contrast, we aim to simplify the procedure and provide a link between existing methods by dividing the unfolding problem in two steps. First, we consider the fully unfolded bands, extended to the whole reciprocal space, through a Fourier decomposition of the Bloch wave functions of the system Kosugi . This results in a non-periodic description, interesting by itself. Second, by what we call a refolding into the crystal primitive Brillouin zone, we recover the conventional unfolded bands of other authors.

II Method

The basic steps of our unfolding method are schematically summarized in Fig. 1. We start by considering the energy bands not just as dispersion relations, but as the density of states in the (first) simulation Brillouin zone (SBZ), the reciprocal of the real-space simulation (super) cell (SC):

[TABLE]

where $\mathbf{K}$ is a wave vector in the SBZ, and $\epsilon_{\mathbf{K},i}$ is the eigenvalue of its $i$th Bloch state, that is, band energy notation .

The normalization of $n_{SBZ}(\mathbf{K},\epsilon)$, as well as those of other densities defined below, is such that they are densities of states per unit of macroscopic volume (as well as per unit of their respective independent variables), what facilitates the comparison between simulation cells of different volumes.

Next, we split the normalized weight $\delta(\epsilon-\epsilon_{\mathbf{K},i})$ according to the squared Fourier coefficients of the corresponding wave function. Summing over Bloch states, we construct the fully unfolded bands as the spectral density, which can also be considered as the local density of states in reciprocal space (qLDOS):

[TABLE]

where $\psi_{\mathbf{K},i}(\mathbf{r})$ is a Bloch wave function (normalized in the SC), $u_{\mathbf{K},i}(\mathbf{r})\equiv\psi_{\mathbf{K},i}(\mathbf{r})e^{-i\mathbf{K}\mathbf{r}}$ is its periodic part, and $\tilde{\psi}_{\mathbf{K},i}(\mathbf{q}),\tilde{u}_{\mathbf{K},i,\mathbf{G}}$ are their respective Fourier transforms:

[TABLE]

[TABLE]

[TABLE]

with $V_{SC}$ the volume of the SC and $\mathbf{G}$ its reciprocal wave vectors. In Eq. (2), $\mathbf{K}_{\mathbf{q}}$ and $\mathbf{G}_{\mathbf{q}}$ are the unique vectors such that: $\mathbf{K}_{\mathbf{q}}$ is within the SBZ; $\mathbf{G}_{\mathbf{q}}$ is a reciprocal wave vector; and $\mathbf{K}_{\mathbf{q}}+\mathbf{G}_{\mathbf{q}}=\mathbf{q}$. A state $\psi_{\mathbf{K},i}$ will contribute to $n(\mathbf{q},\epsilon)$ at points $\mathbf{q}=\mathbf{K}+\mathbf{G}$ for all $\mathbf{G}=\pi N/a$ vectors, due to Bloch’s theorem. We emphasize that $\mathbf{q}$ extends to infinity and $n(\mathbf{q},\epsilon)$ is not periodic in $\mathbf{q}$: although the energies at which $n(\mathbf{q},\epsilon)$ can be nonzero are periodic, these “bands” have a different weight at each Brillouin zone (Fig. 1(c)).

As can be seen by comparing with Eq. (2), the qLDOS is the Fourier-space equivalent of the real-space local density of states (rLDOS),

[TABLE]

The total density of states (DOS) can be obtained by integration of either $n(\mathbf{r},\epsilon)$ or $n(\mathbf{q},\epsilon)$:

[TABLE]

Since $|\tilde{\psi}_{\mathbf{K},i}(\mathbf{q})|^{2}$ is the probability of measuring momentum $\mathbf{q}$ of a given electron, $n(\mathbf{q},\epsilon)$ is the probability of finding an electron (or an empty state) in the system with energy $\epsilon$ and momentum $\mathbf{q}$, and it can thus be directly related with ARPES results if matrix element effects are taken into account Kosugi ; surface ; blackphos .

The qLDOS, that we call fully unfolded bands is the same as the spectral weight of other references Allen2013 ; Medeiros ; Rubel and the plane-wave unfolded spectra introduced by Kosugi et al. Kosugi . Therefore, our approach is a different description, rather than a new method that yields different results. Our emphasis is to provide a clear and simple link with previous methods through the (L)DOS, as well as to generalise band structure analysis to non periodic systems.

The last step in our method is to refold $n(\mathbf{q},\epsilon)$ into a refolding Brillouin zone (RBZ) as

[TABLE]

where $\mathbf{k}$ is within the RBZ and $\mathbf{g}=\pi n/a$ are its reciprocal lattice vectors. Notice that, since $\sum_{\mathbf{G}}|\tilde{u}_{\mathbf{K},i,\mathbf{G}}|^{2}=1$, then   $\sum_{\mathbf{G}}n(\mathbf{K}+\mathbf{G},\epsilon)=n_{SBZ}(\mathbf{K},\epsilon)$, that is, refolding $n(\mathbf{q},\epsilon)$ back into the SBZ recovers the original bands.

Frequently, the simulation cell will be a supercell of the refolding cell. In these cases, the RBZ will be a supercell of the SBZ, and vectors $\mathbf{g}$ will belong to the set of $\mathbf{G}$s (Fig. 1(d)). Nevertheless, this condition is not necessary in our method, and in fact it will not be true in many cases, as for simulation cells of liquids or amorphous systems, or of defects that induce strong deformations cost .

The above full unfolding/refolding method can be immediately applied in a plane wave DFT code Abinit ; QE ; Vasp , since the Fourier coefficients of the Bloch wave functions are then directly available. The slow decay with momentum of all-electron wavefunctions can be addressed by using pseudopotentials or by introducing a momentum cutoff. For a basis of atomic orbitals, we expand the Bloch states as

[TABLE]

where $c_{\mathbf{K},i,\mu}$ are expansion coefficients and $\phi_{\mu}$ are atomic orbitals centered at position $\mathbf{R}+\mathbf{r}_{\mu}$ ($\mathbf{R}$ being SC lattice vectors). Substituting into Eq. (5) we find

[TABLE]

where $\tilde{\phi}_{\mu}(\mathbf{q})$ is the Fourier transform of the numerical atomic orbital $\phi_{\mu}(\mathbf{r})$, with well defined angular momentum quantum numbers $(l_{\mu},m_{\mu})$, that can be decomposed into radial and angular parts:

[TABLE]

[TABLE]

[TABLE]

with $Y_{l,m}(\hat{\mathbf{r}})$ spherical harmonics and $j_{l}(x)$ spherical Bessel functions.

This algorithm has been implemented in the Siesta package Soler2002 . After a converged DFT Siesta calculation, the hamiltonian and overlap matrices, in the atomic basis set, are calculated and written in a file. This file, as well as those specifying the radial numerical atomic orbitals, are read by an external utility program that calculates the fully unfolded and refolded spectra at the desired $\mathbf{q}$ and $\mathbf{k}$ band lines. Some Siesta subroutines are also used by the unfolding/refolding program to obtain the wave function coefficients at each required $\mathbf{K}$ point of the SBZ, as well as to perform the Fourier transforms in Eq. (13).

III Applications

We apply our previously described Siesta implementation to a Si FCC crystal with a single vacancy, a model of amorphous Si (a-Si), a monolayer of graphene with a (585) divacancy, and a rotated graphene bilayer under pressure. We employed the GGA-PBE PBE functional for exchange and correlation and double-$\zeta$ + polarization (DZP) basis sets (double-$\zeta$ (DZ) for a-Si).

III.1 Vacancy in Si crystal

We model a vacancy in a Si FCC crystal using a 63-atom supercell. Its refolded bands into the crystal PBZ vectors are depicted in Fig. 2(a), and can be compared with the bands of the periodic crystal, in Fig. 2(b). Changes are appreciated at a careful sight. The refolded bands become blurred and widened due to the appearance of small splittings. Some of these are a consequence of the supercell approach and they become smoother in a larger SC.

The most relevant changes occur around the Fermi level. A new state arises within the gap, with higher weight around the $\Gamma$ point. The top of the valence band, around $\Gamma$, decreases in energy and in weight. We show the difference between the crystal and defective refolded spectra in Fig. 2(d) around $\Gamma$, at the gap, to remark these changes.

III.2 Amorphous Si

A clear example of disorder is an amorphous solid. In this case, we cannot talk about a proper band structure, but yet the energy dispersion of the electron states provides interesting results. We studied the fully unfolded and refolded bands of a-Si using supercells of 216, 512 and 1024 atoms, modelled by Igram et al. Igram , obtaining equivalent results for all cases. We present the results for the 512-atom supercell.

Fig. 3(a) show the fully unfolded bands along two directions of reciprocal space, corresponding to the $L-\Gamma$ and $\Gamma-X$ directions of the crystal, up to the second Brillouin zone. As expected from the isotropic character of the amorphous solid, its unfolded bands are essentially identical in both directions, yielding a widened free-electron dispersion with an effective mass $m^{*}=1.1$ and a gap of $1.5$ eV, consistent with existing values Igram ; Barber .

The same bands, refolded to the PBZ of the crystal, are shown in Fig.3(b). At close magnification and inspection, a few localised states appear in the gap due to defects in the a-Si model Igram . Interestingly enough, despite its isotropic and non-periodic structure, and the incommensurability of its simulation cell with that of the crystal, the refolded bands of a-Si appear as a blurred version of the crystalline silicon (c-Si) bands, specially at low energies. Also, higher LDOS of a-Si can be appreciated in regions of c-Si band degeneracies. This similarity is independent of the size of the a-Si simulation cell employed. We attribute it to the similarity of the local structure in a-Si and c-Si in bond distances and angles Igram .

III.3 Divacancy in graphene

Graphene is a material with unique electronic properties, but highly sensitive to structural disorder. The presence of defects leads to significant changes on its bands, specially around the Fermi level. Many types of defects have been studied by theorists, such as adatoms, vacancies or Stone-Wales defects, with the aim of predicting their properties and, potentially, using them to tailor the functionalities of graphene. Hence, unfolding their band stuctures will shed light on how they modify the original graphene bands.

We consider a (585) divacancy in a graphene monolayer, a defect that has been synthesised and characterised by Ugeda et al. Ugeda . The (585) defect is formed by two adjacent C vacancies rearranged into two pentagons and one octagon, as shown in Fig. 4(a), with no dangling $\sigma$ bonds. We modelled the divacancy in simulation cells of different sizes, from $7\times 7$ to $14\times 14$. Here we report the $11\times 11$, $12\times 12$ and $13\times 13$ as representative.

The conventional bands of the supercell calculations within their SBZs, Fig. 5, reveal that the former six-fold symmetry of the Brillouin zone is broken, leading to inequivalent $K$ and $M$ points and two emerging mirror planes, Fig. 4(b), as happens in the atomic structure. The $K$ and $K^{\prime}$ points no longer present a Dirac cone, although all of them present a band crossing of $E_{F}$ at different points in their Brillouin zones: near $K^{\prime}$ for the $11\times 11$ SC, around $\Gamma$ for the $12\times 12$ SC and between $M^{\prime}$ and $K$ for the $13\times 13$ SC. Besides this, it is hard to find similarities between them. To relate these overcrowded spectra with the band structure of graphene, we compute the refolded bands into the PBZ in the surroundings of $K$ and $K^{\prime}$ points of graphene.

Fig. 6(a) and (c) show the $11\times 11$ and $13\times 13$ SCs bands refolded into the PBZ. The Fermi level crossing is located around $K^{\prime}$ of graphene, which corresponds to a $K^{\prime}$ point in both SBZs. One of the bands conforming the former Dirac cone stays almost unaltered, while the other one is split, forming two cone tips shifted in the $k_{x}$ direction and connected by a state. In the 2D bands plot of Fig. 7(a) and (c) this is seen with more clarity. The Fermi level decreases from its original value, being coincident with the tip of the lower cone.

In the case of the $12\times 12$ SC, Fig. 6(b), we see a similar behaviour, with slight differences. One band of the cone remains almost unaltered, with a small gap of $0.04$ eV opening at $K^{\prime}$. The other one is split, with one of its fragments conforming a flat state at the Fermi level, leading to a single band crossing. In Fig. 7(b) this flat state and the cone tips can be appreciated. Here, the $K^{\prime}$ of graphene coincides with a $\Gamma$ point of the $12\times 12$ SBZ.

Despite the differences between the three sizes, a general trend is clearly identified after refolding into the PZB of graphene, which is not the case in the traditional bands description of Fig. 5. All cases present a splitting and a shift in $k_{x}$ of the Dirac cone, as well as a single band between the two cone tips as the lowest unoccupied electronic state. The Fermi level is located at the tip of the lower cone. We remark that, in $3n\times 3n$ SCs, the $K^{\prime}$ of the PBZ corresponds to a $\Gamma$ point in the SBZ, unlike in the $(3n+1)\times(3n+1)$ and $(3n+2)\times(3n+2)$ cases. This is consistent with the formal differences existing between both kinds of supercells Lastra . Finally, the refolded bands obtained around $K$ are inverted in $k_{y}$ compared to the refolded bands around $K^{\prime}$ shown above.

III.4 Pressure in rotated graphene bilayers

Another path to modify the electronic behaviour of graphene is given by rotated graphene bilayers. These are defect-free systems with a Moiré pattern, involving large periodicities. At small angles, the interaction between the two layers induces two saddle points in the band structure, along with two logarithmic van Hove singularities in the DOS. As the angle decreases, the singularities approach and, eventually, at the so called magic angles ($\theta=1.1^{\circ}$), they collapse into a single peak at the Fermi level Cao . The same effect has been recently reported to occur for larger angles, when external pressure is applied Felix ; Bezanilla ; Carr ; Yankowitz .

We study a bilayer rotated $5.08^{\circ}$, at equilibrium under increasing pressures up to $1.63$ GPa. We employ a GGA exchange-correlation functional including Van der Waals interactions Berland ; Roman . As none of the monolayers lattice orientations has a prevalence over the other, refolding into the PBZ of one of them is not particularly informative. Therefore, we consider the fully unfolded bands as the adequate tool to analyze this system. We first show the fully unfolded bands of a graphene monolayer, Fig. 8, up its second Brillouin zone, for comparison purposes. Two paraboloids with gaps opening along them and a six-fold symmetry and are clearly distinguished, and conform the dispersion relations of the $\sigma$ and $\pi$ orbitals. In the case of the bilayer, we will restrict the unfolding region to the surroundings of a $K$ (and $K^{5.08^{\circ}}$) point (red square).

In Fig. 9 we show the conventional band structures of the bilayer at the equilibrium distance, at a middle stage and under a pressure of $1.63$ GPa, next to their corresponding DOS. The saddle points and van Hove singularities can be appreciated.

Fig. 10 depicts the evolution of the fully unfolded spectra under increasing pressures. At equilibrium configuration, Fig. 10(a) shows a neat picture of the interaction between the cones of both monolayers, as well as the saddle point emerging in between. The relative maximum intensities of the LDOS are homogeneous in energies, and the states present high dispersion in energies. As pressure is applied, Fig. 10(b),(c),(d) and (e), the cones flatten and the saddle points move towards the Fermi level. We appreciate as well that the two pairs of bands immediately over and below the cones lower their energies and start to merge. Also, in Fig. 10(d) and (e), higher weights correspond to these merging bands, whereas the cone states around the Fermi level tend to be delocalized in many $\mathbf{q}$s. In Fig. 10(f) the cones have collapsed into flat bands. We remark how their weight in this area of reciprocal space is small compared to that of the merging bands, despite the sharp peak on the DOS of Fig. 9(c). This is an indicator of delocalization in $\mathbf{q}$, and is not unexpected, as these states are known to be well localized in the AA stacking region in real space Ugeda ; Felix .

IV Conclusions

We have presented a simple formulation of the band unfolding problem, a tool necessary to extract useful information from the band structure of large supercell calculations. The idea of a full unfolding that expands the bands not only to the primitive cell, but to the full reciprocal space, allows to treat this problem as a decomposition of the wave functions into its Fourier coefficients. A refolding recovers the conventional unfolded bands in the PBZ of the crystal. It is feasible for any eigenstate, regardless of the basis used. In the case of plane wave codes this implementation shall be almost immediate.

We have successfully applied our algorithm to obtain new characterizations of non-trivial physical systems. The fully unfolded bands provide a distribution of the states as a function of their energy and momenta, allowing a direct comparison with experimental photoemission spectra, as well as a way to determine a value of the effective mass of the system under study in a chosen direction of reciprocal space. Refolding into the primitive cell yields clear band spectra that allow comparison with the crystal bands, even identifying crystal-like patterns in an amorphous solid.

The outcomes of this work prove that the underlying state distribution in reciprocal space is much richer than what conventional band structures can evince, transcending the existence of any real or imposed periodicity.

V Akcnowledgements

We are grateful to B. Battarai and D. Drabold for sharing the coordinates of their a-Si model with us. JMS thanks Prof. P. B. Allen for introducing him into the unfolding problem. Funded by Spain’s MINECO grant FIS2015-64886-C5-5-P.