{'fraction_non_alphanumeric': 0.03511670041564532, 'fraction_numerical': 0.018970478524992007, 'mean_word_length': 4.307268099547511, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Long-range interactions between substitutional nitrogen dopants in graphene: electronic properties calculations

Being a true two-dimensional crystal, graphene has special properties. In particular, a point-like defect in graphene may have effects in the long range. This peculiarity questions the validity of using a supercell geometry in an attempt to explore the properties of an isolated defect. Still, this approach is often used in ab-initio electronic structure calculations, for instance. How does this approach converge with the size of the supercell is generally not tackled for the obvious reason of keeping the computational load to an affordable level. The present paper addresses the problem of substitutional nitrogen doping of graphene. DFT calculations have been performed for 9x9 and 10x10 supercells. Although these calculations correspond to N concentrations that differ by about 10%, the local densities of states on and around the defects are found to depend significantly on the supercell size. Fitting the DFT results by a tight-binding Hamiltonian makes it possible to explore the effects of a random distribution of the substitutional N atoms, in the case of finite concentrations, and to approach the case of an isolated impurity when the concentration vanishes. The tight-binding Hamiltonian is used to calculate the STM image of graphene around an isolated N atom. STM images are also calculated for graphene doped with 0.5 % concentration of nitrogen. The results are discussed in the light of recent experimental data and the conclusions of the calculations are extended to other point defects in graphene.


I. INTRODUCTION
Local defects and chemical doping is a welldocumented way to tune the electronic properties of graphene. 1 More specifically, the benefits of doping have been underlined in the context of (bio)sensing, 2 lithium incorporation battery, 3,4 and in other fields. Nitrogen (N) is a natural substitute for carbon in the honeycomb structure due to both its ability to form sp 2 bonds and its pentavalent character. It is not a surprise, then, that many publications deal with the production of N-doped graphene realized either by direct growth of modified layers 5,6 or by post-synthesis treatments. 7,8 Doping singlewall carbon nanotubes has also been considered. [9][10][11] Recent STM-STS experiments of N-doped graphene 5,8 have demonstrated the occurrence of chemical substitution. These experiments provide us with a detailed and local analysis with sub-atomic resolution of the electronic properties of the doped material. The STM images display a pattern having three-fold symmetry around the N atoms and having a strong signal on the C atoms bonded to the dopant, which ab-initio simulations reproduce well. 12 STS measurements have revealed a broad resonant electronic state around the dopant position and located at an energy of 0.5 eV above the Fermi level. 8 It is important to understand the effects of local defects and chemical doping on the global electronic properties of graphene, on its quantum transport properties and on its local chemical reactivity. It is the reason why the electronic properties of graphene doped with nitrogen have been investigated by several groups, mainly on the basis of ab-initio DFT techniques. 5,[12][13][14][15][16][17] The central advantage of the ab-initio approach is to be parameter free. A disadvantage is to be restricted to periodic systems as long as fast Fourier transforms need to be used to link direct and reciprocal spaces. In most instances, doping is therefore addressed in a supercell geometry that precludes the study of low doping concentration or disorder. In the case of single-wall carbon nanotubes, however, the electronic properties of isolated defects have been studied by first-principle methods based on scattering theory. 18 For a nitrogen impurity in the (10,10) armchair nanotube, the local density of states on the N atom shows a broad peak centered at 0.53 eV above the Fermi level. The charge density associated with that quasi-bound state (donor level in the language of semiconductor physics) extends up to ∼1 nm away from the defect. This means that N atoms located ∼2 nm apart can interact, which requires examining with care the intrinsic validity of a supercell method.
The present work, based on both ab-initio DFT and semi-empirical tight-binding electronic structure calculations, aims at looking for interference effects generated by a distribution of N dopants in graphene as compared to the case of an individual impurity. We resort to different tight-binding parametrization strategies. The simplest one, based on just one adjustable parameter related to the defect, permits analytical calculations of the impurity levels. This model is described in Appendix A. We favor a more realistic approach, in which the perturbation induced by the defect is allowed to extend on the neighboring sites. This latter model is used to study the effects of disorder on the local and global electronic properties of doped graphene and to calculate STM im-ages that we compare with available experimental data. The calculations and discussions are developed here for substitutional nitrogen. The conclusions would be qualitatively the same for boron doping and for other types of point-like defects.

II. METHODOLOGY
The SIESTA package 19 was used for the ab-initio DFT calculations. The description of the valence electrons was based on localized pseudo-atomic orbitals with a double-ζ polarization (DZP). 20 Exchange-correlation effects were handled within local density approximation (LDA) as proposed by Perdew and Zunger. 21 Core electrons were replaced by nonlocal norm-conserving pseudopotentials. 22 Following previous studies, 8,12 the first Brillouin zone (BZ) was sampled with a 15 × 15 grid generated according to the Monkhorst-Pack scheme 23 in order to ensure a good convergence of the self-consistent electronic density calculations. Real-space integration was performed on a regular grid corresponding to a planewave cutoff around 300 Ry. All the atomic structures of self-supported doped graphene have been relaxed.
As mentioned above, the DFT calculations are based on a code that requires a periodic system. As a consequence, a supercell scheme was adopted to handle substitutionaly doped graphene. The atomic concentration of dopants is then directly related to the size of the supercell: 0.6 at% for a 9 × 9 supercell and 0.5 at% for a 10 × 10 supercell.
For the tight-binding parametrization, the following, well-established procedure was followed. [13][14][15] The Hamiltonian included the π electrons only, the hopping parameter between nearest-neighbor atoms was set to the ab-initio DFT value γ 0 = -2.72 eV. The same hopping was used between the N atom and its C first neighbors, which is not a severe approximation. It is indeed shown in Appendix A that a non-diagonal perturbation of the Hamiltonian can be offset by a renormalization of the N on-site energy, and this energy is a free parameter to be adjusted to the results of the ab-initio calculations. The C on-site energies were assumed to vary with the distance d to the impurity according to a Gaussian law: where |U | is the depth of the potential well induced by the nitrogen (boron induces a potential hump, instead 13,24 ) and ε C is the asymptotic bulk value of the on-site parameter of carbon. ε C defines the Dirac energy of pure graphene and has no other influence on the density of states (DOS) than a systematic shift of all the electronic levels. The standard deviation σ in eq. (1) was set to 0.15 nm, according to data of Ref. [24], whereas U was chosen so as to best fit DFT local DOS calculated on the N atom in the 9 × 9 supercell (see below). In order to avoid the constraints of a periodic tightbinding Hamiltonian, the local DOS have been computed with the recursion method. 25 150 pairs of recursion coefficients were calculated and extrapolated to their asymptotic values related to the edges ε C ± 3γ 0 of the π bands of pure graphene. It is worth emphasizing that this procedure avoids to introduce any broadening of the energy levels. A drawback is the presence of wiggles that may be observed in some densities of states. They are sorts of Gibbs oscillations generated by the Van Hove singularities, in particular by the abrupt discontinuities of the graphene π DOS at both band edges. Fig. 1(a) shows DFT local densities of states on the N atom, on the first-neighbor atoms (C1) and on carbons located farther (0.5 -0.7 nm) away (bulk), in a 9 × 9 supercell. A broadening of 0.05 eV on a 45 × 45 grid was used for the calculations. The important observation is the existence of two peaks, located at 0.55 and 0.92 eV above the minimum of the density of states in the N local DOS. 8 The present DFT calculations are in overall good agreement with results published for 4 × 4, 13,16 5 × 5, 17 7 × 7, 5 and 10 × 10 supercells. 12 The double-peak structure in the unoccupied part of the local N DOS plotted in Fig. 1(a)) was not observed in previous calculations due to the large energy broadening that was used there, 5,12 except in Ref. [8]. It will be demonstrated below that the observed double-peak structure is a consequence of longrange interactions between the N impurities that were reproduced periodically in the graphene sheet because of the supercell geometry. Fig. 1(b) shows tight-binding densities of states of the same 9 × 9 supercell after Lorentzian broadening aimed at facilitating the comparison with the DFT results. The latter are qualitatively well reproduced by the π-electron tight-binding Hamiltonian by taking |U | = 4 eV in eq. (1). The present work aims at tailoring a simple tool for exploring the effects of the defect configurations. For that reason, achieving the best possible fit of DFT calculations for a specific supercell is superfluous. By comparison, the on-site energies in Ref. [14] can be approximated by a sum of two Gaussians, ε(d) = ε C − 2.95 exp(−0.5d 2 /σ 2 1 ) − 0.59 exp(−0.5d 2 /σ 2 2 ) (in eV), with σ 1 = 0.10 nm and σ 2 = 0.39 nm. The second term has a longer range than the single Gaussian law used in the present calculations, and the potential well at the N location is about 10% shallower than here.

III. PERIODIC DISTRIBUTION OF N
The tight-binding densities of states of the same 9 × 9 superstructure, now calculated without energy broadening, are displayed in Fig. 1(d). A substructure clearly emerges where the broadened local DOS only showed the double-peak feature discussed here above. In addition, Fig. 1(d) reveals that there are much less electronic states in the region between 0 and 1.5 eV on the secondneighbor atoms (C2) than on the first-(C1) and thirdneighbor (C3) carbons (see  The same observation as for C2 can be made for the local DOS on the fourth neighbors (not shown). This alternation between lower and higher densities of states brings out the fact that the two sublattices of graphene are differently affected by a substitutional defect (see Appendix A). The reliability of the tight-binding calculations compared to DFT can be appreciated from Figs. 1(e,f) obtained for a 10 × 10 supercell (0.5 at% concentration of N). The good agreement between the two approaches reinforces the fact that the parametrization of the tightbinding on-site energies used for the 9 × 9 supercell describes well the doped graphene systems and does not require further adjustment depending on the actual concentration. Electronic states localized on and near the N dopant are observed within 1.5 eV from the Dirac energy, as for the 9 × 9 supercell. However, the multi-peak structure in the local N DOS of 10 × 10 ( Fig. 1(f)) is totally different from that of 9 × 9 ( Fig. 1(d)), despite similar N concentrations (0.5 at% and 0.6 at%, respectively). This finding is a clear indication of the importance of interference effects among the supercells.
The band structures of the DFT calculations shown in Fig. 2 provide us with another difference between the two superstructures: the 10 × 10 N-doped superstructure has a direct band gap of ∼0.05 eV at the K point of its Brillouin zone, whereas the 9 × 9 N-doped superstructure has no gap. Symmetry considerations and arguments from perturbation theory developed in Appendix B explain why it is so. The K and K ′ points of graphene move to the K and K ′ points of the folded Brillouin zone of a p × p supercell when p is not divisible by 3, whereas they are both mapped onto the Γ point when p is an integer multiple of 3. 26 In the latter case, the fourfold degeneracy of the Dirac energy at the Γ point of a supercell of pure graphene is partly lifted by the perturbation brought about by the N atoms. As demonstrated in Appendix B, there remain π and π * bands that cross the Dirac energy at Γ. This crossing as clearly visible in Fig. 2(a) for the 9×9 superstructure, which therefore has no gap. In direct space, a doped superstructure has point group symmetry D 3h . The symmetry of the Γ point of the superstructure, D 3h , which is the same as that of the K and K ′ points of pure graphene, 27,28 has irreducible representations of dimension one and two. A twofold degeneracy of some electron energy bands is therefore allowed at the zone center of a doped superstructure. By contrast, this degeneracy is forbidden at the K and K ′ points of the folded zone because the symmetry of theses points, from D 3h in pure graphene, is reduced to C 3h in the doped superstructure. The latter point group has only one-dimensional irreducible representations. As a consequence, the crossing of π and π * bands is symmetry forbidden at K and K ′ . A periodic substitution of N for C therefore opens a gap in these p × p supercells for which p is not divisible by 3, as shown in Fig. 2(b) for the 10 × 10 superstructure. It is demonstrated in Appendix B that the band gap in those superstructures scales with p like E g ≈ V /p 2 , where V is of the order of 10 eV. Other local gaps of the order of V /p 2 appear in the band structures of Figs. 2(a,b) at the edges of the folded Brillouin zone. They produce pseudo-gaps and sub-peaks in the densities of states clearly discernible in Fig. 1 for both the 9 × 9 and 10 × 10 superstructures. As a result, supercells of the order of 30 × 30 at least would be necessary to reproduce the characteristic features of an isolated impurity with an energy resolution better than 0.01 eV.
In the tight-binding calculations, the band gap of the 10 × 10 superstructure is located 0.14 eV below the Dirac point of pure graphene (ε C ). The Fermi level E F of the doped superstructures was found to lie 0.27 eV and 0.25 eV above ε C for the 9 × 9 and 10 × 10 systems, respectively, which means 0.43 eV and 0.39 eV above the minimum of the DOS (the corresponding DFT values are 0.42 and 0.36 eV, respectively). The σ bands below the Fermi energy are completely filled, the π bands of the supercell must accommodate an extra electron brought about by the nitrogen. The occupied local π DOS of the N atom contains 1.36 electron in both the 9 × 9 and 10 × 10 supercells. The remaining 0.64 excess electron is distributed on the surrounding carbons, of which a total of 0.56 sits on the three first neighbors (C1), which are therefore negatively charged.
For both 9 × 9 and 10 × 10 doped supercells, there is a localized state expelled from the π band and located 9.35 eV below E F , to be compared with 8.65 eV in DFT calculations. This localized state weights more than 25% (0.57 state of the N local DOS, which can accommodate 2 electrons including the spin degeneracy). The large value of the weight can be understood from the arguments developed in Appendix A for a simplified tightbinding Hamiltonian.
We conclude here that the supercell technique would require very large supercells (up to 30 × 30) if it were desired to describe the electronic states of an isolated defect. In other words, the supercell technique is badly convergent when the size of the cell increases, particularly in low dimensions. The corresponding long range of a local perturbation produces interferences between the images of the defect generated by the periodic boundary conditions (see also Ref. [29]).

IV. ISOLATED N IMPURITY
We now turn our attention to the case of an isolated N impurity studied by tight-binding. The local DOS on and around an isolated N impurity are presented in Fig. 3. The calculations were performed on a 150x150 cell of graphene (45,000 atoms) with periodic boundary conditions. By restricting the number of pairs of recursion coefficients to 150, the impurities located in adjacent cells do not feel each other, which means that the N atoms are virtually isolated. In the N local DOS, there is a tall asymmetric peak located 0.5 eV above the Dirac energy (ε C ) of graphene which coincides here with the Fermi level. The peak broadens and shifts to 1 eV on the first neighbors and moves up further to 1.5 eV on the second neighbors. The local DOS on the third-neighbor atoms reproduces the resonance peak of the impurity, with a smaller amplitude. The local DOS on the fourth neighbors (not shown) bears resemblance with that of pure graphene. It is interesting to observe in Fig. 3 that there remains almost nothing of the Van Hove singularities of graphene at ε C ± γ 0 on the N and C1 sites.
The problem of an isolated impurity in graphene can be addressed analytically if one simplifies the Hamiltonian to a point-like defect by ignoring the Gaussian delocalization of the perturbation (σ → 0 in eq. (1)). Details are presented in Appendix A. If this simple model captures the essential physics of the problem, it fails in producing all the details set up by the delocalized nature of the perturbation. As shown in the Appendix A, perturbing the first-neighbor C atoms, in addition of course to the substitution site, leads to a better picture.

V. RANDOMLY DISTRIBUTED N ATOMS
Experimentally, nitrogen doping of graphene leads to randomly distributed substitutional sites, 5,8 with perhaps some preference of the dopant atoms to sit on the same sublattice, at least locally. 5 The calculations illustrated in the present Section were performed for nitrogen distributed randomly with an atomic concentration of 0.5%, identical to that realized with the 10 × 10 superstructure. The selection of the substitutional sites was constrained by the requirement that the distance between two nitrogens remains larger than 12σ = 1.8 nm (see eq. (1)). The reason for that was to avoid any overlapping of the potential wells generated by the dopants. Fig. 4 shows a configurational averaging of the local densities of states when the N atoms have equal probabilities to sit on one or the other of the graphene sublattices. This case will be referred to as "unpolarized" (or "uncompensated"). 30 By configurational averaging, it is meant the arithmetic average of local DOS on 50 nitrogens selected randomly -among the 225 dopant atoms that the 150x150 sample box contains -and the arithmetic average of local DOS on 150 randomly-selected C1, C2 and C3 sites. The N, C1, C2, and C3 average local DOS for the unpolarized disordered distribution all have a remarkable similarity with the ones obtained for an isolated impurity (Fig. 3). The result are in agreement with N local DOS for randomly distributed impurities obtained by Lherbier et al. 15 with a small broadening. When the N atoms are put, still randomly, on the same sublattice ("polarized" case"), the multi-peak structure characteristics of the superlattice (see Fig. 1(e,f)) reappear on the average local DOS shown in Fig. 5. Interestingly, the local DOS on the second-neighbor atoms (C2), which belongs to the same sublattice as the N atom they surround, are virtually the same in Figs. 4 and 5, insensitive to the polarization of the N distribution. It is also interesting to remark that, like with the 10 × 10 superstructure, there is a tiny gap of states 0.1 eV below the Dirac energy. The gap looks narrower in case of the unpolarized distribution compared to the polarized one. In the polarized case, the appearance of a gap is natural since the average diagonal potential breaks the symmetry between the two sublattices. Similar effects have been observed in the case of vacancies. 30 Averaging the local densities of states makes visible the distinction between unpolarized and polarized distributions of N. However, for a given configuration of the dopant atoms, the local DOS varies from site to site. How much that variation can be is illustrated in Fig. 6, which compares local DOS calculated on three different N atoms for both the unpolarized and the polarized distributions. The variations from site to site are so important that the local DOS is not a reliable indicator of the global distribution of the nitrogen atoms among the two sublattices. If it is true that the N local DOS displayed in Fig. 6 are less peaked that the N local DOS of the 10 × 10 superstructure (Fig 1(e,f)), identifying which is which would be difficult on the experimental side if one had only STS spectra to deal with.
Finally, it may be important to remark that the minimum of the π density of states around the N dopants does not coincide with the atomic level ǫ C of carbon in graphene (the Dirac point energy; see also Refs. [30,31,32]). Depending on the concentration, the DOS minimum shifts slightly below ǫ C (see e.g. Figs. 4 and 5).

VI. STM IMAGES
A reliable information on the electronic structure of doped graphene can be obtained by STM imaging. A tip polarized negatively compared to the graphene layer probes the unoccupied states where a N substitutional impurity and the adjacent carbon atoms have peaks in their local DOS. The increase of electronic density of states, compared to graphene, in this energy region gives rise to a bright triangular spot in the STM image, 5,8,12,17 with two possible orientations with respect to the honeycomb lattice depending on the sublattice on which the N dopant sits. These two orientations are rotated by 180 • from each other, as observed experimentally. 5,8 Fig. 7 is a tight-binding STM image 33 computed for graphene with a single N impurity. The prominent triangular arrow-head at the center of the image is located on the N, C1 and C3 sites. The second-neighbor carbons (C2) atoms do not participate much to the STM signature of the dopant. This is so because the local DOS in the energy window probed by the STM current (0 to 0.5 eV) is smaller on the C2 atoms than on the N, C1 and C3 atoms (see Fig. 3).

FIG. 7.
Tight-binding STM constant-current image of graphene with an isolated N impurity located near the center, computed for a negatively-polarized tip (Vtip = -0.5 V). The vertical scale is the tip height (Å) above sample. The image size is 2 nm along both sides. Fig. 8 shows two STM images computed for graphene doped with N at 0.5 at% concentration, the dopants being randomly distributed on the two sublattices. In the configuration shown in image (a), there are two impurities located 1.9 nm apart sitting on the same sublattice. The two nitrogens captured in images (b) do not sit on the same sublattice, which explains why the related triangular STM patterns are oriented differently. The scale used for the STM signal (tip height at constant current) is identical for both images. For improving the contrast, the scale has been saturated below 4.0Å and above 5.8Å. Two dopants interfere more when they are located on the two sublattices than when they sit on the same sublattice.
In tight-binding STM theory, 33 there is an empirical parameter coupling the tip apex to the sample atoms, which was chosen independent of the chemical nature of the probed atoms. In other words, the STM calculations differentiate a N atom from the C atoms only through intrinsic effects that the former has on the electronic structure of the doped graphene. In DFT calculations, the standard way to generate an STM image is via Tersoff-Hamann's theory. 34 In addition to local variations of the DOS, the tunnel current depends on the spatial extension of the π orbitals in the direction perpendicular to the atomic plane. For doped graphene, the 2p z orbital of N decreases more rapidly than the C 2p z orbital does, 12 which means that the contrast of the computed image depends on the tip-sample distance. This dependence is illustrated in Figs. 9(a) and (b) that represent DFT constant-current STM images computed for a small current (large tip-sample distance) and for a large current (small tip-sample distance), respectively. The local DOS is larger on top of the N atom than on the top of all the other atoms, but it decreases much more rapidly in the normal direction. This explains the small brightness of the N site for the large tip-sample distance. Experimentally, variations of the STM contrast have also been observed as a function of the current/voltage conditions, 8 where image of the dopant atom protrudes in some occasions and does not in other occasions. It is tempting to attribute this observation to variations of the tip-sample distances depending on the current setpoint, as in Fig. 9. However, the experimental distance is much larger than the one that can be achieved numerically in DFT, due to the fast decay of the localized basis set when moving away from the sample surface.

VII. CONCLUSIONS
Tight-binding and DFT electronic calculations reveal that the perturbation induced by N dopants in graphene is long ranged. This is not a surprise when one remembers that the Green function of the free-electron problem in two dimensions is the Bessel-Hankel function of zeroth order (i/4)H k = 2 /2mE. For the π-electron tight-binding Hamiltonian on the honeycomb lattice, the Green function elements G 0 nm can also be approximated, at low energy, by Hankel functions of δk| r m − r n |, now of order zero or one, depending on whether they refer to sites m and n located on the same sublattice or not. 35 In this regime, the energy E is linearly related to δk by the dispersion relation E = v F δk, where v F is the Fermi velocity, and the Hankel functions are also multiplied by the energy and by periodic functions describing √ 3 × √ 3 modulations related to the Dirac wave vectors K and K ′ . For large separation distances, the Green function elements decay slowly, like 1/ E | r m − r n |/ v F , except for very small excitation energies where the Hankel functions diverge. A more detailed discussion will be given elsewhere. 36 The scattering formalism developed in Appendix A for the case of a point-like defect emphasizes the central role of the Green function in the understanding of the perturbation induced by the defect. Hence the long-range interaction between dopants in graphene. It is important to realize that this long-range interaction is not the consequence of using an extended perturbation of the onsite energies around the defect (eq. (1)). Things actually go the other way round: a defect perturbs the crystal potential far away, because of the slow decay of the scattering it produces. As a consequence, the parameters of the tight-binding Hamiltonian that mimic ab-initio calculations are affected in a sizable neighborhood of the defect. 14 This long-range interaction has several effects. First, it produces interferences between the duplicates of the defect generated by periodic boundary conditions when a supercell approach is used. The latter must therefore be used with caution. Second, substitutional impurities dispersed in graphene cannot be treated as independent defects as soon as the distance between two of them becomes too small, about 2 nm in the case of nitrogen. Increasing the strength of the local potential also increases the range of the perturbation it produces. This unusual behavior is actually related to the properties of the Green functions at low energy (see Appendix A and Refs. [31,37,38]). As mentioned in Appendix A, the deepest defect is the vacancy. Nitrogen, with its 10 eV effective perturbation parameter, is not a shallow defect. Boron would not be a shallow impurity either. B substitution can be described by a potential hump, instead of a well, with a similar |U | as nitrogen and a similar long-range perturbation. 13,24 In the case of more complex defects, such as N plus vacancy, 39 P plus N, 40 or O 3 , 41 it would be more difficult to define local parameters and to assess the importance of long-range interactions along the same way as in this paper.
As demonstrated experimentally, the partition of the honeycomb network in two sublattices has subtle effects on STM image of graphene with substitutional impurities, as exemplified by Fig. 8. The STM image of a dopant in graphene may be influenced by the proximity of another dopant. This is a direct consequence of the defect interactions mediated by graphene. What happens when two impurities come very close to each other remains to be clarified. A new parametrization of the tight-binding Hamiltonian would be required if one had to obtain, from calculations, the STM topography of neighboring N dopants.
Acknowledgments This research has used resources of the Interuniversity Scientific Computing Facility located at the University of Namur, Belgium, which is supported by the F.R.S.-FNRS under convention No. 2.4617.07 and the FUNDP.

APPENDIX A. THE IMPURITY PROBLEM IN GRAPHENE
In this Appendix, a simple but illustrative perturbation model is developed. It is assumed that the on-site energies all take the unperturbed bulk value ε C = 0, except on the nitrogen atom where ε N = U . The tight-binding Hamiltonian can be set in the form: The states |n denote π orbitals centered on sites n: ϕ π (r − n) = r|n . H 0 is the usual tight-binding Hamiltonian where only hopping integrals t nm ≡ γ 0 between first neighbors n and m are kept, and V is the localized potential of the N atom at site 0. The local density of states n(r, E) is given by: where n nm (E) is obtained from the Green function or resolvent G(z) = (z − H) −1 : G(z) can then be calculated in terms of the unperturbed Green function G 0 (z) = (z − H 0 ) −1 ; this is the so-called Koster-Slater-Lifshitz problem. If the matrix elements n|G(E + iε)|m are written G nm (z), we have: In particular, on the impurity site, we have G 00 = G 0 00 /(1 − U G 0 00 ), 42,43 and the local density of states on the N site, n N (E) ≡ n 00 (E) is given by: where F 0 (E) is the real part of the Green function G 0 00 (E + iε), i.e. the Hilbert transform of the unperturbed density of states of graphene n 0 (E).
Close to the Dirac point, n 0 (E) ∼ |E| and F 0 (E) ∼ E ln |E|, and provided U is large enough, n N (E) shows a resonant behavior around the energy E r such that (1−U F 0 (E r )) = 0 (see the open-circle symbol in Fig. 10). This has been discussed in many places. 30,32,35,37,[44][45][46] It turns out that a value for U of the order of −10 eV reproduces a resonance at the position obtained with the previous model (compare Fig. 11(a) with Fig. 3) instead of the value U = −4 eV in the previous approach. However, the shapes of the local DOS differ to some extent between the local-perturbation model developed in this Appendix and the more delocalized potential well considered throughout this paper. The resonance here (Fig. 11(a)) is weaker on the N site than on the sites C1 and C3 located on the other sublattice, which can be explained from the symmetry properties of the unperturbed Green functions: G nm (z) is and odd or even function of z depending on whether the sites n and m belonging to the same sublattice of the graphene structure or not. As a consequence, for example, G 0 00 (z = 0) = 0 whereas the real part of G 0 0n (z = 0) does not vanish when n is a first neighbor, which has a big impact on the local densities of states (see eq. (5)). 36 Also the Van Hove singularities at 2.7 eV are still present here. In addition to the resonant state located above the Dirac point, Fig. 10 reveals the existence of a localized state lying below the π-electron DOS, where the onsite perturbation level 1/U intersects the curve F 0 (E) at −4.5|γ 0 | (solid square symbol). The weight (residue) of this localized state is inversely proportional to the slope of F 0 (E) at the intersection. This is the localized state found at 9.35 eV below E F in tight-binding and 8.65 eV below E F in DFT calculations for the doped superstructures (Section III).
This one-parameter model can be generalized to other type of local defects or chemical doping. B doping, with a positive U , will yield results symmetrical to those obtained for N: it can indeed be anticipated from Fig. 10 that the resonance state will now appear near the top of the occupied states, as confirmed by DFT calculations. 12,13,24 The limit U → ∞ corresponds to the introduction of a vacancy instead of a substitutional impurity. Because of the vanishing of the density of states and of the logarithmic divergence of F 0 (E), this limit is singular and must been studied with special care. 30,32,[45][46][47] The resonant state becomes very narrow, but does not become a genuine bound state. The total density of states varies as 1/([ln |E|] 2 |E|) and the perturbation of the electronic density is concentrated on the first neighbors and on the sites belonging to the sublattice different from that of the vacancy. 36 The simple model developed here above allows one to address the question of the role of N-C hopping interactions on the local DOS on the impurity. When the N-C hopping takes a value γ different from the one γ 0 between the C atoms, the Green function element on the perturbed site can be set in the form where Σ 0 (z) is the self energy of the unperturbed graphene. In the notations of the above formalism, the self energy can be identified to Σ 0 (z) = z − 1/G 0 00 (z). Inserting this expression in the right-hand side of eq. (7) leads to a generalization of eq. (6) valid for γ = γ 0 . In particular, the resonance condition in the vicinity of E = 0 becomes By comparison with the situation depicted in Fig. 10, the intersection of F 0 (E) needs now to be search for with a curve whose ordinate and slope at the origin are −γ 2 /(γ 2 0 |U |) and γ 2 (γ 2 0 − γ 2 )/(γ 2 0 U ) 2 , respectively. If |γ| < |γ 0 |, which should be the case here since the N atom is smaller than the C one, the ordinate at the origin of the curve moves up and its slopes becomes positive. These two effects pull the resonant energy E r closer to the origin compared with the simplest situation where |γ| = |γ 0 |. However, if the hopping perturbation is small, renormalizing the on-site level U to a larger absolute value U eff would produce the same effect. This conclusion validates the approach used so far to not modifying the N-C hopping.
We have finally considered an intermediate model between the Gaussian distribution of the on-site energies and the one-parameter model just described. We took U 0 on the N site and a second perturbation potential U 1 on the first neighbors, with values fixed by eq. (1), U 0 = −4 eV; U 1 = −2.57 eV. The agreement with the full Gaussian model is much better than with the one parameter model (Fig. 11(b)). It is fairly remarkable to realize that the intensity of the resonance state on the N site increases when the potential is delocalized on the first neighbors.

APPENDIX B. BAND GAP OPENING IN DOPED SUPERSTRUCTURES
The following perturbation theory is based on the same model Hamiltonian as in eq. (2), except that the potential V refers now to the perturbation brought about by the periodic array of nitrogen atoms in substitution for C in where K.n in the exponential represents the dot product of the wave vector of the K point of graphene and the position vector of the site n in real space.
It will be assumed that the sublattice 1 is an A sublattice. Then, the states |K B and |K ′B , having no amplitude on sublattice 1, are eigenstates of zero energy for any value of the on-site perturbation U . Very close to zero energy, furthermore, only the ((|K A , |K ′A )) subspace needs to be considered. In particular, where the Kronecker symbol imposes that K ′ − K be a vector of the reciprocal lattice of the supercell, i.e. it is one when p = 3q is a multiple integer of 3 and is 0 if p is not divisible by 3. As a consequence, in the case p = 3q, the four states with zero energy are mapped onto the zone center. Two of them, with non-zero amplitudes on the B sites, remain degenerate with zero energy; the other two states have energy E = 0 and E = −2|U |/p 2 to lowest order in perturbation. There are then two plus one states at zero energy, but this is an accidental degeneracy due to the simplified tight-binding model. When p is not divisible by 3, the degeneracy at the K and K ′ points of the folded zone is lifted: at each point, one state remains at zero energy (|K B or |K ′B , respectively) whereas the other state moves down at energy −|U |/p 2 to lowest order in perturbation.
All these features are clearly visible in the band structures of the 9×9 and 10×10 superlattices shown in Fig. 2. Since a realistic value for |U | is about 10 eV, the gap for the 10 × 10 case should be of the order of 10/p 2 = 0.1 eV, not far from the actual value (see Section III). For the 9 × 9 superstructure, the expected crossing between the two linear π and π * bands is observed at the Γ point close to two parabolic branches separated from each other by a gap about twice this value.