{'fraction_non_alphanumeric': 0.027602478821991017, 'fraction_numerical': 0.010375803058729888, 'mean_word_length': 4.680445664459874, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Unfolding method for the first-principles LCAO electronic structure calculations

Unfolding the band structure of a supercell to a normal cell enables us to investigate how symmetry breakers such as surfaces and impurities perturb the band structure of the normal cell. We generalize the unfolding method, originally developed based on Wannier functions, to the linear combination of atomic orbitals (LCAO) method, and present a general formula to calculate the unfolded spectral weight. The LCAO basis set is ideal for the unfolding method because of the invariance that basis functions allocated to each atomic species are invariant regardless of existence of surface and impurity. The unfolded spectral weight is well defined by the property of the LCAO basis functions. In exchange for the property, the non-orthogonality of the LCAO basis functions has to be taken into account. We show how the non-orthogonality can be properly incorporated in the general formula. As an illustration of the method, we calculate the dispersive quantized spectral weight of ZrB2 slab and show strong spectral broadening in the out-of-plane direction, demonstrating the usefulness of the unfolding method.


I. INTRODUCTION
The Kohn-Sham (KS) framework 1 within the density functional theory (DFT) allows us to investigate a wide variety of imperfect materials such as surface, impurities, and vacancies. 2,3 A widely used method to perform first-principles calculations of such systems is to introduce a supercell which makes studies of various forms of imperfections possible. 2,3 Not restricted by the periodic boundary condition, the Bloch theorem can also be applied to study a non-periodic system by introducing a large supercell which simulates a system where the translational symmetry is highly broken, e.g., the presence of a surface. 4 However, there are at least two drawbacks introduced by the large supercell in analyzing the electronic structure. First, the bands folded heavily in the small first Brillouin zone (BZ) corresponding to the large supercell makes it difficult to analyze how symmetry breakers such as surfaces and impurities perturb the band structure of the normal cell, where by the normal cell we mean a unit cell which is smaller than the supercell, less imperfect, and gives a reference of band structure. Second, it would be difficult to directly compare the heavily folded bands with experimental results. For example, the band structure calculated for the supercell cannot be directly compared with the spectra measured by the angle resolved photoemission spectroscopy (ARPES) without further considering the proper spectral weight, the imaginary part of the retarded one-particle Green function. 5,6 It would be desirable to develop a method which represents the band structure or spectral function of the supercell in terms of eigenstates of a chosen normal cell in order to relieve the two drawbacks in the supercell calculations. Unfolding methods have been proposed as an idea of realizing the change of representation of spec-tral function in terms of the eigenstates of the normal cell, and implemented in a wide variety of ways. For example, one can efficiently and rigorously calculate the band structure of the normal cell from a Γ-point calculation of the large supercell without imperfection by using maximally localized Wannier functions (WFs). 7 Another exact method for unfolding the band structure of a perfect supercell into a bulk dispersion relation, which can be compared to experiments, has been proposed in tightbinding calculations, and was further applied to an imperfect supercell via an averaged Hamiltonian. 8 Unfolding methods have also been introduced to the plane-wave basis sets. 9,10 Recently, a new unfolding approach by using energy-resolved symmetry-respecting WFs has been proposed by representing the spectral function from the supercell calculation in the basis of a conceptual normal cell instead of the representation by the eigenstates of the supercell. 11 The approach allows us to uniquely determine the unfolded spectral weight for a chosen normal cell via WFs and the geometrical structure of the supercell.
While the unfolding methods based on the WFs have been successful in providing detailed physical insights on various systems, 7,11-14 the WFs need to be constructed in the actual implementation. 15,16 For large-scale systems the construction of the WFs can be time-consuming, which may hamper the applicability of the unfolding methods to large-scale systems. In addition to this, one may encounter a difficulty that one to one correspondence between WFs defined in the supercell and normal cell could not be well defined owing to the gauge freedom 15 during the construction of the Wannier functions. Since the most essential step in the unfolding method based on the WFs is to determine one to one correspondence between WFs defined in the supercell and normal cell, arXiv:1211.4690v3 [cond-mat.mtrl-sci] 17 May 2013 the gauge freedom has to be utilized so that the symmetry and shape of WFs can be approximately identical between the two cells. 11 The second issue may also pose a difficulty in applying the unfolding method to systems with strong perturbation. Thus, it would be physically more preferable if the localized basis functions are identical for the same atomic species in all the normal cells arranged in the supercell by following the proposed unfolding method. 11 For the unfolding method, the linear combination of atomic orbitals (LCAO) method [17][18][19][20][21][22] can be regarded as an ideal framework in which the same basis functions are allocated for each atomic species regardless of existence of surface and impurity. It is apparent that one can easily establish one to one correspondence between AOs located in the supercell and normal cell without any ambiguity. In this paper we extend the unfolding method, 11 originally developed based on Wannier functions, to the LCAO method, and present a general formula to calculate the unfolded spectral weight for a chosen normal cell. In exchange for the invariance of the LCAO basis functions, the non-orthogonality between the LCAO basis functions has to be taken into account. We show how the non-orthogonality can be properly incorporated in the general formula. In addition to the ideal property of the LCAO method to the unfolding method, efficient computational methods have been developed by making use of the locality of the LCAO basis functions. Thus, it can be expected that the generalization to the LCAO method extends the applicability of the unfolding method to large-scale systems. [17][18][19][20][21][23][24][25] As an illustration of the method, we apply the method for recovering the bulk dispersion relation in the out-ofplane direction from a DFT electronic structure calculation for a ZrB 2 slab with the (0001) surface, where only one k-point sampling is needed along the direction. A graphene counterpart, silicene, has recently been epitaxially grown on the ZrB 2 (0001) surface. 26 Therefore, it would be important to investigate the surface and slab states of ZrB 2 in order to deeply understand the newfound silicene which could be constituent of future devices as well as extensively studied graphenes. [26][27][28][29] It is shown that the unfolding method reveals the dispersive quantized spectral weight of the ZrB 2 slab and strong spectral broadening in the out-of-plane direction, which is expected to be measured by experiments.
The paper is organized as follows: In Sec. II, the concept of the unfolding method is introduced. After introducing the concept, the unfolding formula is derived in Sec. III. In Sec. IV, we illustrate how to calculate the unfolded spectral weight along the out-of-plane direction in the bulk BZ from a slab calculation for ZrB 2 where a dispersive quantized spectral weight with strong broadening can be revealed. Finally, we summarize our study in Sec. V.

II. THE CONCEPT OF UNFOLDING BAND STRUCTURES
The concept of the unfolding method is introduced in this section. The physical meaning of the conceptual normal cell is also discussed here, which is the key element for the unfolding method. With the introduction of the conceptual normal cell, the method can uniquely determine the spectral weight in the BZ of the conceptual normal cell, which directly presents the strength of the translational symmetry breaking to the band structure of the normal cell. The concept is illustrated with a two-dimensional model system of which primitive cell contains one atom per unit cell as shown in Fig. 1 (a). It is assumed that the Fermi point is located at the Γ point and only a single eigenstate, whose spectral weight is exactly one, is involved at the Fermi point as shown in Fig. 1 (b), where we used the Fermi point instead of the Fermi surface due to the involvement of only the Γ point. To simplify the illustration, we focus only on the spectral weight on the Fermi point, but the idea can be applied for all the k-dependent eigenstates.
To study the same system, one can perform a perfect supercell calculation, for example, the 4×4 supercell as shown in Fig. 1 (c). As long as there is no additional symmetry breaking, the physical properties should be exactly the same as those of the primitive cell. Therefore, the primitive unit cell is referred to as the normal cell. However, it is confirmed that the Fermi points obtained by the supercell calculation as shown in Fig. 1 (d) are different from those of the normal cell. Because the first BZ of the supercell is smaller and is translationally symmetric upon arbitrary integer shifts of reciprocal lattice vectors (G), the repeated counterparts make the comparison difficult. In order to have a direct comparison between these two calculations, one needs to unfold the band structure from the first BZ of the supercell (the reduced-zone representation) to the BZ of the normal cell. In other words, we would like to recover a proper spectral weight from the periodic-zone representation to the extended-zone representation. Obviously, the unfolded band structure is invariant against any arbitrary choice of the supercell without perturbation . In this case, the unfolded band structure should exhibit the same Fermi point as shown in Fig. 1 (b).
In case that the system undergoes translational symmetry breaking via the imperfections as shown in Fig. 1 (e), it is apparent that first-principles calculations have to be performed for the supercell instead of a smaller unit cell. However, it would be expected that a similar spectral weight to Fig. 1 (b) as illustrated in Fig. 1 (f) is recovered by introducing the normal cell if the symmetry breaking is not strong just like the case of the 4×4 perfect supercell. It is also clear that the spectral weight of the supercell eigenstate (exactly one) may not be directly used to compare with ARPES. The measured weight cannot be suddenly modified by the introduction of a weak perturbation. That is, it is possible that the periodicity of the supercell BZ is not the periodicity of the observed spectral weight. Obviously, how to obtain the proper spectral weight via unfolding is an important issue. A solution to do that is simple at least conceptually. We only have to represent the spectral function obtained by the supercell calculation in the eigenstates of the normal cell. However, the definition of the normal cell may not be so obvious in general, since the translational symmetry is broken in the supercell. In the case as shown in Fig. 1 (e), a rather obvious choice is the unit cell plotted with the dashed lines as the same as the primitive unit cell as shown in Fig. 1 (a). However, that is not the only choice. In fact, the normal cell in the unfolding method is regarded as a conceptual unit cell which defines one to one correspondence between localized functions such as WFs and AOs in the supercell and normal cell. If a proper unit cell corresponding to the periodicity of the observed spectral weight is chosen as the normal cell, the resultant spectral weight in the corresponding BZ can be compared to ARPES after the polarization dependent dipole matrix element is taken into account within the sudden approximation. 11 If a reference unit cell is chosen as the normal cell, the resultant spectral weight will provide information of the translational symmetry breaking to the reference system. Therefore, it would be anticipated that the unfolding method can be utilized for a wide variety of studies.
In the current approach, once a unit cell is chosen as the normal cell, we only have to consider one to one correspondence between AOs in the supercell and normal cell in representing the spectral function of the supercell in the eigenstates of the normal cell. Such an assignment is straightforwardly performed in the LCAO method since the same LCAO basis functions are allocated to each atomic species. Then, the information of the translational symmetry breaking in the supercell is stored in two aspects. One is that the symmetry breaking is recorded in the LCAO coefficients of the supercell eigenstates. The other is that the symmetry breaking is built in the position of the basis functions in the supercell. The latter is taken into account by overlap integrals between the LCAO basis functions. Since the eigenstates of the supercell and the overlap integrals are obtained by a conventional first-principles calculation, we discuss in Sec. III how the proper spectral weight can be calculated from the two quantities in detail.

III. THE UNFOLDING FORMULA
In the section, we present a general formula to calculate the unfolded spectral weight. As discussed in the previous section, the proper spectral weight defined in the BZ of a chosen conceptual normal cell is evaluated by using the eigenstates of the supercell and the overlap matrix elements of the LCAO basis functions. The strength of each band's coupling to the symmetry breaker (e.g., impurities, vacancies, dopants, and lattice distortions) can be observed via the unfolded spectral weight.
Let us start to introduce the spectral functionÂ(ω) defined as the imaginary part of one-particle Kohn-Sham Green's function: where 0 + is a positive infinitesimal, and z an energy in the complex plane (z = ω + iη). |KJ denotes the Bloch state with the crystal momentum K and the band index J obtained by a DFT calculation for the supercell. Though the non-spin polarized case is considered in the paper for simplicity of notation, the generalization is straightforward. In the LCAO method, atomic basis functions {|RN } are placed in every unit cell specified with a translational lattice vector R, where N is a symbolic orbital index consisting of atomic position relative to R, a multiplicity index for radial functions, an angular momentum quantum number, and a magnetic quantum number. With the idea of the LCAO method, |KJ is expanded in a form of linear combination of atomic basis functions as with the definition: where L is the number of the unit cells introduced in the Born-von Karman boundary condition, and C KJ N the LCAO coefficient. Although we use the bracket notation to denote the Bloch state |KJ , |KN defined by Eq. (4), and the atomic basis function |RN , the distinction among them is made by alphabet, i.e., the first alphabet is K or R, and the second is J or N (or M ). Then, one of them is uniquely distinguished by the combination of two alphabets. The LCAO coefficients and the eigenvalue ε KJ of |KJ are calculated in the KS framework by solving a generalized eigenvalue problem: with the Kohn-Sham Hamiltonian H given by and the overlap matrix S(K) given by where the overlap matrix element S 0M,RN ≡ 0M |RN reflects the strength of the non-orthogonality between the LCAO basis functions. Since the Bloch state |KJ and the corresponding eigenvalue ε KJ are obtained by the conventional supercell calculation, the spectral function A(ω) is well defined in terms of the representation of the supercell.
We now consider to represent the spectral function in the eigenstates of a normal cell, which is what we mean by unfolding. As discussed in the Sec. II, the normal cell can be chosen as either a unit cell corresponding to the periodicity of spectral weight observed in experiments or a reference unit cell depending on the purpose under consideration. The choice of the normal cell is equivalent to the introduction of a conceptual system of which periodicity is the same as that of the normal cell. If no symmetry breaker is introduced in the supercell, there is no ambiguity for the introduction of the conceptual system. Thus, the unfolding is performed in a precise mathematical sense. For general cases with symmetry breakers, however, such a periodicity of the normal cell is not hold anymore. Nevertheless, we assume that such a conceptual system can be defined and the corresponding Bloch state |kj can be given by where |kn is the counterpart to Eq. (4). In the following discussion, uppercase and lowercase letters will be used for indices associated with the supercell and the normal cell, respectively. The usage of alphabets for the bracket notation in the normal cell follows that of the supercell. The substance of the conceptual system that we introduced here might be regarded as an average of all the normal cells in the supercell. Once the assumption is accepted, the spectral function can be expressed in the representation of the Bloch state |kj of the conceptual system as shown below. By noting that the closure relation in the non-orthogonal LCAO basis functions is given by with the corresponding dual function |kn defined by the spectral function A is expressed as where A kj,kj ≡ kj|Â|kj and Ak n,kn ≡ k n|Â|kn . It is worth mentioning that similar closure relations to Eq. (9) are hold for |KN , |RN , |kn , and |rn . By inserting a closure relation KJ |KJ KJ| in two places of the last formula of Eq. (11), and noting thatÂ is defined by {|KJ }, we obtain the following expression for each k: where the summation over j and J for the left-and righthand sides were dropped, respectively. The omission of the summations is possible by redefining j so that the state specified with j can correspond to the unfolded counterpart to the states specified with J. The overlap integrals km|KJ and KJ|kn appearing in Eq. (12) bridge two representations of the supercell and the normal cell. The integral km|KJ can be written by using Eq. (3) as We further rewrite km|KM in Eq. (13) by making use of two closure relations r |r m rm| and R |RM R M | as follow: where we used the dual orthonormality relation between the original and its dual functions, and l is the number of the unit cells introduced in the Born-von Karman boundary condition for the normal cell. As well, it is easy to see that the overlap integral KJ|kn is written by After we assumed the conceptual system and the Bloch states by Eq. (8), we have never introduced any approximation. All the derivations up to here are rigorous in a precise mathematical sense under the assumption of existence of the conceptual system. However, here we introduce two approximations in evaluating the overlap integrals rm|RM and R N |r n in Eqs. (14) and (16), respectively. The first integral rm|RM is evaluated by assuming that the position of |RM in real space is the same as that of the counterpart defined in the conceptual system, while the second integral R N |r n is evaluated by assuming that the position of |r n in real space is the same as that of the counterpart defined in the supercell system under investigation. The introduction of the approximations may be justified by noting that the true value of each integral is approximately given by the mean of the two approximated values if the two approximations are independently applied to the single integral. In fact, the conceptual system can be introduced by the approximations without addressing a specific system for the conceptual system as shown later on. In order to introduce the approximations, we need to establish one to one correspondence between AOs in the supercell and normal cell. In our case, the establishment of the one to one correspondence is easily realized due to the invariance of the LCAO basis functions. With the introduction of the two approximations, the second integral R N |r n is nothing but the overlap integral used in the supercell calculation. Thus, we only have to focus on the first integral rm|RM , where |RM needs to be relabeled in the representation of the normal cell. By relabeling R → R + r 0 (M ) and M → m (M ), we replace as Although the quantity is nothing but the overlap matrix between AOs in the representation of the normal cell, we further rewrite the integral to simplify the last formula of the spectral function A that we have been pursuing. By inserting the Fourier representation of |rm : and that of |R+r 0 (M ), m (M ) into Eq. (17), and noting that we obtain We are now ready to evaluate A kj,kj given by Eq. (12). By using Eqs. (13)- (16), and (20), the spectral function A kj,kj can be rewritten as It should be noted in the derivation that S −1 nm (k) in Eq. (12) and S mm (M ) (k) in Eq. (20) are replaced by a relation m S −1 nm (k)S mm (M ) (k) = δ nm (M ) . The replacement is a crucial step in our derivation because Eq. (21) does not require any information about the position of the LCAO basis functions defined in the conceptual system, which is the reason why we call the conceptual system. The formula of Eq. (21) still has room to be simplified. By noting that we finally obtain with (25) where A KJ,KJ (ω) is just a delta function at the eigenvalue, δ(ω − KJ ). In Eq. (25) the overlap integral S 0N,rm(M ) is evaluated by assuming the positions of LCAO basis functions used in the actual supercell calculation, which is one of our approximations that we made before. With the two approximations, one can calculate the spectral weight in terms of the normal cell representation without relying on the substance of the conceptual system. Physically, the use of the conceptual system is just to have a guidance in obtaining the proper spectral weight in the desired representation. As can be confirmed, the sum of weight of the unfolded counterparts for each supercell eigenstate is exactly one. It is also noted from Eq. (22) that the spectral function at k is contributed from only the folded counterpart A (k−G)J,(k−G)J . The formula of Eq. (24) allows us to separately calculate the unfolded spectral weight of each supercell eigenstate within a given energy window by using the eigenvalue, the LCAO coefficients, and the overlap matrix elements obtained from the supercell calculation. It becomes clear that the phase factor e ik(r−r (M )) plays an important role to determine a spectral weight of the extended-zone representation in unfolding the band structure of the supercell, and that the effect of translational symmetry breaking has been recorded in not only the LCAO coefficients, but also the overlap matrix elements in the supercell calculation. The formula is general for the localized orbitals and can be applied for both orthogonal and non-orthogonal basis functions. If Eq. (24) is applied for a supercell without any symmetry breakers, the band dispersion of the chosen normal cell is recovered in a precise mathematical sense.
In addition, it is very interesting to see how each basis function contributes to the unfolded spectral weight. The decomposition to each contribution may be possible by defining W k KJM as follows: Then, the spectral function can be written as By decomposing the spectral weight with Eq. (26) into the LCAO basis functions, one can analyze which localized basis functions compose each band dispersion, and which bands correspond to spatially localized states such as surface states of a slab and vacancy states in a bulk. Finally, we summarize the practical procedure to unfold the band structure obtained by the supercell calculation. First, one needs to introduce a supercell and a normal cell, and define a rule of the relabeling for r (M ) and m(M ) in Eq. (25). The second step is to perform the band structure calculation of the supercell, where the LCAO coefficients and the overlap matrix elements are obtained. The third step is to calculate the spectral weight by using Eq. (25) or (26) with the LCAO coefficients and the overlap matrix elements calculated at the second step.

IV. ELECTRONIC STATE OF ZIRCONIUM DIBORIDE SLAB
A common way to identify the surface state is to compare the band structure of the slab model with that of the bulk and/or to highlight the atomic contribution of the terminated layer. The idea can hardly work for a large unit cell of studied surfaces, where the heavily folded bands are involved. The large supercell is especially needed due to a spontaneous reconstruction 30 or the interaction with an added layer. 26,29,31 The unfolding method can provide a way to the direct comparison of the band structures in the same BZ, which enables the visualization of the symmetry breaking and the surface states. In the section, we demonstrate how our unfolding method can be applied to study the electronic states of a ZrB 2 slab, consisting of 17 Zr and 16 B layers, with the Zr-terminated (0001) surface, where the unit cell of bulk ZrB 2 is used as the normal unit cell. Although the electronic states of Zr-terminated ZrB 2 (0001) surface have FIG. 2: (a) The hexagonal unit cell of bulk ZrB2 is presented with one Zr atom at (0,0,0) and two B atoms at (1/3,2/3,1/2) and (2/3,1/3,1/2). (b) The first BZ of the bulk ZrB2. (c) The electronic band structure of bulk ZrB2 along the path indicated in (b). (d) The calculated 33-layer Zr-terminated slab is shown with the (1×1) unit cell of ZrB2(0001) surface. The dashed rectangle gives an idea how thick the vacuum is. For the unfolding method, we choose the unit cell of bulk ZrB2 as the normal cell. Since only the difference between lattice vectors are needed in the unfolding formula, we show one example for relabeling supercell atoms by the normal cell lattice vectors and atoms. The basis functions will follow the assignment of the corresponding atoms. been well studied, 32,33 we are motivated by a recent fabrication of silicene on the ZrB 2 (0001) surface to reinvestigate the surface states in detail. 26 The DFT calculations were performed by the OpenMX code which is based on norm-conserving pseudopotentials generated with multiple reference energies 34 and linear combination of optimized pseudoatomic basis functions. 22,35 The cutoff radius of 7.0 Bohr was used for all the basis functions. For each Zr atom, three, three, and two optimized radial functions were allocated for the s-, p-, and d-orbitals, respectively, as denoted by Zr-s3p3d2, while B-s4p2d1 were allocated for each B atom. The regular mesh of 270 Ry in real space was used for the numerical integrations and solution of the Poisson equation. 19 The exchange-correlation energy functional was treated with a generalized gradient approximation (GGA) proposed by Perdew, Burke and, Ernzerhof. 36 Geometrical structures were optimized until the maximum force on atom becomes less than 3 × 10 −4 Hartree/Bohr. In order to validate the accuracy of the DFT methods for the system, we optimized the lattice constants of the bulk ZrB 2 , where the primitive cell of the bulk ZrB 2 is hexagonal with one Zr atom at (0,0,0) and two B atoms at (1/3,2/3,1/2) and (2/3,1/3,1/2). The optimized lattice constants of bulk ZrB 2 were found to be a = 3.174Å and c = 3.550Å with a k-point sampling of 8×8×5 mesh, which are in good agreement with experimental data (a = 3.170Å and c = 3.533Å) taken from Ref. 37. The geometrical structure and the band structure of the bulk ZrB 2 are shown in Fig. 2.
As a next step, we performed geometry optimization of the ZrB 2 slab with a k-point sampling of 8×8×1 mesh, where the optimized lattice constants obtained by the bulk calculation were used for the a-and b-axes, and the vacuum of 14.292Å were inserted to avoid the interaction between the slabs. Since the unfolding formula can be applied to the calculation with or without adding basis functions in the vacant region, no additional basis functions are added into the vacuum in this study and one can think the corresponding LCAO coefficients and the overlap matrix elements are zero. The optimized structure exhibits that there is no surface reconstruction for the ZrB 2 (0001) surface that has been checked with the (2×2) unit cell of ZrB 2 (0001) surface, and the outermost Zr-Zr distance is 3.492Å which is slightly shorter than the bulk value (3.550Å). Up to this point, the overlap matrix elements between LCAO basis functions had already been obtained. In order to study the electronic band structure in the bulk representation, we need to employ the unfolding method by performing the band structure calculation along the path shown in Fig. 2 (b) where the LCAO coefficients are obtained. An additional effort in evaluating Eq. (25) or (26) is to assign the supercell basis functions to the normal cell units. One example is provided for the relabelings as shown in Fig. 2 (d). Note that only the difference between lattice vectors are needed in the implementation.
Before showing the unfolded band structure, we re-mark several differences of our unfolding method in understanding the surface state from the conventional study which shows only the atomic weight of the folded band structure. First, the unfolded spectral function itself is physically meaningful and the developed broadening is related to the quasiparticle lifetime of an electron in the unfolded |kj state. 5,6 Our method not only shows the total weight but also allows a detailed analysis of the contribution of each atomic basis function to the unfolded spectral function. Since a thicker vacuum would break the translational symmetry of the bulk more, the unfolded weight reflects the relative lifetime for each reference |kj electron that can survive in the slab system. Secondly, even without the surface reconstruction, the folded bands due to the slab thickness still exist, which could result in many shadow bands and/or prominent splittings in the original bulk bands. The unfolding method recovers the proper spectral weight and make the visualization of the effect of symmetry breaking straightforwardly. One more interesting advantage is the capability of recovering the dispersion along the out-of-plane direction which is almost impossible for the conventional electronic band structure calculation to reveal such a dispersion. Recall that only one k-point sampling is needed along that direction in giving rise to horizontal bands. We now present the unfolded band structure in Fig. 3. The result clearly shows that the overall bulk band structure is well maintained, accompanied with some shadow bands and several additional interesting features. Without calculating the proper spectral weight, it is difficult to identify the shadow band since each supercell eigenstate presents an equal amount of spectral weight. One interesting feature is the newly revealed surface bands in comparison to the bulk bands [c.f. Fig. 2 (c)], for example, the bands labeled by S1 and S2 which consist of large contribution of the outermost Zr d orbital. The electrons of the S1 band show diminishing spectral weight while the path approaches to the Γ point, indicating negligible quasiparticle lifetimes and their short mean free paths. Note that the outermost atomic contribution would show a diminishing behavior while approaching to the Γ point for the S1 band without unfolding. 32 However, our unfolding result demonstrates that the negligible spectral weight is expected even with the consideration of the contribution of all the layers.
For the dispersion in the out-of-plane direction, horizontal bands stemmed from the Γ-and M -point calculations are found, and can be identified as the slab states which are the quantized states due to electron confinement caused by the presence of the surface. 38,39 Our unfolding result shows that the discrete energy levels reveal a dispersive behavior as found in Fig. 3. The feature of dispersive quantized spectral weight can be understood in such a way that the thickness of the slab must construct the bulk property to some degrees. As a result, both the bulk dispersion and the quantized levels are observed simultaneously. One can also see strong spectral broadenings in the out-of-plane direction from the The outermost Zr d contribution is colored red and the sum of all Zr and B orbital contributions is colored black. The radii of the red circles are multiplied by ten for a better visualization. The right panel shows the sum of all the unfolded total weight by a contour plot. The bands labeled by S1 and S2 can be easily identified as surface bands. The originally horizontal bands along the out-of-plane directions show a dispersive behavior via unfolding and are found to follow the bulk dispersion, for example, the band labeled by S3.
center of the reference bulk spectral weight [c.f. Fig. 2  (c)]. This shorter quasiparticle lifetime at large deviation from the bulk band is expected, which is missing in the conventional supercell band structure calculation. Of course, the lifetime is expected to be much longer at the bulk correspondence for a thicker slab, where the slab eigenstate is closer to the bulk eigenstate. All the interesting findings shown in Fig. 3 are easily obtained by the unfolding method. Note that the quantum confinement effect has already been observed in multilayer graphene by ARPES. 27

V. SUMMARY
By making use of a conceptual normal cell, we have derived a general formula in unfolding first-principles band structures, which can be applied to the LCAO basis functions being an ideal basis set for the unfolding concept. The translational symmetry breaking is recorded in the LCAO coefficients of the supercell eigenstates and built in the real-space locations of the basis functions. Once a conceptual normal cell is defined, the unfolded spectral weight is properly delivered via the phases of differences between normal cell lattice vectors, the LCAO coefficients, and the overlap integrals between basis functions. The key effort for the unfolding procedure is just to define the conceptual normal cell in real space. We have applied the method for studying the electronic states of ZrB 2 (0001) slab along both the in-plane and the out-ofplane directions. A dispersive quantized spectral weight of the slab states with strong spectral broadenings in the out-of-plane direction can be easily revealed via unfolding. This interesting behavior should be found in all slabs in general, especially the thinner ones with layered structures. Since the LCAO basis functions are efficient in computation and the obtained supercell LCAO coefficients can be directly used in unfolding the band structure, various applications including large-scale systems can be expected by using our unfolding method for studying the translational symmetry breaking.