{'fraction_non_alphanumeric': 0.03596044351213665, 'fraction_numerical': 0.01584923251090467, 'mean_word_length': 4.051976450798991, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Electronic structure of vacancy resonant states in graphene: a critical review of the single vacancy case

The resonant behaviour of vacancy states in graphene is well-known but some ambiguities remain concerning in particular the nature of the so-called zero energy modes. Other points are not completely elucidated in the case of low but finite vacancy concentration. In this article we concentrate on the case of vacancies described within the usual tight-binding approximation. More precisely we discuss the case of a single vacancy or of a finite number of vacancies in a finite or infinite system.


I. INTRODUCTION
Vacancies and various point defects in graphene have been the object of numerous studies both theoretical and experimental. It is of course of utmost interest to understand the role of impurities on the electronic structure of graphene. From the experimental side, let us just mention the problem of doping with nitrogen (or boron) impurities, 1-3 the influence of hydrogen adatoms, 4,5 or the controversial problem of magnetism induced by defects. 6,7 From a theoretical point of view it has been recognized fairly early that the particular electronic structure of graphene with its vanishing density of states at the Dirac points provokes important resonant effects. 8 The case of vacancies is particularly interesting since in the simplest tight-binding model with electron-hole symmetry, resonances occur just at these points. This has been related to the occurrence of so-called zero energy modes where the corresponding wave functions are concentrated on the sublattice different from that of the vacancy. [9][10][11][12] The case of a single impurity is rather well understood, 9,13-15 but some ambiguities remain 16 and other points are not completely elucidated in the case of low but finite vacancy concentration, since the usual tools such as the average t-matrix approximation (ATA) or the coherent potential approximation (CPA) are not accurate enough. [17][18][19] From the numerical side huge systems (more than 10 6 atoms) are necessary to hope to obtain convergence, 20,21 which is out of reach of ab initio calculations. 22 In this article we concentrate on the case of vacancies described within the usual tight-binding approximation.
We first recall this model as well as the properties of the Green functions which will be used at length afterwards. Then, we discuss the case of a single impurity in a large but finite system. In the presence of a vacancy, one state is substracted from the total density of states. A simple algebraic argument shows on the other hand that a "quasi-localized" mode of zero energy should necessary appear simultaneously. Both effects are clearly associated but we will show that this statement hides some subtleties related to the boundary conditions, which accounts for the apparently different results obtained in the literature. In the limit of infinite systems, scattering boundary conditions and the Green function formalism are appropriate, which allows us to discuss in detail the behaviour of the local and total densities of states.
The case of a finite number of vacancies or impurities can be studied quite similarly, but the thermodynamic limit where the number of impurities is infinite and the concentration c is finite is much more difficult to handle, even in the c → 0 limit. This will be discussed elsewhere.

A. The tight-binding model
We use the simplest tight-binding hamiltonian with single transfer integrals −t connecting the first neighbour π orbitals of the graphene structure: 23,24 where |n denotes the π state at site n and where the prime indicates a sum over first neighbours n and m.
The notations used to describe the unit cell and the Brillouin zone are shown in Fig. 1 . The corresponding Bloch functions and eigenvalues are:  where N is the number of units cells. 25 At low energy E = ±νq ; ν = 3a cc t/2, where a cc is the first neighbour carbon-carbon distance and k = K + q or K + q. With our conventions for the δ and K vectors, it can be checked that K α .δ β = (β − α) 2π 3 mod 2π, and then exp(iK α .δ β ) = j β−α where j = e 2iπ/3 .

B. Green functions
The Green function or resolvent is defined from G 0 (z) = (z − H 0 ) −1 , where z is a complex variable, and its matrix elements in the tight-binding basis G nm (z) = n|G 0 (z)|m are related to the local densities of states n nm (E) through: Formally G 0 nm (z) can be expanded in successive powers of 1/z. Actually this expansion is convergent for large enough |z| : The matrix elements of (H 0 ) p are related to paths involving p jumps between first neighbours (moments) and are real. G 0 nm (z) only depends on m − n and is complex through z. As a consequence: Furthermore, as in all alternant lattices with first neighbour coupling, we know that in the moment expansion, there are only even or odd paths depending on sites n and m belonging to the same sublattice or not. As usual we introduce the notations G 0 i j nm , i, j = A, B to specify, when necessary, the sublattice A or B of the sites. Finally this implies the following symmetry properties: where F nm (E) is the real part of G 0 nm (E), so that using Eq. (5), we obtain: and opposite relations for the densities of states:

Explicit expressions
The subject is well documented, 11,12,[25][26][27][28][29] but it is useful here to summarize the results. 30 Explicit expressions for the Green functions are obtained using the Bloch basis: Hence the familiar 2 × 2 matrix representation: Actually the corresponding integral in Eq. (9) can be expressed exactly in terms of elliptic functions and furthermore Horiguchi 31 has shown that all matrix elements can be obtained in terms of a few of them through recurrence relations. This is not necessarily the best procedure from a numerical point of view 28,32 and standard integration or recursion methods are generally more efficient. In practice, since we are principally interested in the low energy regime, it is very fruitful to obtain simplified expressions which allow us to understand general trends. The diagonal matrix element, independent of n plays an important part and is given by: The imaginary part provides us with the density of states per atom, equal to n 0 (E) ≡ n 0 00 (E) = |E|/W 2 ; W 2 = √ 3πt 2 , which is an exact result to first order in energy. The calculation of the full Green function is more delicate since its real part is the Hilbert transform of the density of states and involves high energy states where the linear approximation for the dispersion relation is no longer valid. The recipe is known and consists in introducing a cut-off Λ in k-space. Equivalently we introduce a cut-off in energy which turns out to be equal to W if we want to keep a normalized density of states. Finally: where the proper determination of the logarithm in the complex plane should be taken for Im z > 0, with a cut along the real axis. The approximation is known to be valid for energies lower than t (see e.g. Ref. [19] and Fig. 2). We will denote it in the following as the Debye-like approximation. The other matrix elements can be obtained with similar techniques: the imaginary parts can be calculated exactly in the limit E → 0 in n (E) terms of Bessel functions, and the full functions depend in principle on a cut-off in k-space Λ = W/ν around K and K . Here however, this cut-off is no longer necessary to achieve convergence of the integrals which finally reduce to Hankel functions in the limit Λ → ∞. For example G 0AA 0n (z), n = 0, is obtained from: 1 (nz/ν) cos(K.n − ω n ) .
Hankel functions are defined in the upper part of the complex plane as usual, from analytic continuation of their values on the positive real axis: H 0,1 (z); ω n is the angle between n and the Ox axis which is, according to our convention, the armchair direction. Notice that this precise form of the argument of the cosine depends on our choice of the K point because in the case of the G 0 AB Green function, n is not a lattice vector.

Low energy limit E → 0
Since the Bessel functions depend on the product E n, the E → 0 and n → ∞ limits do not commute. We consider here the E → 0 limit. Using the standard properties of Bessel functions in the limit x → 0 : where γ is the Euler constant, we obtain: Notice finally that one cannot simply set n = 0 in G 0AA 0n (z) to obtain the diagonal Green fonction. The high energy cut-off corresponds to a low distance cut-off about 0.7 a cc . All these results are in full agreement with previous estimates. 12 The behaviour of the Green functions at the origin is therefore different from that at the usual van Hove singularities (discontinuity or logarithmic singularity). 33 In particular all G 0AA 0n (E) vanish, whereas all G 0AB 0n (E) behave as cos(K.n − ω n )/n and are finite, because the cosine never vanishes (see below).

Large distance limit n → ∞
The large distance behaviour is governed by the Hankel functions: where here a = √ 3a cc is the length of unit cell vectors a 1,2 . In fact this limit can be obtained directly using stationary phase arguments 34 and is practically exact at low -but finite -energy.

Symmetry properties
The approximate expressions for the Green functions have the advantage to factorize the dependence on distance. The Bessel functions depend only on the modulus n of n whereas the spatial symmetry is governed by the cosine factors.
a. G 0AA Green functions.
For G 0AA 0n , this factor is equal to cos K.n which has the symmetry of the socalled √ 3× √ 3 superstructure. It takes values 1 and −1/2 depending on the sites belonging to the superlattice or to its motif, respectively. Globally, G 0AA 0n has the full C 6v point symmetry of the hexagonal lattice (Fig. 3). b. G 0AB Green functions.
The case of the Green functions connecting different A and B sublattices is more subtle. Any vector n in this case can be decomposed into a vector of the √ 3 × √ 3 structure and a vector of type δ α . Consider first the atoms obtained from δ 1 along Ox in the negative direction. For these atoms K.n = 0 mod 2π and the geometrical factor is simply cos ω n which takes values ±1 on the 0x axis. Its absolute value decreases when deviating from this axis. In the zigzag direction the angle is never equal to π/2 but reaches this value asymptotically (Fig. 4). The other atoms are simply obtained by applying a three-fold symmetry about the origin, which yields the map shown in Fig. 5.
It is interesting to notice that all atoms in the armchair direction along Ox are already obtained with the single family shown in Fig. 4 so that |Re G 0AB 0n (E)| ∼ 1/n in this case without any oscillatory behaviour at low energy. On the other hand, in the Oy zig-zag direction the Green function is very weak for the family of sites shown in Fig.  4 contrary to the two other ones obtained from ±2π/3 is equal to cos ωn for the atoms shifted from the √ 3 × √ 3 superstructure by the vector δ1. Its value is indicated close to these atoms. The size of the atoms is approximately proportional to the absolute value of this factor.
rotations. This has been studied in more detail in Refs [12] and [29].

Discussion
As pointed above, an important result is that the only non-vanishing Green functions when E → 0 are the real parts of the off-diagonal matrix elements connecting the two different A and B sublattices. Since the real parts are less accurate than the imaginary parts, it is useful to test the accuracy of the previous formulae. Precisely at E = 0 it is in fact possible to account analytically for the cut-off W with the result: Since W/ν 1.56/a cc the correction decreases rapidly with n but is non negligible for the first neighbours for which n = a cc . Let G 0 01 denotes the first neighbour Green function, the previous formula gives G 0 01 (E = 0) 0.29/t whereas the asymptotic formula for an infinite cut-off gives 0.55/t. The exact result is 1/3t and can be deduced from the exact relations obtained from the identity (z − H 0 ) G 0 = 1. Taking the diagonal matrix elements of both sides yields indeed: G 0 01 = (1 − zG 0 00 )/3t = 1/3t when z = 0. In Fig. 6 we compare the "exact" result obtained from the recursion method with the Bessel function approximation; the result is qualitatively good if not quantitatively, but would be improved using a cut-off as mentioned above (see also Ref. [26]). 35 Taking then the matrix elements between first neighbours we obtain also an exact expression for the second neighbour Green function G 0 02 = [(z 2 − 3t 2 )G 0 00 − z]/6t 2 . In Fig. 7 we show the real parts of the first G 0AB 0n Green functions using Bessel functions. By comparison with exact numerical calculations, it can be checked that the larger n, the better the accuracy of the Bessel function approximation.
It should be pointed out also that at finite energy, the Bessel functions begin to oscillate as soon as nE/ν is larger than about 2. If energies are measured in units of t and distances in units of a cc , this means that the low energy expansion of the Green functions G 0 0n (E) are valid when E < 1/n.

III. ELECTRONIC STRCTURE OF THE VACANCY
A. Zero modes or not ?
The existence of particular states of zero energy, the so-called zero modes, has been predicted using fairly simple algebraic arguments which apply to alternant lattices described within a first neighbour tight-binding approximation. 9 The vacancy here is modelled by removing the corresponding state at the vacancy site taken at the origin (A sublattice). Although the single vacancy is known to display a slight reconstruction 36,37 the model remains useful and applies to other situations (hydrogen adatom, pyridine-like nitrogen defect). Removing an atom can be made either explicitly by introducing the hamiltonian H in the new basis: or by introducing a local potential V = |0 U 0| and taking the limit U → ∞, which forces the wave function to vanish on site 0: The only difference is that in the latter case, one bound state is repelled to infinity. Otherwise the two approaches are equivalent, but the latter is known to be more easy to handle. It will be clear in the context what precise hamiltonian is used. Now, starting from pristine graphene with an equal number N of A and B sites, the energy spectrum is symmetric about the origin. In the presence of a vacancy the hamiltonian applies the space of B states of dimension N on the space of A states of dimension N − 1, and there is necessarily a non trivial linear combinaison of B states whose image vanishes. This supplementary state is supposed to be at the origin of the resonance found in the density of states. Actually it is always dangerous to apply theorems valid for finite systems to infinite systems or systems with special boundary conditions, and this is exactly what happens here as we show now.

Vacancy in a finite system
The simplest way to build a perfect finite graphene structure is to introduce periodic boundary conditions. This is possible in many different ways as we have learnt in the case of carbon nanotubes where any non trivial lattice vector can define a wrapping vector C = pa 1 + qa 2 where p and q are integers. A second independent vector defines then a torus. Nanotubes are metallic or semiconductor depending on p − q being a multiple of 3 or not. This depends on K being an allowed wave vector for the Bloch functions or not. This can be extended to a torus. For simplicity we will consider p × p tori in the following. Then K is an allowed wave vector if p is a multiple of 3. It is instructive to look more precisely at the Bloch fonctions |ψ at the K and K points. We have four possibilities: |K A(B) and |K A(B) . The corresponding values of the amplitudes n|ψ are the cubic roots of unity, 1, j = exp(2iπ/3) and j 2 = exp(−2iπ/3) which gives four possibilities indeed. A given state has its amplitudes on a single sublattice and amplitudes rotating with the sequence (1, j, j 2 ) or (1, j 2 , j) on every triangle of the sublattices (Fig. 8). In real space the Bloch states live therefore on the √ 3 × √ 3 structure. The wrapping vectors must therefore be consistent with this periodicity and we recover the fact that p should be a multiple of 3.
a. Metallic case To summarize, in the case of metallic tori, there are exactly four states of zero energy. Otherwise (semiconductor tori) there is no state of vanishing energy. This contradicts the argument sometimes put forward that zero modes only appear in the case where the numbers of A and B sites are different.
What happens then in the presence of a vacancy? Of course in the case of a metallic torus, states on the B sublattice are not perturbed and remain at zero energy. But there is another possibility frequently forgotten, which is to consider the linear combination |K A − |K A , which also vanishes on the vacancy site. Therefore three states survive at zero energy. And there is no other possibility.
In the presence of a localized potential actually, the tight-binding problem is basically exactly solvable for finite as well as for infinite systems. 38,39 The eigenenergies are obtained from the equation: where here G 0 00 (E) is the real Green function for the finite system considered: where the |α and ε α are the (discrete) eigenstates and eigenenergies of the unperturbed system. In the present case the unperturbed eigenstates are Bloch functions and all matrix elements 0|α are equal to 1/ √ N and the solution of Eq. (21) is obtained through the usual procedure of finding the intersections of G 0 00 (E) with 1/U . Because of the divergences at each ε α the new eigenvalues are found to alternate with the unperturbed ones. There are problems in the case of degeneracies but then it is sufficient to imagine that weak appropriate perturbations lift the degeneracies, so that a level with degeneracy g will give rise to g − 1 levels at the same energy, as we found above. Otherwise, the levels alternate strictly. In the case of the vacancy, U → ∞ and the new eigenvalues are the zeros of G 0 00 (E). The conclusion is that in the case of metallic tori, there are three and just three delocalized states of zero energy.
b. Semiconductor case The case of semiconductor tori is different. The K points are no longer allowed and there is a small but finite gap of order 1/L where L is the linear size of the system. By symmetry however G 0 00 (E = 0) = 0 so that there is a single zero mode in that case. Its wave function has been derived by Pereira et al. 40 by combining localized edge states. There is actually a direct and simple proof (see also Refs. [10 and 12]). Consider the state G 0 (z = 0)|0 and let us apply the vacancy hamiltonian H = Q 0 H 0 Q 0 on this state, using the fact that, by definition H 0 G 0 (z = 0) = −1. We have: The first term in the right hand side vanishes indeed by definition. As for the second one, it is proportional to G 0 00 (E) which vanishes also when E = 0. The origin being on a A sublattice, the considered state vanishes on this sublattice because of the symmetry properties of the Green functions. The components ψ n = n|ψ of this state are just the matrix elements of the bare Green function G 0 0n (E = 0). For a large system the sum defining it can be safely replaced by an integral. We can use the results of the previous section so that: when n belongs to the B sublattice, and vanishes on the A sublattice. In the continuous limit, this becomes ψ(r) ∼ e iK.r /(x + iy) + complex conjugate, which is in complete agreement with previous estimates if appropriate correspondences for the coordinates and the K points are made. 9,40 To normalize this state, we calculate n |ψ n | 2 = 0|(G 0 ) 2 |0 = −dG 0 00 (E)/dE for E = 0, which is of order ln N for a (finite) large system, as noticed by many authors. In the present case this can be shown using 0| c. Discussion One conclusion of our discussion is that there is well defined quasi localized zero mode associated to the vacancy if the initial graphene system has no zero mode already, and this depends on the precise boundary conditions used. Quasi localized and itinerant states as well as more or less hybridized states can coexist close to zero energy. This accounts for some previous puzzling results. For example Pereira et al.. 9 reported a first localized state at zero energy and a second state completely itinerant, which corresponds to our "semiconductor" case. Conversely, Huang et al.. used "metallic" boundary conditions and did found itinerant zero modes. 16 Evidently, there are cases where a genuine gap opens: nanoribbons, 41 semiconducting nanotubes, in which case the occurence of a true bound state is clear.

Finite number of vacancies in a finite system
We consider here semiconducting tori. Since the vacancy state G 0 (E = 0)|0 has no amplitude on the A sites, it is not sensitive to other vacancies on this sublattice. Finally with a finite set of vacancies at sites n v is associated a finite number of zero modes G 0 (E = 0)|n v . This can also be seen by treating the case of a second vacancy at site A2 as a perturbation of the case of a single vacancy at site A1. We have again to solve an equation of type α | A2|α | 2 /(E − ε α ) = 0 where now α labels states in the presence of the vacancy A1. There is then a state of zero energy, but its weight | A2|α | 2 vanishes so that this function has still a zero (and not a divergence) at the origin. If instead the second vacancy is on a B site the weight no longer vanishes and the zero mode disappears.

3.
Vacancy in an infinite system From the above discussion it is clear that it is preferable to use a formalism that does not depend on the boundary conditions. The usual method for that introduced in Section II B is to add imaginary parts E → E + iε and calculate density of states variations by taking the limit ε → 0 after taking the limit N → ∞. In the case of a localized potential the exact expression of the Green function is known: where τ is the t-matrix. In the case of a vacancy, in the limit U → ∞ this reduces to: The variation δn n (E) of the local density of states at site n is then given by: Using the results of Section II B for the Green functions , δG nn is then found to vanish at E = 0 on A sites whereas it diverges on the B sites. For these sites the variation of the density of states is given by: and is shown in Fig. 9 in the case of first neighbours, n = a cc . The full density of states n 1 (E) is shown in Fig. 10. The number of states between E = −t and E = +t is about 0.2, but is only equal to 0.08 when the integration is made between -0.1 and 0.1. These states are substracted from the wings of the unperturbed density of states. The corresponding δn n (E) for A sites is quasi negligible, of the order of 0.05 t −1 for the second neighbours. Up to a multiplicative factor the curves are similar for all B sites, for the close neighbours at least. At larger distances the oscillations due to G 0AB 0n appear. As a result interferences occur when summing all contributions. We have actually an exact formula for the total variation: For the integrated density of states N (E) = E dE n(E ), this gives: This is a typical phase shift formula which can be derived quite generally by starting from the identity: 39 N (E) = Im π Tr ln G = Im π ln det G .
The behaviour of arg(G 0 00 ) is particularly simple to study if we use an Argand plot in the plane (Re G 0 00 , Im G 0 00 ) and look at the trajectory of G 0 00 as a function of E. It is more convenient to work with −G 0 00 which has a positive imaginary part so that the angle θ = arg(−G 0 00 ) belongs to the interval [0, π] and δN (E) = −θ/π. The usual definition for the phase shift η is η = arg(1/G 0 00 ) = −θ. It is clear from Fig. 11 that since G 0 00 vanishes at the origin, θ should jump from π to zero. Actually, because of the weakness of the logarithmic singularity of the real part of the Green function, the jump is limited by the value of ε, here equal to 10 −4 and by the resolution of the figure.
The corresponding variations of δN (E) and of its derivative δn(E) are shown in Fig. 12. The latter variation has therefore a δ-function contribution of weight one and also a singularity close to E = 0. This can be shown directly since d/dz(G 0 00 (z)) = 1/z +2/(W 2 G 0 00 (z)) so that, close to E = 0: δN (E) has also particular behaviours close to van Hove singularities which are not well reproduced by our model but which can be qualitatively understood using Argand plots. The striking point is that since δN (E) → −1 when E → ∞ two "states" are substracted from the initial system, one state when E < 0 and another one when E > 0. One state corresponds to the vacancy missing state (in the model with a local potential U → ∞ this state has been pushed to infinity) and the other one to the δ-function. What is surprising here is that this δfunction does not correspond to a single well-defined zero mode but is the result of interferences between many resonant states close to E = 0. The fact that the δfunction has weight one exactly is related to the fact that G 0 00 (z) vanishes at z = 0 and that its real part has a logarithmic singularity (infinite slope). Actually as is clear in Fig. 12 the practical weight is lower than one as soon as we integrate in a small window around E = 0. Introducing a small imaginary part amounts in collision theory to consider wave packets of extension ∼ 1/ε (the Green functions decay exponentially on this scale), very large, but much smaller that the overall size of the system. The wave functions can then be calculated using the Lippmann-Schwinger equation as noticed by Nanda et al. 12 Although not valid strictly at E = 0 it can be used in the limit E → 0 to show that all states there are strongly hybridized with the zero mode G 0 |0 .

Application: spatial variation of the electronic density
The electronic density at point r is given by n(r, E) = r| δ(E − H)|r . Projecting on the π atomic basis, we obtain: where φ(r − n) is the π orbital at site n. The main contribution to the sum is provided by the diagonal terms but interference terms related to the overlap between neighbouring orbitals are not always negligible.
In the case of the vacancy, we have an exact expression for G nm (E) (see Eq. (26)). Close to zero energy the unperturbed electronic density vanishes so that the main contribution to n(r, E) is the diagonal term proportional to (G 0AB 0n (E = 0)) 2 (see Eq. (28)). As a result the spatial variation of n(r, E 0) is driven by the behaviour of the Green functions. As mentioned previously, in the "Bessel" approximation these functions contain a geometrical factor shown in Fig. 5 and a spatial dependence ∼ 1/n. The result is shown in Fig. 13. For clarity the atomc orbitals have just been replaced by s-like functions. The resulting electronic density then mimics what can be observed for example in STM observations at a finite height above the graphene plane. The vacancy (A site) is at the center of the central triangle. The electronic density is concentrated on B sites. The image is characterized by the presence of two bright "triangles" and three arms formed by losanges. Three isolated dots can also be viewed as well as extinctions due to the geometrical factor cos(K.n − ω n ). This type of image is actually typical of that produced by a local resonant defect. 15,42 From the experimental side such images have been observed in nitrogen doped samples. 1,3 They can be due to substitutional nitrogen atoms and/or to so-called pyridine defects. They have also been clearly identified in STM observations. Notice that this typical contrast is not what is frequently mentioned as a √ 3 × √ 3 contrast, which has been observed in carbon nanotubes [43][44][45] as well as in graphene in the presence of extended defects. 46,47 FIG. 13. Local density of states in the presence of a vacancy (courtesy of H.Amara).

Finite number of vacancies
Vacancies on the same sublattice do not interact so that the difficult situations are those of "compensated" systems where the number of vacancies is the same on both sublattices. The t-matrix becomes here a genuine matrix in the space defined by the set of states on the vacancy sites. The label c will denote sites of the vacancy cluster C and P C is the projector on these sites, P C = c |c c|. Introducing potentials U → ∞ on each vacancy site, the t-matrix becomes −(P C G 0 P C ) −1 . Consider first a pair of vacancies on two sites 1 and 2 belonging to sublattices A and B. P C G 0 P C becomes a 2 × 2 matrix whose diagonal matrix elements are equal to G 0 00 . In the limit E → 0, the off-diagonal elements G 0AB 12 are real and independent of the energy, so that the resonance at E = 0 of a single vacancy is replaced by a pair of bonding-antibonding states at energies given by: Notice here that the equation x ln |x| = y has the asymtotic solution x = y/ ln |y| + . . . . Actually for first neighbours, G 0 12 is not really small, the imaginary parts can no longer be neglected and the low energy limit is no longer satisfied. Resonances close to zero energy will only occur for large separations of the vacancies. Consider now a general cluster containing N v vacancies of each type. The condition for resonances becomes det(P C G 0 P C ) = 0. Close to zero energy the only nonvanishing elements of P C G 0 P C are the real part of the G 0AB Green functions. We have therefore to study the eigenvalues of P C G 0 P C considered (at zero energy) as an effective hamiltonian between vacancies on different sublattices, with long range 1/n interactions. 48 This is obviously a difficult problem in the case of a finite concentration of vacancies similar to the impurity band problem. 39 To summarize, we have clarified some points related to the electronic structure of a finite number of vacancies in graphene. The case of a finite vacancy concentration will be dealt with elsewhere.

ACKNOWLEDGMENTS
Numerous discussions with H. Amara are gratefully acknowledged.