Electronic structure of the substitutional vacancy in graphene: Density-functional and Green's function studies

We study the electronic structure of graphene with a single substitutional vacancy using a combination of the density-functional, tight-binding, and impurity Green's function approaches. Density functional studies are performed with the all-electron spin-polarized linear augmented plane wave (LAPW) method. The three $sp^2 \sigma$ dangling bonds adjacent to the vacancy introduce localized states (V$\sigma$) in the mid-gap region, which split due to the crystal field and a Jahn-Teller distortion, while the $p_z \pi$ states introduce a sharp resonance state (V$\pi$) in the band structure. For a planar structure, symmetry strictly forbids hybridization between the $\sigma$ and the $\pi$ states, so that these bands are clearly identifiable in the calculated band structure. As for the magnetic moment of the vacancy, the Hund's-rule coupling aligns the spins of the four localized V$\sigma_1 \uparrow \downarrow$, V$\sigma_2 \uparrow $, and the V$\pi \uparrow$ electrons resulting in a S=1 state, with a magnetic moment of $2 \mu_B$, which is reduced by about $0.3 \mu_B$ due to the anti-ferromagnetic spin-polarization of the $\pi$ band itinerant states in the vicinity of the vacancy. This results in the net magnetic moment of $1.7 \mu_B$. Using the Lippmann-Schwinger equation, we reproduce the well-known $\sim 1/r$ decay of the localized V$\pi$ wave function with distance and in addition find an interference term coming from the two Dirac points, previously unnoticed in the literature. The long-range nature of the V$\pi$ wave function is a unique feature of the graphene vacancy and we suggest that this may be one of the reasons for the widely varying relaxed structures and magnetic moments reported from the supercell band calculations in the literature.

. Sketch of the electronic structure for an isolated substitutional vacancy in graphene. The continuum π and σ bands are shown as dashed curves, while the vacancy-induced localized states, Vσ and Vπ, are denoted by straight lines. E F is the Fermi energy. The occupied vacancy states are indicated by solid circles with a corresponding net magnetic moment of 2µ B . The circular density-of-states (DOS) in the midgap region, labelled π l ↑↓, indicates schematically the antiferromagnetic spin-polarization of the π electron states in the local neighborhood of the vacancy. This spin polarization is responsible for the reduction of the localized magnetic moment from 2µ B (S = 1) to about 1.7µ B in our densityfunctional calculation.

Introduction
Graphene is a material of considerable interest on account of its unusual linearlydispersive Dirac band structure and particle-hole symmetry. [1,2] Vacancy constitutes an important defect center, the electronic structure of which forms the basic foundation for the understanding of the behavior of more complex defects including impurities.
Recently it has been suggested that transition-metal doped graphene with vacancies may have potential application in hydrogen storage. [3] Experimentally, vacancies in graphene have been created intentionally by irradiating materials with electrons and ions [4,5,6,7] and they may also occur in small concentration during the growth process. [8] While an ideal graphene sheet is non-magnetic, experimental observation of magnetism in carbon systems has been long explained in terms of a variety of defects including isolated vacancies, vacancy clusters, or presence of internal or external boundaries as in nanoribbons. [7,9,10,11] There have been several theoretical studies of the isolated vacancy in graphene from first-principles density-functional theory (DFT) [12,13,14,15,16,17,18,19,20,21,22] or Hartree-Fock calculations [23] as well as from tight-binding models [25,26,27]. There is also an enormous amount of related work on the chemisorbed defects such as the hydrogen defects and chemisorbed magnetic atoms. [28] Most of the tight-binding Electronic structure of the substitutional vacancy in graphene: Density-functional and Green's function studies3 Exchange Jahn− Teller field Figure 2. Splitting of the three dangling bond sp 2 σ states of the carbon triangle, denoted by Vσ, and the vacancy-induced zero-mode Vπ state originating from the π band. The splitting of the Vσ states is discussed in detail in Section 3. models have focused on the π bands only, which is clearly inadequate due to the formation of the sp 2 σ dangling bonds in the mid-gap region. The first-principles calculations do include all relevant states in the band structure including the sp 2 σ states; however, in spite of all these works, a clear picture of the vacancy states has not emerged.
In this paper, we study the electronic structure of the vacancy in graphene using the all-electron density functional linear augmented plane waves (LAPW) method along with tight-binding studies as well as the impurity Green's function (GF) approach to interpret the band structure. To our knowledge, this is the first allelectron density functional calculation for the vacancy in graphene reported in the literature. We have already reported the electronic structure for the mono and bilayer graphenes using the same method. [29] In addition to the DFT calculations, the nature of the vacancy-induced states is modeled from the tight-binding and Green's function studies, which help interpret the DFT results.
The basic overall picture of the electronic structure that emerges from our work is summarized in Fig. 1. It shows the standard σ and π bands of graphene plus the vacancy-induced states, denoted by Vπ and Vσ, which are split due to the crystal field, Jahn-Teller, and the Hund's-rule interactions. The Vσ states are made out of the three sp 2 σ dangling bond states, which are located on the three carbon atoms adjacent to the vacancy with their lobes directed towards the vacancy site. With their bonding partners missing, they occur in the midgap region. At the same time, a localized state Vπ gets introduced in the π bands in the midgap region as well, the so called "zero mode" state, whose energy is exactly zero in the nearest-neighbor tight binding approximation. These four states, localized around the vacancy center, can hold eight electrons in total taking into account the spin degeneracy. The level structure of the vacancy-induced states is shown in Fig. 2.
At the same time, electron counting arguments show that the vacancy releases four electrons to be occupied among the above localized states. These electrons include the three orphan sp 2 σ electrons, one from each of the three carbon atoms adjacent to the vacancy, plus one orphan π electron, whose origin may be understood in the following way. Focusing on the π states now and considering a vacancy on the A sublattice, the majority sublattice B has one extra atom, N B − N A = 1, so that the total number of π orbitals is N A + N B , which is the same as 2N A + 1. Out of these, there is one zero-mode state and the electron-hole symmetry of the graphene lattice results in N A band states below E = 0 and the same number above it. (See Fig. 7 for the π band structure). So, of the 2N A + 1 π electrons (one per atom), 2N A fill up the lower bands, leaving a lone orphan π electron. These four orphan electrons (three σ and one π) occupy the vacancy-induced states as indicated in Fig. 1.
The remaining sections are organized as follows. In Section II, we discuss the results of our DFT calculations. Section III discusses the crystal-field and Jahn-Teller splitting of the vacancy-induced localized σ states and Section IV is devoted to the vacancy-induced π states. In Section IV A, we revisit the zero-mode theorem and in Section IV B, we present numerical results for the π states from a numerical diagonalization of the tight-binding Hamiltonian before discussing the vacancy-induced π states using the Green's function approach in Sections IV C and D. Finally the results are summarized in Section V.

Density-Functional Calculations
For the density-functional calculations, we used the all-electron spin-polarized linear augmented plane wave (LAPW) method [30] with the general gradient approximation (GGA) [31] for the exchange-correlation functional. A 72-atom 6×6 supercell was used which included one vacancy site. The LAPW basis functions included the carbon 2s and 2p valence orbitals and a momentum cutoff of RK max = 5.2 was used, with approximately 3500 basis functions and about 50,000 plane waves at each k point. All atomic sphere radii were taken as 0.66Å. The maximum angular momentum for the wave function expansion inside the atomic sphere was kept at l max = 6. Thirty k points in the irreducible Brillouin zone were found to be sufficient for converged results in the self-consistent calculations.
Relaxed structure -First we performed a structural optimization of the lattice constant for pure graphene which yielded about the same lattice constant as the experimental value. For the vacancy calculation, the lattice constant was held fixed at the experimental value and a structural relaxation was performed for the entire structure. The optimization yielded a planar Jahn-Teller (JT) distorted carbon triangle around the vacancy with the carbon atoms outside the triangle relaxed by a much smaller amount. For the carbon triangle, we found two long bonds of length 2.66Å each and a short bond of length 2.40Å (Fig. 3), as compared to the 2.48 A for the undistorted structure. In terms of the standard Jahn-Teller modes of the equilateral triangle, the magnitudes of the distortion are: Q 0 = 0.08Å (breathing mode), Q 1 = 0.166Å (symmetric bond-bending mode), and Q 2 = 0 (asymmetric mode). [32] The relaxed structure for the vacancy reported in the literature varies widely. While some have reported planar structures, [13,17,15,16,23,14] others have found non-planar structures with out-of-plane displacements varying from δz ≈ 0.12 − 0.47 A. [12,19,20,21,22] We find that a paramagnetic relaxation (less accurate for the present problem) yields a non-planar structure δz ≈ 0.27Å, while a spin-polarized calculation yields a planar structure, an observation made by Faccio et al. [17] from their calculations as well using the SIESTA code. We attribute this wide variation in the calculated relaxed structure in the literature partly to the unusual nature of the Vπ bound state, which falls off only as 1/r, leading to a larger width of the Vπ band  Figure 3. The Jahn-Teller distorted planar carbon triangle obtained from the structural relaxation using the all-electron spin-polarized LAPW-GGA method.
in the supercell calculations than is expected from resonance broadening due to the π band continuum.
The calculated vacancy formation energies agree much better between different calculations. Our result for E V = E(graphene + vacancy) − N −1 (N − 1)E(graphene), N being the number of atoms in the graphene supercell, is 7.87 eV, which compares well with the previous calculations [12,20,15] of 7.4 − 7.8 eV as well as with the experimental value of 7.0 ± 0.5 eV. [33] Electronic structure - Fig. 4 shows the band structure, where the vacancy induced Vσ and Vπ states are clearly seen. The momentum points in the Brillouin zone for the band structure plot are defined as: K =x/ √ 3 +ŷ and M =ŷ in units of 2π3 −1 a −1 /n with n = 6 for the 6 × 6 supercell used in the calculation. For this supercell, it can be easily seen by drawing both Brillouin zones that the Dirac points K and K ′ of graphene get folded into the Γ point of the supercell Brillouin zone, so that remnants of the Dirac bands are seen at the Γ point in Fig. 4 just above E F (See also Fig.  7 for the folded graphene tight-binding π bands for the same supercell and note the similarity between the tight-binding π bands and the DFT bands, Fig. 4). Due to symmetry, the σ and π bands don't mix (strictly forbidden for the planar geometry, but also weakened significantly if the relaxed geometry is non-planar), which leads to clearly identifiable vacancy-induced Vσ bands. The Vσ states originating from the dangling bonds are split due to the crystal field, Jahn-Teller, and exchange coupling as indicated in Fig. 2 and discussed in more detail in Section 3. The dispersion of the Vσ bands in the band structure comes from the vacancy-vacancy interactions in different supercells or from the k-dependent interaction with the bonding and the anti-bonding σ bands, both effects being small. For non-planar relaxed structure, they should have a small resonance broadening due to the interaction with the π band continuum. Three electrons occupy these states leading to the occupation Vσ 1 ↑↓, Vσ 2 ↑, with the remaining fourth electron occupying the Vπ ↑ state.
We now turn to a description of the effect of the vacancy on the π states. Basically, the vacancy introduces a sharp, resonance state Vπ in the midgap region. The following summarizes the discussions in Section 4, which are important to keep in mind here: (i) If only NN tight-binding hoppings are kept, then the vacancy introduces a single localized state Vπ at E = 0 and of zero width called the zero-mode state, and its wave function decays as ∼ 1/r with distance in the linear-band approximation.  supercell, all vacancies are located on the same, minority sublattice). (iii) However, due to the 2NN hopping as well as the exchange splitting, the energy of Vπ is different from zero, so that it now has a small but finite width due to resonance broadening with the linear π band continuum consistent with the STM experiments. [10] (iv) In the supercell calculations, the Vπ state acquires an extra broadening due to the slow 1/r decay of the Vπ wave function, because of the interaction between the supercells. Dirac point -In Fig. 4, the Dirac point occurs above the E F (see the bands just above E F at the Γ point, to which the standard Dirac points K and K ′ have become folded). For the truly isolated vacancy, the location of the Dirac point above E F would mean that an infinite number of electrons are transferred from the unfilled part of the Dirac cones to the lone vacancy site, which is impossible. Another way of seeing this is to consider first an infinite lattice without the vacancy. Obviously, the E F occurs at the Dirac point with zero density-of-states as usual. Now, if we introduce a single vacancy into the structure it can only affect the position of E F by ∼ 1/N , where N is the total number of atoms in the lattice, so that E F remains unchanged for the infinite lattice. Of course, the electron states in the local neighborhood of the vacancy are modified, e. g., due to the resonance interaction with the vacancy states or due to the vacancy potential. The Dirac-like bands seen just above E F at Γ in Fig. 4 represent the effect of the vacancy on the electronic structure in the local neighborhood of the vacancy in the supercell calculation.  Figure 5. Sketch of the magnetic moment for an isolated vacancy, emphasizing the spatial extent of the various electronic states. The Vσ electrons are highly localized on the carbon triangle surrounding the vacancy, while the Vπ electron is only "quasi-localized" with its wave function decaying only as 1/r. Hund'srule exchange aligns the Vσ and Vπ electrons, producing a S = 1 state with the nominal magnetic moment of 2µ B . This moment is however reduced by polarization of the π band spins in the vicinity of the vacancy, described by an antiferromagnetic Kondo-like coupling t kπ between the π bands and the localized Vπ and Vσ moments. The π band polarization is about 0.3µ B in our DFT calculations, leading to the net magnetic moment of 1.7µ B .
Magnetic moment -The vacancy magnetic moment consists of two parts as shown schematically in Fig. 5: (i) the localized moment coming from the vacancy states Vπ and Vσ and (ii) the induced moment on the band electrons in the vicinity of the vacancy. One can argue on general grounds that the first contribution should be 2µ B (S = 1), while the second contribution should reduce this value somewhat due to the antiferromagnetic Kondo-like coupling between the localized and the itinerant band spins. Turning to the localized states, the vacancy leaves four electrons to be occupied among the Vσ dangling bond states and the Vπ zero-mode state. Of these, three electrons will occupy the Vσ states, so that one electron resides on each of the three dangling bonds of the carbon triangle. The Coulomb interaction U would prevent the occupation of a fourth Vσ state, so that the remaining electron is energetically favored to occupy the π states. The Hund's coupling between the Vσ and Vπ electrons leads then to a S = 1 state with a magnetic moment of 2µ B . This basic picture is illustrated in Figs. 2 and 5 and it is fully supported by the DFT bands, Fig. 4. This localized magnetic moment of 2µ B is reduced due to the spin-polarization of the π bands in the vicinity of the vacancy.
The spin polarization of the π bands can occur due to two factors: (i) the resonance coupling with the Vπ ↑ electron with the π continuum bands and (ii) the Kondo-like antiferromagnetic interaction between the localized vacancy states and the continuum π states. The first is not well described in a supercell calculation due to the long-range nature of the Vπ state and the second effect is intrinsically not well described within the band theory.
Our DFT calculations yield a magnetic moment of about 1.7µ B . This can be seen by estimating the number of holes in the small hole pocket in the two bands just above E F at the Γ point in the spin-up bands of Fig. 4. The spin-down bands must contain exactly the same number of extra electrons missing from the spin-up bands. Without this pocket of holes, which represents the band polarization in the immediate neighborhood of the vacancy, the magnetic moment would be exactly 2µ B , corresponding to the full occupancy of Vσ 1 ↑↓, Vσ 2 ↑, and Vπ ↑. The existence of the hole pocket reduces this number. We can estimate the number n in the hole pocket by computing the total area of the two hole Fermi surfaces and comparing it to the area of the supercell Brillouin zone, which yields the value n ≈ 0.15. Since the same number of electrons must be accommodated in the spin down bands, this would cause a net reduction of N ↑ − N ↓ by 0.30 leading to a net magnetic moment of 1.7µ B .
In the literature, the calculated magnetic moment varies widely, anywhere between 1.04 − 1.84µ B [13,19,18,12,20,15,14,16,21,22,17,23,24]. Typically, the lower values come from calculations, where the Vπ ↑ and Vπ ↓ bands overlap significantly. We suggest that the variation of the calculated magnetic moment in the literature is due to the intrinsic deficiency of the supercell method in estimating the π magnetic moment due to the slow 1/r decay of the Vπ state, which produces an extra broadening of the Vπ state due to the supercell interaction and does not take into account the full anti-ferromagnetic polarization of the itinerant π band states.
The exchange splitting ∆ of the Vπ state is due to its overlap with the Vσ states which are localized on the three carbon atoms adjacent to the vacancy. It may be estimated from the expression where the Hund's-rule energy is typically J H ∼ 0.9 − 1.0 eV for the atoms and |Ψ 0 | 2 ∼ 0.4 is the combined total density of the Vπ state on the carbon triangle as obtained from the DFT results. The estimated exchange splitting is in agreement with the splitting seen in the DFT bands, Fig. 4.
Relation to the Lieb's Theorem -The Lieb's theorem [34] states that for the repulsive one-band Hubbard model on a bipartite lattice and half-filled band, the ground state has spin S = (1/2)|N B − N A |, N A (N B ) being the number of sites on the two sublattices. It is important to point out that the theorem holds if we consider only the π-band system and also neglect the small second-neighbor interactions that couples the two sublattices. Thus, with a single vacancy present, |N B − N A | = 1 so that according to the Lieb's Theorem we should have a net spin of S = 1/2. However, in addition to the π, we also have the σ electrons. The Lieb result of S = 1/2 for the π electrons is now coupled to the spins of the three σ electrons localized near the vacancy, leading to the net spin S = 1 as indicated in the summary figure, Fig. 1. We have already argued that the magnetic moment of 2µ B corresponding to S = 1 will be reduced due to the polarization of the band electrons in the local neighborhood of the vacancy.

Vacancy-induced Vσ states
The essential features of the density-functional results may be understood by simple tight-binding considerations of the effect of the vacancy on the σ and the π bands. We study the σ states in this section followed by the π states in the next section.
The description of the vacancy-induced V σ states for graphene is rather simple. In graphene, the sp 2 σ states are removed away from E F due to strong interaction with neighbouring orbitals along the C-C bonds. However, with a vacancy present, the three sp 2 σ orbitals of the three NN carbon atoms with their lobes pointed towards the vacancy have their usual bonding partners missing, so that they occur near E F , with their on-site energies ǫ σ slightly below the π orbital energies because of the s orbital component present in the σ states.

Electronic structure of the substitutional vacancy in graphene: Density-functional and Green's function studies9
The crystal-field splitting however will lift the three-fold degeneracy. The main feature can be described by taking into account the 2NN hopping T between the three dangling bonds in the undistorted triangle, leading to the 3 × 3 Hamiltonian: diagonalization of which yields a double degenerate state at E = T and a single degenerate state at E = −2T as shown in Fig. 2, where we call this splitting the crystal-field splitting. The Jahn-Teller distortion of the triangle splits the double degenerate state further, which is described by the unequal hopping T = T ′ . Taking the isosceles-triangle relaxation found in our DFT results, two of the three hopping terms are modified into T ′ as indicated in Fig. 6. From the DFT band structure, we find that T ≈ 1.6 eV, while T ′ ≈ 1.2 eV. The new eigenstates are: ). This simple model suggests a Jahn-Teller distortion of the carbon triangle surrounding the vacancy.
The Jahn-Teller interaction is of the type E ⊗ e (both electronic and vibrational states are doubly degenerate) in a trigonal (D 3h ) symmetry. With this lattice distortion, the trigonal symmetry is broken. The distortion removes the double degeneracy and the two states (shown in Fig. 2 Since there are only three electrons available to the Vσ states, V σ 1 is occupied with two electrons, while the lone remaining electron occupies the V σ 2 state. The spin-degeneracy is removed by the Hund's coupling with the electron occupying the Vπ state, producing the spin structure indicated in Fig. 2. The wave function Ψ 2 corresponding to the V σ 2 state shows that the maximum weight (∼ 66%) comes from the sp 2 σ dangling orbital of the apical atom of the carbon triangle, which is consistent with the spin density plotted in Fig. 6. The Jahn-Teller distortion is actually dynamic, with the carbon triangle tunneling between three equivalent minima on the adiabatic potential surface of the E ⊗ e Jahn-Teller problem, an issue we discuss elsewhere. [35]

Vacancy-induced Vπ states
In this Section, we discuss the origin of the localized state -the so-called "zeromode" state -near the band center of the π bands. Understanding of the origin and the "quasi-localized" nature of the zero-mode state is an essential part of the interpretation of the full band calculation using the density-functional theory.
If only the NN interactions are present, the "zero-mode" state is a localized state with energy exactly at the band center. If in addition the higher-neighbour interactions are also present but not too large, as is the case for graphene, [29] then the localized state turns into a sharply-peaked resonance state owing to its overlap with the π bands and occurs not too far from the band center.

The existence of the zero mode state
According to the zero-mode theorem, [25] which is in fact valid for any bipartite lattice with NN interactions, whenever there is an imbalance in the number of atoms in the two sublattices in a bipartite lattice, viz., n = N B − N A > 0, there are n number Figure 6. Spin density n ↑ − n ↓ at different sites in graphene with a vacancy as obtained from the density-functional calculations. Green (blue) circles indicate positive (negative) values, with the area of the circle being proportional to the spin density. The spin moments on the carbon atoms other than the vacancy triangle are due to the π electrons, which are long-ranged due to the slow 1/r decay of the Vπ state. The hopping integrals T and T ′ between the sp 2 σ orbitals on the carbon triangle adjacent to the vacancy has reference to the model discussed in Section 3.
of degenerate solutions with the eigenvalue ǫ B (the on-site energy of the majority sublattice), with the wave functions residing entirely on this sublattice. This can be seen from the following simple considerations as an alternative to Pereira et al.'s proof which used the rank-nullity theorem in linear algebra. [25] We begin with the conjecture that there are some solutions where the wave functions live completely on the majority sublattice (B) and proceed to find them. Thus we have where Ψ B is a vector in the B sublattice of dimension N B and there is null contribution from the A sublattice to the wave function. It will be clear from the following discussion that for the theorem to hold, the B sublattice N B × N B Hamiltonian is restricted to the diagonal form and there are no restrictions on the remaining parts of the Hamiltonian. The specific form of H BB means that there is no site disorder, nor is there any interactions between the atoms on the B sublattice (hence it will fail if interactions beyond the NN are present, which will produce a non-diagonal H BB ). However, such restrictions need not apply to the A sublattice, so that the N A × N A Hamiltonian H AA for the minority sublattice can have diagonal disorder and also there is no restriction on the form of H BA as well. This means that the A sublattice atoms can interact between themselves and with the B sublattice atoms as well without invalidating the theorem. The wave function Ψ B thus satisfies The first of these equations tells us that if the conjectured solutions of the form (Ψ B , 0) exist, then they must have the enegy E = ǫ B and there would be at most N B

Electronic structure of the substitutional vacancy in graphene: Density-functional and Green's function studies11
number of such degenerate solutions; the equation does not place any constraints on the individual components of Ψ B . Turning to Eq. 6, there are N B components of Ψ B to be determined but only N A < N B equations are there to determine them. This means that the solutions can not be fully determined. However if we specify N B − N A components of Ψ B , then the remaining components are uniquely determined as linear functions of these components. These solutions are therefore of the form where we can choose the φ i 's arbitrarily and the f i 's are then just linear combinations c ij φ j , where the expansion coefficients are determined by H † BA in Eq. 6). Thus the number of linearly independent solutions is given by the number of ways we can choose linearly independent solutions in the subspace (φ 1 , φ 2 , ..., φ NB −NA ), which is clearly N B − N A . This proves the conjecture and the theorem.
It is easy to see why the theorem is not valid if there is on-site disorder on the majority sublattice or interactions beyond the NN, which introduces off diagonal terms in H BB . So, Eq. 4 is not true anymore. This means that Eq. 5 puts constraints on the components of Ψ B in order to satisfy the eigenvalue problem and as a result Eqs. 5 and 6 can not both be satisfied simultaneously. For example, if we use the form Eq. 7 which satisfies Eq. 6, we are only left with the freedom to choose φ 1 , φ 2 , ..., φ NB −NA and this is not enough to satisfy the eigenvalue problem of Eq. 5. There is no such problem The theorem has an important bearing on the results of the supercell calculations, both tight-binding and density functional. In these calculations, the vacancies are repeated in each supercell, connected by the superlattice translational vectors, and hence are all located on the same sublattice, which forms the minority sublattice. If n is the number of supercells in the crystal, then this is also the imbalance in the number of atoms in the two sublattices n = N B − N A . According to the theorem, there should be n zero-modes in the Brillouin zone, which is also precisely the number of Bloch momentum points in the Brillouin zone. These states thus show up in the form of a dispersionless band in the tight-binding supercell calculations at E = 0.
If hopping beyond the NN is present or if the on-site energies of the different atoms are different, then the theorem does not hold. However, the hopping beyond the NN in graphene is small [29] and the on-site energies are negligibly different on sites close to the vacancy as obtained from the DFT calculations. Because these effects are small, a clearly identifiable, nearly-dispersionless zero-mode band is found in the DFT calculations as seen from Fig. 4 as well as in the higher-neighbour tight-binding results [ Fig. 7], where the zero-mode band is marked by the red dots.

Tight-binding results: Direct diagonalization of the Hamiltonian
In order to further understand the formation of the zero-mode states, we have studied the vacancy π bands with the standard tight-binding model of the p z orbitals containing both the nearest neighbour (NN) and the second-neighbour (2NN) interactions. In particular, we look for the vacancy-induced zero-mode states discussed in the previous Subsection.

The tight-binding Hamiltonian is
where −t and t ′ are the NN and the 2NN interactions with the signs chosen such that t, t ′ > 0 (t ≈ 2.91 eV and t ′ ≈ 0.16 eV for graphene [29]). The band structures and the densities-of-states are shown in Figs. 7 and 8. The electron counting in the band structure Fig. 7 is as follows. Both the lower and the upper bands contain in total (integrated over the Brillouin zone) N A states each, while the zero-mode band contains N B − N A states, making a total of N A + N B states, as it must be the case. We have one π electron per site present in the system, so that taking spin into account, the entire lower subband is full, while the zero-mode states are half full. For the single vacancy (N B − N A = 1), this leads to a single occupied electron in the zero-mode states, resulting in S = 1/2 in agreement with the Lieb's Theorem [34].
As discussed in the previous Section, if only the NN interactions are present, we should have N B − N A number of zero-mode states at E = 0 exactly. This is why the zero-mode band in the middle panel of Fig. 7 is completely flat. However, if the 2NN interactions are also present, then the energies of the zero-mode states are not guaranteed to be the same and we see a spread in the energy of these states, which  shows up as a dispersion in the zero-mode band, as seen in the right panel of Fig. 7.
Here, the vacancy site was modeled by simply removing a lattice site, corresponding to the vacancy potential U 0 = ∞. In a real material, however, U 0 is large but finite. The effect of a finite U 0 is that (a) It causes the zero-mode state to occur slightly below the mid-gap (E = 0) and (b) The sharp zero-mode state turns into a resonance state due to interaction with the continuum π bands. This is best described with the Green's function approach discussed in the next Subsection.

Impurity Green's Function and the zero-mode state in the π bands
In this Section, we study the effect of a single impurity on the π electron states by studying the Dyson's equation and show the emergence of the zero-mode state as the strength of the impurity potential U 0 is gradually increased. For the vacancy, this potential is large but finite, so that the results obtained in this Section are helpful in understanding the nature of the zero-mode state in the actual structure with a finite vacancy potential.
The wave function of the zero-mode state is obtained from the Lippmann-Schwinger equation. Consistent with the previous results, [25,36] we find that (a) The zero-mode state consists of wave functions from the majority sublattice only and (b) It is quasi-localized decaying as 1/r as a function of distance from the vacancy in the limit of the linear-band approximation. In addition to these known results, our analysis allows us to (a) obtain the oscillatory phase factors in the zero-mode wave function due to the interference of the two Dirac points and (b) compare the lineardispersion results with the full tight-binding band result by computing the GFs for large distances in both cases.
The vacancy is modeled by adding an on-site perturbation V to the unperturbed NN (NN) tight-binding (TB) Hamiltonian, so that where H 0 = −t c † iα c jβ + H.c., iα being the site-sublattice index, and Electronic structure of the substitutional vacancy in graphene: Density-functional and Green's function studies14  where U 0 is the strength of the potential due to the impurity on the A-sublattice in the central cell. The vacancy corresponds to the value of U 0 → ∞.
The key quantity of interest here is the unperturbed GF, G 0 (E) = (E+iη−H 0 ) −1 , the calculation of which we have discussed in our earlier paper where we studied the RKKY interaction in graphene. [37] As usual, the imaginary part of the GF contains the information about the density-of-states: ρ 0 (E) = −π −1 Im G 0 (E). The GF G(E) in the presence of the perturbation will be obtained from the Dyson's equation.
Since we will be interested in the local density-of-states (LDOS) on the various carbon sites and how they are modified by the presence of the impurity, we will need to calculate the real-space matrix elements G 0 iα,jβ (E) ≡ iα|G 0 (E)|jβ . This may be done by going to the momentum space and defining the Bloch functions for the electrons | kα = N −1/2 i e i k. riα |iα with r iα = R i + τ α being the position vector of the α-th atom in the i-th unit cell. The unperturbed Hamiltonian H 0 in this and d 1 , d 2 , and d 3 are the positions of the three nearest neighbors. Diagonalization of the Hamiltonian yields the eigenenergies E( k) = ±|f ( k)|, which when expanded around the Dirac points lead to the usual linear band structure E( q) = ±v F | q|, where q = k − K D is the deviation from the Dirac point in the Brillouin zone. Here the Fermi velocity v F = 3ta/2, where 'a' is the carbon-carbon bond length. Note that unlike our previous work, [37] v F here is defined to be a positive quantity, since 't' is positive.
The real-space GFs are conveniently obtained by first calculating the momentumspace GF, which can be easily shown to be A Fourier transform then yields the real-space unperturbed GF, viz., which can be calculated by simply a brute-force summation over the Brillouin zone. It can also be computed by a second method using the Horiguchi recursive technique. [38,37] However, the latter, although computationally fast, has convergence problems [39] for distances | R i − R j | ≥ 7a or so, so that this is a better method to use only for smaller distances. The perturbed GF is related to the unperturbed GF through the Dyson's equation: G = G 0 + G 0 V G. Using the localized form of the impurity potential, Eq. 10, and taking the matrix elements, we immediately get for the real-space GF, the result We are specifically interested in the on-site GFs with α = β and R i = R j , which give the LDOS on the A and B sites at distance r iα = R i + τ α from the impurity site. Eq. (12) is easily inverted to yield the perturbed G(E) in terms of the unperturbed G 0 (E), viz., The LDOS at different sites may be obtained by taking the imaginary parts of the diagonal elements of the GF: ρ iα (E) = −π −1 Im G iα,iα (E). It immediately follows from Eq. 13 that the LDOS at the impurity site has the much simpler form where ρ 0 (E) = −π −1 Im G 0 0A,0A (E) is the unperturbed LDOS at the central site, which is of course the same for every site in unperturbed graphene, and F 0 (E) = Re G 0 0A,0A (E). Note that we have defined here ρ 0 (E) to be the unperturbed densityof-state per sublattice per spin (which is independent of the sublattice or the cell index) and ρ iα (E) is the corresponding perturbed quantity for the iα site.
The resonance condition follows from Eq. 14, viz., where E 0 is the resonance energy. The graphical solution of this equation is shown in Fig. 9. There are four solutions for E 0 : The two solutions at E 0 = ±t do not produce much change in the DOS, as may be inferred from Eq. 14, due to the diverging density-of-states ρ 0 (E) there, and the bound state with E 0 → U 0 for large U 0 is inconsequential because it is removed to ∞, which then leaves the sole resonance state indicated by the black dot in Fig. 9. Its energy goes to zero in the limit U 0 → ∞ and the oscillator strength to one, producing the zero-mode state for the vacancy.
The total DOS in the presence of the perturbation may be computed by taking the trace of Eq. (13) for the entire lattice. Using the identity iα G 0 iα,0A (E)G 0 0A,iα (E) = −dG 0 0A,0A (E)/dE and some tedious algebra, the result is Similarly, by taking the trace of Eq. (13) over the cell index only, the individual sublattice DOS may be found, which for the A sublattice reads as where ξ(E) = (1/N ) k [G 0 AA (kE)] 2 and the densities of states are, again, per sublattice and per spin. A similar expression for ρ B (E) reads . (18) It can be verified from Eqs. (16) -(18) that ρ tot (E) = ρ A (E) + ρ B (E), so that these equations are consistent. The numerical results are summarized in Figs. 10, 11, and 12. The factors multiplying the 1/N in Eqs. 17 and 18 are the changes in the DOS ∆ρ A (E) and ∆ρ B (E) introduced by the impurity potential, which are shown in Fig. 10. Fig. 11 shows the emergence of the zero-mode in the density-of-states with E = 0 and that this state resides entirely on the B sublattice in the limit U 0 → ∞. Fig. 12 shows the LDOS on the impurity site (ρ 0A ) and on the nearest (ρ 0B ) and the next nearest sites (ρ 1A ).
The width of the resonance peak increases with the resonance energy E 0 . Keeping the linear term in the expansion of F 0 (E), viz., where the resonance width is given by The width is zero if E 0 = 0 and increases with energy as shown in Fig. (13).

Impurity state wave function
In this section, we study the impurity state Vπ wave function from the Lippmann-Schwinger equation. The analysis allows us to obtain the well known 1/r decay of the vacancy state; however, in addition we also obtain the oscillatory behavior of the wave function due to the interference between the two Dirac cones.
Our starting point is the Lippmann-Schwinger equation |Ψ = |Ψ 0 + G 0 V |Ψ , where |Ψ 0 is the unperturbed state. For the localized impurity potential on the central A site, V = U 0 |0A 0A|, the Lippmann-Schwinger equation leads to the wave function We are interested in the low-energy behavior, since that's the region where the resonance state gets introduced by the impurity as seen from Fig. 9. The GFs for the full tight-binding band structure as well as for the linear bands were computed in our previous work. [37] For the linear bands and in the low energy limit, the results are Electronic structure of the substitutional vacancy in graphene: Density-functional and Green's function studies18 (E) E/t Figure 12. Local density-of-states at the impurity site ρ 1A (top), the NN site ρ 0B (middle), and the next NN site ρ 1A (bottom) obtained from Eqs. 13 and 14 for different strengths of the impurity potential U 0 /t = 0, 2, and 5, denoted by black dashed, black solid, and red dashed lines respectively. As U 0 → ∞, the top LDOS goes to zero (except for the bound state beyond the top of the band whose energy goes to ∞), and the zero-mode state lives only on the B sublattice, as indicated from the middle and the bottom panels. The prominent zero-mode peak in the middle panel for U 0 /t = 5 will develop into a δ-function peak at E = 0 as the impurity potential U 0 → ∞.
where A c is the unit cell area in graphene, K 0 and K 1 are the modified Bessel functions of the second kind and r is the distance vector between the two atoms: r = r iA − r 0A for the first GF and r = r iB − r 0A for the second. The multiplicative factors are β = e i K· r + e i K ′ · r , which is a real number for the graphene lattice and which is purely imaginary and the polar angle θ r = tan −1 (y/x) is defined with the x  [29] if only the NN hopping is kept. direction taken to be along the vector K ′ − K, which are two adjacent Dirac points in the Brillouin zone.
Using the small z expansion for the Bessel functions: [40] K 0 (z) = − ln(z/2)−γ and K 1 (z) = z −1 + 2 −1 z ln(z/2) + (γ − 1/2)z/2, where γ ≈ 0.577 is the Euler-Mascheroni constant, we find the Bessel functions entering the expressions for the GFs (Eq. 22) to become, in the low energy limit: Plugging these into Eq. 22, taking the E → 0 limit , and retaining the lowest-order terms in energy, we find the results There is no guarantee that these results, calculated for the linear bands extending to infinite energy, should agree even for low energies with the GFs for the actual bands, e.g., as obtained with the tight-binding band structure. Certain elements must exactly agree at low energies, for example, the imaginary part of G 0 0A,0A , which yields the density-of-states, since at low energies, it is controlled by the Fermi velocity v F alone. We nevertheless find that the expressions Eq. 25 do agree quite well with the tight-binding GFs, the agreement becoming better with increasing distances. A comparison between the low-energy GFs Eq. 25 and the full GFs for a typical case is shown in Fig. 14, which also illustrates the symmetry of the GFs. A notable exception is the real part of the on-site GF G 0 0A,0A (E), where the substitution of r = 0 in Eq. 25 yields a divergent result. However, we find that the tight-binding GF in this case can be fitted to the expression Eq. 25 for G 0 0A,0A (E), provided we use the value r ≈ 0.6 a instead of r = 0.
The nature of the impurity state immediately follows from the Lippmann-Schwinger expression Eq. 21. First of all, notice an important point from the expression for the GF Eq. 25, viz., that all GFs vanish at E = 0 except for the real part of G 0 iB,0A , which is finite and decays as 1/r. This is precisely what leads to the property that the zero-mode state resides on the majority sublattice B only and its wave function decays inversely with distance. These features are true if only the NN interactions are present on the graphene lattice. If second NN interactions are present, then there is no electron-hole symmetry and the behavior of the GFs near the resonance energy differs from Eq. 25. The form of the GFs for the latter case is such that both sublattices contribute to the resonance state near E = 0, an issue that is discussed in detail elsewhere. [41] Returning to the Lippmann-Schwinger equation Eq. 21 and inserting into it the low-energy expansion for the GFs (Eq. 25) and then taking the limit of the resonance energy E 0 = 0, it can be easily seen that as E 0 → 0 in the limit U 0 → ∞, the impurity wave function follows the behavior This is an important result, which states that in the NN approximation, only the B sublattice component survives for the zero-mode state, it being the stronger infinity. The surviving component is found to be simply proportional to the real part of the inter-sublattice GF, since its imaginary part vanishes. Using Eq. 25 and evaluating Im α from Eq. 23, we finally get the desired result Ψ B (r) = N r sin[( K − K ′ ) · r/2 − θ r ] cos[( K + K ′ ) · r/2 − π/3], (28) where we have suppressed the cell index i, N is a constant, r is again the actual distance vector of the B site with respect to the impurity position, and the two Dirac points in the Brillouin zone may be taken as: K = 2πa −1 3 −3/2 (−1, √ 3) and K ′ = 2πa −1 3 −3/2 (1, Eq. 28 is the central result of this Subsection that describes the 1/r decay of the vacancy-induced Vπ state along with the phase factors. The long-range nature 1/r of the wave function (28) is well-known, [36] but the oscillatory factor due to the interference effect of the two Dirac points is new. The same kind of interference is also present in the oscillations of the RKKY interactions. [42,37] The wave function is not square integrable because we used the linear band structure, but it will be if we take the full band structure into account. Eq. 28 nevertheless describes the gross features of the zero-mode state. The wave function changes sign along different directions, e. g., it changes sign along the zigzag direction but not along the armchair direction. The kinetic energy gained by the delocalization of the wave function is exactly cancelled by anti-bonding components present in the wave function, so that its energy still equals the on-site energy in spite of the delocalization. The calculated wave function for the zero-mode state is shown in Fig. 15. We note that a recent study has shown that the 1/r decay of the vacancy state remains unchanged even when a repulsive Coulomb interaction is included in the tight-binding Hamiltonian. [43]

Summary
In summary, we have studied the electronic structure of graphene with a single substitutional vacancy from density-functional calculations using the all-electron LAPW method and interpreted the results with the help of the tight-binding model and the impurity Green's Function approach. We find that the vacancy induces four localized states, viz., three Vσ dangling bond states on the carbon triangle and one Vπ resonance state. The dangling bond states cause a Jahn-Teller distortion, which we found to be a planar distortion of the carbon triangle. Hund's coupling between these electrons would then produce the S = 1 state at the vacancy center as indicated in the summary figure Fig. 1. The magnetic moment has two components: (i) The component 2µ B coming from the localized vacancy states Vσ and Vπ and (ii) An opposite component of several tenths of µ B coming from the spin-polarization of the continuum π band states in the vicinity of the vacancy. The second part is not well described in the supercell band calculations due to the slow 1/r decay of the "quasi-localized" Vπ wave function. This long-range decay also means that in an experimental sample it is only for the extremely low vacancy concentration that the truly isolated vacancy limit is reached and as a result the magnetic moment is expected to be dependent on the vacancy concentration.
In addition to the density-functional calculations, we also studied the formation of the Vπ state in detail from the impurity Green's function approach for the isolated vacancy, which provided important insight in the interpretation of the results of the band calculations and the formation of the zero-mode states in the π bands. This zero mode state is a slowly-decaying localized state that lives mostly on the majority sublattice. It spreads into the minority sublattice (the one containing the vacancy) and becomes a resonance state due to the second and the higher-neighbor interactions as well as the finite strength of the vacancy potential. The Green's function approach provided a sinusoidal phase factor associated with the Vπ wave function described by Eq. 28. In addition to the understanding of the vacancy electronic structure, our work provides important insight necessary for the understanding of impurities in general such as iron and cobalt dopants and other complex defects.