Unveiling oxidation mechanism of bulk ZrS2

Transition metal dichalcogenides have shown great potential for next-generation electronic and optoelectronic devices. However, native oxidation remains a major issue in achieving their long-term stability, especially for Zr-containing materials such as ZrS2. Here, we develop a first principles-informed reactive forcefield for Zr/O/S to study oxidation dynamics of ZrS2. Simulation results reveal anisotropic oxidation rates between (210) and (001) surfaces. The oxidation rate is highly dependent on the initial adsorption of oxygen molecules on the surface. Simulation results also provide reaction mechanism for native oxide formation with atomistic details.


Introduction
Transition metal dichalcogenides (TMDCs) are layered materials with promising electronic and optoelectronic applications such as field-effect transistors (FETs). Among TMDCs, ZrS 2 exhibits superior electrical properties [1]. However, ZrS 2 is known to oxidize under ambient conditions [2]. Formation of the native oxide in TMDCs and their properties dictate device processing and their applicability. For example, oxidation of TMDCs results in a reduction of on-state current in FET [3] and changes in work function [4]. While the long-term stability of TMDCs has been studied in different environments to slow down their degradation [5], oxidation mechanisms of ZrS 2 remain less understood.
In this work, we performed first principles-informed reactive molecular dynamics (RMD) simulations [6] to study atomistic oxidation processes in ZrS 2 . We first developed new reactive forcefield (ReaxFF) parameters for Zr/O/S using multi-objective genetic algorithm (MOGA) [7][8][9]. The optimized forcefield is able to reproduce quantummechanically computed charge values and bond-population dynamics. Using this forcefield, we performed RMD simulations and compared the oxidation processes at (210) and (001) surfaces. We observed a higher oxidation rate at (210) surface compared to (001) surface. On both surfaces, the oxidation process was found to be diffusion controlled.

Methods
ReaxFF parameters were optimized in two steps. In the first stage, we optimized the atomic charges of Zr, O and S atoms. The ground truth for charge optimization consisted of Mulliken charges in 53 small clusters, which were computed using Q-Chem program [10].
To refine the forcefield to better reproduce the groundtruth Zr-S and Zr-O bond-population dynamics in small quantum molecular dynamics (QMD) simulation of ZrS 2 oxidation, RMD simulations were performed with the same schedule for the QMD. QMD simulation was performed using the Vienna Ab initio Simulation Package (VASP) [11], while RMD simulation was performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) package [12]. We simulated a Zr 48 S 96 slab immersed in 72 O 2 molecules in the canonical (or NVT) ensemble at 1,200 K for 11,155 fs. The numbers of Zr-S, Zr-O, S-S and S-O bonds were calculated as a function of time. The root mean-square deviations in the numbers of these bonds between QMD and RMD simulations were tabulated as a 4-column list, for MOGA training of the ReaxFF parameters.
In the second stage of optimization, we trained Zr-S and Zr-O 2-body forcefield terms to correctly describe the bond breaking and bond formation during the oxidation reaction of ZrS 2 . For MOGA training, we used the NSGAII algorithm [13] implemented in our EZFF forcefield training software [14]. The mutation probability was set to be 0.2, and the population size was 1,024. The evolution process consisted of four steps in each iteration ( Fig. 1): (1) select parents which are fit for reproduction; (2) perform crossover and mutation on the selected parents to generate feasible solutions, i.e., offspring; (3) combine the original parents and the offspring, select the next population using fast non-dominated sorting; (4) truncate the population size to maintain a constant population size. Out of the parent and offspring generations, we maintain the same population size to form the next parent population. Selection is based on rank, and if individuals with same rank are encountered, crowding distance is compared. The criterion is to choose first by lower rank and then by higher crowding distance. The stopping rule here was a maximum number of 200 iterations or no improvement in fitness value for some fixed number of generations.

Results
To assess the accuracy of the optimized forcefield on reproducing atomic charges, we compared the ReaxFF charges with the density functional theory (DFT) charges computed by Q-Chem. Figure 2 shows small differences between the two calculations, which confirms that the charge values obtained by optimized forcefield are consistent with DFT charges.
In MOGA training to improve bond-population dynamics, we considered three genetic-algorithm variants: NSGAIII, epsNSGAII and NSGAII. Figure 3 compares time evolution of the numbers of Zr-S and Zr-O bonds between QMD (blue) and RMD (red) simulations, both before and after the corresponding genetic-algorithm refinement of ReaxFF parameters. After the parameter tuning, Fig. 3b-d, f-h exhibits much better agreement between QMD and RMD simulations compared with those before the training in Fig. 3a, e. Thus, RMD simulations with the refined ReaxFF parameters provide quantitatively correct bond-population dynamics consistent with the ground-truth QMD. We decided to use NSGAII because it showed best accuracy. After 130 iterations, we stopped the training process because error did not decrease further. From Fig. 3, we see that the forcefield trained by NSGAII algorithm is accurate enough to be used in RMD simulations to study the oxidation of ZrS 2 .
Using the refined ReaxFF parameters, we performed RMD simulations in order to compare the oxidation processes of ZrS 2 slabs on (210) and (001) surfaces. Periodic boundary conditions have been applied to both systems. To accelerate the oxidation process within accessible simulation time, the oxygen partial pressure was set to be 1.26 kbar and reaction temperature to be 1200 K. Figure 4a shows the initial simulation setups for both simulations. Each system consisted of a six-layer ZrS 2 slab with (210) or (001) surface, which was immersed in O 2 atmosphere. Figure 4b shows snapshots of both (210) and (001) simulations at different times, t = 0.02, 0.2, 0.3 and 0.4 ns. At 0.02 ns, the oxide grows much faster in the right slab surface in (210) simulation, where zirconium atoms are directly exposed to O 2 atmosphere. The anisotropic rate of oxidation is thus highly dependent on the initial adsorption. Detailed analysis reveals that the subsequent oxidation pathway in both directions involves adsorption of oxygen on ZrS 2 surface, followed by amorphization and oxygen transport into bulk, leading to a continuous oxidation. Detailed oxidation mechanisms are explained elsewhere [15].
As a closer look into the oxidation mechanism, Fig. 5 shows time evolution of the oxide depth at different surfaces. The general growth trend is similar for the right and left (210) surfaces, which oxidize much faster than both top and bottom (001) surfaces. As discussed before, the zirconium atoms are directly exposed to O 2 atmosphere on the right (210) surface, which are easier to be oxidized. According to Li et al. [16], the adsorption energy of O 2 on (210) edges are generally more negative than the pristine (001) surface for TMDCs, though they considered Zr edge with S coverage of Since Zr edge coordinated with less sulfur atoms have more negative adsorption energy and hence more energy is released by oxidation reaction, (210) surface is oxidized much faster, where the zirconium atoms on the right (210) surface are coordinated with 3 sulfur atoms. On (001) surface, the oxide grows very slow before 0.2 ns. A sudden increase in the oxide depth at 0.2 ns in the top (001) surface is caused by the formation of a disordered ZrS 2 region between ZrS 2 layers, and extensive oxidation assisted by the O transport therein, as shown in Fig. 4b. Both top and bottom (001) surfaces show a similar oxidation trend later from 0.3 to 0.4 ns. While a similar fast oxide growth period is observed for (001) and (210) surfaces, the oxidation rate is slightly faster on (210) surface in that period. This is likely a geometric effect. As can be seen in Fig. 4a, spacing between ZrS 2 layers provide channels for oxygen to diffuse in the (210) case, which facilitates faster oxidation.

Summary
In summary, we have optimized reactive forcefield parameters for Zr/O/S using MOGA. Using the optimized forcefield, we have performed RMD simulations to study oxidation of ZrS 2 . We found faster oxidation dynamics on (210) surface than on (001) surface. The calculated oxide growth dynamics indicates a diffusion-controlled oxidation mechanism on both surfaces.