Cu - Nitrogen doped graphene (Cu–N/Gr) nanocomposite as cathode catalyst in fuel cells – DFT study

Novel Cu-nitrogen doped graphene nanocomposite catalysts are developed to investigate the Cu-nitrogen doped fuel cell cathode catalyst. Density functional theory calculations are performed using Gaussian 09w software to study the oxygen reduction reaction (ORR) on Cu-nitrogen doped graphene nanocomposite cathode catalyst in low-temperature fuel cells. Three different nanocomposite structures Cu2–N6/Gr, Cu2–N8/Gr and Cu–N4/Gr were considered in the acidic medium under standard conditions (298.15 K, 1 atm) in order to explore the properties of the fuel cell. The results showed that all structures are stable at the potential range 0–5.87 V. Formation energy, Mulliken charge and HOMO-LUMO energy calculations showed that Cu2–N6/Gr and Cu2–N8/Gr are more stable structure-wise, while free energy calculations showed that only Cu2–N8/Gr and Cu–N4/Gr structures support spontaneous ORR. The maximum cell potential under standard conditions was shown at 0.28 V and 0.49 V for Cu2–N8/Gr and Cu–N4/Gr respectively. From the calculations, the Cu2–N6/Gr and Cu2–N8/Gr structures are less favorable in H2O2 generation; however, Cu–N4/Gr showed the potential for H2O2 generation. In conclusion, Cu2–N8/Gr and Cu–N4/Gr are more favorable to ORR than Cu2–N6/Gr.


Introduction
Energy security and sustainability are the impetus for modern economics. Transformation of the energy ecosystem requires the search for cost-effective, secure, and sustainable energy technologies to replace incumbent energy sources. Fuel cells contribute well to this transformation as they evidently have their merits, particularly for compactness and zero emissions, such as proton exchange membrane fuel cells (PEMFCs) [1].
Low-temperature fuel cells operate at temperatures below 200 • C and a carbon support catalyst is deployed to improve the oxygen reduction reaction (ORR). Carbon support catalysts are not conducive to high temperature operations, leading to an increased rate of carbon corrosion and, therefore, degrading the carbon support [2]. Corrosion of the carbon support will directly affect the catalyst

Methodology
DFT calculations were performed with Gaussian 09w software using the B3LYP/3-21G basis set with non-periodic and nondispersion interaction as Bhatt et al. [21]. Defect structure is visualized using GaussView 6.0. All atoms in every structure were relaxed by optimization. A common single pristine graphene layer was developed, and it is used to develop other structures (Cu 2 -N 6 /Gr, Cu 2 -N 8 /Gr, and Cu-N 4 /Gr) to investigate ORR considering the defect structure at the catalyst surface. All structures were developed to contain only pyridinic nitrogen. After optimization, zero imaginary frequency for dual atom structures and one (1) imaginary frequency for the Cu-N 4 /Gr structure. Additionally, the optimized charges of all three catalysts were zero and singlet spin for all dual atom structures. The molecules used in the ORR steps were individually optimized. Bond lengths and bond angles were measured and compared with those of previous research.
Catalyst structures based on Cu confirmed their stability with the cathode potential. Stability was inspected using the formation energy (ΔE) defined accordingly in the Eq. (1) [16,21], Here, E graphene+(Ma− N b ) is energy of optimized graphene layer with Cu-N defect. The 'a' and 'b' are positive integers defining the selected Cu-N defects (a = 1,2 and b = 4, 6,8). M and N are Cu and nitrogen, respectively. The μ C and μ N are chemical potential of carbon defined as total energy per carbon atom for defect-free graphene, and the chemical potential of nitrogen defined as half of the total energy of N 2 molecule, respectively. x and y are the number of nitrogen atoms added and the carbon atoms removed during the defect formation, respectively. E graphene is the energy of optimized pristine graphene layer. E M is the total energy of M n + defined as Eq. (2), where, E(M) is the total energy of isolated M (M = Cu) in the gas phase and, n, e, U are the number of electron transfer (+2), electron charge, and external potential in order. Binding energies (BE) were calculated from the ORR steps in an acidic medium, as defined by Eq. (3) [16,20,21], Here, E defect+molecule is the total energy of molecules adsorbed by defect graphene. E defect is the total energy of defect graphene configuration, and E molecule is isolated molecule species (O 2 , O, H 2 O, OOH, OH, H 2 O 2 ). Negative signed binding energies (E defect+molecule < (E defect + E molecule )) indicate it is more favorable for molecules to be attached to the defect configuration. The formation of H 2 O 2 during ORR is also considered. ORR steps in H + medium are defined as follows [16], * indicates the configuration of the defect.
The free energies were calculated for each ORR step of all defects configurations as defined by Eq. (4) [16], ΔE is the energy from the DFT calculation to the relevant reaction step, and ΔZPE is the correction of zero-point energy is obtained from the NIST database and DFT calculations. T and S are the absolute temperature and entropy, respectively. Here, T = 298.15 K, and the entropy value was obtained from the NIST database. ΔG U = -eU where U and e are the electrode potential and the charge transferred, respectively. ΔG pH = K B T × ln 10 × pH, where K B is Boltzmann's constant and T = 298.15 K ΔG field is a contribution of the interaction of adsorbate with the local electric field in the electric double layer formed in the vicinity of the cathode, which is negligible according to a Nørskov et al. [24]. Corrections of zero-point energy and entropy values were obtained from frequency calculations for the molecule adsorb by defect graphene. The free energy vs. reaction coordinate graphs were drawn with different potentials (U) for each structure. Additionally, the free energy of O 2 , 4.92 eV, was obtained from the literature under standard conditions with reaction  [24,28]. The Mulliken charge population was applied to molecules as defined by Eq. (5) [29], The Q after and Q before are charges of the molecule (X) after and before adsorption, respectively. The energy gap (E g ), chemical hardness(η), chemical potential (μ), and electrophilicity index (ω) were calculated using the Koopman's principle for optimized structures as defined by Eqs. (6)-(9), respectively [30][31][32][33], Here, I and A are the ionization potential (≅ − E HOMO ) and the electron affinity (≅ − E LUMO ), respectively.

Formation energy
The formation energies of the Cu 2 -N 6 /Gr, Cu 2 -N 8 /Gr, and Cu-N 4 /Gr structures were − 23.60, − 23.49 and − 13.22 eV in order at zero potential (U = 0). These values were only applied under open circuit conditions, and they changed with different cathode potentials according to Eq. (1) [16,21]. Table 1 shows that the formation energy increased with increasing external potentials. A negative sign means energy value gain when a structure forms. Therefore, the negative formation energies are favorable for the formation of stable structures, and the formation energies of all the structures studied remain negative between the potential range of 0-5.87 V. The critical cell potential of negative formation energies was observed at 5.89, 5.87, and 6.60 V for Cu 2 -N 6 /Gr, Cu 2 -N 8 /Gr, and Cu-N 4 /Gr, respectively, see Table 1. The critical U value at zero formation energy demarcates the stability/instability of the cathode. The structures forming stability follow the increasing order of Cu-N 4 /Gr < Cu 2 -N 8 /Gr < Cu 2 -N 6 /Gr, as indicated by the external potential, U, at zero.

Table 2
Binding energies (BE) of Cu 2 -N 6 /Gr, Cu 2 -N 8 /Gr, Cu-N 4 /Gr structures with and without solvent energies (SE) [34] and shortest distance (d) between Cu-O, O-O atoms in angstroms (Å). The Mulliken charge of the adsorbate (Q x ) and the Mulliken charge of the Cu site (Q Cu ) for the shortest distance of the optimized structure. (01) and (02) in the Cu 2 -N 6 /Gr defect represent the different ORR paths 03-04-06-08 and 01-05-08 in Fig. 1 (Table 2), and with and without the water solvent energies obtained from the work of Sha et al. [34]. Kattel et al. and Bhatt et al. reported that the binding energies of O 2 to the structure of Co-N 4 /Gr and Cu-N 2 /Gr were − 0.67 and − 2.90 eV, respectively. As a first ORR step, O 2 binding to the defect shown in this study is significantly high, except for the Cu-N 4 /Gr structure in relation to Cu-N 2 /Gr [16,21]. Structure-wise Cu-N 2 /Gr has a higher O 2 binding energy than Cu-N 4 /Gr, also, Cu 2 -N 6 /Gr has a higher binding energy compared to the Cu 2 -N 8 /Gr structure. This indicates that adding more nitrogen atoms to Cu did not create a good O 2 adsorbent defect. The efficiency of the O 2 adsorbent in decreasing order is as follows, The binding energies of OH and OOH on the Cu 2 -N 6 /Gr structure are shown to be − 5.12 and − 4.13 eV and − 3.81 and − 3.27 eV, respectively, for binding onto the Cu 2 -N 8 /Gr structure. The binding energies for OH and OOH to the Cu-N 2 /Gr structure were reported to be − 2.68 and − 1.81 eV, respectively, while for the Cu-N 4 /Gr structure were − 1.99 and − 1.70 eV [21]. The OH and OOH binding energies also signify the above N atoms and the idea of a stable defect. The energy values, shown in Table 2, also suggest that the Cu 2 -N 6 /Gr defect supports high ORR molecule binding.
There are two possibilities of ORR pathways for Cu 2 -N 6 /Gr as shown in Fig. 1 and Table 2. Both pathways show that the O 2 adsorbent is high relative to the other structures considered in this study. For the structure of Cu 2 -N 6 /Gr, the 1-5-8 pathway has a higher O 2 adsorbent affinity than the 3-4-6-8 pathway. For both the 1-5-8 and 3-4-6-8 pathways, the initial step is the same, but having different binding energies may cause the ORR to take two different paths. Fig. 2

Effects of the Mulliken charge
The Mulliken charge for the adsorbate of the optimized structures is shown in Table 2 Table 2 indicate that the bonding capability is dependent on the difference of Q x and Q Cu and the bond length d Cu-O . Strong bonding is favored by a large difference in Q x and Q Cu and a short bond length d Cu-O .
Separate generation of the OH molecule may be due to the attraction of electrons by the Cu atoms. According to Mulliken charges, N atoms lose fewer electrons to Cu atoms of Cu 2 -N 6 /Gr on demand of electrons compared to Cu-N 4 /Gr. Therefore, the Cu atom of Cu 2 -N 6 /Gr may try to replenish its electron demand from O atom charges by forming a strong Cu-O bond (Fig. 5(b) and (c)) unlike Cu-N 4 /Gr which form H 2 O 2 . Furthermore, this scenario could be explained as π bonding electron on O atoms that attracts to each Cu Nallathambi et al. reported that H 2 O 2 generation could stop the ORR four-electron path halfway [36]. The ORR steps of an acidic medium show that a two-electron pathway is needed for the generation of H 2 O 2 . However, the optimization did not produce a two-electron pathway for the Cu 2 -N 6 /Gr structure. Cu 2 -N 8 /Gr structures only follow the ORR 02-04-05-08 path (Fig. 1). It is a four-electron pathway without H 2 O 2 generation. The *OOH intermediate on its ORR path is more favorable to form dual OH over H 2 O (Figs. 4(c) and 6(a)-(c)). Binding of the H atom (01) to the second O atom (Fig. 6(a)) causes the breakage of the O-O bond due to π electron deviation. This situation creates the H 2 O molecule short time. However, the formation of *2OH is made possible by the attraction of electrons from the O atom of H 2 O by the Cu atom caused by weak Cu-N bonding. This in turn would weaken one of the O-H bonds and lead to the release of H to react to form OH ( Fig. 6(b)). Since the *OOH of Cu 2 -N 8 /Gr structure creates the *2OH configuration by H atoms binding to both O atoms (Fig. 6(c)), the possibility of H 2 O 2 formation is reduced.
The ORR associated with the Cu-N 4 /Gr structure follows only the 02-04-06-08 path ( Fig. 1) with *O 2 , *OOH, *O, *OH, *H 2 O intermediates respectively (Fig. 7(a)-(e)). The optimized structures show the possibility of H 2 O 2 generation (Fig. 7(f)). With a binding energy of − 0.92 eV and a Cu-O bond length of 2.18 Å, H 2 O 2 can easily break free from the Cu-N 4 /Gr defect (Fig. 7(f)). The O-O bond in the simulation is shown to be stable, as suggested by Bhatt et al. [21]. The Mulliken charge population value of H 2 O 2 , 0.026 e, indicates that H 2 O 2 has the ability to break free from the catalyst surface as H 2 O. Since H 2 O 2 directly participates in the degeneration of the fuel cell membrane [37,38], the structure of Cu-N 4 /Gr can significantly reduce the efficiency of the fuel cell. Additionally, the two-electron pathway of H 2 O 2 generation can reduce the current density of the fuel cell.

HOMO and LUMO energy calculations
The HOMO and LUMO energies were calculated for every step of the ORR and structural defects to further investigate the binding energy. The energy gap (E g ), chemical hardness (η), chemical potential (μ), and electrophilicity index (ω) of *O 2 , *H 2 O, *H 2 O 2, intermediates and defect structures were calculated as shown in Table 3 using Eqs. (6)- (9) [29,31,32].
According to Eqs. (6) and (7) η and E g are directly associated with each other. The chemical hardness elucidates the stability and reactivity of the structure. Therefore, stable structures are less reactive, and high-reactive structures are less stable [8,29,31]. Defects can be ranked in decreasing order of chemical hardness as Cu 2 -N 8 /Gr > Cu 2 -N 6 /Gr > Cu-N 4 /Gr. This suggests that the reactivity of the Cu-N 4 /Gr structure is higher than that of other structures, but it is also less stable. This reactivity could cause electrons in the valence band to be easily excited to occupy the conduction band. The reactivity order is confirmed by the formation energy shown in Table 1. The Cu 2 -N 8 /Gr structure is more stable than others and is the least reactive of the three defect structures, and Cu 2 -N 6 /Gr has medium structural stability and reactivity according to chemical hardness values.
The chemical hardness of the oxygen bonding decreasing order Cu 2 -N 8 /Gr = Cu-N 4 /Gr > Cu 2 -N 6 /Gr indicates that Cu 2 -N 6 /Gr is less stable and more reactive than the other two structures. If the ORR first step is less stable and has high reactivity, the probability of  achieving the next ORR step is high. Therefore, the Cu 2 -N 6 /Gr defect is more suitable to complete the ORR first step. From previous calculations, the binding energy of the second step for ORR of Cu 2 -N 6 /Gr structure is higher than that of the first step, and this collaborates well with the chemical hardness data. The chemical hardness of the H 2 O bond in decreasing order is as follows: Cu 2 -N 8 / Gr > Cu-N 4 /Gr > Cu 2 -N 6 /Gr. This indicates that the H 2 O bonded to the Cu 2 -N 8 /Gr defects is more stable than the other defects. DFT calculations also show that the longest bond length of Cu-O and the highest binding energy would increase the probability of breaking H 2 O from the Cu 2 -N 8 /Gr defect sites. In the case of H 2 O bonding to the highly reactive Cu 2 -N 6 /Gr defect, this would result in an oxygen evolution reaction (OER) [20,39].
The chemical potential is directly related to Mulliken's electronegativity [31]. The chemical reactivity of the structure increases with decreasing chemical potential. The chemical potential decreases in the following order of defects: Cu 2 -N 8 /Gr > Cu 2 -N 6 /Gr = Cu-N 4 /Gr. This shows that the chemical reactivity of the Cu 2 -N 6 /Gr and Cu-N 4 /Gr structures is slightly higher than that of Cu 2 -N 8 /Gr. Furthermore, the chemical reactivity of the oxygen-bonded structures is the highest for Cu 2 -N 6 /Gr, confirming the previous chemical hardness evaluation of Cu 2 -N 6 /Gr.
Since electrophilicity is the electron accepting capacity from the surrounding atoms to form a stable energy state, and this relates directly to the structural stability. The electrophilicity of the studied defects has the following decreasing order, Cu-N 4 /Gr > Cu 2 -N 6 / Gr > Cu 2 -N 8 /Gr, and this also confirms with the reactivity order of the structures in terms of chemical hardness. The oxygen-bonded Cu 2 -N 6 /Gr shows a higher electron acceptation than the other structures, hence, indicating that the Cu 2 -N 6 /Gr structure favors electron adsorption to form a stable structure. Two ORR steps of Cu 2 -N 6 /Gr that generate different paths could occur to this reactivity. The Cu 2 -N 8 /Gr and Cu-N 4 /Gr structures show similar electrophilicity values for the O 2 bonding. The stability of the Cu-N 4 /Gr structure is attributed to the four nitrogen atoms that donate electrons to the Cu atom. Compared to the Cu 2 -N 6 /Gr structure, which has six nitrogen atoms, the Cu-N 4 /Gr structure is more stable ( Table 1). The Cu 2 -N 8 /Gr structure has twice as many nitrogen atoms as Cu-N 4 /Gr, but its electrophilicity value is likely to be low, similar to that of Cu-N 4 /Gr. Furthermore, H 2 O-bonded Cu 2 -N 8 /Gr showed the lowest electrophilicity value for the last ORR step, indicating a low probability that a further reaction will take place.

Free energy
Free energy values were calculated for the three defect structures with possible ORR pathways in an acidic medium as shown in Fig. 8(a-d). Cu 2 -N 6 /Gr shows a spontaneous reaction (downhill) for the first two steps of ORR with cell potential (U) 0-1.00 V (Fig. 8  (a) and (b)). The second to the third steps of the ORR for all potentials show uphill. From the third step to the last downhill, the configuration was found again for 0 V-1.00 V. This indicates that a large energy barrier is present between the second step and the third step, which prevents the spontaneous reaction of ORR from being completed. The highest reaction step barrier, 5.34 eV and 5.11 eV, are shown by Cu 2 -N 6 /Gr, for both (01) and (02) paths at U = 0 V, respectively, in *O-O hydrogenate by attaching the H atom to form *OOH and *O-OH. In summary, even with good binding energies, both ORR paths of Cu 2 -N 6 /Gr structure did not show spontaneity.
The Cu 2 -N 8 /Gr structure shows spontaneous reaction for all ORR steps for cell potential ranges from 0 to 0.28 V. With cell potential >0.28 V, the fifth step of ORR becomes uphill, showing the *2OH → *OH + H 2 O reaction process is not spontaneous (Fig. 8(c)). A further increase in the cell potential caused an increase in the uphill processes. The Cu-N 4 /Gr structure shows a spontaneous process for all ORR steps until the cell potential reaches 0.49 V. After the cell potential exceeds 0.49 V (U > 0.49 V), the ORR step of forming *O + H 2 O becomes uphill (Fig. 8(d)). Then, with the enhancement of the potential, the uphill process also increased. The maximum cell potentials for spontaneous reactions are shown to be 0.49 and 0.28 V for Cu-N 4 /Gr and Cu 2 -N 8 /Gr, respectively. These results show that the open circuit voltage (OCV) is relatively better than the Pt catalyst cathode, which initially shows 0.9 V, as reported by Kim et al. [40].
A comparison of the formation energies and chemical hardness indicates that the Cu-N 4 /Gr structure has a lower probability of forming and is less stable than other structures. Therefore, generating 0.49 V as cell voltage is a low possibility. Furthermore, the Cu-N 4 /Gr defect structure has a greater possibility of forming H 2 O 2 , as indicated in Table 2. The free energy diagram (Fig. 9) shows that the formation of H 2 O 2 was spontaneous within the range of 0-0.49 V. After exceeding 0.49 V, the H 2 O 2 formation step becomes uphill. Both the free energy and the maximum potential of the downhill process to form H 2 O 2 and H 2 O from *OOH are the same. Furthermore, the Mulliken charges of the O atoms of *OOH are approximately equal (− 0.312 e and − 0.331 e, respectively), indicating that the final step of ORR of Cu-N 4 /Gr is solely dependent on the binding of the H atom to the O atom of *OOH (Fig. 7(b)    However, this study has several limitations that should be acknowledged. First, the 3-21G basis set used in this study can be considered small. However, it was selected based on Bhatt et al.'s study [21], which compared Cu-N 2 /Gr results with the current study. Furthermore, this study did not include periodic boundary conditions or dispersion interactions, unlike the work by Xiao et al. [20] on Cu-N 4 /Gr. However, these studies cannot be used as a basis for a comparison with the current study due to their use of different DFT simulation software. The use of Gaussian09 instead of VASP-like simulation software was due to financial constraints. Although VASP-like software is better suited for surface simulation and PBC approach calculations, this study employed a single layer of graphene to examine the impact of defect structures on the surface of the catalyst in relation to ORR, similar to the approaches taken by Bhatt et al. [21] and Kattel et al. [16]. Future studies related to this research will employ dispersion correction with the PBC approach to provide a basis for comparison with the current study.

Conclusion
Computational calculations predicted the Cu 2 -N 8 /Gr and Cu-N 4 /Gr motifs, showing promising spontaneous ORR with maximum cell potential of 0.28 and 0.49 V under standard conditions. The Cu 2 -N 6 /Gr motif is not favorable for spontaneous ORR. From the formation energy calculations, Cu 2 -N 8 /Gr, Cu-N 4 /Gr, and Cu 2 -N 6 /Gr are stable for potential range (U) 0-5.87, 0-6.6, and 0-5.89 V, respectively. The maximum cell potentials of Cu 2 -N 8 /Gr and Cu-N 4 /Gr are within the potential range of formation energy, showing structural stability. Based on the chemical hardness, chemical potential, electrophilicity, and Mulliken charge calculations, the Cu 2 -N 8 /Gr shows a higher stability and favors ORR than Cu 2 -N 6 /Gr. Cu 2 -N 6 /Gr is more favorable for ORR than Cu-N 4 /Gr by binding energy values. The Cu 2 -N 6 /Gr motif is not recognized as a spontaneous catalyst by free energy calculations. Therefore, considering these three motifs, only Cu 2 -N 8 /Gr and Cu-N 4 /Gr show promising ORR ability. Consequently, the possibility of cell voltage solely depends on Cu 2 -N 8 /Gr could occur due to its high stability. According to DFT optimization, H 2 O 2 is generated only from Cu-N 4 /Gr. The other structures follow four-electron transfer pathways to generate high electron density. The Cu-Nitrogen doped non-PGM catalyst with the considered motifs has a relatively high current density even with low cell voltage when comparing to Pt-like catalysts. The amount of H 2 O 2 generated by Cu-N 4 /Gr could be limited by its low stability. Bhatt et al. have also shown that the Cu-N 2 /Gr motif has a high probability of generating H 2 O 2, but the − 5.68 eV formation energy [21] indicates that it has lower stability than Cu-N 4 /Gr. In this study, a novel Cu 2 -N 8 /Gr catalyst is a good candidate for fuel cell application as it has good cell potential and non-H 2 O 2 formation. Furthermore, the three structures presented in this study show a better O 2 binding energy than the Pt catalyst and PGM alloys [41].

Author contribution statement
Yashas Balasooriya: Conceived and designed the experiments; Performed the experiments; Analyzed and interpreted the data; Wrote the paper.
Pubudu Samarasekara, Muhammad Raziq Rahimi Kooh: Conceived and designed the experiments; Contributed reagents, materials, analysis tools or data; Wrote the paper.
Chee Ming Lim, Yuan-Fong Chou Chau: Contributed reagents, materials, analysis tools or data; Wrote the paper. Roshan Thotagamuge: Conceived and designed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper.

Data availability statement
Data included in article/supplementary material/referenced in article.

Funding
The work described in this article was supported by the Universiti Brunei Darussalam Research Grant UBD/RSCH/1.9/FICBF(b)/ 2022/017.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.