Electronic state unfolding for plane waves: energy bands, Fermi surfaces and spectral functions

Modern computing facilities grant access to first-principles density-functional theory study of complex physical and chemical phenomena in materials, that require large supercell to properly model the system. However, supercells are associated to small Brillouin zones in the reciprocal space, leading to folded electronic eigenstates that make the analysis and interpretation extremely challenging. Various techniques have been proposed and developed in order to reconstruct the electronic band structures of super cells, unfolded into the reciprocal space of an ideal primitive cell. Here, we propose an efficient unfolding scheme embedded directly in the Vienna Ab-initio Simulation Package (VASP), that requires modest computational resources and allows for an automatized mapping from the reciprocal space of the supercell to primitive cell Brillouin zone. This algorithm can computes band structures, Fermi surfaces and spectral functions, by using an integrated post-processing tool (bands4vasp). The method is here applied to a selected variety of complex physical situations: the effect of doping on the band dispersion in the BaFe$_{\rm 2(1-x)}$Ru$_{\rm 2x}$As$_2$ superconductor, the interaction between adsorbates and polaronic states on the TiO$_2$(110) surface, and the band splitting induced by non-collinear spin fluctuations in EuCd$_2$As$_2$.


Introduction
Material science simulations adopting periodic boundary conditions in the framework of density functional theory (DFT) may require large unit cells in order to model long or broken periodicity in crystals. Supercells (i.e., large unit cells built by stacking of smaller primitive cells forming an ideal Bravais lattice), are used to study the effects of lattice impurities (e.g., local dislocations, defects, doping), but also to investigate domain boundaries, magnetic orders, surface reactivity, and structural reconstructions, just to name a few common applications. 1,2 While well developed facilities and efficient DFT packages are capable to deal with hundreds and even thousands of atoms in large cells, the analysis of the electronic properties (such as the energy band structure or the Fermi surface) gets complicated by the shrinking of the Brillouin zone (Bz) and the consequent folding of the eigenstates in the reciprocal space, 3 that prevent also a genuine comparison with photo-emission spectroscopy experiments. 4,5 Figure 1 shows an example of the intricate band structure typically obtained by using a supercell: Clearly, the bands calculated by using the primitive cell allow instead for a more Figure 1: Example of eigenstate folding for the pristine BaFe 2 As 2 compound. Band structures obtained by using a supercell (a) and a primitive cell (b). straightforward analysis. The intricate supercell states can be unfolded back into the larger Brillouin zone of the primitive cell by applying the unfolding technique. 3,[6][7][8][9][10] This technique is based on the projection P Km of the supercell eigenstates |Km on the primitive cell eigenstates |kn : where m and n denote energy band indices at vectors K and k i in the reciprocal space of the supercell and primitive cell, respectively. This projection represents the amount of Bloch character of the states |k i n contributing to |Km , which allows for a direct connection between the reciprocal space of the supercell and the primitive cell. By weighting the contributions of all single states with P Km , it is indeed possible to obtain an effective band structure (EBS) of the supercell unfolded in the larger Brillouin zone of the primitive cell. Equation 1 can be rewritten in terms of states of the supercell only, if the eigenstates |Km are expanded in terms of a plane-wave basis set with coefficients C m,K : By expressing P Km in this form, all information required from the primitive cell are purely geometric and collected by the reciprocal lattice vector g applied to the k i vectors in the reciprocal space. This alternative formulation brings the advantage to avoid any calculation on the primitive cell as well as any direct comparison between the two spaces, which could turn out to be technically challenging. The Vienna Ab-initio Simulation Package (VASP) is an optimal candidate for the implementation of unfolding calculations: This code can deal with large cells efficiently, and plane waves are used as basis functions. Moreover, a basic implementation of Eq. 2 is already available in the recent VASP releases. 11 In this work we extend the original unfolding scheme aiming to reduce the memory requirements, and to simplify the user interface for both input parameters and extraction of output data. Specifically, we implemented an automatized scheme for generating the supercell reciprocal-space vectors K starting from given k i vectors in the primitive cell space. The calculation of the P Km projection can be limited only to the automatically determined (K, k i ) pairs of interest, saving memory resources in the calculation. The user interface has been simplified, including an automatic initialization of the primitive cell; moreover, the user is provided with a post-processing package for the analysis of the results, "bands4vasp". 12 This updated implementation of the unfolding technique in VASP can be efficiently applied to a wide range of physical problems: in the following we show examples of such applications, including a benchmark on the superconductivity emerging on BaFe 2 As 2 upon Ru doping, the interplay between polarons and adsorbates on semiconducting TiO 2 (110) surface, and the non-collinear magnetic ordering on EuCd 2 As 2 .

Methodology
Electronic eigenstates calculated for supercells, can be unfolded into the reciprocal space of the primitive cell by using external packages based on different methods, [13][14][15][16][17] or directly in VASP that implements Eq. 2. 11,[18][19][20] We have optimized this unfolding algorithm by carefully considering the relation between the K and k i vectors, as discussed below and sketched in Fig The determinant |M | of the transformation matrix defines the ratio between the super-cell and primitive cell volumes in the direct (V and v) and reciprocal (W and w) spaces: |M | = V /v = w/W . The eigenstates of a single K point in the supercell Bz correspond to eigenstates of the primitive cell folded from different primitive cell points k i , with i running from 1 to |M | (see Fig. 2); these points are connected by linear combinations of supercell reciprocal lattice vectors {G} i : The folding problem can be equivalently expressed by considering that the eigenstates of one k i vector fold to one unique K point in the supercell first Brillouin zone, determined by one specific combination of supercell reciprocal lattice vectors {G} 0 : 7 see also the straight arrow in Fig. 2. The two equations above allow for an efficient mapping of supercell and primitive cell reciprocal spaces. We implemented the possibility to limit the calculation of Bloch character to specific (K, k i ) pairs fulfilling Eq. 4 for any supercell K vector defined in input, excluding all other pairs which would trivially result in P Km (k i ) = 0. This restriction reduces considerably the computational effort, as fewer Bloch characters need to be evaluated: the number of evaluated characters for states on any K is given by |M |.
Additionally, the calculation of the Bloch character can be further limited to selected k i vectors of interest. In fact, the user is typically interested in retrieving the eigenstates for selected k i vectors from the folded supercell states, rather than exploring all contributions to the supercell K vectors. Therefore, the k i vectors can be initialized by the user, then they get automatically translated in the supercell reciprocal space by the transformation where K B and k b i represent the vector coordinates expressed in the supercell and primitive cell reciprocal spaces, respectively. The calcu-lation of the Bloch character can be then executed as in the original implementation, but it is limited to (K, k i ) pairs satisfying Eq. 5: This approach drastically reduces the computational effort of the algorithm, since the P Km (k i ) character needs to be evaluated on only one single K for any given k i (see also discussion in the Benchmark Section).
The implementation of these automatized features simplifies the initialization of the unfolding calculation for the user. Moreover, the primitive cell lattice vectors can also be automatically determined from the supercell by the program, by simply inverting Eq. 3, if the transformation matrix M is specified in input: The extraction of the output data is quite straightforward as well. In order to further facilitate the analysis of unfolding calculations, we make available a post-processing package for band structure analysis (bands4vasp), 12 that can also be used for the construction of unfolded band structures, Fermi surfaces and spectral functions, as well as the automatic calculation of Fermi vectors (i.e., the k i vector of eigenstates at the Fermi level). We recall that, in the framework of unfolding calculations, the spectral function A is approximated as where δ(E m − E) are Dirac delta functions centered around E m energies.

Benchmark and Results
We tested our implementation of the unfolding algorithm embedded in VASP by considering the electronic properties of BaFe 2(1-x) Ru 2x As 2 (with x = 0 and 0.25 for the undoped and doped cases, respectively). This material is indeed a good testbed, since many reference data are available in literature. [19][20][21] We take this opportunity also to describe features included in the post-processing bands4vasp package, such as the visualization of band structures, Fermi surfaces and spectral function, and the calculation of Fermi vectors (Sec. 3.1). Finally, we show novel and more challenging applications of our machinery: by performing calculations on large cells modeling the TiO 2 (110) surface, the unfolding analysis reveals formation of flat bands originating from in-gap polaronic states that are perturbed by the interaction with CO adsorbates deposited on the material surface (Sec. 3.2). Moreover, we describe the band splitting occurring in the transition process from non-collinear paramagnetic to ferromagnetic ordering in EuCd 2 As 2 (Sec. 3.3).

Analysis of metallic states in
The metal-to-superconductor transition driven by Ru doping in the BaFe 2 As 2 pnictide has attracted wide interest in both theoretical and experimental communities, [22][23][24][25][26][27][28][29] and the unfolding technique has proven itself largely useful to support density-functional theory investigations: One of the most evident results observed from the effective band structures is the progressively closure of hole pockets upon doping, due to a coupling with structural distortions. 20, 21 We performed spin-unpolarized DFT calculations on BaFe 2 As 2 by maintaining a similar computational setup as in Ref. 19, but applying the updated version of the unfolding algorithm. We modeled the BaFe 2 As 2 structure adopting supercells of different size (A 2 ,A 8 ,A 16 ), constructed by applying the following transformation matrices (as described in Eq. 3) The M 2 matrix represents the transformation from primitive cell to conventional unit cell for compounds with the body-centered tetragonal I4/mmm space group such as BaFe 2 As 2 : 30 We remark that the M matrix elements are required to be integer in the unfolding formal-ism, but this requirement does not prevent us to model rotations or non trivial transformations. 9 Figure 3: Effective band structure of the supercell constructed by the M 2 transformation matrix (gray-to-green color gradient), compared to the band structure of primitive cell (blue).
We performed a preliminary test by considering the pure (undoped) BaFe 2 As 2 crystal and comparing the effective band structure obtained from the supercell with the bands calculated directly for the primitive cell. In fact, no difference should be observed between the two band structures, once the supercell states are unfolded into the reciprocal space of the primitive cell. Figure 3 shows the unfolded bands of the supercell A 2 constructed by the transformation matrix M 2 . As expected (due to |M 2 | = 2), we obtained for the supercell double as many of the bands obtained directly from the primitive cell calculation, shown in the figure for comparison. The unfolding algorithm is able to correctly identify the band folding, and to assign the Bloch character P Km accordingly: the gray color in the gradient palette identifies the folded bands with P Km = 0 that belong to different points in the primitive cell reciprocal space. Usually, the folded bands should not be considered in the analysis of the electronic properties of the material, and can be omitted from the graph (e.g., by setting the gradient palette with P Km = 0 to the same color as the background of the image). Figure 4 shows the effects of 25% of Ru doping on the electronic properties of BaFe 2(1-x) Ru 2x As 2 (x = 0.25). The band structure in Fig. 4a correctly reproduces the closure of the hole pockets around the Γ and Z high symmetry points. 20 The intersection of the bands with the Fermi level can be identified automatically, by using the bands4vasp tool (see inset in Fig. 4a): the eigenstates around Fermi are assigned to different bands by considering the Bloch character and the orbital symmetry, and the intersection of every band with Fermi is found by an interpolation. This feature is extremely useful for the analysis of the Fermi vectors, that get progressively shortened for the hole pockets around Γ upon Ru doping in this material. Moreover, by collecting all Fermi vectors in the Brillouin zone, it is possible to construct 3-dimensional Fermi surfaces (see right panel in Fig. 4a). The 2dimensional cut of the basal plane around Γ (see inset in Fig. 4a) highlights the importance of an accurate sampling of the reciprocal space, by comparing the Fermi surface obtained using two different approaches: The right side of the Fermi surface is constructed by using a radial distribution of the k points, leading to a better resolved description of the states around Γ, as compared to the resolution obtained by adopting a conventional rectangular grid (left side).
The bands4vasp tool can also efficiently extract the atomic orbital character of the eigenstates from unfolding calculations. Figures 4bd show the effective band structures with Bloch character represented by the size of the circles, and the color gradient representing the projection of the states on the d yz (Fig. 4b) and d x 2 −y 2 (Fig. 4c) orbitals (essentially due to Fe atoms), and the overall contribution from As atoms (Fig. 4d). The contribution of Fe states around Γ (d yz , and d xz not shown here), Z (d x 2 −y 2 ) and on the Bz border stems out clearly from both the band structure and the corresponding 3-dimensional Fermi surfaces. Similarly, the hybridization with As atoms leads to states at the Fermi level around Z and X.
By looking carefully at the d x 2 −y 2 orbital (in Fig. 4c), we note a band crossing the Fermi level in the Γ-Z direction. This feature is highlighted in Figure 5, that shows 2-dimensional Fermi surfaces for the basal plane and the parallel planes along Γ-Z, obtained for the spectral  Fig. 5a). Conversely, the spectral function for the inner state (with d x 2 −y 2 orbital symmetry) is absent, and becomes progressively better defined when moving towards Z (Fig. 5b,c). The spectral function allows us also to easily identify band degeneracy and crossing points between different bands, that are revealed by a more in-tense value: In Fig. 5b the crossing of the d xz and d yz states determines four points with high value for the spectral function around the center of the plane. This crossing corresponds to a progressive band switching between d xz and d yz states, clearly identifiable by looking at the evolution of the orbital symmetry of the internal and external rings moving from Γ to Z.
We conclude our benchmark by commenting on memory requirements of the unfolding algo- rithm. The automatic determination of (K, k i ) pairs (as in Eqs. 5 and 6) reduces the computational effort, by limiting the unfolding procedure only to points of interest. For small systems, as the A 2 supercell, we counted a memory gain of about 20% when using 300 k points, which increases to approximately 30% and 50% for the larger A 8 and A 16 cells, respectively. The lower computational requirements are very useful for performing unfolding calculations that model complex systems: some original application of the unfolding method are presented in the following sections.

In-gap polaronic states and adsorbates on TiO 2 (110) surface
The unfolding algorithm can be applied also for DFT calculations on systems with reduced dimensionality, such as the surface slab shown in Figure 6, modeling the pristine rutile TiO 2 (110) termination. Unit cells modeling surfaces of solids in VASP contain a vacuum region in order to interrupt the periodicity of the system along the surface normal. 31 Typically, several atomic layers are also included in the model in order to mimic the properties of the bulk below the surface (see Fig. 6a). Primitive cells and supercells share the same vector perpendicular to the surface as in Fig. 6a; in order to model sur-face reconstructions or defects, supercells are constructed by enlarging the lateral size (see Fig. 6b), leading to folding of electronic states in analogy with the bulk. Rutile TiO 2 (110) supercells can be used to model the formation of oxygen vacancies on the surface, that lead to stabilization of small electron polarons, i.e., electrons strongly localized on Ti ions and coupled with the phonon field. [32][33][34] Small polarons are typically associated to eigenstates appearing in the energy band gap of semiconductors; moreover, polaron localization can occur on different sites, with different formation energy (in rutile, subsurface Ti ions are preferred, over surface sites). 35 Polarons are known to drastically affect the electronic and chemical properties of the hosting material, with substantial impact on the applications: we focus here on the chemical activity of rutile surface, by considering the interplay between polarons and CO adsorbates, and the effects on the eigenstates. 36 Figures 6c-e collect the results obtained for large 6 × 2 supercells containing 363 atoms, including two CO molecules and two polarons (technical details of the calculation in Ref. 36). As shown in Figs. 6c,d, the CO can adsorb on Ti sites at different distance from the polarons. The corresponding effective band structure (unfolded on the surface primitive cell) shows the appearance of the strongly localized polaronic states, revealed by two flat in-gap bands (one per polaron) in the majority spin channel (Fig. 6e). By looking closer at these in-gap bands (inset in Fig. 6e), we note that the two polaronic states are not degenerate, due to the interaction with the CO molecules. The band appears more perturbed for the polaronic state closer to the CO molecule, as manifested by the increased band width. Perturbations of the polaronic in-gap states may originate also from the repulsive interaction of polarons at small distance, as described in Ref. 37 for polarons and bipolarons.

Non-collinear ferromagnetic fluctuations in EuCd 2 As 2
We describe here the application of the unfolding algorithm on systems with noncollinear magnetic ordering. 15 We consider the paramagnetic-to-ferromagnetic transition in EuCd 2 As 2 , an interesting semimetal showing emergence of Weyl fermions in the para-magnetic phase due to spin fluctuations of Eu magnetic moments. 38 Figure 7 compares the effective spectral functions calculated for EuCd 2 As 2 with different magnetic orderings (technical details of calculations are described in Ref. 38). The paramagnetic phase was modeled by a large supercell including 16 Eu atoms with magnetic moments fixed to random orientations, resulting in a vanishing total magnetization: the corresponding spectral function unfolded in the reciprocal space of the primitive cell is shown in Fig. 7a. The flat f bands of Eu atoms appearing around −1.5 eV show an evident incoherence, due to the random orientation of the magnetic moments. The three p bands of As atoms appear strongly spin-degenerate: the spectral function is very effective in capturing the band degeneracy, as degenerate bands result in higher values of the spectral character integrating the contribution from every state, as described in Eq. 7 (at variance with band structures, showing instead the Bloch character of every state The ferromagnetic phase shows interesting changes (Fig. 7c). First, we note that the f bands are more coherent, as expected, due to the ferromagnetic alignment of all Eu magnetic moments. Remarkably, the ferromagnetic order induces a splitting of the p states of As atoms: we observe indeed six bands, lifting the spin degeneracy of the three p bands in the paramagnetic phase (note also the lower spectral function value, as compared to the paramagnetic case).
Although the study of ferromagnetic systems could be done directly in the primitive cell, the supercell approach allowed us to study the paramagnetic-to-ferromagnetic transition by considering ferromagnetic domains embedded in a paramagnetic environment. In Fig. 7b, we show the spectral function of the system including a large ferromagnetic domain (consisting of 10 Eu atoms with aligned magnetic moments), and a smaller region (6 atoms) with Eu magnetic moments constrained to random directions. The splitting of p orbital persists in this transition state: This is an example of the effect of spin fluctuations on the param-agnetic phase of the compounds. In smaller ferromagnetic domains, the band splitting is gradually reduced, progressively converging towards the paramagnetic degeneracy (results obtained by modeling different size of the ferromagnetic domain are available in Ref. 38).

Conclusions
In summary, we report here our optimization of the unfolding scheme embedded in VASP, characterized by a simplified user interface, and reduced memory requirements, thanks to an efficient mapping between the reciprocal spaces of the supercell and primitive cells. The construction of effective band structures, spectral functions, Fermi surfaces and projections of electronic states on orbitals and ions, is further facilitated by the bands4vasp post-processing package.
The unfolding scheme is extremely useful in the interpretation of the results obtained by supercell approaches, and it facilitates the comparison with the experimental observations, especially in the field of spectroscopy. The application range is very broad: We take here the BaFe 2(1-x) Ru 2x As 2 superconductor as a benchmark, given the large amount of data available in literature. Moreover, we considered the adsorption of CO molecules on the rutile TiO 2 (110) surface, in order to show the suitability of the algorithm for very large supercells, such as those required in surface science calculations. In this case study, the effective band structures highlight the interactions of adsorbates with strongly localized polarons, revealed by perturbation of the flat polaronic bands. Finally, we performed non-collinear calculations for the EuCd 2 As 2 semimetal: The supercell approach allowed us to model the paramagnetic phase, by constraining magnetic moments along random directions, resulting in a vanishing total magnetization. The corresponding spectral function reveals the spin degeneracy of shallow states below Fermi, that is lifted by including spin-fluctuations via formation of ferromagnetic domains.
The unfolding algorithm proposed here rep-resent a useful computational tool for a wide range of physical and chemical phenomena requiring very large supercells, thanks to the improved interface, the reduced computational requirements and the powerful analysis tool.
graphene: Retaining an effective primitive cell band structure by band unfolding.