Unfolding method for first-principles LCAO electronic structure calculations

The Kohn–Sham (KS) framework [1] within the density functional theory (DFT) allows us to investigate a wide variety of imperfections in materials such as surfaces, 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., by the presence of a surface [4]. However, there are at least two drawbacks introduced by the large supercell in analysing the electronic structure. First, as the bands are folded heavily in the small first Brillouin zone (BZ) corresponding to the large supercell it is difficult to analyse 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 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 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 alleviate the two drawbacks in supercell calculations. Unfolding methods have been proposed as an idea for realizing a change of representation of the spectral 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 tight-binding calculations, and has been further applied to an imperfect supercell via an averaged Hamiltonian [8, 9]. Unfolding methods have also been introduced to plane-wave basis sets [10, 11]. 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 [12]. This 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 insight into various systems [7, 12–15], the WFs need to be constructed in the actual implementation [16, 17]. 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 [16] 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, the gauge freedom has to be utilized so that the symmetry and shape of the WFs can be approximately identical between the two cells [12]. 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 were identical for the same atomic species in all the normal cells arranged in the supercell by following the proposed unfolding method [12]. For the unfolding method, the linear combination of atomic orbitals (LCAO) method [18–23] can be regarded as an ideal framework in which the same basis functions are allocated for each atomic species regardless of the 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 [12], 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 for the unfolding method, efficient computational methods have been developed by making use of the locality of the LCAO basis functions [18–26]. Thus, it can be expected that the generalization to the LCAO method will extend the applicability of the unfolding method to large-scale systems.

As an illustration of the method, we apply the method for recovering the bulk dispersion relation in the out-of-plane direction from a DFT electronic structure calculation for a ZrB₂ slab with a (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₂(0001) surface [27]. Therefore, it would be important to investigate the surface and slab states of ZrB₂ in order to deeply understand the newfound silicene which could be a constituent of future devices as well as extensively studied graphenes [27–30]. It is shown that the unfolding method reveals the dispersive quantized spectral weight of the ZrB₂ 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 section 2, the concept of the unfolding method is introduced. After introducing the concept, the unfolding formula is derived in section 3. In section 4, 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₂, where a dispersive quantized spectral weight with strong broadening can be revealed. Finally, we summarize our study in section 5.

The concept of the unfolding method is introduced in this section. The physical meaning of the conceptual normal cell is also discussed here, as this is the key element of 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 of the band structure of the normal cell. The concept is illustrated with a two-dimensional model system in which the primitive cell contains one atom per unit cell, as shown in figure 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 figure 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.

Zoom In Zoom Out Reset image size

To study the same system, one can perform a perfect supercell calculation, for example, the 4 × 4 supercell as shown in figure 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 figure 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 figure 1(b).

In the case where the system undergoes translational symmetry breaking via imperfections, as shown in figure 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 figure 1(b), as illustrated in figure 1(f), would be recovered by introducing the normal cell if the symmetry breaking were 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, the means of obtaining the proper spectral weight via unfolding is an important issue. A solution to do this 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 shown in figure 1(e), a rather obvious choice is the unit cell plotted with the dashed lines, as it is the same as the primitive unit cell shown in figure 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 [12]. If a reference unit cell is chosen as the normal cell, the resultant spectral weight will provide information on the translational symmetry breaking in 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 on the translational symmetry breaking in the supercell is stored in two ways. 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 into 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 section 3 how the proper spectral weight can be calculated from the two quantities in detail.

3.1. Overview

3.2. LCAO method

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 the 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, the Bloch state |KJ〉 corresponding to the crystal momentum K and the band index J is expanded in the form of a linear combination of atomic basis functions as

$$\begin{eqnarray} &&\vert K J\rangle =\mathop{\sum }\limits _{N}C_{N}^{K J}\vert K N\rangle \end{eqnarray} \tag{ 1 }$$

with the definition

$$\begin{eqnarray} &&\vert K N\rangle =\frac{1}{\sqrt{L}}\mathop{\sum }\limits _{R}{\mathrm{e}}^{\mathrm{i}K R}\vert R N\rangle , \end{eqnarray} \tag{ 2 }$$

where L is the number of unit cells introduced into the Born–von Karman boundary condition, and ${C}_{N}^{K J}$ is the LCAO coefficient. Although we use the bracket notation to denote the Bloch state |KJ〉,|KN〉 defined by equation (2), and the atomic basis function |RN〉, the distinction among them is made by letter, i.e., the first letter is K or R, and the second is J or N (or M). Then, each of them is uniquely distinguished by the combination of two letters. The LCAO coefficients and the eigenvalue ε_KJ of |KJ〉 are calculated in the KS framework by solving a generalized eigenvalue problem

$$\begin{eqnarray} &&{H}^{K}{C}^{K}={\epsilon }^{K}S(K){C}^{K} \end{eqnarray} \tag{ 3 }$$

with the Kohn–Sham Hamiltonian H given by

$$\begin{eqnarray} &&{H}_{M N}^{K}=\mathop{\sum }\limits _{R}{\mathrm{e}}^{\mathrm{i}K R}\langle 0 M\vert \hat {H}\vert R N\rangle \end{eqnarray} \tag{ 4 }$$

and the overlap matrix S(K) given by

$$\begin{eqnarray} &&{S}_{M N}(K)=\mathop{\sum }\limits _{R}{\mathrm{e}}^{\mathrm{i}K R}{S}_{0 M,R N}, \end{eqnarray} \tag{ 5 }$$

where the overlap matrix element S_0M,RN ≡ 〈0M|RN〉 reflects the strength of the non-orthogonality between the LCAO basis functions. Note that the closure relation in the non-orthogonal LCAO basis functions is given by

$$\begin{eqnarray} &&\mathop{\sum }\limits _{K N}\vert \tilde {K N}\rangle \langle K N\vert =\mathop{\sum }\limits _{K N}\vert K N\rangle \langle \tilde {K N}\vert =I \end{eqnarray} \tag{ 6 }$$

with the corresponding dual function $\vert \tilde {K N}\rangle$ defined by

$$\begin{eqnarray} &&\vert \tilde {K N}\rangle =\mathop{\sum }\limits _{M}\vert K M\rangle S_{M N}^{-1}(K). \end{eqnarray} \tag{ 7 }$$

3.3. The change of representation of the spectral function

Let us start by introducing the spectral function $\hat {A}(\omega )$ defined as the imaginary part of the one-particle Kohn–Sham Green function

$$\begin{eqnarray} &&\hat {A}(\omega )=-\frac{1}{\pi }\mathrm{Im}\hat {G}(\omega +\mathrm{i}{0}^{+}) \end{eqnarray} \tag{ 8 }$$

with

$$\begin{eqnarray} &&\hat {G}(z)=\mathop{\sum }\limits _{K J}\frac{\vert K J\rangle \langle K J\vert }{z-{\varepsilon }_{K J}}, \end{eqnarray} \tag{ 9 }$$

where 0⁺ is a positive infinitesimal, and z an energy in the complex plane (z = ω + iη). Since the Bloch state |KJ〉 and the corresponding eigenvalue ε_KJ are obtained by the conventional supercell calculation within DFT, the spectral function $\hat {A}(\omega )$ is well defined in terms of the representation of the supercell. Although the non-spin polarized case is considered in this paper for simplicity of notation, the generalization is straightforward.

We now consider representing the spectral function in the eigenstates of a normal cell, which is what we mean by unfolding. As discussed in section 2, the normal cell can be chosen as either a unit cell corresponding to the periodicity of the 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 for which the periodicity is the same as that of the normal cell. If no symmetry breaker is introduced into the supercell, there is no ambiguity in 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 does not hold anymore. Nevertheless, we assume that such a conceptual system can be defined and the corresponding Bloch state |kj〉 can be given by

$$\begin{eqnarray} &&\vert k j\rangle =\mathop{\sum }\limits _{n}C_{n}^{k j}\vert k n\rangle , \end{eqnarray} \tag{ 10 }$$

where |kn〉 is the counterpart to equation (2). 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 letters 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 using the closure relation in the non-orthogonal LCAO basis functions as given by equation (6) for the normal cell, the spectral function A is expressed as

$$\begin{eqnarray} &&A=\mathop{\sum }\limits _{k j}{A}_{k j,k j}=\mathop{\sum }\limits _{k n}{A}_{\tilde {k n},k n}, \end{eqnarray} \tag{ 11 }$$

where ${A}_{k j,k j}\equiv \langle k j\vert \hat {A}\vert k j\rangle$ and ${A}_{\tilde {k n},k n}\equiv \langle \tilde {k n}\vert \hat {A}\vert k n\rangle$. In the concept of unfolding, |kj〉 simply represents one of the |KJ〉 counterparts in the first BZ of the conceptual unit cell. By inserting a closure relation ∑_KJ|KJ〉〈KJ| in two places in the last formula of equation (11), and noting that $\hat {A}$ is defined by {|KJ〉}, we obtain the following expression for each k:

$$\begin{eqnarray} &&{A}_{k j,k j}=\mathop{\sum }\limits _{m n K}S_{n m}^{-1}(k)\langle k m\vert K J\rangle {A}_{K J,K J}\langle K J\vert k n\rangle , \end{eqnarray} \tag{ 12 }$$

where the summations over j and J for the left- and right-hand sides have been dropped, respectively. Omission of the summations is mathematically possible if A_kj,kj does not rely on the substance of the conceptual system, which will be demonstrated in section 3.4. Then, one can redefine j so that the state specified with j can correspond to the unfolded counterpart to the state specified with J.

3.4. Simplification

In section 3.3, the spectral function A involves the information on the conceptual system, for example, |km〉,|kn〉, and the overlap matrix in the inverse form. Here, we will introduce two approximations to simplify the unfolding formula. The overlap integrals 〈km|KJ〉 and 〈KJ|kn〉 appearing in equation (12) bridge two representations of the supercell and the normal cell. The integral 〈km|KJ〉 can be written by using equation (1) as

$$\begin{eqnarray} &&\langle k m\vert K J\rangle =\mathop{\sum }\limits _{M}C_{M}^{K J}\langle k m\vert K M\rangle . \end{eqnarray} \tag{ 13 }$$

We further rewrite 〈km|KM〉 in equation (13) by making use of equation (2) for both 〈km| and |KM〉 as follows:

$$\begin{eqnarray} &&\langle k m\vert K M\rangle =\mathop{\sum }\limits _{r R}\frac{{\mathrm{e}}^{-\mathrm{i}k r}}{\sqrt{l}}\langle r m\vert R M\rangle \frac{{\mathrm{e}}^{\mathrm{i}K R}}{\sqrt{L}}, \end{eqnarray} \tag{ 14 }$$

where l is the number of unit cells introduced in the Born–von Karman boundary condition for the normal cell. Moreover, it is easy to see that the overlap integral 〈KJ|kn〉 is written by

$$\begin{eqnarray} &&\langle K J\vert k n\rangle =\mathop{\sum }\limits _{N}{C_{N}^{K J}}^{\ast }\langle K N\vert k n\rangle \end{eqnarray} \tag{ 15 }$$

with

$$\begin{eqnarray} &&\langle K N\vert k n\rangle =\mathop{\sum }\limits _{{r}^{\prime}{R}^{\prime}}\frac{{\mathrm{e}}^{-\mathrm{i}K{R}^{\prime}}}{\sqrt{L}}\langle {R}^{\prime}N\vert {r}^{\prime}n\rangle \frac{{\mathrm{e}}^{\mathrm{i}k{r}^{\prime}}}{\sqrt{l}}. \end{eqnarray} \tag{ 16 }$$

After assuming the conceptual system and the Bloch states by equation (10), we have not 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 equations (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 its 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 its 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 a 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 the 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 relabelled in the representation of the normal cell. By relabelling R → R + r₀(M) and M → m'(M), we replace as

$$\begin{eqnarray} &&\langle r m\vert R M\rangle \rightarrow \langle r m\vert R+{r}_{0}(M),{m}^{\prime}(M)\rangle . \end{eqnarray} \tag{ 17 }$$

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〉,

$$\begin{eqnarray} &&\vert r m\rangle =\frac{1}{\sqrt{l}}\mathop{\sum }\limits _{k}{\mathrm{e}}^{-\mathrm{i}k r}\vert k m\rangle , \end{eqnarray} \tag{ 18 }$$

and that of |R + r₀(M),m'(M)〉 into equation (17), and noting that

$$\begin{eqnarray} \langle k m\vert {k}^{\prime}n\rangle &=&\frac{1}{l}\mathop{\sum }\limits _{r}{\mathrm{e}}^{\mathrm{i}({k}^{\prime}-k)r}\mathop{\sum }\limits _{{r}^{\prime}}\langle r m\vert {r}^{\prime}n\rangle {\mathrm{e}}^{\mathrm{i}{k}^{\prime}({r}^{\prime}-r)}\nonumber\\ &=&{\delta }_{k{k}^{\prime}}{S}_{m n}(k), \end{eqnarray} \tag{ 19 }$$

we obtain

$$\begin{eqnarray} &&\langle r m\vert R M\rangle =\mathop{\sum }\limits _{k}\frac{{\mathrm{e}}^{\mathrm{i}k(r-R-{r}_{0}(M))}}{l}{S}_{m{m}^{\prime}(M)}(k). \end{eqnarray} \tag{ 20 }$$

We are now ready to evaluate A_kj,kj given by equation (12). By using equations (13)–(16) and (20), the spectral function A_kj,kj can be rewritten as

$$\begin{eqnarray} {A}_{k j,k j}&=&\mathop{\sum }\limits _{K M N{r}^{\prime}R{R}^{\prime}}\frac{{\mathrm{e}}^{-\mathrm{i}k(R+{r}_{0}(M)-{r}^{\prime})}}{l L}{\mathrm{e}}^{\mathrm{i}K(R-{R}^{\prime})}\nonumber\\ &&\times ~C_{M}^{K J}{{C}_{N}^{K J}}^{\ast }\langle {R}^{\prime}N\vert {r}^{\prime}{m}^{\prime}(M)\rangle {A}_{K J,K J}. \end{eqnarray} \tag{ 21 }$$

It should be noted in the derivation that ${S}_{n m}^{-1}(k)$ in equation (12) and S_mm'(M)(k) in equation (20) are replaced by a relation ${\sum }_{m}{S}_{n m}^{-1}(k){S}_{m{m}^{\prime}(M)}(k)={\delta }_{n{m}^{\prime}(M)}$. The replacement is a crucial step in our derivation because equation (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 it the conceptual system. The formula of equation (21) still has room to be simplified. By noting that

$$\begin{eqnarray} &&\frac{1}{L}\mathop{\sum }\limits _{R}{\mathrm{e}}^{\mathrm{i}(K-k)R}=\mathop{\sum }\limits _{G}{\delta }_{k-G,K}, \end{eqnarray} \tag{ 22 }$$

and

$$\begin{eqnarray} &&\mathop{\sum }\limits _{{r}^{\prime}{R}^{\prime}}{\mathrm{e}}^{-\mathrm{i}K{R}^{\prime}}{\mathrm{e}}^{\mathrm{i}k{r}^{\prime}}\langle {R}^{\prime}N\vert {r}^{\prime}{m}^{\prime}(M)\rangle \nonumber\\ &&\quad \quad =~\mathop{\sum }\limits _{{r}^{\prime}{R}^{\prime}}{\mathrm{e}}^{-\mathrm{i}(k-G){R}^{\prime}}{\mathrm{e}}^{\mathrm{i}k{r}^{\prime}}\langle {R}^{\prime}N\vert {r}^{\prime}{m}^{\prime}(M)\rangle \nonumber\\ &&\quad \quad =~\mathop{\sum }\limits _{{R}^{\prime}}{\mathrm{e}}^{\mathrm{i}G{R}^{\prime}}\mathop{\sum }\limits _{{r}^{\prime}}{\mathrm{e}}^{\mathrm{i}k({r}^{\prime}-{R}^{\prime})}\langle 0 N\vert ({r}^{\prime}-{R}^{\prime}){m}^{\prime}(M)\rangle \nonumber\\ &&\quad \quad =~L\mathop{\sum }\limits _{{r}^{\prime}}{\mathrm{e}}^{\mathrm{i}k{r}^{\prime}}\langle 0 N\vert {r}^{\prime}{m}^{\prime}(M)\rangle , \end{eqnarray} \tag{ 23 }$$

we finally obtain

$$\begin{eqnarray} &&{A}_{k j,k j}(\omega )=\frac{L}{l}\mathop{\sum }\limits _{K G}{\delta }_{k-G,K}{W}_{K J}{A}_{K J,K J}(\omega ) \end{eqnarray} \tag{ 24 }$$

with

$$\begin{eqnarray} &&{W}_{K J}=\mathop{\sum }\limits _{M N r}{\mathrm{e}}^{\mathrm{i}k(r-{r}^{\prime}(M))}C_{M}^{K J}{{C}_{N}^{K J}}^{\ast }{S}_{0 N,r m(M)}, \end{eqnarray} \tag{ 25 }$$

where A_KJ,KJ(ω) is just a delta function at the eigenvalue, δ(ω − _KJ). In equation (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 the 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. The proper spectral weight in the desired representation is given by the guidance of the conceptual system. The connection between A_kj,kj and its counterpart A_{(k−G)J,(k−G)J} can be further discussed. First, it can be found in equation (22) that the spectral function at k is contributed only from the folded counterpart A_KJ,KJ with k − G = K. Secondly, j can be defined as one of the unfolded counterparts of J without ambiguity. It can be confirmed that the sum of weights of the unfolded counterparts for each supercell eigenstate |KJ〉 is exactly one, which further advocates that the two approximations are proper.

The formula of equation (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 in determining the 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 equation (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.

3.5. Orbital contribution

With the obtained unfolding formula, it is also 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 J M}^{k}$ as follows:

$$\begin{eqnarray} {W}_{K J M}^{k}&=&\frac{L}{l}\mathop{\sum }\limits _{G}{\delta }_{k-G,K}C_{M}^{K J}\mathop{\sum }\limits _{r N}{\mathrm{e}}^{\mathrm{i}k(r-{r}^{\prime}(M))}\nonumber\\ &&\times \, {C_{N}^{K J}}^{\ast }{S}_{0 N,r m(M)}. \end{eqnarray} \tag{ 26 }$$

Then, the spectral function can be written as

$$\begin{eqnarray} &&{A}_{k j,k j}(\omega )=\mathop{\sum }\limits _{K}{A}_{K J,K J}(\omega )\mathop{\sum }\limits _{M}W_{K J M}^{k}. \end{eqnarray} \tag{ 27 }$$

By decomposing the spectral weight with equation (26) into the LCAO basis functions, one can analyse 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 for the relabelling for r'(M) and m(M) in equation (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 equation (25) or (26) with the LCAO coefficients and the overlap matrix elements calculated at the second step.

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 the studied surface, where heavily folded bands are involved. The large supercell is especially needed due to a spontaneous reconstruction [31] or an interaction with an added layer [27, 30, 32]. The unfolding method can provide a means of direct comparison of the band structures in the same BZ, which enables the visualization of the symmetry breaking and the surface states. In this section, we demonstrate how our unfolding method can be applied to study the electronic states of a ZrB₂ slab, consisting of 17 Zr and 16 B layers, with a Zr-terminated (0001) surface, where the unit cell of bulk ZrB₂ is used as the normal unit cell. Although the electronic states of the Zr-terminated ZrB₂(0001) surface have been well studied [33, 34], we are motivated by a recent fabrication of silicene on a ZrB₂(0001) surface to reinvestigate the surface states in detail [27].

The DFT calculations were performed by the OpenMX code which is based on norm-conserving pseudopotentials generated with multiple reference energies [35] and linear combination of optimized pseudoatomic basis functions [23, 36]. A 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. A regular mesh of 270 Ryd in real space was used for the numerical integrations and solution of the Poisson equation [20]. The exchange–correlation energy functional was treated with the generalized gradient approximation (GGA) proposed by Perdew et al [37]. Geometrical structures were optimized until the maximum force on an atom became less than 3 × 10⁻⁴ Hartree Bohr⁻¹. In order to validate the accuracy of the DFT method for the system, we optimized the lattice constants of bulk ZrB₂, where the primitive cell of bulk ZrB₂ 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₂ were found to be a = 3.174 Å and c = 3.550 Å with a k-point sampling of an 8 × 8 × 5 mesh, which are in good agreement with experimental data (a = 3.170 Å and c = 3.533 Å) taken from [38]. The geometrical structure and the band structure of bulk ZrB₂ are shown in figure 2.

Zoom In Zoom Out Reset image size

As a next step, we performed geometry optimization of the ZrB₂ slab with a k-point sampling of an 8 × 8 × 1 mesh, where the optimized lattice constants obtained by the bulk calculation were used for the a- and b-axes, and a vacuum of 14.292 Å was inserted to avoid interaction between the slabs. Since the unfolding formula can be applied to the calculation with or without the addition of basis functions in the vacant region, no additional basis functions are added into the vacuum in this study and one can think that the corresponding LCAO coefficients and the overlap matrix elements are zero. The optimized structure exhibits that there is no surface reconstruction for the ZrB₂(0001) surface, which has been checked with the (2 × 2) unit cell of the ZrB₂(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 have 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 figure 2(b), where the LCAO coefficients are obtained. An additional effort in evaluating equation (25) or (26) is to assign the supercell basis functions to the normal cell units. One example is provided for the relabellings, as shown in figure 2(d). Note that only the differences between lattice vectors are needed in the implementation.

Before showing the unfolded band structure, we remark on several differences of our unfolding method for understanding the surface state from conventional study which shows only the atomic weight of the folded band structure. Firstly, 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 makes the visualization of the effect of symmetry breaking straightforward. One more interesting advantage is the capability of recovering the dispersion along the out-of-plane direction, while it is almost impossible for a 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 figure 3. The result clearly shows that the overall bulk band structure is well maintained, accompanied by 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 (cf figure 2(c)); for example, the bands labelled by S1 and S2 which consist of large contributions from the outermost Zr d orbital. The electrons of the S1 band show diminishing spectral weight while the path approaches the Γ point, indicating negligible quasiparticle lifetimes and their short mean free paths. Note that the outermost atomic contribution would show a diminishing behaviour while approaching the Γ point for the S1 band without unfolding [33]. However, our unfolding result demonstrates that a negligible spectral weight is expected even with consideration of the contribution of all the layers.

Zoom In Zoom Out Reset image size

For the dispersion in the out-of-plane direction, horizontal bands stemming 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 [39, 40]. Our unfolding result shows that the discrete energy levels reveal a dispersive behaviour, as found in figure 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 degree. 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 centre of the reference bulk spectral weight (cf figure 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 figure 3 are easily obtained by the unfolding method. Note that the quantum confinement effect has already been observed in multilayer graphene by ARPES [28].

By making use of a conceptual normal cell, we have derived a general formula for unfolding first-principles band structures, which can be applied to the LCAO basis functions as an ideal basis set for the unfolding concept. The translational symmetry breaking is recorded in the LCAO coefficients of the supercell eigenstates and built into 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 the 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 study of the electronic states of a ZrB₂(0001) slab along both the in-plane and the out-of-plane 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 behaviour should be found in all slabs in general, especially 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.

This work was supported by the Strategic Programs for Innovative Research (SPIRE), MEXT, the Computational Materials Science Initiative (CMSI), and Materials Design through Computics: Complex Correlation and Non-Equilibrium Dynamics, A Grant in Aid for Scientific Research on Innovative Areas, MEXT, Japan. This work was also supported by the Funding Program for Next Generation World-Leading Researchers (GR046). The benchmark calculations were performed using the Cray XC30 machine at Japan Advanced Institute of Science and Technology (JAIST).