Selective nitrogen doping of graphene due to preferential healing of plasma-generated defects near grain boundaries

Hyperspectral Raman IMAging (RIMA) is used to study spatially inhomogeneous polycrystalline monolayer graphene films grown by chemical vapor deposition. Based on principal component analysis clustering, distinct regions are differentiated and probed after subsequent exposures to the late afterglow of a microwave nitrogen plasma at a reduced pressure of 6 Torr (800 Pa). The 90 × 90 µm2 RIMA mapping shows differentiation between graphene domains (GDs), grain boundaries (GBs), as well as contaminants adsorbed over and under the graphene layer. Through an analysis of a few relevant band parameters, the mapping further provides a statistical assessment of damage, strain, and doping levels in plasma-treated graphene. It is found that GBs exhibit lower levels of damage and N-incorporation than GDs. The selectivity at GBs is ascribed to (i) a low migration barrier of C adatoms compared to N-adatoms and vacancies and (ii) an anisotropic transport of C adatoms along GBs, which enhances adatom-vacancy recombination at GBs. This preferential self-healing at GBs of plasma-induced damage ensures selective incorporation of N-dopants at plasma-generated defect sites within GDs. This surprising selectivity vanishes, however, as the graphene approaches an amorphous state.


INTRODUCTION
Numerous industrial applications of graphene rely on the largearea synthesis, which is commonly achieved by chemical vapor deposition (CVD) 1 . This method leads to the formation of numerous grain boundaries (GBs) between graphene domains (GDs). The former has been extensively studied since they are extremely different from GDs with respect to their electronic 2-4 , mechanical 5 , magnetic 6 , and chemical 7,8 properties. GBs also have shown both etching enhancement 9,10 or self-healing 11 , depending on the irradiation conditions. The understanding of extended defect topology 12 at GBs is limited and quickly vanishes when various growth conditions are considered during irradiation of graphene 13,14 . As pointed out by Malola et al. 7 "Grain boundaries […] are like snowflakes-there is no flake like another"; this clearly highlights their inherent complexity due to a sheer number of different possible GB configurations, each exhibiting its own intrinsic properties 7,[15][16][17] . This thus makes difficult the study of damage formation by ion or plasma irradiation of polycrystalline graphene films in presence of GBs.
Recently, advances in a non-intrusive spectroscopic technique called hyperspectral Raman imaging (RIMA TM -Photon ETC.) have been applied to study the inhomogeneity of plasma-induced disorders in the graphene lattice 18,19 . This method has been able to detect strong areal differences in polycrystalline monolayer graphene samples grown by CVD that were supposedly highly uniform. Being able to study spatial discrepancies is fundamental since it is well-known that the functionalization of graphene differs according to local defects initially present in the pristine lattice, be it GBs, dislocations, or impurities such as PMMA residues.
Inspired by ref. 18 , the present study capitalizes on RIMA capacities to examine the respective Raman signatures of GDs and GBs following the exposure of polycrystalline monolayer graphene to the late afterglow of a microwave nitrogen plasma at reduced pressure. More precisely, a method based on principal component analysis (PCA) filtering of hyperspectral Raman mappings 19 is developed and used to analyze behaviors observed in different regions of the CVD-grown graphene samples. The method is able to differentiate clusters that present instances of similar Raman signatures. Thus, it can probe the late-afterglow nitrogen plasma treatment effect on these regions 20,21 . Coupled with X-ray photoelectron spectroscopy, the doping level of the plasmatreated sample is assessed. By monitoring local variations in the initial pristine state and in the resulting plasma-treated state of the same graphene film, the study sheds light on the dynamical behaviors of N and C adatoms during plasma treatment. Using recent literature, this analysis brings a deeper understanding of doping selectivity in nitrogen plasma treatments.

Imaging and clustering
Graphene samples (grown on copper foil by CVD 22 and then transferred on SiO 2 /Si substrates using a standard PMMA procedure 23 ) were exposed to the late-afterglow region of a reduced-pressure plasma sustained by microwave electromagnetic fields 20,21 . A schematic of the plasma reactor is presented in Supplementary Fig. 1. In this study, the sample undergoes subsequent plasma treatment steps between which the sample is probed by X-Ray photoelectron spectroscopy (XPS) (spot size of 400 microns) and over a 90 × 90 μm 2 area using a Raman Imager (RIMA TM ) from Photon ETC 19 . The spatial resolution of Raman measurements is 390 nm. The RIMA allows the acquisition of a high number of spectra (here 116,281) from which parameters of the different bands are extracted. The G (~1580 cm −1 ) and 2D (~2700 cm −1 ) bands are prominent features of the untreated sample, while the D (~1350 cm −1 ) and D' (~1600 cm −1 ) bands rise with the generation of plasma-induced disorders 24 . Typical RIMA measurements, as well as examples of data preparation and processing, are presented in Supplementary Fig. 2 and Supplementary Figs. 3 and 4, respectively. PCA-based clustering is performed to highlight the differences between the regions probed inside the RIMA probed area. Scanning electron microscopy (SEM) is carried out after the last plasma treatment to assess the topographical state of the different clusters over the polycrystalline monolayer graphene film. Figure 1a, c presents SEM images at different magnifications over specific regions of the graphene sample after the final plasma treatment. The clusters identified by PCA analysis over the same regions are shown in Fig. 1b, d. We notice the dominance of the green cluster (72%) matching the graphene domains (GDs). The red (2%) and black (7%) clusters appear to be prevalent at linear defects, seemingly grain boundaries (GBs), and defects in their vicinity. Magenta (5%) and blue (15%) regions are mostly local defective graphene spots (GSs) at the center of the graphene domains. Such comparison thus demonstrates the capability of Raman spectroscopy with RIMA when coupled with a clustering method to highlight the local variation of supposedly uniform monolayer graphene films. Since some of the GSs appear to be indistinguishable from GDs in the SEM images, the clustering technique permits not only to reveal these defective regions but also allows their detailed characterization before and after plasma irradiation. Similar features can be seen by optical imaging over the whole area probed by RIMA analysis. In order to avoid content duplication, the results are presented in Supplementary Fig. 5.
SEM images presented in Fig. 1a-d reveal several topographical features often coinciding with the Raman mapping, in particular the presence of darker lines linked to the boundaries of the various growth domains in polycrystalline graphene films. These lines linked to GBs can also be seen in the comparison between optical image and RIMA cluster mapping presented in Supplementary Fig. 5. Note that, due to the width discrepancy between the Raman pixels (390 nm × 390 nm) and the GBs (typically 2-3 nm wide 25 ), details seen with SEM (electron beam diameter~1 nm) may not always be observed in the Raman mapping. Another interesting feature of the SEM images displayed in Fig. 1a, c is the presence of darker circular areas (circled in orange) throughout the whole graphene surface. Typically organized following straight lines, these have a consistent size of~1 µm in diameter. They also exist in two distinct shades: dark or light gray. As can be seen in Fig. 1c, these regions tend to be crossed by crack-like fractal linear defects emanating from GBs. While they resemble holes in SEM images, their Raman signature (red and black clusters, see all details below) is similar to GBs. These circular features were also examined by Energy-Dispersive X-ray Spectroscopy; the results are shown in Supplementary Note 1. On these darker circular areas, a rise of carbon by about a factor of two with respect to GDs was observed. No significant change in the other elements was seen. This indicates that these dark and light gray spots represent bilayer graphene domains.
Another distinct feature is the presence of numerous isolated dots (highlighted in a yellow circle) spread in a rather uniform way throughout the graphene surface. Considering the surface sensitivity of SEM at 3-kV acceleration voltage (electron penetration depth of about 150 nm in SiO 2 26 ), this could correspond to contaminants introduced either during the growth by CVD and/or through the transfer process of the graphene sample onto SiO 2 /Si substrate. Most likely, these dots are mostly contaminants commonly seen in CVD-graphene grown on copper in quartz furnace due to a devitrification of the quartz tube 27 . When untreated, these carbon-coated SiO x -based particles appear as small white dots uniformly spread on the untreated surface of the graphene. This aspect was confirmed by Energy-Dispersive X-ray Spectroscopy; the results are also shown in Supplementary Note 1. After plasma treatments and laser exposure (from Raman imaging) of the region of the monolayer graphene film that becomes increasingly amorphous, we noted a change in the morphology of these contaminants (white dots become darker). Graphene surrounding these defects should exhibit modified strain and/or doping, which signals are different in Raman with respect to GDs. Here, these contaminants are linked to the GSs component of the Raman mapping (blue and magenta clusters).

Damage generation
To easily assess relevant properties of CVD-grown graphene films before and after plasma treatment, Raman band parameters can be plotted using normalized 2D histograms such as in Fig. 2a. A maximum count normalization is made, and the observed contours are set at 10% of the maximum for each cluster. This allows to easily highlight how the properties of the graphene differ from a region to another. The red and black clusters present lower G band frequency (ω G ) and higher compressive strain (i.e. higher 2D frequency; ω 2D ) than the green regions. This is typical of CVD-grown GBs 28 and strengthens the aforementioned observation deduced from Fig. 1. The GSs (blue and magenta clusters) and GDs (green) present similar ω G and ω 2D behaviors with a tail at slightly lower values for the former. However, multiple factors can influence band frequency, including doping 29 and/or damage 30 . As discussed below, further investigation is required to identify the nature of these clusters. For comparison purposes, Fig. 2b presents a similar 2D histogram for the full region probed without any cluster differentiation (standard method). In such case, the data do not permit the study of low-count regions. While Raman mappings could help distinguish two regions for a given peak parameter, it cannot do so without parametrization along with two different band parameters. Mean values of Raman band parameters for each cluster are provided in Supplementary Table 1I. A well-established method of characterization of damage generation and damage type in graphene is to plot the ratio of the area of the D and G bands (A D :A G ) as a function of the bandwidth of the G band (Γ G ) 31 . Figure 3 presents such analysis for different successive plasma treatment times of the same graphene sample: t = 0 (a), 60 s (b), 120 s = 2 × 60 s (c), 180 s = 3 × 60 s (d), and 240 s = 4 × 60 s (e). The top orange curve is related to pure 0D defect generation (point defects), while the bottom one corresponds to pure 1D defect generation (line defects and reduction of the crystallite sizes) 31 . Due to the change in the coverage ratio of structurally damaged area and activated area, both damage generation pattern (0D or 1D) reveal different Raman signal evolution as the disorder increases. For 0D defects, a sharp rise is followed by a decrease of A D :A G as a function of Γ G as structurally damaged areas start covering activated areas of nearby defects. This occurs for low inter-defect distances (L D ). The 1D defects show a slower and monotonous increase of A D :A G as a function of Γ G as the size of nonacristaline (L a ) domains decreases.
In Fig. 3a, it can be seen that the GDs (green) of the untreated graphene sample present a very narrow distribution at low Γ G and A D :A G , whereas GBs (red and black) show much larger Γ G while keeping a low A D :A G . The GS (magenta and blue) data exhibit a rather large A D :A G ratio with a rather high and broadly-distributed Γ G . After the first plasma treatment, Fig. 3b reveals a strong increase of A D :A G for the GDs, while Γ G is rather constant. This corresponds to a strong 0D damage generation to reach about halfway of the 0D-type defect curve. Regarding the GSs, they show a lower A D :A G with a slightly higher Γ G . Most notably, the GBs seem to present a strong sturdiness to the plasma treatment as the A D : A G values of the red and black clusters are much lower than the ones of the GDs. This conclusion cannot necessarily be extrapolated to any type of sample as growth conditions can alter the properties of GBs 12 . Since each pixel probes, an area far greater than the actual width of GBs (390 × 390 nm 2 versus 2-3 nm 25 ), the black and red clusters correspond to areas containing various densities and/or types of GBs 7 .
After two subsequent 60-s treatments (Fig. 3c), the A D :A G ratio from GDs reaches the maximum of the empirical 0D defect generation curve (y~110 eV 4 , x~25 cm −1 ). At this point, the sample has reached the end of the first stage of damage generation (stage 1 [30][31][32] ). Further damage progressively causes the amorphization of the graphene (stage 2 [30][31][32], which is characterized by global disorder with an increasing amount of sp 3 C-C bonds. At the end of stage 1, GBs still show good resistance to the plasma-induced damage, exhibiting both smaller values of Γ G and A D :A G ratios. The last two plasma treatments (Fig.  3d, e) reveal very similar results. As the sample is brought closer to an amorphous state, additional plasma treatments have seemingly less effect on the Raman band parameters 30,32 . The GDs and GSs show similar signals even though a slight decrease in their A D :A G ratios and an increase in their Γ G (up to 45 cm −1 ) are observed. On the other hand, GBs show a lower increase in their Γ G , which is particularly noticeable in the red cluster data. Figure 3f presents a summary of the clusters' damage generation over all the treatment steps where the mean value for each step and cluster type is plotted (error bars correspond to the standard deviation). The maximum A D :A G value for GBs is much lower than for GDs, whose damaging behavior is close to the empirical 0D damage generation curve (top green curve). This highlights the differences in the damage generation: for GBs, a change in both crystalline sizes (L a ) and inter-defect length (L D ) is present, while GDs essentially see a change in L D 31 . The assessment of the damage with the D:G band ratio for monolayer graphene can be misleading due to its simultaneous dependence on charge carrier density. Indeed, values derived from D:G plots decrease for increasing doping (both positive or negative) 33,34 . In contrast to the D:G ratio, the D:2D ratio steadily increases with rising lattice disorder 35 , and its value remains independent upon increasing strain and doping 34,36 . This explains the recent use of this marker in damage assessment studies in graphene. In this work, the band ratio D:2D is therefore used as a direct indicator of the disordered state of the graphene lattice. Furthermore, the 2D:G ratio is known to be very sensitive to perturbation of the pristine graphene lattice; its value sharply decreases with increasing disorders of any types 37,38 , including GBs 39 . The 2D:G ratio also decreases for bi-and multilayer graphene 40 . While 2D:G is also known to be influenced by the doping level 34 , this ratio (unlike D:G) decreases monotonously with increasing damage and thus facilitates the interpretation of the results. In this framework, Fig. 4 presents the evolution of the 2D:G ratio as a function of the corresponding D:2D ratio.
As can be seen in Fig. 4, almost all clusters present large 2D:G for the untreated sample. This implies the growth of good quality graphene. As expected, a decrease of 2D:G of similar amplitude is found all over the sample after the plasma treatment, but the increase in D:2D is larger for clusters linked to GD and GS regions. This difference of behavior hints towards inhomogeneities in the variation of strain and or doping. Thus GBs seems to undergo a much less aggressive disorder generation: the increase of D:2D from GB-related spectra is two times lower than the one of GDs. Regarding 2D:G, each cluster presents a sharp drop with increasing plasma treatment time with black and red clusters linked to GBs experiencing lower decreases.
Optical strain and doping decoupling The band frequency variations of G (ω G ) and 2D (ω 2D ) can provide significant insights into the doping and strain levels of the graphene film. Indeed, the 2D band shows a large shift with a change in lattice strain while its counterpart varies strongly with the graphene doping level 29,34 . However, the method as proposed by Lee et al. for separating strain and doping 29 does not necessarily translate to damaged samples such as those  investigated here since ω G is also related to defect density for amorphous, nanostructured, and diamond-like carbon sp 2 materials 30 . Based on the work of Bruna et al. 34 , it is possible, however, to extract the behavior of ω G and ω 2D as a function of charge carrier density for damaged graphene (Supplementary Note 2). The linear behavior of ω 2D and ω G highlighted by Lee et al. 29 is lost for pdoped damaged graphene. A complete characterization of the band frequency shifts would therefore be required to decouple strain from doping for graphene undergoing increasing disorder. In the current work, the correlation extracted from Bruna et al. 34 is used instead since the plasma-induced disorder is relatively important. In such conditions, for n-doped graphene, ω 2D becomes proportional to ω G (Supplementary Note 2). Since we expect the nitrogen plasma to induce n-type doping, this modified method is better adapted, as discussed below. Moreover, Bisset et al. 28 revealed the importance of the crystallinity of the sample in the analysis of the ω G and ω 2D variations in response to strain. Indeed, highly polycrystalline CVD-grown samples (crystallite size L a~1 µm) show an increase of ω G for compressive strain at GBs, as opposed to the decrease observed in exfoliated graphene (L a~2 0 µm). Since their laser beam diameters are of the order of the CVD-grown graphene domain size, this implies that the x-scale proposed by Lee et al. should be reversed for GB-related data, i.e., for red and black clusters.
In line with these studies, Fig. 5a presents the variation of ω 2D and ω G mean values for all the regions and between the subsequent plasma treatment steps. For ω G , all regions initially show a small redshift followed by a large blueshift. For ω 2D , two distinct behaviors are observed: GBs show a monotonous increase while all other regions undergo a redshift followed by a large blueshift. Figure 5b, c presents the change in doping and strain of each cluster as a function of the plasma treatment time. Such values are obtained by projecting the values of Fig. 5a along the charge carrier density axis introduced by Bruna et al. 34 , and the strain axis introduced by Bissett et al. 28 . Only the shifts of the mean values are considered (to best distinguish the clusters), and the error bars result from the propagation of the errors of the linear fits performed on the relation between strain or charge carrier concentration on Raman band positions (as extracted from the literature 29,33,34 ). In order for any discussion about variations in doping and strain to be relevant, one must first remove the expected effect of damage generation for these conditions. As demonstrated by Eckmann et al. 41 , no clear change in ω G and ω 2D should be observed for damaged graphene before the amorphization stage. Entering stage 2 of the amorphization trajectory, the 2D and G bands show large redshifts and blueshifts, respectively. In addition, the effect is much stronger for the 2D band and thus the projection may bias the calculation of the strain from Fig. 5b, but only for the last two treatments and for the heavily damaged clusters (GDs and GSs).
Based on this analysis, the strain and doping behaviors presented in Fig. 5b, c reveal two distinct regimes. The first stage corresponds to a transition from graphene to nanocrystalline graphene (stage 1) and the second to the transition towards amorphous carbon (stage 2) 30,32,41 . The same transition is highlighted in Fig. 3 (transition towards amorphous carbon is responsible for a decrease of D:G at high disorder levels 31 ). More specifically, the first three graphene states along the amorphization trajectory (plasma treatment times: 0, 1, and 2 min) fall into stage 1, in which doping and strain behaviors differ among the different clusters. Most noticeably, the doping of GBs remains constant, while the concentration of charge carriers decreases for the other regions of the sample. This feature provides the first hint of more selective nitrogen incorporation at GDs (so-called segregation of nitrogen dopants at GDs 42 ). In addition, since both p-and n-type dopings increase ω 2D and ω G , a decrease along the axis used for the so-called strain-doping decoupling is related to either a decrease of p-or n-type doping. In Fig. 5c, this doping reduction for GDs and GSs reaches a maximum between 1 and 2 min of plasma treatment. It is usually expected that untreated graphene films grown by CVD show impurity-based (unintentional) p-type doping 43,44 . Therefore, it can readily be expected that such initial p-type doping first disappears before n-type doping from nitrogen incorporation can appear.
As for the strain evolution during stage 1, it can be seen in Fig. 5b that all clusters sustain a compressive strain with GBs suffering most of that strain. It is worth mentioning that nitrogen incorporation is expected to induce a mild compressive strain to the graphene due to its inferior bond length (C-C: 1.42 Å versus C-N:1.37 Å 45 ). Thus, a doping-related strain is expected to be relatively small for low-to-medium nitrogen incorporation levels 46 . Therefore, the strain features observed in Fig. 5b are most likely due to plasma-induced disorder than to plasma-induced doping in monolayer graphene films.
As discussed above, during the transition toward amorphization (stage 2), ω 2D and ω G are expected to redshift and blueshift, respectively 41 . This effect is not considered for the optical separation of strain and doping as it would require criteria to access the transition towards amorphization. As detailed previously and according to Vinchon et al. and Bruna et al. 11,34 , D:2D is independent of doping, and thus it represents an ideal candidate to discuss qualitatively the effects of an increase of damage on the expected change of the calculated strain (Lee et al.'s model 29 ). According to Eckmann et al. 41 , a rise in lattice disorder (in stage 2) leads to a decrease of ω 2D and an increase of ω G with a slope of around Δω 2D /Δω G~− 1.5. This would translate in an overall rise of tensile strain (see Fig. 3b of ref. 29 ), which is consistent with the presence of defects (such as vacancies) allowing strain relaxation. In our case, Fig. 5b does reveal a rise in tensile strain prior to the amorphization between t = 1 and 3 min. In addition, such an increase is not present for the GB-derived clusters (red and black): only a reduction of the decreasing slope is observed. This may be explained by their lower degree of disorder (lower D:2D ratios, Fig. 4), allowing for less strain relaxation than GDs and GSs. Moreover, such an increase of compressive strain might arise due to numerous adatom incorporations leading to inverse Stone-Wales defects 47 . As for the doping in stage 2, it can be seen that the estimated doping level for all clusters follows a similar rise. This rather spatially-uniform change in doping over the entire graphene surface is linked to the transition towards amorphization. More details are provided below.
X-ray photoelectron spectroscopy (XPS) is performed to further support the incorporation dynamics of nitrogen atoms extracted from the decoupling of strain and doping. For each plasma treatment step, XPS is used to extract the N:C ratio and the various components of the N1s high-resolution spectrum. Survey spectra reveal a moderate amount of nitrogen incorporation following subsequent exposure of the graphene films to the late afterglow of microwave N 2 plasmas (N1s/(N1s + C1s) = 3.7%, 5.7%, and 7.6% at 1 min, 2 min, and 3 min of plasma treatment, respectively). XPS surveys are provided in Supplementary Fig. 6. High-resolution spectra of C1s and N1s are presented at Fig. 6a, b. The C1s region presents the main sp 2 C-C bond typical of graphene at 284.6 eV 48-53 as well as additional features at higher binding energies related to N bonding (e.g., sp 2 C = N at 285.5 ± 0.1 eV 53-55 and sp 3 C-N at 287.5 ± 0.1 eV [55][56][57] or O (C-OH at 286.5 ± 0.2 eV 50-52,57 , C = O and/ or O = C-O at 288.5 ± 0.2 eV 48,51,55,56,58 and O = C-OH at 289.4 ± 0.2 eV 57 ). The N1s region reveals a strong contribution of the pyrrole-nitroso type of inclusions at 400.1 eV 53,59-62 with a reasonable amount of pyridine (399.2 eV 53,54,62,63 ) and graphitic incorporations (401.3 eV 53,54,[61][62][63]. Oxygen inclusion arises due to the exposition to atmospheric conditions of the plasma-related sample between each measurement. Its presence is limited for the untreated sample and is most likely caused by surface contamination 21 . Plasma-induced disorder creates dangling bonds responsible for oxygen inclusion after the treatment and further exposition to ambient air. This contamination is not believed to play a major role in the change in charge carrier concentration or its homogeneities across the graphene surface. Each nitrogen inclusion is expected to contribute differently to the modification of the Fermi level in monolayer graphene films. While graphitic inclusions induce n-type doping (0.54 e/N), pyridine provides p-type doping (0.45 p/N) 64 . Pyrrole inclusions are less discussed in the literature, but one expects their effects to be limited for low-to-medium nitrogen contents 65 . The product of the nitrogen atomic fraction (N1s/(N1s + C1s)) from the survey scans by the respective percent contribution of each nitrogen component from the high-resolution N1s scans is used to estimate the corresponding charge carrier concentration for each plasma treatment step. These results are compared in Fig. 6c to those extracted from Raman optical decoupling of strain and doping, for GDs only. Here, the absolute change in charge carrier density obtained from Raman analysis and presented in Fig. 5c was divided by the atomic density of carbon atoms in graphene domains (3.85 × 10 15 cm −2 ) to obtain percent variations. (i) Because XPS measurements are spatially averaged, it is impossible to distinguish in XPS the different regions or clusters (as in RIMA), and (ii) GDs represent more than 70% of the graphene surface (see Fig. 1 and Supplementary Fig. 5), it is therefore expected that XPS data are mostly due to GDs. As can be seen in Fig. 6c, similar trends and values are obtained from both XPS and RIMA. However, it is worth highlighting that due to the similar nitrogen content of both graphitic and pyridinic inclusions and their relatively low percent contribution in the overall N1s signals, the corresponding errors in XPS data analysis can be large.

DISCUSSION
Further analysis of the Raman band parameters displayed in the previous section can be used to better understand how each region of the graphene film is altered by the plasma treatment. Despite the use of a reactive plasma source, the analysis proposed by Cancado et al. 31 still holds and reveals the creation of 0D defects at the GDs (Fig. 3). In most plasma irradiation studies, such defects are linked to ion bombardment of the graphene sample. In the late-afterglow region of the microwave nitrogen plasma, however, the kinetic energy of positive ions interacting with the graphene sample (of around 0.5 eV, Supplementary Note 3) is too low to cause any relevant damage. Since pure knock-on collisions can be ruled out in such conditions, other damage formation mechanisms involving local potential energy transfers must be involved 66 , including the surface neutralization of N 2 + −16 eV, the surface deexcitation of N 2 (A) metastable species −6 eV, and the surface recombination of nitrogen atoms −10 eV 67 . Similar energy exchange processes were proposed to reduce the energy barrier by a factor of 4 (from 9 eV to 2.3 eV) for Frenkel pair formation in graphene-like structures under a carbon adatom flux 68 . Note that, over the range of experimental conditions examined in this study, high-energy photons emanating from the main microwave plasma zone cannot reach the lateafterglow region due to the use of an adequately-shaped knee in the discharge tube ( Supplementary Fig. 1).
In addition to energy transfer processes, residual oxygen species present in the microwave plasma and flowing afterglow (as a result of base pressure impurities and outgassing of plasma reactor walls) could also induce mild chemical etching of carbon atoms in monolayer graphene films. Such a phenomenon is expected to be mostly specific to GBs, as previously observed in O-and H-bearing plasmas 9,[69][70][71] . Surprisingly, GBs reveal less damage generation than all other regions: this is evidenced by a much weaker D:2D ratio at GBs than in GDs. In addition, D:D' ratios are constant for all clusters ( Supplementary Fig. 7), suggesting that there is no preferential defect type generation 32 . Thus, damage formation must either be slower at GBs due to their intrinsic organization or there is a preferential self-healing of plasma-induced damage in these regions. In plasma irradiation conditions leading to a large density of carbon adatoms, Vinchon et al. 11 recently revealed that the GBs were more resilient than GDs due to a more efficient Frenkel pair recombination in their vicinity. Briefly, such preferential healing of plasmagenerated defects near grain boundaries results from (i) the difference between the migration energy of carbon adatoms (0.4 eV) and vacancies (1.6-3.0 eV) 72 and (ii) the anisotropic transport of these 0D defects along the axis of GBs 17,73 . This induces an accumulation of carbon adatoms at the GBs, which enhances the probability for adatom-vacancy recombination in the vicinities of GBs. Such a mechanism would obviously result in a net loss of vacancies close to the defect sites of black and red clusters, which explains the low D:2D ratios observed for the latter. This mechanism is expected to vanish as the sample gets more damaged due to the limited transport of C adatoms towards GBs (carbon adatoms now become mostly trapped by disorders in the GDs). Anisotropic transport of carbon adatoms along GBs is also strongly reduced as the graphene approaches the amorphous state. In such conditions, carbon adatoms coverage is expected to become much more uniform. Due to the limited amount of data, however, this transition toward the amorphous state cannot be studied in detail.
As mentioned above, the change in doping extracted from the method presented by Lee et al. 29 (updated with measurements of Bissett et al. 28 and Bruna et al. 34 ) reveals much larger doping at GDs than at GBs (Fig. 5c). Such dopant segregation at GDs is revealed, while some disorder is induced at GBs. This selective incorporation of nitrogen atoms at GDs only occurs for the first two plasma treatments, i.e. when graphene is still in stage 1 along the amorphization path. The explanation is twofold. First, as evidenced by the observed variations in D:2D ratios (Fig. 4), GDs sustain more damage overall than GBs, the latter being subject to a preferential self-healing phenomenon 11 . This results in a larger population of vacancies at GDs, which are favored sites for N-atoms incorporation. Second, a crucial aspect of dopant selectivity is the large difference between the migration barrier of carbon and nitrogen adatoms on the graphene surface: 0.4 eV for C 72 and 1.1 eV for N 74 . This means that C adatoms are highly mobile on monolayer graphene films over the range of experimental conditions investigated here (~room temperature 75 ), while N-adatoms are much less mobile. Thus, considering the C adatoms and vacancies behaviors in the study of Vinchon et al. 11 , a strong population imbalance is expected to arise throughout the sample between C-and N-adatoms with a greater density of carbon adatoms at the GBs. Eventually, as the sample engages its transition towards amorphization (stage 2), the incorporation becomes more homogeneous throughout the whole graphene surface and affects all clusters in a very similar way (Fig.  5c). This inhibits the preferential self-healing mechanism such that the N-incorporation dynamics in stage 2 becomes less selective.
It is worth highlighting that the amplitude of the nitrogen content (3-8%) is similar to what is typically observed for other treatments of CVD-grown monolayer graphene films (1-16%) [76][77][78] . The assessment of the electronic doping is not always performed, specially its spatial distribution and its behavior for increasingly disordered graphene. Zhao et al. 42 have shown the synthesis of Ndoped graphene by CVD and have revealed localization of nitrogen incorporation at graphene domains, where a much lower 0.4% graphitic content is shown responsible for n-type doping. Here, plasma-generated nitrogen atoms are found in many more inclusion configurations, but specific incorporation in the domains is still present even for the much larger nitrogen content.
A combination of Hyperspectral Raman IMAger (RIMA) and PCA clustering was used to distinguish different regions of a polycrystalline graphene film grown by CVD and exposed to the late afterglow of a microwave N 2 plasma at low pressure. The precision of the technique, verified by optical microscopy and scanning electron microscopy, is effective in highlighting the evolution of different regions of the sample: graphene domains (GDs), over-and subgraphene contaminants (GSs) as well as grain boundaries (GBs). Through careful decoupling of strain and doping effects in Raman spectroscopy, selective doping by plasma-generated nitrogen species is revealed at the GDs for the first two plasma treatments. Such dopant segregation is believed to be made possible by a combination of two phenomena. First, a preferential self-healing mechanism occurring at GBs leads to a decrease of the population of vacancies in their vicinity, the latter being a favored site for Natoms incorporation. Second, an imbalance is expected to occur resulting in a greater density of carbon and nitrogen adatoms at the GBs and GDs, respectively. After the third plasma treatment, the higher defect density throughout the entire graphene surface reduces the mobility of C adatoms toward the GBs, which stop the preferential self-healing mechanism, giving more homogenous Nincorporation across the different regions. Over the range of experimental conditions investigated, n-type doping is evaluated at 0.24%, a result confirmed by X-ray photoelectron spectroscopy. This powerful technique demonstrated the importance of the migration barrier of each plasma-generated defect; accumulating and incorporating more mobile defects in the vicinity of the GBs while defects with higher migration barriers segregate within GDs. These results and the characterization method presented in this study are of particular relevance for the understanding of surface processes in 2D materials and paves the way toward complete tailoring of the doping levels in graphene films.

Scanning electron microscopy
Scanning electron microscopy (SEM) were deliberately chosen not to use in between the different steps of plasma irradiation due to the possible electron beam-induced contamination effects 79 . Moreover, the increasingly damaged state of the film after plasma treatment might cause important surface charging effects as the sample gets more amorphous and thus perturbs subsequent Raman analysis. The SEM analysis was performed using a JEOL JSM-7600F in secondary emission mode at 3.0 kV; this allows for a good spatial resolution while ensuring a relatively small electron penetration depth. X-ray photoelectron spectroscopy X-ray photoelectron spectroscopy (XPS) is performed to assess the doping level of the polycrystalline monolayer graphene films. The setup is a Thermo Scientific K-Alpha (CAE detector, 180°double-focusing hemispherical analyzer, 128 channel detector) operating with pass energies of 200 eV and 50 eV for a survey and high-resolution scans, respectively. The step energy is 1 eV for the surveys and 0.1 eV for the high-resolution spectra. A flood gun is used to ensure minimal shifts due to charging. The spot size of 400 microns is centered on the area probed by the RIMA system.

Plasma treatment
A gap-type wave launcher (surfatron) maintains a 2.45 GHz surface wave along an 8-mm external (6-mm internal) diameter fused silica discharge tube. The whole apparatus is extensively described in the previous publications 20,21,80 . In this study, the pressure is set to 6 Torr (800 Pa), the gas flow to 100 standard cube centimeter per minute (sccm), and the injected power to 32 W. The resulting plasma length is about 2 cm. From the surfatron gap, the early afterglow peaks at 4 cm while the sample position is set at 27 cm. The sample is far enough from the energetic species of the plasma discharge and thus interacts with a much less damaging environment. Still, three different reactive species (plasmagenerated N atoms, N 2 (A) metastable species, and N 2 + ) are present in significant quantities in the afterglow and can induce a significant change of the graphene properties 20 . Their respective densities in the flowing afterglow region were previously determined by optical emission spectroscopy and Langmuir probe measurements 80 . The setup is presented in Supplementary Information S-I.

RIMA measurements
The setup is detailed elsewhere 18 . Briefly, a 3.25-W laser at 532 nm uniformly irradiates a 130 × 130 μm 2 area via a beam shaping device and a ×100 objective. The resulting Raman emission is then spectrally separated using a volumetric Bragg tunable filter. Despite the beam shaping device, the absolute value mappings of Raman bands still show a large difference between the side and the center of the detector. The line ratio is not affected but the signal to noise varies. In this study, to reduce the effect of a large variation of the signal-to-noise ratio between the center and the sides of the hyperspectral cube, 20-μm bands are cropped from each side resulting in a 90 × 90 μm 2 final probed area. The RIMA experiments can be time-consuming and thus a 3 × 3 binning is used in order to reduce by a factor 9 the length of each experimental run. Typical unfiltered spectra are presented in Supplementary Information S-II.

Data preparation and processing
The core of the processing methods employed in this article is detailed elsewhere 18 . The noise removing algorithm and baseline subtraction was improved. The methods are inspired by the work of Antonelli et al. 81 . The following points detail each step.

PCA-based noise filtering
First data centering and normalization are carried out. PCA decomposition (centering) is then applied to extract the first 30 components maps of the measurements. A Gaussian unmixing is used to differentiate the ten clusters core to this study. This is done to ensure that the following noise removal algorithm is able to properly remove the noise of the data of each cluster. Indeed, PCA is intrinsically biased by the intensity of each dimension and therefore mappings containing a low amount of highly distinguishable spectrum induce large error to the noise reduction of these outliers. Moreover, the number of components required to reconstruct all the data would be increased and could possibly reduce the quality of the noise reduction. Base on the work of Antonelli et al. 81 , we thus try to remove a uniform simulated noise contribution for the data while ensuring that the noise is similar for all spectra and is of adequate amplitude. We must first estimate the intensity of the noise to be removed. We do so using the modified criteria proposed in another work 81 This criterium is verified for each cluster in order to select the components up to the change in the slope of the average spectrum error. Extreme precision for this first evaluation of experimental noise is not mandatory as a correction of experimental noise is later assessed. We then perform a weighted average of the estimated experimental noise values. To converge toward a more precise value of the experimental noise, we then evaluate a minimum delta value for all clusters simultaneously (adjusting the root mean square of the error by the normalization factor of each cluster). (Eq. (26) from ref. 81 ) This allows to obtain the corrected experimental noise value. Through delta minimization again, we find the number of components in each cluster that allows to remove a constant noise of the appropriate amplitude. Filtered spectra are presented in Supplementary Information S-III.

PCA assisted artifact subtraction
The nature of the artifact of the Raman Imager (RIMA) is yet unknown. Its main cause is believed to be the fluorescence of the system objectives. The shape of the artifact can drastically change from one position to another on the CCD. Measurement is taken on a clean substrate of SiO 2 . Through PCA decomposition (centering, no normalization), four main components G. Robert Bigras et al.
are extracted. A linear combination of those can be used to describe the behavior of all the artifacts on the CCD.

Spectrum fitting
It results that four parameters are fitted simultaneously to the band parameters to ensure proper baseline substation. The constraints imposed by the 4 core-artifacts enable better fitting of the baseline on the extreme regions of the spectra, i.e., regions where polynomial fitting tends to diverge and somewhat induces error on band fittings. The fitting method is least-square fitting with centering, normalization, and dynamic bounds. Line shapes are set to Lorentzian for all bands. Spectra fitting and artifact subtraction can be assessed for typical spectra presented in Supplementary Information S-III.

PCA clustering
In order to highlight the differences between the regions probed inside the selected area of polycrystalline monolayer graphene films, a processing method was developed. All spectra from the Raman map are first centered using their mean values and normalized according to their standard deviations. PCA is performed using all spectral values as dimensions and after subtraction of the mean values of the spectra. The first 30 components are then separated into 10 clusters using Gaussian unmixing. The number of clusters is arbitrary and is taken high enough to make sure that the most distinguishable regions over the polycrystalline monolayer graphene films are separated into different clusters. The number of components taken for noise filtering is typically lower than 30 18 and thus the number of considered PCA components is high enough to ensure adequate distinction of the relevant areas. Since the number of components is arbitrary, some clusters are combined by comparing various band parameters (mainly D over G band intensity ratio as well as width and position of the G band). This allows for the definition of final clusters that highlight regions of distinct Raman signatures.
To investigate the evolution of each cluster as the plasma treatment steps are performed, an image registration (allowing only translation and rotation) is carried out to align every plasma-treated measurements relative to the initial region of the untreated graphene sample. This enables a precise study of the evolution of different graphene areas as a function of their initial pristine state (grain boundary, local defects, contaminants, etc.) 18 .