First-principles calculation of the electronic properties of graphene clusters doped with nitrogen and boron: Analysis of catalytic activity for the oxygen reduction reaction

Recent studies suggest that the carbon-alloy catalyst with doped nitrogen may be a powerful candidate for cathode catalyst of fuel cell. In this paper, we aim to clarify the microscopic mechanisms of the enhancement in the catalyst activity caused by nitrogen doping using a simple graphene cluster model. Our analysis is based on the density-functional electronic-structure calculations. We analyze modiﬁcations in the electronic structures and the energetical stability for some different conﬁgurations of N doping. We extend the analysis to the case of codoping of nitrogen and boron and propose two possible scenarios explaining the further enhancement of catalytic activity by N and B codoping.


I. INTRODUCTION
Due to serious global environmental problems and also the limitation in natural resources, the polymer electrolyte fuel cell ͑PEFC͒ has been regarded as one of the promising candidates for the future energy source, particularly for automobiles. It has many advantages such as high efficiency, low operation temperature, and, in principle, no environment load. However, as the reaction rate of oxygen reduction reaction ͑ORR͒ at the cathode is much slower than that of the hydrogen oxidation reaction at the anode, plenty amount of Pt is currently used as cathode catalyst for ORR. This hinders the wide use of PEFC because of the high cost and limited natural resources of Pt. Although there are some other technical issues to enhance the efficiency and durability of PEFC, search for cathode catalysts alternative to Pt is one of the most important tasks in the research of PEFC.
Several groups have reported that high ORR activities are found in the carbon-based materials doped with a certain amount of nitrogen ͑and sometimes also boron͒. 1- 18 We call the carbon-based materials containing some different elements carbon alloys. 19 Although the research on catalytic activities of carbon alloys with N doping started in 1926 ͑Ref. 1͒ and the catalytic activity enhancement by N-doped carbon in the fuel-cell cathode was reported in 1960s, 2 intensive study started only recently aiming at replacing the Ptbased catalyst with carbon-alloy catalyst ͑CAC͒. It is indeed rather surprising that the performance of CAC is getting close to that of Pt-based catalyst. Because carbon and nitrogen are abundant in nature, CAC is now a promising candidate for the Pt-free catalyst.
However, in order to design and produce CAC with higher performance, the detailed knowledge about the reaction site and the reaction mechanism is needed. Some of CACs have been synthesized by using carbonization processes of carbon polymers and nitrogen precursors such as phthalocyanine containing transition-metal elements ͑Fe, Co, Ni, and Cu͒. 4,11 Due to this synthesis process, one may suspect that the high performance of CAC is due to the presence of those transition-metal elements. Although there may be cases where the transition-metal elements play crucial roles, some analyses show that the CAC activity is not at all affected with acid washing which clearly decreases the amount of transition-metal elements in the sample. 11,18,20 The result implies that N-doped CAC without transition-metal elements can exhibit high catalytic ORR activity. 20 Transition-metal elements in the starting materials may help other elements, C and N, to form some particular structure suitable for the ORR activity. Therefore, in parallel with the study of the role of transition-metal elements, intensive studies have been performed to obtain insight into the characteristic features of nitrogen doping into CAC.
X-ray photoelectron spectroscopy ͑XPS͒ ͑Refs. 13, 15, and 18͒ and x-ray absorption spectroscopy ͑XAS͒ ͑Ref. 17͒ measurements suggest that nitrogen atoms in the carbon networks could occupy three or four chemically different sites, i.e., pyridine ͑or pyridiniumlike͒, pyrrolelike, graphitelike, and oxide. The variation in the relative intensity of these signals was found to be correlated with the ORR activity. These experimental analyses 17,18 suggest that the graphitelike nitrogen has a strong tendency of enhancing the ORR activity while the pyridinic nitrogen 21 ͑coordinated with two car-bon atoms͒ may have a tendency to suppress the ORR activity. With the present resolution of XPS and XAS, it is not possible to distinguish between the graphitelike nitrogen near the zigzag edge and the one inside the graphite sheet.
We have also performed the first-principles moleculardynamics simulations to analyze the reaction site and reaction process. 22 This analysis supports the experimental results. 17,18 In the presence of pyridinic nitrogen, the activation barrier becomes much higher for the O 2 -molecule adsorption to the carbon next to the nitrogen and also to the nitrogen itself. As for the graphitelike nitrogen, our simulation suggests that only the one next to the zigzag edge ͑for example, N2 in Fig. 8͒ enhances the O 2 adsorption and that the subsequent ORR process leads to complete reduction in O 2 to two H 2 O with a small activation barrier of about 5 kcal/mol. On the other hand, the graphitelike nitrogen inside the graphene sheet does not contribute to ORR activity.
In the present work we study the relative stability of various configurations of single N dopant in a graphene cluster to supplement the analysis in our previous paper. 22 We give detailed discussion on the modification of the electronic structure caused by N doping and particularly the dopant-site dependence of the electronic structure. We then study the interaction between two doped N atoms and obtain results suggesting that two doped N atoms repel each other. We extend our study to B dopant, being stimulated by the work by Ozaki et al., 6,7 which showed that the ORR activity of CAC is enhanced further by the coexistence of N and B. Starting with single B dopant cases, we study the interaction between N and B and also the electronic structure with N and B coexisting in a graphene cluster.
Here we make some remarks on our work on structural stability. First, we do not take into account possible structural modifications of zigzag edge recently discussed by Koskinen et al. 23 and Wassmann et al. 24 assuming that the chemical potential of hydrogen molecule is fairly low and also considering the situation described below. Second, the actual sample can have different structures depending on the starting materials and the synthesis processes. The structure, particularly the dopant structure, seems to be strongly affected by the transition-metal ions in the starting materials. 11,17,18 After synthesizing the sample, it is sometimes washed by acid solution in order to remove transition-metal ions. This implies that the structures in the actual sample may not be in the thermal equilibrium. In the present work, the transitionmetal ions are not taken into account. Taking account of these facts, we admit that our argument on the structural stability of dopants may have some limited meaning in the actual samples and we therefore study the electronic structures even for the atomic configurations which may not be so stable within the present model of graphene clusters.

II. COMPUTATIONAL METHOD
The first-principles calculations based on densityfunctional theories 25,26 and norm-conserving pseudopotential method 27,28 are performed by OPENMX ͑Ref. 29͒ code. The generalized gradient approximation proposed by Perdew, Burke, and Ernzerhof 30 is adopted. The basis functions con-sist of pseudoatomic orbitals 31,32 are specified as s 2 p 2 d 1 for carbon, nitrogen and boron, and s 2 p 1 for hydrogen atoms. The cutoff radii of all the orbitals are 4.5 a.u. The real-space grid technique is used with an energy cutoff of 150 Ry in numerical integrations of the Kohn-Sham equation and the Poisson equation is solved with the fast Fouriertransformation algorithm. The graphene cluster as shown in Fig. 1 of about 20 Å ϫ 20 Å with all the edge carbons terminated with monohydrogen is used for the simulation. The cluster consists of ten zigzag chains each of which contains 19 carbon atoms. As a whole, the cluster contains 190 carbon atoms and 38 hydrogen atoms. A periodic boundary condition is used and the separation between neighboring clusters is not less than 12 Å in all three dimensions. The atomic structures are fully relaxed in all calculations.

III. RESULTS AND DISCUSSIONS
A. Electronic structure and magnetic state of pure graphene cluster First we studied the electronic structure of a pure graphene cluster of Fig. 1. It is well known that an infinite graphene sheet is a zero gap semiconductor and that the band energy disperses linearly with the wave vector near the Fermi level. In a graphene cluster or a ribbon, special edge states exist along the zigzag edges, forming the flat bands at the Fermi level in the zigzag-type graphitic ribbon with H-terminated edges. [33][34][35][36] The presence of such a flat band is one of the typical conditions for the appearance of ferromagnetism and in fact the single zigzag edge becomes ferromagnetic. For a zigzag-type ribbon, two zigzag edges exist and the ferromagnetic moment on one edge couples antiferromagnetically with the one on another edge. The situation is nearly the same in the graphene cluster and the spin-density distribution in our cluster is shown in Fig. 1, where the con- tour for a given spin-density value ͑0.003 B / Bohr 3 ͒ is plotted. Clearly the magnetic moment is well localized along the zigzag edge reflecting the character of zigzag edge states. The antiferromagnetic configuration is more stable than the ferromagnetic one between the two zigzag edges by 0.045 eV per cluster. According to the calculation 37 using ribbon configurations, the corresponding energy difference is 0.004 eV per a pair of edge carbons on two sides of a ribbon which consists of eight zigzag chains across the ribbon. Our cluster has a corresponding width of ten zigzag chains and nine pairs of edge carbons. Therefore, our energy difference per cluster may correspond to 0.005 eV per a pair of edge carbons. As our cluster has a slightly larger width, our energy difference is a bit larger than that of Ref. 37.
In the cluster model, all the carbon atoms along the edge are not equivalent and the magnetic moment varies along the zigzag edge as shown in Fig. 2. Except on the edge carbons near the corner, the magnetic moment is about 0.26Ϯ 0.2 B . This value is consistent with the result of other calculations. 33,38,39 In the following sections, we will discuss the modifications in the electronic structure caused by doping of N and B. As a preparation for such discussions, we present in Fig. 3 detailed information about orbitals near the Fermi level. Due to the antiferromagnetic spin configuration shown in Fig. 1, a given orbital in the up-spin state has its counterpart in the  down-spin state with the same energy and they are mutually connected by the mirror-symmetry operation with respect to the bisecting plane denoted by m y in Fig. 1. The unoccupied orbitals a, b, and c in the up-spin state and orbitals h, i, and j in the down-spin state are zigzag edge states. The orbitals with higher energies are extended over the whole cluster. Similarly, the occupied orbitals d, e, and f in the up-spin state and orbitals k, l, and m in the down-spin state are zigzag edge states. The deeper occupied orbitals are extended over the whole cluster. Only these zigzag edge orbitals are spin polarized. The occupied orbitals d, e, and f and the unoccupied orbitals h, i, and j have weight on the atomic sites with net up-spin moment in Fig. 1. The orbital h is a counterpart of the orbital f and their energy separation ͑0.75 eV͒ comes from the exchange splitting. Similarly the exchange splitting is 0.54 eV between e and j and 0.48 eV between d and i. The same situation is seen between occupied down-spin orbitals and unoccupied up-spin orbitals which are localized in the right side of the cluster.
The total and local density of states ͑TDOS and LDOS͒ on some specific carbon atoms indicated in the right panel of Fig. 4 are shown in the left panel. The atom denoted with Cc is nearly at the center of the cluster and C3Ј is a counter atom of C3 located at the opposite side of the cluster. Only atoms along zigzag edges have significant peaks in the LDOS near the Fermi level. Clearly C3Ј and C3 have antiparallel magnetic moments. The number of electrons at each site for each spin state and the sum over spin states are also shown in Fig. 4. These are the result of Mulliken population analysis. At the center of the cluster ͑Cc͒, the number of electrons is nearly 4.0 being equal to the valence number of a C atom. For C atoms close to edges, the number of elec-trons is generally smaller than 4.0 by 0.1 or 0.2. This does not necessarily imply that the edge C atoms are positively charged. There is always some ambiguity in the assignment of charge states. In the Mulliken analysis, the overlap charge is evenly distributed to the neighboring atoms and this procedure produces smaller number of electrons at the edge C atoms and larger number of electrons at the hydrogen atoms.

Energetical stability of single N dopant
As mentioned in Sec. I, CAC with N doping is known to exhibit high catalytic activities for ORR. We have so far analyzed the role of N doping with the first-principles molecular-dynamics simulations 22 and also tried to identify experimentally 17,18 the particular configuration of doped N which promotes ORR. These analyses suggest that among some different configurations of doped N, those N located at the sites next to the zigzag edge activate the neighboring two carbons along the zigzag edge. On the other hand, the pyridinic N located along the zigzag edge may negatively affect ORR activity though it is abundant in most samples.
In this section, we analyze the stability of some possible configurations of graphene cluster with a single N dopant and also their electronic structures. As was already pointed out in Sec. I, the atomic structure of a sample depends on the starting materials and the synthetic process and may not necessarily correspond to the thermal equilibrium. Nevertheless, the structural stability analysis based on the total-energy calculation and the detailed study of the electronic structure may provide us with hints how to improve the CAC performance.
First we checked the stability of hydrogen termination for an N atom sitting along the zigzag edge as well as the armchair edge. For a pyridine molecule, N is not terminated with hydrogen but pyridine can be easily protonated to form pyridinium ion. Therefore, whether the N atom along the zigzag edge is hydrogen terminated or not may be rather subtle. Table I shows the result of total-energy comparison for hydrogen termination for some cases, for all of which only the    charge-neutral case is considered. In the second column of the table, the total-energy difference between hydrogenterminated N case and a case with a hydrogen molecule and nonhydrogen-terminated N. The negative ͑positive͒ value in the second column means that hydrogen termination is energetically stable ͑unstable͒. The result of Table I implies that the N atom along the edge of a graphene cluster, no matter along the zigzag edge or armchair edge or the corner site, has a tendency to hydrogen termination while pyridine has an N atom without hydrogen termination. This is consistent with the other first-principles calculation and experimental data. 40,41 In the following discussions, whenever N is located along graphene edges, we assume that it is hydrogen terminated. Note that the basis-set superposition error ͑BSSE͒ is taken into account in this calculation. The variation in the total energy with a C atom replaced with an N atom at different sites in the graphene cluster is listed in Fig. 5. The total energy for an N atom located near the center of the graphene cluster is taken as the reference of energy. The results show that the most favorable site for a doped N atom should be along the zigzag edge. The site dependence of the formation energy of N dopant in the carbon nanoribbons was studied by Yu et al. 42 and the present result is consistent with theirs. The strong stabilization of N along the zigzag edge ͑pyridinic N͒ has close relation with the existence of zigzag edge states near the Fermi level. Although the carbons along the zigzag edge are partially stabilized by the spin polarization of the zigzag edge state, they still have high-energy levels near the Fermi level and are energetically rather unstable. Therefore, the replacement of the zigzag edge C with N is energetically favorable because the high-energy levels are pulled down to a lower-energy region due to stronger attractive potential of N. The plot of the stabilization energy of N against the edge position in Fig.   2 clearly indicates such close correlation between the presence of edge state and the stabilization of N. The larger the magnetic moment on the edge carbon in a pure graphene cluster, the stronger the stabilization of N substitution.
In our previous paper, 22 we showed that an N atom located at a site next to the zigzag edge ͑for example, N2 in Fig. 8͒ activates the neighboring C atoms ͑C1 and C3 in Fig.  8͒ along the zigzag edge while an N atom along the zigzag edge ͑for example, N1 in Fig. 6͒ suppresses the activity of the neighboring edge C atoms. The present result shows that the chemically inactive configuration is much more stable than the chemically active configuration. Anyway, this may be an example of a general trend in the relation between energetical stability and the chemical reactivity.

Modifications in the electronic structure by a single N dopant
In this section, we give a detailed analysis on the electronic structures of N1-type and N2-type configurations of doped N. We will point out that the perturbation in the electronic structure by N doping is quite different between N1 and N2 and give a clear picture to understand the underlying mechanism producing the difference.   ber of electrons at C2 decreases by about 0.05. These results suggest that the charge state of C atoms around N1 is virtually unchanged from the one in the pure graphene.
In order to understand the mechanism of electronicstructure modification caused by N1, we give the detailed information about orbitals near the Fermi level in Fig. 7. As there are no mirror symmetries in the system, the symmetry between up-spin and down-spin states is also broken. First, we pay attention to the up-spin states. As orbitals a, b, and c for the unoccupied up-spin states are mostly localized in the right side of the cluster, they are not affected by N1 whose position is denoted by a filled circle in the right middle panel. As for the occupied up-spin states, the original three orbitals d, e, and f of Fig. 3 are reconstructed to form two orbitals d and e of Fig. 7 which have vanishing amplitude at N1. These two orbitals have only minor energy shifts from the corresponding original orbitals d and e of Fig. 3. The third orbital f of Fig. 7 which has significant amplitude at N1 has a relatively large energy shift. The orbital g is of extended state origin and modified from the corresponding one in Fig. 3 by the attractive potential of N1. Anyway, the number of occupied states does not change in the up-spin state. The situation is different in the down-spin state. As the unoccupied down-spin orbitals are localized in the left side of the cluster ͑see orbitals h, i, and j in Fig. 3͒, they are affected significantly by N1. As pointed out already for the pure graphene case in Fig. 3, the orbitals h, i, and j are the counter part of orbitals d, e, and f. Therefore the same story for the orbitals d, e, and f is applied to orbitals h, i, and j. Now the two orbitals h and i which have vanishing amplitude at N1 are little affected by N1 and stay above the Fermi level. On the other hand, the third one must have significant amplitude at N1 and is pulled down to be occupied as the orbital m in Fig. 7. The orbitals k, l, and m of the pure graphene cluster in Fig. 3 have no amplitude in the left side of the cluster and are little affected by N1 and stay nearly in the same state in Fig.  7 as orbitals j, k, and l. The orbital n has the same character as the orbital g. As explained already, the deep potential of N1 brings one of the originally unoccupied down-spin states into a lower energy to be occupied and the total charge neutrality is maintained. Figure 8 shows TDOS and LDOS for a cluster containing one N atom as shown as N2 in the right panel of Fig. 8. By comparing the LDOS in Fig. 8 with that in Fig. 6, we note significant difference in the electronic structure between N2 and N1 cases. In the energy range for the occupied states in   other hand at C1 and C3, the originally unoccupied downspin edge states in Fig. 4 are now pulled down to be occupied. The increased LDOS right below the Fermi level may be the origin of enhanced ORR catalytic activity of C1 and C3. 22 This result was pointed out already by Yu et al. 42 Naively one may think that an additional valence electron of N2 is donated to neighboring C atoms at the zigzag edge. However, we show in the following that this is not a correct picture. According to Mulliken population analysis, the numbers of electrons at C1 and C3 are slightly reduced from those of pure graphene shown in Fig. 4. Moreover, the number of electrons at N2 is slightly larger than that at N1. Below we give an interpretation to the mechanism leading to the change in the electronic structure from Fig. 4 to Fig. 8.
The most important difference between N1 and N2 is that while C1 ͑N1͒ site has significant contribution from edge states, C2 ͑N2͒ site has vanishingly small contribution from them. Therefore, while N1 can directly modify the edge states, N2 cannot. The problem is basically the same as the effective filling of d bands of Ni, Pd, and Pt by nontransition element impurities and a two-step screening mechanism was proposed by Terakura and Kanamori. 43 The same screening mechanism is applicable to the case of N2 and the basic idea is schematically illustrated in Fig. 9. In the first step from ͑a͒ to ͑b͒ in Fig. 9, even the deeper potential of N2 than that of C2 ͑in Fig. 4͒ cannot pull down the unoccupied edge states but the wave functions in the occupied states are modified so as to increase the number of electrons at N2 to attain its local charge neutrality. This leads to the electron deficiency at C1 and C3 because the total number of electrons accommodated in the states below E F remains unchanged unless any new additional states are introduced below E F . The situation is schematically illustrated in ͑b͒ and at the right bottom corner in Fig. 9. In the second step from ͑b͒ to ͑c͒, the electron deficiency at C1 and C3 causes attractive potential there, which pulls down the unoccupied edge state below E F . With these two steps, the neutrality of the whole system as well as at each site is achieved. Figure 10 shows the orbitals close to the Fermi level in the case of N2. A remarkable aspect is that all the wave functions shown here have vanishingly small amplitude at N2, being consistent with the LDOS of Fig. 8. Therefore, the modification of wave functions is not controlled directly by the potential of N2 but by the screening mechanism of Fig. 9. In the up-spin state, the three unoccupied orbitals a, b, and c are the same as those of pure graphene and the N1 cases. Next three orbitals d, e, and f below the Fermi level are to some extent affected by the presence of N2 and may be obtained by reconstruction of the three corresponding orbitals in the pure graphene. The orbital g is little affected by N2. In the down-spin state, the orbitals h, i, and j are the counterpart of orbitals d, e, and f but are modified slightly in a different way. Among three orbitals h, i, and j, the orbital j with large amplitude at C1 and C3 are pulled down below the Fermi level. Orbitals k, l, m, and n are basically the same as those of pure graphene.

Interaction between two doped N atoms
In Sec. III B 1, we showed that a single-doped N atom prefers the site along the zigzag edge. In this section, we study the interaction between two doped N atoms. As there are too many possible configurations for two doped atoms, we study only the cases where one of the two doped atoms is always located at a particular site along the zigzag edge, the site 11 in Fig. 2. The sites are renumbered as in Fig. 11 for convenience in the following discussion and the site 11 in Fig. 2 is now site 0. The total energy of the configuration in which the first N atom is at the site 0 and the second N atom is nearly at the center of the cluster is taken as the reference. The dependence of the relative total energy on the site of the second N atom is shown in Fig. 11. The result seems to suggest that the second N atom also prefers the zigzag edge except the close neighbors of the first N atom. The interaction energy E int ͑i , j͒ between an N atom at site i and another N atom at some other site j is obtained with the following equation: The first term in the right-hand side is the total energy of the pure graphene cluster and the second term is that of the case with two doped N atoms at site i and site j. The third and fourth terms denote the total energy of the case with a single N atom. With i fixed to be the site 0, E int ͑i =0, j͒ is shown in the inset of Fig. 11 for some different j. Clearly, if two doped N atoms are close, they interact repulsively and for farther interatomic distance beyond j = 5, they become weakly at-  9. ͑Color online͒ A screening mechanism for N2 in Fig. 8. ͑a͒ The original LDOS and the carbon potential in pure graphene cluster are schematically shown. "edge C" denotes a carbon atom along a zigzag edge, for example, C1 and C3, and "edge-1 C" denotes a carbon atom located at a site next to a zigzag edge, for example, C2 in Fig. 4. ͑b͒ "edge-1 N" denotes N2. A deeper potential of N than C cannot pull down the unoccupied edge state but increases the amplitude of wave functions of occupied states at N2 so that the charge neutrality is satisfied at N2. This causes the deficiency in the electron number at the neighboring carbons C1 and C3. The situation is illustrated in the right bottom panel. ͑c͒ The electron deficiency at C1 and C3 produces attractive potential which pulls down the unoccupied edge states to be filled. The charge neutrality at every site is thus approximately achieved. tractive or almost noninteracting. Therefore the probability of having two doped N atoms at two neighboring zigzag sites will be quite low and the model used in the analysis of oxygen adsorption process may be unlikely. 44

C. Effects of codoping of boron and nitrogen
As was mentioned in Sec. I, CAC with codoping of N and B exhibits further enhancement of ORR activity. 6,7 Below we propose two possible scenarios for the enhancement of ORR by codoping of N and B. One is that B may help N to be located at the favorable sites for ORR, for example, N2.
Another is compensation of the effect of the pyridinic N by forming an N-B complex. In this section we study the role of B doping and explore these possibilities. We start the discussion with a single B dopant case.

Energetical stability and electronic structure of single B dopant
First we studied the stability of hydrogen ͑H͒ termination for B sitting at the edge site and found that similarly to the N case, B also prefers H termination. Therefore in the following discussion, we always assume that B is terminated by H whenever it occupies the edge site. We then calculated the relative total energy of a system with a single B dopant at some different sites and the result is shown in Fig. 12. A doped B atom also prefers the zigzag edge sites with slightly reduced degree compared with the case of doped N. Figures 13 and 14 show TDOS and LDOS for graphene cluster with a single B dopant at two different sites. These figures correspond to Figs. 6 and 8 for a doped N atom and the effects of doped B can be understood just by changing the sign of the perturbing potential. With reference to C, the perturbing potential by B ͑N͒ dopant is repulsive ͑attractive͒. Fig. 13 with Fig. 6, we can clearly see that the peak structures in LDOS of N1 around −0.6 to −0.8 eV are shifted to around +0.6 to +0.8 eV for B1 case. In both cases, at the neighboring edge carbons, the edge-state contributions are significantly suppressed. The number of electrons at B1 seems to be appreciably less than 3.0. This is due to more diffuse character of 2p orbitals of B and an artifact of Mullikin population analysis. For B2 dopant, the neighboring edge carbons ͑C1 and C3͒ have enhanced unoccupied edge FIG. 11. ͑Color online͒ Graphene cluster with two N dopants. One of the two N dopants is always located at the site numbered "0" along the zigzag edge. The total energy of the system with two N dopants is measured with reference to the case where the second N dopant is located at the center of the cluster. The number shown at each site denotes the relative total energy ͑in eV͒ of the system where the second N dopant is located at the site. The inset plot denotes the interaction energy ͑in eV͒ between two N dopants, the first one being at the site labeled as 0 and the second one at another site. Our previous analysis suggests that the dopant N2 can activate the neighboring carbons, C1 and C3, for ORR. A possible reason for this higher catalytic activity is due to the larger number of electrons with higher energy at C1 and C3, which may make the electron transfer from C1 and C3 to the O 2 molecule easier. The effects of B2 and N2 have the electron-hole symmetry and B2 may activate reactions which induce electron transfer from reactant molecules to CAC, i.e., oxidation process.

Codoping of B and N
In the actual sample, the atomic configuration may not necessarily be in thermal equilibrium. Nevertheless energetical stability of N2 may enhance the catalytic activity further.
Here we first propose a possibility that the codoping of B and N may stabilize the N2-type configuration.
As B prefers the zigzag edge site, we first put a B atom at a particular edge site as marked with open circles in Figs. 15͑a͒-15͑d͒ and then put an N atom ͑filled circle͒ at the neighboring sites. Following the procedure described in Sec. III B 3, we calculate the relative total energy of a system with doped B and N for different sites for N dopant. The numbers outside parentheses in Fig. 15 denote the result. The numbers inside parentheses denote the interaction energy as evaluated with the equation similar to Eq. ͑1͒. It is important to note that strong attractive interaction works between neighboring B and N and that an N atom can be located at the site next to the zigzag edge ͑such as N2͒ fairly stably forming a B-N complex ͓Fig. 15͑c͔͒ though both B and N occupying separate zigzag edge sites ͓Fig. 15͑b͔͒ are slightly more stable. Figure 16 shows TDOS and LDOS for a B-N complex of Fig. 15͑c͒. At the edge carbon C3, the weight of zigzag edge states just below the Fermi level is enhanced compared with that for pure graphene shown in Fig. 4. Some weight of edge states can also be observed at B1. Our recent study of oxygen adsorption suggests that the reactivity of C3 in the present case may be not so strong as that in Fig. 8. This may be due to some compensation effect between B1 and N2.
Here we assume that by doping B first, the edge site will be preferentially occupied by B and then the subsequent N doping will keep the original B location by some kinematical reasons. Note, however, that the configuration of Fig. 15͑e͒ is more stable than those of Figs. 15͑b͒ and 15͑c͒. Therefore, in thermal equilibrium, the configuration of Fig. 15͑e͒ will be more favored. Unfortunately, this configuration whose elec-    tronic structure is shown in Fig. 17 seems to be not favorable for activating ORR of CAC. This is another example showing that the stable configuration will not exhibit high catalytic activity. Nevertheless, comparison between Figs. 6 and 17 suggests that the existence of B next to the edge N may recover the amplitude of zigzag edge states at the neighboring edge C atoms. Therefore, the negative effect of the edge N atoms on the ORR activity will be compensated to some extent by forming the N-B complex. This is another possible reason for the enhancement of ORR activity by N and B codoping.
The local structures of N and B in CAC were studied by XPS for N 1s and B 1s states. 5,6 The peaks were decomposed into four components corresponding to different environments. Particularly, in the analysis of N 1s XPS, it was concluded that the peak assigned as N-II is a new component associated with N-B codoping and assigned to B-N-C type local structures. 6,7 The relative weight of this component is about 0.2 irrespective of the concentration of B and always much larger than what we expect from the random distribution of B dopants, implying that N and B attract each other as suggested by the present analysis.
In order to go beyond these B-N complexes, we further studied the possibility of N-B-N complex such as the configuration ͑d͒ in Fig. 18. In this complex, we can expect that one of the combinations N-B or B-N in N-B-N may be responsible to the electronic compensation and that the remaining N will act as N2 in Fig. 8. This expectation is supported to some extent by the LDOS analysis shown in Fig. 19 and also by our separate study on oxygen adsorption at edge carbons in the presence of N-B-N complex. For example, Fig. 19 shows that at C3 and C5, the weight of edge states just below the Fermi level is significantly enhanced. Another interesting feature in Fig. 19 is that the LDOS of B1 has now significant contribution from the edge states suggesting that B1 itself may be chemically active. Despite these nice features expected for the N-B-N complex of Fig. 18͑d͒, the configuration of Fig. 18͑e͒ is energetically much more stable. As expected from the atomic configuration, this complex is rather harmful to the CAC activity for ORR, although the electronic structures are not shown here. Again the stable configuration will not have high chemical activity.
Detailed study of reactivity for B-N complexes will be published in a separate paper. For ͑a͒ to ͑d͒, the B dopant ͑open circle͒ is first located along the zigzag edge and the N dopant ͑filled circle͒ is put at some different sites. In the configuration ͑e͒, N is located along a zigzag edge site. The number outside parenthesis is the relative total energy of each configuration with the configuration ͑a͒ taken as the reference. The interaction energy between two dopants is estimated by an equation similar to Eq. ͑1͒ and is listed in the parenthesis.  FIG. 16. ͑Color online͒ TDOS and LDOS for a B-N complex in a graphene cluster. C3Ј is along the other zigzag edge. At the edge carbon C3 next to N, the weight of zigzag edge states just below the Fermi level is enhanced compared with that for pure graphene. Some weight of edge states can also be observed at B1.  FIG. 17. ͑Color online͒ TDOS and LDOS for a N-B complex in a graphene cluster. C3Ј is along the other zigzag edge. At the edge carbon C3 next to B, the weight of zigzag edge states just below the Fermi level is nearly the same as that for pure graphene.

IV. CONCLUSIONS
In this paper, we tried to clarify the microscopic mechanisms of the enhancement of ORR activity of CAC by N doping using a simple graphene cluster. We analyzed the energetical stability for some different configurations of doped N and interdopant interactions. Detailed analysis was also made for modifications in the electronic structure caused by the dopant atoms.
We found that the pyridinic N is the most stable configuration for doped N in a graphene cluster. In other words, the pyridinic N will stabilize the zigzag edges of graphene and will help the zigzag edges grow in the synthesis process. In this sense, pyridinic N may contribute to ORR activity of CAC. At the same time, we found that for a given zigzag edge, the presence of pyridinic N will have a negative effect on the ORR activity. On the other hand, while the doped N at the site next to the zigzag edge can activate the neighboring edge carbons, this configuration is energetically not so stable. We admit that the reactivity and atomic configuration of actual samples strongly depend on the starting materials and synthesis processes and that the final atomic configuration of a sample may not necessarily be in the thermal equilibrium. Nevertheless, it will be helpful to improve the catalytic activity of CAC if the configuration favorable to reactivity en-hancement is also energetically favorable. In order to realize such situation, we proposed codoping of B and N. Our analysis suggests that codoping will help N to be located at a particular site suitable for enhancing the catalytic activity of CAC or will compensate the negative effect for ORR caused by the pyridinic N.
A first-principles molecular-dynamics study is now underway to clarify the activity of some dopant configurations and also the subsequent ORR processes.  FIG. 19. ͑Color online͒ TDOS and LDOS for a ͑N-B, N͒ complex in a graphene cluster. The B dopant is located along the zigzag edge, and two N dopants are next to it. C3Ј is along the other zigzag edge. At the edge carbons C3 and C5 next to N, the weight of zigzag edge states near the Fermi level is enhanced compared with that of pure graphene cluster. Some weight of zigzag edge states can also be observed at B1.