Surface oxidation mechanism of a refractory high-entropy alloy

High-entropy alloys (HEAs) synthesized using refractory elements are being strongly considered as candidates for high temperature structural applications. The role of compositional changes of HEA surfaces due to oxidation is crucial to sustain the material properties, but a detailed description of the thermodynamic mechanism driving the adsorption of oxygen on such complex surfaces is absent. We examine and explain the reaction process of oxygen on a representative refractory HEA surface using first principles and atomistic thermodynamic models. The HEA surface is highly reactive to oxygen yielding a full monolayer coverage at temperatures between 300 and 1500 K. The preferential adsorption of oxygen to specific sites of the HEA surface is attributed to the electronic configuration of the bonding shells of the constituent surface atoms. On further oxygen addition, the oxygen atoms diffuse into the bulk regions of the alloy. Manipulation of temperature and oxygen pressure reveals that it is difficult to rid the alloy surface of oxygen even at extremely low pressures of 10−9 bar at 2000 K.


INTRODUCTION
High-entropy alloys (HEAs), a subset of multi-principal element alloys (MPEAs), contain N principal elements typically in near equiatomic concentrations with N ≥ 5. 1 For the majority of the HEAs investigated to date, the predominant phases formed are disordered solid solutions primarily as face-centered cubic (FCC) or body-centered cubic (BCC) crystals, and some hexagonal-closed packed (HPC) phases. These disordered phases are favored over the ordered phases due to many competing parameters such as the enthalpy of mixing, the valence electron concentration, atomic radii, lattice structure and electronegativity of the constituent elements as well as the high ideal configurational entropy of mixing, ΔS conf = − R∑X i ln(X i ) that is equal to or greater than 1.5R, 2,3 where X i represents the atomic fraction of element i and R is the molar gas constant. This high entropy of mixing reduces the driving forces for the formation of intermetallic phases. 4 However, the intermetallic phases formed in some cases are due to a considerable atomic size mismatch and large negative enthalpies of mixing of the alloying elements. 5,6 Ordered solid-solutions have, nevertheless, been observed to exist in some HEAs though they are typically less stable, which affects their properties and microstructures, e.g., through increase in density, higher melting temperatures and increased brittleness of the resulting alloys. 3,[7][8][9] New structural alloys that can survive extreme temperatures with good oxidation resistance and mechanical properties at elevated temperatures are required in applications such as highspeed turbines and in the aerospace industry. Refractory HEAs offer promising prospects for such applications. Several potential refractory HEAs 7,10-14 have been recently identified based on contemporary approaches that involve manipulation of valence electron concentrations, entropy of mixing, the Hume-Rothery rules among others. 1,15,16 These HEAs have been reported to provide promising alternatives to Ni-based superalloys. Typically, these alloys are composed of three or less phases, with the BCC phase being predominant as the primary phase. 4 Limited literature on the degradation mechanisms of some HEA surfaces due to oxidation and other physical effects are available, 17 while the corrosion resistance of several HEAs has been determined both experimentally and computationally primarily focusing on anodic and cathodic reactions. 18,19 However, to the best of our knowledge, this report is the first study that describes the fundamental thermodynamic mechanisms analyzing how oxygen adsorbs on complex surfaces, such as those of HEAs.
Recently, Singh et al. 20 proposed refractory Mo-W-Ta-Ti-Zr HEAs with compositions that exhibited greatly enhanced modulus of elasticity (3 × at 300 K) over the near equiatomic counterparts. Two compositions were found to exhibit mechanical properties better than commercial Mo-based alloys (such as TZM) with a 2.3 × higher modulus of elasticity at 2000 K: (Mo 1-z W z ) 0.85 Ta 0.10 (TiZr) 0.05 where the alloy identified as c10 (here referred to as HEA1) has z = 0.5 and another alloy denoted as c10 has z = 0.05 (denoted here as HEA2). The 0.85 signifies that Mo and W together occupy 85% by mole fraction of the alloy. The z value of 0.5 indicates each of Mo and W is 0.425 mole fraction in that HEA, while z = 0.05 indicates Mo is 0.8075 and W is 0.0425 by mole fraction in the corresponding alloy. As these refractory alloys are potentially useful for operations in high temperature environments, knowledge of their oxidation behavior is required to develop oxidationresistant HEAs. The data on the oxidation mechanisms of HEAs is sparse, and more so for the Mo-W-Ta-Ti-Zr refractory HEA.
Here, we employ density functional theory (DFT) calculations and thermodynamic modeling to understand and explain the oxidation behavior of HEA1 and HEA2 alloys in oxidizing environments. The special quasi-random structures (SQS) method is used to generate random solid solutions for the chosen HEA1 and HEA2 compositions. Different surface models are generated from the constructed random solid solution and the model with the lowest surface energy used. Sequential adsorption of oxygen on the alloy surface is carried out and characterized by the Gibbs free energy of adsorption until full coverage of the surface is obtained. Subsequently, we examine the diffusion of oxygen into the bulk regions of these alloys.

RESULTS AND DISCUSSION
Alloy structure properties Upon relaxing the unit cell of the generated SQS structures, we note that the changes in cell vectors relative to the ideal unrelaxed SQS structures are minimal indicating stability of the produced random solid solution. The pair correlation functions for the alloys are reproduced perfectly for the SQS-40, SQS-80, SQS-120, and SQS-160 structures, where SQS-N denotes structures consisting of N atoms. All the generated structures are monoclinic. Figure 1 shows the ideal unrelaxed structures created using the mcsqs code.
The calculated formation enthalpies for all the generated SQS-N structures are presented in Fig. 2. As the modeled structures contain considerable number of atoms, the calculated formation enthalpies are all negative or closer to zero highlighting the stability in the generated structures. Thus, we do not expect clustering or segregation of the atoms in either the HEA1 or the HEA2 alloy.
The lattice parameter a calculated for all the SQS-N structures for both HEA1 and HEA2 alloys are listed in Table 1. The estimated lattice parameters decrease as the SQS size increases with convergence achieved when 80 or more atoms are used in generating the random solid solution. The calculated values are in good agreement with the computed value of 3.171 Å by Singh et al. 20 using Korringa-Kohn-Rostoker Coherent-Potential Approximation. As the calculated lattice parameter converges at SQS-120, we used the 120-atom structure for further calculations. This selection also ensures that we have a good representation of all atoms on the surface slab for the oxidation studies due to the low concentration of Ti and Zr atoms in the alloys.
The total (TDOS) and projected (PDOS) density of states for SQS-120 are illustrated in Fig. 3. The Fermi level lies in the pseudo-gap of HEA1, which reflects the stable nature of the bulk alloy obtained through SQS. However, from the DOS, we note that a considerable number of states are present at the Fermi level, most of which are contributed by Mo and W -d states due to the high relative concentrations of Mo and W atoms in the alloy. Additionally, there are contributions from Mo -p, Ta -d, W -p states at the Fermi level. These predictions suggest that Mo and W atoms are reactive towards oxygen upon adsorption due to the higher percentage of these atoms compared to the other elements. Thus, the presence of more Mo and W active sites compared to Ta, Ti and Zr atoms will lead to higher oxygen coverage at Mo and W rich regions.
Generation of appropriate surface structures The calculated cleavage energies for the (001), (110), and (111) surfaces of HEA1 are listed in Table 2. The (001) face that has the lowest surface energy is selected for the oxidation studies. Due to the random nature of the alloy structure, different surface terminations are possible and hence surface energies are calculated for the different (001) terminations ( Table 2). The    Fig. S1. All four surfaces yield lower surface energies than the (110) and (111) surfaces, while surface I assumes the lowest. Note that the unit for surface energies is meV.Å −2 , and hence the effective difference in stability between the four surfaces is minimal. In order to observe the effect of oxygen on all the constituent atoms of the HEA, surface termination IV that contains all five atoms of the alloy is used for the adsorption studies. As the difference between HEA1 and HEA2 is the relative fraction of Mo and W atoms (the total number of Mo and W atoms remain the same in both alloys), we assume that the compositions of the most stable HEA1 and HEA2 (001) surfaces are similar (minor difference in the number of Mo and W atoms) and the difference in their surface energies is negligible (in meV.Å −2 ). Hence, the subsequent oxidation studies are performed on HEA1 with similar observations expected for HEA2. Sequential oxidation of the (001) surface While oxygen is adsorbed on the (001) surface of HEA1 in both associative (O 2 ) and dissociative (O) modes, our results suggest a preference for dissociative adsorption of oxygen atoms. A strong interaction of the first O atom with the surface is characterized by a high adsorption energy. The O atoms show preferential bonding at the threefold fcc hollow sites between three surface metal atoms. The first O atom adsorbs to two Zr and one Ti surface atoms. The calculated bond distances between the O atom and the two Zr atoms are 2.108 Å and 2.175 Å, while that between the O and Ti atoms is 1.960 Å. Figure 4 illustrates the adsorption of oxygen on the HEA1 (001) surface at various coverages.
The calculated adsorption energy (E ads ) is 5.45 eV. An earlier investigation on oxidation of ZrC found O to adsorb strongly at fcc three-fold hollow sites mostly between surface Zr atoms with E ads = 5.83 eV on the (111) surface that primarily had only Zr atoms. 21 The high adsorption energy predicted from our computations indicates a high reactivity of oxygen towards the HEA1 alloy surface.
We find that the first four added oxygen atoms prefer to bond at sites with available Zr atoms with successive adsorption energies of 5.45, 5.01, 4.84 and 5.10 eV respectively. The first O atom with the highest E ads involves Zr and Ti while the next two highest E ads (5.01 and 5.10 eV) involve Zr and Ta atoms. The lowest E ads of 4.84 eV involves Zr with two Mo atoms. Thus, a high affinity of oxygen to both Zr atoms is noted followed by Ti, Ta, and Mo though the percentage of Zr, Ti, and Ta active sites are lesser than the Mo and W active sites present on the surface.
Both Zr and Ti have the highest number of unfilled d-orbitals (five unfilled d-orbitals), including two partially filled orbitals compared to Mo, W, Ta. This electronic structure results in the There is further rearrangement upon adsorption of the 18th O atom, which is located at a fourfold hollow site between a Zr, two Mo and a W atom with E ads = 3.78 eV. The 19th O atom is at a site between two Mo atoms and one W atom with E ads = 3.63 eV. The 20th O atom adsorbs between two W atoms and one Mo atom with E ads = 3.11 eV. When the 21st O atom adsorbs at a site between a Ti atom and two Mo atoms with E ads of 3.49 eV, there is rearrangement of the 5th O atom (between Ti, W, and Mo atoms) as well as the 1st O atom (which moves further beneath the surface plane by forming a bond with a Ti atom in the second layer). There is further rearrangement when the 22nd O atom adsorbs between three Mo atoms with E ads = 3.73 eV. The 23rd O atom adsorbs between two Mo atoms and one W atom with E ads = 3.47 eV. Upon adsorption of the 24th O atom, the 1st O atom moves further inwards towards the second atomic layer whiles the 24th atom is located atop the 1st O atom. Thus, as the full monolayer or ML (here, defined as occupation of all surfaceactive sites) coverage is achieved, oxygen begins to move to the bulk propagating through the region rich in Zr and Ti atoms.
The mean adsorption energies E ads as well as the successive adsorption energies for all other coverages are provided in Table 3. In general, the average adsorption energy as well as the successive adsorption energies for the O atoms decrease with increasing coverage. The strongest O interactions are observed in Zr rich regions of the alloy.
The lowest successive adsorption energy of 3.11 eV is noted for the adsorption of the 20th O atom. This step involves the oxygen atom occupying a three-fold fcc hollow site between a Mo and two W surface atoms. Mo has all 4d orbitals partially filled and energetically stable compared to the other elements, which leads to low affinity for oxygen and results in the lowest successive adsorption energy. Also, it is important to note that W has more d orbitals partially filled than Ta. Thus, the weakest interactions of O atoms are with W and Mo atoms, but given their higher mole fraction relative to other elements, the corresponding regions contribute to enhanced oxygen coverage on the HEA surface. The next lowest successive adsorption energy of 3.24 eV is recorded during the adsorption of the 24th O atom above the site of the 1st O atom by pushing the latter into the sub-surface plane. Thus, the initial region at which surface O atoms begin to diffuse from the surface layer into the bulk is the subdomain rich in Zr and Ti atoms as shown in Fig. 4.
We further examine the diffusion of oxygen into the bulk structure of the HEA1 alloy by adding a 25th O atom to the fully oxidized surface layer. We observed an exothermic reaction for such process with E ads = 2.60 eV, which is lower than the adsorption energies for all surface processes. We corroborate that complete surface coverage is achieved before diffusion of any oxygen atoms into the bulk structure of the alloy. The 25th O atom adsorbs at the region where the 1st O atom initially adsorbed. Upon addition of the 25th O atom, the 1st O atom rearranges to occupy a spot on the plane of the first surface layer of atoms (where the 13th O atom is initially adsorbed). The 25th O atom is pushed further down into the inter-planar region between the first and second sub-surface layers and diffuses into the bulk (shown in Fig. 5). This atom sits in the planar region of the second subsurface layer at a hollow site between Mo, W, Ta, Ti atoms and directly atop a Ta atom in the third sub-surface layer.
Atomistic thermodynamic modeling We investigate the effects of temperature and pressure on the alloys by thermodynamic modeling using stability plots for the HEA when exposed surfaces are oxidized. The Gibbs free energy of adsorption per surface area is determined as a function of pressure for different temperature regimes. As shown in Fig. 6 for representative oxygen coverages at lower to mid-range E. Osei-Agyemang and G. Balasubramanian temperatures (300-1500 K), the surface is covered with a monolayer of oxygen. The variations of temperature and pressure are not sufficient to eliminate oxygen from the alloy surface at low to medium temperatures. To remove some of the adsorbed oxygen atoms from the HEA1 surface, a higher temperature of~2000 K is selected. The distribution of the Gibbs free energy against pressure at 2000 K is presented in Fig. 7 (Figure inset  Between the pressures of 10 2 -10 −4 bar, the (001) surface of HEA1 is fully covered (1 ML) with oxygen. As the oxygen pressure is reduced from 10 −4 to 10 −5 bar at 2000 K, there is the removal of one O atom with the 0.958 ML coverage becoming the most stable. Between 10 −5 and 10 −5.8 bar, there is further removal of four oxygen atoms, yielding a coverage of 0.792 ML on the surface. As the oxygen pressure is further reduced from 10 −5.8 to 10 −6.6 bar, one more oxygen atom is removed from the surface resulting in 0.750 ML coverage. Reducing from 10 −6.6 to 10 −7 bar of oxygen pressure, there is further removal of one oxygen atom to yield 0.708 ML coverage. With further lower oxygen pressure regimes, another oxygen atom is removed from 10 −7 to 10 −7.6 bar to yield 0.667 ML coverage. Between 10 −7.6 and 10 −7.8 bar, four oxygen atoms are further removed from the surface, yielding approximately. 0.500 ML coverage. From 10 −7.8 to 10 −8.4 bar oxygen pressure, another O atom is eliminated from the surface to yield 0.458 ML coverage, while from 10 −8.4 to 10 −8.6 bar pressure reduction, two oxygen atoms are released from the surface, which yields a 0.375 ML coverage. Between the extremely low oxygen pressures of 10 −8.6 -10 −9 bar, an additional oxygen atom gets removed to form 0.333 ML coverage. From this analysis, it is evident that even at extremely low oxygen pressures of 10 −9 bar at 2000 K, a significant number of oxygen atoms are present on the surface (0.333 ML). Thus, temperature and pressure variations are not sufficient to retrieve the bare alloy surface once oxidation sets in.
The oxidation process on the (001) surface of the refractory [(Mo 1−z W z ) 0.85 Ta 0.10 (TiZr) 0.05 (z = 0.5)] HEA alloy is examined. The SQS method is used to generate random solid solution structures with the desired composition while surface energies of different terminations are used to select an appropriate alloy surface for oxygen adsorption. Sequential adsorption of oxygen atoms at the available alloy active sites suggests a very strong interaction of the HEA with oxygen. We note a high affinity of oxygen to Zr atoms followed by Ti, Ta and Mo. The high affinity of oxygen to Zr and Ti atoms compared to the other elements is due to the electronic configuration of the bonding shells. A complete coverage with oxygen occurs for the HEA surface as long as there is the Fig. 5 Side view of the 25th O atom adsorption on HEA1 (001) surface. Yellow sphere (broken circle area) represents the 25th O atom in the sub-surface region Fig. 6 The variation in the Gibbs free energy with pressure for the adsorption of oxygen atoms on the (001) surface of HEA1 alloy for 300-1500 K from 10 −9 to 10 2 bar O 2 pressure. The distribution shows the HEA1 surface to be completely covered with a monolayer of oxygen at all pressures between 300 and 1500 K E. Osei-Agyemang and G. Balasubramanian availability of active sites. Once the exposed surface layer is completely covered with oxygen, further addition of oxygen promotes diffusion of O atoms into the second and third subsurface layers (bulk regions) of the alloy. As the full monolayer coverage is achieved, oxygen begins to move into the bulk propagating through the region rich in Zr and Ti atoms.
We characterize the thermodynamic properties of the oxidation process as a function of temperature and pressure using the Gibbs free energy of adsorption per surface area. We observe that the surface is fully covered (a single monolayer i.e., 1 ML) with oxygen at temperatures between 300 and 1500 K. At extremely high temperatures (2000 K) we observe the removal of oxygen atoms with the manipulation of the oxygen pressure. Even at 10 −9 bar, there is still 0.333 ML oxygen coverage on the (001) alloy surface. The bare alloy surface could not be completely recovered through variations of temperature and pressure.
Further efforts are required to understand the effect of different surface compositions with varying concentrations of the refractory elements to analyze the effect of such compositional changes on the oxidation mechanism. Future efforts will also involve understanding the diffusion of oxygen into the bulk alloy and measuring the kinetic rate of the associated processes. Additionally, the effect of oxidation on mechanical properties such as coefficient of thermal expansion will be an area of future investigation.

SQS models
In order to generate the random solid solutions that represent the HEA1 and HEA2 alloys, the SQS method is employed. For ordered systems, DFT can be used to generate cells with periodic boundary conditions. However, one encounters complications when treating disordered solid solutions for which the use of DFT in generating structures is inappropriate. For fivecomponent HEA1 or HEA2 alloys, one strategy to obtain a random solid solution structure is to construct large supercells where the atoms A, B, C, D, and E are used to randomly decorate the host lattices. 22 DFT calculations are however limited by the number of atoms that can be considered in the simulated system and since such large supercells for random solid solutions will involve many atoms, this method is computationally prohibitive. Zunger et al. [23][24][25] developed the SQS method to overcome the complexities associated with the large supercells as well as the meanfield theories such as coherent particle approximation (CPA) among others. SQS are built with small-unit-cell periodic structures that mimic the most relevant information such as pair and multisite correlation function of the alloy. 22 In generating the SQS structure, a distribution of distinct local environments is maintained. The average of those distinct local environments corresponds to the random solid solution. Many important properties of an alloy can be obtained from a single DFT calculation on an SQS structure based on the existence of distinct local environments. The errors introduced by representing random solid solutions with small periodically repeated unit cells can be mitigated by constructing SQS structures to reproduce the correlation functions between the first few nearest neighbors. Also, the errors introduced by periodicity are shifted to the relatively distant neighbors since the near neighbor interactions are more important.
The mcsqs code implemented in the Atom-Theoretic Automated Toolkit (ATAT) 26 is used to generate various SQS-N structures (N = 40, 80, 120, and 160) for the bcc HEA1 and HEA2 random solid solutions. Initially, all structures with N atoms per unit cell based on the bcc lattice are generated exhaustively with the mcsqs code. The pair and multi-site correlation functions for each structure are then constructed, and the structures that best fit the correlation functions of the alloys for specified sets of pair and multisite parameters are explored. The pair correlation functions are constrained to match with those of a random solid solution up to the seventh-nearest neighbor.

DFT calculations
The Vienna ab initio Simulation Package (VASP) 27 29 are used to represent the core electrons as well as the core part of the valence electron wavefunctions that are kept frozen. The PAW representation reduces the number of planewaves required to effectively describe the electrons close to the nuclei. The generalized gradient approximation (GGA) exchange correlation functional as parametrized by Perdew, Burke and Ernzerhof (PBE) 30 is used while employing the Methfessel-Paxton 31 smearing scheme by setting the gamma parameter to 0.1 eV. An energy cut-off of 600 eV is used for the planewaves expansion. A Monkhorst-Pack 32 special grid sampling of the k-points for integration of the Brillouin zone yields 4 × 4 × 4 k-points representing 36 irreducible number of sampling points for all bulk calculations and 4 × 4 × 1 k-points for surface calculations. The self-consistent field procedure is used for resolution of the Kohn-Sham equations by setting energy changes for each cycle at 10 −4 eV as the convergence criterion between two successive iterations. Both volume and shape of the unit cell as well as all the atomic positions for the generated SQS structures are thoroughly relaxed. However, the positions of all the ions in the two top most layers Fig. 7 The variation in the Gibbs free energy with O 2 pressure (from 10 −9 to 10 2 bar) for the adsorption of oxygen atoms on the (001) surface of HEA1 alloy at 2000 K. The inset shows the stability plots within the region of 10 −9 -10 −7 bar O 2 pressure. Stability at different coverages change as the pressure is varied are relaxed until the net forces acting on them are smaller than 10 −2 eV/Å for all surface calculations.
In order to analyze the phase stability of the HEA1 and HEA2 alloys, the formation enthalpies of the SQS-N structures are calculated using the DFT bulk total energies of the constituent elements in their ground state. In selecting the most appropriate surface for oxygen adsorption studies, surface energies for the (100) Due to the random nature of the HEA, cleavage of the bulk yields two inequivalent surfaces with different compositions. In order to limit the number of surfaces to be dealt with, cleavage energy is computed for the as-cleaved (001), (110), and (111) facets from and the surface with the lowest cleavage energy is used for further calculations. Here, E bulk is the energy per atom in the bulk (formula unit), E i slab is the slab energy for the i slab termination, n is the number of formula units in the bulk and A is the exposed surface area.
As different surface terminations can be obtained for either the (001), (110) or (111) surfaces, we select the appropriate surface by generating symmetric slabs with similar terminating layers (i) on both sides due to the random nature of the alloy structure. Their surface energies are computed using Eq. 2 as: Here, E slab represents the total energy of the i terminated symmetric slab, E bulk is the bulk energy per unit formula of HEA1 or HEA2, m is the total number of bulk formula units in both slabs, and A is the exposed surface area. The surface energy in Eq. 2 is for the rigid surface without any relaxations. For the relaxed surface, the relaxation energy for the slab is added to that of the rigid surface and the surface energy is calculated from Eq. 2.
(1 × 1) surface unit cells generated from the selected SQS structure are used for all surface calculations. A 15 Å vacuum is included as separation between two periodically repeated slabs to avoid surface-surface interactions. For calculations on adsorption of oxygen on the selected surface, four atomic layers are used. The topmost two layers of the surface slabs are relaxed while the remaining layers are kept fixed to mimic bulk properties.
Eleven Mo atoms, eight W atoms, two Ta atoms, one Ti atom and two Zr atoms are exposed on the selected surface as shown in Fig. 8. The (1 × 1) supercell used for oxygen adsorption on the selected surface has both a and b parameters equal to 13.1558 Å, while α and β are 90°and ϒ is 93.3723°with an exposed surface area of 172.7754 Å 2 .
Oxygen adsorbs on the surfaces in both associative and dissociative modes, and the adsorption energies are computed from Eq. 3 as: E O/Surf is the energy of the oxidized surface, E Surface is the energy of the clean surface and E O2 is the energy of the gas phase O 2 molecule. The binding energy for O 2 is predicted to be -6.04 eV which compares well with −6.08 eV calculated in earlier reports. 33 The binding energy is overestimated by 0.87 eV as compared to the experimental value of -5.17 eV 34 which is typical for GGA DFT calculations. To rectify the O 2 overbinding, the adsorption energies are corrected by adding half of the overbinding energy (~0.44 eV) to each O atom adsorbed.
The effect of different coverages of oxygen are evaluated by using the various active metal sites with increasing number of oxygen atoms or molecules at the adsorption sites after identifying the preferred mode of O 2 adsorption. As there are 24 active metal sites on the exposed surface, we evaluate 0.042 monolayer (ML), 0.083 ML, 0.125 ML, up to 1 ML full coverage. The number of available metal sites on the exposed surface defines the different coverages. We observe that spin polarized calculations have negligible effect on the O 2 adsorption so non-spin polarized calculations are used.

Atomistic thermodynamic model
We use the well-established thermodynamic model 35 to provide a relationship between the calculated parameters and experimental working conditions. Here, we assume that the adsorbed atoms or molecules on the surfaces are in thermodynamic equilibrium with the gas phase that serves as a reservoir for gas molecules. This presumption permits for a definition of the Gibbs free energy of adsorption (Δ r G) as a function of thermodynamic parameters such as temperature, pressure and chemical potential using Eq. 4: Here, E ZPE(surface/O2) is the zero point energy (ZPE) contribution of the O 2 molecule when adsorbed on the surface. ΔE o is the difference in the electronic energies of the oxidized surface and the gaseous O 2 molecule as well as any released species, while Δμ(T, p) is the difference in chemical potential of the gas phase O 2 molecules and all other released species of the reaction. ΔE o can be computed according to the surface reaction as described by Eq. 5: HEA1 Surface þ n O 2 ð Þ $ HEA1 oxidizedÀsurface þ released species (5) ΔE o is determined as [E el (oxidized-surface) + E released-specie -E el (Surface)n E el (O2) ]. This approach provides the corresponding reaction energies of the small molecules on the oxidized surface. Changes in the chemical potential are brought about by the thermal contributions of the gas phase molecules and it comprises of the temperature dependent terms Δμ°(T) and RT, as in Eq. 6.
Δμ T; P ð Þ¼Δμ o T ð Þ þ RT ln P P o Changes in chemical potential of the gas phase molecules can be calculated using statistical thermodynamics (Eq. 7) as: Thermal contributions of the small molecules upon adsorption on the surface are estimated from the changes in the vibrational, rotational and translational degrees of freedom. The values of Δμ°(T) at different temperatures are computed from standard statistical thermodynamic formulas combined with calculated frequencies from the equilibrium geometry. The total pressure P°in Eq. 6 is set at 1 bar. The distribution of the calculated Gibbs free energies (Δ r G) as a function of pressure p at different temperatures with the different coverages of O 2 molecules exemplify the thermodynamic stability of these surfaces. The Bader code 36 optimized for VASP is used to analyze the charge densities obtained for both the clean as well as the oxidized surfaces.

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