Ordered Semiconducting Nitrogen-Graphene Alloys

The interaction between substitutional nitrogen atoms in graphene is studied by performing first principles calculations. The nearest neighbor interaction between nitrogen dopants is highly repulsive because of the strong electrostatic repulsion between nitrogen atoms, which prevents the full phase separation in nitrogen doped graphene. Interestingly, there are two relatively stable nitrogen-nitrogen pairs due to the anisotropy charge redistribution induced by nitrogen doping. We reveal two stable semiconducting ordered N doped graphene structures C3N and C12N through the cluster expansion technique and particle swarm optimization method. In particular, C12N has a direct band gap of 0.98 eV. The heterojunctions between C12N and graphene nanoribbons might be promising organic solar cells.


Introduction
Graphene, a single layer of carbon atoms arranged in a honeycomb lattice, has been the focus of recent research efforts [1][2][3], due to its unique zero-gap electronic structure and the massless Dirac fermion behavior. The unusual electronic and structural properties make graphene a promising material for the next generation of faster and smaller electronic devices.
However, graphene lacks an essential feature for controlled and reliable transistor operation, namely, a band gap around the Fermi level. Several schemes have been proposed to open a band gap in graphene. One is the substrate-induced gap for graphene supported on SiC [4], but the experimental realization of this idea turned out to be very difficult and controversial [5]. Another approach is the creation of gaps through confinement [6], such as in narrow graphene nanoribbons (GNRs). However, a large scale production of narrow GNRs is still very challenging. In addition, it was demonstrated that the band gap of a graphene bilayer can be controlled externally by applying a gate bias [7,8]. Unfortunately, bandgap tailoring by external electrostatic gate has limited tunability for the transistor on-off ratio. Chemical functionalization [9,10] is an alternative way to manipulate the electronic properties of graphene. The hydrogenation of graphene [9,[11][12][13][14][15][16][17], as a prototype of covalent chemical functionalization, was studied extensively. However, the fully hydrogenated graphene, i.e., graphane CH, is an insulator with a very large band gap (larger than 5 eV) [18], which prevents from using graphane in a field-effect transistor. There were some efforts to tune the band gap of two-dimensional (2D) graphene by varying the hydrogen concentration.
Unfortunately, the phase separation between bare graphene and graphane will take place spontaneously, leading to an essential zero band gap for the whole system [12]. Similar phenomenon was found to also occur in the oxidized graphene [19].
Previous studies [20][21][22][23][24][25] revealed that doping graphene and related materials with nitrogen is effective to tailor its electronic property and chemical reactivity because of the stronger electronegativity of nitrogen compared to that of carbon and conjugation between the nitrogen lone pair electrons and the graphene π-system. This could create novel nanomaterials and expand its already widely explored potential applications. Experiments [26] suggest that nitrogen species have been incorporated into the graphene structure with content in the range of 4.5-16.4%. Usually, the nitrogen atoms are substitutionally incorporated into the basal plane of graphene in the forms of pyridinic, pyrolic and graphitic nitrogen 2 bonding configurations. X-ray photoelectron spectroscopy showed that the pyridinic and pyrolic nitrogen may lie at the edge or the defect sites and the graphitic nitrogen (the nitrogen replacing carbon in graphene plane) is dominant [27]. A recent scanning tunneling microscopy study also showed that individual nitrogen atoms were incorporated as graphitic dopants [28]. Theoretically, the isolated substitutional nitrogen doping in graphene was studied before [28]. However, the nature of the interaction between substitutional nitrogen atoms is far from clear. How the electronic structures of graphene is modified by multiple nitrogen dopants is an interesting open question to be addressed.
In this paper, we examine the collective behavior of nitrogen doping in graphene on the basis of density functional calculations. We find that the nearest neighbor interaction between nitrogen dopants is highly repulsive because of the strong electrostatic repulsion between negatively charged nitrogen atoms. Our calculations reveal two stable ordered N doped graphene structures C 3 N and C 12 N. Surprisingly, nitrogen atoms do not simply doping electrons to the graphene π * state. Instead, both ordered structures are found to be semiconducting. In particular, C 12 N has a direct band gap of 0.98 eV. The substitutional nitrogen raises the valence band maximum (VBM) and conduction band minimum (CBM) level in N doped graphene, leading to a type-II band alignment between C 12 N and 6-AGNR.
We propose that C 12 N/6-AGNR is a promising organic solar cell with C 12 N as donor and 6-AGNR as acceptor, respectively.

Computational Method
In our DFT calculations, the local density approximation (LDA) was adopted unless otherwise stated. The plane-wave cutoff energy for wavefunction was set to 500 eV. The ionelectron interaction was treated using the projector augmented wave (PAW) [29] technique as implemented in the Vienna ab initio simulation package [30]. For the Brillouin zone sampling, we use 2×2×1 Monkhorst-Pack k-mesh for the 12×12×1 graphene supercell. And for C 3 N (C 12 N), we used the 12 × 12 × 1 (6 × 6 × 1 ) Monkhorst-Pack k-mesh. The constant temperature first-principles molecular dynamics simulations were performed to check the thermal stability of C 3 N and C 12 N. The "cluster expansion" [31] of the alloy Hamiltonian was carried out by the ATAT package [32]. The global minimum optimization was performed using the two-dimensional (2D) particle swarm optimization (PSO) technique, which we proposed recently to predict the most stable 2D crystals [33]. For the PSO simulations, we use the Crystal structure AnaLYsis by Particle Swarm Optimization (CALYPSO) code [34].

Interaction Between Nitrogen Dopants
To investigate the interaction between two substitutional nitrogen atoms, we use a 12×12 supercell of graphene. The use of such a large supercell is essential due to the relative long range interaction between nitrogen dopants. And test calculation shows that using larger supercells gives similar results. As is well known, the honeycomb lattice of graphene consists of two interpenetrating triangular sublattices, namely sublattice A and sublattice B. We consider the interaction between a nitrogen atom at 0A site and those at other sites [see Fig. 1(a)]: All interaction with the interaction range less than 6.5Å are included, and we also include several other longer range interaction N-N pairs (0A-11A, 0A-12B, 0A-13B). For each N-N pair, the internal atomic structure is fully relaxed. With respect to the energy of longest N-N pairs 0A-13B considered in this work (15.60Å), We plot the interaction energy of other N-N pairs in Fig. 1(b). First, we can see that all the interaction energy are positive, i.e., the interactions are repulsive. This can be understood because nitrogen has a high electronegativity than carbon thus nitrogen atoms are surrounded by more electrons which leads to the repulsive electrostatic interaction. Second, the nearest neighbor interaction is extremely repulsive. This is because the two nitrogen atoms are directly connected along the bond direction leading to a very strong electrostatic repulsion.
Third, generally the interaction energy decreases with the pair distance because of the nature of coulombic electrostatic interaction. Fourth, for N-N pairs with short pair distances (less than 6.5Å), the interaction energy curve has two local minima at N-N pairs 0A-3B and 0A-7B. The interaction between 0A and 3B and that between 0A and 7B are only slightly repulsive with the interaction energies of 0.08 eV and 0.03 eV, respectively.
The nontrivial non-monotonic dependence of the pair interaction energy on the N-N distance and the presence of local minima at N-N pairs 0A-3B and 0A-7B are interesting findings to be understood. Our calculations show that the strain relaxation is not responsible for the stability of the 0A-3B and 0A-7B N-N pairs because the energies of the unrelaxed pair structures display the same trend. Instead, we find that the stability of the 0A-3B and 0A-7B N-N pairs is due to the anisotropic electron charge density redistribution induced by the nitrogen substitution. Fig. 1(b) shows the charge difference between the isolated N-doped graphene and the undoped graphene. We can see that the charge redistribution in graphene is mainly along the nearest neighbor bonding direction. Since sites 3B and 7B are nearly opposite to the bonding direction, there is much less charge redistribution around site 3B (7B) than that around 2A and 4B (6A and 8B). Therefore, in the case of the 0A-3B and 0A-7B N-N pairs, the coulombic repulsion will be weaker than that between the other neighboring pairs.
Our calculations show that the interaction between substitutional nitrogen atoms in graphene differs fundamentally from that between adsorbed hydrogen (oxygen) atoms on graphene. In the hydrogenated (oxidized) graphene case, hydrogen (oxygen) tends to be close to each other in order to lower the kenetic energy of the carbon π electrons, leading to the phase separation between fully hydrogenated (oxidized) graphene and bare graphene [12,19].
The repulsive interaction between substitutional nitrogen atoms will prevent the phase separation of nitrogen doped graphene into undoped graphene part and highly nitrogen doped graphene part, in contrast to the case of the hydrogenation and oxidization of graphene [12,19].

Ordered Semiconducting Nitrogen-Graphene Alloys
In our above discussions, 0A-3B and 0A-7B are two relatively stable configurations of N-N pairs. For each type of two pairs, we can construct a graphene superstructure with a uniform distribution of nitrogen dopants. For the 0A-3B N-N pair, the ordered nitrogen doped graphene structure is a 2 × 2 superstructure (C 3 N) of the graphene, as shown in Fig. 2(a). We can see that the doped nitrogen atoms themselves form a honeycomb lattice, similar to graphene but with a doubled lattice constant. As a consequence, the carbon atoms are separated into isolated six-membered rings, as shown in the right panel of Fig. 2(a). For the 0A-7B N-N pair, the ordered nitrogen doped graphene structure is a √ 13 × √ 13-R46.1 • superstructure (C 12 N) of the graphene, as shown in Fig. 2(b). In this superstructure, the carbon atoms are connected by carbon-carbon bonds to form a single carbon domain, in contrast to the C 3 N case.
To examine the stability of the C 3 N and C 12 N superstructures, we use "cluster expansion" method [31] and direct comparison of total energies between different superstructures. In the "cluster expansion" method, the alloy Hamiltonian is mapped onto a generalized Ising Hamiltonian. In brief, for some nitrogen doped graphene configurations, we perform spin polarized DFT calculations to relax the cell and internal atomic coordinates. The energies of the relaxed structures are used to extract the interaction parameters of the alloy Hamiltonian. After obtaining the interaction parameters (See Supplemental Material [35]) of the alloy Hamiltonian, we find that the C 3 N superstructure shown in Fig. 2(a) is the most stable structure for the C x N alloy with the 25% nitrogen concentration (at least for cells with no more than 32 atoms). For the C 12 N superstructure with the 0A-7B N-N pairs, the relatively long range (5.11Å) pair interaction is not included in the alloy Hamiltonian. Therefore, we generate all superstructures with no more than 26 atoms and the same nitrogen concentration (7.7%) as C 12 N using the linear scaling algorithm recently proposed by Hart and Forcade [36]. Our calculations show that the C 12 N superstructure shown in Fig. 2(b) has the lowest energy among all considered structures. The lattice dynamics calculation shows that both C 3 N and C 12 N have no imaginary frequency phonon and thus are stable [35].
First-principles molecular dynamic simulations up to 500 K were performed to check the stability of C 3 N and C 12 N. We find that the N atoms only vibrate around the equilibrium position, which confirms the thermal stability of C 3 N and C 12 N.
Above we proposed two possible stable ordered nitrogen-carbon alloys (C 3 N and C 12 N shown in Fig. 2) by considering only the graphene-based structures. Here, we perform PSO global minimum optimization to confirm the proposed structures. In the 2D PSO simulations [33], the graphene lattice structure is not assumed. Instead, we generate random structures (both atomic positions and cell parameters) to initialize the simulations. The initial structures are subsequently relaxed before performing the PSO operations. The population size is set to 30. We consider all possible cell sizes with the total number of atoms no more than 26. The number of generations is fixed to 30. Our simulations show that the most stable 2D nitrogen-carbon alloy structure with 25% N concentration is indeed the structure (p6mm plane group) shown in Fig. 2(a) [see also Fig. 3(a)]. The second lowest energy 2D structure with the plane group p2mg [ Fig. 3(b)] has a higher energy by 0.07 eV/N. It should be noted that both p2mg and p6mm structures can be viewed as nitrogen doped graphene structures although the graphene honeycomb lattice is not assumed in the PSO simulations. 6 The PSO simulations for the 2D nitrogen-carbon alloy structure with 7.7% nitrogen concentration do not give structures more stable than the C 12 N shown in Fig. 2(b), in support of the stability of the proposed C 12 N structure. We note that the formation energies of C 3 N and C 12 N are positive if graphene and N 2 molecule are taken as the references. However, many experiments [20,28] showed that N atoms can be incorporated into graphene. Here we show that C 3 N and C 12 N are the most stable structures once N atoms are incorporated into graphene, as can be seen from the convex hull plot in the new supporting material [35].
The band structures for the C 3 N and C 12 N superstructures from the LDA calculations are shown in Fig. 4(a)  And C 3 N has one unique Clar formula, thus belongs to the class 1CF of the (pseudo)-allbenzenoid structure. In contrast, C 12 N has no Clar sextets, and thus are nonbenzenoid.
Because each nitrogen atom forms three single bonds with the neighboring carbon atoms, the nitrogen atoms are excluded when constructing the Clar formula. We note that in the Clar formula for the -NH terminated graphene nanoribbon, the two N-C bonds were also treated as single bonds [38].  [39] showed that all graphene related structures has a relative large band gap except for the (pseudo)-all-benzenoid structure with class nCF (n >= 3) such as graphene which has small or zero band gap. Therefore, we expect that there is a large band gap between the occupied carbon π states and the unoccupied carbon π states in both C 3 N and C 12 N. Our analysis for the wavefunction characters indicates that the highest occupied band near Γ is mainly contributed by nitrogen p z orbitals and there are some contributions from carbon p z orbitals when approaching the Brillouin zone boundary. The lowest unoccupied band is mostly contributed by the carbon p z orbitals. The gap between the highest occupied carbon π state and the CBM is larger than 1.5 eV (in LDA) in both C 3 N and C 12 N, which is in accord with the Clar's theory. The small band gaps for C 3 N and C 12 N are due to the presence of the nitrogen p z state (VBM state) above the highest occupied carbon π state.
C 3 N has a smaller band gap because the direct hopping between N p z orbitals in C 3 N leads to the larger dispersion of the highest occupied band. We note that the dilute limit of nitrogen substitutional doping in graphene was experimentally studied by Zhao et al. [28], where the extra electron of the nitrogen dopant was found to delocalized into the neighboring graphene lattice, resulting in the n-type behavior. In this work, we predict that collective nitrogen doping can induce exotic semiconducting behavior in graphene. Experimentally, some disorder in the nitrogen distribution might occur due to the thermal fluctuation. To see the effect of the disorder on the band gap, we have simulated a disordered configuration with 12.5% N concentration using a 128-atom special quasirandom structure [40]. The LDA band gap is found to be 0.08 eV. The real band gap should be larger because the LDA functional underestimates the band gap. This suggests that the disorder might change the magnitude the band gap, but shall not close it. Moreover, because the ordered C 3 N and C 12 N is more stable at low temperature than the disordered one, we believe reasonably large band gap should exist in nitrogen doped graphene synthesized at low temperature.
It is well-known that LDA seriously underestimates the band gap of semiconductors, thus we calculate the electronic structures of C 3 N and C 12 N by employing the screened Heyd-Scuseria-Ernzerhof 06 (HSE06) hybrid functional [41][42][43], which was shown to give a good band gap for many semiconductors including graphene related π systems [44]. Our

Potential Organic Solar Cell
The fact that C 12 N has a direct band gap of 0.98 eV suggests that C 12 N might be a promising solar cell absorption material. To investigate this possibility, we calculate the imaginary part of the frequency dependent dielectric function via summation over pairs of occupied and empty states without considering the local field effects using the HSE06 screened hybrid function. Fig. 5(a) shows the calculated result (ǫ 2 xx ) for one of the diagonal parts of the in-plane components. We can see that there is a very strong peak at around This suggests that C 12 N can act as a donor. If an appropriate acceptor material can be found, we can construct a photovoltaic system with a similar architecture as the P3HT/C 60 -PCBM bulk heterojunction solar cell. We find from HSE06 calculations that the the VBM and CBM levels of 6-AGNR are −4.67 eV and −3.15 eV, respectively. Therefore, the band alignment between C 12 N and 6-AGNR is of type-II, as shown in Fig. 5(b). In this proposed C 12 N/6-AGNR solar cell, C 12 N can absorb the solar energy between 0.98 eV and ∼1.50 eV.
And 6-AGNR can absorb the solar energy above the band gap 1.52 eV. The eletron-hole pair excited by photo can be separated by the field due to the type-II band alignment.

Conclusion
In summary, we have studied the interaction between substitutional nitrogen atoms in graphene and find that the nearest neighbor interaction between nitrogen dopants is highly repulsive because of the strong electrostatic repulsion between nitrogen atoms. Our calculations reveal two stable ordered N doped graphene structures C 3 N and C 12 N. Both structures are semiconducting, in contrast to the common belief that nitrogen atoms simply dope electrons to the graphene Dirac cone. In particular, C 12 N has a direct band gap of 0.98 eV.
The substitutional nitrogen raises the VBM and CBM level in N doped graphene, leading to a type-II band alignment between C 12 N and 6-AGNR. We propose that C 12 N/6-AGNR is a promising organic solar cell with C 12 N as donor and 6-AGNR as acceptor, respectively.
Our study suggests that nitrogen doping might be a promising way to open a band gap in grapene for the application in electronics and photovoltaics.