Atomistic Mechanisms of Codoping-Induced p- to n-Type Conversion in Nitrogen-Doped Graphene

It was recently shown that nitrogen-doped graphene (NG) can exhibit both p- and n-type characters depending on the C-N bonding nature, which represents a significant bottleneck for the development of graphene-based electronics. Based on first-principles calculations, we herein scrutinise the correlations between atomic and electronic structures of NG, and particularly explore the feasibility of converting p-type NG with pyridinic, pyrrolic, and nitrilic N into n- or bipolar type by introducing an additional dopant atom. Out of the nine candidates, B, C, O, F, Al, Si, P, S, and Cl, we find that the B-, Al-, and P-codoping can anneal even relatively large vacancy defects in p-type NG. It will be also shown that, while the NG with pyridinic N can be converted into n-type via codoping, only the bipolar type conversion can be achieved for the NG with nitrilic or pyrrolic N. The amount of work function reduction was up to 0.64 eV for the pyridinic N next to a monovacancy. The atomistic origin of such diverse type changes is analysed based on Mulliken and crystal orbital Hamiltonian populations, which provide us with a framework to connect the local bonding chemistry with macroscopic electronic structure in doped and/or defective graphene. Moreover, we demonstrate that the proposed codoping scheme can recover the excellent charge transport properties of pristine graphene. Both the type conversion and conductance recovery in codoped NG should have significant implications for the electronic and energy device applications.


I. INTRODUCTION
Because of its unique structural, electronic, and transport properties, [1][2][3][4] graphene is regarded as one of the best candidate materials to be incorporated into the next-generation electronic, energy, and bio devices. [5][6][7] To realize such potential for device applications, reliable methods to tailor the atomic and electronic structures of graphene is required. [8][9][10] Doping of graphene has been extensively investigated in this context, and among various options nitrogen doping emerged as one of the most effective schemes to improve diverse functionalities of graphene [10][11][12][13][14][15][16][17][18][19][20] and especially achieve n-type doping that is crucial for electronics applications. [21][22][23][24][25] Unfortunately, recent experimental and theoretical studies have revealed that N-doped graphene (NG) can assume both n-and p-type characters depending on the bonding nature of N atoms, namely, n-type for the graphitic N and p-type for the pyridinic, pyrrolic, and nitrilic Ns. [25][26][27] The question then arises as to how one can achieve robust n-type NG, and at the more fundamental level which factors determine the structure-property relationships of NG. The latter should have significant implications for practical applications. For example, the bonding state of the N atom was found to critically affect the oxygen Carrying out first-principles density functional theory (DFT) and DFT-based non-equilibrium Green's function (NEGF) calculations, we herein systematically investigate the atomistic origins of p-type character in NG with pyridinic, pyrrolic, and nitrilic N and the feasibility of achieving robust n-type graphene by incorporating a codopant atom. Graphene codoped with N and B or N and P has been already synthesised and shown to improve the ORR performance, [28,29] but the atomistic details and mechanism of the synergetic effects associated with the B, N-and P, N-codoping are not understood yet. We first consider the energetic feasibility of introducing an additional atom into various p-type NG defect sites and show that B, Al, Si and P atoms can structurally anneal even relatively large vacancy defects next to the pyridinic, nitrilic, and pyrrolic N atoms. We find that, except Si, the B, Al, and P codoping can convert the pyridinic NG into n-type and the nitrilic and pyrrolic NG into bipolar type, and thus effectively eliminate the p-type character of NG. Based on the Mulliken and crystal overlap Hamilton population (COHP) analyses, we further establish the basis to understand how the macroscopic electronic type change in NG is induced in the atomistic bonding viewpoint. Finally, we will demonstrate that the additional benefit of the codoping approach is the recovery of excellent charge transport capacity of pristine graphene, in line with the recovery of the sp 2 bonding network upon the healing of vacancy defects.

II. CALCULATION METHODS
We performed the spin-polarized DFT calculations within the Perdew-Burke-Ernzerhof parameterization of generalized gradient approximation [30] using the SIESTA software.
[31] The atomic cores were replaced by the Troulier-Martins-type norm-conserving pseudopotentials, [32] and the double-ζ-plus-polarization level numerical atomic basis sets defined by the confinement energy of 80 meV were adopted. For the supercell of 9 × 9 graphene unit cells (see Supporting Information Fig. S1), we used a mesh cut-off energy of 200 Ry and the 3 × 3 × 1 k-point sampling in the Monkhost-Park scheme. [33] A finer 30 × 30 × 1 k -point mesh was sampled for the calculation of density of states (DOS) and COHP. The vacuum region of the supercell along the direction perpendicular to the graphene plane was set to be 25Å. We have calculated the work function Φ using the equation: where E F is the Fermi energy and V vac is the macroscopic average potential in the vacuum, defined as the midpoint between the graphene layer and its neighbouring images. [34] The doping charge density was calculated according to where E D is the Dirac point energy and v F is the Fermi velocity, 1 × 10 6 m/s. [3,34] Transmission functions T (E) were calculated using the DFT-based NEGF method, [35] as implemented in Tran-SIESTA. [36] We used the periodic cell composed of six dimer lines along the transport normal direction (14.80 A) and sampled 66 k ⊥ points. Along the transport direction, we used eight and two zigzag chains to model the channel and electrode regions, respectively, and the surface Green's functions were obtained for the corresponding electrode models sampled with the 25 k points (see Supporting Information Fig. S1). In calculating T (E), the energy was scanned from −1.0 eV to 1.0 eV with respect to E F with the 0.001 eV resolution. In many cases, we have crosschecked the validity of our results using Se-qQuest and our in-house NEGF code. [37,38]

III. RESULTS AND DISCUSSION
In this work, in addition to graphitic N (N gr (Fig. 1a), we considered four other representative NG conformations: pyridinic N next to a monovacancy (V 1 -N py , Fig. 1b) and a divacancy (V 2 -N py , Fig. 1c), nitrilic N next to a divacancy (N nit , Fig. 1d), and pyrrolic N next to a trivacancy (N pyrr , Fig. 1e). [25,39,40] Recent experimental and theoretical reports showed that, despite high formation energies, vacancy defects can form within the graphene sheets by, e.g., ion or electron irradiation. [41] In the presence of vacancy defects, the pyridinic N configurations were shown to become energetically more favourable than the graphitic N one. [26,42]  In Table 1, we show the zero-temperature formation energies of the V 1 -N py , V 2 -N py , N nit , and N pyrr structures calculated according to where E N −doped/pristine are the total energies of Ndoped/pristine graphenes, µ C/N/H are the chemical potentials (total energies per atom in their elemental reference phases) of C/N/H, x is the number of C atoms removed from the graphene sheet during the vacancy formation, y is the number of N atoms, and z is the number of H atoms (2 for N nit and 0 otherwise). The chemical potentials of C, N, and H were extracted from graphene, N 2 molecule, and H 2 molecule, respectively. Note that we include two pyridinic conformations, because the V 2 -N py configuration is energetically even more favorable than the N nit and N pyrr counterparts. For the V 1 -N py and V 2 -N py cases, we additionally considered the substitution of different numbers of N atoms around the vacancy sites, V 1 -N py α (α = 1 ∼ 3) and V 2 -N py β (β = 1 ∼ 4), and found that the maximum concentration of pyridinic N atoms around the mono-and di-vacancies (trimerized V 1 -N py 3 , tetramerized pyridine-like V 2 -N py 4 ) is energetically preferred (Table 1). For V 1 -N py α and V 2 -N py β , we found that the maximum concentration of pyridinic N atoms around the mono-and di-vacancies (trimerized V 1 -N py 3 , tetramerized pyridine-like V 2 -N py 4 ) is energetically preferred (Table 1). [26] An important objective of the present work is to provide the atomistic understanding on the nitrogen bonding nature in NG and its modification upon codoping. For this purpose, we will demonstrate that the combination of bandstructure, DOS, Mulliken population, and COHP [43] provides a systematic and detailed information that guides the route toward the desired n-type or bipolar conversion. The COHP is defined by the DOS multiplied by the Hamiltonian of the corresponding element, and the negative COHP values (−COHP) give positive and minus signs for the bonding and antibonding states, respectively. The DOS, COHP normalised by the number of bonds, band structure data of N gr , V 1 -N py 3 , V 2 -N py 4 , N nit , and N pyrr are shown in Fig. 1, and their electronic structures are summarised in Table 2. The doping type of graphitic NG is n-type, i.e., E D is located below E F (Fig. 1a middle and bottom panels). This results from the electron transfer from the N gr atoms to the π* states of graphene. Visualising the Mulliken charge populations as in Fig. 2a, we can observe that the electrons donated from N gr are distributed throughout the graphene basal plane or become mobile charges. The amount of donated charges extracted from the population analysis is 0.402 electron per nitrogen atom, which is in excellent agreement with the experimental estimation. [44] It results in the upward shift of E F from E D of the pristine graphene case into the graphene conduction band, or the lowering of the work function by 0.37 eV. The calculated charge-carrier density of 11.31 × 10 12 cm −2 for our model with a N atom doping ratio of 0.62 % N atoms per C atom (N atom density of 2.31×10 13 cm −2 ) is in good agreement with the experimental estimation of 5.42×10 12 electrons per cm 2 for the 0.34 % N atoms per C atom doping. The COHP curve shows that these impurity resonant states have antibonding characters, being identified as the strong negative COHP peaks right above E F (thus E D ). The energetic locations of the antibonding COHP peaks are E F + 0.07 eV and E F + 0.45 eV. These peak locations depend on the supercell size, or the doping ratio, and is expected to converge to an experimentally observed single peak at about E F + 0.14 eV. [27,45] On the other hand, the electronic structures of the other three V 1 -N py 3 , N nit , and N pyrr NG are p-type, or E F has been shifted downward into the valence bands (Figs. 1b, 1d, and 1e middle and bottom panels). This should mainly result from the presence of vacancy defect states (whose LDOS are shown in Figs. 1b, 1d, and 1e top panels), which behave as acceptor states with missing π-electrons. [24,25,27] Note that in these cases the impurity states appear in the COHP plots as strong antibond peaks right below E F (thus E D ). Interestingly, the V 2 -N py 4 case shows a bipolar character (E D nearly coincides with E F ). [27] To understand such discrepancies, we analyzed Mulliken populations as shown in Figs. 2b, 2c, and 2d for V 1 -N py 3 , V 2 -N py 4 and N pyrr , respectively (The N nit diagram is similar to N pyrr and is not shown). They show that the spatial range of charge redistribution around the nitrogen-vacancy complexes discriminate the V 1 -N py 3 and N pyrr cases from the V 2 -N py 4 counterpart: We observe the depletion of electrons throughout the entire graphene basal plane in the former, but rather localised charge depletion around the nitrogen-vacancy site in the latter (Compare Fig. 2c with Figs. 2b and 2d). Related with these charge transfer characters, the impurity states strongly hybridize with the carbon 2p z orbitals or graphene π states (existence of an antibonding C 2p z -N COHP peak) in V 1 -N py 3 (Fig. 1b COHP panel, shaded region near E F ), whereas the hybridisation is very weak (absence of an antibonding C 2p z -N COHP peak) in V 2 -N py 4 ( Fig. 1c COHP panel, shaded region at E − E F ≈ 0.5 eV).
We now move on to consider the possibility of codoping the p-type NG. Since the presence of a vacancy defect is the precondition of the existence of V 1 -N py 3 , V 2 -N py 4 , N nit , and N pyrr sites, we focused on whether the vacancy defects can be healed by introducing a codopant. In addition to C, eight light elements, B, O, F, Al, Si, P, S, and Cl, were considered as codoping candidates. For the fully optimised structures, X+V 1 -N py 1∼3 , X+V 2 -N py 1∼4 , X+N nit , and X+N pyrr (X = B, C, O, F, Al, Si, P, S, and Cl), we have calculated the formation energies, (4) where E codoped is the total energy of codoped NG, µ C/N/H/X are the chemical potentials of C, N, H, and the codopant X, x/y are the numbers of removed/added C/N atoms with respect to the pristine graphene, and z is the number of added H atoms (2 for N nit and 0 otherwise). The chemical potentials of B, O, F, Al, Si, P, S, and Cl were obtained from the α-rhombohedral boron, O 2 molecule, F 2 molecule, face-centred cubic aluminum, diamond cubic Si, black phosphorus, orthorhombic sulphur (α-S 8 ), and Cl 2 molecule, respectively. The forma-tion energies relative to those of NG, which represent the energetic feasibility of annealing the vacancy defects of p-or bipolar-type NG by the codoping approach, are summarised in Fig. 3  First, we find that, except for O+N pyrr , both N nit and N pyrr defects can be annealed by all of the above codoping elements (negative E rf ), because the formation energies of N nit and N pyrr are rather large to begin with (Table 1) and their relatively large vacancies can easily accommodate an additional dopant atom. However, incorporating an additional dopant atom into the energetically more favourable pyridinic NG is found to be rather difficult: It was determined that the V 1 -N py 3 and V 2 -N py 4 defects can be annealed only with a B, F, Al, Si, or P atom. It was determined that the energetically less favourable pyridinic NG with lower N concentrations can also accommodate these elements (see Supporting Information Figs. S2 and S3). The results, e.g., indicates that a N py to N gr conversion scheme solely based on the carbon source [24] will be limited in that di-or larger vacancies cannot be annealed. Although the F codoping of p-type NG might be energetically feasible, it is found that the F atom cannot structurally passivate the vacancy defects in the V 1 -N py 3 and V 2 -N py 4 cases and accordingly the electronic structures of F+V 1 -N py 3 and F+V 2 -N py 4 still maintain the p-type and bipolar characters, respectively (see Supporting Information Fig. S4). This leaves B, Al, Si, and P as the codoping candidates for the n-type conversion of p-type NG.
We first focus on the pyridinic NG codoped with B, Al, Si, and P, whose fully relaxed geometries and electronic structures are shown in Fig. 4. The C-N bond lengths of pyridinic NG are elongated when the vacancy is annealed by the codopant (Table 2). We further find that the trend of energetic stability of codoped NG (Fig. 3) is strongly related with the planarity of their geometries. For example, B can anneal the V 1 -N py 3 defect without protrusion (Fig. 4a), which makes its E rf larger than those of Al, Si, and P. Similarly, Al can anneal the V 2 -N py 4 defect in a fourfold-coordinated planar configuration (Fig. 4f), which makes it an energetically more favourable codopant than B, Si, and P.
In terms of electronic structure (Fig. 4 bottom panels and Table 2), we find that all the X+V 1 -N py 3 and X+V 2 -N py 4 cases (X = B/Al/Si/P) become n-type, or E D is now located below E F . The passivation of V 1 -N py 3 and V 2 -N py 4 vacancy defects by the codoping of B/Al/Si/P allows the charge transfer from the N atoms to the graphene π* states (Figs. 2e and 2f) as in the case of N gr (Fig. 2a). The alteration of the electronic type change is also easily identified in the COHP data, in which we observe strong antibonding C-N peaks that are pinned near E F or the downward shift of E D below E F . The reduction of work function for V 1 -N py 3 and V 2 -N py 4 was up to 0.62 eV and 0.52 eV upon codoping of B/Si and Al, respectively.
To further understand the microscopic mechanisms of electronic type conversion in X+V 1 -N py 3 andX+V 2 -N py 4 , we have decomposed the C-N COHPs into different orbital contributions (Fig. S5) and found that the C-N antibonding COHP peaks mostly come from the C 2p z -N antibonds. Namely, the insertion of a B/Al/Si/P atom into V 1 -N py 3 or V 2 -N py 4 and the subsequent annealing of vacancy allows the N atoms to directly couple with the sp 2 C network and behave like the graphitic Ns. Whereas the amount of downward shift of E D tends to decrease with the lowering of N concentrations, we have confirmed that the B, Al, Si, and P codoping of V 1 -N py 1∼2 and V 2 -N py 1∼3 in general changes their electronic type to n-or at least bipolar ones (see Supporting Information Figs. S6-S10).
Compared with the pyridinic cases, the nitrilic and pyrrolic NG show a different type conversion behavior and thus become interesting comparative systems ( Fig. 5 and Table 2). The C-N bond lengths in nitrilic and pyrrolic NG are more elongated than in the pyridinic counterparts. We also observe that, except the case of Si+N nit (Fig. 5c bottom panels), all of them show the bipolar property (Table 2). The different behaviour between the N py family (can be converted to the n-type) and N nit or N pyrr (can be converted only to the bipolartype) upon the introduction of a codoping element can be understood from their band structure plots. While the pyridinic defect states can energetically mix with the π states of graphene ( Figs. 1b and 1c bottom panels), the nitrilic and pyrrolic defects states are energetically located further away from E F and cannot strongly hybridise with the graphene π states (Figs. 1d and 1e bottom panels).
In spite of the difference, we can generally conclude that the p-type character of the structurally more distorted and energetically less probable N nit and N pyrr defects (Table 1) can be also eliminated by the B, Al, and P codoping. The COHPs of C-N in these cases are nearly zero around E D , or the N atoms in the codoped N nit and N pyrr do not affect the sp 2 carbon network. The nature of charge distribution in X+N nit and X+N pyrr can be again visualised in the Mulliken population analysis plot as shown for the representative B+N pyrr case in Fig. 2g, which shows that the net induced charge is localised near the defect sites as in the (bipolar) V 2 -N py 4 (Fig. 2c). It remains to explain why the Si-codoped N nit maintains the p-type character (Fig. 5c bottom panel). To understand its origin, we have additionally analysed the COHPs of various bonds (See Supporting Information  Fig. S11), and found that Si+N nit can be discriminated from the B/Al/P+N nit counterparts by a strong antibonding Si-C COHP peak that appears below E D and is pinned near E F (Fig. 5c bottom panel). On the other hand, e.g., the B-C in B+N nit has the bonding character around E F (Fig. 5a bottom panel).
We finally demonstrate a very desirable feature of the codoping approach in terms of charge transport characteristics. We show in Figs. 6a and 6b the transmission functions of (codoped) NG for the representative cases of (B+)V 1 -N py 1 and (P+)V 2 -N py 1 . Note that in these NEGF calculations, unlike the DFT results that have been discussed so far, we are considering isolated doping sites sandwiched by two semi-infinite pristine graphene electrodes (See Supporting Information Fig. S1). The work function of the entire junction system is thus determined by the pristine graphene region, which results in the Dirac point located at E F . The resulting currentbias voltage curves calculated according to the Landauer-Büttiker formula, where µ 1 = E F − 0.5eV and µ 2 = E F + 0.5eV , are shown together in Figs. 6c and 6d, respectively. For the untreated p-type NG models (V 1 -N py 1 and V 2 -N py 4 ), we find that the localised defect states around vacancy sites result in currents much decreased from that of pristine graphene. However, upon introducing codopant atoms and healing the vacancy defects (B+V 1 -N py 1 and P+V 2 -N py 1 ), we almost completely recover the currents of pristine graphene. We have recently obtained a similar behaviour for the B-N edge-doped graphene nanoribbons [46].
The availability of a simple yet effective method to convert NG with mixed n-and p-type characters into nearuniform n-type NG will improve the performance and reliability of various graphene-based electronic devices such as all-graphene p-n junctions, field-effect transistors, integrated circuits, etc. [47][48][49][50] We additionally note that the codoping-induced type change in NG may have significant implications for energy applications. [10][11][12][13][14][15][16][17][18][19][20] In the effort to identify the catalytically active sites in NG that result in significantly improved ORR performance, it was recently shown that the ORR activity in NG-based fuel cell cathodes is proportional to the graphitic N content. [17] On the other hand, it was experimentally demonstrated that the codoping of B or P into NG improves the ORR performance. [28,29] While the precise origin of the enhanced ORR activity in NG and codoped NG is not clear yet, we can propose that the codopinginduced conversion of mixed-type NG into uniformly ntype NG (i.e. generating more "effectively" graphitic N atoms) and the enhanced charge transport capacity will play an important role in enhancing the ORR activity.

IV. CONCLUSION
In summary, based on fully first-principles computations, we have predicted that the post-synthetic codoping of p-type NG with B, Al, and P can anneal even relatively large vacancy defects associated with pyridinic, pyrrolic, and nitrilic Ns. We also found that, while the pyridinic NG can be converted into n-type, only the bipolar type conversion is allowed for the nitrilic and pyrrolic counterparts. More importantly, employing novel analysis schemes, we have revealed the origin of p-type or bipolar behaviour of NG, and showed that the codoping and consequent healing of vacancy defects strongly modify the C-N antibonding character, its energetic position, and the spatial distribution of electrons donated from N atoms. We expect this will provide in the future a framework to systematically characterise doped and/or defective graphene. Furthermore, we demonstrated that the codoping of p-type NG can recover the excellent charge transport capacity of pristine graphene, and suggested that together with the electronic type change it should have significant implications for energy and electronic device applications. In view of experimental advances already achieved in the control of doping characters of NG, [15,17] we envisage the facile post-synthetic codoping scheme proposed here will lead to new pathways for tailoring and enhancing the properties of graphene at the atomic level.
Note added in proof : After submission of this work, we became aware of Ref. 51, in which the authors experimentally synthesised P, N-codoped graphene by chemical vapour deposition method and found much improved airstable n-type characteristics. In addition, in Ref. 52, the authors theoretically predicted that the dissociative adsorption of H 2 molecule on the trimerized and tetramerized pyridine-type defects is energetically favourable, and the adsorption of the two H atoms change their electronic properties from p-type to n-type doping.