# Hydroxyl-Induced Modification of Oxygen Activation and Desorption Free Energy on Defective Tetragonal Zirconia Catalysts

Sara Fazeli, Pascal Brault, Amaël Caillard, Eric Millon

► To cite this version:

Sara Fazeli, Pascal Brault, Amaël Caillard, Eric Millon. Hydroxyl-Induced Modification of Oxygen Activation and Desorption Free Energy on Defective Tetragonal Zirconia Catalysts. Journal of Physical Chemistry C, 2024, 128 (47), pp.20066-20080. 10.1021/acs.jpcc.4c04866. hal-04931089

HAL Id: hal-04931089

https://hal.science/hal-04931089v1

Submitted on 5 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Hydroxyl-Induced Modification of Oxygen Activation and Desorption
Free Energy on Defective Tetragonal Zirconia Catalysts

Sara Fazeli¹,²*, Pascal Brault¹,²*, Amaël Caillard¹, Eric Millon¹

¹GREMI, CNRS - Université d'Orléans, 45067 ORLEANS Cedex 2, France

²Simulation moléculaire pour tous - MS4ALL, 1 Avenue du champ de Mars, 45100 Orléans, France

seyedehsara.fazeli@univ-orleans.fr, pascal.brault@univ-orleans.fr

## Abstract

Describing the activation of O₂ on metal surfaces is crucial for understanding fundamental electrochemical processes, such as the oxygen reduction reaction (ORR) in hydrogen fuel cells. This study explores how defects influence O₂ adsorption mechanisms on a zirconia-based cathode. In the first step, we model O₂ adsorption on two defective surfaces: oxygen-deficient t-ZrO₂-x and oxynitride t-ZrO₂-xNₓ, in an aqueous solution. We describe various O₂ adsorption states by analyzing charge transfer and cohesive energy changes in O₂ molecules, Zr active sites, and defects. The results suggest that O₂ adsorption mechanisms on the surfaces of t-ZrO₂-x and t-ZrO₂-xNₓ occur through dissociative and associative pathways, respectively. Additionally, O₂ adsorption on t-ZrO₂-xNₓ leads to the departure of N dopants from the surface, which is unfavorable for catalytic activity. In the second step, we modify the surfaces of t-ZrO₂-x and t-ZrO₂-xNₓ with the hydroxyl (OH) group. Afterward, we simulate the O₂ activation process on these modified surfaces and identify the most probable active sites. Our findings reveal that OH groups stabilize N dopants on hydroxylated t-ZrO₂-xNₓ, preventing their loss. Moreover, OH groups influence the O₂ adsorption mechanism on t-ZrO₂-x, shifting towards associative O-O bond breaking. Conversely, O₂ adsorption on hydroxylated t-ZrO₂-xNₓ remains molecularly associative. Overall, on hydroxylated surfaces, O₂ adsorption involves stronger charge transfer among oxygen, defects, and Zr active sites. In the third step, we explore O₂ desorption trends from these surfaces. This entails analyzing O₂ desorption using steered molecular dynamics (SMD) to generate potential mean force (PMF) profiles and applying Jarzynski's equality to calculate the free energy of desorption. Herein, we find that the free energy of O₂ desorption from hydroxylated surfaces is lower, indicating a more spontaneous process compared to t-ZrO₂-x and t-ZrO₂-xNₓ. Moreover, we discover that oxygen has the highest tendency to desorb from the hydroxylated-ZrO₂-x surface, attributed to the lowest free energy involved in pulling oxygen from the surface, potentially influencing ORR acceleration.
These findings offer valuable guidance for developing efficient non-platinum-based cathode materials, particularly in catalysis applications.

**Keywords:** tetragonal zirconia, defects, catalysts, adsorption, free energy, potential mean force, jarzynski's equality, molecular dynamics simulation

# 1. Introduction

Motivated by the constant reduction of finite fossil resources¹, the discovery of renewable², and affordable energy³ alternatives has become one of the most significant scientific challenges. This has prompted the development of generating energy conversion devices utilizing the proton exchange membrane fuel cells (PEMFCs)⁴,⁵,⁶. The oxygen reduction reaction (ORR)⁷,⁸ is crucial in PEMFC devices, but its slow rate at the cathode poses a key challenge in developing efficient and affordable fuel cells. In reality, the ORR reaction involves multiple components such as a metal surface, solvent water molecules, an oxygen molecule, electrolyte anions or cations, and a hydroxide ion⁹,¹⁰ in alkaline media¹¹ or a proton in acidic media¹²,¹³. Providing insights into ORR mechanisms on various catalysts is crucial for enhancing fuel cell efficiency and driving research to discover advanced catalysts. Platinum-based catalysts¹⁴,¹⁵ exhibit promising ORR activity, but their high cost requires substitution with cost-effective, high-performance¹⁶ alternatives such as transition metal oxide (TMOs)-based catalysts¹⁷,¹⁸,¹⁹. TMOs exhibit enhanced specific activity for ORR due to their unique electronic structures and surface properties. These materials would perform comparable or superior to platinum-based catalysts under particular conditions²⁰.

Among different TMOs, tetragonal zirconia (t-ZrO₂)¹⁷,²¹ is a stable electrocatalyst in advanced electrochemical devices operating in aqueous environments. The catalytic effectiveness of t-ZrO₂ relies on its stability in the tetragonal phase²². This stability is achieved through methods such as doping²³,²⁴ and inducing oxygen vacancy (V_o)²⁵,²⁶ defects. We specifically introduced oxygen vacancies and nitrogen doping into t-ZrO₂ to enhance its catalytic properties for ORR by modifying its electronic structure, reducing the band gap, and increasing active sites for oxygen interaction²². Nitrogen doping was selected based on its known ability to improve catalytic activity in similar systems by altering electronic properties and promoting O₂ adsorption²⁷, ²⁸. Muhammady et al.²⁹ reported the impact of V_o and N-doping on ORR catalytic activities in t-ZrO₂ (101), concluding a minor role of defects in enhancing ORR catalytic
activity on both defective surfaces. Hence, additional research is required to design and optimize targeted t-ZrO₂ catalysts for improved ORR acceleration. For instance, introducing a hydroxyl (-OH) group on a defective t-ZrO₂ surface alters its composition, charge transfer, reactivity, and generates new catalytic sites for ORR³⁰,³¹. These active sites affect vital ORR steps, including O₂ adsorption mechanism, reduction, and ultimately altering catalytic activities. We previously³² showed that hydroxylation of oxygen-deficient and nitrogen-doped zirconia improves their ORR performance. However, the specific roles of these surface modifications differ: nitrogen doping promotes the formation of key intermediates in the ORR, while oxygen vacancies help stabilize these intermediates during the reaction. Therefore, more investigation is needed to determine which modifications-oxygen vacancies, nitrogen doping, or hydroxylation- are most effective for accelerating ORR.

A major motivation for the investigation of O₂ activation on transition metal surfaces is crucial for predicting oxidation catalysis, pivotal in the chemical industry and especially impacting ORR kinetics³³. The O₂ adsorption mechanisms take place through associative or dissociative pathways³⁴. Indeed, the dissociation of O₂ into oxygen is anticipated to be the crucial, rate-limiting step in the ORR³⁵. The strength of O₂ binding to the surface relies on the types of defects and functional groups. This bond strength is associated with the degree of charge transfer to O₂ across the surface, influencing the kinetics of O₂ dissociation and, consequently, the efficiency of ORR³⁶,³⁷. Li et al.¹¹ investigated the impact of OH⁻ ions on O₂ adsorption on Pt surfaces in alkaline solutions. Their findings suggest that the co-adsorption of OH⁻ and H₂O directly affects O₂ adsorption. However, the roles of OH groups, as well as V>O and N-doping defects on t-ZrO₂ surfaces in O₂ adsorption mechanisms, are not fully understood.

Computer simulations complement state-of-the-art experimental techniques, providing an alternative method to explore the precise O₂ adsorption and ORR mechanisms³⁸,³⁹. Numerous theoretical studies have focused on O₂ adsorption mechanisms on metal surfaces and the performance of ORR³⁹,⁴⁰,⁴¹,⁴². Sokolovic et al.⁴³ used DFT to identify O₂ adsorption mechanisms on rutile TiO₂ (110), revealing a preference for dissociative adsorption on surfaces with (V>O) at room temperature. Xia et al.²⁴ investigated the adsorption and dissociation of O₂ on yttria-stabilized zirconia (YSZ) surfaces utilizing DFT calculations. Their findings revealed different mechanisms for O₂ adsorption on pure ZrO₂ and YSZ surfaces. They emphasized the dissociative activation of O₂ at oxygen vacancy sites induced by yttrium.

While the adsorption of O₂ on catalyst surfaces has received significant attention, numerous fundamental questions remain regarding the O₂ desorption mechanism from the surface to interact with other species involved in ORR. Understanding how adsorbed O₂ on catalyst
surfaces, whether associative or dissociative, contributes to ORR is intricately linked to the
crucial parameter of the required free energy for O₂ desorption from the surface⁴⁴,⁴⁵. Certainly,
a lower free energy value directly relates to easier O₂ separation from the surface, thereby
enabling control over the acceleration or deceleration of the ORR process.

Molecular dynamics (MD) simulation<sup>46,47</sup> is a powerful tool to uncover the role of solvent and free energies of O<sub>2</sub> adsorption and desorption on metal surfaces. This is particularly evident through the potential of mean force (PMF) calculated by steered molecular dynamics (SMD) simulations<sup>48</sup>. SMD based on Jarzynski's equality<sup>49,50</sup> has become popular for its simplicity and effectiveness in calculating free energy. Wang et al.<sup>7</sup> used SMD simulations to determine the free energy profiles of oxygen adsorption on Pt surfaces. They found that O<sub>2</sub> adsorption on Pt (110) and Pt (111) is thermodynamically favorable due to negative adsorption free energies. It is crucial to highlight that the mechanisms of O<sub>2</sub> desorption from the surface of t-ZrO<sub>2</sub>, with coexisting defects and hydroxyl groups, are still unclear. These mechanisms govern the required free energy for O<sub>2</sub> desorption and significantly influence ORR performance, remaining a topic of controversy. In this work, we initially used MD simulations to investigate the influence of oxygen vacancy and nitrogen doping on O<sub>2</sub> adsorption mechanisms at the t-ZrO<sub>2</sub>/water interface. Following that, we explored how an OH group affect the charge transfer between O<sub>2</sub> and the surfaces of two defective catalysts, t-ZrO<sub>2-x</sub> and t-ZrO<sub>2-x</sub>N<sub>x</sub>, during associative or dissociative adsorption. Furthermore, O<sub>2</sub> desorption process from defective t-ZrO<sub>2-x</sub> and t-ZrO<sub>2-x</sub>N<sub>x</sub> surfaces, with and without OH, was analyzed via free energy calculation using PMF based on Jarzynski's equality. Since, free energy values are crucial for predicting the ease of oxygen desorption from surfaces, directly affecting the difficulty of the ORR, this study may offer valuable insights for researchers developing non-platinum-based cathodes to accelerate ORR.

## 2. Research Methodology

Molecular dynamics (MD) simulations were performed for the t-ZrO₂/water system using COMB3⁵¹ potential in the LAMMPS Package⁵². The details of COMB3 potential are available in Supporting Information. The simulation cell consists of 96 formula units of tetragonal phase ZrO₂. The structural detail of the unit cell of t-ZrO₂ (001) is illustrated in Figure S1. The pristine t-ZrO₂ was initially relaxed by conducting a canonical NVT⁵³ ensemble simulation for 300 ps. A timestep of 0.1 fs is used. The equation of motion is advanced in time by the Verlet algorithm. The canonical NVT ensemble, wherein the number (N) of atoms, volume (V), and temperature
(T) remain constant, was selected. To maintain a constant temperature at 300K, the Nosé-Hoover⁵⁴ thermostat is employed on the translational degrees of freedom of atoms, utilizing a temperature damping parameter set at 100 fs. Then, we obtained under-stoichiometric oxide t-ZrO₂-x by removing one oxygen atom from the optimized pristine t-ZrO₂. The defective surface of oxynitride t-ZrO₂-xNₓ was achieved with a 2N₀ + V₀ configuration. Details of the process used to identify the most stable structure of t-ZrO₂-xNₓ are provided in the supporting information. Both defective structures underwent relaxation under the same conditions as the pristine t-ZrO₂. Subsequently, to attain hydroxylation of the surfaces, we placed one OH group on both the under-stoichiometric oxide t-ZrO₂-x and oxynitride t-ZrO₂-xNₓ. This allowed for the adsorption of OH on the surfaces under identical simulation conditions. To investigate the O₂ adsorption process on two defective t-ZrO₂ samples, a solution containing 82 water molecules, one O₂ molecule was prepared and placed above the catalysts. Water molecules were randomly inserted to fill the space within the box, resulting in a water slab with an approximate density of 1 g/cm³. The MD simulations were continued until the O₂ adsorption process on each catalyst was completed, ensuring the system reached equilibrium in terms of temperature and pressure. Additionally, there were no changes in the position or bond length of the adsorbed O₂. Free boundary condition (FBC) is applied in the Z direction, and periodic boundary conditions (PBC) are imposed in the X and Y directions.

After completing the O₂ adsorption process, it is essential to investigate the O₂ desorption process from different surfaces. The strength of O₂ binding to the surface directly affects the ease of O₂ desorption, subsequently influencing the next step of the ORR involving interactions with other species. The energy required for O₂ desorption can be linked to the free energy change associated with separating O₂ from the surface. The potential of mean force (PMF) can represent the free energy profiles of the O₂ desorption process. To determine the PMF of O₂ desorption from various surfaces, steered molecular dynamics (SMD) simulations were applied. The SMD simulations were conducted utilizing the “fix SMD” feature within the LAMMPS software. In SMD simulations, time-dependent external detaching forces were applied to an oxygen molecule adsorbed on a catalyst surface to facilitate its separation during the desorption process. In addition, in the SMD simulation, the final configuration of O₂ adsorbed on all surfaces was used as the initial position for the constant velocity pulling to obtain the PMF profile. The distance in the z-direction between the center of mass of the two oxygen atoms in O₂ and a point within the solution was chosen as the reaction coordinate. An external force was applied to the oxygen atoms (in cases of dissociative adsorption) or oxygen molecules (in cases of associative adsorption) to pull them out of the defective and hydroxylated t-ZrO₂ surfaces.
This force was applied in the opposite direction to their initial adsorption process. The pulling velocity in SMD simulations balance two key factors: minimizing nonequilibrium effects and maintaining computational feasibility. A very low pulling velocity would improve the accuracy by reducing irreversible work contributions, but it would also significantly increase simulation time. To investigate the sensitivity of the PMF to pulling velocity, a range of constant velocities from 0.001 to 0.005 Å/ps was tested. This range, as demonstrated in the literature<sup>7</sup>, provides an effective balance between accuracy and computational efficiency. The resulting PMF profiles showed minor shifts, but the overall trends and conclusions remained consistent. Importantly, the free energy difference for O<sub>2</sub> desorption across the surfaces varied by less than 5%, well within the acceptable range for nonequilibrium simulations. Based on this, a pulling velocity of 0.005 Å/ps was selected for all subsequent SMD simulations. The pulling distance was set to 4Å above the surfaces and within the solution. The spring constant, k, was set to 20 kcal/mol·Å<sup>2</sup>, which is consistent with the value reported in the literature<sup>7,11</sup>, providing sufficient stiffness to maintain the reaction pathway without over-constraining the system. This choice ensures a smooth detachment of the O<sub>2</sub> molecule while minimizing oscillations and deviations during the pulling process. While the chosen setup balances accuracy and efficiency, we acknowledge potential limitations, particularly related to nonequilibrium effects at higher pulling velocities. Although our relatively slow pulling velocity mitigates this, future studies could explore slower velocities or alternative methods like umbrella sampling to better capture equilibrium conditions. Additionally, the constant spring force used was effective, but future investigations could explore adaptive spring constants that adjust dynamically based on system response for further accuracy. To acquire the PMF profile for O<sub>2</sub> desorption from defective t-ZrO<sub>2</sub>, SMD simulations were conducted ten times, and the free energy was derived from the average work across these ten trajectories using Jarzynski's equality. Jarzynski's equality defines the connection between equilibrium-free energy differences and the work done during nonequilibrium processes. Consider a procedure where a parameter of the system, γ, shifts from its initial value γ<sub>0</sub> at time zero to its value γ<sub>t</sub> at time t. The second principle of thermodynamics states that the average work performed on the system cannot be less than the disparity between the free energies associated with the initial and final values of the parameter, denoted by γ<sup>50</sup>.

$$ \Delta F = F(\gamma_t) - F(\gamma_0) \le \langle W \rangle \quad (1) $$

Where the equality holds exclusively under conditions of quasi-static processes. In accordance with this inequality, a nonequilibrium process sets only an upper boundary for the difference in
free energy. However, Jarzynski's equality (as documented in the literature) is a relationship that remains applicable regardless of the speed of the process⁵⁰:

$$e^{-\beta\Delta F} = \langle e^{-\beta W} \rangle \qquad (2)$$

where $\beta = 1/k_B T$ is the inverse temperature and $k_B$ is the Boltzmann constant. $\Delta F$ indicates the free energy difference, illustrated by PMF profiles in this study, while $w$ signifies the nonequilibrium work performed in SMD simulations. It can be computed by integrating the force exerted on the target molecule's center of mass over the pulling distance for each pulling process:

$$W(t) = f(t) * \vartheta \int_{0}^{t} dt \qquad (3)$$

### 3. Results and Discussions

#### 3.1. $O_2$ adsorption on t-ZrO$_{2-x}$ and the hydroxylated t-ZrO$_{2-x}$

Understanding $O_2$ adsorption mechanisms on various surfaces is crucial for evaluating the effectiveness of materials in ORR. Our focus was on the charge transfer between $O_2$ molecule and zirconia surfaces in a solution. A simulation was conducted to model the catalyst/water interface containing of a t-ZrO$_{2-x}$ surface with an $O_2$ molecule, and water layers taking into account solvent effects at 300 K over 0.5ns. Figure 1a illustrates the change in charge over time for each oxygen atom ($O_A$ and $O_B$, as shown in the Figure 2) within the $O_2$ molecule on the t-ZrO$_{2-x}$ surface. Additionally, Figure 1b illustrates the computed cohesive energy of each oxygen atom in $O_2$ during its adsorption onto the t-ZrO$_{2-x}$ surface, providing insights into their respective stabilities. Moreover, Figure 2 depicts the interaction between $O_A$ and $O_B$ at the t-ZrO$_{2-x}$ surface-solution interface, enhancing our understanding of $O_2$ activation mechanisms. The charges on $O_A$ and $O_B$, as well as the bond length between $O_A$ and $O_B$, were calculated to be in the ranges of 0.07 to 0.01, -0.01 to -0.02, and 1.1Å to 1.2Å, respectively, over a duration of 50 ps. After 100 ps, these values were reported to be 0.004, 0.04, and 1.2Å, respectively. Analysis of Figures 1a-b reveal fluctuations in the charge and cohesive energy of $O_A$ and $O_B$ before 500 ps, with the charge increasing and the cohesive energy decreasing. The bond length remains approximately constant. The charges on the Zr$_1$ and Zr$_2$ sites interacting with the $O_2$
molecule change from initial values of 1.72 and 1.80 (Figure 2-I) to 1.94 and 1.88 at 500 ps (Figure 2-III), respectively. This indicates that during O₂ adsorption, charge transfer occurs from the Zr atoms to form bonds with the O₂ molecules. At 500ps, there is a great gap between the charges and cohesive energies of O_A and O_B is observed (Figures 1a-b).

By 500 ps, the bond length between O_A and O_B increases, and the atoms become adsorbed onto two different Zr sites. The O_A-O_B bond length at this step is 1.8Å (It is worth to note that the bond length of O_A-O_B over time is detailed in Table S1). In addition, the charges on O_A and O_B are reported to be 0.12 and -0.17, respectively. The cohesive energies are -1.15 eV for O_A and -5.18 eV for O_B. In addition, the charges on Zr₁ at the site of O_A adsorption and on Zr₂ at the site of O_B adsorption are 1.91 and 1.88, respectively. This discrepancy suggests a dissociative pathway for O₂ adsorption on t-ZrO₂-x, where the bond between O_A and O_B breaks, and each atom adsorbs onto different Zr atoms of t-ZrO₂-x (as shown in Figure 2-III). After 500 ps, until the end of the simulation, the charge and cohesive energy values of both O_A and O_B exhibit minor fluctuations. However, while the charge and cohesive energy of O_A remain relatively constant, those of O_B experience a slight decrease, as depicted in Figure 1a-b. Essentially, the adsorption of O_B occurs at a deeper level within the t-ZrO₂-x structure, rendering it more susceptible to the chemical environment, whereas the adsorption of O_A is closer to the surface. Consequently, O_B exhibits greater stability at the adsorption site compared to O_A (Figure 2-IV). To confirm that O₂ adsorption on the surface of t-ZrO₂-x is dissociative, we calculated the spin densities of oxygen atoms O_A and O_B at two key stages: the initial state (t = 0 ps) and state (t = 3000 ps), where the dissociation and full adsorption of the O atoms on the surface have occurred. In the initial state, the spin densities of O_A and O_B were calculated as 3.29 electron/Å³ and 3.31 electron/Å³, respectively. After 3000 ps, following the completion of the dissociation process, the spin densities of O_A and O_B increased to 4.64 electron/Å³ and 4.39 electron/Å³, respectively. This substantial increase in spin density (from 3.29–3.31 to 4.64–4.39) indicates a stronger interaction between the oxygen atoms and the oxygen-deficient surface, which is consistent with the dissociation of the O₂ molecule. The small spin density difference in the initial state (Δ spin density = 0.02) aligns with the expected values for an isolated O₂ molecule, where the two oxygen atoms have unpaired, spin-polarized electrons. Upon dissociation, the interaction between the individual oxygen atoms and the surface becomes more pronounced, as evidenced by the larger spin density difference in the final state (Δ spin density = 0.25). This increase in spin density suggests electron transfer between the surface and the O₂ molecule, facilitating the bond-breaking process that leads to dissociative adsorption.
Analysis of O₂ interaction with hydroxylated t-ZrO₂-x under identical simulation conditions as for t-ZrO₂-x indicates O₂ movement within water, atop OH groups, with its axis slanting to the surface, indicating strong interaction. It is worth noting that O₂ remains on the surface of hydroxylated t-ZrO₂-x for a longer time compared to its adsorption on the t-ZrO₂-x surface. This prolonged time is attributed to the increased interaction with the OH group before the molecule finally adsorbs onto the surface. As illustrated in Figure 1c and Figure 1d, fluctuations in both the charge and energy levels of atoms Oₐ and Oₓ were observed within the timeframe of zero to 5000ps. After 100 ps, the O₂ molecule is positioned above the OH group, with charges of 0.05 for Oₐ and 0.16 for Oₓ, and cohesive energies of -2.14 eV for Oₐ and -1.71 eV for Oₓ (Figure 3-II). The charges on H and O of the OH group are 0.05 and -0.59, respectively. At 706 ps, O₂ undergoes a configuration change relative to OH, leading to a significant decrease in charges (qOₐ=0.02, qOB=0.06) and cohesive energies (-2.90 eV for Oₐ and -2.20 eV for Oₓ), as shown in the Figure 1c-1d. Notably, the charges on H and O of OH change to 0.14 and -0.58, respectively. By 1407 ps, O₂ and OH experience a configuration change, with O₂ aligning atop the OH group. The charges are 0.09 for Oₐ and -0.02 for Oₓ, with cohesive energies of -2.50 eV for Oₐ and -2.92 eV for Oₓ. At 3000 ps and 3500 ps, there is a significant decrease in the charges and cohesive energies for Oₐ and Oₓ, respectively (Figure 3-III). This indicates increased interaction with the OH group, leading to changes in the charge of H to 0.03 and the charge of O in OH to -0.87 at 3500 ps. It is worth noting that through the simulation up to 4000 ps, the bond lengths of Oₐ and Oₓ remain almost constant, fluctuating within the range of 1.1-1.2 Å. The changes in the Oₐ-Oₓ bond length during O₂ adsorption on hydroxylated t-ZrO₂-x are listed in Table S2.

Notably, at 4000 ps, the O₂ molecule moves towards the nearby Zr atom, adjacent to the OH group. At this point, a sudden drop in charges (qOₐ = -0.08, qOB = -0.47) is observed (Figure 3-IV). As well as, following the adsorption on the Zr site, there was a sudden decrease in the cohesive energy values of Oₐ from -3.40 to -4.29 eV and Oₓ from -2.92 eV to -5.30 eV, respectively. The changes in charge and cohesive energy at the Zr site have been observed in Figure 3. The charge increases from 1.87 before 4000 ps to 1.97 at 4000 ps. Additionally, the cohesive energy shows a sudden decrease from 1.70 eV at 3500 ps to 0.55 eV at 4000 ps. This drop is attributed to the stabilization of both Oₐ and Oₓ sites and the Zr site due to their connection with the Zr atom, typically accompanied by substantial charge transfer from the metal surface to O₂. Interestingly, despite both Oₐ and Oₓ atoms being adsorbed on Zrᵢ (as depicted in the Figure 3-IV, and 3-V), the bond length increases from 1.2 Å to 1.8 Å, leading to the breakage of the Oₐ-Oₓ bond. This kind of associative adsorption process leads to a final
$O_2$ adsorbed structure with higher stability, demonstrated by lower cohesive energy values for $O_A$ and $O_B$, and a more pronounced charge transfer, resulting in more negative charges on $O_A$ and $O_B$ compared to t-ZrO$_{2-x}$. However, in the case of hydroxylated t-ZrO$_{2-x}$, both $O_A$ and $O_B$ are adsorbed on the same Zr site with similar stability. These mechanisms can directly influence the contribution of each oxygen atom in the ORR. These findings suggest that hydroxylation of the t-ZrO$_{2-x}$ surface can significantly improve the efficiency and effectiveness of ORR, making them more promising candidates for catalytic applications in ORR processes compared to the t-ZrO$_{2-x}$ surface. The presence of OH groups facilitates a more favorable $O_2$ adsorption pathway, leading to enhanced stability, prolonged adsorption duration, and increased charge transfer, all of which are desirable traits for ORR catalysts.

<figure><img src="image_2.png" /><figcaption>Figure 1. (a) The variation in charge on each oxygen atom ($O_A$ and $O_B$) within the $O_2$ molecule on the t-ZrO$_{2-x}$ surfaces. (b) The cohesive energy changes of each oxygen atom in $O_2$ during adsorption onto the surface of the t-ZrO$_{2-x}$ over time (c) the variation in charge on each oxygen atom ($O_A$ and $O_B$) within the $O_2$ molecule on the hydroxylated t-ZrO$_{2-x}$ surfaces (d) the cohesive</figcaption></figure>
energy changes of each oxygen atom in O₂ during adsorption onto the surface of the
hydroxylated t-ZrO₂-x

<figure><img src="image_3.png" /><figcaption>Figure 2. Evolution of the O₂ adsorption process on the surface of t-ZrO₂-x in aqueous solution.
The labeled snapshots correspond to specific points in Figure 1a.</figcaption></figure>
<figure><img src="image_4.png" /></figure>
Figure 3. Evolution of the O₂ adsorption process on the surface of the hydroxylated t-ZrO₂-x in aqueous solution. The labeled snapshots correspond to specific points in Figure 1c.

## 3.2. O₂ adsorption on t-ZrO₂-xNₓ and the hydroxylated t-ZrO₂-xNₓ

A crucial aspect to consider in studying N-doped ZrO₂ structures is ensuring that the chosen structures are representative of the most stable N-doped configurations. Prior to investigating the O₂ adsorption mechanism on the t-ZrO₂-xNₓ surface, we initially examined the most stable structures of adsorption of OH on various sites of t-ZrO₂-xNₓ. The hydroxylated t-ZrO₂-xNₓ structures are shown in Figure S2. Subsequently, the most stable structure of both t-ZrO₂-xNₓ and hydroxylated t-ZrO₂-xNₓ is chosen, and a model simulating the catalyst/water interface enclosing an O₂ molecule in a solution, under identical conditions to those of the t-ZrO₂-x catalysts, undergoes annealing at 300K for 0.5ns. Initially, we investigated the O₂ adsorption mechanism on the surface of t-ZrO₂-xNₓ. The charge and cohesive energy of individual oxygen atoms in O₂ during their adsorption onto t-ZrO₂-xNₓ surface over time are illustrated in Figures 4a and Figure 4b, respectively. Moreover, Figure 5 illustrates the schematic depiction of the movement and subsequent adsorption of the O₂ molecule on the surface of t-ZrO₂-xNₓ at different steps throughout the simulation.

As depicted in Figure 4a, the initial charges on O<sub>A</sub> and O<sub>B</sub> are 0.01 and 0.03, respectively (Figure 5-I). Additionally, the initial charge on the Zr atom adjacent to the oxygen vacancy (Zr<sub>OV</sub>) is 1.08. After 10 ps, the oxygen molecule forms a bond with Zr<sub>OV</sub>, resulting in electron transfer primarily from the transition metal atom to the oxygen molecule. Consequently, the charge on O<sub>A</sub>, O<sub>B</sub>, and Zr<sub>OV</sub> shifts to -0.21, -0.14, and 1.34, respectively. The bond length between O<sub>A</sub> and O<sub>B</sub> is reported as 1.2Å. Additionally, the cohesive energy of O<sub>A</sub>, O<sub>B</sub>, and Zr<sub>OV</sub> is calculated as -4.74 eV, -3.80 eV, and -2.44 eV, respectively. In this step, four water molecules are simultaneously adsorbed on different Zr sites adjacent to both the N-doped and oxygen vacancy (V<sub>O</sub>). Up to 100 ps, the charges on both O<sub>A</sub> and O<sub>B</sub> fluctuated, eventually changing to -0.28 and -0.12, respectively, by 100 ps, while the charge on Zr<sub>VO</sub> remained relatively stable at 1.35. Moreover, a similar trend is observed in the alteration of cohesive energy for the mentioned atoms, with cohesive energy values recorded as -6.18 eV, -4.83 eV, and -3.20 eV for O<sub>A</sub>, O<sub>B</sub>, and Zr<sub>OV</sub>, respectively (Figure 5-II). In this step, more water molecules are adsorbed on Zr sites, which are opposite site of V<sub>O</sub> and N-doped. As shown in Figures 2a-b, After 100 ps, the charges and cohesive energies of O<sub>A</sub> and O<sub>B</sub> exhibit fluctuations around a constant value. Furthermore, in subsequent simulation steps (200ps), it is shown that more water molecules
tend to adsorb on Zr sites distant from both defects (Figure 5-III). This suggests that the initial adsorption of O₂ molecule on ZrOv site and nearby nitrogen doped may affect surface chemistry and energetics, thereby influencing water molecule adsorption. This suggests that the initial adsorption of an O₂ molecule on the ZrOv site and the nearby nitrogen doping may affect the surface chemistry and energetics, thereby influencing the adsorption of water molecules. The presence of adsorbed O₂ molecule could obstruct water molecule adsorption on nearby Zr site by altering surface coverage and competing for available adsorption sites. Interestingly, at 220 ps, it is observed that due to the absence of water blockage on the side of defects, one of the nitrogen atoms migrates to the surface (Figure 5-IV). After 250ps, it completely exits the surface and begins interacting with the solvent phase. This process is not favorable for catalytic behavior. Hence, it is vital to improve the surface of t-ZrO₂-xNx to prevent the instability of nitrogen atoms. The investigation of O₂ interaction with hydroxylated t-ZrO₂-xNx surface reveals a significantly different behavior compared to t-ZrO₂-xNx one.

Figure 4c and Figure 4d show the variation in charge and cohesive energy of each oxygen atom in the O₂ molecule during its interaction and adsorption on hydroxylated t-ZrO₂-xNx, respectively. As depicted in the Figure 6, the O₂ molecule interacts with the hydroxylated t-ZrO₂-xNx surface, mainly engaging with sites located on the opposite side of the N-doped and VO regions. It is interesting to note that the O₂ molecule is positioned near the OH group, while simultaneously, the configuration of OH shifts horizontally towards the vicinity of the N-doped sites at t = 100 ps (Figure 6-II). The charges on Oₐ and O_b are -0.02 and -0.04, and the corresponding cohesive energies are -2.82 eV and -4.23 eV, respectively. At 403 ps, a sudden rise in the charge of Oₐ and a drop in the charge of O_b are observed (Figure 4c). In this step, the configuration of the O₂ molecule is parallel to the surface, stabilizing on Zr site (as shown in the Figure 6-III). The cohesive energies of the stabilized Oₐ and O_b atoms reach -4.01 eV and -5.35 eV, respectively (Figure 4d). From 403 ps until the end of the simulation, the fluctuations in the charge and cohesive energies of Oₐ and O_b stabilize around a constant value. In addition, the phenomenon of associative adsorption of O₂ is identified at the end of the simulation, characterized by Oₐ-O_b bond length 1.22Å.

It is important to note that during the whole of simulation, both N-doped positions shift closer to the surface, interacting with adsorbed OH group. This leads to a considerable change in the charges on both N₁ and N₂. As shown in Figure 6-III, the charges on N₁ and N₂ changed from -0.24 and -0.17 in the initial state to -0.24 and 0.16 at 403 ps, corresponding to the shift of the N atoms toward the surface. Additionally, the cohesive energies of N₁ and N₂ changed from -1.48 eV and -1.51 eV in the initial state to -2.24 eV and -1.57 eV at 403 ps, respectively. Despite
these movements, they remain stable within the structure until the end of the simulation of t= 3800ps (Figure 6-IV, and 6-V). Furthermore, prior to the transition in N configurations, the charges on the H and O atoms of OH fluctuated between 0.13 to 0.1 and -0.72 to -0.46, respectively. As the N-doped site adjusts its configuration to achieve a more stabilized position through interaction with the OH group, the charges and cohesive energies of the O and H atoms of OH change until 403 ps. It is important to note that before 403 ps, alongside O₂ adsorption, the majority of water molecules are also adsorbed on the surface opposite the defect. Due to the presence of OH group near the defects, N atoms interact with this OH group for an extended period, preventing the site from being blocked by water molecules. This interaction favors the stability of nitrogen doped on the surface. These results represent a significant achievement in surface modification. Hydroxylation of the t-ZrO₂-xNx surface stabilizes N atoms and prevents their expulsion from the surface.

<figure><img src="image_5.png" /><figcaption>Figure 4. (a) The variation in charge on each oxygen atom (O<sub>A</sub> and O<sub>B</sub>) within the O<sub>2</sub> molecule on the t-ZrO<sub>2-x</sub>N<sub>x</sub> surfaces. (b) The cohesive energy changes of each oxygen atom in O<sub>2</sub> during adsorption onto the surface of the t-ZrO<sub>2-x</sub>N<sub>x</sub> over time (c) the variation in charge on each oxygen atom (O<sub>A</sub> and O<sub>B</sub>) within the O<sub>2</sub> molecule on the hydroxylated t-ZrO<sub>2-x</sub>N<sub>x</sub> surfaces (d)</figcaption></figure>
the cohesive energy changes of each oxygen atom in O₂ during adsorption onto the surface of
the hydroxylated t-ZrO₂-xNx.

<figure><img src="image_6.png" /><figcaption>Figure 5. Evolution of the O₂ adsorption process on the surface of the t-ZrO₂-xNx in aqueous solution. The labeled snapshots correspond to specific points in Figure 4a.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6. Evolution of the O₂ adsorption process on the surface of the hydroxylated t-ZrO₂-xN_x in aqueous solution. The labeled snapshots correspond to specific points in Figure 4c</figcaption></figure>
### 3.3 Free energy calculation from SMD simulation

In the previous sections, we illustrated how altering defect types and surface modification via hydroxylation can change the $O_2$ adsorption behavior on defective t-ZrO$_2$, a crucial initial stage in ORR. Indeed, the $O_2$ adsorption mechanism can significantly influence subsequent ORR steps. Additionally, the ease of $O_2$ desorption from the surface, enabling its interaction with other species involved in ORR, influences the catalyst's reactivity. In this section, the $O_2$ desorption processes were studied using the PMF method. The PMF profile represents the free energy landscape of a system along pulling distance. After completing the $O_2$ adsorption on various defective structures of t-ZrO$_2$, the O atoms were pulled out through SMD simulations to obtain the PMF.

PMF profile enable a detailed understanding of underlying atomic-level processes and allow the exploration of system states that are inaccessible through experimental methods. In most cases, free energy calculations focus on relative free energies, with the primary interest being how the free energy changes as a function of an external parameter or an internal coordinate. In this section, we calculate the free energy relative to an internal coordinate using Jarzynski's equality method.

#### 3.3.1 PMF profile of $O_2$ desorption from the surface of t-ZrO$_{2-x}$

As discussed in Section 3.1, during the adsorption of an $O_2$ molecule on t-ZrO$_{2-x}$, the two oxygen atoms dissociate and adsorb onto separate zirconium sites. Subsequently, the desorption of one oxygen atom ($O_A$) could indeed influence the charge state of the remaining atom ($O_B$) due to charge redistribution, and this could potentially affect the adsorption energy of the second oxygen atom. In this case, two individual springs for pulling the oxygen atoms ($O_A$ and $O_B$) from the surface are employed, with the aim of providing a clearer picture of the free energy landscape for each oxygen atom's desorption. This approach allowed us to compute the free energy barrier for each atom separately, which can be useful in understanding the strength of binding for both $O_A$ and $O_B$. Consequently, we applied an external force to each $O_A$ and $O_B$ to detach them from the surface of t-ZrO$_{2-x}$, in the direction opposite to the adsorption process. Figure 7 shows the free energy of $O_A$ and $O_B$ desorption from t-ZrO$_{2-x}$, calculated using Jarzynski's equality, versus pulling distance. The schematic depicting the corresponding $O_A$ and $O_B$ desorption, is illustrated in Figure 8. Additionally, the resulting PMF profile for the $O_A$ and $O_B$ moving out of the mentioned surface is shown in Figure S3.
As shown in Figure 7, the ΔF values for O<sub>A</sub> and O<sub>B</sub> desorption from the surface of t-ZrO<sub>2-x</sub> are negative at pulling distances up to 0.17Å and 0.23Å, respectively. Beyond these distances, the free energy values for both O<sub>A</sub> and O<sub>B</sub> desorption become positive. From a thermodynamic perspective, if the value of the free energy change (ΔF) is negative, the process is spontaneous. A more negative free energy change indicates a more favorable spontaneous process. Conversely, if ΔF is positive, the process is non-spontaneous.

<figure><img src="image_8.png" /><figcaption>Figure 7. Free energy profile based on Jarzynski's equality of O<sub>2</sub> desorption from the surfaces of the t-ZrO<sub>2-x</sub>.</figcaption></figure>

The free energy (and PMF curve) for O<sub>A</sub> rises until 0.68 Å, reaching a value of 3.1 kcal/mol (Figure 7 and Figure 8-II). The ΔF value then continues to increase gradually, with a gentle slope, reaching 3.8 kcal/mol at a pulling distance of 0.95 Å (Figure 8-III). Beyond this point, the fluctuation in ΔF values between 0.95 Å and 2.13 Å is not significant. Up to 2.12 Å, the ΔF value increases to approximately 5.03 kcal/mol. As the pulling distance extends to 3.15 Å (Figure 8-IV), the ΔF value slightly decreases, ultimately reaching 4.12 kcal/mol at 3.99 Å (Figure 8-V).

In contrast, the PMF and free energy curves for O<sub>B</sub> desorption exhibit notable disparities compared to those of O<sub>A</sub> desorption from the surface of t-ZrO<sub>2-x</sub>. Initially, the free energy values for O<sub>B</sub> desorption steadily rise, reaching a ΔF value of 8.2 kcal/mol at 0.86 Å pulling distance.
Subsequently, between 0.86 Å and 2.22 Å, the ΔF values increase with a more gradual slope, peaking at 14.69 kcal/mol at 2.22 Å. Beyond this point, up to 3.99 Å, the ΔF values fluctuate around a relatively constant value of approximately 14.39 kcal/mol. The significant contrast in free energy values between O<sub>A</sub> and O<sub>B</sub> desorption can be attributed to the differing strengths of adsorption of the oxygen atoms on the surface. According to the Figure 2-IV, O<sub>B</sub> is adsorbed at a deeper site, forming a stronger bond with the surrounding atoms, while O<sub>A</sub> is adsorbed at a greater distance from the surface. As discussed in the section 3.1, the charge transfer between Zr and O<sub>B</sub> is greater than that of Zr-O<sub>A</sub>, resulting in stronger adsorption of O<sub>B</sub>. Consequently, more work is required to desorb O<sub>B</sub> from the surface, as illustrated in the PMF profile in Figure S3. The greater work required for non-equilibrium oxygen desorption corresponds to larger values of free energy for the desorption process, rendering it less spontaneous. Therefore, it can be concluded that the desorption of O<sub>B</sub> is more non-spontaneous compared to O<sub>A</sub> due to the higher free energy values associated with its desorption process.
<figure><img src="image_9.png" /><figcaption>Figure 8. Evolution of the O₂ desorption process on the surface of t-ZrO₂-x in aqueous solution. The labeled snapshots correspond to specific points in Figure 7.</figcaption></figure>
### 3.3.2 PMF profile of O₂ desorption from the surface of t-ZrO₂-x and the hydroxylated t-ZrO₂-xNₓ

As detailed in Section 3.1, the O₂ adsorption mechanism on hydroxylated t-ZrO₂-x follows an associative pathway involving the breaking of the Oₐ-Oₛ bond. Therefore, two external force were applied to the each adsorbed Oₐ and Oₛ to facilitate its removal from the surface of the hydroxylated t-ZrO₂-x. The free energy curve and PMF profile for O₂ desorption from the hydroxylated t-ZrO₂-x surface are depicted in Figure 9 and Figure S3, respectively. In addition, the steps of Oₐ and Oₛ desorption from the surface of hydroxylated t-ZrO₂-x are illustrated in Figure 10.

As depicted in Figure 9, at small pulling distances (0-0.2 Å for Oₐ, and 0-0.3Å for Oₛ), the free energy values are negative, indicating an energetically favorable state. During these small displacements, the interaction between Oₐ, Oₛ, and the surface remains significant, keeping the free energy negative and reflecting the stabilization provided by the surface. Beyond these small distances, the free energy for both Oₐ and Oₛ steadily increases, reaching ΔF values of 3.74 kcal/mol and 1.30 kcal/mol, respectively, at a pulling distance of 0.6Å. As shown in Figure 4, the free energy curve of Oₐ is greater than that of Oₛ. This pattern indicates that more water molecules are adsorbed near Oₐ compared to Oₛ (Figure 10-I). Therefore, Oₐ interacts more with the surrounding water molecules as it moves through the solution. These interactions create more barriers for the pulling of Oₐ under external force, resulting in larger values of free energy and work needed for separating it from the surface. The ΔF value then slightly increases, reaching 5.22 kcal/mol for Oₐ and 3.09 kcal/mol for Oₛ at a pulling distance of 1.4 Å. At this pulling distance, Oₐ and Oₛ form a bond and start to desorb from the surface (Figure 10-III). Beyond this point, O₂ desorbed completely from the surface, and moved outside of the solution (Figure 10-IV). Therefore, the ΔF values start to be constant until 3.9Å, with a corresponding free energy of 7.50 kcal/mol for Oₐ and 4.50 kcal/mol for Oₛ.
<figure><img src="image_10.png" /><figcaption>Figure 9. Free energy profile based on Jarzynski's equality of O₂ desorption from the surfaces of the hydroxylated t-ZrO₂-x and the hydroxylated t-ZrO₂-xNₓ.</figcaption></figure>

In the case of t-ZrO₂-xNₓ, as discussed in Section 3.2, nitrogen-doped atoms are not stable during O₂ adsorption on the surface and tend to migrate into the solution. Consequently, this instability poses a challenge for catalytic behavior. Therefore, we chose not to investigate the PMF profile of this type of catalysts. We compared the PMF profile and subsequently free energy curve of O₂ desorption on the surface of hydroxylated t-ZrO₂-xNₓ with the hydroxylated t-ZrO₂-x, and the result is depicted in Figure S3, Figure 9, respectively. Furthermore, the evolution of O₂ desorption from the surface of hydroxylated t-ZrO₂-xNₓ over time is depicted in Figure 11.

As previously described, the O₂ molecule adsorbs associatively with a weak bond to the Zr site on the surface of hydroxylated t-ZrO₂-xNₓ (Figure 11-I'). Consequently, an external force was applied to the associatively adsorbed O₂ molecules on the surface of hydroxylated t-ZrO₂-xNₓ. First, the free energy of pulling an O₂ molecule from the surface of hydroxylated t-ZrO₂-xNₓ is negative, ranging from -0.33 to -0.15 kcal/mol at pulling distances of 0.005 to 0.1Å. As the pulling distance increases, the free energy curve and PMF continue to rise, reaching a ΔF value
of 1.55 kcal/mol at a pulling distance of 0.44 Å. Up to this point, the ΔF values for O₂ desorption from the surface of hydroxylated t-ZrO₂-xNₓ exhibit a similar trend to those from hydroxylated t-ZrO₂-x. Beyond 0.44 Å, the free energy profiles of O₂ desorption from hydroxylated t-ZrO₂-x and t-ZrO₂-xNₓ surfaces deviate significantly. The ΔF values for hydroxylated t-ZrO₂-xNₓ are lower and larger than the ΔF values for Oₐ and Oₓ desorption from hydroxylated t-ZrO₂-x surfaces. The reason for this behavior is that when force is applied to pull the associatively adsorbed O₂ in hydroxylated t-ZrO₂-xNₓ, both oxygen atoms of the O₂ molecule are subjected to the pulling force and are simultaneously pulled during the desorption process. As a result, the bond between the Oₐ atoms remains intact, making it easier to pull the molecule and requiring less work. In the range of 0.44-0.92 Å, the ΔF values remain relatively constant (Figure 11-II'). However, beyond this point, they sharply increase until reaching 2.81 Å (Figure 11-III'), with a ΔF value of 11.86 kcal/mol. This increase is attributed to the interactions with surrounding adsorbed water molecules, which hinder the pulling of the O₂ molecule, necessitating additional force for it to move in the opposite direction. This is reflected in the continuous rise in the free energy values. Upon reaching a pulling distance of 2.8 Å, there is a slight change and decrease in the PMF and free energy curve, indicating the movement of the O₂ molecule out of the solution(Figure 11-IV').
<figure><img src="image_11.png" /><figcaption>Figure 10. Evolution of the O₂ desorption process on the surface of hydroxylated t-ZrO₂-x in aqueous solution. The labeled snapshots correspond to specific points in Figure 9.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 11. Evolution of the O₂ desorption process on the surface of hydroxylated t-ZrO₂-xN_x in aqueous solution. The labeled snapshots correspond to specific points in Figure 9.</figcaption></figure>

4. Conclusion

We conducted molecular dynamics simulations to investigate the oxygen adsorption mechanisms on two catalyst surfaces, t-ZrO2-x and t-ZrO2-xNx. Our study revealed that the presence of hydroxyl groups significantly influenced the stability and pathways of O2 adsorption, determining whether the process occurs through associative or dissociative mechanisms. Our results demonstrated that O2 molecules adsorbed dissociatively on t-ZrO2-x at
two different Zr sites. In contrast, on hydroxylated t-ZrO2-x, prolonged interaction between O2 and the OH group led to associative O2 adsorption near the oxygen vacancy and hydroxyl group, with the O-O bond breaking. Furthermore, O2 adsorption on t-ZrO2-xNx occurred associatively. However, during O2 adsorption on t-ZrO2-xNx, N dopants tend to leave the surface, which is unfavorable for catalytic activity. The hydroxyl group played a crucial role in modifying the t-ZrO2-xNx surface, stabilizing N dopants, and preventing their loss. Consequently, the O2 adsorption mechanism on both t-ZrO2-xNx and hydroxylated t-ZrO2-xNx occurred associatively. Furthermore, a series of MD simulations were employed to utilize the steered molecular dynamics (SMD) method, facilitating the construction of a model for oxygen detachment from catalyst surfaces. Subsequently, we performed PMF analysis using Jarzynski's equality, derived from the SMD simulations, to assess the impact of defect nature and the presence of hydroxyl groups on the free energy values associated with oxygen desorption from different catalysts. The findings indicated that the free energy linked to oxygen desorption from the examined surfaces is the primary driving force behind the spontaneous process. Moreover, the presence of the hydroxyl group was found to notably influence the free energy of the O2 desorption process. Oxygen desorption from hydroxylated surfaces was observed to be lower than from surfaces lacking hydroxyl groups, indicating a more spontaneous process. Besides, the PMF profile of the hydroxylated t-ZrO2-x showed that less work is required to pull O2 from its surface compared to hydroxylated t-ZrO2-xNx. This implies that the increased spontaneity may lead to higher reactivity and an accelerated ORR compared to hydroxylated t-ZrO2-xNx.

## ASSOCIATED CONTENT

### Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### Supporting Information

The supporting information includes the following details: Third-generation COMB (COMB3) potential: A comprehensive explanation of the COMB3 potential used in the simulations. Structural details of a unit cell of t-ZrO2(001): Description and specifications of the tetragonal t-ZrO2(001) unit cell structure. Exhaustive search of optimized defective t-ZrO2-xNx: Detailed procedures and results from the optimization search for defective nitrogen-doped t-ZrO2.
Comparative PMF profiles of O₂ desorption: A figure showing the PMF profiles of O₂ desorption from the surfaces of t-ZrO₂-x, t-ZrO₂-xNₓ, hydroxylated t-ZrO₂-x, and hydroxylated t-ZrO₂-xNₓ. Changes in Oₐ-Oₛ bond length: Two tables detailing the changes in the Oₐ-Oₛ bond length during O₂ adsorption on: t-ZrO₂-x, and hydroxylated t-ZrO₂-x

# AUTHOR INFORMATION

## Corresponding Authors

**Sara Fazeli**- GREMI, CNRS - Université d'Orléans, 45067 ORLEANS Cedex 2, France,
Simulation moléculaire pour tous - MS4ALL, 1 Avenue du champ de Mars, 45100 Orléans,
France

Email: [seyedehsara.fazeli@univ-orleans.fr](mailto:seyedehsara.fazeli@univ-orleans.fr), [sara.fazeli@ms4all.eu](mailto:sara.fazeli@ms4all.eu)

**Pascal Brault**- GREMI, CNRS - Université d'Orléans, 45067 ORLEANS Cedex 2, France,
Simulation moléculaire pour tous - MS4ALL, 1 Avenue du champ de Mars, 45100 Orléans,
France

Email: [pascal.brault@univ-orleans.fr](mailto:pascal.brault@univ-orleans.fr), [pascal.brault@ms4all.eu](mailto:pascal.brault@ms4all.eu)

## Authors

**Amaël Caillard**- GREMI, CNRS - Université d'Orléans, 45067 ORLEANS Cedex 2, France,
Email: [amael.caillard@univ-orleans.fr](mailto:amael.caillard@univ-orleans.fr)

**Eric Million**- GREMI, CNRS - Université d'Orléans, 45067 ORLEANS Cedex 2, France,
Email: [eric.millon@univ-orleans.fr](mailto:eric.millon@univ-orleans.fr)

## Acknowledgements

This work has been funded by the French National Research Agency (ANR) through the InnOxiCat project ANR-20-CE05-0010.

## References

(1) Wang, J.; Azam, W. Natural Resource Scarcity, Fossil Fuel Energy Consumption, and Total Greenhouse Gas Emissions in Top Emitting Countries. *Geosci. Front.* **2024**, 15
(2), 101757.

(2) Gong, J.; Luque, R. Catalysis for Production of Renewable Energy. *Chem. Soc. Rev.* **2014**, *43* (22), 7466–7468.

(3) Chen, Z.; Wei, W.; Ni, B.-J. Cost-Effective Catalysts for Renewable Hydrogen Production via Electrochemical Water Splitting: Recent Advances. *Curr. Opin. Green Sustain. Chem.* **2021**, *27*, 100398.
https://doi.org/https://doi.org/10.1016/j.cogsc.2020.100398.

(4) Ronovsky, M.; Pan, L.; Klingenhof, M.; Martens, I.; Fusek, L.; Kus, P.; Chattot, R.; Mirolo, M.; Dionigi, F.; Burdett, H. Assessing Utilization Boundaries for Pt-Based Catalysts in an Operating Proton-Exchange Membrane Fuel Cell. *ACS Appl. Energy Mater.* **2023**, *6* (17), 8660–8665.

(5) Zhao, W.; Xu, G.; Dong, W.; Zhang, Y.; Zhao, Z.; Qiu, L.; Dong, J. Progress and Perspective for In Situ Studies of Oxygen Reduction Reaction in Proton Exchange Membrane Fuel Cells. *Adv. Sci.* **2023**, 2300550.

(6) Tripkovic, V.; Vegge, T. Potential-and Rate-Determining Step for Oxygen Reduction on Pt (111). *J. Phys. Chem. C* **2017**, *121* (48), 26785–26793.

(7) Wang, S.; Zhu, E.; Huang, Y.; Heinz, H. Direct Correlation of Oxygen Adsorption on Platinum-Electrolyte Interfaces with the Activity in the Oxygen Reduction Reaction. *Sci. Adv.* **2021**, *7* (24), eabb1435.

(8) He, C.; Sankarasubramanian, S.; Matanovic, I.; Atanassov, P.; Ramani, V. Understanding the Oxygen Reduction Reaction Activity and Oxidative Stability of Pt Supported on Nb-Doped TiO2. *ChemSusChem* **2019**, *12* (15), 3468–3480.

(9) Viswanathan, V.; Hansen, H. A.; Rossmeisl, J.; Nørskov, J. K. Unifying the 2e- and 4e-Reduction of Oxygen on Metal Surfaces. *J. Phys. Chem. Lett.* **2012**, *3* (20), 2948–2951.

(10) Yamamoto, Y.; Kasamatsu, S.; Sugino, O. Scaling Relation of Oxygen Reduction Reaction Intermediates at Defective TiO2 Surfaces. *J. Phys. Chem. C* **2019**, *123* (32), 19486–19492.

(11) Li, R.; Li, H.; Xu, S.; Liu, J. Theoretical Investigation of the Effect of OH- Ions on O2 Adsorption on Low-Index Pt Surfaces in Alkaline Solution. *Appl. Surf. Sci.* **2015**, *351*, 853–861. https://doi.org/https://doi.org/10.1016/j.apsusc.2015.06.042.

(12) Lankiang, S.; Chiwata, M.; Baranton, S.; Uchida, H.; Coutanceau, C. Oxygen Reduction Reaction at Binary and Ternary Nanocatalysts Based on Pt, Pd and Au. *Electrochim. Acta* **2015**, *182*, 131–142.
https://doi.org/10.1016/j.electacta.2015.09.061.

(13) Ishihara, A.; Nagai, T.; Ukita, K.; Arao, M.; Matsumoto, M.; Yu, L.; Nakamura, T.; Sekizawa, O.; Takagi, Y.; Matsuzawa, K.; Napporn, T. W.; Mitsushima, S.; Uruga, T.; Yokoyama, T.; Iwasawa, Y.; Imai, H.; Ota, K. Emergence of Oxygen Reduction Activity in Zirconium Oxide-Based Compounds in Acidic Media: Creation of Active Sites for the Oxygen Reduction Reaction. *J. Phys. Chem. C* **2019**, *123* (30), 18150–18159. https://doi.org/10.1021/acs.jpcc.9b02393.

(14) Zorko, M.; Farinazzo Bergamo Dias Martins, P.; Connell, J. G.; Lopes, P. P.; Markovic, N. M.; Stamenkovic, V. R.; Strmcnik, D. Improved Rate for the Oxygen Reduction Reaction in a Sulfuric Acid Electrolyte Using a Pt (111) Surface Modified with Melamine. *ACS Appl. Mater. Interfaces* **2021**, *13* (2), 3369–3376.

(15) Li, H.; Zhao, H.; Tao, B.; Xu, G.; Gu, S.; Wang, G.; Chang, H. Pt-Based Oxygen Reduction Reaction Catalysts in Proton Exchange Membrane Fuel Cells: Controllable Preparation and Structural Design of Catalytic Layer. *Nanomaterials* **2022**, *12* (23), 4173.

(16) Huang, Y.; Li, X.; Shen, P. K.; Yang, C.; Zhu, J. High-Performance MoP-Mo2C/C Heterogeneous Nanoparticle Catalysts for Alkaline Hydrogen Evolution and Oxidation Reactions. *ACS Mater. Lett.* **2024**, *6* (5), 1678–1685.

(17) Liu, J.; Fung, V.; Wang, Y.; Du, K.; Zhang, S.; Nguyen, L.; Tang, Y.; Fan, J.; Jiang, D.; Tao, F. F. Promotion of Catalytic Selectivity on Transition Metal Oxide through Restructuring Surface Lattice. *Appl. Catal. B Environ.* **2018**, *237*, 957–969.

(18) Mladenović, D.; Mladenović, A.; Santos, D. M. F.; Yurtcan, A. B.; Miljanić, Š.; Mentus, S.; Šljukić, B. Transition Metal Oxides for Bifunctional ORR/OER Electrocatalysis in Unitized Regenerative Fuel Cells. *J. Electroanal. Chem.* **2023**, *946*, 117709. https://doi.org/10.1016/j.jelechem.2023.117709.

(19) Wolf, S.; Roschger, M.; Genorio, B.; Garstenauer, D.; Hacker, V. Mixed Transition-Metal Oxides on Reduced Graphene Oxide as a Selective Catalyst for Alkaline Oxygen Reduction. *ACS omega* **2023**, *8* (12), 11536–11543.

(20) Suntivich, J.; Gasteiger, H. A.; Yabuuchi, N.; Nakanishi, H.; Goodenough, J. B.; Shao-Horn, Y. Design Principles for Oxygen-Reduction Activity on Perovskite Oxide Catalysts for Fuel Cells and Metal–Air Batteries. *Nat. Chem.* **2011**, *3* (7), 546–550.

(21) Hofmann, A.; Clark, S. J.; Oppel, M.; Hahndorf, I. Hydrogen Adsorption on the Tetragonal ZrO 2 (101) Surface: A Theoretical Study of an Important Catalytic Reactant. *Phys. Chem. Chem. Phys.* **2002**, *4* (14), 3500–3508.
(22) Fazeli, S.; Brault, P.; Caillard, A.; Thomann, A.-L.; Millon, E.; Atmane, S.; Coutanceau, C. Atomistic Insight into the Defect-Induced Tunable Plasticity and Electronic Properties of Tetragonal Zirconia. *Mater. Today Commun.* **2023**, *36*, 106731.

(23) Qiu, H.; Zhang, Y.; Huang, W.; Peng, J.; Chen, J.; Gao, L.; Omran, M.; Li, N.; Chen, G. Sintering Properties of Tetragonal Zirconia Nanopowder Preparation of the NaCl + KCl Binary System by the Sol–Gel–Flux Method. *ACS Sustain. Chem. Eng.* **2023**, *11* (3), 1067–1077. https://doi.org/10.1021/acssuschemeng.2c05908.

(24) Xia, X.; Oldman, R. J.; Catlow, C. R. A. Oxygen Adsorption and Dissociation on Yttria Stabilized Zirconia Surfaces. *J. Mater. Chem.* **2012**, *22* (17), 8594–8612.

(25) Zheng, X.; Yu, P.; Liu, Y.; Ma, Y.; Cao, Y.; Cai, Z.; Zhou, L.; Huang, K.; Zheng, S.; Jiang, L. Efficient Hydrogenation of Methyl Palmitate to Hexadecanol over Cu/m-ZrO2 Catalysts: Synergistic Effect of Cu Species and Oxygen Vacancies. *ACS Catal.* **2023**, *13* (3), 2047–2060.

(26) Yang, X.; Chen, L.; Yu, X.; Yu, J.; Han, D.; Xiao, M.; Ge, M.; Yang, X. Exploration of Atomic Interfaces with Inherent Oxygen Vacancies in Zirconia for Toluene Oxidation. *J. Mater. Chem. A* **2023**, *11* (1), 287–296.

(27) Livraghi, S.; Chierotti, M. R.; Giamello, E.; Magnacca, G.; Paganini, M. C.; Cappelletti, G.; Bianchi, C. L. Nitrogen-Doped Titanium Dioxide Active in Photocatalytic Reactions with Visible Light: A Multi-Technique Characterization of Differently Prepared Materials. *J. Phys. Chem. C* **2008**, *112* (44), 17244–17252.

(28) Cao, X.; Zheng, S.; Wang, T.; Lin, F.; Li, J.; Jiao, L. N-Doped ZrO2 Nanoparticles Embedded in a N-Doped Carbon Matrix as a Highly Active and Durable Electrocatalyst for Oxygen Reduction. *Fundam. Res.* **2022**, *2* (4), 604–610.

(29) Muhammady, S.; Haruyama, J.; Kasamatsu, S.; Sugino, O. Effect of Nitrogen Doping and Oxygen Vacancy on the Oxygen Reduction Reaction on the Tetragonal Zirconia (101) Surface. *J. Phys. Chem. C* **2022**, *126* (37), 15662–15670.

(30) Jian, M.; Liu, J.-X.; Li, W.-X. Hydroxyl Improving the Activity, Selectivity and Stability of Supported Ni Single Atoms for Selective Semi-Hydrogenation. *Chem. Sci.* **2021**, *12* (30), 10290–10298.

(31) Añez, R.; Sierraalta, A.; Martorell, G.; Sautet, P. Stabilization of the (1 1 0) Tetragonal Zirconia Surface by Hydroxyl Chemical Transformation. *Surf. Sci.* **2009**, *603* (16), 2526–2531.

(32) Fazeli, S.; Brault, P.; Caillard, A.; Thomann, A.-L.; Millon, E.; Coutanceau, C.;
Atmane, S. Influence of Defects on Barrier Energy Formation for OOH\* Intermediate in ORR on Tetragonal-ZrO2 with Adsorbed-Hydroxyl. *J. Phys. Chem. C* **2024**, 128 (28), 11572–11582. https://doi.org/10.1021/acs.jpcc.4c01144.

(33) Yang, F.; Liao, L.; Li, M.; Mei, D.; Chen, Y. Kinetics Study on O2 Adsorption and OHad Desorption at Pt (111), Its Implication to Oxygen Reduction Reaction Kinetics. *Chinese J. Chem. Phys.* **2014**, 27 (4), 479–484.

(34) Montemore, M. M.; Van Spronsen, M. A.; Madix, R. J.; Friend, C. M. O2 Activation by Metal Surfaces: Implications for Bonding and Reactivity on Heterogeneous Catalysts. *Chem. Rev.* **2017**, 118 (5), 2816–2862.

(35) Yim, W.-L.; Klüner, T. Promoting O2 Activation on Noble Metal Surfaces. *J. Catal.* **2008**, 254 (2), 349–354. https://doi.org/10.1016/j.jcat.2008.01.018.

(36) Tamtaji, M.; Li, Y.; Cai, Y.; Liu, H.; Goddard III, W. A.; Chen, G. Reaction Mechanism and Kinetics of Oxygen Reduction Reaction on the Iron–Nickel Dual Atom Catalyst. *J. Mater. Chem. A* **2023**, 11 (46), 25410–25421.

(37) Song, K.; Yang, B.; Zou, X.; Zhang, W.; Zheng, W. T. Unified ORR Mechanism Criteria via Charge–Spin–Coordination of Fe Functional Units. *Energy Environ. Sci.* **2023**.

(38) Hu, Y.; Sun, J.; Wei, H.; Ai, M.; Li, Z. Adsorption Characteristics and Oxygen Reduction Reactions on Pristine and Pt-, Co-Decorated Antimonenes: A DFT-D Study. *New J. Chem.* **2020**, 44 (3), 1138–1146.

(39) Li, B.; Holby, E. F.; Wang, G. Mechanistic Insights into Metal, Nitrogen Doped Carbon Catalysts for Oxygen Reduction: Progress in Computational Modeling. *J. Mater. Chem. A* **2022**, 10 (45), 23959–23972.

(40) Fang, C.; Zhou, J.; Zhang, L.; Wan, W.; Ding, Y.; Sun, X. Synergy of Dual-Atom Catalysts Deviated from the Scaling Relationship for Oxygen Evolution Reaction. *Nat. Commun.* **2023**, 14 (1), 4449.

(41) Singh, S. K.; Takeyasu, K.; Nakamura, J. Active Sites and Mechanism of Oxygen Reduction Reaction Electrocatalysis on Nitrogen-doped Carbon Materials. *Adv. Mater.* **2019**, 31 (13), 1804297.

(42) Li, R.; Wang, L.; Yue, Q.; Li, H.; Xu, S.; Liu, J. Insights into the Adsorption of Oxygen and Water on Low-Index Pt Surfaces by Molecular Dynamics Simulations. *New J. Chem.* **2014**, 38 (2), 683–692.

(43) Sokolović, I.; Reticcioli, M.; Čalkovský, M.; Wagner, M.; Schmid, M.; Franchini, C.; Diebold, U.; Setvín, M. Resolving the Adsorption of Molecular O2 on the Rutile TiO2
(110) Surface by Noncontact Atomic Force Microscopy. *Proc. Natl. Acad. Sci. 2020*, 117 (26), 14827-14837.

(44) Binninger, T.; Kowalski, P. M.; Eikerling, M. H. Oxygen Desorption – Critical Step for the Oxygen Evolution Reaction. *Curr. Opin. Electrochem. 2023*, *42*, 101382.
https://doi.org/https://doi.org/10.1016/j.coelec.2023.101382.

(45) Deng, X.; Min, B. K.; Guloy, A.; Friend, C. M. Enhancement of O2 Dissociation on Au (111) by Adsorbed Oxygen: Implications for Oxidation Catalysis. *J. Am. Chem. Soc.* **2005**, *127* (25), 9267-9270.

(46) Badar, M. S.; Shamsi, S.; Ahmed, J.; Alam, M. A. Molecular Dynamics Simulations: Concept, Methods, and Applications. In *Transdisciplinarity*; Springer, 2022; pp 131–151.

(47) Wang, Y.; Yan, H.-M.; Cao, H.; Chen, J.; Yang, H.; Zhu, J.; Sun, J. Atomic-Scale Understanding of Electrified Interfacial Structures and Dynamics during the Oxygen Reduction Reaction on the Fe-N4/C Electrocatalyst. *ACS Catal.* **2023**, *13* (16), 11080-11090.

(48) You, W.; Tang, Z.; Chang, C.-E. A. Potential Mean Force from Umbrella Sampling Simulations: What Can We Learn and What Is Missed? *J. Chem. Theory Comput.* **2019**, *15* (4), 2433-2443.

(49) Raman, S.; Utzig, T.; Baimpos, T.; Ratna Shrestha, B.; Valtiner, M. Deciphering the Scaling of Single-Molecule Interactions Using Jarzynski's Equality. *Nat. Commun.* **2014**, *5* (1), 5539.

(50) Park, S.; Khalili-Araghi, F.; Tajkhorshid, E.; Schulten, K. Free Energy Calculation from Steered Molecular Dynamics Simulations Using Jarzynski's Equality. *J. Chem. Phys.* **2003**. https://doi.org/10.1063/1.1590311.

(51) Liang, T.; Shin, Y. K.; Cheng, Y.-T.; Yilmaz, D. E.; Vishnu, K. G.; Verners, O.; Zou, C.; Phillpot, S. R.; Sinnott, S. B.; Van Duin, A. C. T. Reactive Potentials for Advanced Atomistic Simulations. *Annu. Rev. Mater. Res.* **2013**, *43*, 109-129.

(52) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. *J. Comput. Phys.* **1995**. https://doi.org/10.1006/jcph.1995.1039.

(53) Fraternali, F. Restrained and Unrestrained Molecular Dynamics Simulations in the NVT Ensemble of Alamethicin. *Biopolym. Orig. Res. Biomol. 1990*, *30* (11-12), 1083-1099.

(54) Evans, D. J.; Holian, B. L. The Nose-Hoover Thermostat. *J. Chem. Phys.* **1985**, *83* (8), 4069-4074.
<table><tr><td>Abstract</td><td>1</td></tr><tr><td>1. Introduction</td><td>2</td></tr><tr><td>2. Research Methodology</td><td>4</td></tr><tr><td>3. Results and Discussions</td><td>7</td></tr><tr><td>3.1. O<sub>2</sub> adsorption on t-ZrO<sub>2-x</sub> and the hydroxylated t-ZrO<sub>2-x</sub></td><td>7</td></tr><tr><td>3.2. O<sub>2</sub> adsorption on t-ZrO<sub>2-x</sub>N<sub>x</sub> and the hydroxylated t-ZrO<sub>2-x</sub>N<sub>x</sub></td><td>13</td></tr><tr><td>3.3. Free energy calculation from SMD simulation</td><td>18</td></tr><tr><td>3.3.1. PMF profile of O<sub>2</sub> desorption from the surface of t-ZrO<sub>2-x</sub></td><td>18</td></tr><tr><td>3.3.2. PMF profile of O<sub>2</sub> desorption from the surface of t-ZrO<sub>2-x</sub> and the hydroxylated t-ZrO<sub>2-x</sub>N<sub>x</sub></td><td>22</td></tr><tr><td>4. Conclusion</td><td>26</td></tr><tr><td>ASSOCIATED CONTENT</td><td>27</td></tr><tr><td>Acknowledgements</td><td>28</td></tr><tr><td>References</td><td>28</td></tr></table>