Density Functional Theory Investigation of Temperature-Dependent Properties of Cu-Nitrogen-Doped Graphene as a Cathode Material in Fuel Cell Applications

In this study, density functional theory (DFT) was used to investigate the influence of temperature on the performance of a novel Cu-nitrogen-doped graphene Cu2-N8/Gr nanocomposite as a catalyst for the oxygen reduction reaction (ORR) in fuel cell applications. Our DFT calculations, conducted using Gaussian 09w with the 3–21G/B3LYP basis set, focus on the Cu-nitrogen-doped graphene nanocomposite cathode catalyst, exploring its behavior at three distinct temperatures: 298.15 K, 353.15 K, and 393.15 K, under acidic conditions. Our analysis of formation energies indicates that the structural stability of the catalyst remains unaffected as the temperature varies within the potential range of 0–7.21 V. Notably, the stability of the ORR steps experiences a marginal decrease with increasing temperature, with the exception of the intermediate OH + H2O (*OH + H + *OH). Interestingly, the optimization reveals the absence of single OH and H2O intermediates during the reactions. Furthermore, the OH + H2O step is optimized to form the OH + H + OH intermediate, featuring the sharing of a hydrogen atom between dual OH intermediates. Free energy calculations elucidate that the catalyst supports spontaneous ORR at all temperatures. The highest recorded maximum cell potential, 0.69 V, is observed at 393.15 K, while the lowest, 0.61 V, is recorded at 353.15 K. In particular, the Cu2-N8/Gr catalyst structure demonstrates a reduced favorability for the H2O2 generation at all temperatures, resulting in the formation of dual OH intermediates rather than H2O2. In conclusion, at 393.15 K, Cu2-N8/Gr exhibits enhanced catalyst performance compared to 353.15 K and 298.15 K, making it a promising candidate for ORR catalysis in fuel cell applications.


Introduction
Fuel cells operating at low temperatures have attracted considerable attention due to their enhanced durability, reliability, and efficiency with low weight, making them more portable compared to their counterparts.These cells typically operate at temperatures below 473.15 K, often around 373.15 K, and commonly rely on carbon-supported catalysts to enhance the oxygen reduction reaction (ORR).However, these carbon supports are susceptible to rapid corrosion at higher temperatures, resulting in catalyst agglomeration, material durability, and subsequent reductions in efficiency [1].
Walch et al. revealed the effectiveness of the DFT method in studying ORR performance and catalyst properties [2].Commercially available low-temperature fuel cells predominantly employ platinum group metals (PGMs), such as platinum (Pt), palladium (Pd), and iridium (Ir), as catalysts.These costly materials effectively facilitate hydrogen ionization at the anode [3,4] and provide an efficient ORR surface at the cathode [5].However, as the carbon support undergoes oxidation, the catalyst layered upon it tends to agglomerate, resulting in reduced fuel cell efficiency over short operational cycles.Catalyst materials made of transition-metal-based perovskite oxides have also been used as effective catalysts compared to PGM.However, the process of creating those materials is relatively costly due to the need for La, Sr-like alkaline-earth/rare-earth materials, and transition metals [6].To address this issue and lower the costs, non-platinum group metals (non-PGM) have been explored as potential catalysts.These alternatives, typically composed of nitrogen atoms doped onto graphite layers [7][8][9] or transition metal layers [10], offer advantages such as ease of fabrication, high electron conductivity, and cost-effectiveness [11][12][13].
Low-temperature fuel cells, such as proton exchange membrane (PEM) fuel cells, have shown improved efficiency with increasing temperature [14].A temperature range of 353.15-403.15K has been identified as particularly promising for low-temperature fuel cells [15].Some studies have reported significant performance enhancements in PEM fuel cells at a temperature of 393.15K [15,16].Raising the temperature has the benefits of increased voltage and efficiency, but it also introduces the drawback of reduced durability [14,16,17].
The output voltage, as per the Nernst equation, is directly proportional to the temperature, with a higher temperature resulting in enhanced charge particle kinetics and an increased fuel cell voltage [18].However, the lowest performance of the fuel cell is observed at around 343.15 K [19].As Guan et al. explained, the reason for the lowest performance at 343.15 K and the highest performance at 393.15 K may be a phase change that occurs with increasing temperature.Furthermore, as temperature increases, the oxygen evolution reaction (OER) takes place, leading to the reconstruction of the catalyst surface [20].To achieve a voltage increase, the enhancement in voltage must exceed the voltage loss due to the negative thermodynamic correlation between the open-circuit voltage and temperature [18,21].
The rate of catalyst decay is influenced by several factors, including temperature, ORR, hydrogen oxidation reaction (HOR), high potential, and the pH of the medium [14].The formation of carbon monoxide (CO) at the cathode is due to the reaction between water and the carbon support of the catalyst that can cover the catalyst, and this significantly reduces the electrochemical reaction [22].This process takes place mainly at temperatures below 373.15K [23].Above 373.15 K, the water generated from the reaction turns into vapor and is rapidly removed.This dehydration process causes a reduction in CO generation from the catalyst [24].
Numerous non-PGM catalysts, such as Fe/N/Gr, Co/N/Gr, and Mn/N/Gr, have been extensively studied, and some have even reached commercialization [25,26].For instance, Mn/N/Gr catalysts have been demonstrated to be unstable for potential above 0.53 V [27], while Co/N/Gr-based catalysts tend to produce hydrogen peroxide (H 2 O 2 ) during ORR, which reduces the current density and can damage the fuel cell membrane [28].In contrast, Fe-based non-PGM catalysts have been shown to be promising, especially with dual metal structures such as Fe 2 -N 6 /Gr and Fe 2 -N 8 /Gr [29,30].Increasing the Fe content to 1.0 wt % can significantly reduce H 2 O 2 generation [31].However, a drawback of Febased catalysts is the formation of Fe(OH) 2 during ORR, which can lead to the removal of Fe atoms from the catalyst layer [32].Furthermore, Fe oxidation can contribute to the catalyst layer degradation [33].
Previous research on copper-based non-PGM single metal catalysts has shown them to have good spontaneous ORR [34,35].In our previous work, we performed a DFT investigation on Cu-based novel duel metal non-PGM catalyst structures (Cu 2 -N 6 /Gr and Cu 2 -N 8 /Gr) and a single metal structure (Cu-N 4 /Gr) [36], and our study showed that all three structures exhibit higher stability than the previously studied Cu-N 2 /Gr structure.The study also showed that the possibility of H 2 O 2 formation is much less with dual metal structures compared to single metal structures.It is noted that the study was carried out at an initial temperature of 298.15K and pressure of 1 atm without applying van der Waals force correction.In conclusion, the Cu 2 -N 8 /Gr structure showed promising results, confirming its potential for high cell potential with spontaneous ORR and the absence of H 2 O 2 formation.Additionally, Cu-N 4 /Gr was also shown to have a high cell potential with spontaneous ORR but has a higher probability of forming H 2 O 2 .
In this study, DFT calculations were performed to determine the effect of temperature on the Cu 2 -N 8 /Gr catalyst in an acidic medium.To avoid the limitations of previous studies, a dispersion correction was applied to the DFT calculation.[14,15,19].Therefore, for this study, 298.15, 353.15, and 393.15K temperature values were selected to investigate the impact of temperature on the fuel cell performance of Cu 2 -N 8 /Gr as a cathode catalyst in an acidic medium.At these three temperature values, the structural stability, binding ability of ORR steps, the break-free ability of H 2 O, evidence of H 2 O 2 generation, and maximum cell potential were predicted as assessments of objectives.

Stability and Reactivity
Stability and reactivity were predicted using formation energy data [28,35], and the HOMO-LUMO energy calculation [37][38][39][40] is elaborated in Equation ( 6) and Equations ( 15)- (18), respectively.Under open circuit conditions, the Cu 2 -N 8 /Gr structure provided formation energy values of −28.856, −28.865, and −28.864 eV for the three different temperatures (298.15,353.15, and 393.15 K) at zero potential (U = 0), as presented in Table 1.According to Equation (6), an increase in the potential value leads to an enhancement of the formation energy, resulting in a decrease in the structure's stability.Minor changes in formation energy, at the same potential, generated due to thermal correction were small in comparison to the optimized energy.Negative values in the formation energies indicate that the energy gained favors structural stability.Notably, the formation energy of the structure at all the selected temperatures remained negative for external potentials ranging between 0 and 7.21 V.According to free energy data, the highest maximum cell potential of 0.69 V is observed under an external potential of 7.21 V, which favors this situation.Table 1 illustrates the results showing that the catalyst structure remains stable at all three temperatures as long as the external potential values do not exceed 7.21 V.
In the case of the dual metal structure, Cu 2 -N 8 /Gr involves multiple copper atoms that coordinate with eight nitrogen atoms.This configuration results in a significantly more extensive bonding network.It enhances coordination stability when compared to single metal structures, such as Cu-N 4 /Gr and Cu-N 2 /Gr, in which only one metal atom coordinates with a smaller number of nitrogen atoms.This increased stability is compatible with findings reported by Yang et al. [41].Additionally, the dual metal structure is likely to facilitate a more efficient distribution and delocalization of electrons due to the presence of multiple Cu and N atoms.The electrostatic potential (ESP) diagram (Figure 1), which supports the electron distribution theory, for all three temperatures, shows no significant differences in electronic structure.The red area represents the distribution of negative charge (electrons), and the blue and white areas represent positive charge and depleted charge density, respectively.In the case of the dual metal structure, Cu2-N8/Gr involves multiple copper atoms that coordinate with eight nitrogen atoms.This configuration results in a significantly more extensive bonding network.It enhances coordination stability when compared to single metal structures, such as Cu-N4/Gr and Cu-N2/Gr, in which only one metal atom coordinates with a smaller number of nitrogen atoms.This increased stability is compatible with findings reported by Yang et al. [41].Additionally, the dual metal structure is likely to facilitate a more efficient distribution and delocalization of electrons due to the presence of multiple Cu and N atoms.The electrostatic potential (ESP) diagram (Figure 1), which supports the electron distribution theory, for all three temperatures, shows no significant differences in electronic structure.The red area represents the distribution of negative charge (electrons), and the blue and white areas represent positive charge and depleted charge density, respectively.HOMO-LUMO energy calculations, as explained in Equations ( 15)- (18), were carried out to investigate the stability and reactivity of the Cu2-N8/Gr and *O2, *H2O, and OH + H2O intermediates [37][38][39][40].Specifically, the energy gap (Eg), chemical hardness (η), chemical potential (µ), and electrophilicity index (ω) were calculated for these structures.Here, the *O2 and 2*H2O intermediates were considered for the initial and final ORR steps, respectively.These intermediates are important when compared to other ORR steps due to the breaking of the O2 bond with the catalyst that initiates the ORR, and 2H2O is released from the catalyst to conclude the ORR.The OH + H2O intermediate was also considered due to its abnormality.Within all the temperature ranges considered, the HOMO-LUMO calculations remained constant, except for the OH + H2O intermediate, indicating that the catalyst structure did not change for the given temperatures.
Chemical hardness is a measure that reflects the stability and reactivity of a structure, where stable structures tend to be less reactive, and highly reactive structures are typically less stable [9,38,42].Table 2 shows that the OH + H2O intermediate's stability increases in the following temperature order: 353.15K < 298.15K < 393.15 K.This indicates that at 393.15 K, the intermediate becomes more stable, while at 353.15 K it becomes highly reactive compared to the other temperatures.The ESP diagram in Figure 2 reveals that as the temperature increases from 298.15 to 353.15 K, the electron charge density decreases HOMO-LUMO energy calculations, as explained in Equations ( 15)-( 18), were carried out to investigate the stability and reactivity of the Cu 2 -N 8 /Gr and *O 2 , *H 2 O, and OH + H 2 O intermediates [37][38][39][40].Specifically, the energy gap (E g ), chemical hardness (η), chemical potential (µ), and electrophilicity index (ω) were calculated for these structures.Here, the *O 2 and 2*H 2 O intermediates were considered for the initial and final ORR steps, respectively.These intermediates are important when compared to other ORR steps due to the breaking of the O 2 bond with the catalyst that initiates the ORR, and 2H 2 O is released from the catalyst to conclude the ORR.The OH + H 2 O intermediate was also considered due to its abnormality.Within all the temperature ranges considered, the HOMO-LUMO calculations remained constant, except for the OH + H 2 O intermediate, indicating that the catalyst structure did not change for the given temperatures.
Chemical hardness is a measure that reflects the stability and reactivity of a structure, where stable structures tend to be less reactive, and highly reactive structures are typically less stable [9,38,42].2 reveals that as the temperature increases from 298.15 to 353.15 K, the electron charge density decreases around the Cu-N structure.The pointed green area in Figure 2b indicates the potential depletion of the electron density.This causes the lower stability of the OH + H 2 O intermediate at 353.15 K.However, when the temperature reaches 393.15 K, this depletion area is significantly reduced.Furthermore, at 393.15 K, the electron distribution becomes more symmetric (Figure 2c) around the Cu-N structure, as indicated by the similarly sized green area on the opposite side of the blue area.The high stability of the OH + H 2 O intermediate at 393.15 K is attributed to this symmetric electron charge distribution.According to Watanabe et al., vibrational motion can significantly change the electron density [43], and as the temperature increases, the vibration frequency and direction of the atoms change rapidly.This is attributed to variations in the ESP at different temperatures.To confirm the vibration effects, the vibration modes at different frequencies were examined at a temperature of 353.15 K.A frequency around 289.32 cm −1 indicates the displacement of H atoms, as shown in Figure 3.This mode also has a significantly high infrared value of 25.30 compared to other frequencies.Other vibrational modes were considered but then neglected because their vibration direction did not align with the charge distribution in the ESP diagram.
Here, the H atoms represent the most positive charge in the ESP diagram.Therefore, the long displacement vector of the H atom (left H atom) compared to the abnormal electron depletion area is highlighted in Figure 2b.The displacement vector of the hydrogen atom located on the right-hand side in Figure 3 exhibits a smaller magnitude in comparison to the hydrogen atom situated on the left-hand side, an observation consistent with the electrostatic potential (ESP) depicted in Figure 2b.The middle H atom, with the smallest displacement vector, creates a bottleneck in the dumbbell-shaped electron-depletion area in the ESP diagram.around the Cu-N structure.The pointed green area in Figure 2b indicates the potential depletion of the electron density.This causes the lower stability of the OH + H2O intermediate at 353.15 K.However, when the temperature reaches 393.15 K, this depletion area is significantly reduced.Furthermore, at 393.15 K, the electron distribution becomes more symmetric (Figure 2c) around the Cu-N structure, as indicated by the similarly sized green area on the opposite side of the blue area.The high stability of the OH + H2O intermediate at 393.15 K is attributed to this symmetric electron charge distribution.According to Watanabe et al., vibrational motion can significantly change the electron density [43], and as the temperature increases, the vibration frequency and direction of the atoms change rapidly.This is attributed to variations in the ESP at different temperatures.To confirm the vibration effects, the vibration modes at different frequencies were examined at a temperature of 353.15 K.A frequency around 289.32 cm −1 indicates the displacement of H atoms, as shown in Figure 3.This mode also has a significantly high infrared value of 25.30 compared to other frequencies.Other vibrational modes were considered but then neglected because their vibration direction did not align with the charge distribution in the ESP diagram.Here, the H atoms represent the most positive charge in the ESP diagram.Therefore, the long displacement vector of the H atom (left H atom) compared to the abnormal electron depletion area is highlighted in Figure 2b.The displacement vector of the hydrogen atom located on the right-hand side in Figure 3 exhibits a smaller magnitude in comparison to the hydrogen atom situated on the left-hand side, an observation consistent with the electrostatic potential (ESP) depicted in Figure 2b.The middle H atom, with the smallest displacement vector, creates a bottleneck in the dumbbell-shaped electron-depletion area in the ESP diagram.Furthermore, at 353.15 K, certain peaks of the OH + H2O intermediate exhibit a shift towards lower wavenumbers within the range of 1200-1600 cm −1 .This observation implies thermal denaturation, as explained by Panick et al. [44], or a phase change associated with temperature, as elucidated by Guan et al. [20].These findings contribute valuable insights into the temperature-dependent dynamics of the studied OH + H2O intermediate.

ORR Steps and Binding Energy
In this study, the ORR steps did not indicate deviations at different temperatures.However, when the intermediate 2OH is bonded with an H atom to proceed to the next step, some abnormalities are observed.The expected outcome after *2OH − + H + was *OH with a separated H2O molecule.However, after optimization, the *OH + H + *OH intermediate was formed, where both OH molecules shared a single H atom (Figure 5).With the exception of *OH + H + *OH (intermediates after *2OH − + H + ), the other intermediates did not indicate any deviations in terms of bond lengths, angles, and binding strength with different temperatures (Table 3).

ORR Steps and Binding Energy
In this study, the ORR steps did not indicate deviations at different temperatures.However, when the intermediate 2OH is bonded with an H atom to proceed to the next step, some abnormalities are observed.The expected outcome after *2OH − + H + was *OH with a separated H 2 O molecule.However, after optimization, the *OH + H + *OH intermediate was formed, where both OH molecules shared a single H atom (Figure 5).With the exception of *OH + H + *OH (intermediates after *2OH − + H + ), the other intermediates did not indicate any deviations in terms of bond lengths, angles, and binding strength with different temperatures (Table 3).In the first step of ORR, the oxygen molecule exhibits a significant binding strength value compared with that reported by Kattel et al. and Bhatt et al. [28,35], and it is slightly mitigated with increasing temperatures.This trend has been observed for all other intermediates, as recorded in Table 3, except for the OH + H2O intermediate (*OH + H + *OH).Initially, the *OH + H + *OH intermediate shows a reduction in its binding strength that mirrors the behavior of the other intermediates.However, at 393.15 K, its binding strength increases compared to the lower temperature of 353.15 K.In the first step of ORR, the oxygen molecule exhibits a significant binding strength value compared with that reported by Kattel et al. and Bhatt et al. [28,35], and it is slightly mitigated with increasing temperatures.This trend has been observed for all other intermediates, as recorded in Table 3, except for the OH + H 2 O intermediate (*OH + H + *OH).Initially, the *OH + H + *OH intermediate shows a reduction in its binding strength that mirrors the behavior of the other intermediates.However, at 393.15 K, its binding strength increases compared to the lower temperature of 353.15 K.
Binding energy calculations reveal that the highest binding energy for the OH intermediate was at −6.45 eV, which is observed only at an initial temperature of 298.15 K.At other temperatures, the binding energy significantly decreases to −0.25 eV at 353.15 K and −0.24 eV at 393.15 K.These binding energy values at higher temperatures are not favorable for the ORR, as they represent weaker binding strength for all intermediates.In a typical case, the OH intermediate is generated after the water molecules are rendered and moved away, and subsequently, an H atom attaches itself to the O intermediate, as shown in the ORR reaction Equations ( 11)-( 13).The optimized reactions show that the H atom is attached to the 2OH intermediate, instead of forming a water molecule and OH intermediate.This results in the formation of the OH + H + OH intermediate, as previously explained.Therefore, the weakest bond strength of the OH intermediate is not considered a step in the ORR.
In the final step of the ORR, it is expected that the water molecule will form and easily break away from the catalyst.The calculations confirmed that H 2 O has the lowest binding strength, which favors the described situation.Furthermore, the longest bond length of As the OOH intermediate emerges during the ORR (reaction Equation ( 2)), there is a significant possibility of H 2 O 2 generation in the middle of the reaction.To investigate this, the H atom is introduced to the first oxygen atom to determine whether H 2 O 2 is formed.Figure 6 illustrates this reaction at an initial temperature of 298.15 K.As reported in our previous study [36], under the initial temperature and the other two temperatures, the formation of the 2OH intermediate occurs instead of H 2 O 2 molecules.This phenomenon can be attributed to the electronegativity of nitrogen (N) being higher than that of copper (Cu).Binding energy calculations reveal that the highest binding energy for the OH intermediate was at −6.45 eV, which is observed only at an initial temperature of 298.15 K.At other temperatures, the binding energy significantly decreases to −0.25 eV at 353.15 K and −0.24 eV at 393.15 K.These binding energy values at higher temperatures are not favorable for the ORR, as they represent weaker binding strength for all intermediates.In a typical case, the OH intermediate is generated after the water molecules are rendered and moved away, and subsequently, an H atom attaches itself to the O intermediate, as shown in the ORR reaction Equations ( 11)-( 13).The optimized reactions show that the H atom is attached to the 2OH intermediate, instead of forming a water molecule and OH intermediate.This results in the formation of the OH + H + OH intermediate, as previously explained.Therefore, the weakest bond strength of the OH intermediate is not considered a step in the ORR.
In the final step of the ORR, it is expected that the water molecule will form and easily break away from the catalyst.The calculations confirmed that H2O has the lowest binding strength, which favors the described situation.Furthermore, the longest bond length of As the OOH intermediate emerges during the ORR (reaction Equation ( 2)), there is a significant possibility of H2O2 generation in the middle of the reaction.To investigate this, the H atom is introduced to the first oxygen atom to determine whether H2O2 is formed.Figure 6 illustrates this reaction at an initial temperature of 298.15 K.As reported in our previous study [36], under the initial temperature and the other two temperatures, the formation of the 2OH intermediate occurs instead of H2O2 molecules.This phenomenon can be attributed to the electronegativity of nitrogen (N) being higher than that of copper (Cu).In the Cu2-N8 structure, where one Cu atom is attached to four N atoms, this arrangement is responsible for the electron cloud depletion around the Cu atoms.According to the studies by Xiao et al. [34] and Noh et al. [45], the H atom could bind to the terminal O atom (Figure 7a) of the OOH intermediate to form an H 2 O molecule.In this study, the possibility of this reaction was determined via the optimization process.For all three temperatures, the final optimization resulted in the formation of a dual OH intermediate (Figures 7c and 8c) instead of the formation of the water molecule.Figure 7 depicts the reaction at a temperature of 298.15 K.When the introduced H atom attempts to bond with the second O atom of OOH to form H 2 O, the previously bonded H atom breaks its bond with the second O atom and moves toward the first O atom (Figure 7b).During this process, a water molecule forms in the initial optimization steps, and then the Cu atom near the second O atom attracts it to establish a strong bond.This phenomenon is the same as that explained in previous H dual OH molecules are formed instead of H2O2 (white, grey, blue, orange, and red spheres represent hydrogen, carbon, nitrogen, copper, and oxygen atoms, respectively).
According to the studies by Xiao et al. [34] and Noh et al. [45], the H atom could bind to the terminal O atom (Figure 7a) of the OOH intermediate to form an H2O molecule.In this study, the possibility of this reaction was determined via the optimization process.For all three temperatures, the final optimization resulted in the formation of a dual OH intermediate (Figures 7c and 8c) instead of the formation of the water molecule.Figure 7 depicts the reaction at a temperature of 298.15 K.When the introduced H atom attempts to bond with the second O atom of OOH to form H2O, the previously bonded H atom breaks its bond with the second O atom and moves toward the first O atom (Figure 7b).During this process, a water molecule forms in the initial optimization steps, and then the Cu atom near the second O atom attracts it to establish a strong bond.This phenomenon is the same as that explained in previous H2O2 generation by the OOH intermediate.Due to electron depletion around the terminal O of *OOH, the O atom could not bear two H atoms as H2O.This could cause the initially bonded H atom to break its bond with the second O atom.The changes observed in the 2OH and OH + H + HO intermediates, as opposed to O + H2O and OH +H2O in the ORR, may be attributed to the dual metal catalyst structure, as explained before.The following reaction Equations ( 1)-( 5) depict the optimized ORR steps from the current study in an acidic medium.The asterisk (*) indicates the defect's configuration.

Reaction Spontaneity
All free energy calculations were performed for the ground state of the ORR steps, as shown in Figure 9.After conducting free energy calculations for all three temperatures, unlike in our previous studies [36], it was found that the maximum cell potential was greater than 0.60 V.At the initial temperature of 298.15 K, all reaction coordinates exhibited a downhill process (spontaneous process) until the cell potential reached 0.66 V (maximum cell potential).Beyond the maximum cell potential, the last step of the ORR became an uphill process (Figure 9a), indicating an increase in the energy barrier between the reaction coordinates.Further enhancement of cell potential is attributed to the initiation of the uphill process from the 2OH intermediate step.The difference in maximum cell potentials between the previous and current findings may be attributed to the inclusion of a dispersion correction for van der Waals forces in the Gaussian calculations of the current study.

Reaction Spontaneity
All free energy calculations were performed for the ground state of the ORR steps, as shown in Figure 9.After conducting free energy calculations for all three temperatures, unlike in our previous studies [36], it was found that the maximum cell potential was greater than 0.60 V.At the initial temperature of 298.15 K, all reaction coordinates exhibited a downhill process (spontaneous process) until the cell potential reached 0.66 V (maximum cell potential).Beyond the maximum cell potential, the last step of the ORR became an uphill process (Figure 9a), indicating an increase in the energy barrier between the reaction coordinates.Further enhancement of cell potential is attributed to the initiation of the uphill process from the 2OH intermediate step.The difference in maximum cell potentials between the previous and current findings may be attributed to the inclusion of a dispersion correction for van der Waals forces in the Gaussian calculations of the current study.By increasing the temperature to 353.15 K, the maximum cell potential decreased to 0.61 V (Figure 9b).Beyond 0.61 V, the chemical reactions become uphill processes from the 2OH intermediate step to the OH + H 2 O intermediate step.This adverse trend intensifies with higher temperatures.When the temperature was raised to 393.15 K, the maximum cell potential reached the highest value of 0.69 V (Figure 9c).For cell potentials higher than 0.69 V, the OH + H 2 O intermediate step acted as a barrier to the downhill process of the ORR.According to the free energy data, further increases in the cell potential augmented the uphill process between the 2OH and OH + H 2 O intermediates.The order of maximum cell potential at different temperatures was as follows: 353.15K < 298.15K < 393.15 K, indicating that temperature and maximum cell potential were not directly proportional.The decrease in maximum cell voltage at around 353.15 K suggests a decline in cell performance, which aligns with the findings of Ferraris et al. [19], who reported the lowest cell potential around the operating temperature of 343.15 K, supporting our observation.Furthermore, other investigations [15,16] suggest that at 393.15 K PEM fuel cell performance improved dramatically in congruence with cell voltage escalation due to high temperature.By increasing the temperature to 353.15 K, the maximum cell potential decreased t 0.61 V (Figure 9b).Beyond 0.61 V, the chemical reactions become uphill processes from the 2OH intermediate step to the OH + H2O intermediate step.This adverse trend inten sifies with higher temperatures.When the temperature was raised to 393.15 K, the max mum cell potential reached the highest value of 0.69 V (Figure 9c).For cell potentia higher than 0.69 V, the OH + H2O intermediate step acted as a barrier to the downhi process of the ORR.According to the free energy data, further increases in the cell poten tial augmented the uphill process between the 2OH and OH + H2O intermediates.Th order of maximum cell potential at different temperatures was as follows: 353.15K 298.15K < 393.15 K, indicating that temperature and maximum cell potential were no directly proportional.The decrease in maximum cell voltage at around 353.15 K sugges a decline in cell performance, which aligns with the findings of Ferraris et al. [19], wh reported the lowest cell potential around the operating temperature of 343.15 K, suppor ing our observation.Furthermore, other investigations [15,16] suggest that at 393.15 PEM fuel cell performance improved dramatically in congruence with cell voltage escal tion due to high temperature.
In the free energy calculations, the observed up-down behavior of the maximum ce potential with temperature is mainly caused by OH + H2O intermediate step.The ES diagram of the OH + H2O intermediate (Figure 2c) delineates that at 393.15 K, this inte In the free energy calculations, the observed up-down behavior of the maximum cell potential with temperature is mainly caused by OH + H 2 O intermediate step.The ESP diagram of the OH + H 2 O intermediate (Figure 2c) delineates that at 393.15 K, this intermediate gains high stability due to a dense and symmetric electron density distribution.This stability, along with higher electron mobility at higher temperatures, could contribute to a higher maximum cell potential observed at 395.15K compared to 353.15 K, a conclusion supported by HOMO-LUMO energy calculations [23,24].
The free energy analysis of the first and last steps in the ORR revealed a 15 eV gap, showing a departure from the typical values observed in conventional ORR catalysts.In a work by Kattel et al. [28], the free energy of O 2 was determined through the reaction O 2 + 2H 2 → 2H 2 O, resulting in a free energy change of 4.92 eV.This value contrasts with most other density functional theory (DFT) calculations, where such values typically fall below 8 eV [27,46].While a high free energy gap might suggest reduced effectiveness of the catalyst, insights from Santisouk et al. [47] challenge this notion.Binding energy, bond length, and HOMO-LUMO energy calculations indicate the effectiveness of the catalyst with intermediates.Notably, the introduction of dispersion correction (GD3BJ) in the current study contributes significantly to this observed change.In contrast, a previous study [36] of ours in the same research field did not exhibit such a substantial free energy gap, with the O 2 to 2H 2 O free energy change hovering around 6 eV.This prompts consid-eration of alternative dispersion correction methods that could potentially enhance the accuracy of the free energy results.
Furthermore, given the novelty of the catalyst's structure, additional studies are warranted to validate its effectiveness.It is worth noting that some ORR free energy calculations in the literature indicate differences between the first and last ORR steps ranging from 8 to 14 eV [34,47,48] emphasizing the variability in observed outcomes across different studies.
Table 4 presents the binding energy, formation energy, and bond length of molecules on the catalyst surface at the initial temperature for various catalyst types.All the catalysts under investigation are non-platinum group metal (non-PGM) types, featuring transition metals and nitrogen doping.In terms of formation energy, our current study yields the most promising results, with the cobalt (Co)-based catalyst displaying the least promising outcome.Specifically, the titanium (Ti)-based catalyst exhibits the strongest binding capability for O 2 and H 2 O intermediates, suggesting its limited effectiveness as a catalyst material.Conversely, the Co-based catalyst demonstrates the weakest binding capability for the O 2 intermediate, indicating its diminished efficacy as a catalyst.In the context of H 2 O 2 formation at the conclusion of the ORR, our current study reveals no H 2 O 2 formation, whereas all other catalysts exhibit such potential.Notably, the Co-based catalyst, characterized by its weakest binding strength and longest bond length, implies the likelihood of H 2 O 2 formation and easy detachment from the catalyst layer.Furthermore, all dual molecules in our current study exhibit dual bonding, contributing to their strong attachment to the catalyst.These findings provide valuable insights into the diverse binding characteristics of the catalysts, shedding light on their potential efficacy in the context of the ORR.
Nonetheless, it is important to acknowledge several limitations of this study.Firstly, the use of a small 3-21G basis set may be considered a limitation.However, this choice was made to compare the results with previous investigations by Bhatt et al. [35], as well as to maintain consistency with our earlier study [36] and to compare with the Cu-N 2 /Gr results.Additionally, insufficient computational power meant that the study could not use a higher basis set.Furthermore, this study did not incorporate periodic boundary conditions (PBC), which distinguishes it from the approach employed by Xiao et al. [34].The utilization of different DFT simulation software in these referenced studies precludes direct comparison with our findings.It should also be noted that the preference for Gaussian software over VASP or ab initio molecular dynamics simulation software was influenced by financial constraints and a lack of specialists.In this study, a single layer of graphene was employed to examine the impact of the defect structures on the surface of the catalyst in the context of ORR, a methodology similar to those of Bhatt et al. [35] and Kattel et al. [28].In future studies related to this research, the use of the PBC approach will be explored to facilitate more direct comparisons with the current study and potentially expand the scope of investigation.

Methodology
The DFT calculations were performed using the Gaussian 09w software, applying the B3LYP method with a 3-21G basis set.The GD3BJ van der Waals correction (dispersion correction) was applied for every calculation.The form of the defect structures was visualized with GaussView 6.0.All atoms in each structure were relaxed using an optimization function in Gaussview 6 with the previously stated basis set and method.A single pristine graphene layer was created as the base, and from it, the Cu 2 -N 8 /Gr structure was developed to investigate the ORR under various temperatures (298.15,353.15, and 393.15 K).The pristine graphene layer design at the initial temperature has C-C and C-H bond lengths of 1.43 and 1.08 Å, respectively, and a C-C-C angle of 119.92°which is compatible with the studies by Bhatt et al. and Hernandez et al. [35,49].The structure that features the pyridinic nitrogen was constructed.Following the initial optimization, structures with a 1-unit imaginary frequency (=1) were reoptimized by using a stability function in GaussView 6.Individual optimization of the molecules involved in the ORR steps was carried out, with subsequent measurement and comparison of bond lengths and bond angles with previous studies by Xiao et al., Bhatt et al., on this subject.
The stability of the catalyst structure with the cathode potential at different temperatures was confirmed.Stability was inspected using the formation energy (∆E) defined as [28,35], Here, E graphene+(Cu 2 −N 8 ) is the energy of the optimized Cu 2 -N 8 graphene layer with the thermal correction to energy.µ C and µ N are the chemical potential of carbon, which is defined as the total energy per carbon atom for defect-free graphene, and the chemical potential of nitrogen, defined as half of the total energy of the N 2 molecule with the thermal correction to energy, respectively.x represents the number of nitrogen atoms added, and y represents the carbon atoms removed during the defect formation.M is the metal atom (Cu) doped in the system.E graphene is the energy of the optimized pristine graphene layer with the thermal correction to energy.E M is the total energy of M n+ defined as where E(M) is the total energy with the thermal correction to energy of isolated M (M = Cu) in the gas phase and n, e, and U are the number of electron transfer (+2), electron charge, and the external potential, respectively.Binding energies (BEs) were calculated for the ORR steps in an acidic medium, as defined by [28,34,35], Here, E de f ect+molecule is the total energy with the thermal correction to the energy of molecules adsorbed by the defective graphene.E de f ect is the total energy with the thermal correction to energy of the defective graphene configuration, and The free energies were calculated for each ORR step of all defect configurations as defined by [28], ∆G = ∆E + ∆ZPE − T∆S + ∆G U + ∆G pH + ∆G f ield (14) ∆E represents the energy calculated from DFT related to the relevant reaction step with the thermal correction to energies, while ∆ZPE is the correction obtained from DFT calculations, accounting for the zero-point energy.The value of entropy (S) was obtained from the DFT calculations for the three temperatures (T = 298.15,353.15, and 393.15 K) as a variable in this study, depending on temperature (T).In the context of this research, ∆G U = −eU where U and e are the electrode potential and the charge transferred, respectively.Additionally, ∆G pH = k B T × ln10 × pH, where k B is the Boltzmann's constant, and T is the temperature.Since this study considered acidic media, pH = 0. Therefore, normally the ∆G pH value is considered to be zero.∆G f ield is the 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 Nørskov et al. [50].Corrections of zero-point energy and entropy values were obtained after frequency calculations.Free energy vs. reaction coordinate graphs were plotted with different potentials (U) for all temperature values.

Figure 1 .
Figure 1.Electrostatic potential diagram of the defect structure at all temperatures (red and blue areas represent the negative and positive charges, respectively, and the increasing density of the color represents the increasing charge density).

Figure 1 .
Figure 1.Electrostatic potential diagram of the defect structure at all temperatures (red and blue areas represent the negative and positive charges, respectively, and the increasing density of the color represents the increasing charge density).

Figure 2 .
Figure 2. Electrostatic potential diagram of the OH + H 2 O intermediate at the temperatures of (a) 298.15K (b) 353.15 K, and (c) 393.15K (red and blue areas represent the negative and positive charges, respectively, and increasing density of color represents increasing charge density).

Figure 2 .Figure 3 .
Figure 2. Electrostatic potential diagram of the OH + H2O intermediate at the temperatu 298.15K (b) 353.15 K, and (c) 393.15K (red and blue areas represent the negative and charges, respectively, and increasing density of color represents increasing charge density

Figure 3 .
Figure 3. Top and side displacements of the *OH + H + *OH intermediate around the 289.32 cm −1 frequency according to vibrational data at 353.15 K (white, grey, blue, orange, and red spheres represent hydrogen, carbon, nitrogen, copper, and oxygen atoms, respectively).The chemical reactivity of the structure increases as the chemical potential decreases.Chattaraj et al. [38] established a direct connection between Mulliken charge negativity and chemical potential.For all temperatures, the chemical potentials of the O 2 and 2H 2 O intermediates, and the Cu 2 -N 8 /Gr structure remain unchanged.The chemical potentials of O 2 and 2H 2 O, which are −3.615 and −3.322 eV, respectively, indicate that O 2 has higher reactivity than H 2 O. Therefore, the optimized 2H 2 O intermediate is more stable, suggesting that it can easily detach from the catalyst structure due to its weak bond strength and long bond length.The chemical potential of the OH + H 2 O intermediate decreases in the following order for different temperatures: 298.15K > 393.15K > 353.15 K.This order aligns with the high reactivity of the OH + H 2 O intermediate observed in the chemical hardness data.The electrophilicity indices of O 2 and 2H 2 O at all temperatures are 8.854 and 8.048 eV, respectively.Since the electrophilicity index represents the ability to attract electrons from the surroundings to form a stable structure, the O 2 intermediate is more reactive than 2H 2 O, as explained in the chemical potential data.For the OH + H 2 O intermediate, the lowest and highest electrophilicity indices are observed at temperatures of 393.15 and 353.15 K, respectively, which is in agreement with the chemical hardness data.The order of reactivity for the OH + H 2 O intermediate is as follows for different temperatures: 353.15K > 298.15K > 393.15 K, which confirms the chemical potential data, indicating its highest reactivity at 353.15 K.The infrared (IR) spectra, illustrated in Figure 4a-d, were examined for the catalyst (*), *O 2 , *OH + H 2 O, and *H 2 O + H 2 O. Notably, the peaks remained unchanged with variations in temperature, except for *OH + H 2 O, as depicted in Figure 4c.This consistent behavior in the IR spectrum substantiates stability across temperature variations and aligns with the results obtained from other calculations, such as HOMO-LUMO energy calculations assessing stability and reactivity.

Figure 4 .
Figure 4. IR spectra of different temperatures for (a) the catalyst (*), (b) *O 2 , (c) OH + H 2 O, and (d) H 2 O + H 2 O.Moreover, Figure 4c highlights that at temperatures of 298.15K and 393.15 K, similar peaks are observed.However, at 353.15 K, shifts or variations in peaks are evident, indicating distinct behaviors of the *OH + H 2 O intermediate at this temperature.Simulation results suggest vibrational activity around the frequency of 289.32 cm −1 for the OH + H 2 O intermediate, while higher peaks signify vibrations of C and H atoms at the catalyst layer's edge.Furthermore, at 353.15 K, certain peaks of the OH + H 2 O intermediate exhibit a shift towards lower wavenumbers within the range of 1200-1600 cm −1 .This observation implies thermal denaturation, as explained by Panick et al. [44], or a phase change associated with temperature, as elucidated by Guan et al. [20].These findings contribute valuable insights into the temperature-dependent dynamics of the studied OH + H 2 O intermediate.
2.144 Å between the Cu and O atoms of the catalyst intermediate confirms the weak bond strength.However, according to the optimized ORR steps, the possibility of a single H 2 O intermediate forming at the end of the reaction is not significant.Due to the presence of the OH + H + OH intermediate, the possibility of the formation of a 2H 2 O intermediate is substantial at the end of ORR.The binding energies of the 2H 2 O intermediates are −2.12,−2.11, and −2.10 eV at temperatures of 298.15, 353.15, and 393.15 K, respectively.When compared to the other ORR steps, the 2H 2 O intermediate possesses the weakest bond strength and the longest bond length at all temperatures, favoring the situation.
In the Cu 2 -N 8 structure, where one Cu atom is attached to four N atoms, this arrangement is responsible for the electron cloud depletion around the Cu atoms.Consequently, this can lead to the electrons around the oxygen of the OOH intermediate being highly attracted to Cu atoms.The high binding strengths of the 2OH intermediate, which are −10.97,−8.51, and −8.50 eV at the temperatures of 298.15, 353.15, and 393.15 K, respectively, support this explanation.When a H atom is attached to the first oxygen atom, it can draw electrons away from the first oxygen atom.This effect weakens the O-O bond of the OOH intermediate and leads to the generation of the 2OH intermediate.
2.144 Å between the Cu and O atoms of the catalyst intermediate confirms the weak bond strength.However, according to the optimized ORR steps, the possibility of a single H2O intermediate forming at the end of the reaction is not significant.Due to the presence of the OH + H + OH intermediate, the possibility of the formation of a 2H2O intermediate is substantial at the end of ORR.The binding energies of the 2H2O intermediates are −2.12,−2.11, and −2.10 eV at temperatures of 298.15, 353.15, and 393.15 K, respectively.When compared to the other ORR steps, the 2H2O intermediate possesses the weakest bond strength and the longest bond length at all temperatures, favoring the situation.
Consequently, this can lead to the electrons around the oxygen of the OOH intermediate being highly attracted to Cu atoms.The high binding strengths of the 2OH intermediate, which are −10.97,−8.51, and −8.50 eV at the temperatures of 298.15, 353.15, and 393.15 K, respectively, support this explanation.When a H atom is attached to the first oxygen atom, it can draw electrons away from the first oxygen atom.This effect weakens the O-O bond of the OOH intermediate and leads to the generation of the 2OH intermediate.

Figure 6 .
Figure 6.Optimization process when the H atom is introduced to the first O atom of the OOH intermediate at a temperature of 298.15 K. (a) Initially, the H atom binds to the first O atom of OOH to form H2O2, (b) in the middle state of optimization, and (c) in the final optimization step where

Figure 6 .
Figure 6.Optimization process when the H atom is introduced to the first O atom of the OOH intermediate at a temperature of 298.15 K. (a) Initially, the H atom binds to the first O atom of OOH to form H 2 O 2 , (b) in the middle state of optimization, and (c) in the final optimization step where dual OH molecules are formed instead of H 2 O 2 (white, grey, blue, orange, and red spheres represent hydrogen, carbon, nitrogen, copper, and oxygen atoms, respectively).
2 O 2 generation by the OOH intermediate.Due to electron depletion around the terminal O of *OOH, the O atom could not bear two H atoms as H 2 O.This could cause the initially bonded H atom to break its bond with the second O atom.

Figure 7 .
Figure 7. Optimization process when a H atom is introduced to the second O atom of the OOH intermediate at the temperature of 298.15 K: (a) initially, the H atom binds with the second O atom of the OOH intermediate to form H2O; (b) in the next step, the H atom binds to the second O atom, causing the previously bonded H atom to break its bond and migrate to the first O atom; and (c) in the final optimization step, a dual OH molecule is formed instead of H2O (white, grey, blue, orange, and red spheres represent hydrogen, carbon, nitrogen, copper, and oxygen atoms, respectively).

Figure 8
Figure 8 illustrates the optimized ORR steps, which differ from previous studies [28].The changes observed in the 2OH and OH + H + HO intermediates, as opposed to O + H2O and OH +H2O in the ORR, may be attributed to the dual metal catalyst structure, as explained before.The following reaction Equations (1)-(5) depict the optimized ORR steps from the current study in an acidic medium.The asterisk (*) indicates the defect's configuration.

Figure 7 .
Figure 7. Optimization process when a H atom is introduced to the second O atom of the OOH intermediate at the temperature of 298.15 K: (a) initially, the H atom binds with the second O atom of the OOH intermediate to form H 2 O; (b) in the next step, the H atom binds to the second O atom, causing the previously bonded H atom to break its bond and migrate to the first O atom; and (c) in the final optimization step, a dual OH molecule is formed instead of H 2 O (white, grey, blue, orange, and red spheres represent hydrogen, carbon, nitrogen, copper, and oxygen atoms, respectively).

Figure 8
Figure 8 illustrates the optimized ORR steps, which differ from previous studies [28].The changes observed in the 2OH and OH + H + HO intermediates, as opposed to O + H 2 O and OH +H 2 O in the ORR, may be attributed to the dual metal catalyst structure, as explained before.The following reaction Equations (1)-(5) depict the optimized ORR steps from the current study in an acidic medium.The asterisk (*) indicates the defect's configuration.

Figure 9 .
Figure 9. Free energy diagrams of Cu2-N8/Gr at different temperatures under various cell potentia in an acidic medium: (a) 298.15 K, (b) 353.15 K, and (c) 393.15K (note: the OH + H2O intermedia is the same as the OH + H +OH intermediate, as explained previously).

Figure 9 .
Figure 9. Free energy diagrams of Cu 2 -N 8 /Gr at different temperatures under various cell potentials in an acidic medium: (a) 298.15 K, (b) 353.15 K, and (c) 393.15K (note: the OH + H 2 O intermediate is the same as the OH + H +OH intermediate, as explained previously).
Salam et al., Amrit et al., and Ferraris et al. have reported that the performance of low-temperature fuel cells changes dramatically at temperatures around 298.15, 353.15, and 393.15K due to changes in the particle kinetics with temperature using experimental methods, as previously stated

Table 1 .
Formation energies of the Cu 2 -N 8 /Gr structure determined at the temperatures (T) of 298.15, 353.15, and 393.15K under various external potentials (U).

Table 2
shows that the OH + H 2 O intermediate's stability increases in the following temperature order: 353.15K < 298.15K < 393.15 K.This indicates that at 393.15 K, the intermediate becomes more stable, while at 353.15 K it becomes highly reactive compared to the other temperatures.The ESP diagram in Figure

Table 3 .
The binding energies (BEs) of the Cu2-N8/Gr structure at different temperatures and the shortest distance (d) between Cu-O and O-O atoms in angstroms (Å).

Table 3 .
The binding energies (BEs) of the Cu 2 -N 8 /Gr structure at different temperatures and the shortest distance (d) between Cu-O and O-O atoms in angstroms (Å).Temperature/K Molecule BE/eV d O-O /Å d Cu-O /Å

Table 4 .
Comparison of binding energy (BE), bond length of molecule to catalyst (d M-O ), and formation energy at zero potential (E F(U = 0) ) of different catalysts in this study with the published studies, at initial temperature.M-O is defined as M = metal atom in catalyst surface and O = oxygen atom of molecule.