Molecular dynamics of liquid Zirconia: Effect of Pressure

: This paper carefully examined the influence of pressure within the range of (0 – 20) GPa on the self- diffusion coefficients of zirconium and oxygen in cubic zirconia by a molecular dynamics (MD) simulations at a temperature of 3000 K using 96 -1500 particles. Data obtained showed that, at 96 particles the self-diffusion coefficients decreased from 1.60 x 10 -4 (Å2/ps) to 0.76 x 10 -4 (Å2/ps) between 0.0 GPa and 20.0 GPa. As the pressure is increased, the self-diffusion coefficient for zirconium and oxygen atoms decreases. Particularly, in the region of high pressure, the decrease in self- diffusion coefficient for oxygen is nonlinear. This might be as a result of instability associated with large systems where one oxygen atom moves over the maximum potential of another oxygen atom. This is consistent with previous MD studies for oxygen in liquid SiO 2 . The pair correlation function was found to increase with pressure. These thermo physical properties find practical applications in the electro ceramic industry.

Zirconia (ZrO2) has been found to be useful in the manufacture of electro ceramics, ceramic tiles, ceramic sensors, medical devices, refractory materials and paint pigments. This has attracted the attention of both experimental and theoretical investigations for the past few decades. In recent times, a lot of work has been carried out in the production of ZrO2 (Hrovat et al., 2003;Blumenthal, 1998;Agrawal and Sudhakar, 2002;Hannink et al., 2000;Reynen et al., 1981;Dodd et al., 2001;Paola et al., 2007).Other important features of zirconia include the fact that at low pressure, it crystallizes in three phases: monoclinic (P21/c) below 1400K,tetragonal (P42/nmc) between 1400 -2570K and cubic (Fm3m) between 2570 -2980K structures (Oliveira et al., 2001;Murty et al., 1994). It is also a very hard and stable material.
It is important to note that the dynamics in liquid materials is governed by atomic diffusion which in turn describes the single-particle transport property in molten state. To have better understanding of this property, series of studies have been carried out. For example, the diffusion property of zirconia was studied by Joshua Pelleg(2016). Here, the selfdiffusion of oxygen in zirconia melt was reported. Also, self-diffusion coefficients of zirconium were measured between 1773K and 1373K using radioactive Zr 95 as a tracer (Kidson and McGurn, 1961). Suarez et al. (2008) performed a sintering kinetic process on cubic zirconium and reported the cation diffusion of Zr between 1473K and 1673K. The oxygen-18 gas-solid exchange technique was also used to study the diffusion of oxygen in monoclinic zirconia, and self-diffusion of oxygen was calculated (Keneshea and Douglass, 1971). However, to the best of our knowledge, there is no information on the study of pressure dependence of self-diffusion coefficients for zirconium and oxygen in liquid zirconia. Therefore, in this study we used the molecular dynamics simulations to study the influence of pressure within the range of (0 -20)GPa on the selfdiffusion coefficients of zirconium and oxygen in cubic zirconia.

MATERIALS AND METHODS
The face-centered cubic crystal structure of zirconia (ZrO2) was used for our simulations. We started by first converting the primitive cell of bulk zirconia (consisting of one zirconium (Zr) atom and two oxygen (O) atoms) to a unit cell of four Zr atoms and eight O atoms. This conventional unit cell is repeated twice along the A, B and C axes to create a 2X2X2 cubic supercell of bulk ZrO2. The resulting supercell now contains 96 particles (32 Zr atoms and 64 oxygen atoms). Similarly, the unit cell is repeated three, four and five times along A, B and C axes to respectively generate 3X3X3, 4X4X4, and 5X5X5 cubic supercells consisting of 324, 768 and 1500 particles.

NENUWE, ON; AZI, OJ
We investigated the dynamics of liquid cubic ZrO2 for systems consisting of 96, 324, 768 and 1500 particles at 3000K, and at different pressure of 0, 5, 10, 15 and 20 GPa. The calculations are done using the molecular dynamics (MD) technique as implemented in the ATK-VNL (Atomistix Toolkit, 2017;Griebelet al., 2007;Griebel and Hamaekers 2004). In our MD method the interactions between the atoms are described by the Wang_2012 potential function (Wang et al., 2012) as shown below: Where V(rij) is the potential energy between atoms i and j separated by rij, C, A and ρ are adjustable parameters, qiand qj charges at a separation of rij, and 0 ε is permittivity of free space. The MD process is divided into three steps: (i) heating the solid ZrO2 above its melting point (2988K), (ii) annealing the resulting liquid ZrO2 at 3000K to bring the system to equilibrium, (iii) further annealing the equilibrated liquid ZrO2 to obtain more physical data for accurate calculations of the atomic diffusion coefficient, by using the mean-square displacement extracted from the respective MD trajectory.
To achieve the first step, we use the atk-forcefield code together with the Wang_2012 (Wang et al., 2012) potential for the simulation. The first MD object is added to heat the crystalline ZrO2 from 1500K to 3000K, using an NPT Martyna Tobias Klein (NPT-MTK) barostat (Martyna, et al., 1994). The NPT barostat is used to equilibrate the ZrO2 supercell under constant temperature and pressure. The NPT-MTK was allowed to run for 200000 steps at log intervals of 2000 and Maxwell-Boltzmann type initial velocity was used at 1500K. The NPT reservoir temperature and final temperature is respectively taken as 1500K and 3000K. This allows us to increase the temperature of the supercell gradually above the melting point of ZrO2. It implies that the ordered solid ZrO2 will eventually become a liquid, which usually exhibits a short-range order as illustration by the radial distribution function g(r). Since our system has cubic symmetry, we apply isotropic pressure of 0 GPa for the simulation. This enables all cell vectors to be changed by the same factor, which preserves the shape of the cell during the simulation.
For the second step of the MD simulation, this is annealing liquid ZrO2 to bring it to equilibrium. We set the MD object to equilibrate the supercell structure of liquid ZrO2 at constant temperature 3000K, using the same NPT barostat. This helps us to eliminate any memory effect of the initial solid zirconia structure on the physical properties of liquid zirconia obtained from the initial structure. Here, we also used 200000 steps at log intervals of 2000 for the NPT-MTK. At this stage the Maxwell-Boltzmann temperature, reservoir and final temperatures are all set to 3000K. The isotropic pressure is also taken as 0GPa.
Finally, we used the NPT-MTK to collect sufficient physical data required for calculation of self-diffusion coefficient of zirconium and oxygen in liquid zirconia. The same number of steps and log intervals are used. This time the initial velocity is taken as the configuration velocity. This allows the system to use the initial velocities from the annealing stage (second MD) for the third step of the MD simulation. The reservoir and final temperature are also taken as 3000K. Then, the system is allowed to run the entire simulation. When the calculation is done, we repeated the entire process using isotropic pressure of 5, 10, 15 and 20 GPa. For each pressure, the calculation is allowed to run successfully and the results are analyzed by extracting the mean-square displacement (MSD). The slope of the MSD gives the self-diffusion coefficients as discussed by Nenuwe and Osafile (2018) in their calculation of self-diffusion coefficients of liquid aluminum at different temperatures.

RESULT AND DISCUSSION
Curves representing the mean-square displacement with time as a function of pressure and the number of atoms for different systems are displayed in Figures 1  and 2 for zirconium and oxygen atoms, respectively. Figures 1 and 2 show that in the region of high pressure the MSD versus simulation time graphs for all systems converge to a well-defined slope throughout the simulation time. The convergence of the results with increasing number of particles in Figures 1 and 2 suggests that even 96 particles are sufficient to determine the diffusion property of ZrO2 under high pressure. We observed that both at low and high pressure the diffusion behaviour of the system do not strongly depend on the number of particles used. This is in agreement with James et al. (1990).  The self-diffusion coefficients for zirconium and oxygen atoms in liquid ZrO2 were calculated by taking the slopes of each curve (using a Python script) at different pressures, and the results are displayed in Table 1 and Figure 3. These results strongly suggest that as we increase the pressure of the system, the amplitude of self-diffusion maximum for zirconium decreases, and the position of the diffusivity is shifted to higher pressure for all system sizes. This trend of pressure dependence has been reported earlier by Rudd et al. (2012). For oxygen, the diffusivity maximum decreased and later increased and decreased again with increasing pressure. This might be as a result of an instability associated with large systems where one oxygen atom moves over the potential maximum of another oxygen atom. This anomalous behaviour for oxygen was also reported by James et al., (1990) in their molecular dynamics study of .
The radial distribution function (RDF) or pair correlation function g(r) being an effective tool for describing the structure of disordered molecular systems was used to determine the convergence of the simulation at low and high pressure as the size of the system is increased from 96 to 1500 particles. In solid, g(r) is known to have infinite number of pointed peaks whose heights and separations are characteristic of the lattice structure. While in the liquids phase, the pair correlation function has small number of peaks at short distances that steadily decays to unity at longer distances r(Å) (Andrew, 2001). We demonstrated in Figures 4 and 5 that at both 0 and 20 GPa the pair correlation function g(r) versus position curves for all systems converged very well. This is an indication that the simulations reached an equilibrium state and the diffusion coefficients are well defined. Also, as the position r(Å) increases, g(r) decays to a constant value.

Fig 3. Self-diffusion coefficients as a function of pressure for Zirconium and Oxygen
It is clear that, at short distances less than atomic diameter, g(r) is zero. This is as a result of strong repulsive forces. In Figure 4, with the 96 particles system, at zero pressure (0GPa), the first and largest peak for O-O occurs at 2.625 Å, with g(r) having a value of 2.348. This is an indication that it is two times more likely that two oxygen atoms would be found at this separation. Then, the pair correlation function curve falls rapidly and passes through a minimum value around 3.125 Å. The probability of finding two oxygen atoms within this separation is less. Similarly, the first peak for Zr-O occurs at 2.125 Å, with the pair correlation g(r) having a value of 4.534. This is an indication that it is four times more likely that one zirconium atom and one oxygen atom would be found at this position. Then, g(r) falls to a minimum value of 0.0378 around 3.375 Å. In this case the chances of finding zirconium-oxygen atoms around this separation are almost zero. Also, the first and largest peak for Zr-Zr curve occurs around 3.675 Å, with the pair correlation g(r) having a value of 5.583. This suggests that it is five times more likely to find two zirconium atoms at this position. Beyond this point the radial distribution function falls rapidly and passes through a minimum value (0.0034) around 4.425 Å.
Here, the probability of finding two zirconium atoms within this separation goes to zero. At long distances, g(r) approaches 1 which implies the absence of long-range order. It is important to state that absence of long-range order is an indication that solid ZrO2 has completely melted.For the system with 1500 particles, at zero pressure (0 GPa), the first and largest peak for the O-O curve occurs at 2.575 Å, with g(r) having a value of 2.3156. This is an indication that it is two times more likely that two oxygen atoms would be found at this separation. Then, the RDF falls quickly and passes through a minimum value of 0.4863 around 3.075 Å. The probability of finding two oxygen atoms within this separation is less than the probability at 2.575 Å. For the Zr-O curve, the first peak occurs at 2.125 Å, with the correlation function g(r) having a value of 4.8019. This is an indication that it is four times more likely that one zirconium atom and one oxygen atom would be found at this separation. Then, g(r) falls to a minimum value of 0.0161 around 3.325 Å. In this case the probability of finding zirconiumoxygen atoms within this separation goes to zero. For the Zr-Zr curve, the first peak occurs at 3.675 Å, with g(r) having a value of 5.7619. This is an indication that it is five times more likely that two zirconium atoms would be found at this separation. Then, g(r) falls to a minimum value of 0.0017 around 4.425 Å. In this case the probability of finding two zirconium atoms within this separation goes to zero. At a pressure of 20GPa, in the 96 particles system, the first and largest peak for the O-O curve occurs at 2.525 Å, with g(r) having a value of 2.553. This is an indication that it is two times more likely that two oxygen atoms would be found around this point. Then, the RDF falls quickly and passes through a minimum value of 0.2833 around 3.025 Å. Therefore, the probability of finding two oxygen atoms within this separation is less. For the Zr-O curve, the first peak occurs at 2.125 Å, with g(r) having a value of 5.3577. This is an indication that it is five times more likely that one zirconium atom and one oxygen atom would be found at this separation. Then, g(r) falls to a minimum value of 0.0045 around 3.175 Å. In this case the probability of finding zirconium-oxygen atoms with this separation goes to zero. For the Zr-Zr curve, the first peak occurs at 3.575 Å, with g(r) having a value of 5.8494. This is an indication that it is five times more likely that two zirconium atoms would be found at this separation. Then, g(r) falls to a minimum value of 0.000 around 4.325 Å. In this case the probability of finding zirconium-oxygen atoms within this separation is zero. The effect of pressure on the radial distribution function is also examined for all systems simulated. The pressure dependence of pair correlation function g(r) is displayed in Figure 6. It is clear that the pair correlation function is sensitive to pressure. The first peak in the pair correlation function g(r) at each simulated pressure increases with pressure and this is shown in Figure 6. This signifies increase in the number of Zr-Zr, Zr-O and O-O bonds. Particularly, for the 768 and 1500 particles systems, it was observed that the radial distribution function curves for Zr-Zr show a non-linear increase. The might be associated with the size of the system. Conclusion: This study has described the effect of pressure in liquid zirconia at 3000K as the pressure is increased from 0.0 GPa to 20 GPa. Results obtained suggest that self-diffusion coefficients decrease with increasing pressure. Also, the pair correlation function is observed to be sensitive to pressure. The behavior of self-diffusivity and pair correlation function of liquid materials in response to effects of pressure has practical applications in the electro ceramic industry.