Coulomb charging energy of vacancy-induced states in graphene
Vladimir G. Miranda, Luis G. G. V. Dias da Silva, Caio H. Lewenkopf
CCoulomb charging energy of vacancy-induced states in graphene
V. G. Miranda, Luis G. G. V. Dias da Silva, and C. H. Lewenkopf Instituto de F´ısica, Universidade Federal Fluminense, 24210-346 Niter´oi, RJ, Brazil Instituto de F´ısica, Universidade de S˜ao Paulo, C.P. 66318, 05315–970 S˜ao Paulo, SP, Brazil (Dated: August 27, 2018)Vacancies in graphene have been proposed to give rise to π -like magnetism in carbon materials,a conjecture which has been supported by recent experimental evidence. A key element in this“vacancy magnetism” is the formation of magnetic moments in vacancy-induced electronic states.In this work we compute the charging energy U of a single-vacancy generated localized state for bulkgraphene and graphene ribbons. We use a tight-binding model to calculate the dependency of thecharging energy U on the amplitudes of the localized wave function on the graphene lattice sites. Weshow that for bulk graphene U scales with the system size L as (ln L ) − , confirming the predictionsin the literature, based on heuristic arguments. In contrast, we find that for realistic system sizes U is of the order of eV, a value that is orders of magnitude higher than the previously reportedestimates. Finally, when edges are considered, we show that U is very sensitive to the vacancyposition with respect to the graphene flake boundaries. In the case of armchair nanoribbons, wefind a strong enhancement of U in certain vacancy positions as compared to the value for vacanciesin bulk graphene. PACS numbers: 73.22.Pr, 75.75.-c,73.20.Hb
I. INTRODUCTION
Lattice defects have long been considered as an unde-sired presence in micro- and nanostructured devices. Inmany cases, they were deliberately avoided in the manu-facturing processes since they modify the electronic andstructural properties and are detrimental to electronictransport. More recently, this scenario has been changingas it has been shown that defects themselves can give riseto interesting physical phenomena in many condensedmatter systems.Graphene is one of the main platforms that is con-tributing to the growing interest in defective nanostruc-tures. In addition to the remarkable and already thor-oughly explored properties of clean graphene flakes, experiments in graphene with vacancies reveal a newroute to extend the plethora of fascinating aspects of thismaterial. .Some of the important discoveries pointed out by ex-periments in graphene with vacancies include the onset ofmagnetic behavior in a p -block system, signaturesof the Kondo effect in the transport properties andthe recently reported atomic charge collapse. Theory predicts that a vacancy originates a midgapstate pinned at E = 0 and localized around the “vacancysite”. This result is not exclusive to graphene andgeneralizes to all class of systems that can be representedby a bipartite-like hamiltonian with nearest neighbor in-teractions only . For such systems, a counting rulestates that whenever an imbalance N I = | N A − N B | (cid:54) = 0between the number of constituents of the sublattices A and B exists, at least N I states pinned at zero energyrise. The picture described above relies on a simple modelfor the vacancy considering graphene π -band electronsonly. This “ π -like” magnetism in graphene with va- cancies has been reported experimentally recently, though its relevance for the vacancy-induced magnetismis a matter of a long debate in the DFT community. As we show next, our results reinforce the arguments infavor of the relevance of π -like magnetism in graphenewith vacancies. Moreover, our arguments should be alsovalid for H adatoms, which can be described by an onsitescalar potential model. A crucial aspect of the vacancy state and that is keyfor the understanding of the onset of magnetism, Kondoeffect, and atomic collapse within this model is the charg-ing energy of the vacancy state U . This parameter,also referred to as the “Hubbard-U”, encodes the lo-cal electron-electron Coulomb interaction of the vacancystate. As such, it leads to the spin splitting of thevacancy midgap state which is essential for the onsetof magnetic behavior. In addition, the interplay be-tween U and the hopping between the midgap state intothe band states sets the condition for the onset of theKondo effect. Moreover, it has been argued thatthe strength of U controls the vacancy charge and thusthe critical coupling that determines the appearance ofthe atomic collapse .The central focus of this work is the determination ofthe charging energy U of vacancies in graphene. We ad-dress this issue for bulk graphene as well as for graphenewith edges. Previous analytical results on the wave func-tion of a vacancy in bulk graphene and graphenewith edges are extremely relevant for the results wepresent.We recall that, for a single vacancy in bulk graphene,the vacancy wave function decays as ∼ /R with R beingthe distance from the vacancy site. This unusualbehavior produces a non-normalizable wave function andhence the state is said to be quasi-localized with the de-gree of localization of the system decaying as the system a r X i v : . [ c ond - m a t . m e s - h a ll ] M a y enlarges. This feature has led some authors to claim thiswould make U negligible for real samples. It is well known that the presence of edges stronglymodifies the electronic properties of graphene.
Whenedges are introduced, important changes occur also onthe vacancy state.
Analytical accounts on theinfluence of edges in the vacancy-generated state ingraphene have been recently put forward in Refs. 35and 36. In these works, Deng and Wakabayashi stud-ied the influence of edges on the vacancy localized stateof a semi-infinite sheet and graphene nanoribbons. Interestingly, for the semi-infinite systems in the pres-ence of an armchair edge, the vacancy state decays as1 /R and is thus normalizable. When a zigzag edge isconsidered, the influence on the vacancy state dependson which sublattice the vacancy is created. Deng andWakabayashi find that if the vacancy and the edge sitesbelong to different sublattices, the vacancy causes no ef-fect on the system zero energy states. In contrast, whenboth, edge and vacancy belong to the same sublattice,the vacancy yields a strongly distorted, non-normalizablezero-mode wave function around the vacancy site. Thecase of vacancies in graphene nanoribbons was stud-ied by the same authors in Ref. 36. They claimthat a single-vacancy has no effect on the zero-energystates of zigzag-terminated edge ribbons, quantum dots,and armchair-terminated metallic ribbons. However, forarmchair-terminated semiconductor ribbons, the vacancygenerates a square-normalizable wave function pinned at E = 0.Although the main focus of the present paper is theanalysis of the charging energy of the vacancy-generatedstate in bulk graphene and graphene with edges, wealso perform a systematic study of the vacancy-generatedwave functions for graphene with edges. We combine an-alytical and numerical techniques to study how the prop-erties of the vacancy state and its respective U is alteredfor different system sizes, edges and vacancy-edge-sitesdistance.One of the most important results of this work is thederivation of an analytical expression for the computa-tion of U . This result reveals that in addition to the inter-site contribution to the Coulomb charging at thevacancy, there is also a sizable intra-site contributionwhich has been so far overlooked. We note that, forbulk graphene, Ref. 33 finds that the 1 /R decay of themidgap state leads to U ∼ (2 π ln L ) − , with L being thelinear system size. Based on this scaling behavior, the au-thors estimate U ∼ U scaling predicted in Ref. 33, how-ever our U estimates are orders of magnitude larger thantheir predictions and in line with RPA/Hubbard modelcalculations. Very recently, scanning tunneling spectroscopy (STS)experiments in graphene deposited on a Rh foil reportthe appearance of spin-split states near vacancy sites,which is consistent with the presence of strong onsiteinteractions. The STS data show splittings of about 20 - 60 meV, which can, in principle, be directly comparedto the effective U in their system (which depends on theeffective dielectric constant at the vacancy site).For graphene ribbons, U calculations have been per-formed for the case of armchair ribbons in Ref. 37. Theauthors find out that U is related to the inverse par-ticipation ratio (IPR) of the vacancy state, in agreementwith the results we present below. However, their predic-tion that U vanishes for increasing ribbon widths is notsupported by our study. In this paper, we also perform asystematic study of the of the IPR and U of the vacancystates for different ribbon widths and edge-vacancy-sitesdistances, a study lacking in the literature to the best ofour knowledge. Our results point that for some vacancy-edge-sites distances, U decreases with increasing ribbonwidths and approaches the bulk estimates, however re-maining close to ∼ U predicted elsewhere. We alsofind that there are some vacancy-edge-sites distance con-figurations for which a directionality effect of the vacancywave function makes the IPR and U system-size inde-pendent, being a truly localized state. This is the secondmain contribution of our work.The paper is structured as follows. In Sec. II, wepresent the tight-binding formalism used to obtain thevacancy wave functions and to derive the analytical ex-pression of U . In Sec. III, we use the results of the pre-vious section to obtain U estimates in bulk graphene forvarying system sizes. In Sec. IV, we make a thoroughstudy of the vacancy state in graphene armchair ribbonsfor different configurations. We calculate the IPR of themidgap states to quantitatively evaluate the degree of lo-calization of such states. We finish the section addressingthe issue of vacancies in the presence of zigzag and quan-tum dots with both zigzag and armchair edges and showthat our main findings for the armchair edges seems toremain robust to the additional presence of zigzag edges.In Sec. V we evaluate U for armchair ribbons and showthat U mimics the IPR behavior. Also, we show thatour results hold irrespective if the ribbons are semicon-ducting or metallic. Finally, we present our concludingremarks in Sec. VI. II. VACANCY-INDUCED MIDGAP STATES:WAVE FUNCTION AND CHARGING ENERGY
In this section, we present a single-orbital tight-bindingmodel description of the midgap states due to a singlevacancy in graphene monolayer systems. We then de-rive analytical expressions for the Coulomb charging en-ergy for these localized states. For notation compact-ness, we consider the charging energy U in vacuum. Ifthe graphene layer is in contact with another medium,the calculated U should be divided by the correspondingdielectric constant (cid:15) .We describe single-particle spectrum and wave func-tion using the model Hamiltonian given by H = H + V, (1)where V accounts for a single monovacancy and H is thepristine graphene tight-binding Hamiltonian, namely H = − t (cid:88) (cid:104) i,j (cid:105) ( | i (cid:105)(cid:104) j | + H . c . ) (2)where t ≈ . | i (cid:105) corresponds to a p z orbital placed at the i th site of thegraphene honeycomb lattice with interatomic separation a = 1 .
41 ˚A, (cid:104)· · · (cid:105) restricts the sums to nearest neighborssites.There are several equivalent ways to account for thevacancy.
For analytical calculations it is convenientto model a monovacancy placed at the site v by an on-sitepotential term, namely, V = V | v (cid:105)(cid:104) v | (3)and take the limit | V /t | (cid:29)
1. Alternatively, one can alsouse V (cid:48) = t (cid:88) (cid:104) v,j (cid:105) ( | v (cid:105)(cid:104) j | + H . c . ) , (4)which corresponds to turning off the hopping terms thatconnect the vacancy site v to its nearest neighbors. Fromthe numerical point of view both models, V and V (cid:48) , givethe same results for the low-energy single-particle proper-ties of the system. (This statement is also corroboratedby the good agreement between our results and those ofRef. 35.)The single-particle electronic properties of a systemwith an impurity can be analytically obtained from a T -matrix analysis. The T -matrix for a Hamiltonian ofthe form H = H + V reads T = V (1 − V G ) − , where G = ( E + iη − H ) − is the Green’s function of thepristine graphene system. For | V /t | (cid:29)
1, the T matrixreduces to T ( E ) = − | v (cid:105)(cid:104) v |(cid:104) v | G ( E ) | v (cid:105) . (5)In general, the system wave functions are obtainedfrom the Lippmann-Schwinger equation, | ψ (cid:105) = (1 + G T ) | ψ (0) (cid:105) , namely ψ E ( i ) = ψ (0) E ( i ) − (cid:104) i | G ( E ) | v (cid:105)(cid:104) v | G ( E ) | v (cid:105) ψ (0) E ( v ) (6)where ψ (0) E ( i ) = (cid:104) i | ψ (0) E (cid:105) is the unperturbed wave functionamplitude at the site i , solution of H | ψ (0) E (cid:105) = E | ψ (0) E (cid:105) .The knowledge of (cid:104) i | G ( E ) | j (cid:105) allows one to analyticallycalculate ψ E ( i ).Of particular interest is the vacancy-induced midgapstate | ψ (cid:105) at E = 0. However, for these states, a care-ful study of the T ( E →
0) behaviour should be taken into account since, depending on the particular geome-try under interest, T ( E →
0) diverges and the strategyoutlined above to obtain the wave functions is not alwaysstraightforward.
In addition, it has been found thatif the system has a finite gap, as in a armchair semicon-ducting ribbon, the vacancy-induced midgap state be-comes a truly bound state and the wave function is nolonger given by the expression above (see, for instance,Refs. 35 and 36 for a thorough discussion of this case).We also analyze the midgap state problem numerically.The vacancy breaks down the lattice translational invari-ance. For bulk systems, we find it convenient to intro-duce a supercell of dimension N tot = N × M sites withperiodic boundary conditions at its edges (the descrip-tion of the unit cell for graphene nanoribbons is pre-sented in Sec. IV). Figure 1 shows the superlattice geom-etry: the atomic lattice sites i are denoted by the labels( m, n, S ), where n identifies the zigzag chain (withrespect to the vacancy site) to which the site i belongsand m is the position of the i − th site within the chain,while S = A or B is the sublattice index. m .... N z i g z ag c ha i n s x(m,n) .... − − − .... n M · · · m M · · · x ( m, n ) ······ n N ··· z i g za g c h a i n s N /2 -N /2+1 FIG. 1. Schematics of the supercell used in this work. Thesupercell is defined by M/ N “horizontal” zigzag ones. m is the discrete horizontal positionwith respect to the edge and n labels the horizontal zigzagchains starting from the one containing the vacancy, indicatedby the cross. The single-orbital graphene tight-binding Bloch basisis given by χ k j ( r ) = 1 √ N sc (cid:88) R (cid:48) e i k · R (cid:48) φ ( r − t j − R (cid:48) ) , (7)where the sum runs over all N sc supercells centered at R (cid:48) , t j gives the position of the j th carbon atom in thesupercell, and φ ( r ) is the p z orbital atomic wave function.We take the limit of N tot = N × M (cid:29)
1. Due toband folding, the Γ point gives a good representationof the first Brillouin zone of the supercell.
For thisreason, we restrict our calculation to k = 0 and drop the k label from now on. The crystal electronic single-particleeigenstates are the solutions of H Ψ ν ( r ) = E ν Ψ ν ( r ) andread Ψ ν ( r ) = (cid:88) i ψ ν ( i ) χ i ( r ) , (8)where ψ ν ( i ) is the ν th tight-binding wave function am-plitude at the i th site and χ i ≡ χ k =0 ,i . The midgapvacancy-induced state corresponds to ν = 0. For laterconvenience, we introduce the envelope wave functions ψ ν ( r i ) ≡ √A ψ ν ( i ) , (9)where A is the supercell area. This approach allows usto calculate ψ ( r i ) by direct diagonalization.Let us now analyze the charging energy U correspond-ing to the Coulomb energy associated with a double oc-cupation of the midgap state | ψ (cid:105) . It is rather temptingto use the envelope wave function ψ ( r ) to evaluate U , namely U = e (cid:90) d r (cid:90) d r (cid:48) | ψ ( r ) | | ψ ( r (cid:48) ) | | r − r (cid:48) | . (10)As shown below, this is the main contribution to thecharging energy U due to the π orbitals. However wefind another important contribution, which has been ne-glected so far.The Coulomb energy associated to the double occupa-tion of the vacancy-induced midgap state reads U = e (cid:90) d r (cid:90) d r (cid:48) | Ψ ( r ) | | Ψ ( r (cid:48) ) | | r − r (cid:48) | , (11)where Ψ ( r ) is the three dimensional wave function of themidgap state given by Eq. (8). U can then be expressedin terms of the tight-binding amplitudes as U = (cid:88) i,j,i (cid:48) ,j (cid:48) ψ ∗ ( i ) ψ ( j ) ψ ∗ ( i (cid:48) ) ψ ( j (cid:48) ) W ij,i (cid:48) j (cid:48) , (12)where the sums run over all N tot atomic sites of the su-percell and W ij,i (cid:48) j (cid:48) = e (cid:90) d r (cid:90) d r (cid:48) φ ∗ ( r − t i ) φ ( r − t j ) 1 | r − r (cid:48) | φ ∗ ( r (cid:48) − t i (cid:48) ) φ ( r (cid:48) − t j (cid:48) ) . (13)Equation (12) does not include contributions from atomic orbitals located at different supercells, namely, R (cid:54) = R (cid:48) .Since the Coulomb integral decreases rapidly as the orbital centers are separated, one only expects significant inter-supercell contributions from orbitals located at the edges of neighboring supercells. Those correspond roughly to afraction 1 / √ N tot of the supercell sites and can be safely ignored in the limit of N tot (cid:29) W ij,i (cid:48) j (cid:48) in the context of generalizing the tight-binding ideas to obtain an atomistic total-energy method(see, for instance, Refs. 54–57 and references therein).The leading matrix elements W ij,i (cid:48) j (cid:48) correspond toan intra-atomic (on-site) Coulomb repulsion matrix el-ements, where all orbitals belong to the same atom, andto inter-atomic (non-local) terms, where i = j is in oneatom and i (cid:48) = j (cid:48) on another.Accordingly, we decompose U as U = U + U , (14)where U consists of intra-atomic Coulomb repulsionterms, while U contains the inter-atomic ones.Let us first consider U , namely U = e (cid:90) d r (cid:90) d r (cid:48) | φ ( r ) | | φ ( r (cid:48) ) | | r − r (cid:48) | (cid:88) i | ψ ( i ) | , (15)where the Coulomb integral U orbital = e (cid:90) d r (cid:90) d r (cid:48) | φ ( r ) | | φ ( r (cid:48) ) | | r − r (cid:48) | , (16) can be evaluated, for instance, by an expansion of | r − r (cid:48) | − in spherical harmonics. Equation (16) rep-resents the Hartree contribution to U . In the literature U orbital was estimated to be ∼
17 eV for free standinggraphene and, if screening effects from electrons ofbands other than the π are taken into account, U orbital reduces to ∼ . Hence U = U orbital (cid:88) i | ψ ( i ) | . (17)Let us now address U , which is more convenientlyexpressed by changing the integration variables as r → r − t i and r (cid:48) → r (cid:48) − t j , namely U = e (cid:88) i (cid:54) = j | ψ ( i ) ψ ( j ) | (cid:90) d r (cid:90) d r (cid:48) | φ ( r ) | | φ ( r (cid:48) ) | | r − r (cid:48) + t i − t j | . (18)We write U = e (cid:88) i (cid:54) = j | ψ ( i ) ψ ( j ) | | t i − t j |× (cid:90) d r (cid:90) d r (cid:48) | φ ( r ) | | φ ( r (cid:48) ) | (cid:115) | δ r | | t i − t j | + 2 δ r · ( t i − t j ) | t i − t j | , (19)where δ r = r − r (cid:48) . The orbital wave functions am-plitudes φ ( r ) decay quickly for r/a (cid:38) | r − r (cid:48) | /a (cid:38)
1. These observations suggest that U canbe approximated by the lowest order Taylor expansionin powers of | δ r | / | δ t | of the square root at the r.h.s. ofEq. (19). This can be checked quantitatively by compar-ing with the exact values of the integral in Eq. (18). Wefind that our approximation overestimates U by about ≈ Hence, U ≈ e (cid:88) i (cid:54) = j | ψ ( i ) | | ψ ( j ) | | t i − t j | (cid:90) d r (cid:90) d r (cid:48) | φ ( r ) | | φ ( r (cid:48) ) | ≈ e (cid:88) i (cid:54) = j | ψ ( i ) | | ψ ( j ) | | t i − t j | , (20)since the orbitals are normalized.In the continuum limit, the sums in Eq. (20) can bechanged to integrals and the site amplitudes can be re-placed by an envelope wave function ψ ( r ), leading to U ≈ e (cid:90) d r (cid:90) d r (cid:48) | ψ ( r ) | | ψ ( r (cid:48) ) | | r − r (cid:48) | . (21)Many-body corrections to intra and inter atomicCoulomb repulsion terms of Eq. (13) have been cal-culated in Ref. 59. The latter uses Wannier-like or-bitals projected in the p z bands. This offers the ad-vantage of separating the p z contribution to the chargescreening thereby accounting for the effective partialtwo-dimensional screening expected for the electrons ingraphene.In summary, the general expression for the Coulombrepulsion term U is given in terms of a three dimen-sional integral given by Eq. (11). We show how the lattercan be cast into a two dimensional form proposed in theliterature. We also find an additional significant contri-bution to U that is proportional to U orbital , given by Eq.(17). III. U OF VACANCY-INDUCED LOCALIZEDSTATES IN BULK GRAPHENE
In this Section, we study the dependence of U with thesystem size. In order to reach our goal, we recall that the wave function ψ ( r ) of a localized state due to a single-vacancy in a clean bulk graphene monolayer reads ψ ( r ) = N r sin (cid:104) ( K − K (cid:48) ) · r − θ r (cid:105) × cos (cid:104) ( K + K (cid:48) ) · r − π (cid:105) , (22)where r = ( x, y ) is a coordinate vector with origin atthe vacancy site, N is the normalization constant, K =2 π/ (3 √ a )( − , √
3) and K (cid:48) = 2 π/ (3 √ a )(1 , √
3) denotethe two inequivalent Dirac points in the first Brillouinzone, and θ r = arctan( x/y ). Here, we consider finitesystems with N tot = L × L sites. Having an analyticalexpression for ψ ( r ) is key for the study of systems with L (cid:38) (or larger than 10 sites), since the computationtime for exact diagonalization scales with L .Pereira and collaborators showed that for bulkgraphene (cid:88) i | ψ ( i ) | ∝ L ) , (23)and thus U ∝ U orbital (ln L ) . (24)Similarly, a rough estimate of U can be obtained byusing the approximation ψ ( r ) ≈ N /r . By normalizingthe envelope wave function, the charging energy given byEq. (21) reads U ∝ e /(cid:15) (ln L ) − .We use ψ ( r ) given by Eq. (22) to estimate the charg-ing energy U of the vacancy-induced state as a functionof L . We insert the lattice wave function amplitudes ψ ( i ) = √A ψ ( r i ) in Eqs. (17) and (20) to numericallyobtain U and U , respectively. To account for the effectof the substrate, we use (cid:15) = 4, which is consistent withthe value measured for graphene deposited on SiO . Figure 2 gives U , U , and U approx = U + U as a func-tion of L . We find that the dependence of the chargingenergy (in eV) with system size is accurately fitted by U approx = 0 .
32 + 82(ln L ) − . Figure 2 shows that the U estimates are almost an order of magnitude larger thanthose of U .The extrapolation of our results to graphene sheetswith areas of about 1 µ m (corresponding to L ≈ or10 sites), gives U ≈ . as well as from DFT calcula-tions for small lattice sizes. Therefore, our results showthat the π − like magnetism is not negligible in model sys-tems with realistic sample sizes, in line with experimentalevidences. IV. GRAPHENE WITH EDGES
In this section we study the influence of the edges onthe vacancy-induced states. More specifically, we present (log L ) ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● . . . U ( e V ) ( ln L ) − ● U off + U diag U off U diag U ( e V ) U +U U U U ( e V ) ( ln L ) − U + U U U (log L ) ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● . . . U ( e V ) ( ln L ) − ● U off + U diag U off U diag U ( e V ) U +U U U (log L ) ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● . . . U ( e V ) ( ln L ) − ● U off + U diag U off U diag U ( e V ) U +U U U (ln L ) FIG. 2. (Color online) Scaling of the diagonal and off-diagonalcharging terms U and U as a function of system size. Valuesare in eV. a systematic study of the degree of localization and thecharging energy U of the vacancy-induced states in arm-chair nanoribbons. We also provide a comparative analy- sis with the cases of vacancies in zigzag nanoribbons andquantum dots. The main observation is that both local-ization and U are strongly dependent on the vacancy-edge distance. A. Midgap state wave function in the presence ofarmchair edges
We follow two alternative approaches to study thecharacteristics of the vacancy-induced state: one thatconsiders the analytical expression for the vacancy statederived in Ref. 36 and a numerical one which is based onthe numerical diagonalization of the tight binding Hamil-tonian described in Sec. II with armchair boundary con-ditions.Let us consider, without loss of generality, a vacancycreated at sublattice A. The wave function of | ψ (cid:105) re-sides solely on sublattice B. In the notation presented inSec. II, the expression for ψ ( m, n, S ) presented in Ref. 36is written as: ψ ( m, n, S ) = − (cid:88) r I r ( n ) (cid:20) cos (cid:18) πr ( x B ( m, n ) − x ( m , M + 1 (cid:19) − cos (cid:18) πr ( x B ( m, n ) + x ( m , M + 1 (cid:19)(cid:21) δ S,B M + 1 . (25)We recall that m and n label a given site i and thesubindex 0 stands for the vacancy site (see Fig. 1). Here,armchair edges correspond to sites with m = 1 , m = M − , M in Fig. 1. Following Ref. 36, n = 0stands for the zigzag row where the vacancy is posi-tioned and the index increases (decreases) as one movesto rows on top (bottom), see left axis of Fig. 1. Fi-nally, x B ( m, n ) = m/ B lattice site labeled by( m, n ).The function I r ( n ) is given by I r ( n ) = 2( − n (cid:26) Θ( π − k r )[2 cos( k r )] − ( n +1) , n ≥ k r − π )[2 cos( k r )] (1 − n ) , n < , (26)where the “wave number” k r obeys the quantizationrule dictated by M , the finite number of sites along thetransversal direction (see lower panel in Fig. 1): k r = 2 πM + 1 r, r = 1 , , · · · , M/ , (27)where r is the band index. For the case where mod( M +1 , (cid:54) = 0, the nanoribbon is semiconducting, oth-erwise the system is metallic. In Sec. V, we use theseanalytical expressions for the computation of the charg-ing energy.We also analyze the effects of vacancies in armchairgraphene nanoribbons by numerical diagonalization ofthe tight binding Hamiltonian H = H + V (cid:48) , see Eqs. (2) and (4). As before, due to the lack of translationalsymmetry, we consider supercell with a N tot = M × N sites. We use periodic boundary conditions at the zigzagchains characterized at n = N/ n = − N/ and is located solely in thesublattice opposite to that of the vacancy. In Fig. 3, theblack dot indicates the vacancy site. The radii of the bub-bles are proportional to the amplitude of the wave func-tions and the blue (red) color stands for positive (nega-tive) sign of the wave function. Note that the sites thatcarry the larger weights are those closer to the vacancysite.As expected, the vacancy originates a midgap boundstate pinned at zero energy. We show in the nextsection that this numerical approach provides resultswhich are in very good agreement with the analyticaltreatment discussed above.
B. Midgap state degree of localization versusvacancy-edge distance.
Although Ref. 37 presents an extensive study of thebehavior of vacancies created in semiconducting armchairribbons, one important aspect that was not properly ex-
IPR=0.121 P = 0 . (a) IPR=0.0151 P = 0 . (b) FIG. 3. (Color online) Wave function amplitudes for a va-cancy at two different sites (black dots) in a armchair ribbon:(a) vacancy two sites away from the edge and (b) vacancy onesite away from the edge. Note that the degree of localizationdrops by an order of magnitude as the vacancy is moved onesite towards the edge. plored so far is the dependency of the degree of localiza-tion of the vacancy state when the position of vacancysite is changed. Before addressing the calculation of theCoulomb charging energy for the vacancy-localized state,we discuss the nature of the localized wave functions asa function of the vacancy-edge distance D .We characterize the degree of localization of the wavefunctions by the Inverse Participation Ratio (IPR) which, for a given state of energy E ν , reads: P ν = (cid:88) i | ψ ν ( i ) | , (28)The IPR is contained in the interval (0 ,
1] and the closerto the upper (lower) limit, the higher (lower) is the degreeof localization of the wave function. Note that, for thevacancy state, the P = U /U orbital . Hence the behaviorof P mimics the one followed by U .We find that the vacancy state is extremely sensitiveto the vacancy-edge distance as shown in Fig. 3. Fig-ure 3a shows the midgap state for a vacancy placed twosites away from the edge ( D = 2, in units of a √ / D = 1, as shown in Fig. 3b), the vacancy state changesabruptly and the wave function extends much more thanin the previous configuration. This behavior is qualita-tively seen in Fig. 3 and quantified by the IPR. Noticethat the IPRs of the two configurations differ by an orderof magnitude.Our numerical calculations show a non-monotonicaldecrease of the P as a function of the vacancy-edge dis-tance. For very narrow ribbons, this effect is very subtle(see, for instance, the behavior of the IPR for the ribbonwith width M = 6 in Fig. 4). For small M , the confine-ment due to a finite width competes with the localizationdue to the vacancy. Figure 4 indicates that the IPR doesnot depend on M for some specific sites, while for othersites the IPR decreases with increasing M . The size inde-pendent IPRs correspond to truly localized states due tothe vacancy, while the others behave as quasi-localizedstates as those due to vacancies in bulk graphene .This interpretation plays a key role in our analysis and,to the best of our knowledge, has been unnoticed so far. ●●●●●●●●●●●●●●●●●●●● D I P R ● M = = = = = P D FIG. 4. (Color online) Midgap state P as function as a func-tion of D , the distance of the vacancy to its closest nanoribbonedge (in units of a √ / The results of Fig. 3 are obtained from numerical di-agonalization with M = N = 40 (largest size in Fig. 4),and are consistent with the analytical approach. Figure 4shows that P becomes increasingly independent of sys-tem size as the ribbon width is increased. This behavioris further confirmed by the analysis of U presented inFig. 8.All these results are obtained for semiconducting arm-chair graphene nanoribbons. For metallic armchair rib-bons we also find vacancy-induced states which alsopresent the behavior observed in Figs. 3 and 4. Thisfinding is at odds with the analysis of Ref. 36. We willreturn to this discussion on Sec. V where we study thebehavior of U in armchair nanoribbons.For semiconducting graphene nanoribbons, our find-ings for P can be qualitatively understood as follows.The combination of the bipartite nature of the honey-comb lattice and the presence of edges gives origin toa peculiar modulation of the degree of localization as afunction of vacancy-edge distance. Some insight is pro-vided from recent studies of the Kitaev model in thegapped phase, where it has been shown that a vacancyinduces a zero mode state with a specific directionalitysuch that the wave function of this state is nonzero only ina wedge emanating from the vacancy position and zeroelsewhere. This directionality is only sublattice de-pendent, but independent from the site chosen withinthe same sublattice. We recall here that semiconduct-ing armchair ribbons constitute a realization of a gappedhoneycomb model. Hence, we also observe a direction-ality pattern in the vacancy zero mode states in theseribbons. We note that in distinction to the Kitaev model,the systems we study have edges that strongly influencethe behavior of the midgap states, particularly, their di-rectionality and degree of localization as a function of thevacancy site position.Further insight is obtained by inspecting the wavefunction with energy closest to zero of a clean armchairribbon, illustrated in Fig. 5. By inspecting the probabil-ity of finding the states at the sites along any of the “hor-izontal” zigzag chains starting from the armchair edge,see Fig. 5, one clearly observes a pattern of two maximafollowed by a minimum. We note that this pattern iscomplementary to the one we obtain for P as we movethe vacancy across a zigzag chain. Hence, if a vacancy isplaced at a site where the state closest to the zero modeof the clean ribbon has a maximum, the vacancy statesuffers a larger repulsion and has to extend more, theopposite occurs when the vacancy is placed at sites withvanishing weights in the clean system.This view is confirmed by the study of midgap statesdue to hydrogen adatoms (not shown here) instead of va-cancies. Interestingly, the manner in which the vacancyzero mode is more localized/extended is consequence ofthe directionality of the midgap state. In the sites where P has a maximum, the vacancy state is directed towardsthe closest edge becoming more concentrated. In thecases where a minimum in the IPR occurs, the vacancystate is directed to the farther edges or to the the ribbon’slongitudinal direction and hence becomes more extended.This is an interesting property which could allow thetunability of a site-site entanglement of this system andalso tune the transport properties of a specific edge dueto interactions with the vacancy state. The transportalong a chosen edge could be blocked or reduced just by selecting the vacancy/adatom position respective thechosen edge. As we show in Sec. V, all these findingshave a strong impact on the charging energy U . A r m c ha i r edge ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ...... ... a r m c h a i r e d g e FIG. 5. (Color online) Site dependency of the amplitudes ofthe state with energy closest to zero energy of a pristine arm-chair ribbon. The circles radii are proportional to the wavefunction amplitude, red and blue correspond the negative andpositive amplitudes respectively.
C. Vacancies with other kinds of edges: zigzag andquantum dots
For other kinds of edges, namely zigzag, chiral andthose found in quantum dots, we also observe vacancy-induced states | ψ (cid:105) .For zigzag ribbons, Ref. 36 predicts that a vacancydoes not affect the zero energy states of the system. Ourstudy agrees with this result. However, we find thatthe vacancy gives rise to pairs of eletron-hole states withsame energy and amplitude (but opposite phases), as ex-pected from chiral symmetry preservation. In addition, | ψ (cid:105) clearly hybridizes with the edge-localized states ofthe edge with sites belonging to the opposite sublatticeof the vacancy site (see Fig. 6a). These observations arein agreement with transport simulations for this kind ofsystem. The situation becomes more complicated in the caseof graphene quantum dots, whose edges are (in general)a combination of zigzag and armchair chains. We findthat the midgap states hybridize with the states localizedat the system edges. In addition, we observe an energyshift of the vacancy state, similarly to the case of zigzagribbons. Here the presence of the armchair-like edges alsoplays a role: We find a modulated behavior of the IPR ofthe vacancy state as the distance of the vacancy site fromthe armchair edge is varied. As in the case of armchairribbons, we also observe increasing degree of localizationof the vacancy state as the distance from the armchairedge is decreased. The IPR has its maximum for D = 2as for armchair ribbons, see Fig. 6b.Although we did not make an extensive numerical anal-ysis of vacancies in zigzag ribbons or graphene quantumdots, the discussion above is useful to emphasize that themain results of our paper for the armchair ribbons shouldremain valid for realistic samples in which other kinds ofedges and/or disorder appear. IPR=0.0093 (b) P = 0 . (a) IPR=0.03 (a) P = 0 . (b) FIG. 6. Wave function amplitudes for the vacancy state in azigzag nanoribbon (a) and for a quantum dot (b).
V. CHARGING ENERGY ESTIMATES FORARMCHAIR RIBBONS
We focus our attention on the analysis of vacancy-induced states in armchair graphene nanoribbons sincein this case the degree of localization of the | ψ (cid:105) statesis much more pronounced, as discussed in the previoussections. More specifically, we calculate the charging en- ergy U and U , given respectively by Eqs. (17) and (20),as a function of the vacancy-edge distance D and thenanoribbon width.Both the analytical and the tight-binding approachesgive the envelope wave functions ψ ( i ) which allows oneto numerically evaluate U and U . We compare the re-sults and find good qualitative agreement, as shown inFig. 7. For the cases we considered (with M = 220), thetight-binding values for U are about 30% smaller (for thesites with the larger U ) than those obtained using theanalytical wave functions. For U , which is one order ofmagnitude larger than U , both approaches agree within5%. In Appendix A, we provide evidence that the dis-crepancy is not related to the size of the supercell. Wespeculate that the larger discrepancy in the case of U arises from the fact that it scales with the numerically-obtained amplitudes as | ψ ( i ) | (see, e.g., Eq. (17)) while U scales roughly as | ψ ( i ) ψ ( j ) | with i (cid:54) = j (Eq. (20)).Thus, the numerical values of U tend to be more sen-sitive to small numerical errors in the calculation of theenvelope wave functions ψ ( i ) than those obtained for U . (a) (b) FIG. 7. (Color online) Charging energy U as a function ofthe vacancy-edge distance D for an armchair nanoribbon ofwidth M = 220. U = U + U is calculated by using numericaland analytical wave functions. Panel (a) corresponds to theintra-atomic U contribution, while (b) to the inter-atomicCoulomb term U . The insets show maxima of U and U versus D in a log-log scale. . The calculation of U using the analytical expressionfor the ψ ( i ) is significantly faster and demands much0less memory than the tight-binding approach, which re-quires storage of the Hamiltonian matrix elements. Thisallows us to address ribbons with micron size widths.The results presented next are obtained using analyticalwave functions.Figure 8 shows U and U as a function of the vacancy-edge distance D . The results clearly indicate that U isalmost independent of the ribbon width M for very largesystems. Moreover, for the sites for which U ( U ) hasa minimum, the computed charging energies approachthe bulk values obtained in Fig. 2. The dependence of U with D follows the same pattern as the IPR (for U this is expected, since the latter is proportional to P ),namely, an oscillatory behavior as the vacancy is movedaway from the edge with a fast decay in the modulation(see insets of Fig. 7). When the vacancy is far from theribbon edges, U is suppressed with respect to its largestvalue by an order of magnitude and approaches the bulkestimate. (a) (b) FIG. 8. (Color online) Charging energy U as a function of thevacancy-edge distance D for armchair nanoribbons of differentwidths, namely, M = 40 , , and 1000. The charging energyis split in an intra-atomic U (a) and an inter-atomic Coulombterm U (b). The insets show the tight-binding result for ametallic ribbon with M = 44. . We note that the ribbon with M = 44 is metallic. Inthis case, the analytical treatment of Ref. 36 predicts thatthe vacancy does not produce a localized state pinned atzero energy. Following the notation introduced in Ref. 36,we find that, when the vacancy is placed at the so-callednodal line (sites where the wave function of the lowest energy state of the clean system vanishes, which corre-spond to the sites with the smallest probabilities in Fig.5), a localized state pinned at zero energy arises, orig-inating the maxima in the charging energy seen in theinset of Fig. 8 for the width M = 44. Interestingly, ifthe vacancy is placed at a site off the nodal line, we donot observe a localized state pinned at zero energy, aspredicted in Ref. 36. Instead, we find that the vacancyoriginates an electron-hole pair of states close to zero en-ergy with a localized character. The charging energy forthese states correspond to the “minima” in the modu-lated pattern shown in the insets of Fig. 8. We notethat these estimates are in excellent agreement with thebehavior found for the semiconducting ribbons, indicat-ing that the values for U are robust irrespective to themetallic or semiconducting character of the ribbons.The appearance of localized states is metallic ribbons isimportant for the discussion of the observation of boundstates immersed in the continuum (BIC) in graphene as the localized states here occur within a region of finitedensity of states. VI. CONCLUSIONS AND OUTLOOK
In this paper we study the charging energy U , a keyelement to understand the magnetic properties of the sys-tem, of a localized state due to a single-vacancy in mono-layer graphene bulk and in graphene nanoribbons.We find that U can be expressed in terms of two maincontributions, U and U , corresponding respectively tointrasite and intersite electron-electron interactions. Weshow that U can be identified with the effective low en-ergy expression for the Coulomb energy associated withtwo-dimensional electronic wave functions. Although U is the dominating term, there are several scenarios wherethe U contribution to the charging energy can becomeimportant. For instance, a recent study suggests theuse of impurities to design a lattice structure that cangive rise to a graphene-based spin-liquid system. Thereit is argued that the onset of the spin-liquid regime de-pends on the ratio between the charging energy of theimpurity bound state and the hopping between these im-purity states, that makes the precise assessment of U very critical to infer the system behavior.Our systematic study of the charging energy confirmsthe heuristic prediction that U scales with the samplelength L as (ln L ) − . Our estimates for U support thepicture of π magnetism in realistic sample sizes, contraryto previous results. Edges change significantly this simple bulk scaling. Inthis paper, we perform a systematic investigation of themidgap states in armchair ribbons. We establish the de-pendence of U and the IPR with the vacancy-edge dis-tance and discuss how this is related to the directionalityof the wave functions. Some midgap states are truly lo-calized (unnoticed so far) with IPR and U an order ofmagnitude larger than the bulk value. Similar behav-1ior is also found for metallic armchair ribbons, whichcould be a manifestation of a BIC. This particular caseis at odds with the wave function analysis presented inRef. 36. For the remaining cases the overall agreementis very good. Moreover, we show that other kinds ofedges also affect the vacancy-induced state. In particular,for zig-zag edges and quantum dots, the vacancy-inducedstates have energies shifted from the Dirac point due to astrong hybridization with edge-induced localized states.This observation can be an indication that our resultsfor armchair edges can be robust in samples with edgesother than zigzag and armchair.Using the tight-binding model, we have checked thatan hydrogen adatom gives raise to localized states thatshare some of the features of those caused by a vacancy.Their IPR, for instance, have the same oscillatory behav-ior we find for the armchair ribbons with vacancies. Thisopens the possibility of having a controlled way to ver-ify our results by the precise manipulation of H adatomsusing a STM tip and to observe the generated states pat-tern by a STS measurement near the adatom, as recentlyshown in Ref. 31. The results we present can also bechecked experimentally in other platforms since they arevalid for other systems with bipartite lattices. Artificialgraphene could also be used to verify experimentally ourfindings.We stress that a correct assessment of the substratedielectric constant (cid:15) is key for a quantitative comparisonof our U estimates with the experimental values. We re-call that the estimates we present here were obtainedconsidering (cid:15) = 4 which is characteristic of grapheneon top of SiO . Depending on the dielectric environ-ment to which graphene is exposed (cid:15) can vary by twoorders of magnitude. Hence, the dielectric media towhich graphene is submitted in the recents experimentsof Refs. 14 and 31, which use graphene on top of SiCand Rh substrates, respectively, can have an importantinfluence in the discrepancies between our estimates andthe observed U . A systematic experimental study of therole of the substrate and proper characterization of thedielectric constant that should be used in the U estimatesis necessary to clarify this issue.Our U estimates can be directly comparable to thatof Ref. 33, which also estimates U for graphene on SiO substrates. We recall that our results are 2 − U , whose value iscomputed using the (numerically precise) vacancy-statewave function amplitudes. By contrast, the values re-ported in Ref. 33 are based on scaling arguments forthe vacancy-induced wave function, which do not takeproperly into account the behavior of the vacancy wavefunction amplitudes along the sites of the system.Finally, we note that though the results presented for U are for single vacancies, the expression we derive isalso valid for the multivacancy case. Though an analyt-ical expression for the multivacancy wave function case is not available in the literature, those can be obtainednumerically and used to obtain U in a similar fashion asthe one we use here for the single vacancy.In summary, we present a full derivation of U forvacancy-induced localized states in graphene systems.Our results help the understanding of defect-induced car-bon magnetism, allowing contact with the typical exper-iments that use samples with billions of atoms, whereedges, multivacancies and other kinds of disorder arepresent, which are beyond the reach of other methods,such as DFT ACKNOWLEDGMENTS
This work has been supported by the Brazilian fundingagencies CAPES, CNPq, and FAPERJ.
Appendix A: Convergence analysis
The tight-binding model calculations are implementedusing a supercell of size M × N . In the case of graphenenanoribbons, M is determined by the ribbon width, whilethe longitudinal length N has to be conveniently chosensince the system is periodic in this direction. We optimizeour calculations by choosing the smallest value of N forwhich U becomes independent (within less than 1%) of N . U ( e V ) D N = = = = = (a) U ( e V ) D N = = = = = (a) (b) FIG. 9. (Color online) Charging energy as a function of dis-tance for armchair ribbons with length M = 40 and withvarying ”infinite length” N . (a) Results obtained from theanalytical expression of the vacancy wave function. (b) Re-sults from the tight binding model. In Fig. 9 we display the behavior of the U as a func-2tion of edge distance D for M = 40 and varying N . Theresults are obtained by using the analytical wave func-tion (Fig. 9a) and the tight binding model (Fig. 9b).For M = 40 we perform tight-binding simulations for10 ≤ N ≤ N andshow that the results already converge for modest valuesof N . These calculations also indicate that the discrep-ancies observed in the estimates of U obtained from thetwo approaches is not an artifact of the supercell sizes weuse. A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S.Novoselov, and A. K. Geim, Rev. Mod. Phys. , 109(2009). A. K. Geim and K. S. Novoselov, Nat. Mater. , 183 (2007). K. S. Novoselov, V. I. Fal’ko, L. Colombo, P. R. Gellert,M. G. Schwab, and K. Kim, Nature , 192 (2012). J.-H. Chen, W. G. Cullen, C. Jang, M. S. Fuhrer, andE. D. Williams, Phys. Rev. Lett. , 236805 (2009). J.-H. Chen, L. Li, W. G. Cullen, E. D. Williams, and M. S.Fuhrer, Nat. Phys. , 535 (2011). R. R. Nair, M. Sepioni, I.-L. Tsai, O. Lehtinen,J. Keinonen, A. V. Krasheninnikov, T. Thomson, A. K.Geim, and I. V. Grigorieva, Nat. Phys. , 199 (2012). R. R. Nair, I.-L. Tsai, M. Sepioni, O. Lehtinen,J. Keinonen, A. V. Krasheninnikov, A. H. Castro Neto,M. I. Katsnelson, a. K. Geim, and I. V. Grigorieva, Nat.Commun. , 2010 (2013). S. Just, S. Zimmermann, V. Kataev, B. B¨uchner,M. Pratzer, and M. Morgenstern, Phys. Rev. B , 125449(2014). V. W. Brar, R. Decker, H.-M. Solowan, Y. Wang,L. Maserati, K. T. Chan, H. Lee, O. Girit, A. Zettl, S. G.Louie, M. L. Cohen, and M. F. Crommie, Nat. Phys. ,43 (2011). M. M. Ugeda, I. Brihuega, F. Guinea, and J. M. G´omez-Rodr´ıguez, Phys. Rev. Lett. , 096804 (2010). M. M. Ugeda, D. Fern´andez-Torre, I. Brihuega, P. Pou,A. J. Mart´ınez-Galera, R. P´erez, and J. M. G´omez-Rodr´ıguez, Phys. Rev. Lett. , 116803 (2011). J. Mao, Y. Jiang, D. Moldovan, G. Li, K. Watanabe,T. Tanigushi, M. R. Masir, F. M. Peeters, and E. Y. An-drei, Nat. Phys. , 1745 (2016). G. Lopez-Polin, C. Gomez-Navarro, V. Parente, F. Guinea,M. I. Katsnelson, F. Perez-Murano, and J. Gomez-Herrero, Nat. Phys. , 26 (2015). Y. Zhang, S.-Y. Li, W.-T. Li, J.-B. Qiao, W.-X. Wang,L.-J. Yin, and L. He, arXiv:1604.06542 (2016). O. V. Yazyev, Rep. Prog. Phys. , 56501 (2010). V. G. Miranda, L. G. G. V. Dias da Silva, and C. H.Lewenkopf, Phys. Rev. B , 201101 (2014). L. Fritz and M. Vojta, Rep. Prog. Phys. , 032501 (2012). V. M. Pereira, F. Guinea, J. M. B. Lopes dos Santos,N. M. R. Peres, and A. H. Castro Neto, Phys. Rev. Lett. , 036801 (2006). V. M. Pereira, J. M. B. Lopes dos Santos, and A. H. CastroNeto, Phys. Rev. B , 115109 (2008). S. Yuan, H. De Raedt, and M. I. Katsnelson, Phys. Rev.B , 115448 (2010). M. Inui, S. A. Trugman, and E. Abrahams, Phys. Rev. B , 3190 (1994). G. Santhosh, V. Sreenath, A. Lakshminarayan, andR. Narayanan, Phys. Rev. B , 054204 (2012). A. J. Willans, J. T. Chalker, and R. Moessner, Phys. Rev.B , 115146 (2011). This result holds if the A and B sublattice constituentsare of the same kind. If this is not the case, the states arenot pinned at E = 0 . It worth noting that vacancy zeromodes in graphene can also be generated even if N A = N B as long as the system is composed of clusters for which N A (cid:54) = N B inside the cluster . J. J. Palacios and F. Yndur´ain, Phys. Rev. B , 245443(2012). C.-C. Lee, Y. Yamada-Takamura, and T. Ozaki, Phys.Rev. B , 014401 (2014). B. R. K. Nanda, M. Sherafati, Z. S. Popovi´c, and S. Sat-pathy, New J. Phys. , 83004 (2012). M. Casartelli, S. Casolo, G. F. Tantardini, and R. Marti-nazzo, Phys. Rev. B , 195424 (2013). P. Francis, C. Majumder, and S. Ghaisas, Carbon , 358(2015). W. S. Paz, W. L. Scopel, and J. C. C. Freitas, Solid StateCommun.
175 - 176 , 71 (2013). H. Gonzalez-Herrero, J. M. Gomez-Rodriguez, P. Mallet,M. Moaied, J. J. Palacios, C. Salgado, M. M. Ugeda, J.-Y.Veuillen, F. Yndurain, and I. Brihuega, Science , 437(2016). E. H. Lieb, Phys. Rev. Lett. , 1201 (1989). M. A. Cazalilla, A. Iucci, F. Guinea, and A. H. CastroNeto, cond-mat/1207.3135 (2012). X. Dou, V. N. Kotov, and B. Uchoa, arXiv:1602.01477(2016). H.-Y. Deng and K. Wakabayashi, Phys. Rev. B , 115413(2014). H.-Y. Deng and K. Wakabayashi, Phys. Rev. B , 035425(2015). J. J. Palacios, J. Fern´andez-Rossier, and L. Brey, Phys.Rev. B , 195428 (2008). L. Brey and H. A. Fertig, Phys. Rev. B , 235411 (2006). H. Zheng, Z. F. Wang, T. Luo, Q. W. Shi, and J. Chen,Phys. Rev. B , 165414 (2007). E. R. Mucciolo, A. H. Castro Neto, and C. H. Lewenkopf,Phys. Rev. B , 075407 (2009). A. R. Carvalho, J. H. Warnes, and C. H. Lewenkopf, Phys.Rev. B , 245444 (2014). C. Tao, L. Jiao, O. V. Yazyev, Y.-C. Chen, J. Feng,X. Zhang, R. B. Capaz, J. M. Tour, A. Zettl, S. G. Louie,H. Dai, and M. F. Crommie, Nat. Phys. , 616 (2011). G. Z. Magda, X. Jin, I. Hagymasi, P. Vancso, Z. Osvath,P. Nemes-Incze, C. Hwang, L. P. Biro, and L. Tapaszto,Nature , 608 (2014). S. Wang, L. Talirz, C. A. Pignedoli, X. Feng, K. Muellen,R. Fasel, and P. Ruffieux, ArXiv e-prints (2015),arXiv:1511.04940 [cond-mat.mtrl-sci]. D. A. Bahamon, A. L. C. Pereira, and P. A. Schulz, Phys. Rev. B , 165438 (2010). A. Orlof, J. Ruseckas, and I. V. Zozoulenko, Phys. Rev.B , 125409 (2013). N. Gorjizadeh, A. A. Farajian, and Y. Kawazoe, Nan-otechnology , 15201 (2009). D. Midtvedt and A. Croy, J. Phys.: Condens. Matter ,045302 (2016). M. Sch¨uler, M. R¨osner, T. O. Wehling, A. I. Lichtenstein,and M. I. Katsnelson, Phys. Rev. Lett. , 036601 (2013). E. N. Economou, Green’s Functions in Quantum Physics,3rd ed. (Springer, Berlin, 2006). E. Kaxiras, Atomic and Electronic Structure of Solids, 1sted. (Cambridge University Press, Cambridge, 2003). A. Zunger and R. Englman, Phys. Rev. B , 642 (1978). W. Ku, T. Berlijn, and C.-C. Lee, Phys. Rev. Lett. ,216401 (2010). M. E. A. Coury, S. L. Dudarev, W. M. C. Foulkes, A. P.Horsfield, P.-W. Ma, and J. S. Spencer, Phys. Rev. B ,075101 (2016). W. A. Harrison, Phys. Rev. B , 2121 (1985). C. M. Goringe, D. R. Bowler, and E. Hern´andez, Rep.Prog. Phys. , 1447 (1997). M. Elstner, D. Porezag, G. Jungnickel, J. Elsner,M. Haugk, T. Frauenheim, S. Suhai, and G. Seifert, Phys.Rev. B , 7260 (1998). R. G. Parr, D. P. Craig, and I. G. Ross, J. Chem. Phys. , 1561 (1950). T. O. Wehling, E. S¸a¸sioˇglu, C. Friedrich, A. I. Lichtenstein, M. I. Katsnelson, and S. Bl¨ugel, Phys. Rev. Lett. ,236805 (2011). Our approximation is also consistent with the ratios be-tween U , U , and U (bare) reported in Ref. 59. C.-P. Lu, G. Li, K. Watanabe, T. Taniguchi, and E. Y.Andrei, Phys. Rev. Lett. , 156804 (2014). D. Soriano, D. Van Tuan, S. M-M Dubois, M. Gmitra,A. W. Cummings, D. Kochan, F. Ortmann, J.-C. Charlier,J. Fabian, and S. Roche, 2D Materials , 022002 (2015). On the other hand, we find that when the vacancy stateis moved in the direction parallel the edges,the degree oflocalization is unchanged, which is quite reasonable sincethe system is “infinite” in this direction. The latter states were probably the states observed inRefs. 37 and 48 which led them to conclude that goingfrom a very narrow ribbon to bulk graphene would turnthe localized state into a non-normalizable one with van-ishing spin splitting for the midgap state. Hence, a vanish-ing magnetization. J. W. Gonz´alez, M. Pacheco, L. Rosales, and P. A. Orel-lana, EPL (Europhysics Letters) , 66001 (2010). N. Cort´es, L. Chico, M. Pacheco, L. Rosales, and P. A.Orellana, EPL (Europhysics Letters) , 46008 (2014). V. N. Kotov, B. Uchoa, V. M. Pereira, F. Guinea, andA. H. Castro Neto, Rev. Mod. Phys. , 1067 (2012).68