Characterizing the ZrC(111)/c-ZrO2(111) Hetero-Ceramic Interface: First Principles DFT and Atomistic Thermodynamic Modeling

The mechanical and physical properties of zirconium carbide (ZrC) are limited to its ability to deteriorate in oxidizing environments. Low refractory oxides are typically formed as layers on ZrC surfaces when exposed to the slightest concentrations of oxygen. However, this carbide has a wide range of applications in nuclear reactor lines and nozzle flaps in the aerospace industry, just to name a few. To develop mechanically strong and oxygen-resistant ZrC materials, the need for studying and characterizing the oxidized layers, with emphasis on the interfacial structure between ZrC and the oxidized phases, cannot be understated. In this paper, the ZrC(111)//c-ZrO2 (111) interface was studied by both finite temperature molecular dynamic simulation and DFT. The interfacial mechanical properties were characterized by the work of adhesion which revealed a Zr|OO|Zr|OO//ZrC(111) interface model as the most stable with an oxygen layer from ZrO2 being deposited on the ZrC(111) surface. Further structural analysis at the interface showed a crack in the first ZrO2 layer at the interfacial region. Investigations of the electronic structure using the density of state calculations and Bader charge analysis revealed the interfacial properties as local effects with no significant impacts in the bulk regions of the interface slab.


Introduction
Zirconium Carbide (ZrC), being a non-oxide ultra-high temperature ceramic, is an interesting material for several applications in severe environments. It is mostly used in environmentally harsh and demanding conditions including cutting tools, nuclear plant inner coatings, turbine components in the aerospace industry, and as refractory ceramics in the steel industry [1,2]. It fulfills all these application requirements due to its excellent mechanical and physical properties with a high melting point of 3430 • C.
However, a serious problem with the oxidation of ZrC was encountered with this material when used in harsh conditions. Indeed, ZrC forms low refractory oxides at 500-600 • C [3]. The oxide causes deterioration of the physical and mechanical properties and defeats the purpose of its general applications. Thus, there is a need to study and properly understand the oxidation process and mechanism on the low index surfaces of this nanocrystallite material. ZrC, being cubic, has three distinct low index surfaces: (100), (110), and (111) surfaces, with (100) being the most stable [4,5]. Studies on the oxidation process have been carried out on the (100) surface [6][7][8][9] and all provide similar results: the ZrC(100) Since ZrC was used as the substrate, all calculation parameters were based on the optimized values for ZrC. Density functional theory (DFT) was performed using the Vienna ab initio Simulation Package (VASP) [26] based on Mermin's finite temperature DFT [27]. The following electronic configurations were used for Zr, C, and O atoms, respectively: [Kr]4d 2 5s 2 , [He]2s 2 2p 2 , and [He]2s 2 2p 4 . We used Projected Augmented Wavefunction (PAW) [28] to describe the core electrons and the core part of the valence electrons wavefunctions and this helped in reducing the number of plane waves required to describe the electrons close to the nuclei. The Kohn-Sham valence states were expanded in a plane wave basis set with a kinetic energy cutoff of 500 eV. The generalized gradient approximation (GGA) parametrized by Perdew, Burke, and Ernzerhof (PBE) [29] was used for the exchange correlation part. Dispersion corrections have not been added to the DFT energies as, at the interfaces, covalent bonds are formed and the atomic density of the two solids are comparable. The interface interaction may be slightly overestimated but the stability order and the conclusions will not be influenced. The Methfessel-Paxton [30] smearing scheme was utilized with the gamma parameter set to 0.1 eV. For all bulk calculations, a k-sampling of a 9 × 9 × 9 mesh using the standard Monkhorst-Pack [31] special grid was employed. However, 9 × 9 × 1 k-points sampling was selected for all surface and interface slab calculations. The Kohn-Sham equations were resolved using the self-consistent field (SCF) procedure and assumed to be converged when energy changes of 1 × 10 −4 eV were obtained between two successive iterations.

Finite Temperature Molecular Dynamics (MD)
To confirm the experimental findings on the analysis of ZrC nanocrystals, finite temperature ab initio molecular dynamic simulation was performed to provide a first approximation on the nature of ZrO 2 formed on the ZrC surfaces. A (2 × 2) supercell was used for all MD simulations (13.40 × 13.40 × 30.00 Å, 60 atoms). It started with a nine-layer thick ZrC(111) substrate (terminating with four Zr atoms on both sides of the surface slab) by depositing Zr and O atoms onto the exposed ZrC(111) surface to form about two layers of ZrO 2 . The ions were initially kept at T = 100 K within the micro canonical ensemble (NVE) and the velocities were scaled upwards at different steps until a final temperature of 1000 K was reached. This temperature was selected to provide an allowance for the possible formation of the m-ZrO 2 phase which is stable at temperatures below 1450 K. A 1 fs time step was used. The resulting equilibrium structure after 20 ps was then quenched from 1000 K to 500 K. Geometries at minima on the potential energy surface were selected and optimized at a higher precision of calculation to obtain a final structure.

Bulk ZrC and c-ZrO 2 Phases
To obtain parameters optimized for the system, bulk calculations were performed on ZrC and c-ZrO 2 phases.
Energy versus volume data were obtained for both ZrC and c-ZrO 2 and finally fitted with Murnaghan's equation of state. The optimized lattice parameters were calculated from this fitting.
An optimized k-points of 5 × 5 × 5 Monkhorst-Pack grid producing 63 irreducible k-points was used for bulk ZrO 2 calculations and the same kinetic energy cut-off of 500 eV for the ZrC bulk was used for the ZrO 2 bulk. All the c-ZrO 2 bulk used contained four formula units as shown in Figure 1. the two solids are comparable. The interface interaction may be slightly overestimated but the stability order and the conclusions will not be influenced. The Methfessel-Paxton [30] smearing scheme was utilized with the gamma parameter set to 0.1 eV. For all bulk calculations, a k-sampling of a 9 × 9 × 9 mesh using the standard Monkhorst-Pack [31] special grid was employed. However, 9 × 9 × 1 k-points sampling was selected for all surface and interface slab calculations. The Kohn-Sham equations were resolved using the self-consistent field (SCF) procedure and assumed to be converged when energy changes of 1 × 10 −4 eV were obtained between two successive iterations.

Finite Temperature Molecular Dynamics (MD)
To confirm the experimental findings on the analysis of ZrC nanocrystals, finite temperature ab initio molecular dynamic simulation was performed to provide a first approximation on the nature of ZrO2 formed on the ZrC surfaces. A (2 × 2) supercell was used for all MD simulations (13.40 × 13.40 × 30.00 Å , 60 atoms). It started with a nine-layer thick ZrC(111) substrate (terminating with four Zr atoms on both sides of the surface slab) by depositing Zr and O atoms onto the exposed ZrC(111) surface to form about two layers of ZrO2. The ions were initially kept at T = 100 K within the micro canonical ensemble (NVE) and the velocities were scaled upwards at different steps until a final temperature of 1000 K was reached. This temperature was selected to provide an allowance for the possible formation of the m-ZrO2 phase which is stable at temperatures below 1450 K. A 1 fs time step was used. The resulting equilibrium structure after 20 ps was then quenched from 1000 K to 500 K. Geometries at minima on the potential energy surface were selected and optimized at a higher precision of calculation to obtain a final structure.

Bulk ZrC and c-ZrO2 Phases
To obtain parameters optimized for the system, bulk calculations were performed on ZrC and c-ZrO2 phases.
Energy versus volume data were obtained for both ZrC and c-ZrO2 and finally fitted with Murnaghan's equation of state. The optimized lattice parameters were calculated from this fitting.
An optimized k-points of 5 × 5 × 5 Monkhorst-Pack grid producing 63 irreducible kpoints was used for bulk ZrO2 calculations and the same kinetic energy cut-off of 500 eV for the ZrC bulk was used for the ZrO2 bulk. All the c-ZrO2 bulk used contained four formula units as shown in Figure 1.

Interface Model Construction
To construct the interface model, the stacking direction at the interface is initially selected and there should be a proper commensurability factor between the two bulk phases with respect to the interface plane [32]. The required surfaces are subsequently cleaved from the two bulk phases along the selected surface normal. Each of the revealed surfaces will have different atomic arrangements and configurations. The interface is finally created by bringing the two surfaces in contact with each other and then fully relaxed to obtain a final optimized interfacial geometry.

Commensurate Phases and Surface Structures
Among the low index ZrC surfaces, the (100) stoichiometric and non-polar surface is found to be the most stable [4,5]. However, even though the (111) surface is polar upon cleaving from the bulk phase by terminating on one side with a carbon layer and the other side with a zirconium layer, surface reconstruction revealed a more stable surface terminated on both sides with Zr atom layers [4]. With a lattice parameter of 6.698 Å, b of 5.801 Å, and β = 60 • , the exposed surface area of the ZrC(111) surface is 38.854 Å 2 . On accounts of several studies made on the c-ZrO 2 surfaces, the (111) surface is found to be the most stable [33,34]. Surface energies are calculated for one layer up to six layers of ZrC to ascertain the effect of the layer thickness on the surface energy. The surface energies are calculated as where E slab is the total energy of the surface slab, Ebulk is the energy per formula unit of ZrC in the corresponding bulk, A is the surface area, and n is the number of formula units in the surface slab.
Surface energies are also computed for the (111) terminations of c-ZrO 2 . The surface energies were calculated for different numbers of layers, starting from one to six layers of ZrO 2 . Upon cleaving the c-ZrO 2 along the [111] direction, a polar slab is obtained with an OO layer terminating on one side and a Zr layer terminating on the other side. The slab can, however, be terminated in three different arrangements as OO|Zr|OO|Zr|OO-, Zr|OO|Zr|OO|Zr-, and O|Zr|OO|Zr|O- (Figure 1). Only symmetric slabs were used (slabs with mirror symmetry) in the calculation of the surface energy to eliminate the net dipole moment. The calculation of the interface tension defined in a subsequent section requires these surface energies. Thus, surface energies of three different terminations were calculated: Zr-termination, O-termination, and OO-termination.
To access surface energies of both stoichiometric and non-stoichiometric slabs, the surface grand potential is defined as Ω i which implies contact of the Zr and O reservoirs with the surface: N Zr and N O are the numbers of Zr and O atoms in the slab with µ Zr and µ O being the chemical potential of Zr and O, respectively. E i slab is the total energy of the surface slab and A is the surface area. The chemical potentials of Zr and O are related by bulk ZrO 2 in the expression: µ ZrO 2 = E bulk ZrO 2 = µ Zr + 2µ O with E bulk ZrO 2 being the total energy per bulk ZrO 2 unit. Rearranging this expression and substituting it in Equation (1), we obtain the following: Defining the chemical potential of O in relation to the chemical potential of the reference state, O 2 which is defined as half the total energy of O 2 gas as and substituting it in Equation (2) with further rearrangements, we obtain: If we make the following definition: where γ i is the surface energy of the stoichiometric part of the selected slab. Substituting Equation (4) into Equation (3), we come to the following expression for the surface grand potential: Thus, the surface grand potential is defined in terms of the surface energy arising from the stoichiometric part of the slab, and another part correcting for the extra number of Zr or O atoms.
From Equation (5), a range of ∆µ O values can be accessed if we define the lower and upper limits. In defining the upper limit of the O chemical potential, we assume that the chemical potential of O must be lower than the energy of O in its reference stable gaseous state. Thus comes the upper limit of the O chemical potential: For the lower limit of the O chemical potential, if we combine the expressions µ ZrO 2 = E bulk O 2 /2 , and make rearrangements, the lower limit of the O chemical potential is obtained: A plot of the surface grand potential Ω i against the accessible range of O chemical potential is obtained for both stoichiometric and non-stoichiometric slabs for easy comparison of surface energies.
Phase commensurability is one major problem that is encountered when forming interfaces. The two surfaces used in forming the interface must be coherent due to the periodic boundary condition imposed in the calculation. The surface misfit parameter Υ can be used to obtain highly coherent interfaces [35]. This parameter is defined as: Thus, a unit cell of c-ZrO 2 with a surface area of S B is forced into coherency onto a substrate ZrC(111) with surface area S A , and the resulting overlap area between the two surfaces is S A-B . The misfit parameter measures the average length scale misfit between the two-unit cells [21] rather than an area misfit. In Table 1, the calculated misfit parameters between ZrC(111) substrate surface and all c-ZrO 2 surfaces are summarized. It is apparent from this table that the ZrC(111)||c-ZrO 2 (111) interface combination has the lowest misfit parameter of 8.2% and is acceptable. The resulting interface unit cell defined by the substrate ZrC(111) is 6.698 Å × 6.698 Å which is small and can be easily managed by the DFT calculation. The misfit parameter, being a geometrical measure, cannot be used alone in building the interface. It must be combined with other models. Two models are widely known for ensuring the commensurability of two different phases when forming an interface. Within the first approach, the unit cells of the two phases are multiplied by a factor corresponding to the other unit cell until both cells are commensurate with each other. The resulting supercell is usually large and unbearable for ab initio calculations. However, the resulting interface is coherent with very small mismatch parameters [36].
The second method is widely used [37][38][39][40][41][42] as it results in small and manageable interface supercells (a single unit cell), suitable for ab initio calculations. In this model, the lattice parameters of the phase considered as the substrate are considered for the interface, with the lattice parameter of the other phase scaled until a perfect match with the substrate lattice is obtained.

Geometrical Models for Interface
Within the slab model described for studying the interface, a thickness of 10.945 Å of ZrC (nine layers) was used. This thickness was enough to mimic the electronic structure when ionic positions in the bulk are relaxed. The c-ZrO 2 (111) units were then pinned into the registry, layer by layer on the exposed ZrC(111) surface. Thus, in straining the c-ZrO 2 to match the dimensions of the ZrC surface, coherent interfaces are ensured. The interface unit cell is therefore determined by the bulk and surface parameters of the ZrC(111). In this manner, the unit cell lattice parameter of the c-ZrO 2 (111) is shrunk by about 8%. After fixing the geometries of the two surfaces at the interface, the remaining degrees of freedom in the resulting interface structure are perpendicular to the interface and the interface chemical composition [43]. From one to five layers of the c-ZrO 2 (111) units were built on the ZrC(111) surface.
In Figure 2, we provide side views of the interface models used with the different numbers of ZrO 2 layers. Each c-ZrO 2 bilayer is approximately 3.5 Å thick. All interface models used were symmetric with respect to the center of the interface slab to remove any long-range dipole-dipole interaction between exposed surfaces. A total of 14 Å of vacuum was applied between two subsequent interface slabs to avoid any physical interactions between the slabs. The interface slab used thus has a configuration of: -ZrC(111)|c-ZrO 2 (111)|vacuum|ZrC(111)|c-ZrO 2 (111)|vacuum|ZrC(111)|c-ZrO 2 (111)|vacuum|ZrC(111)-Since the ZrC(111) slab is a reconstructed structure with four extra Zr atoms, the interface chemical composition is dependent on the number of Zr atoms (ZrC side) and the terminating layer of the c-ZrO 2 (111) phase. In constructing the interface, three different terminations along the c-ZrO 2 [111] direction were considered: Zr|OO|Zr|OO|Zr|OO-, O|Zr|OO|Zr|OO|Zr|O-, and OO|Zr|OO|Zr|OO|Zr-. We also considered the OO|Zr|OO|Zr|OO|Zr-on a ZrC(111) surface with an oxidized layer. A total of four different interface models were built as shown in Figure 2.

Mechanics and Cohesion at the Interface
In defining the interface cohesion and stability, one important parameter used is the interface tension G int , defined as the reversible work needed to separate the interface into two free surfaces [44]. According to this definition, an assumption made is that both diffusional and plastic degrees of freedom are suppressed and hence negligible. The greater the G int value, the higher the energy needed to separate the interface into two surfaces.
The interface tension can be defined according to the Dupre equation in terms of the interface and free surface energies as [45,46]:

Mechanics and Cohesion at the Interface
In defining the interface cohesion and stability, one important parameter used is the interface tension ɣint, defined as the reversible work needed to separate the interface into two free surfaces [44]. According to this definition, an assumption made is that both diffusional and plastic degrees of freedom are suppressed and hence negligible. The greater the ɣint value, the higher the energy needed to separate the interface into two surfaces.
The interface tension can be defined according to the Dupre equation in terms of the interface and free surface energies as [45,46]: || − 2 is the interface energy also known as the adiabatic work of adhesion, Wad > 0, , and − 2 are the relaxed surface energies of the ZrC(111) and c-ZrO2(111) surfaces, respectively. In this definition, the relative strength of the interface versus the bulk bonds decides the preference for the formation of either the interface or the open surfaces [22].
A measure of whether there is the formation of an interface or the free surfaces can be determined by the interface tension. The magnitude and sign of ɣint (Equation (10)) provide a measure of whether the interface bonds are stronger than the internal bonds in the separate phases [22]. The criteria are that 0 < ɣint < + − corresponds to σ ZrC||c−ZrO 2 is the interface energy also known as the adiabatic work of adhesion, W ad > 0, σ ZrC , and σ c−ZrO 2 are the relaxed surface energies of the ZrC(111) and c-ZrO 2 (111) surfaces, respectively. In this definition, the relative strength of the interface versus the

Mechanics and Cohesion at the Interface
In defining the interface cohesion and stability, one important parameter used is the interface tension ɣint, defined as the reversible work needed to separate the interface into two free surfaces [44]. According to this definition, an assumption made is that both diffusional and plastic degrees of freedom are suppressed and hence negligible. The greater the ɣint value, the higher the energy needed to separate the interface into two surfaces.
The interface tension can be defined according to the Dupre equation in terms of the interface and free surface energies as [45,46]: || − 2 is the interface energy also known as the adiabatic work of adhesion, Wad > 0, , and − 2 are the relaxed surface energies of the ZrC(111) and c-ZrO2(111) surfaces, respectively. In this definition, the relative strength of the interface versus the bulk bonds decides the preference for the formation of either the interface or the open surfaces [22].
A measure of whether there is the formation of an interface or the free surfaces can be determined by the interface tension. The magnitude and sign of ɣint (Equation (10)) provide a measure of whether the interface bonds are stronger than the internal bonds in the separate phases [22]. The criteria are that 0 < ɣint < + − 2 corresponds to weakly coupled interfaces and ɣint < 0 to strongly coupled interfaces. The calculated values of and − 2 used here are obtained from their respective relaxed bulk equilibrium phases (i.e., strain-free surface slabs).
In parallel, the adiabatic work of adhesion Wad is defined as: A measure of whether there is the formation of an interface or the free surfaces can be determined by the interface tension. The magnitude and sign of G int (Equation (10)) provide a measure of whether the interface bonds are stronger than the internal bonds in the separate phases [22]. The criteria are that 0 < G int < σ ZrC + σ c−ZrO 2 corresponds to weakly coupled interfaces and G int < 0 to strongly coupled interfaces. The calculated values of σ ZrC and σ c−ZrO 2 used here are obtained from their respective relaxed bulk equilibrium phases (i.e., strain-free surface slabs).
In parallel, the adiabatic work of adhesion W ad is defined as: is the total energy of the fully relaxed interface slab, A is the interface area, and E tot ZrC and E tot c−ZrO 2 are the total energies of the fully relaxed isolated ZrC(111) and c-ZrO 2 (111) slabs, respectively. Usually, the calculated Wad value is a lower bound as compared to values obtained in cleavage experiments due to dissipative processes in physically separating the interface [44] There is no relation between characterizing the interfacial strength and the bulk strain when depositing the c-ZrO 2 . Hence, the E tot c−ZrO 2 value used is the total energy of the strained c-ZrO 2 for commensurability with the ZrC surface. In this manner, the strain energy component between E tot c−ZrO 2 and E tot ZrC||c−ZrO 2 is canceled out since the c-ZrO 2 is in the strain state [22].
Aside from the relaxed work of adhesion, the rigid work of adhesion W rigid ad can be used in characterizing the interface cohesion and stability. In this definition, the same strained state is ensured to exist in both the interface and the free surfaces. This provides maximum cancelation for the strain energy in the calculated interface energy [43]. This quantity gives information purely on the bonds formed at the interface irrespective of the free surfaces. It is calculated by separating the optimized interface structure into the different phases and rigidly calculating their energies without allowing the phases to fully relax. Equation (11) is finally applied in calculating the rigid work of adhesion.

Interfacial Thermodynamics
For each couple of interface models used, the most stable chemical composition of the interface is determined by a thermodynamic approach. A thermodynamic grand canonical ensemble treatment is used to compare the relative stability of the models with different chemical compositions. Within such an ensemble, all models are assumed to be in chemical and thermal equilibria with bulk phases and the relevant thermodynamic quantity is the grand potential. An assumption made is that the entropic and volumetric contributions to the grand potential are negligible. For an interface of ZrC and c-ZrO 2 , the interface grand potential is defined as: where Ω i/j int is the interface grand potential, Ω ZrO 2 sur f is the surface grand potential of the exposed c-ZrO 2 side of the interface slab and Ω i/j slab , Ω ZrC , and Ω ZrO 2 are the grand potential of the interface slab, ZrC slab, and ZrO 2 slabs, respectively. N ZrC and N ZrO 2 are the number of ZrC and ZrO 2 units in the respective slabs. Substituting the following: (12), we obtain: where µ Zr , µ C , and µ O are the chemical potentials of Zr, C, and O, respectively. N k is the number of that species, E i/j slab is the total energy of the interface slab, E ZrC bulk and E ZrO 2 bulk are the bulk energies of ZrC and ZrO 2 , respectively. Upon rearrangement of Equation (13), we obtain the following: If we make the following definition, (14) with subsequent rearrangement, we obtain: If we define another quantity: Substituting Equation (16) into Equation (15), we obtain an expression for the interface grand potential as: Thus, for each overlayer termination i, the interface grand potential Ω i/j int depends on ∆µ O and ∆µ Zr . A derivation of the upper and lower boundaries of the O and Zr chemical potentials is provided here as: - For ∆µ Zr , the upper boundary is ∆µ Zr = µ Zr − E bulk Zr < 0 and lower boundary is ZrC and E f ZrO 2 are the formation energies of ZrC and c-ZrO 2 , respectively. Thus, the Zr chemical potential range is defined as: ZrO 2 /2 , and the range of chemical potentials for O is:

ZrC-ZrO 2 : From Experiments to Theoretical Studies
As already described in previous papers [23][24][25] and according to TEM-ED experiments, the presence of a different phase from ZrC bulk was noticed at the surface of the particles. An EDX analysis revealed this phase to be zirconium oxide with an estimated thickness of 5 nm. Two different orientations were observed according to high image resolution (polycrystalline oxide layer). The dhkl indexation of the diffraction patterns of the different orientations showed the presence of cubic ZrO 2 , mainly the (111) phase with some traces of the tetragonal ZrO 2 (101 and 011) phases (see Supplementary Materials, Figures S1 and S2). To more precisely investigate the interface, a theoretical point of view may therefore be necessary.

Finite Temperature Molecular Dynamic Simulation
This section discusses the results obtained during the MD simulation. A haphazard ZrO 2 structure was observed to grow on the ZrC(111) surface at 1000 K. This structure shows O atoms forming three-fold hollow bonds between three Zr atoms of the ZrC(111) surface at the interface. Upon quenching to a T = 500 K, a more ordered structure was obtained ( Figure 3). The observed pattern of the ZrO 2 atomic arrangements matches the crystal structure of c-ZrO 2 (111) and t-ZrO 2 (101) structures. Thus, the MD simulation confirms the formation of ZrO 2 on the ZrC(111) surface from atomic depositions. This further complements the experimental results.

Surface and Bulk Properties of ZrC and ZrO2
Details of the optimized lattice parameter, bulk modulus as well as the pressure derivative of the bulk modulus of ZrC, are provided in a previous paper [47]. The optimized lattice parameter for the ZrC bulk is 4.736 Å . The fitted lattice parameter for c-ZrO2 is 5.143 Å . For t-ZrO2, the calculated a parameter is 3.649 Å and c = 5.257 Å and the tetragonal distortion dz = ∆z/c = 0.050. For m-ZrO2, the calculated values are a = 5.243 Å , b = 5.307 Å , c = 5.412 Å , and β = 99.20o. All these bulk parameters are well reproduced and are in excellent agreement with both experimental and other calculated values.
In a previous paper [4], we calculated the surface energy for a ZrC(111) surface terminating with four Zr atoms on both sides of the exposed surface at the different chemical potentials of C. The surface energy was 0.169 eV Å −2 at = . In Figure 4, a stability plot is provided for the surface grand potential of each of the surface terminations of c-ZrO2(111). The surface termination O|Zr|OO|Zr|O-is observed to be the most stable. The same stable termination is found in a different theoretical work [31].

Surface and Bulk Properties of ZrC and ZrO 2
Details of the optimized lattice parameter, bulk modulus as well as the pressure derivative of the bulk modulus of ZrC, are provided in a previous paper [47]. The optimized lattice parameter for the ZrC bulk is 4.736 Å. The fitted lattice parameter for c-ZrO 2 is 5.143 Å. For t-ZrO 2 , the calculated a parameter is 3.649 Å and c = 5.257 Å and the tetragonal distortion dz = ∆z/c = 0.050. For m-ZrO 2 , the calculated values are a = 5.243 Å, b = 5.307 Å, c = 5.412 Å, and β = 99.20o. All these bulk parameters are well reproduced and are in excellent agreement with both experimental and other calculated values.
In a previous paper [4], we calculated the surface energy for a ZrC(111) surface terminating with four Zr atoms on both sides of the exposed surface at the different chemical potentials of C. The surface energy was 0.169 eV Å −2 at µ C = E bulk C . In Figure 4, a stability plot is provided for the surface grand potential of each of the surface terminations of c-ZrO 2 (111). The surface termination O|Zr|OO|Zr|O-is observed to be the most stable. The same stable termination is found in a different theoretical work [31].
For six layers of c-ZrO 2 (111) with the O|Zr|OO|Zr|O-termination, the calculated surface grand potential at µ O = E Gas O 2 /2 is 0.054 eV Å −2 which agrees very well with 0.048 eV Å −2 in different works [35]. This is not surprising as this surface termination is the only stoichiometric structure used. Focusing on the most stable surface termination, surface energies were calculated with different numbers of ZrO 2 layers. Table 2 provides a summary of the surface energies with different numbers of layers. The surface energies of the c-ZrO 2 (111) surface with O|Zr|OO|Zr|O-termination converges after three layers. Thus, the calculated surface energy of 0.054 eV Å −2 was used to compute the interface tension in Equation (10).  For six layers of c-ZrO2(111) with the O|Zr|OO|Zr|O-termination, the calculated surface grand potential at = 2 /2 is 0.054 eV Å −2 which agrees very well with 0.048 eV Å −2 in different works [35]. This is not surprising as this surface termination is the only stoichiometric structure used. Focusing on the most stable surface termination, surface energies were calculated with different numbers of ZrO2 layers. Table 2 provides a summary of the surface energies with different numbers of layers. The surface energies of the c-ZrO2(111) surface with O|Zr|OO|Zr|O-termination converges after three layers. Thus, the calculated surface energy of 0.054 eV Å −2 was used to compute the interface tension in Equation (10). Analysis of the interfacial cohesion begins with the rigid work of adhesion, previously defined in Section 2.4.3. With the rigid work of adhesion, the bulk properties of ZrC and c-ZrO2 are canceled out and the resulting parameter depends solely on interfacial properties. Table 3 provides a summary of the rigid work of adhesion values. In this table, it is apparent that the rigid work of adhesion is always lower than the relaxed work of adhesion for nearly all the interface models considered. Thus, the relaxation of the interface structure in the relaxed work of adhesion contributes significantly to the interface properties and results in a higher value than the rigid calculation. This relaxation is involved in releasing the strain imposed in the c-ZrO2 over-layer when it is forced into the registry with the ZrC substrate.

Rigid Work of Adhesion
Analysis of the interfacial cohesion begins with the rigid work of adhesion, previously defined in Section 2.4.3. With the rigid work of adhesion, the bulk properties of ZrC and c-ZrO 2 are canceled out and the resulting parameter depends solely on interfacial properties. Table 3 provides a summary of the rigid work of adhesion values. In this table, it is apparent that the rigid work of adhesion is always lower than the relaxed work of adhesion for nearly all the interface models considered. Thus, the relaxation of the interface structure in the relaxed work of adhesion contributes significantly to the interface properties and results in a higher value than the rigid calculation. This relaxation is involved in releasing the strain imposed in the c-ZrO 2 over-layer when it is forced into the registry with the ZrC substrate.  Table 3 and Figure 5 show that, using the rigid work of adhesion, the most stable interface model is the c-ZrO 2 (Zr|OO|Zr|OO-) interface model compared to the others. With this interface model, the convergence of the W rigid ad rigid parameter begins at three c-ZrO 2 layers.  Table 3 and Figure 5 show that, using the rigid work of adhesion, the most stable interface model is the c-ZrO2 (Zr|OO|Zr|OO-) interface model compared to the others. With this interface model, the convergence of the rigid parameter begins at three c-ZrO2 layers.

Relaxed Work of Adhesion
The relaxed work of adhesion is calculated for all the fully relaxed interfacial systems by allowing the separated ZrC and c-ZrO2 slabs to fully relax. This parameter characterizes the interfacial bond strengths. In Table 3, the relaxed work of adhesion calculated for the four interface models are summarized. According to Table 3, the interface model involving two layers of oxygen from the ZrO2 side (Zr|OO|Zr|OO-) is the most stable in terms of interfacial strength as observed in the calculated relaxed work of adhesion. Figure 6 provides a pictorial view of the relaxed work of adhesion for all the interface models considered. There is a convergence of the work of adhesion after three layers of ZrO2 are deposited on the ZrC. The Wad values initially decrease from one layer to two layers of ZrO2 and sharply rise at three layers of ZrO2, from which point it converges. Looking at the most stable model (Zr|OO|Zr|OO-), the Wad value sharply increases

Relaxed Work of Adhesion
The relaxed work of adhesion is calculated for all the fully relaxed interfacial systems by allowing the separated ZrC and c-ZrO 2 slabs to fully relax. This parameter characterizes the interfacial bond strengths. In Table 3, the relaxed work of adhesion calculated for the four interface models are summarized. According to Table 3, the interface model involving two layers of oxygen from the ZrO 2 side (Zr|OO|Zr|OO-) is the most stable in terms of interfacial strength as observed in the calculated relaxed work of adhesion. Figure 6 provides a pictorial view of the relaxed work of adhesion for all the interface models considered. There is a convergence of the work of adhesion after three layers of ZrO 2 are deposited on the ZrC. The Wad values initially decrease from one layer to two layers of ZrO 2 and sharply rise at three layers of ZrO 2 , from which point it converges. Looking at the most stable model (Zr|OO|Zr|OO-), the W ad value sharply increases from 0.251 eV Å −2 at two layers of ZrO 2 to 0.965 eV A −2 at three ZrO 2 layers for the on-top Zr mode of adsorption. This phenomenon is exactly the opposite of what is observed in ceramics deposited on metals where the first one and two layers of deposition are rather stronger than the deposition of three or more layers of ZrO 2 [22]. Thus, the interfacial strength depends on the first three ZrO 2 layers deposited. The same feature has been observed for metals deposited on ceramics where the metals are predicted to wet the ceramic surface but then ball up when more than one monolayer of metal is added [48][49][50]. The on-top interface models form stronger interface structures than the fcc models. Weak interfaces are formed when one and two layers of ZrO 2 are deposited but strong interfaces are obtained when three or more layers of ZrO 2 are added.
The interface tension defined in Section 2.4.3 is a good parameter for assessing interfacial mechanics and strength. It provides a measure for comparing the strength of bonds at the interface and in the corresponding bulk phases. According to the criteria defined, 0 < G int < σ ZrC + σ c−ZrO 2 is linked to a weakly coupled interface and G int < 0 to strongly coupled interfaces. With an asymptotic value of σ ZrC(111) = 0.169 eV Å −2 and σ c−ZrO 2 = 0.054 eV Å −2 combined with W ad = 0.965 eV Å −2 for the most stable interface (Zr|OO|Zr|OO-), the calculated interface tension is G int = −0.742 eV Å−2. This shows that the interfacial bonds are stronger than the internal bonds in each ceramic bulk phase. When we consider the surface energy for the c-ZrO 2 (111) with OO|Zr|OO|Zr|OO termination in an oxygen-rich environment (∆µ O = 0) as 0.292 eV Å −2 , the calculated interface tension is −0.504 eV Å −2 which is still less than zero. Moreover, G int ≤ 0 corresponds to a layer-by-layer growth of the ceramic known as the Frank-van-der-Merwe (FM) mode and the mixed-mode is also known as the Stranski-Krastanov (SK) growth mode. from 0.251 eV Å at two layers of ZrO2 to 0.965 eV A at three ZrO2 layers for the on-top Zr mode of adsorption. This phenomenon is exactly the opposite of what is observed in ceramics deposited on metals where the first one and two layers of deposition are rather stronger than the deposition of three or more layers of ZrO2 [22]. Thus, the interfacial strength depends on the first three ZrO2 layers deposited. The same feature has been observed for metals deposited on ceramics where the metals are predicted to wet the ceramic surface but then ball up when more than one monolayer of metal is added [48][49][50]. The on-top interface models form stronger interface structures than the fcc models. Weak interfaces are formed when one and two layers of ZrO2 are deposited but strong interfaces are obtained when three or more layers of ZrO2 are added. The interface tension defined in Section 2.4.3 is a good parameter for assessing interfacial mechanics and strength. It provides a measure for comparing the strength of bonds at the interface and in the corresponding bulk phases. According to the criteria defined, 0 < ɣint < is linked to a weakly coupled interface and ɣint < 0 to strongly coupled interfaces. With an asymptotic value of = 0.169 eV Å −2 and = 0.054 eV Å −2 combined with Wad = 0.965 eV Å −2 for the most stable interface (Zr|OO|Zr|OO--), the calculated interface tension is ɣint = −0.742 eV Å−2. This shows that the interfacial bonds are stronger than the internal bonds in each ceramic bulk phase. When we consider the surface energy for the c-ZrO2(111) with OO|Zr|OO|Zr|OO termination in an oxygen-rich environment (∆μO = 0) as 0.292 eV Å −2 , the calculated interface tension is −0.504 eV Å −2 which is still less than zero. Moreover, ɣint 0 corresponds to a layer-by-layer growth of the ceramic known as the Frank-van-der-Merwe (FM) mode and the mixed-mode is also known as the Stranski-Krastanov (SK) growth mode.

Relaxed Work of Adhesion
A description of the structure and properties of the most stable interface models is provided here. The relaxed stable structures for the Zr|OO|Zr|OO//ZrC(111) interface model using 1, 2, 3, 4, and 5 layers of c-ZrO2(111) are shown in Figure 7. The interface structure appears to depend somehow on the number of ceramic layers. For this stable interface model, even though the starting geometry was O atoms from the ZrO2 side adsorbed directly on top of Zr atoms of the ZrC, the final geometry was O atoms adsorbing at three-fold hollow fcc sites between three Zr (ZrC) atoms. For all layers of ZrO2, there is a rearrangement of the Zr|OO|Zr|OO-atoms upon forming the interface into a more stable O|Zr|O|Zr|OO-structure. Thus, the exposed surface of the slab terminates with an O layer.

Relaxed Work of Adhesion
A description of the structure and properties of the most stable interface models is provided here. The relaxed stable structures for the Zr|OO|Zr|OO//ZrC(111) interface model using 1, 2, 3, 4, and 5 layers of c-ZrO 2 (111) are shown in Figure 7. The interface structure appears to depend somehow on the number of ceramic layers. For this stable interface model, even though the starting geometry was O atoms from the ZrO 2 side adsorbed directly on top of Zr atoms of the ZrC, the final geometry was O atoms adsorbing at three-fold hollow fcc sites between three Zr (ZrC) atoms. For all layers of ZrO 2 , there is a rearrangement of the Zr|OO|Zr|OO-atoms upon forming the interface into a more stable O|Zr|O|Zr|OO-structure. Thus, the exposed surface of the slab terminates with an O layer. The crystal shape of the ZrC phase is maintained. However, there is a transformation of the c-ZrO2 phase into m-ZrO2. The phase transition is highly evident in the middle of the third, fourth, and fifth layers of the ZrO2 interface structures with 3-fold and 4-fold O atoms as well as 6-fold and 7-fold Zr atoms. The (c→m) transformation is, however, not a complete one due to the restraint imposed on periodic boundary conditions. At T = 0 K, the m-phase is about 7% larger in volume than the c-phase and hence this transformation reduces the misfit of about 8% already calculated for this interface model. A similar transformation pattern is found in ceramics deposited on metals [22].
Two types of oxygen bonds are observed at the interface: O1 atoms (ZrO2) closer to the interface plane, bonding at fcc sites on the ZrC surface, and O2 atoms (ZrO2) above the interface plane, bonding directly on top of Zr (ZrC) atoms. The fcc bonds are exactly the same found in a previous study when the ZrC(111) surface was completely oxidized with The crystal shape of the ZrC phase is maintained. However, there is a transformation of the c-ZrO 2 phase into m-ZrO 2 . The phase transition is highly evident in the middle of the third, fourth, and fifth layers of the ZrO 2 interface structures with 3-fold and 4-fold O atoms as well as 6-fold and 7-fold Zr atoms. The (c→m) transformation is, however, not a complete one due to the restraint imposed on periodic boundary conditions. At T = 0 K, the m-phase is about 7% larger in volume than the c-phase and hence this transformation reduces the misfit of about 8% already calculated for this interface model. A similar transformation pattern is found in ceramics deposited on metals [22].
Two types of oxygen bonds are observed at the interface: O1 atoms (ZrO 2 ) closer to the interface plane, bonding at fcc sites on the ZrC surface, and O2 atoms (ZrO 2 ) above the interface plane, bonding directly on top of Zr (ZrC) atoms. The fcc bonds are exactly the same found in a previous study when the ZrC(111) surface was completely oxidized with a monolayer of oxygen [10]. It was also noticed to passivate the ZrC(111) surface with no further diffusion of oxygen into the bulk. In the 1-layer ZrO 2 interface, no fcc bonds of the O1 atoms were observed due to the O2 atoms pushing outwards from their bulk positions. These O2 atoms do not show any bonding with the Zr (ZrC) atoms in all layers of ZrO 2 deposited. The bond distances of the O2 type atoms with the Zr (ZrC) atoms have a minimum of 4 Å. Consequently, at more than three ZrO 2 layers, there is a structural failure of ZrO 2 coating on the ZrC substrate. This crack is highly visible in the 4-ZrO 2 layer interface structure in Figure 7. The crack leaves a monolayer of oxygen deposited on the ZrC. Thus, there is a mixed-mode of deposition of ZrO 2 on ZrC. A monolayer of oxygen is formed on the ZrC surface and the remaining ZrO 2 layers ball up. This is not surprising as the very negative interface tension calculated in Section 3.4.2 suggests a mixed-mode of deposition. Additionally, the high work of adhesion, showing an over-adhered interface, results in failure in the ceramic layer, a few angström distances from the interface plane. Using

Thermodynamic Stability of the Interface
In this section, an analysis of the thermodynamic stability of the different interface models used is considered. The stability of the interface is calculated with respect to the bulk phases and not the surface slabs forming the interface. The interface grand potential Ω i/j int provides a measure of stability for the different interface models in different terminations. According to Equation (17), the interface grand potential has an interface-dependent term ∅ i/j int which only differs from the corresponding surface term γ i in Equation (4). As such, for the purpose of brevity and clarity, Table 4 provides the ∅ i/j int values. These terms do not depend on the chemical potentials of excess Zr and O species. In Table 4, it is apparent that the interface is readily formed for all models (due to the negative values of the grand interface dependent terms) with the OO|Zr|OO|Zr-ZrC(111) model being the least stable amongst them. It is, however, evident from Table 4 that the most stable model is the Zr|OO|Zr|OO//ZrC(111) fcc model, contradicting the top Zr|OO|Zr|OO//ZrC(111) model calculated to be the most stable using the relaxed work of adhesion parameters. However, it is worth pointing out that even though the fcc and top models used different starting points of interfacial atom adsorption, the two models resulted in the same configuration with interface O atoms bonding at fcc sites. In addition, the difference between the calculated relaxed work of adhesion values for the two models is only 0.042 eV Å −2 using four layers of c-ZrO 2 . The interface grand potential thus corroborates the stability criteria established in Section 3.4.2, with Zr|OO|Zr|OO//ZrC(111) being the most stable. When the chemical potentials of excess O and Zr atoms are considered, the calculated Ω i/j int value, when minimized for every (∆µ Zr , ∆µ O ) pair, resulted in Zr|OO|Zr|OO//ZrC(111) being the most stable model formed. The electronic features at the interface were analyzed by first considering the density of states when the free surfaces form the interface.
As an initial description of the density of states (DOS) at the interface, a total DOS (TDOS) was obtained by projecting the density of states onto all atoms at the interface and compared with the DOS of the individual separate surfaces. In Figure 8, the TDOS for the most stable interface model, Zr|OO|Zr|OO//ZrC(111) using three layers of ZrO 2 is shown. Figure 8  To understand the shift in bands when the surface atoms come together to form the interface, the DOS are projected onto the atoms at the interface and the corresponding atoms in the surface slabs in Figure 9. Upon forming the interface, the Zr (ZrC) conduction bands shift to higher energies as they are filled with electrons from O (ZrO2) with a dip at the Fermi level, further stabilizing the interface formed compared to the high level of states at the Fermi level of the corresponding surface. The Zr (ZrC) atoms form a new sharp core state at −18 eV corresponding to the Zr-O bonds formed at the interface. Thus, the Zr (ZrC)-O (ZrO2) bond closest to the interface plane is highly localized. There are no significant changes in the C (ZrC) states upon forming the interface with only the valence bands shifting slightly to higher energies. For Zr (ZrO2) bands, they are shifted to To understand the shift in bands when the surface atoms come together to form the interface, the DOS are projected onto the atoms at the interface and the corresponding atoms in the surface slabs in Figure 9. Upon forming the interface, the Zr (ZrC) conduction bands shift to higher energies as they are filled with electrons from O (ZrO2) with a dip at the Fermi level, further stabilizing the interface formed compared to the high level of states at the Fermi level of the corresponding surface. The Zr (ZrC) atoms form a new sharp core state at −18 eV corresponding to the Zr-O bonds formed at the interface. Thus, the Zr (ZrC)-O (ZrO2) bond closest to the interface plane is highly localized. There are no significant changes in the C (ZrC) states upon forming the interface with only the valence bands shifting slightly to higher energies. For Zr (ZrO2) bands, they are shifted to lower energies upon forming the interface. The sharp Zr (ZrO2) band at −15.6 eV which was initially mixed O (ZrO2) closest to the interface plane is now lost. This further explains the breakage of the c-ZrO2 (111) phase after depositing an oxygen layer on the ZrC(111) surface. These Zr (ZrO2) bands become diffuse when forming the interface as their bonds with the closest O (ZrO2) to the interface plane are broken and the electrons are delocalized. The high surface states of the O (ZrO2) closest to the interface plane are quenched and their sharp peaks at −15.6 eV are also lost as they break off from the c-ZrO2 phase and deposit on the ZrC(111) phase. In the valence band, the O (ZrO2) closest to the interface plane interacts strongly with the Zr (ZrC) bands due to their sharp peaks while the O (ZrO2) far from the interface plane interacts weakly with its diffuse bands and delocalized electrons.  Figure 10 provides a spatial profile of the electronic structure upon moving from the stable bulk faces to form the exposed surfaces and subsequently the interface. Thus, in Figure 10, the DOS are projected onto the atoms parallel to the interface plane, moving from bulk ZrC to the interface region, then bulk c-ZrO 2, and finally to the exposed c-ZrO 2 (111) surface. The C bands in bulk ZrC remain fairly the same upon forming the interface. In addition, the Zr valence band of ZrC are shifted to lower energies when they form the interface with the introduction of a sharp core state at −18.1 eV which mixes with the O core states of ZrO 2 closest to the interface plane. At the interface region, the valence band is a mixture of Zr (ZrC) d-bands and O (ZrO 2 ) p-bands with a covalent bonding nature. There is a clearly diffuse interaction of the O (ZrO 2 ) p-bands far from the interface plane, further explaining the weak interaction of the second O layer (ZrO 2 ) far from the interface plane when the interface is formed. It is also apparent that the O (ZrO 2 ) p-bands far from the interface maintain the bulk nature of such bands with a similar feature for the Zr (ZrO 2 ) d-bands far from the interface plane. The conduction band in this bulk ZrO 2 region is mainly Zr-d states. Moving to the exposed surface on the ZrO 2 side, the Zr-d states in the conduction band are reduced, shifted to lower energies, and new states appear close to the Fermi level while the valence band becomes less diffuse. The exposed ZrO 2 surface is stabilized by the oxygen termination as there are relatively no states around the Fermi level. This termination is the same foundation for the most stable exposed surface of c-ZrO 2 (111).
plane, further explaining the weak interaction of the second O layer (ZrO2) far from the interface plane when the interface is formed. It is also apparent that the O (ZrO2) p-bands far from the interface maintain the bulk nature of such bands with a similar feature for the Zr (ZrO2) d-bands far from the interface plane. The conduction band in this bulk ZrO2 region is mainly Zr-d states. Moving to the exposed surface on the ZrO2 side, the Zr-d states in the conduction band are reduced, shifted to lower energies, and new states appear close to the Fermi level while the valence band becomes less diffuse. The exposed ZrO2 surface is stabilized by the oxygen termination as there are relatively no states around the Fermi level. This termination is the same foundation for the most stable exposed surface of c-ZrO2(111). Figure 10. PDOS of atoms at the interface, in the bulk, and the exposed surface of the four c-ZrO2 layer interface slab (Zr|OO|Zr|OO-ZrC (111)). The states are aligned parallel to the interface plane with each region marked on the top layer. The Fermi level is aligned at the zero energy position.

Charge Transfer Analysis
A Bader charge analysis is provided in this section to understand how charges are transferred when forming the interface from the cleaved surfaces. The charge transfer on the atoms is arranged in a spatial profile of the atoms from the surface to respective layers in the different phases forming the interface as described by Christensen and Carter [21]. As a first step, we provide a charge transfer analysis for the most stable interface with different layers of c-ZrO2(111) starting from one to five in a spatial profile and finally analyzing the four c-ZrO2 (111) Figure 10. PDOS of atoms at the interface, in the bulk, and the exposed surface of the four c-ZrO 2 layer interface slab (Zr|OO|Zr|OO-ZrC (111)). The states are aligned parallel to the interface plane with each region marked on the top layer. The Fermi level is aligned at the zero energy position.

Charge Transfer Analysis
A Bader charge analysis is provided in this section to understand how charges are transferred when forming the interface from the cleaved surfaces. The charge transfer on the atoms is arranged in a spatial profile of the atoms from the surface to respective layers in the different phases forming the interface as described by Christensen and Carter [21]. As a first step, we provide a charge transfer analysis for the most stable interface with different layers of c-ZrO 2 (111) starting from one to five in a spatial profile and finally analyzing the four c-ZrO 2 (111) layer stable interface model, Zr|OO|Zr|OO//ZrC(111) top site. Table 5 shows that upon forming the interface from the free surfaces, there is a significant amount of charge transfer in the interfacial region. ) closest to the interface from the remaining part of the ZrO 2 phase. In all the layers of ZrO 2 deposited, there is essentially no significant charge redistribution in the atomic layers farther away from the interface moving toward the bulk region of the respective phases. With three or more ZrO 2 layers deposited, the high charge transferred from interface cations to anions releases part of the axial strain imposed on the ZrO 2 phase as it is brought into registry with the ZrC substrate upon forming the interface. The magnitude of the charges transferred from Zr (ZrC) to O (ZrO 2 ) for the different layers of ZrO 2 is in the same pattern as the relaxed work of adhesion computed for the different numbers of layers. In Table 6, a spatial profile of the charge transfer is provided to obtain a further understanding of the electronic structure of the interface formed and how it affects the corresponding bulk and surface properties. The most stable Zr|OO|Zr|OO-ZrC(111) top site model with four layers of ZrO 2 is used in this analysis. In this table, the average charge distribution per layer is arranged in a profile moving from the interface plane towards the bulk and then the exposed surface regions. In this analysis, it is essential to compare the ∆Q s values of Table 6 when analyzing the surface regions (with respect to atoms in the corresponding surface slabs), while ∆Q b values are used when making bulk comparisons (with respect to corresponding atoms in the bulk phases). In the interfacial region, there is a significant amount of charge transfer, compared to the corresponding free surface. This high amount of charge distribution mainly originates from Zr (ZrC) and O atoms (ZrO 2 ) at the interface. However, moving away from the interface into the bulk regions of both phases, and finally, to the exposed surface area of ZrO 2 , there is virtually no charge redistribution. The ∆Q b values for the bulk regions also show no charge redistribution when forming the interface. This further ascertains the interface to be affected only by atomic layers closer to the interface plane. Table 6. The spatial profile of charge transfer is arranged along the normal direction of the interface. ∆Q s is the difference in the charge of the ion with the corresponding ion in the isolated surface slab used to create the interface, and ∆Q b is the difference in charge between the ion and the corresponding ion in the bulk structure.

Summary and Conclusions
Different theoretical approaches have been used to characterize the oxide layer formed on ZrC nanoparticles. Finite temperature ab initio molecular dynamic was used to build two layers of ZrO 2 on a ZrC(111) surface. An ordered layer of ZrO 2 was observed to form on the ZrC(111) surface. Then, periodic DFT was used to characterize the interface between ZrC(111) and c-ZrO 2 (111). The preferred interface observed was formed between Zr terminated ZrC(111) and an OO-terminated c-ZrO 2 (111) leading to a final interface model Zr|OO|Zr|OO//ZrC(111). The main mechanical property used to characterize the interface was the relaxed work of adhesion W ad . This value reveals the Zr|OO|Zr|OO//ZrC(111) model as the most stable interface with a high W ad value compared to the other models. A thermodynamic analysis investigating the interface grand potential also confirmed the Zr|OO|Zr|OO//ZrC(111) model as the most stable interface. A close examination of the structural properties revealed the deposition of an oxygen layer from the ZrO 2 phase onto the ZrC(111) surface with the remaining part of ZrO 2 breaking off from the interface, suggesting a crack at this interface. This also was in very good agreement with a previous study on the oxidation of ZrC(111) surfaces which revealed the formation of an oxygen layer on the surface with the further addition of O atoms leading to an endothermic reaction. In addition, the electronic structure of the interface was analyzed using the density of states calculated for the interface models. The DOS showed no major changes induced by the formed interface but only features of the deposited oxygen layer. In using the Bader charge analysis, an enormous amount of charge transfer at the interface was discovered, originating mainly from cations on the ZrC side (Zr) to the O atoms from ZrO 2 at the interfacial region. Moreover, there was virtually no charge redistribution in the bulk regions of the interface slab. Thus, the interfacial properties were governed by local effects, only confined to the first two atomic layers around the interface plane. To conclude, this study brought new insights into the structure and composition of the surface/interface of ZrC. This will pave the way for new opportunities in the field of hybrid materials with high-performance properties in harsh environments.