Bending the ORR Scaling Relations on Zirconium Oxynitride for Enhanced Oxygen Electrocatalysis

Technologies like polymer electrolyte membrane fuel cells play an important role in environmentally friendly energy conversion. Essential for their commercialization is the development of cheap and efficient electrocatalyst for the oxygen reduction reaction (ORR). Non‐platinum group metal (PGM) based catalysts has exhibited favourable activity in acidic electrolytes. In this study, we computationally explore the catalytic sites on the (111) surface of hydrated zirconium oxynitride using periodic DFT calculations. The thermodynamically limiting step is determined to be the O2 activation step, (O2→OOH*) for the majority of the sites, hence strengthening the OOH* intermediate binding will improve ORR activity. The calculations also reveal that the determined scaling relation for GOH‐GOOH ( ΔGOOH=0.78ΔGOH+ ${{\rm{\Delta }}G_{OOH} = 0.78{\rm{\Delta }}G_{OH} + }$ 3.52 eV) is comparable to the standard ORR scaling relation ( ΔGOOH=ΔGOH+ ${{\rm{\Delta }}G_{OOH} = {\rm{\Delta }}G_{OH} + }$ 3.2 eV) however GO*‐GOH scaling ( ΔGO*=1.04ΔGOH+ ${{\rm{\Delta }}G_{O^{\ast} } = 1.04{\rm{\Delta }}G_{OH} + }$ 2.16 eV) deviates dramatically from the typical ORR scaling relation ( ΔGO*=2ΔGOH ${{\rm{\Delta }}G_{O^{\ast} } = 2{\rm{\Delta }}G_{OH} }$ ).


Introduction
The increasing fossil fuel consumption has led to increased levels of atmospheric carbon dioxide. The combustion of fossil fuels has been linked to anthropogenic global warming, which might result in a 4°C increase in the global mean temperature during the 21st century under certain "business as usual" fossilrich scenarios. [1,2] Apart from carbon dioxide emission, fossil fuel combustion is also responsible for the emission of many other toxic air pollutants. The synergistic effect between climate change and air pollution can magnify health hazards. [3] Thus technologies for more environmentally friendly energy conversion are actively being pursued. One such technology is the polymer electrolyte membrane fuel cell (PEMFC).
In PEMFCs, a fuel such as methanol or hydrogen is oxidized at the anode and oxygen is reduced to water at the cathode. [4] Possible applications of PEMFC ranges from stationary to mobile and portable applications, in e. g. auxiliary power units (APU) and fuel cell electric vehicles (FCEV). [5] However the widespread use of hydrogen fuel cells continue to face several challenges the most significant of which is the high cost of the platinum-based catalyst used in the electrodes, which contributes to almost 55 % of the total cost. [6] Ideally to decrease the overall cost of PEMFCs, catalysts for both the electrodes should be replaced but the slow oxygen reduction reaction (ORR) at the cathode requires much more platinum than the much faster hydrogen oxidation reaction at anode. The sluggish kinetics of the ORR at the PEMFC cathode can lead to nearly 30 % efficiency decrease of the PEMFCs. [7] Thus, the development of non-precious metal catalyst for the oxygen reduction reaction, that would be cheap and effective, has received much attention. [8][9][10][11] Some of the possible alternatives can be catalysts based on transition metal in nitrogen doped graphitic cathodes, [11] oxides and oxynitrides of transition metals e. g. tantalum, [12,13] zirconium oxide based catalysts with multiwalled carbon nanotube (MWCNTs) support [14] and niobium-titanium complex oxides, [15] among others. Exploring the ORR mechanisms on various catalyst surface, and consequently the active sites, using Density Functional Theory (DFT), helps us gain fundamental understanding regarding the same. This knowledge can be used in developing rational design strategies for improved ORR catalysts. [16][17][18][19][20] Ideal electrocatalysts for hydrogen fuel cells should allow reversible reaction. For ORR to be reversible the adsorption Gibbs energy DG ads ð Þ profile plotted along the reduction steps should be flat. [21,22] None of the currently identified catalysts, however, satisfy the ideal ORR condition. Even on a platinum surface, the ORR intermediates HOO* is adsorbed too weakly and HO* too strongly.
Nevertheless, finding a material with larger DG OOH ð Þ and smaller DG OH ð Þ is restricted by the linear scaling existing within the DG ads ð Þ of the reaction intermediates which is observed almost universally for many materials. [22,23] Scaling relations refer to the correlation between binding energies of different reaction intermediates. One can explain universal scaling of ORR as strong adsorption of HOO* on a metal surface being associated with similarly strong adsorption of HO* on the same when both adsorbates have similar bonding configuration with the surface. However, this might not hold true for complex oxides and oxynitrides of transition metals since they allow various bonding configurations. These oxides and oxynitrides can provide different adsorption sites like the lattice oxygen, metal sites, oxygen vacancy or impurity sites. [24] Computational investigations of ORR on defective titanium oxide surfaces indicate possible deviations from the standard scaling relations. [24] Experimental investigations also reveal that several transition metal complex oxides and oxynitrides, like that of zirconium, titanium and tantalum, show promising activity as ORR catalysts in terms of lower overpotential than platinum ORR catalysts. [14,24,25] Hence, we would like to explore the possibility of oxynitride of another transition metal, i. e. zirconium, exhibiting deviations from standard scaling relations while showing lower overpotential than the standard platinum catalysts.

Structure of Zirconium Oxynitride
Zirconium Oxynitride are compounds with structure similar to fluorite type crystal in which part of oxygen of zirconia is replaced by nitrogen and anion vacancies. Zr 2 ON 2 crystal belongs to space group Ia-3 (206) or according to the Pearson notation has the Pearson symbol of cI80. [26,27] The Pauling electronegativity difference between for oxygen (3.50) and nitrogen (3.07) gives the oxynitride covalent characteristics, which different from oxides with ionic characteristics. [28] Bulk optimisation of the representative structure of Zr 2 ON 2 gives us a cubic lattice crystal structure with lattice constant of 1.01 nm which agrees very well with the experimental value of 1.0135 nm. [26] In our study we are exploring ORR on the (111) surface of Zr 2 ON 2 slab. [29][30][31] It has a hexagonal lattice structure with lattice parameter a = 1.42 nm and a total vacuum of 1.6 nm (with 0.8 nm on both side of the slab). Water, being one of the products of the ORR, is expected to be present in the fuel cell, hence we investigate the possibility of water molecule(s) interacting with and dissociating on the catalyst surface. From Pourbaix stability analysis, [32] we conclude that water molecules interacts with the catalyst surface, dissociates over it and gets adsorbed on the surface. Various dissociation and the subsequent adsorption mechanism is studied as illustrated in Figure 1c. We presumed the water molecule could, potentially, split in different ways, either into an O atom and an OH ion or an O atom and two H atoms. These species, so formed, could further be adsorbed on the surface in many different ways. For example, the OH ion could be adsorbed into an intrinsic vacancy site and H on a surface N. Further, we have also tested low coverage (1 site) and high coverage (4 sites) configuration. The entries in the legend of Figure 1c represent the different configurations tested. The subscript V denotes an interstitial vacancy site, and the subscript L denotes a lattice site; for example, N L and O L represents nitrogen and oxygen lattice adsorption site. O V in the label stands for an oxygen atom getting adsorbed in an intrinsic vacancy site. In contrast, an N L and O L associated with H is a hydrogen atom getting adsorbed over an oxygen or nitrogen atom lattice. For example, the (O V + N L H + N L H 2 ) would stand for the configuration where, after the water split, the O atom has adsorbed into an intrinsic vacancy site. In contrast, one H atom has adsorbed onto a surface N while two have adsorbed onto another N. According to the stability plot for water splitting over Zr 2 ON 2 slab, the energetically most favorable configuration is water dissociating into an O atom and two H atoms with the O getting adsorbed in the intrinsic anion vacancy sites and the H on surface N, the configuration of the same is illustrated in Figure 1b. The structural illustration for the rest of the configurations of water splitting has been provided in the Supporting Information. The catalytic activity of the hydrated Zr 2 ON 2 towards ORR is evaluated systematically by studying all the Zr-ontop site for their ORR activity.

Free Energy Diagram
It is assumed the ORR catalysed by Zr 2 ON 2 proceeds through an associative mechanism, which involves hydrogenation of adsorbed molecular oxygen on the catalytic surface [33] which is evidenced from the free energy diagrams (FED) in Figure 2. The possibility of ORR on stoichiometric Zr 2 ON 2 surface taking place through dissociative mechanism was also investigated. Since the dissociative mechanism involves direct dissociation of O 2 without the formation of the OOH* intermediate, the energy barrier of O 2 dissociation is relevant in this scenario. Our test calculations found a high O 2 dissociation barrier in the range of 1.2 to 1.55 eV. A barrier of more than 0.8 eV is generally considered unfavorable for the dissociative mechanism at low temperature; therefore, we concluded that the Zr 2 ON 2 surface favours the associative mechanism for ORR.
We calculated the Gibb's free energy of adsorption for the three ORR intermediates i. e. hydroxyl (OH*), hydroperoxyl (OOH*) and oxide (O*) using the computational hydrogen electrode (CHE) approach and have investigated all the Zrontop sites as potential catalytic sites and the free energy diagrams (FED) show that majority of them provide stable binding sites for these ORR intermediates. We had also sampled various orientations of the adsorbates over the Zr 2 ON 2 slab and the reported adsorption energies correspond to the geometric configurations that give us the lowest adsorption energy. With respect to adsorption of different intermediates on the catalytic sites of the hydrated Zr 2 ON 2 (111) slab, we can make certain observations. During adsorption of OH and OOH intermediates; they adsorb directly on all of the ontop Zr sites. Although for most of the ontop-Zr catalytic sites, direct adsorption of the O* intermediate also takes place, for certain sites they form NHO* or OÀ O* peroxide-like intermediates upon adsorption. Even in the cases where there is no complex formation, we can see some of the O* adsorbs on the Zr-ontop sites while others move towards a bridge site.
The FEDs for each adsorption sites are initially constructed at the equilibrium potential. (All the FEDs are provided in the Supporting Information). The favourable ORR catalytic sites are the ones for which the constructed FEDs have reaction steps have potential energy that are downhill at the highest possible potential. Within the CHE model the highest potential at which all the reaction steps are downhill in free energy is called the thermodynamic limiting potential, U L ð Þ. The theoretical overpotential (η) is defined as the difference between the equilibrium potential and the values of U L . The η as well as the U L gives us the measure of the activity of a catalyst. Higher the U L and correspondingly lower the η implies higher catalytic activity.
In Figure 2, the FEDs have been drawn at an overpotential of 0.4 V. A FED constructed at equilibrium potential is provided in the Supporting Information for reference. In Figure 2  When we observe all the FEDs together, we can observe that there are different values of U L corresponding to every Zrontop sites. The values for the U L and η corresponding to every Zr-ontop sites is given in Table S1 of the Supporting Information. This brings us to the conclusion that different catalytic sites possess different catalytic activity, where the sites with higher U L are more active. We observe the presence of two catalytic sites where the U L is negative; this means those two sites are unfavorable for the ORR in an H 2 fuel cell. These two sites, site A and C, along with another site, G, with very low limiting potential (U L = 0.1), are the ones where formation of a NHO* complex takes place upon O* adsorption. The FED corresponding to site A is illustrated in Figure 2a. Hence this leads us to believe that the formation of this NHO* complex leads to deactivation or very low activity of the catalytic sites. We also notice the formation of a peroxide-like complex upon O* adsorption on certain sites. For sites B, F, and K the value of U L is high (< 0.5 V) while for other site E it is low (> 0.15 V). Hence it can be concluded that the OÀ O peroxide-like complex formation does not seem to have any effect on the catalytic activity of those sites. For most of the sites, however, the adsorption of O* does not give rise to formation of any complex.
We have also, very briefly, explored the possibility of Zr 2 ON 2 surface as a catalyst for oxygen evolution reaction (OER). The reader can refer to the Supporting Information regarding the same.
We had, also, investigated the possibility of the Zr 2 ON 2 surface undergoing a Mars van-Krevelen mechanism [34] involving N or O vacancies and computed the O and N vacancy formation energies for the same. However the O and N vacancy formation energies were quite high, in the range of more than 2 eV for N-vacancy formation and around 3.5 eV for O-vacancy formation. These values lead us to conclude that such a mechanism was unlikely to take place on the Zr 2 ON 2 surface. A detailed table with all the vacancy formation energies is given in the Supporting Information. Additionally, vacancy formation energy values indicate the stability of the Zr 2 ON 2 surface. This conclusion is also corroborated by various experimental studies. [35,36]

Scaling Relations and Volcano Plots
Observing Figure 3a and 3b, we can conclude that scaling relations exist between the ORR intermediates on Zr 2 ON 2 catalyst. The OOH* or O* vs OH* scaling relations imply that there exists a single independent variable that can be used to describe the binding energies of all the ORR intermediates on the catalyst surface. We have used OH adsorption free energy as the descriptor for our scaling relations. (Note: The O* adsorption energies corresponding sites where NHO* complex formation takes place are not included in the G O* -G OH scaling relation since they form outliers.) Since the theoretical overpotential is a function of free energy and correspondingly that of binding energy, scaling relations are also used to determine the limiting potential. From the scaling relation plots in Figure 3, it can be seen that there exists a linear correlation between the the adsorption energies of OOH* and O* vs OH*. The scaling relations that are determined for the Zr 2 ON 2 catalyst from the adsorption energies of the intermediates are as follows: The determined scaling relation for G OH -G OOH is comparable to that of the standard ORR scaling relation: however the one for G O* -G OOH deviates dramatically from the typical ORR scaling relation: The standard scaling relation is typically observed for a wide range of materials like porphyrins, [37] boron nitride, [38] MoS 2 , [38] ultrathin metal oxide sheets, [38] MN 4 motifs in graphene, [38] carbons and HAB-based coordination polymers. [22,39] A previous study [40] reports scaling relations for platinum slab with explicit solvation with no significant deviation from the standard scaling relation. This leads makes us believe hydration is not the reason for the dramatic deviation of the G O* -G OH scaling that we observe here.
As we have discussed earlier, the OH and OOH intermediates adsorb directly on all of the ontop Zr sites without exception. However this is not true for the O* intermediate. On few sites it forms complexes like the NHO* complex or OÀ O* peroxide-like complex. Even in absence of complex formation, adsorption on bridge sites is observed for some sites while for others absorption on ontop sites is observed. We refer to this phenomena as the atypical adsorption of O* intermediate on the Zr 2 ON 2 surface and believe this is the cause for deviation in OOH vs O* scaling relation. This is evident from Figure 3b where we can see the outliers circled in blue corresponding to sites where NHO* formation takes place. Hence we believe the atypical adsorption of O* intermediate also causes the deviation in ΔG O* = 2ΔG OH scaling relation.
While discussing the binding energies of the ORR intermediates to the catalytic sites on the Zr 2 ON 2 surface it is worth mentioning that most sites bind the ORR intermediates too weakly. Hence, contrary to many transition metal ORR catalysts where *OH!H 2 O is the potential limiting step, [32,41] in case of Zr 2 ON 2 catalyst O 2 activation, O 2 !*OOH, is the potential limiting step for majority of the catalytic sites. To maximise the descriptive power of the volcano plot, we therefore use *OOH adsorption free energy as the descriptor instead of using *OH adsorption free energy. The volcano plot with the commonly used *OH descriptor is shown in Supporting Information. The volcano indicates that the minimum theoretical overpotential of the Zr 2 ON 2 catalyst is around 0.45 eV. This value of overpotential is in agreement with the experimentally determined values for ZrO x N y -MWCNT catalyst at a current density of 0.1 A cm À 2 . [14] Similar values have been reported by other experimental studies as well. [42][43][44] The weak binding of ORR intermediates further agrees with the experimentally observed O 2 reduction. [17] The weak binding of reaction intermediates also indicates that several sites are also selective for the 2 electron pathway of ORR. [17] The scope of our study is limited to the thermodynamics of reaction intermediates, which have been found to well describe trends in ORR activity on the terraces of metal surfaces well. Incorporating surface kinetics is expected to provide similar trends as the thermodynamic analysis while suppressing the activity around the apex of the volcano. [16]

Density of States
The total and partial density of states (TDOS and pDOS) are plotted to clarify the bond formation. As SCAN functional is better suited for computing band gaps, we use it for our density of states calculations. [45] In our DOS plot the Fermi level (E f ) is at 0 eV. Figure 4 shows the total and partial density of states of stoichiometric Zr 2 ON 2 slab, hydrated Zr 2 ON 2 slab and a hydrated Zr 2 ON 2 slab with an O* adsorbate. From the Figures 4a  and 4b, we can see that the there is a reduction of band gap upon hydration of the Zr 2 ON 2 slab. The stoichiometric slab has a band gap of 0.8 eV which reduces to 0.4 eV when the slab is hydrated. This implies that the hydration of the slab makes it a better conductor. The fermi level is within the band gap, and the top of the valence band is formed mainly by N 2p states and to a lesser extend O 2p and Zr states. The conduction band has mainly Zr 3d character. As we can see from Figures 4b and  4c, there is a increase in DOS and N pDOS intensity around the E f upon adsorption of O* on the hydrated Zr 2 ON 2 slab. Upon adsorption of O*, the fermi level moves below the top of the valence band indicating that N, O, and Zr contributes with charge to bind the O* adsorbate. The magnitude of DOS at E f can serve as indicator for the ability to form bonds with adsorbed species. [46] Along with the DOS and pDOS plots, we have also studied several electronic properties like band centre, band-width and fractional band filling. [47][48][49] We have also tried to find out if relationships between electronic properties and binding energies exist. However none of the analysis carried out could capture the complexity of the Zr 2 ON 2 catalyst.
In charge density difference maps (present in the Supporting Information) we can see two different kinds of spatial distribution of electronic charge density denoted by different colours. The yellow regions signify increase in electronic charge density and blue region signify decrease in electronic charge density in those regions. This helps us to draw the conclusion that the charge transfer occurs at the adsorption site. Most significant charge transfer is near the surface with the charge being transferred from the zirconium atom to the O* adsorbate. Some charge transfer is also observed for the neighbouring nitrogen and oxygen atoms near the ontop-Zr sites.

Conclusions
Our DFT calculations show that Zr 2 ON 2 is a viable candidate as an ORR catalyst with Zr-ontop sites acting as the active catalytic sites. However, the sites that shows NHO* complex formation during O* adsorption are either inactive or very mildly active sites. Ideal electrocatalysts for ORR must allow it to occur reversibly and complex oxides of transition metals are viable candidates for the same. [50] Our calculations also revealed that the universal scaling line between O* and OH* adsorption free energy can be broken by using a complex oxide surface of a transition metal oxynitride. This does, however, not lead to improved ORR activity. On the basis of scaling relations, we also plot the activity volcano which helps us draw the conclusion that the maximum limiting potential for the Zr 2 ON 2 catalyst is approximately 0.70 V which agrees with the experimentally determined onset potential. [14,42,43] From the ORR kinetics and FEDs of individual catalytic sites we also understand that all sites are not equally active.

Computational Details
All density functional theory (DFT) calculations are done using the Vienna Ab Initio Simulation Package (VASP) and the core electrons are described with the projector augmented wave (PAW) method. [51,52] The Perdew-Burke-Ernzerhof functional (PBE) is used to describe the exchange and correlation energy. [53] The Atomic Simulation Environment (ASE) is used to set up and analyse the structures, which are optimised with a plane wave cutoff of 500 eV and Brillouin-zone integration is performed with a Monkhorst-Pack k-point mesh of 2 × 2 × 1. [54,55] The self-consistent electron density loop is converged to 10 À 6 eV and the structures are relaxed until all forces are below 0.02 eV Å À 1 .
We have studied the (111) surface of zirconium oxynitride (Zr 2 ON 2 ) which is modeled by a periodic non-polar slab. The model includes 3 cationic layers with the bottom layer fixed in its position and the upper two layers allowed to relax. The associative mechanism for ORR, given below, is well established, where the reaction proceeds through four electron-proton transfer steps. [40] * þ O 2 þ H þ þ e À ! * OOH (5) where * denotes the active sites on the catalyst surface.

Adsorption Free Energy
The standard computational hydrogen electrode (CHE) approach is used to calculate the reaction free energies of electrochemical reaction steps. [21] Using the CHE, the reaction free energy, ΔG i , of a reduction step can be written as: where DG 0 i is the reaction free energy at 0 V versus the reversible hydrogen electrode (RHE) and U L is the potential measured against the reversible hydrogen electrode. The DFT reaction energy corresponding to the following reactions are defined here as the adsorption energies. * þ 2H 2 O ! * OOH þ 1:5H 2 ; D ads EðOOHÞ (12) For the adsorbed species, the Gibb's free energies are calculated using the ASE package with the assumption that adsorbate degrees of freedom can be approximated by independent quantum mechanical harmonic oscillators and the vibrational frequencies of only the adsorbates are calculated. The free energy calculations for gas molecules are carried out with the ideal gas approximation [56] and those for adsorbates are carried out with harmonic approximation [57][58][59] where finite differences with 2 displacements of 0.01 Å are used to calculate force constants. The Gibb's free energy includes the zero point energy and entropy. For the gas phase molecules, Gibb's free energies are calculated assuming them to be ideal gas molecules. We use the equilibrium potential derived from DFT calculations in this study which is determined to be 1.15 V, slightly lower than the experimentally determined one of 1.23 V. The experimentally determined equilibrium potential corresponds to the value of reaction free energy of the overall reaction O 2 þ 2H 2 ! 2H 2 O ð Þ that is 4.92 eV. However, the calculated value of equilibrium potential of 1.15 V is determined for our calculations of the H 2 O formation free energy of 2.3 eV with the PBE functional. A simple one-to-one correspondence between equilibrium potential and adsorption free energy helps us determine the theoretical equilibrium potential of 1.15 V. Difference is due to DFT error on the reaction energy.