Ab initio study of the processes of nitrogen functionalisation in graphene

Nitrogen functionalisation of graphene is studied with the help of ab initio electronic structure methods. Both static formation energies and energy barriers obtained from nudged elastic band calculations are considered. If carbon defects are present in the graphene structure, low energy barriers on the order of 0.5 eV were obtained to incorporate nitrogen atoms inside the sheet. For defect-free graphene, much larger barriers in the range of 3.70–4.38 eV were found, suggesting an external energy source is required to complete this type of incorporation.


Introduction
Graphene is a two-dimensional material composed of a single layer of carbon atoms in a honeycomb lattice.Since its first synthesis/isolation [1], in 2004, a considerable amount of research activities has been devoted to its study due to its remarkable properties [2][3][4][5][6][7].The sp 2 hybridisation of its carbon orbitals means graphene is a robust yet flexible material with a rather high thermal conductivity [3].Its atomic thickness also makes it nearly transparent to light [4].
This article presents a theoretical study of nitrogen doping of graphene performed within the framework of the density functional theory (DFT) [27].The methodologies of the two main computational strategies used, namely the examination of the relative stability of various defects and nitrogen configurations via formation energy (FE) calculations and the investigation of various reaction paths using the nudged elastic band (NEB) [28] method, are presented in section 2. The results of the FE calculations are presented and analysed in section 3 in order to obtain insights into various doping processes.
In the context of post-growth functionalisation, plasma treatments [29,30] are advantageous because of the presence of dissociated nitrogen atoms that are much more reactive than bonded N 2 molecules.With a correctly tuned treatment, other plasma particles can also provide surface activation, which helps overcome reaction barriers, without damaging significantly the material.NEB computations were used to investigate the energy barriers that a nitrogen atom must overcome to be incorporated in the graphene sheet, both with Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. and without a pre-existing defect.These results are presented in section 4, along with an assessment of the qualitative model used to explain the functionalisation process.Final observations regarding the present study are available in section 5.

Computational methods
Spin-polarised DFT calculations were performed, as implemented in the open-source software BigDFT [31,32].The exchange-correlation energy was evaluated using the generalised gradient approximation (GGA) [33].The core electrons were replaced by relativistic separable dual-space Gaussian pseudopotentials [34], including non-linear core corrections [35].Due to the semimetallic nature of graphene, the Pulay potential mixing procedure [36] was used for the last few occupied orbitals to reach convergence.The other parameters are further explained below.

Formation energy calculations
Formation energy (FE) calculations were used to study the stability of different defects and doped configurations.They aim to compare structures with varying compositions relative to a common reference point.The general formula is where E for is the formation energy of the defective structure, E tot is its total DFT energy and the sum runs over the chemical elements present in the structure and N i and μ i stand respectively for the number of atoms and the chemical potential associated with each element.
In the present study, the structures contain only carbon and nitrogen atoms.As in many previous works on the subject [37][38][39][40][41][42][43], pristine graphene and the N 2 molecule are used as reference points, meaning the chemical potentials in equation (1) are equal to 1/N S and 1/2 of the total energy of these structures, with N S being the number of atoms in the considered supercell.
Ideally, the FE would be obtained for an isolated defect in an otherwise perfect graphene structure.However, since periodic boundary conditions are used, there is a residual interaction between defects of neighbouring cells as the range of the strain field can be quite extended.To assess the importance of this interaction, two different supercell setups were used for the calculations.
The BigDFT code is limited to orthorhombic simulation cells.For this reason, the orthorhombic primitive cell of graphene containing four atoms was used as a starting point.The supercells used were constructed as 6 × 4 and 9 × 5 replicas of this minimal four atoms cell.This amounts to respectively 96 and 180 atoms of pristine graphene.These combinations were chosen for two reasons: they create supercells that are close to being square, and they also ensure that the special points K of the primitive Brillouin Zone, where the Fermi energy is located, are mapped back to Γ.
More details on this last point are available in the supplemental material (SM) [44].
The BigDFT input parameters were converged to obtain a 50 meV precision on the FE.Non-shifted two-dimensional grids were used to sample the momentum space of sizes 4 × 4 for the 96 atoms supercell and 3 × 3 for the 180 atoms supercell.
The cell lattice parameters were determined by minimising the energy of a pristine graphene cell.The cell dimensions were then fixed for all subsequent calculations in order to minimise the impact of defects on the graphene structure of neighbouring cells.For each configuration, internal forces were relaxed to less than 5 × 10 −4 Ha Å −1 , using the stabilised quasi-Newton optimisation algorithm [45].This small residual force value allows for adequate precision on the ionic positions.

Reaction paths calculations
Reaction paths were computed with the NEB module as implemented in BigDFT [46].The Fletcher-Reeves algorithm [47] was used to find saddle points of the potential energy hypersurface.Final converged calculations used the 96-atom supercell and strings of at least 17 images with a 2 × 2 k-point sampling grid.Residual perpendicular forces on images were relaxed until they were smaller than 5 × 10 −3 Ha Å −1 or, in more difficult cases, until the energy barrier value was sufficiently stable between iterations.The Climbing Image (CI-NEB) algorithm [48] was used to ensure one of the images converged to the saddle point and thus corresponds directly to the energy barrier.This algorithm modifies the applied forces on the highest energy image to push it closer to the saddle point.

Studied structures
The most common defects found in pure and nitrogen-doped graphene were studied with the help of 96 and 180-atom supercells.Due to the periodic nature of the simulation cells, a large number of atoms is necessary to reduce the interaction between neighbouring strain fields induced by the defect, which has been shown to impact FE values [49] and properties [50] of impurities.
The orthorhombic nature of the simulation cells means that most studied defects had two nonsymmetrically equivalent orientations possible, as shown in the SM [44].The results presented in table 1 correspond to the lowest energy of these possible configurations.Figure 1 presents the studied planar structures in which no atom is found outside of the plane.Structures with defects containing only carbon atoms are represented in figures 1(a)-(c).The monovacancy is created simply by removing one carbon atom from the pristine mesh of graphene.The divacancy defect is obtained by removing a second carbon atom adjacent to the first one.A partial reconstruction of the dangling bonds around the defects is observed for both structures, as it is favourable for the carbon atoms to keep their sp 2 hybridisation.The third carbon-only defect (figure 1(c)) is the Stone-Wales (SW) defect, in which a carbon-carbon bond rotates 90°from its original orientation.
The next structures correspond to various nitrogen functionalisation configurations.They are presented in figures 1(d)-(f).The first one is known as graphitic or substitutional doping, where a nitrogen atom replaces a carbon atom, creating three bonds with its nearest neighbours and keeping the graphene structure mostly intact.Graphitic doping was found to induce n-type doping [37], or an increase of the Fermi energy.The second type of nitrogen functionalisation is the pyridinic structure.In this case, the nitrogen atom is paired with a vacancy on a neighbouring site and only creates two bonds with neighbouring carbon atoms.It nevertheless stays in a hexagonal structure.
The last planar structure, visible in figure 1(f), is the pyrrolic functionalisation.The nitrogen atom also forms two bonds with neighbouring carbon atoms.The difference is that it sits in a pentagonal cycle rather than a hexagonal.This type of functionalisation is usually found around large vacancies or  grain boundaries, which makes it harder to study in a periodic simulation setup.In figure 1(f), only the nitrogen on the bottom is of pyrrolic type.Without the other two, it would stabilize in pyridinic form as it is more stable in small vacancies.Therefore, the reported FE should not be interpreted directly as the pyrrolic FE but can be compared to previous results from the literature.Both pyridinic and pyrrolic functionalisation are associated with p-type doping [51,52], or a decrease of the Fermi energy.
There is also the possibility of non-planar adsorbed atoms on graphene.For both carbon and nitrogen adatoms, the ground state is found to be around 1.7 Å above a carboncarbon link.This bridge configuration is represented in figure 2. Those states have a role to play in certain experiments and technological applications, for example, graphene growth by chemical vapour deposition (CVD) [53,54].The model proposed in section 4 also expects that adsorbed nitrogen atoms play a key role in post-growth functionalisation treatments.

Formation energy results
A complete list of FE values for all the above-mentioned configurations is presented in table 1 along with similar results from the literature.Energy discrepancies smaller than 0.3 eV are found in most cases.They can be mainly explained by cell sizes and choices of exchange-correlation functionals and pseudopotentials.Another source of discrepancies is the inclusion or omission of the special point K in the momentum space sampling.If not included in the sampling, the FE can be overestimated, as discussed in the SM [44].
Some relevant insight can be obtained by analysing table 1.First, our numerical results are consistent when comparing defects across both supercell sizes.It suggests that defect-defect interaction is acceptably low with the 96-atom cell, and it can be considered as converged in terms of cell size.That being said, the FE of larger defects, such as the divacancy and SW, show larger differences between 96-and 180-atom cells.This is mostly due to the larger strain field they generate in the graphene structure, which can be divided by more atoms in larger supercells.
Another interesting observation is that a divacancy's energy is lower than that of a single vacancy.This can be explained by the partial reconstruction of the bonds when one or two carbon atoms are removed.With two vacancies, all four adjacent carbon atoms have a third neighbour to bond with after reconstruction, meaning that they still are in a favourable sp 2 environment (see figure 1(b)).This is not possible for a monovacancy: two of the three carbon atoms close to the defect do tend to form a pentagon in order to keep the sp 2 environment, but the third one has to remain isolated [49], as shown in figure 1(a).As a consequence, a high defect concentration should increase the presence of double vacancies rather than isolated ones.Similar effects have been both predicted [43,55] and experimentally observed [56] regarding the favoured presence of complex defects rather than a high concentration of simple defects.
Notably, the FE of a SW defect is relatively high for a configuration where no carbon atoms are missing.Previously calculated energy barriers show that the formation of this defect is much harder than the inverse reaction, with respectively 9.2 eV and 4.4 eV barriers [57].Experimental results also seem to support this conclusion [58].
Since all FEs presented in table 1 are positive values, one could think the functionalised structures are unfavourable.However, functionalisation methods such as plasma-assisted doping [29,30] involve atomic nitrogen with an energy of about 5.2 eV compared to a nitrogen molecule.This means, for instance, that an isolated nitrogen would lower its energy by around 0.8 eV upon adsorption on a graphene surface when compared to the value for the bridge N configuration in table 1.This configuration is thus energetically favourable in such conditions.The energy of a graphitic nitrogen atom is also much lower than the combined energies of an adsorbed nitrogen and a vacancy.This suggests that graphitic doping is highly favourable when a nitrogen atom is adjacent to a vacancy.Upon completion of the process, the dopant is stable in the material.The same conclusion holds for pyridinic and a divacancy defect.
As raised earlier, the pyrrolic FE should not be analysed directly.It is nevertheless interesting to observe such a small value for a structure combining three heteroatoms and a vacancy.This hints at the existence of more complex defective structures around grain boundaries and vacancy clusters, especially for a higher concentration of dopants.This phenomenon has been observed experimentally [59].
In light of these results, the 96-atom cell was deemed sufficient to analyse small graphene defects and was used in the next section for reaction path analysis.

Reaction paths results
The present section explores possible transitions between atomic nitrogen in the gas phase and functionalised graphene.This situation is mostly applicable to post-growth functionalisation methods where the dissociation of the nitrogen molecule is already assured.It has been found that the type and concentration of nitrogen defects vary significantly depending on experimental setups and parameters [60], confirming that understanding the dynamics of the process is important to have more control over the final product.A simplified qualitative model will be presented, explaining the complete reaction path shown in figure 3.This model helps obtain a deeper understanding of the complete process leading to nitrogen functionalisation of graphene.It is based on the computation of energetic barriers between each intermediate steps, done with the CI-NEB algorithm.A large range of different scenarios can be explored in this fashion in order to study the functionalisation process from start to finish.Horizontal axis of figures in this section correspond to total atomic displacements.These values are obtained by summing the small displacements between each image, for the atom with the largest total displacement during the reaction.

Nitrogen adsorption
The starting configuration of the model corresponds to an atomic nitrogen far enough from the graphene so that they do not interact with each other.Inside a plasma used for experimental treatment, a significant proportion of dissociated nitrogen is created by way of an electric discharge.It is assumed that they are, by far, the major source of dopants due to the high chemical stability of N 2 molecules.In section 3, the graphene sheet is chosen as the reference energy.An isolated atomic nitrogen in the plasma then has a 5.23 eV energy, which corresponds to the DFT value for μ N from equation (1).When it is adsorbed to a bridge configuration, the energy of the system is lowered to 4.40 eV, as reported in table 1.This first step in the qualitative model is therefore not a limiting factor for the nitrogen functionalisation of graphene.

Nitrogen diffusion
Adsorbed atoms can diffuse from one bridge site to an adjacent one.By repeating this process, an atom can move around on a surface and may eventually reach a more reactive site for the functionalisation reaction.Energy barriers have already been calculated for carbon and nitrogen diffusion on graphene, using GGA.Previously reported values for carbon diffusion are 0.43 eV [10] and 0.55 eV [42] for calculations using the PBE functional and 0.47 eV [61] for the PW91 functional.For the nitrogen diffusion, they are 0.86 eV [10,42](PBE).Complete results from our calculations are presented in figure 4. They yield energy barriers of 0.46 eV and 0.95 eV for these reactions, in good agreement with the literature.
The mobility of carbon and nitrogen atoms on graphene can then be estimated, using the Arrhenius equation, defined as where k is the rate of reaction, σ is the number of symmetrically equivalent reactions (4 for this diffusion), A is a preexponential factor, E a is the calculated energy barrier and k B T is the thermal energy.The pre-exponential factor can be calculated precisely, but as a first-order approximation, it can be set to the highest phonon frequency [62], which has been reported as 1586 cm −1 [63].At room temperature, the order of magnitude of k is therefore 10 6 Hz for carbon and 0.1 Hz for nitrogen.The higher mobility of carbon on graphene is well-known and sufficient for CVD [54].As a comparison, very mobile atoms, such as fluorine have been reported at 10 11 Hz at the same temperature [64].Even though the reaction rate is much lower in the case of the diffusion of nitrogen atoms, the area covered can still be significant for an experimental treatment that lasts a few minutes.Assuming a sufficient exposure to nitrogen, reactive sites will be reached by adatoms.

Main functionalisation reaction
As shown earlier, an adsorbed nitrogen atom can move just above the graphene surface to potentially reach a reaction site.However, to obtain one of the functionalised state presented  in section 3, an incorporation mechanism is needed, for which an energy barrier can be calculated.Various incorporation mechanisms were explored in the present study, both with and without the presence of native defects in the graphene sheet.They are presented and analysed in the next two subsections.
4.3.1.Incorporation with defects.The presence of defects in graphene should make the incorporation of nitrogen atoms into the sheet much easier, as the presence of dangling bonds increases reactivity.Furthermore, higher FEs for a starting configuration compared to the final one usually means an easier reaction to complete.Precise activation energies still need to be calculated in order to verify these assumptions.Many defects and incorporation mechanisms could have been considered, but only two case studies are presented here.The first one concerns the transition from a monovacancy to graphitic doping, that is from figures 1(a)-(d), while the second focuses on the transition from a divacancy to pyridinic doping, i.e. from figures 1(b)-(e).As these mechanisms do not require the evacuation of an extra carbon atom after the incorporation, the last reaction step in figure 3 can be ignored.
The starting configuration for the first transition is an adsorbed nitrogen above one of the six carbon-carbon bonds around a single vacancy.The presence of a nitrogen atom will naturally move the partially reconstructed bond opposite to the nitrogen, as can be seen in the first inset of figure 5(a), which means all starting positions are symmetrically equivalent.It is interesting to note this structure is 3.40 eV lower in energy than the system with an isolated vacancy and an isolated adsorbed nitrogen atom, indicating that the nitrogen atom is saturating some dangling bonds.The second inset shows a metastable state where the nitrogen is bonded to only one carbon.This intermediate state raises the energy barrier compared to the initial state, because it is lower in energy.Still, the energy curve reported in figure 5(a) exhibits a low barrier value of 0.55 eV.This following reaction releases 7.50 eV when the proper incorporation happens, due to the low FE of the graphitic doping.This reaction is thus highly favourable and should happen quickly once its starting conditions are met.
The second case-study corresponds to the nitrogen incorporation inside a divacancy.The initial configuration is the one where the nitrogen atom is adsorbed above a carboncarbon bond, as shown in the first inset of figure 5(b).There is a second non equivalent bond neighbouring the incorporation site, but we found that a nitrogen atom adsorbed there will more favourably diffuse to the shown configuration than incorporate.As discussed earlier, divacancies are more stable than single vacancies and this structure is only 0.92 eV lower in energy than the isolated defect and the isolated adsorbed nitrogen bridge.
After a diffusion towards the defect, the nitrogen can be incorporated.It then creates pyridinic doping, such as in figure 1(e).The curve in figure 5(b) looks rather similar to the one in figure 5(a), as it exhibits a small energy barrier followed by a significant energy decrease.The value of that barrier is 0.46 eV, meaning that this process is also expected to be favourable and easily completed.
Incorporation in these types of defects has been observed qualitatively using empirical potentials [65].Other mechanisms could be studied, with larger defects or a combination of small ones.The above results nevertheless indicate that the reactivity of defects ease the restoration of the original physical structure.It is worth noting that the energy barriers for the incorporation into those two defects are lower than the barrier for nitrogen diffusion shown in figure 4, suggesting that defects exercise a kind of attraction on neighbouring nitrogen adatoms.with an adsorbed nitrogen atom as in figure 2 and end with a nitrogen atom in a graphitic configuration, as in figure 1(d).
The first mechanism serves as a reference point for the next ones.It simply consists of the nitrogen atom in a bridge configuration switching places with one of the two neighbouring carbon atoms beneath it.The nitrogen atom is then in a graphitic position and the carbon atom is in a bridge configuration, adsorbed on the surface.Following this first reaction, the carbon atom can diffuse away so that a pure graphitic doping is obtained.The obtained energy curve and relaxed trajectory for this mechanism are presented in figure 6(a).The incorporation reaction has a barrier value of 4.38 eV, which is much higher than the values reported earlier in the presence of defects.To be completed, this mechanism needs to disturb the graphene structure and temporarily break some bonds.Finally, the ejected carbon atom can do a first diffusion away from the nitrogen, with a 0.71 eV barrier.It is interesting to note that this final structure, an adsorbed carbon one bond away from the incorporated nitrogen, is 0.82 eV more favourable than having those two defects isolated.Therefore, the ejected carbons in all reactions presented in the current section would be more likely to stay near the incorporation site.
Other mechanisms were explored in order to find more favourable ways to incorporate nitrogen without native defects.It was predicted, precedently to the discovery of graphene, that the presence of a carbon adatom could catalyse reactions in carbon materials [66].More recently, a similar effect was demonstrated for the formation and healing of SW defects in graphene [67].A similar mechanism is studied here, with an adsorbed nitrogen instead of carbon.While forming the SW defect, the nitrogen is incorporated in-plane and rotates with a neighbouring carbon inside the sheet.Meanwhile, another carbon atom (intially two bonds away from the nitrogen) moves to an adsorbed position.The energy curve associated with this reaction is presented in figure 6(b).The incorporation mechanism is made of two consecutive energy barriers, of respectively 3.98 and 2.11 eV.They roughly correspond to the bond rotation and the expulsion of the second carbon atom.This reaction is followed by the rotation back in place of the nitrogen-carbon bond, limited by a barrier of 3.34 eV.The maximal and leading barrier for this reaction is thus lower than the 4.25 eV barrier of the first mechanism, but some extra steps are required to complete the reaction.Thus it is not clear which mechanism is the most favourable.
The last mechanism is somewhat similar to the first one: instead of exchanging places with a carbon atom, the adsorbed nitrogen pushes it to the other side of the sheet and takes its place inside the graphene lattice.The energy curve related to this process can be seen in figure 6(c).This mechanism is the most favourable one to incorporate graphene without defects as it features a 3.70 eV barrier, the lowest value reported here.The following carbon diffusion, to isolate the graphitic configuration, is still associated to a 0.71 eV barrier.One must be cautious with this last result, given that the experiments are often performed with a graphene sheet on top of a substrate.This must affect the electronic environment under the sheet, and could affect this mechanism more than the others.Nevertheless, the presence of intercalated atoms between graphene and its substrate has been experimentally observed after nitrogen doping [68], which suggests that this type of reaction might play an important role in nitrogen incorporation.
All these results show that the incorporation of nitrogen atoms without defects in graphene is much less favourable than in the presence of defects.Frequencies of such reactions, according to equation (2), are extremely low, ranging from 10 −48 to 10 −57 Hz at room temperature.Such reactions would therefore need an external source of energy activating the surface and helping overcome the barriers.In the context of a nitrogen plasma treatment, a significant proportion of excited metastable N 2 molecules [69] is observed.Each molecule carries around 6 eV of vibrational energy [70], and could provide the required activation energy during the plasma treatment.The significant presence of dopants after a few minutes [30] could then be explained, even for graphene sheets with a very low initial concentration of defects.

Conclusion
Nitrogen functionalisation is a promising way to tailor graphene for technological applications.This work gave some insights on the functionalisation process at the atomic level, both with FE and NEB calculations.Energy values for various defects were obtained and compared to previous results from the literature, mostly agreeing with them.It was shown that a 96-atom supercell is sufficient to simulate most of the strain relaxation around localized defects.A qualitative model to explain the whole functionalisation process was presented.To quantify this model, NEB calculations were conducted.They showed that nitrogen incorporation inside native defects is a favourable process even at room temperature, but that incorporation without defects is harder and would require an external source of energy.

Figure 1 .
Figure 1.Planar structures for which the FE was computed.Results are available in table 1.

Figure 2 .
Figure 2. Adsorbed nitrogen bridge configuration on graphene.Left side is a vertical projection and right side corresponds to a 60°p rojection.Carbon adatoms occupy approximately the same position.

Figure 3 .
Figure 3. Qualitative model to explain the complete functionalisation reaction.The initial state is a gaseous atomic nitrogen far from a graphene sheet and the final state is a planar functionalised state.Each step of the reaction is limited by a transition state which determines a corresponding energy barrier, represented by red lines.Figure 4. Calculated energy curves for carbon (black) and nitrogen (green) diffusion on graphene from a bridge site to a neighbouring one.Both physical paths are similar, although the nitrogen follows a more curved trajectory.

Figure 4 .
Figure 3. Qualitative model to explain the complete functionalisation reaction.The initial state is a gaseous atomic nitrogen far from a graphene sheet and the final state is a planar functionalised state.Each step of the reaction is limited by a transition state which determines a corresponding energy barrier, represented by red lines.Figure 4. Calculated energy curves for carbon (black) and nitrogen (green) diffusion on graphene from a bridge site to a neighbouring one.Both physical paths are similar, although the nitrogen follows a more curved trajectory.

4. 3 . 2 .
Figure Energy curve for the two mechanisms inside native defects.

Figure
Figure Energy curve for the three incorporation mechanisms without native defects.

Table 1 .
FE values of various configurations, in eV.Results computed during this study are in the top two rows and are compared to results from the literature in the bottom six rows.