Size Optimization of a N-Doped Graphene Nanocluster for the Oxygen Reduction Reaction.

N-Doped graphene nanoclusters (N-GNCs) are promising electrocatalysts for the oxygen reduction reaction (ORR) at the cathode of fuel cells. In this study, the dependence of the ORR activity on the size of N-GNCs was investigated using first-principles calculations based on density functional theory. The maximum electrode potential (U Max) was estimated from the free energy of the reaction intermediates of the ORR. U Max was predicted to show a volcanic trend with respect to the cluster size. The results suggest that C215H36N with a radius of 13.6 Å is the best candidate for ORRs and is better than platinum in terms of U Max. The volcano-shaped plot of U Max is attributed to the switch of the reaction step that determines U Max, which is caused by the destabilization of reaction intermediates. Such changes in the stability of the intermediates can be explained by the decrease in the local density of states at the reaction site, which is due to the development of the so-called edge state at the zigzag edge. The establishment of experimental techniques to control the cluster size and doping position will be the key to superior catalyst preparation in the future.

Introduction

Fuel cells are in great demand as eco-friendly energy systems that do not emit carbon dioxide. In the fuel cells, energy is obtained through the hydrogen oxidation reaction (HOR) at the anode and the oxygen reduction reaction (ORR) at the cathode. In order for fuel cells to spread widely throughout society, it is necessary to solve the lower reaction efficiency of the ORR compared with that of the HOR. Platinum-based alloy materials are widely used as cathode catalysts with high ORR activities. However, the high cost and low durability of platinum-based catalysts inhibit the propagation of fuel cells.[1] While the development of new catalysts to replace platinum is being vigorously pursued, it was recently reported that N-doped graphene exhibits a high ORR activity comparable to that of platinum.[2−16] Various claims have also been made about the local coordination of nitrogen atoms in graphene, but it has been suggested that nitrogen atoms prefer to be located near the zigzag edge of graphene.[17−21] Therefore, many theoretical researchers have been focusing on the structure and the electronic states of the nitrogen atoms at the edges of graphene.

It has been theoretically revealed that a N-doped graphene nanocluster (N-GNC) with a hexagonal shape exhibits a high catalytic activity for the ORR.[13,15,16,22] Recently, Ganyecz and Kállay systematically explored the effect of the position of N on the ORR activity in a GNC of a certain size;[16] scaling relations with regard to the free energy of the intermediates can be derived for N-GNCs. Zhang et al. investigated the ORR on N-GNC where the N atom was doped near the edge and concluded that the size of the N-GNC affects the ORR activity depending on the stability of the reaction intermediates.[23] Furthermore, it has been shown the ORR activity for N-GNCs in the vicinity of the edge depends strongly on doping and reaction sites.[24] In our previous study, the dependence of the ORR on the location of the N atom in the cluster was investigated.[22] When the N atom is doped at the midportion of the cluster, the GNC shows a relatively higher ORR activity, that is, a high maximum electrode potential (UMax), and a selectivity not for the two-electron (2e–) pathway but instead for the direct four-electron (4e–) pathway.[22] It would be interesting to know whether the wide bulk-like basal surfaces of the cluster or the localized electronic states at the edges have more of an effect on the value of UMax for the ORR. Our previous studies provided insight into the expression of high ORR activities. The effects of edge-localized electronic states do not promote the high ORR activity, but the ORR may be enhanced by another factor such as a kind of confinement effect.[22,24] However, only few attempts have so far been made to determine the impact of cluster size effect on the ORR activity of N-GNCs while focusing on the edge and the so-called confinement effect. Furthermore, the origin governing the stability of the reaction intermediates in N-GNC has not been discussed so far from an electronic points of view.

In this study, we investigate the ORR activity of N-GNCs with various sizes and edgeless N-doped graphene with periodic boundary conditions using the first-principles calculations based on density functional theory (DFT). We first investigate the effect of the physical boundary, namely the edge, on the ORR activity. Here, we further show the N-GNC confined by edges shows a higher ORR activity than infinite N-doped graphene. Finally, we predict both that there is an optimal cluster size in terms of the UMax of the ORR and the reaction selectivity for the direct four-electron pathway and discuss the electronic origin of the existence of the optimal cluster size for the ORR.

Computational Methods and Models

Computational Methods

For the finite system, namely N-GNCs, we used the Gaussian 09 code[25] and adopted the hybrid B3LYP functional[26,27] and the 6-31G(d,p) basis set in this code. The atomic structure of the N-GNC was optimized so that each component of the interatomic force was below 0.0003 Ha/Bohr. The ORR in carbon systems is known to proceed by two main pathways. One is the 4e– pathway in which O2 is reduced to H2O via OOH, O, and OH. The other pathway is the 2e– one, where the O2 molecule is reduced to H2O2 via OOH, O, and OH. In the 2e– pathway, the reduction process of H2O2 is eliminated because the H2O2 molecule is less likely to be adsorbed onto the N-GNC. The stability of the intermediates on the N-GNC determines the electrocatalytic activity of the reduction in the 4e– and the 2e– pathways. Unlike Pt-based catalysts, it is known that the energy barrier of the so-called dissociative pathway is very high for graphene-based catalysts,[28−30] so we considered only the so-called associative pathway in this study.

The electrocatalytic activity of the ORR was evaluated based on the so-called computational hydrogen electrode model.[31] The chemical potential for (H+ + e–) corresponds to that of 1/2H2 in the gas phase with respect to the standard hydrogen electrode (SHE). In this study, a pressure of 1 bar, pH 0, and T = 298 K were adopted for the evaluation of the free energy for each ORR process. The reaction free energy (ΔG) was calculated as following formula:where ΔG0 is the Gibbs free energy, ΔG is the electrode potential, ΔGpH the effect of the solvent, ΔGW is the stabilization energy by water, and ΔGfield is the effect of the local electric field near the surface of the electrode. ΔG0 consists of the following terms:where ΔE is the difference in the internal energy between the reaction intermediates and the final products. ΔZPE and TΔS, the zero point energy and the entropy, respectively, were calculated on the basis of the vibrational frequency calculation. ΔGpH was set to 0, i.e., we supposed that the ORR would take place under acidic conditions (pH 0). Here, we did not take ΔGfield into account because the absolute value of ΔGfield was negligibly small (≈10–2 eV).[32] Water is known to play a very important role in the oxygen reduction reaction in graphene systems.[22] Therefore, the effect of water was included as an effect of the solution, and ΔGW was assumed to be the stabilization energy by an H2O molecule. It has been revealed that the present values of ΔGW are quantitatively equivalent to those using other solvent models.[22] The modeling details for the electrochemical reaction process were described in our previous paper.[22] We evaluated UMax by drawing a free energy diagram with varying values of U. UMax is the maximum electrode potential at which all reaction processes are exothermic. The larger the value of UMax, the lower the overpotential.

To extract the effect of edges on the ORR activity for the in-plane of N-doped graphene, we also considered a periodic model of isolated-N-doped graphene without edges. A (8 × 8) supercell of graphene where one C atom in the unit cell was substituted by one N atom was used. The interaction energy between N atoms is considered to be almost negligible because the nitrogen atoms are separated from each other by 20.3 Å in the plane.[33] In the supercell geometry, a vacuum region in the surface vertical direction was set to 20 Å to decouple the periodic images. The Vienna ab initio simulation package (VASP) code[34,35] was used for the supercell system. A plane-wave basis set with an energy cutoff of 600 eV was applied to the wave function for all calculations using the projector-augmented wave (PAW) method.[36,37] The so-called PBE (Perdew, Burke, and Ernzerhof)-type generalized gradient functional was adopted as the exchange-correlation functional.[38] Integration in the reciprocal space over the 2D Brillouin zone was carried out using two independent k-points in the irreducible Brillouin zone of the (8 × 8) supercell. The atomic position was optimized such that the force acting on each atom became less than 1.0 × 10–2 eV/Å.

Computational Models

Since hexagonal-shaped GNCs with various sizes have been fabricated experimentally so far,[39−48] we adopted the zigzag edge N-GNCs with the hexagonal shape. Figure shows the models of the N-GNCs. Nitrogen is known to be doped in graphene in a variety of configurations other than the graphitic configuration.[49] Recently, it has been shown that nitrogen atoms can be doped into in-plane graphitic sites instead of edges using a certain nonequilibrium synthesis method called the solution plasma method.[50,51] Further, since our previous study verified that the ORR performance of graphitic N-GNC is best when N is doped as far away from the edge as possible,[22] the nitrogen atom was located in the central part of the cluster. The reaction sites were assumed to be on the C atom adjacent to the doped N atom[22] on which the reaction intermediates (OOH, O, and OH) for the ORR are adsorbed most stably. The reaction sites are denoted by letters, specifically “a” and b (see Figure ). Calculation models referred to “CHN–a(b)”, where “X” and “Y” are the numbers of C and H atoms, respectively, and “a” or “b” is a reaction site. It is noted that the change in the cluster size corresponds to the change in the nitrogen concentration, since one nitrogen atom is doped in each cluster.

Figure 1

Models of the following N-GNCs: (a) C53H18N, (b) C95H24N, (c) C149H30N, and (d) C215H36N. The white, gray, and blue balls indicate H, C, and N atoms, respectively. The red circles indicate reaction sites “a” and “b”.

Results and Discussion

Figure shows the free-energy diagrams for the model of C53H18N-b as an example of a relatively small cluster. The horizontal axis in this figure shows that the ORR progresses from left to right; the reaction free energy, ΔG, for the ORR intermediates is traced from the initial O2 adsorption toward the final H2O or H2O2 production at zero cell potential (U = 0 V) and the equilibrium potential (Ueq = 1.23 V). Here, UMax is defined as the maximum electrode potential such that all reaction steps are exothermic. It can be seen that for the 4e– pathway the ORR at U = 0 V proceeds while gradually reducing the free energy, that is, all reactions spontaneously proceed toward H2O generation. In this case, the UMax value was calculated as 0.58 V/SHE for the 4e– pathway. On the other hand, it can be seen for the 2e– pathway that the H2O2 generation process from OOH adsorption becomes uphill even at U = 0 V. As a result, UMax is apparently a negative value (−0.19 V). This means that the ORR for the 2e– pathway is aborted at the OOH adsorption step and the reaction does not proceed any further. Therefore, the model C53H18N-b has a high selectivity for the 4e– pathway. In addition, it was found that H2O2 is energetically unstable, at least in the vicinity of graphitic N, and leaves the surface without barriers. Therefore, we expect that corrosion of the GNC is unlikely if graphitic N can be prepared, even if the reaction through the 2e– pathway occurs. All models also have the selectivity for the 4e– pathway, i.e., UMax(4e–) becomes larger than UMax(2e–). Diagrams of the reaction step for all other models are shown in the Supporting Information (Figures S1–S4).

Figure 2

ORR diagrams for the model C53H18N-b under acidic conditions (a) for the 4e– pathway at U = 0 V (zero cell potential), Ueq = 1.23 V (equilibrium potential), and UMax = +0.58 V (maximum potential) at which all reaction steps are exothermic and (b) for the 2e– pathway at U = 0 V, Ueq = 0.68 V, and UMax = −0.19 V. The value of U was estimated with respect to the SHE.

We estimated values of values of UMax for all models, C53H18N–C215H36N. Figure shows the change in UMax as a function of the cluster size. Circles and crosses show the calculated values of UMax for the 4e– and the 2e– pathways, respectively. The UMax values for the 4e– pathway are always higher than those for the 2e– pathway regardless of the cluster size. This means that the N-GNCs with sizes from C53H18N to C215H36N have a reaction selectivity for the 4e– pathway. As the cluster size increases, the values of UMax for the 4e– and the 2e– pathways increase. For the 4e– pathway, a high UMax value means a high capability of ORR catalysis. On the other hand, the increase of the UMax value for the 2e– pathway causes H2O2 generation, leading to the low durability of the electrocatalyst. The theoretical limit of UMax equal to the equilibrium potentials for the ORR under the acidic conditions are 1.23 and 0.68 V for the 4e– and the 2e– pathways, respectively. For the N-GNC larger than C215H36N, it is not obvious that the value of UMax approaches the theoretical limit value.

Figure 3

UMax values of the N-GNCs with various sizes. The circles and the crosses show the UMax values calculated for the 4e– and the 2e– pathways, respectively. The solid and dashed lines show values of UMax calculated from the fitted and extrapolated values for ΔGdiff, respectively.

We also evaluated the values of UMax for the periodic model of an isolated N-doped graphene sheet without edges (shown in the Supporting Information Figure S5). The UMax values of the periodic model for the 4e– and the 2e– pathways are nearly the same (∼0.4 V), showing no selectivity for the 4e– pathway. In both reaction pathways, the process of OOH adsorption determines the value of UMax. What has to be noticed is that the UMax value of the periodic model is lower than those of N-GNCs, even though the values of UMax for N-GNCs increase with the increasing cluster sizes. This suggests that there exists a N-GNC with the maximum value of UMax.

To find the optimal size of N-GNCs, we estimated the dependence of UMax on the cluster radius by extrapolating the free-energy values of reaction intermediates for the models of C53H18N–C215H36N. The adsorption stability of the reaction intermediates was evaluated from the relative free energies of the intermediates, namely ΔGOOH, ΔGO, and ΔGOH. UMax was determined by the difference in ΔG at each reaction step, ΔGdiff. Plots in Figure show the ΔGdiff values for the 4e– and the 2e– pathways. The ΔGdiff value becomes the smallest for the step between OH adsorption and H2O generation (ΔGOH – ΔGH), which determines the value of UMax for the 4e– pathway. Here, we consider a linear interpolation of the change in ΔGdiff. The dashed lines in Figure show the linearly extrapolated values of ΔGdiff for the reaction intermediates. While the value of ΔGOH – ΔGH increases with the increasing cluster size, the value of 4.92 – ΔGOOH decreases. The smallest value of ΔGdiff switches from ΔGOH – ΔGH to 4.92 – ΔGOOH with a radius of about 14 Å for the 4e– pathway. As for the 2e– pathway, the smallest value of ΔGdiff switches from ΔGOOH – ΔGH to 1.36 – ΔGOOH with a radius of about 21 Å. As a consequence, the maximum electrode potential shows a volcano-shaped plot because the process that minimizes ΔGdiff is switched. In addition, these results are consistent with the result for the periodic model: the UMax values for the periodic model are dominated by the adsorption of OOH for both 4e– and 2e– pathways, resulting in the no selectivity for the 4e– pathway and UMax values lower than those of N-GNCs. Furthermore, as can be seen in Figure , the free energies of OH and OOH are linear with a high correlation coefficient, indicating a high linear relationship between them.[16,52,53]

Figure 4

ΔGdiff values of the reaction intermediates for the models for the (a) 4e– and (b) 2e– pathways. In panel a, the squares, rhombuses, triangles, and the circles represent 4.92 – ΔGOOH, ΔGOOH – ΔGO, ΔGO – ΔGOH, and ΔGOH – ΔGH values for the 4e– pathway, respectively. In panel b, the squares and the circles represent 1.36 – ΔGOOH and ΔGOOH – ΔGH values for the 2e– pathway, respectively. The solid and the dashed lines show linearly fitted and extrapolated values of ΔGdiff, respectively. Values for the coefficient of determination (R2) were 0.80, 0.31, 0.21, 0.75, 0.79, and 0.79 for 4.92 – ΔGOOH, ΔGOOH – ΔGO, ΔGO – ΔGOH, and ΔGOH – ΔGH for 4e– and 1.36 – ΔGOOH and ΔGOOH – ΔGH for 2e–, respectively.

The solid and dashed lines in Figure show fitted and extrapolated UMax values based on the smallest ΔGdiff obtained from the fitted and extrapolated ΔG values for the reaction intermediates, respectively. The dependence of UMax on the cluster radius leads to a volcano-shaped trend, which has a maximum value. The maximum UMax values are estimated to be 1.07 and 0.68 V with a radius of 13.8 and 21.0 Å for the 4e– and the 2e– pathways, respectively. For the range of cluster radii less than about 21 Å, the UMax values for the 4e– pathway are always higher than those for the 2e– pathway, showing the reaction selectivity for the 4e– pathway. On the other hand, for the range of cluster radii over about 21 Å, the UMax values for the 4e– and the 2e– pathways are nearly the same. This means that there is no selectivity of the reaction pathway. Thus, it is suggested that C215H36N with a radius of 13.6 Å exhibits the best performance for the ORR from the viewpoint of UMax. Indeed, the value of UMax = 1.05 V for the C215H36N surpasses that for platinum.[54] Furthermore, Figure indicates that the 2e– reaction does not occur in clusters with a radius of about 9 Å or less. Therefore, a cluster with a radius of about 9 Å, i.e., C95H24N, is the best choice for the ORR in terms of the selectivity for the 4e– pathway.

Finally, we discuss the reason for the existence of the optimal cluster size for the ORR. A donor electron from the doped N atom enhances the chemical bonding between reaction intermediates and N-GNC.[24] As shown in Figure S6, as the cluster size increases, the local density of states at the reaction site decreases, which is conducive to chemical bonding with the reaction intermediates. This leads to the destabilization of the reaction intermediates. As a matter of fact, the value of ΔG for the reaction intermediate become higher with the increasing cluster size, as shown in the Figure S7. As the cluster size increases, the length of the straight portion of the edge increases. Infinitely long zigzag edges have been known to have the so-called edge states, which are localized just at the edges.[55] It can be interpreted that as the cluster size increases, the contribution of the edge states in the SOMO increases, resulting in a decrease in the density of states at the reaction site. Indeed, it has been revealed that the edge state develops in GNCs larger than C96H24.[21] If the cluster size is further increased and a complete edge state is generated, the ORR activity is expected to disappear because the extra electron provided by the nitrogen atom contributes completely to the formation of the edge state.[20] Therefore, the ORR activity of a larger cluster is not expected to converge to the one of the periodic system.

Conclusions

The ORR activity for N-GNCs with the various sizes has been investigated using the first-principles calculations within DFT. The UMax value of the N-GNC reaches a maximum with a radius of about 14 Å because the step that determines the UMax values switches from the H2O generation step to the OOH adsorption one. The maximum UMax value of N-GNCs for the 4e– pathway has been estimated to be 1.07 V with respect to the SHE, surpassing that of platinum. The C215H36N with a radius of 13.6 Å is the best size in view of UMax and is expected to be a potential electrocatalyst for a fuel cell. In terms of the selectivity for the 4e– pathway, the C95H24N cluster turned out to be the best choice for the ORR, although UMax was somewhat lower (∼0.8 V). Such a size-dependent ORR activity of the N-GNC is derived from the change in the confinement of a donor electron from the doped N atom. The adsorption energy of the reaction intermediates varies continuously because the spread of the supplied electrons from the doped nitrogen atoms to the surrounding atoms varies with the cluster size. This is due to the development of the so-called edge states as the cluster size increases, causing a decrease in the number of electrons contributing to the chemical bond of the reaction intermediate at the reaction site. As a result, according to the so-called Sabatier principle, C215H36N has the largest UMax, that is, the smallest overpotential in the case of a N-doped GNC.

Since N atoms thermodynamically prefer to be near the edge, it is experimentally difficult to fabricate N-GNCs containing N atoms inside clusters in a thermal equilibrium. Recently, however, attempts have been made to fabricate N-GNCs containing N atoms in the cluster plane using nonequilibrium processes such as the solution plasma processes.[56,57] If it succeeds, N-GNCs will be breakthrough catalysts for the ORR. The reaction when multiple nitrogen atoms are doped in one cluster is also of interest, but we would like to leave that for future research.