{'fraction_non_alphanumeric': 0.030187150201332082, 'fraction_numerical': 0.008863676653075695, 'mean_word_length': 4.522024891623549, '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}}
Friedel Oscillations in Graphene: Sublattice Asymmetry in Doping

Symmetry breaking perturbations in an electronically conducting medium are known to produce Friedel oscillations (FOs) in various physical quantities of an otherwise pristine material. Here we show in a mathematically transparent fashion that FOs in graphene have a strong sublattice asymmetry. As a result, the presence of impurities and/or defects may impact the distinct graphene sublattices very differently. Furthermore, such an asymmetry an be used to explain the recent observations that Nitrogen atoms and dimers are not randomly distributed in graphene but prefer to occupy one of its two distinct sublattices. We argue that this feature is not exclusive of Nitrogen and that it can be seen with other substitutional dopants.


I. INTRODUCTION
Nanoscale characterization techniques are fundamentally based on the changes experienced by otherwise pristine materials in the presence of symmetry-breaking impurities and defects. Such symmetry-breaking, for example in a Fermi gas, induces perturbations in the electronic environment of the gas through the scattering of its electrons 1 . These changes in the electronic scattering manifest as spatial oscillations, called Friedel oscillations (FOs), in quantities like the local density of states (ρ) and the carrier density (n), which radiate away from the location of the symmetry breaking perturbation and decay with the distance from the perturbation, D, with a rate linked directly to the dimensionality of the system and to some extent the resolution of the measurement. Much attention has been focused recently on such symmetry breaking in graphene 2-9 , a hexagonal lattice of sp-2 hybridised carbon with a wide range of unique properties 10 . Fig. 1 shows a schematic of the lattice with several different kinds of impurities, where atomic sites in the two triangular interconnected triangular sublattices composing graphene are represented by black and white symbols. In graphene, the vanishing of the density of states at the Dirac Point affects the decay rate of the change in carrier density (∆n) from D −2 , expected in a 2-D system, to a faster D −3 rate for ungated/undoped graphene and the oscillations disappear due to their commensurability with the lattice spacing 2,4,8 .
Previous studies examining the analytical behaviour of ∆n FOs in graphene have generally relied on a linearisation of the electronic bandstructure near the Dirac points and the introduction of a momentum cutoff [2][3][4]8 . In the current work, we present an alternative framework which removes these assumptions and matches numerical results exactly in the long-distance limit and over large energy ranges, paving the way for applications to other electronic quantities and graphene-like materials. The methodology we use is similar to that used to describe the Ruderman-Kittel-Kasuya-Yosida (RKKY) interaction in graphene 11,12 , a coupling effect between magnetic impurities which has been studied extensively and, like FOs, is another manifestation of symmetry breaking, making this work a natural extension of those techniques.
This methodology is applied to a range of commonly investigated impurity configurations, namely single and double substitutional impurities, vacancies, and realistic instances of the more commonly found top-and bridgeadsorbed atoms. The analytical expressions derived for the fluctuations in electron density for these impurities are corroborated with numerical calculations to confirm the predicted behaviour. We note how important features of the FO are dictated by the bipartite nature of the graphene lattice and furthermore, that different behaviours are observed for adsorbed impurities connecting symmetrically or asymmetrically to the two sublattices in graphene. The framework is extended to consider similar oscillations which occur in the formation energy of two impurities introduced into the graphene lattice in close proximity to each other. Such FOs in formation energy are consistent with recent experimental findings of sublattice-asymmetric doping of nitrogen substitutional impurities in graphene 13,14 . The ability to dope one sublattice of graphene preferentially opens many possibilities, including the opening of a bandgap 15 , and it is very interesting that such behaviour may be the manifestation of the strong sublattice dependence seen in FOs in graphene.
The paper is organised as follows. Section II introduces the relevant mathematical methods required, in particular the Green functions for graphene and the perturbations associated with different impurity configurations. Section III details the analytic approximations for the FOs in ∆n for the range of impurity type shown in Fig.  1 and compares their predictions to fully numerical calculations. As an application of our methodology Section IV investigates the appearance of sublattice-asymmetry in nitrogen-doped graphene and presents a simple tight binding model for the long-range sublattice ordering witnessed in recent experiments.

A. Green Functions
We begin by outlining the Green Functions (GFs) methods which play a central role in our approach as the Friedel oscillations in electron density and local density of states, ∆ρ and ∆n respectively, are directly obtainable using them. The retarded single-body GF, G ij , between two cells i and j in the pristine lattice is calculated through diagonalisation of the nearest-neighbour tightbinding Hamiltonian using Bloch's Theorem and can be expressed as an integral over the Brillouin Zone 11 where the energy, E, includes an infinitesimal positive imaginary part, t = −2.7eV is the pristine graphene lattice hopping integral, r = ma 1 + na 2 (where m, n ∈ Z) is the spatial separation of the unit cells i and j containing the relevant sites which is expressed in terms of the lattice vectors a 1 and a 2 shown in Fig. 1 and the variable f (k) = 1 + e −ik.a1 + e −ik.a2 . In our analytic work, we principally examine armchair direction separations, de-fined by m = n and work in distance units of m+n 2 . The matrix form of Eq. (1) captures the inter-sublattice nature of the GF calculation between the two sites in the graphene unit cell, such that ij is the pristine lattice GF from the s 1 sublattice site in cell i to the s 2 sublattice site in cell j. For conciseness we will omit the sublattice indices from hereon and use g ij to denote the pristine lattice GF between two sites on the same sublattice sites unless specified otherwise.
To aid numerical calculation one of the two integrals in Eq. (1) can be solved analytically, and in some cases it is possible to approximate the second integral using the stationary phase approximation (SPA) 11 . These methods take advantage of the highly oscillatory nature of the integrand for large r and to approximate the integral at values near the stationary points of the phase, ignoring the highly oscillating parts which mostly cancel each other. The advantage of the SPA is that it is applicable across the entire energy band, with no intrinsic momentum or energy cut-offs, although maximum achievable experimental gating (through ion gels) is limited to E F ∼ 0.5|t| currently 16 . Applying the SPA assumption to g ij gives us a sum of terms in the form where the coefficients A(E) are dependent on the sublattice configuration of sites i and j and the direction between them. For two sites on the same sublattice separated in the armchair direction we find valid for |E| < |t| where D = m+n 2 is the separation between the sites i and j and Q(E) is associated with the Fermi wavevector in the armchair direction. This approximation works best for separations beyond a few lattice spacings in the armchair and zig-zag directions as other directions are not always analytically solvable 11 . For an armchair separation of D = 5 the agreement between analytic and numerical calculations of g ij is excellent, with less than 1% deviation over 95% of the energy spectrum.
To calculate the fluctuations in properties of a system when an impurity is introduced, we will usually require the difference between the GFs describing the pristine (ĝ) and perturbed (Ĝ) systems. This can be expressed, using the Dyson equation 17 , as whereV describes the potential applied to the pristine system to introduce the impurity. By finding suitable i1 |τ | 2 g aa 1−g aa g 11 |τ | 2 (g i1 +g i2 ) 2 g aa τ 2 1−(g 11 +g 12 +g 21 +g 22 )g aa τ 2 TABLE I: Perturbation operators (V ) and the corresponding fluctuation in the system GF at a general lattice site i (∆Gii). Our notation follows the rule that a corresponds to an adsorbate atom site and the numerical indices specify the location of impurities -e.g. 0 for a substitutional impurity and 1 or 2 for the carbon host sites of adsorbates.
descriptions and parameterizations for different impurities in the graphene lattice we can find directly the change in the corresponding perturbed lattice GF ∆Ĝ, and also how the ρ and n are altered from the pristine system. Exact parameterization, however, is not that important when considering the qualitative features of the phenomena we investigate here. More precise parameterization can be found, for example, by comparison to DFT calculations 18 .

B. Impurities
In this section we present the types of impurities that will be considered and the respective perturbative potentials used to describe them. Substitutional impurities in graphene, shown schematically by A and B in Fig 1, occur when single carbon atoms are replaced with dopant atoms such as nitrogen.
The simplest way to model them is by introducing a per-turbationV Subs. (Table I) which alters the onsite energy of the site in question by a quantity λ. The onsite energy of the carbon sites neighbouring the impurity and their overlap integrals with the impurity site are presumed to be unchanged in this simple model. However they can be easily incorporated, as can additional orbitals beyond our single orbital approximation for the impurities, if a more accurate parameterization is required. The presence of a substitutional impurity at site 0 induces a fluctuation of the diagonal matrix element GF, for example ∆G ii at site i, which can be found through applyingV Subs. to Eq. (4) the result of which is shown in Table I. Vacancy defects, formed by the removal of a carbon atom from the graphene lattice ( Fig. 1 (C)), can be considered as a phantom substitutional atom in the limit where the onsite energy,λ → ∞. In a physical context this is effectively excluding the state to the electrons in the lattice. Experimentally, vacancies can be induced by ion bombardement 19 . An approximation of the corresponding ∆G Subs.
ii , using the limit λ → ∞, is shown in Table I. Alternatively this can be modelled by removing the hopping between the site and the nearest neighbours.
Adsorbates, which bind to atoms in the graphene lattice, are characterised by an atom with onsite energy a which is initially disconnected from the graphene sheet with an onsite GF g aa = 1/(E − a ). A top-adsorbed adatom (Fig 1 (D)) is connected to the lattice by the perturbationV T op in Table I Table I.

C. Charge Density Perturbations
The effect of a general perturbationV on the charge density on site i (∆n i ) is calculated through an integration of the change in LDOS, ∆ρ i , where f (E) is the Fermi Function and ∆ρ i = − 2 π Im∆G ii relates the change in local density of states to the perturbed diagonal GF ∆G ii . Physically, n i is the number of states below the Fermi energy (E F ) that are filled by electrons on site i, where ρ i describes the energy distribution of these states at the site. When calculating ∆n numerically the integral is evaluated along the imaginary energy axis to avoid discontinuties along the real axis. This transformation is done through forming a contour in the upper half plane, which contains no poles, and evaluating the contour integral via Cauchy's Theorem.

D. Changes in Total System Energy
Green Functions methods can also be used to quantify other phenomena associated with the lattice, for example the change in the total system energy due to the bonding of impurities follows from a sum rule and is given explicitly by the Lloyd Formula 20 This quantity is directly related toV and is useful for finding energetically favourable impurity positions and can be used to investigate the dispersion and clustering of impurities 21,22 . For multiple impurities the expression will contain interference terms which result in changes from the single impurity case and can reveal favoured configurations in the lattice.

A. Weak Substitutional Impurities
To begin, we consider the charge density variations ∆n at all lattice sites surrounding a substitutional impurity of strength λ = |t|, situated on a site in the black sublattice. A numerical evaluation using Eq. (5) and ∆G Subs.
ii from Table I  , where we see FOs in the charge density radiating away from the central impurity on both sublattices with a wavelength determined by E F . There is clear sublattice asymmetry in ∆n between the black (a) and white (b) sublattices with ∆n swapping signs between the sites in the same unit cell. This signature is important when considering multiple impurities, which we will discuss in sections III B and IV. It is possible to approximate ∆ρ and ∆n along the armchair direction (dashed line in Fig. 2 (a)), by applying the SPA approach and the Born Approximation, which is valid for weak scatterers of strength λ |t|, to ∆G Subs. ii , resulting in where g 2 i0 can then be expressed using Eq. (2) as and we can solve the integral via contour integration in the upper-half of the complex energy plane, the poles in the integrand being given by the Matsubara frequencies. Taking the limit of zero temperature gives the sum The sum coefficients γ (E) are related to the SPA coefficients and defined as where B = A 2 and X ( ) denotes the th derivative of the function X with respect to energy. B ( ) , Q (0) and Q (1) are evaluated at the Fermi energy E F . The first order term of ∆n thus decays as D −2 with an oscillation period determined by Q 0 , and thus k F (Fig. 2 (c)). At the Dirac point we find that Q 0 → π in the phase factor which causes the sign-changing oscillations to become commensurate with the lattice spacing and seemingly disappear for all terms. Additionally, the energy dependent term B 0 → 0 and so the leading term of the series γ 0 = 0. Taking the next leading term of the series in Eq. 9 ( = 1) gives ∆n i decaying as D −3 and a comparison between the SPA approximation and the numerical result is shown in Fig. 2 (d). We also note that ∆ρ decays as D −1 as shown in Fig. 2 (e) away from the Dirac Point, which can be inferred directly from Eq. (7). The sublattice asymmetry is quite clear in the cross-section of ∆ρ, which is shown for sites on the same (black) and opposite (red) sublattice as the impurity. The opposing sign of the oscillations is a signature of FOs in ∆ρ, ∆n, ∆E and also of the RKKY interaction 12 . These results agree with previous work 2,4 , but due to our approximation method requiring large D we are limited to the long-range behaviour in the ranges of 5+ unit cells and so miss short wavelength features which are present in the region immediately surrounding an impurity which have been investigated in more depth by Bacsi and Virosztek 2,3 . In addition to what is shown in Fig. 2 we found excellent agreement between numerical and analytic calculations for all energy and sublattice configurations. When considering two or more substitutional impu-rities we can extend the matrixV Subs. (Table I) to include additional sites by addition of extra perturbations at the corresponding locations, for example for two identical impurites at arbitrary sites A and B we havê V = λ|A A| + λ|B B|. Fig. 3 (a) shows a contour plot of ∆n on the black sublattice for two such impurities spaced by D = 80 in the armchair direction, where the SPA can be used to approximate ∆n along the dashed line. Generally, if λ t then we find through Eq. (5) and the Born Approximation that to first order ∆n has the form

B. Multiple Weak Substitutional Impurities
where ∆ A = g iA λg Ai , ∆ A and ∆ B arise from the effects of the individual isolated impurities and the extra term ∆ AB describes the interference effect. The approach for the single impurity can be adapted and applied to all three terms, where the integration of ∆ A and ∆ B will be identical to the single impurity case and the interference term ∆ AB can be approximated well if A and B are at least a few unit cells apart. Clearly ∆ AB is dependent on λ 2 , but also on the separation of A and B. Thus this term decays rapidly when the scatterers are weak and/or separated by several unit cells. Consider a cross-section of the charge density fluctuations indicated by the black dashed line in Fig. 3  a). Applying the SPA derived in the previous section we can match the interference pattern seen, as shown in Fig.  3 b), achieving an excellent match between analytic and numerical methods. By breaking apart the SPA terms we find that the dominant contribution is from the isolated impurities (given by ∆A (black curve) and ∆B (red curve) in Fig. 3 c)) with a very small contribution from the interference term ∆AB regardless of chosen energy (Fig. 3 d)).
Whilst the SPA works well when D is suitably large and λ small, the simple approximation breaks down when the impurities are moved to within a couple of unit cells as the contribution of the interference terms becomes more important, especially for impurities on opposite sublattices. In addition to the interference term ∆ AB , the region of significant overlap between the ∆ A and ∆ B increases and the FO patterns observed become more complex. In Fig. 4 we examine the numerical contour plot on the black sublattices of ∆n for two such configurations, namely two substitutional atoms of the type considered in Fig. 2  first case, we note the strongly sublattice dependent behaviour noted in Section III A is still present, whereas in the nearest-neighbour impurity case it has been mostly washed out due to a superposition of the features observed in panels a) and b) of Fig. 2. The strikingly different interference patterns present clear signatures for the two cases and this general qualitative difference in FOs may make impurity configurations easier to distinguish. The importance of this cross sublattice effect is apparent when considering FOs in other quantities, and will be discussed later in the context of energetically favourable doping configurations for multiple nitrogen substitutional impurities in graphene.

C. Vacancies and Strong Scatterers
Taking the limit λ → ∞ for a substitutional impurity, as shown in Table I, corresponds to placing a vacancy in the lattice and yields We note that the change in charge density on the impurity site, given by g i0 = g 00 in Eq. (12), becomes ∆n 0 = −1 which corresponds to a complete depletion of electrons on this site. The pristine onsite g ii can be approximated very well for energies in the linear regime as 23 and this approximation works at energies up to approximately E F ∼ |t|/2. Eq. (12) can be solved in a similar fashion to that of the single substitutional impurity by observing that g ii is a function of energy only and absorbing it into the usual γ term in Eq. (9) and Eq. (10), then following through with the usual derivation. There is no pole in the upper half plane for this integrand, so the evaluation method remains unchanged from the weak impurity case. Vacancies can be calculated numerically by either inducing a very large λ value on a site, or by disconnecting the site from its neighbours. A comparison of this numerical calculation and the SPA approximation of ∆n V acancy i (Table I) on sites within the same sublattice as the impurity is shown for a finite E F in Fig. 5 (a), and excellent agreement is seen between the two. The general features are similar to those noted for weaker impurities, namely sign-changing oscillations and a D −2 decay.
Interestingly, when we set E F = 0 we find a complete absence of FOs with ∆n = 0 at all sites, corresponding to no change in n from the pristine state. This behaviour is noted whether we examine sites on the same or on the opposite sublattice to the impurity, and is in stark contrast to the case of weaker scatterers where a non-oscillatory D −3 behaviour was noted. To try to understand this unusual behaviour, we examine the imaginary-axis integral that we need to solve to find ∆n, which is of the form We will see that it is the behaviour of the individual GFs along the imaginary axis which determines the presence or absence of the carrier density oscillations. At E F = 0 both g 00 and g i0 , if i is on the black sublattice, become entirely imaginary along the integration range. Similarly if i is on the white sublattice, g i0 is entirely real, ensuring that on both sublattices the integrand becomes entirely imaginary over the integration range and ∆n vanishes since it depends on the real part of the integrand only.
In the limit of large values of λ the perturbations described by Eq.(14) will vanish. In Fig. 5 (b), we examine the charge density fluctuations on the site nearest the impurity on both the same (black curve, A) and opposite (red, dashed curve, B) sublattice as λ is increased at E F = 0. These sites will have the largest change in carrier density due to their proximity to the impurity. We note that λ ≈ 4t (λ ≈ 10t) causes the largest amplitude in ∆n on the black (white) sublattice, with the amplitude of ∆n decreasing as λ increases further . It should be emphasised that this disappearance of ∆n in the λ → ∞ limit occurs only at E F = 0, whereas other energies will show the familiar D −2 oscillatory pattern.
It is straightforward to evaluate ∆ρ i from the ∆G V acancy ii approximation in Table I, where there is a clear asymmetry between the black and white sublattices at the Dirac Point, as g •• i0 → 0 (with • being the minority sublattice -a vacancy removes a carbon atom from the sublattice) and g •• i0 → +i∞ (• is the majority sublattice as it has more carbon atoms). This results in a zero density of states, similar to undoped graphene, on the minority sublattice and leads to divergencies in ∆ρ on the majority sublattice, which correspond to the widely predicted midgap resonance states seen in many previous works [24][25][26] . This is an explanation of the phenomenon of vanishing ∆n for the vacancy case as mentioned previously in this section, where the vacancy introduces a resonance peak at E = 0 in the LDOS of the majority sublattice sites with the bound state outside the band being removed to infinity. As the total number of states must be preserved such that −2 π Im ∞ −∞ dEρ = 2 at all sites the deformation of the LDOS by the vacancy to form the peak at E = 0 ensures that when E F = 0 the filled electron and hole states are symmetrical as in the pristine case leading to ∆n = 0 for these cases. On the minority sublattice the LDOS is symmetrical about E = 0 but without the resonance so again electron-hole symmetry is preserved. Once more the sublattice nature of graphene introduces significant asymmetries in the features of FOs surrounding a defect.

D. Adsorbates
An analytic expression for ∆n T op i , the FO induced by a top-adsorbed impurity (as shown in Fig 1 (d)), can be found using similar methods to that of a single substitutional atom, since ∆G T op ii in Table I is analytic in the upper half plane and there are no additional poles other than those at the Matsubara Frequencies. However, an approximation of g 11 as in Eq. (13) is required as this GF is beyond the scope of the SPA. This approximation restricts our results to the linear regime.
We consider two realistic cases of hydroxyl (OH − ) and hydrogen (H + ) adsorbates, using the parameters OH− = −2.9t , τ OH− = 2.3t for hydroxyl 27 and H+ = −t/16 , τ H+ = 2t for hydrogen 28 . The results for these impurities for E F = 0 and E F = 0.2|t| cases are pre- and H + (middle row) and to bridge-adsorbed carbon (bottom panel). The top two cases show both EF = 0 (left panels) and EF = 0.2|t| (right panels). Numerical data (black) is compared with the analytic expression (green), with corresponding log-log plots as insets. The red lines in the insets correspond to D −2 fits for both EF = 0.2|t| plots (panels (b) and (d)) and D −3 for the other cases. In panel c), we note that due to the parameterization of the H + , the SPA is a poor match to the numerics at EF = 0 away from the asymptotic limit. The charge density perturbations for a bridge-adsorbed carbon adatom, shown here as a numerical calculation (black) in the bottom panel, do not vanish at the Dirac Point, unlike a top adsorbed carbon, due to cross sublattice interferences. sented in Fig. 6. We note an excellent match between the numerical (black circles) and analytic (green lines) approaches for all cases with the exception of Hydrogen at E F = 0. The parameterization of hydrogen leads to divergencies at E F = 0 which requires us to look much further away to see agreement (approx. 100 unit cells). The proximity of the hydrogen onsite term to the Dirac point produces a resonance condition at E F = 0 and we see a clear change in the decay rate in the region close to the adsorbate in Fig. 6, with an eventual D −3 rate far from the defect. This has been studied recently by Mkhitaryan and Mishchenko 29 . In the case E F = 0 a decay rate of D −2 is found for both adsorbates, with D −3 at the Dirac Point, with proviso, matching the behaviour of a substitutional impurity which could be expected due to the possibility of modelling the effect of the adatom through a self-energy term replicating a substitutional atom 17 .
We note that altering the onsite energy to a = E F = 0 forces a resonance condition and confirms previous findings where the adatom behaviour can be similar to that of a vacancy 28 and the ∆n FOs disappear completely, which we may expect, for example, in the case of a top adsorbed carbon. However, carbon prefers a bridge adsorbed configuration 30 and due to this bonding type the interference effects from the two host sites, which are on opposite sublattices, lead to finite charge density perturbations at the Dirac Point. The presence of FOs in ∆n for a bridge adsorbed carbon at E F = 0 can be seen by considering an arbitrary site i and the corresponding ∆G Bridge ii (Table I).
In section III C we noted that, when considering imaginary axis integration of offdiagonal GFs at E F = 0, that g ab is either entirely real (for opposite sublattice propagators) or entirely imaginary (for same sublattice propagators). However from the form of ∆G Bridge ii it is clear that there will exist both real and imaginary terms in the integrand for ∆n i , which was not the case for top-adsorbed or vacancy impurities. This cross sublattice interference is key to the non-vanishing ∆n FOs for all bridge adsorbed atoms, regardless of parameterization, and in Fig. 6(e) we show a numerical plot of the familiar D −3 decay in ∆n at the Dirac Point. The usual D −2 oscillations are recovered with doping.

IV. SUBLATTICE ASYMMETRY IN NITROGEN DOPED GRAPHENE
Recent experimental works involving substitutional nitrogen dopants in graphene have reported a distinct sublattice asymmetry in their distribution, where the impurities are discovered to preferentially occupy one of the two sublattices, instead of being randomly distributed between them. This effect is noted at both long and short ranges 13 and is corroborated with DFT results 14 . A distinct and controllable sublattice asymmetry in doped graphene presents many interesting possibilities, among them the possibility of inducing a band gap by controlling the dopant concentration -an important step in the development of graphene-based field-effect transistors 15 . As remarked in the Introduction, Green Functions methods can be extended beyond FOs in ∆n and ∆ρ to include other quantities such as the change in total system energy (∆E), due to a perturbationV , using Eq. (6). The calculation of ∆E allows the investigation of favourable impurity configurations, and we will now apply this method to study substitutional nitrogen impurities in graphene within a simple tight-binding model where the nitrogen impurities are characterised 18 by λ N = −10eV.

A. Total Energy Change
We will first consider the interaction between two identical substitutional impurities at sites A and B with onsite energies λ, as in the multiple scattering case dis-cussed previously, so that the determinant in Eq. (6) becomes det(Î −ĝV ) = (1 − g AA λ)(1 − g BB λ) − g 2 AB λ 2 . It is possible then to identify two separate contributions to ∆E. The first one is associated with the individual impurities and is independent of the separation of A and B. Using g AA = g BB , we find where the superscript S refers to the substitutional impurities. The second contribution is an interaction term dependent on their separation through the off-diagonal GF g AB , (16) such that ∆E = 2∆E S 1 + ∆E S AB . It is easy to see that ∆E S 1 is both separation and configuration independent and takes a constant value.
To investigate the favourability of different configurations we define the dimensionless configuration energy function (CEF) for substitutional impurities This quantity describes the change in energy of the system due to the interference between the two impurities, relative to the total energy change in the system for two non-interacting (infinitely separated) impurities. Positive values of the CEF correspond to less favourable configurations whereas negative values correspond to favourable configurations which decrease the total energy of the system. By calculating the CEF for different impurity configurations we can establish which are energetically favourable and thus more likely to be realised in experiment. A map of β S (E F = 0) values for a large number of different configurations is shown in Fig.  7 (a). One impurity is fixed at the red circle corresponding to a site on the black sublattice. β S is then calculated with the second impurity located at each of the sites on the map, with the shading of the triangle surrounding each site corresponding to the β S value for that configuration. We note that, with the exception of nearestneighbour site impurities, a general trend is seen where the second impurity prefers to locate on the same sublattice as the initial impurity. This trend gives rise to the chequerboard-like pattern seen in Fig. 7 (a) where same sublattice (black) sites are surrounded by darker triangles, corresponding to lower energy configurations, than the opposite sublattice (white) sites. It is also clearly visible in Fig. 7 (b) where we plot β S (E F = 0) for the two sublattices separately for armchair direction separations. We also note here that the magnitude of β S decays as D −3 . This is the same rate noted for FOs in ∆n for substitutional impurities and is easily explained by examining g AB and hence distance dependence of Eq. (16) to first order in λ, which results in a similar equation to Eq. (8) for FOs in ∆n. Indeed, the same plot for E F = 0 in Fig. 7 (c) reveals an oscillatory behaviour and D −2 decay rate, again matching the ∆n behaviour. Thus the FOs both in ∆n for a single substitutional impurity and in β S for a pair of impurities display the same distance dependent behaviour due to the similar dependence on off-diagonal Green functions that appears in both quantities. An important point to note is the discrepancy for the nearest neighbour impurity cases which Fig. 7 (a) shows to be most favourable. In this case the asymptotic behaviour extracted from our analytic expressions is evidently not yet valid. However we should also expect that our parameterization approach is not adequate to describe such impurities, as it treats the nitrogen impurities separately and neglects, for example, the additional overlap matrix elements required for two neighbouring nitrogen impurities. Another point to note is that the introduction of nitrogen impurities into graphene increases the total number of electrons in the system and leads to a shift in the Fermi energy. Thus for higher concentrations of nitrogen, the range for which same-sublattice doping is preferred is reduced due to the presence of the β S oscillations in Fig. 7 (c). However, the strength of the preference within this range may be increased by the slower rate of decay predicted for doped graphene.

B. Impurity Segregation
Recent experimental observations, corroborated by DFT calculations, seem to suggest that two nitrogen impurities in close proximity to each other prefer to occupy sites on the same sublattice in a quasi-neighbouring configuration 13,14 as shown in Fig. 4 (Left, Inset), hereon referred to as N •• 2 configuration. This configuration is preferred over the nearest-neighbour configuration N •• 2 seen in Fig. 4 (Right, Inset) which was calculated as being the most favourable configuration in the simple model above. Despite the limitations of the simple model for small separations, we note that the experimentally observed configuration fits the general trend of small separation, same sublattice configurations being the most favourable that the asymptotic behaviour of our model suggests.
The real advantage to our FO-based approach to studying such systems becomes clear when we consider multiple N •• 2 type impurities. Experimental evidence suggests that not only do pairs of nitrogen impurities prefer the N •• 2 to the N •• 2 configuration, but that a pair of N •• 2 -like impurities prefer to locate on the same sublattice. In other words, that two N •• 2 or two N •• 2 impurities are formed in preference to one of each. This behaviour leads to domains in nitrogen-doped graphene with a large sublattice asymmetric doping. Such behaviour has been predicted to lead to interesting and useful transport properties 15 . Numerical investigation of such systems using DFT calculations is limited to small separations which makes it difficult to explore the behaviour which emerges in more highly-doped larger scale systems. By extending the model discussed for single substitutional N dopants in the previous section, we can use the Lloyd model to investigate these types of systems. In this setup, the individual impurities are now N •• 2 or N •• 2 defects, shown as insets in Fig 7 (e). We consider a N •• 2 -type impurity at location A and introduce a second N •• 2 or N •• 2 impurity a distance D away at site B. We calculate the CEF for such a configuration analogously to Eq. (17), where now ∆E N 2 1 is the total change in energy in introducing a single N •• 2 or N •• 2 impurity. This quantity is plotted for the case when the second impurity is also a N •• 2 (blue curve) and when it is a N •• 2 (red dashed curve) in Fig. 7 (d). We note that, similar to the substitutional impurity case, same sublattice impurity configurations are preferential.
To benchmark our calculations, it is worth comparing our ∆E N 2 AB results to the DFT calculation performed in Lv et al. 31 for a single value of separation. In this work, calculations were performed for both two N •• 2 -type impurities and for a configuration with one N •• 2 and one N •• 2 . In both cases the impurity pairs had a separation of approximately D = 7. An energy difference of 14meV is reported using the DFT calculation 31 compared to 2.3 meV for the tight-binding model. In both cases the double N •• 2 configuration was energetically favourable. The numerical discrepancy between the results is to be expected due to the overly simple paramterisation of the N •• 2 impurity employed in the tight-binding model. However, the qualitative results for this model are not strongly affected by the local impurity parameterization, indicating also that the long-range sublattice ordered doping behaviour may not be unique to nitrogen. We emphasise that the same-sublattice configuration preference is noted for all separations in our model, explaining the long-ranged ordering seen in experiment 13,14 .
In a similar manner to that discussed for substitutional impurities, a finite concentration of N •• 2 impurities shifts the Fermi energy away from the Dirac point and introduces oscillations in β 2N (D). These oscillations, seen in impurity. However, small increases in E F would preserve the same sublattice preference in local regions. This suggests that larger scaled N-doped systems may have alternating domains where each of the sublattices is dominant, and this prediction is consistent with experimental observations 14 . An extension of the model discussed here to include a more accurate parameterization of the individual N •• 2 impurities would provide a transparent and computationally efficient method to explore the formation and size of such domains, and to determine their dependence on the concentration of N dopants and the resultant Fermi energy shift.

V. CONCLUSIONS
In this work we have derived analytic approximations for change in carrier density (∆n) and local density of states (∆ρ) in the high symmetry directions in graphene by employing a Stationary Phase Approximation of the lattice Green Functions. We obtain excellent agreement with numerical calculations for single and double substitutional atoms, vacancies, top adsorbed and bridge adsorbed impurities in the long-range limit, finding ∆n decays with distance (D) as D −2 for all impurities for E F = 0. At the Dirac Point, due to the disappearing density of states, we find the Friedel Oscillations in ∆n away from substitutional, top and bridge adsorbed impurities decay with as D −3 , but that in the case of vacancies and a resonant top adsorbed carbon ∆n is unchanged from the pristine case on all lattice sites due to the symmetry of electron and hole states in the LDOS profile around the Dirac Point. In the case of top adosrbed carbon, which is less energetically favoured than the more naturally occuring bridge adsorbed carbon, the cross sublattice interference present in the bonding arrangement ensures ∆n does not vanish and can be seen through the Green Functions of the system. Furthermore by expressing the total change in system energy due to the introduction of impurities through the lattice Green Functions we investigated how a sublattice asymmetry of both single and pairs of nitrogen dopants in graphene arises. We demonstrate that the dopant configuration energy is minimised where they share the same sublattice, a result which agrees with recent experiments [13][14][15] where such a distinct sublattice preference was found.
It is possible for our method to be extended further by applying it to strained graphene, as has been done with the SPA approach to the RKKY 32 and the behaviour of FOs in such a system has been studied theoretically only very recently 5 where a change in decay behaviour and sublattice asymmetry of the FOs due to the merging of the two inequivalent Dirac Points in the Brillouin Zone caused by inducing strain was found.

VI. ACKNOWLEDGEMENTS
Authors J.L. and M.F. acknowledge financial support from the Programme for Research in Third Level Institutions (PRTLI). The Center for Nanostructured Graphene (CNG) is sponsored by the Danish National Research Foundation, Project No. DNRF58.