Visualizing influence of point defects on electronic band structure of graphene
VVisualizing influence of point defects on electronic band structure of graphene
M. Farjam
School of Nano-Science, Institute for Research in Fundamental Sciences (IPM), P.O. Box 19395-5531, Tehran, Iran
The supercell approach enables us to treat the electronic structure of defective crystals, butthe calculated energy bands are too complicated to understand or to compare with angle-resolvedphotoemission spectra because of inevitable zone folding. We discuss how to visualize supercellband structures more effectively by incorporating in them unfolded spectral weights and orbitaldecompositions. We then apply these ideas to gain a better understanding of the band structure ofgraphene containing various types of points defects, including nitrogen impurity, hydrogen adsorbate,and vacancy defect, and also the Stone-Wales defect.
I. INTRODUCTION
The electronic band structure of a crystal is an impor-tant tool to understand its electronic, optical, transportand magnetic properties. Applied to graphene, the bandstructure already shows us many interesting facts aboutits electronic properties.
We can see that graphene isa gapless semiconductor, and its low-energy excitationshave linear dispersion, which implies that they have theproperties of massless Dirac fermions. Defects and impurities have a significant influence onthe electronic properties of semiconductors, and can beintroduced deliberately to tailor their electronic struc-ture. In graphene, atomic impurities, point and struc-tural defects, edges and substrates can all modify theelectronic structure in important ways. Two simple im-purities are substitutional dopants formed by carbonneighbors in the periodic table, boron and nitrogen,which can turn graphene into a p- or n-type semiconduc-tor, respectively. Two other simple defects of a differenttype are vacancies and hydrogen adsorbates, which in-duce midgap states and magnetic moments in graphene. A small structural defect is the 57-57 defect known asStone-Wales (SW) defect.
A convenient way of treating point defects is via the su-percell approach within density-functional theory (DFT)and tight-binding method calculations. The supercelldevice allows us to use methods designed for periodicstructures to treat relatively isolated defects. The work-around has a price: the Brillouin zone (BZ) becomessmaller as the supercell gets larger, which results in thefolding of the bands. The complicated bands are difficultto comprehend or to compare with angle-resolved photoe-mission spectroscopy (ARPES) experiments. This oftenleaves the density of states (DOS) as the sole option forspectral analysis.However, there are ways to improve visualization ofband structures. A simplification is made possible by theunfolding method, which allows plotting effective bandsin the larger BZ of the normal system, i.e., one describedin terms of the primitive unit cell.
Furthermore, anenhancement can be made by incorporating orbital con-tributions in both supercell and normal cell band struc-tures. Orbital decomposition is more commonly usedin obtaining partial density of states (PDOS), but it is sometimes seen in band structure plots as well. In thispaper, we develop ways of implementing the above ideasin the visualization of electronic structure data. We thenapply them to the case of graphene containing point de-fects. Our results are based on density-functional basedtight-binding (DFTB) calculations for 5 × II. COMPUTATIONAL DETAILSA. DFTB calculations
We calculate the prerequisite electronic structure datafor the unfolding method by the DFTB+ code, whichis an efficient and accurate implementation of the DFTBmethod. Since the method is based on a small set ofnon-orthogonal atomic orbitals, 2 s , 2 p x , 2 p y , and 2 p z ,for carbon and nitrogen, and 2 s for hydrogen, it is alsoquite suitable for the unfolding formula to be describedbelow.We use n × n supercells, which encompass N u = n normal unit cells, or 2 n honeycomb lattice sites. Thisimplies for n = 5 a 2% concentration of periodically ar-ranged point defects, which is sufficiently small to con-sider them as relatively isolated, but large enough to pro-duce visible effects. We need to run the code severaltimes. Having defined a supercell with our choice of de-fect placed in a tentative position, in a first calculationwe allow the atomic positions to relax. Then using theoptimized coordinates, we calculate the DOS and PDOSusing a finer k -mesh for the desired accuracy. Finally, weuse the well-converged charges saved from the previousself-consistent calculation to obtain the data for the bandstructure along the given paths in supercell and normalcell Brillouin zones, respectively.The basic ingredients needed for the unfolding methodconsist of eigenvector coefficients and overlap integrals.Thus concomittant with the band structure calculations, a r X i v : . [ c ond - m a t . m t r l - s c i ] M a r FIG. 1. (a) Band structure of graphene obtained from super-cell calculations. (b) Unfolded bands in the normal Brillouinzone. The light gray bands are supercell band structure inextended zone. (c) Total DOS. (d) Color map used in bandstructure plots. (e) The 5 × E F ± we instruct the code to save these data in designated files.To estimate the size of the data, which may grow quitelarge, let N b be the total number of orbitals in a supercell,which is also the total number of bands, and let N k bethe number of k points used in the calculation of theband structure. The overlap integrals are represented asa sparse matrix, so their size can only be estimated tobe less than N R × N b × N b × (size of a real number), where N R is the few number of unit cells required to includenonzero overlaps. On the other hand, the total size ofthe complex eigenvectors, which is by far the larger set ofdata, is exactly 2 × N k × N b × N b × (size of a real number). B. Unfolding procedure
We need an unfolding formula that takes into accountnon-orthogonality of atomic orbitals. Such a formula hasbeen derived recently by Lee et al. in Ref. 11. In thissection, we describe their method from the viewpoint ofnumerical implementation.It is customary to use upper and lower cases to refer to supercell and normal cell variables, respectively, so K and k denote the wavevectors of the supercell Brillouinzone (SBZ) and normal cell Brillouin zone (NBZ), re-spectively. A given k in NBZ folds onto exactly one K inSBZ but, in contrast, each K unfolds back to N u normalsystem k vectors in NBZ. The unfolding method assignsthe same energy at a given K to the multiple unfolded k ’s with proper spectral weights. In practice this simpleprocedure can be followed. If we choose a k path in theNBZ, and regard it as a path in the extended zone of thesupercell, then each K needs to be unfolded to only a sin-gle k = K in order to achieve the desired band structurerepresentation. We just calculate the energy eigenvaluesalong this path for the supercell sytem, and determinetheir spectral weights by the unfolding formula.It is useful to note some properties of non-orthogonalorbitals. The normalization of an eigenvector, | K J (cid:105) , interms of its expansion coefficients, is given by (cid:88) MN C K J ∗ N S NM ( K ) C K JM = 1 , (1)where S ( K ) is the overlap matrix of Bloch orbitals, S NM ( K ) = (cid:88) R e i K · R S NM ( R ) . (2)Here the overlap integrals are defined by S NM ( R ) ≡ S N, R M = (cid:104) N | R M (cid:105) . (3)In Eq. (1), J and N (or M ) are band and orbital indices,respectively, and are both in the range of 1 to N b . Break-ing down the contributions in Eq. (1), we can write theorbital populations for an eigenstate as P K JN = Re (cid:40) C K J ∗ N (cid:88) M S NM ( K ) C K JM (cid:41) . (4)The unfolding procedure must uncover the hiddentranslational symmetry inherited by the supercell systemfrom the normal system. A useful connection between thetwo systems is expressed by the Fourier relation, N u (cid:88) G e i G · r = δ r , R , (5)where R and r are the corresponding lattice vectors, re-spectively, and the sum is over the set of N u reciprocallattice vectors of the supercell system that unfold theSBZ onto the NBZ. An orbital described in the supercellsystem by the pair R , M can be described in the normalsystem through the mapping, R , M → R + r M , m M , (6)where r M is one of the N u normal system lattice vectorswithin a supercell. The range of m M is the (possiblyvariable) number of orbitals in each normal unit cell.The unfolded spectral weight is given by W K J ( G ) = 1 N u (cid:88) MN C K J ∗ N U NM ( k ) C K JM , (7)where k = K + G , and U NM ( k ) = (cid:88) r e i k · ( r − r M ) S N, r m M , (8)which is a particular Fourier sum of overlap integrals thatencapsulates the unfolding information. We shall call itthe unfolding matrix . A convenience of this unfoldingformula is that it does not require the definition of avirtual crystal for its implementation.Noting the similarity between U ( k ) and S ( K ), we writepartial unfolded orbital weights as W K JN ( G ) = 1 N u Re (cid:40) C K J ∗ N (cid:88) M U NM ( k ) C K JM (cid:41) , (9)so that W K J ( G ) = (cid:80) N W K JN ( G ). We can verify thatthe unfolding matrix is hermitean, so that the spectralweights are properly real. This can be shown by writing S N, r m M as S r N n N , r m M , and using the symmetry proper-ties of overlap integrals. Another checkpoint is the spe-cial case of an orthonormal basis set. We can easily seethat orthogonality implies that S N, r m M = δ r , r N δ n N ,m M ,from which we obtain the simpler result, U NM ( k ) = e i k · ( r N − r M ) δ n N ,m M . (10)The main numerical task pertains to the unfolding ma-trix in Eq. (8). First, there is a sum over the infinite set oflattice vectors r , but the actual set is sharply restrictedby the range of overlap integrals. Second, it must benoted that the overlap integral S N, r m M is equivalent toone in the standard form of S N, R M (cid:48) , supplied by theelectronic structure code. The supercell and normal celllattice vectors are related by r = R + r i , where r i is anormal lattice vector within a supercell, while the orbitalindex M (cid:48) is determined by the pair of r , m M .Two sum rules can help us verify the results of ourcalculations. By using Eq. 5, we can prove the first one, (cid:88) G W K J ( G ) = 1 . (11)The second one, which is more practical for our proce-dure, is given by (cid:88) J W K J ( G ) = N b N u , (12)where the right hand side is just the average numberof basis orbitals in the primitive unit cell of the normalsystem. This sum rule is a consequence of the propertiesof the spectral function. FIG. 2. (a) Band structure of nitrogen-doped graphene ob-tained from supercell calculations. (b) Unfolded bands in thenormal Brillouin zone. (c) Total DOS and PDOS of nitrogenorbitals. (d) Color map used in band structure plots. (e)The 5 × III. RESULTS AND DISCUSSIONA. Perfect graphene
Band structure of perfect graphene provides a goodtesting ground for our numerical procedure, since weknow exactly what to expect. Furthermore, it allowsus to measure the computational load of the unfoldingprocedure, which is about the same in all cases that weconsider. Using a Fortran serial code, on a computer with
Intel R (cid:13) Core TM i7-3610QM CPU, the runtime for the un-folding of this example was about one minute.We first explain our approach to using additional dataconsisting of unfolded spectral weights and orbital popu-lations in the visualization of band structures. (A some-what different approach is described in the
SupplementalMaterial of Ref. 16.) In the ( matplotlib ) 2D plotting pack-age, used here, we can specify colors and transparency interms of parameters r, g, b and α , respectively, all in the0–1 range. Hence in the rgb color space (0,0,0) and(1,1,1) represent black and white, respectively, which arevisible or invisible on a white background. The α param-eter controls transparency, so that α = 0 is completelytransparent and thus invisible, and α = 1 is completelyopaque. We adopt a simple scheme and always decom-pose the orbital contributions into only two sets, whichare represented by blue (0,0,1) and red (1,0,0), respec-tively. We then use the normalized weighted average ofthe colors to represent mixing of the two sets, and set α equal to the total spectral weight, which results in thecolor map shown in Fig. 1(d). It is possible to partitionthe orbitals into more than two sets but the color mapsbecome more complicated. We must mention that thiscontinuous range of colors is unnecessary for the presentcase of perfect graphene, where the weights turn out tobe zero or unity, but it will be essential for defectivegraphene.Figure 1(a) shows the band structure, using 25 k pointsalong ΓKM path in the SBZ. We have assigned our twopure colors to the two sets consisting of 2 s , 2 p x , 2 p y orbitals and the 2 p z orbital, respectively. This has madethe σ and π bands to appear in these colors, since thesebands are mutually orthogonal in graphene. However,heavy folding makes it hard to recognize the shape of thebands.Figure 1(b) shows the unfolded bands. We make thecalculations using the same 5 × k points. The ex-pected bands of graphene have emerged, in blue and red,out of the complicated supercell bands, shown in lightgrey. Figures 1(f–h) show the same data as in Figs. 1(a–c), but in the important range of low energies. B. Nitrogen impurity
We repeated similar electronic structure calculationsfor a 5 × σ and π bands, because of the pres-ence of nitrogen impurity. In addition, unfolding showsthat the weights of these states have nonuniform distri-butions in the NBZ.In Figs. 2(f–h) we show the data for the low energyspectrum. Because of doping by nitrogen, the Fermi levelhas shifted above the neutrality point. The presence ofimpurity has also caused a small gap to open at the Diracenergy.
FIG. 3. (a) Band structure of graphene with a hydrogenadsorbate obtained from supercell calculations. (b) Unfoldedbands in the normal Brillouin zone. (c) Total DOS and PDOSof hydrogen orbital. (d) Color map used in band structureplots. (e) The 5 × C. Hydrogen adsorbate
Results of calculations for a single hydrogen adsorbateon a 5 × sp defect. Such adefect effectively creates a π orbital vacancy, which man-ifests itself as a midgap state at the Fermi level. Wecan see the midgap state in all the band structure plotsas well as the DOS plots. Moreover, most of the weight ofthe midgap band comes from carbon 2 p z orbitals. Thereis a bound state visible in the gap above the π bandswith contributions from both C and H atoms, Figs. 3(a–c). We note, in Fig. 3(g), that the weight of the midgapband is accumulated near the K point, a fact that hasalso been observed in ARPES spectra of hydrogenatedgraphene. FIG. 4. (a) Band structure of graphene with a vacancyobtained from supercell calculations. (b) Unfolded bands innormal Brillouin zone. (c) Total DOS and the PDOS of threecarbon atoms bordering the vacancy site. (d) Color map usedin band structure plots. (e) The 5 × D. Vacancy defect
Electronic structure calculation for graphene with amissing atom is shown in Fig. 4. Although not visiblein Fig. 4(e), the relaxation of atoms results in a Jahn-Teller distortion. There are four missing orbitals at thevacancy site, which can give rise to midgap bands, andother spectral changes. Three midgap bands are clearlydisplayed in Figs. 4(f–h), where the use of the color maphas also revealed the difference in their orbital decompo-sition. The two pure colors were assigned to two sets ofcarbon atoms, those nearest the vacancy site and thoseon the remaining sites.
E. Stone-Wales defect
Lastly, we performed calculations for the SW defect.A similar example was recently treated with a differentband unfolding methodology in Ref. 24, which could pro-vide an interesting comparison. Although the Dirac coneis generally preserved as pointed out in Ref. 24, a small
FIG. 5. (a) Band structure of graphene with a Stone-Walesdefect obtained from supercell calculations. (b) Unfoldedbands in normal Brillouin zone. (c) Total DOS. (d) Colormap used in band structure plots. (e) The 5 × gap has appeared at the K point of the unfolded bands,shown in Fig. 5(g). A close examination of the foldedbands, Fig. 5(f), shows that there is actually no gap inthe energy spectrum, but the band crossing has slightlyshifted from its SBZ K point toward the Γ point. Also,the SW defect has caused a lot of smearing and breakingup of the bands, and appearance of new bands, especiallyfor the σ bands. It can be expected that the magnitudeof the defect-induced effects will diminish smoothly asthe size of the supercell is gradually increased.In all cases, we checked that the sum rule, (12), holds.For pure graphene, SW defect or nitrogen impurity, thesum adds up to N b /N u = 200 /
25 = 8, while for hydrogenand vacancy defects it adds up to 201 /
25 = 8 .
04, and196 /
25 = 7 .
84, respectively. The variety of the examplesused show the versatility of our procedure.
IV. CONCLUSIONS
We have implemented the unfolding method of Lee et al. , have proposed how to visualize additional datafrom spectral weights in energy dispersion plots, and havedemonstrated these for selected defective 2D graphenes.Although the DFTB+ code was used in this work, theunfolding procedure described can more generally be em-ployed for any electronic structure code that uses a basisof non-orthogonal atomic orbitals. It can be predictedthat a great number of graphene-based systems, involvingdefects, impurities, adsorbates and substrates can takeadvantage of the methods discussed here. The simpli-fied effective band structures could be more easily com-prehended, and supercell calculations may be compared directly with normal cell calculations. The major attrac-tion, however, is the correlation of the unfolded bandstructure with the measurements of angle-resolved pho-toemission spectroscopy (ARPES) experiments. ACKNOWLEDGMENTS
I have greatly benefited from discussions withAlex Gr¨uneis. A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S.Novoselov, and A. K. Geim, Rev. Mod. Phys. , 109(2009). P. R. Wallace, Phys. Rev. , 622 (1947). K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang,Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A.Firsov, Science , 666 (2004); K. S. Novoselov, A. K.Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V.Grigorieva, S. V. Dubonos, and A. A. Firsov, Nature ,197 (2005); Y. Zhang, Y. W. Tan, H. L. Stormer, andP. Kim, ibid . , 201 (2005). D. Usachov, O. Vilkov, A. Gr¨uneis, D. Haberer, A. Fe-dorov, V. K. Adamchuk, A. B. Preobrajenski, P. Dudin,A. Barinov, M. Oehzelt, C. Laubschat, and D. V. Vyalikh,Nano Lett. , 5401 (2011). O. V. Yazyev and L. Helm, Phys. Rev. B , 125408 (2007). A. Stone and D. Wales, Chem. Phys. Lett. , 501 (1986). F. Banhart, J. Kotakoski, and A. V. Krasheninnikov, ACSNano , 26 (2011). P. B. Allen, T. Berlijn, D. A. Casavant, and J. M. Soler,Phys. Rev. B , 085322 (2013). T. B. Boykin and G. Klimeck, Phys. Rev. B , 115215(2005); T. B. Boykin, N. Kharche, G. Klimeck, andM. Korkusinski, J. Phys.: Condens. Matter , 036203(2007). W. Ku, T. Berlijn, and C.-C. Lee, Phys. Rev. Lett. ,2010 (2010). C.-C. Lee, Y. Yamada-Takamura, and T. Ozaki, J. Phys.:Condens. Matter , 345501 (2013). V. Popescu and A. Zunger, Phys. Rev. Lett. , 236403(2010); Phys. Rev. B , 085201 (2012). A. N. Rudenko, F. J. Keil, M. I. Katsnelson, and A. I.Lichtenstein, Phys. Rev. B , 081405(R) (2013). B. Aradi, B. Hourahine, and T. Frauenheim, J. Phys.Chem. A , 5678 (2007). We used the mio-0-1 setof Slater-Koster parameters, obtained from M. Elstner, D. Porezag, G. Jungnickel, J. Elsner,T. Frauenheim, S. Suhai, and G. Seifert, Phys. Rev. B , 7260 (1998). T. Berlijn, D. Volja, and W. Ku, Phys. Rev. Lett. ,077005 (2011). J. D. Hunter, Computing in Science and Engineering , 90(2007). Ph. Lambin, H. Amara, F. Ducastelle, and L. Henrard,Phys. Rev. B , 2012 (2012). T. O. Wehling, M. I. Katsnelson, and A. I. Lichtenstein,Phys. Rev. B , 085428 (2009). M. Mirzadeh and M. Farjam, J. Phys.: Condens. Matter , 235304 (2012). M. Farjam, D. Haberer, and A. Gr¨uneis, Phys. Rev. B. , 193411 (2011). D. Haberer, L. Petaccia, M. Farjam, S. Taioli, S. A. Ja-fari, A. Nefedov, W. Zhang, L. Calliari, G. Scarduelli,B. Dora, D. V. Vyalikh, T. Pichler, C. W¨oll, D. Alf`e, S. Si-monucci, M. S. Dresselhaus, M. Knupfer, B. B¨uchner, andA. Gr¨uneis, Phys. Rev. B , 165433 (2011). B. R. K. Nanda, M. Sherafati, Z. S. Popovi´c, and S. Sat-pathy, New J. Phys. , 083004 (2012). P. V. Medeiros, S. Stafstr¨om, and J. Bj¨ork, Phys. Rev. B89