University of Southern California Dissertations and Theses

Reactive and quantum molecular dynamics study of materials: from oxidation and amorphization to interfacial charge transfer

(USC Thesis Other)

Reactive and quantum molecular dynamics study of materials: from oxidation and amorphization to interfacial charge transfer

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content Reactive and Quantum Molecular Dynamics Study of Materials: From Oxidation and Amorphization to Interfacial Charge Transfer by Liqiu Yang A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY CHEMICAL ENGINEERING December 2023 Copyright 2023 Liqiu Yang ii Acknowledgements I would like to thank Prof. Priya Vashishta, Prof. Aiichiro Nakano, and Prof. Rajiv Kalia from Collaboratory for Advanced Computing and Simulations (CACS). I am grateful to Prof. Priya Vashishta, who has been guiding me through my PhD, and it would be hard for me to finish without his wisdom and advice. I am grateful to Prof. Rajiv Kalia, who also taught me a lot in his class. I would like to thank Prof. Paulo Branicio, who is a good professor and helped me cultivate my research skills. I would like to thank Patricia Wong, who has been giving me lots of advice on life and guiding me get used to the abroad environment. I would like to thank Dr. Ken-ichi Nomura, who also gives lots of good advice during my PhD. I would also like to thank Wei Wu for serving as my committee member. I am grateful to Sungwook Hong, Aravind Krishnamoorthy, Subodh Tiwari, Lindsay Bassman, Kuang Liu, Ankit Mishra, Beibei Wang, Pankaj Rajak. Their guidance helped me a lot to start my PhD journey and helped me gain research skills. I would also like to thank other CACS group members -- Tom Linker, Taufeq Mohammed Razakh, Ruru Ma, Nitish Baradwaj, Jingxin Zhang, and Tian Sang. I would like to thank Prof. Rafael Jaramillo and Seong Soon Jo, for their valuable experimental collaborations. I would also like to thank Prof. Fuyuki Shimojo, Shogo Fukushima, and Kohei Shimamura, who help me to use QXMD software. I also appreciate discussions with Antonina L Nazarova and Hikaru Ibayashi. I am also thankful to Marco Olguin, who provides lots of technical support for using USC’s Center for Advanced Research Computing (CARC). I acknowledge Viterbi school of engineering/graduate school fellowship and teaching/research assistantship from University of Southern California, and TLARGI fellowship. I am thankful to iii CARC at USC, Theta and Polaris from Argonne, where I have been performing most of my simulations and calculations. I would also like to thank fundings from U.S. Department of Energy that supported part of my research. I am thankful to my friends, Suyue Yuan, Jiaoyang Li, Fan Ke, Kewei Wang, and other friends. I always remember the happiness and support they brought me. Last, I must thank my family for their unconditional love and support and encouragement through my PhD life, which give me energy to overcome the barriers when I encounter problems. iv Table of Contents Acknowledgements .........................................................................................................................ii List of Tables .................................................................................................................................vi List of Figures................................................................................................................................vii Abstract .................................................................................................................................xi Chapter 1: Introduction............................................................................................................. 1 1.1 ZrS2 Oxidation ................................................................................................................ 1 1.2 Photoexcitation-induced GeTe Local Disorder............................................................... 2 1.3 Surface Transfer Doping of MoO3-x on Hydrogenated Diamond ................................... 3 Chapter 2: Methods................................................................................................................... 5 2.1 Molecular Dynamics (MD) Simulation .......................................................................... 5 2.2 Interatomic Force Field................................................................................................... 7 2.3 Periodic Boundary Conditions (PBC)............................................................................11 2.4 Reactive Molecular Dynamics (RMD)..........................................................................11 2.5 Density Functional Theory (DFT) ................................................................................ 12 2.6 Quantum Molecular Dynamics (QMD)........................................................................ 13 2.7 Nonadiabatic Quantum Molecular Dynamics (NAQMD)............................................ 16 2.8 Visualization ................................................................................................................. 16 2.9 Machine Learning and High Performance Computing................................................. 16 Chapter 3: ZrS2 Oxidation ...................................................................................................... 18 3.1 Introduction................................................................................................................... 18 3.2 Method .......................................................................................................................... 21 3.3 Force Field Optimization .............................................................................................. 21 3.3.1 Training ReaxFF Parameters .................................................................................... 21 3.3.2 Workflow for Force Field Training........................................................................... 23 3.3.3 Force Field Accuracy ................................................................................................ 24 3.4 Oxidation Mechanism................................................................................................... 27 3.4.1 Comparison of Thermal Oxidation between MoS2 and ZrS2.................................... 27 3.4.2 Comparison of Oxidation on ZrS2 (210) and (001) Surfaces ................................... 31 3.4.3 DFT Calculations of the Gibbs Free Energy Change of Adsorption ........................ 34 3.4.4 Pressure-Controlled Layer-by-Layer to Continuous Oxidation................................ 39 v 3.4.5 Critical Oxidation Moments ..................................................................................... 56 3.5 Summary ....................................................................................................................... 63 Chapter 4: Photoexcitation-Induced GeTe Local Disorder..................................................... 66 4.1 Introduction................................................................................................................... 66 4.2 Method .......................................................................................................................... 68 4.3 Results........................................................................................................................... 68 4.3.1 Structural Changes Mechanism ................................................................................ 69 4.3.2 Electronic Structure Calculations ............................................................................. 74 4.4 Summary ....................................................................................................................... 85 Chapter 5: Surface Transfer Doping of MoO3-x on Hydrogenated Diamond ......................... 87 5.1 Introduction................................................................................................................... 87 5.2 Method .......................................................................................................................... 88 5.3 Results........................................................................................................................... 89 5.3.1 Interfacial Structure .................................................................................................. 90 5.3.2 Charge Transfer......................................................................................................... 91 5.4 Summary ....................................................................................................................... 97 Chapter 6: Conclusions........................................................................................................... 98 References ............................................................................................................................. 101 vi List of Tables Table 3. 1. Chemical potential of oxygen ..................................................................................... 37 Table 3. 2. Pseudo-code for cell-based surface construction method ........................................... 50 Table 5. 1. Charge transfer per elemental atoms due to deposition .............................................. 92 vii List of Figures Figure 2. 1. Dimensionless Lennard-Jones potential. ..................................................................... 9 Figure 2. 2. Energy terms included in ReaxFF where Ebond is the bond energy, Elp is the lone pair energy, Eover and Eunder are the over- and under- coordination energies, Eval is the valence angle energy, Epen is the penalty energy, Ecoa is the three-body conjugation energy, Etors is the torsion angle energy, Econj is the four-body conjugation energy, Ehbond is the hydrogen bond interaction energy, EvdWaals, ECoulomb and Elg are the van der Waals, Coulomb energies and dispersion energy, respectively. .................................................................................................... 10 Figure 3. 1. Training set for optimizing charges. .......................................................................... 22 Figure 3. 2. Workflow for training bonds population. .................................................................. 24 Figure 3. 3. Comparison of QM charges (black color) and ReaxFF charges (red color) for selected clusters. Green, yellow and red spheres represent Zr, S and O atoms, respectively. ...... 25 Figure 3. 4.Reduction of bond error as a function of the epoch in multi-objective genetic algorithm training.......................................................................................................................... 26 Figure 3. 5.Comparison of QMD (blue) and RMD (red) simulation results for time evolution of Zr-S (a, b) and Zr-O (c, d) bonds, both with initial (a, c) and optimized (b, d) RMD force fields.............................................................................................................................................. 26 Figure 3. 6. Comparison for the performance of force fields optimized by different algorithms. Time evolution of Zr-S and Zr-O bonds are plotted for QMD (blue) as well as for RMD (red) with force fields using (a) NSGAIII, (b) epsNSGAII, and (c) NSGAII algorithms..................... 27 Figure 3. 7. RMD simulations of MoS2 and ZrS2 oxidation. Snapshots of (a-b) MoS2 oxidation and (c-d) ZrS2 oxidation at time t = 0 and 2.5 ns, respectively. Spheres represent individual atoms: Blue (Mo), green (Zr), yellow (S) and red (O). Time evolution of corresponding bond count for (e) MoS2 and (f) ZrS2 oxidation. ................................................................................... 29 Figure 3. 8. Atomistic mechanisms of ZrS2 oxidation. Cross section (a) and top-down (b) views at t = 2 ps, showing O2 adsorption and Zr-O bond formation that involves both 2- and 3- fold coordinated oxygen atoms (enclosed by black rectangles). Cross section (c) and close-up (d) view at t = 2.5 ns, showing the formation of an amorphous matrix and the closing of the vdW gaps; S atoms are not shown in (d) for clarity. (e) Oxygen transport mechanism of Zr-O bond switching, including oxygen switching between 2- and 3-fold coordination: (i-ii) Breakage of Zr1-O1 bond and decrease in coordination number (3→2); (ii-iii) formation of Zr4-O1 bond and increase in coordination number (2→3); (iii-iv) breakage of Zr2-O1 bond and reduced coordination number (3→2); (iv-v) formation of Zr5-O1 bond and coordination number recovery (2→3)................................................................................................................. 31 Figure 3. 9. (a) Initial system setups of ZrS2 slabs with (210) and (001) surfaces. Green, yellow and red spheres represent Zr, S and O atoms, respectively. (b) Snapshots of oxidization process in both (210) and (001) simulations at times, t = 0.02, 0.2, 0.3 and 0.4 ns...................... 33 Figure 3. 10. Oxide depth as a function of time on different surfaces: The right and left (210) surfaces (black) as well as top and bottom (001) surfaces (red)................................................... 33 Figure 3. 11. Calculated binding energies of oxygen on Zr sites of ZrS2 (001) surface for various coverages.......................................................................................................................... 36 Figure 3. 12. Pressure dependence of the change of Gibbs free energy upon oxygen-atom adsorption. Calculated change of Gibbs free energy, ∆, upon oxygen adsorption on ZrS2(001) viii surface with varying oxygen chemical potential ∆μO. The corresponding pressure at various temperatures is shown as bars below the plot............................................................................... 38 Figure 3. 13. Initial and final configuration of RMD simulations of ZrS2 oxidation. Initial configuration in the parallel front view (a) and perspective side view (b). Final configurations of ZrS2 oxidation at time 0.6 ns under pressure p = 0.1 GPa (c), 0.08 GPa (d) and 0.06 GPa (e). Green, yellow and red spheres represent Zr, S and O atoms, respectively............................. 40 Figure 3. 14. Snapshots of RMD simulations of ZrS2 oxidation at different pressures and times at temperature 300 K. (a-c) ZrS2 oxidation after time t = 0.2 ns at pressures 0.1 GPa (a), 0.08 GPa (b) and 0.06 GPa (c). (d-f) ZrS2 oxidation after time t = 0.6 ns at pressures 0.1 GPa (d), 0.08 GPa (e) and 0.06 GPa (f). Oxygen atoms are colored yellow. Zr atoms are colored according to the number of coordinated oxygen atoms, i.e., Zr atoms with no coordinated oxygen atoms are colored red, while the color of Zr atoms coordinated to more oxygen atoms is graded toward blue. To better visualize the oxide growth, surfaces of the oxidation front are shown in cyan color, and S atoms are omitted for clarity. ............................. 41 Figure 3. 15. Pressure dependence of ZrS2 oxidation dynamics. Time evolution of the oxide depth (a), number of Zr-O bonds (b), and number of Zr-S bonds (c). Blue, red and black curves represent oxidation at pressures 0.1, 0.08 and 0.06 GPa, respectively.............................. 42 Figure 3. 16. Time evolution of Zr-O pair distribution function. Blue, red and black curves are for pressures 0.1, 0.08 and 0.06 GPa, respectively, at times t = 0.1 ns (a), 0.2 ns (b), 0.3 ns (c), 0.4 ns (d), 0.5 ns (e), and 0.6 ns (f).................................................................................... 44 Figure 3. 17. Snapshot of interface structure at time 0.6 ns under pressure 0.1 GPa. The color scheme is the same as that in Figure 3.1 in the main text. Zr-O bonds are drawn to highlight the progress of oxidation. .............................................................................................. 46 Figure 3. 18. Zr-O bond-length distribution at 0.6 ns under pressures (a) 0.1 GPa, (b) 0.08 GPa, and (c) 0.06 GPa................................................................................................................... 47 Figure 3. 19. O-Zr-O bond-angle distribution at 0.6 ns under pressures (a) 0.1 GPa, (b) 0.08 GPa, and (c) 0.06 GPa................................................................................................................... 47 Figure 3. 20. Snapshots of ZrS2 oxidation under 0.1 GPa at times t = 0.3 ns (a), 0.4 ns (b), and 0.5 ns (c). The color scheme is the same as that in Figure 3.14............................................. 47 Figure 3. 21. Mean square displacement of oxygen atoms in the amorphous oxide after the onset of continuous oxidation at time 0.16 ns under pressure 0.1 GPa. ....................................... 48 Figure 3. 22. Oxygen concentration profile along the slab-normal direction, z, at time 0.47 ns under pressure 0.1 GPa. ................................................................................................................ 49 Figure 3. 23. Fitting of the oxide growth depth vs. time for pressures (a) 0.1 GPa, (b) 0.08 GPa and (c) 0.06 GPa. Simulation results are shown blue, red and black, respectively, for 0.1, 0.08 and 0.06 GPa, while fits are shown as gray dashed lines. All fittings are with 95% confidence bound. The three stages of oxidation (I, IIa and IIb) are indicated by magenta, green and cyan backgrounds, respectively.................................................................................... 52 Figure 3. 24. Snapshots of RMD simulation trajectory for ZrS2 oxidation under 0.1 GPa. Two examples in (a) and (b) show the transport of O atoms (labeled O1 in purple color) assisted by the Zr atoms (labeled Zr1 in blue color) bonded to inner layer. .............................................. 55 Figure 3. 25. Snapshots of O transport through gap closing for ZrS2 oxidation under 0.1 GPa in (a) NVT and (b) NPT ensembles. ............................................................................................. 56 Figure 3. 26. Substitution of a S atom by an O atom. (a) Snapshots of an oxygen atom replacing a sulfur atom. (b) Schematic of the substitution. .......................................................... 58 Figure 3. 27. Movement of a zirconium atom toward the inner layer. ......................................... 59 ix Figure 3. 28. Oxygen transfer in the ZrS2 inner layer through bond switching and structure rearrangement. .............................................................................................................................. 61 Figure 3. 29. Charge redistribution and bonding characteristics during O1-S1-Zr4-S2-Zr1 ring rearrangement. Mulliken charges (a) and COHP (b) analysis for t = 2.341 ps. Mulliken charges (c) and COHP (d) analysis for t = 2.551 ps...................................................................... 62 Figure 3. 30. Oxygen transfer in the ZrS2 inner layer through Zr-O bond rotation...................... 63 Figure 4. 1. Changes in the structure of GeTe with photoexcitation. (a) GeTe structure after thermalization at 300 K. (b-e) GeTe structures after 5 ps of photoexcitation at different excitation levels n: n = 2.6% (b), 4.0% (c), 5.2% (d) and 7.5% (e). Green and yellow spheres represent Ge and Te atoms, respectively. ...................................................................................... 69 Figure 4. 2. Change in Ge-Te bond length distribution after thermalization at 300 K. (a) GeTe bond length distribution in the perfect cubic structure at 0 K. (b) Ge-Te bond length distribution in the GeTe crystal thermalized at 300 K. . ............................................................... 70 Figure 4. 3. Photoexcitation-induced structural change in GeTe. (a) Mean square displacement of Ge as a function of time; (b) Evolution of the fraction of wrong bonds. The black curve corresponds to adiabatic MD simulation, n = 0.0%. Cyan, red, blue, and green curves correspond to excitation of n = 2.6, 4.0, 5.2, and 7.5%, respectively............................... 71 Figure 4. 4. Local disorder triggered by Ge diffusion at n = 4.0%. (a-c) Shift of Ge atom from octahedral to tetrahedral sites for t = 1.833 ps (a), 2.094 ps (b), and 3.942 ps (c). (d-f) Schematic of the corresponding local disorder mechanism.......................................................... 73 Figure 4. 5. GeTe density of states (DOS) at different simulation times for photoexcitation level n = 4.0%. (a) Total DOS at different times after photoexcitation. The black, red, blue, and green curves represent GeTe before photoexcitation (t = -0.1 ps) and after photoexcitation at t = 0.01, 0.1, and 2.5 ps, respectively. (b) Partial DOS of cubic GeTe. (c) Partial DOS of the relaxed structure. (d-f) Partial DOS at excitation level n = 4.0% after t = 0.01 ps (d), 1 ps (e), and 2.5 ps (f). The blue, orange, yellow, purple, green, and cyan curves in Figures 4.5(b-f) represent Ge s, Ge p, Ge d, Te s, Te p, and Te d, respectively. ..................................................... 74 Figure 4. 6. Atom and orbital projected band structure of cubic GeTe (Fm3m). (a) Ge s, (b) Ge p, (c) Te s, (d) Te p. The colors indicate the weight of the projections. The Fermi level energy is set at 0 eV. (e) The Brillouin zone paths used. ............................................................. 75 Figure 4. 7. Average Mulliken charge as a function of time for different photoexcitation levels. (a-d) photoexcitation for n = 2.6% (a), 4.0% (b), 5.2% (c) and 7.5% (d). The blue and red curves represent the average Mulliken charge for Ge and Te atoms, respectively. . .............. 76 Figure 4. 8. Mulliken bond overlap as a function of time for Ge-Ge, Ge-Te, and Te-Te bonds. The black curve corresponds to time t = -0.1 ps (before photoexcitation). Red, blue, and green curves represent results at t = 0.01, 1.0, and 2.5 ps (after photoexcitation)........................ 77 Figure 4. 9. Crystal orbital Hamiltonian population bonding analysis of Ge-Te and Ge-Ge interaction. (a-c) Bond length (a), COHP (b) and (c) before Ge1-Ge2 bond formation for t = - 0.001 ps. (d-f) Bond length (d), COHP (e and f) after Ge-Ge bond formation for t = 2.122 ps. (g-i) Bond length (g), COHP (h and i) after Ge-Ge bond formation for t = 3.942 ps................... 79 Figure 4. 10. COHP bonding analysis of Ge1-Te5 and Ge1-Te6 interaction. (a-c) COHP bonding analysis for t = -0.001 ps (a), 1.833 ps (b) and 2.094 ps (c). The red and blue curves in (a-c) represent the interactions of Ge1-Te5 and Ge1-Te6, respectively.................................... 80 Figure 4. 11. ELF of Ge-Te and Ge-Ge interaction. (a) ELF before Ge1-Ge2 bond formation for t = -0.001 ps. (b) ELF after Ge-Ge bond formation for t = 2.122 ps. (c,d) ELF after Ge-Ge x bond formation for t = 3.942 ps. The scale is from 0 (blue) to 1 (red). ELF = 1 represents perfect covalent bonds and any values between 0.5 and 1 reveal covalent bonds of various bonding strength, although ELF = 0.5 gives a metallic system. ................................................... 81 Figure 4. 12. Evolution of Kohn-Sham energy eigenvalues and the excited orbitals for photoexcitation level n = 4.0%. (a) Kohn-Sham energy eigenvalues as a function of time. The red, yellow, and blue curves represent orbitals that are empty, singly occupied, and doubly occupied, respectively. (b-d) Excited orbitals for t = 1.833 ps (b), 2.094 ps (c) and 3.942 ps (d).................................................................................................................................................. 82 Figure 4. 13. Charge density difference (CDD) of different GeTe structures. (a, d) Bond lengths and CDD analysis for cubic GeTe; (b, e) for distorted rock-salt (rhombohedral) GeTe; (c, f) for structure after relaxation at 300 K. The Ge-Te bonds in (a-c) are colored according to their bond length, the CDD maps in (d-f) show the difference in the electron charge distribution as a result of bonding in the different structures. ...................................................... 83 Figure 5. 1. MoO3-x encapsulated H-diamond systems with x being 0.7 (a), 0.4 (b), and 0.1 (c). The grey, red, blue, green balls represent C, H, Mo, and O atoms, respectively.................... 90 Figure 5. 2. Mo-O pair distribution for MoO3-x encapsulated H-diamond systems with x being 0.7 (a), 0.4 (b), and 0.1 (c). ................................................................................................. 90 Figure 5. 3. Mo-O bond length distribution for MoO3-x encapsulated H-diamond systems with x being 0.7 (a), 0.4 (b), and 0.1 (c). ...................................................................................... 91 Figure 5. 4. Average Bader charges of atoms before and after deposition. Dimond-MoO3-x systems with x being 0.7 (a), 0.4 (b), and 0.1 (c). The grey, red, blue, green balls represent C, H, Mo, and O atoms, respectively. ........................................................................................... 93 Figure 5. 5. Charge density profiles for MoO2.9 deposited hydrogenated diamond. (a) Side view of charge density difference of MoO2.9 deposited hydrogenated diamond. Yellow and brown regions represent electron and hole accumulation (electron depletion), respectively. (b) Side view of iso-surface of charge density at the interface with iso-level of -0.08 signifies. the hole distribution, shown in dark brown. The color scheme for atoms is the same as that in Figure 5.1. ..................................................................................................................................... 94 Figure 5. 6. Density of states for pristine hydrogenated diamond and MoO3-x deposited hydrogenated diamond systems with x being 0.7, 0.4, and 0.1. EF represents the Fermi level. Gray dashed and dash-dotted lines indicate the valence-band top and Fermi energies, respectively. .................................................................................................................................. 95 Figure 5. 7. (a) Density of states for pristine hydrogenated diamond and molten MoO2.9. (b) Density of states for MoO2.9 deposited hydrogenated diamond.................................................... 96 xi Abstract The end of Moore’s law poses enormous scientific and technological challenges. First-principles quantum molecular dynamics (QMD) and reactive molecular dynamics (RMD) simulations could revolutionize future semiconductor technology by providing fundamental mechanism at atomic and electronic scale underlining the synthesis and performance of these future semiconductor devices. This thesis contains results from three projects: (1) Controlled oxidation of 2D transition metal dichalcogenide (TMD) to synthesize semiconductor-oxide heterostructure; (2) nonthermal switching of phase change material (PCM); and (3) performance of MoO3 for hole doping of hydrogenated diamond. We investigate controlled oxidation of 2D transition metal dichalcogenide to synthesize semiconductor-oxide heterostructure: While transition metal dichalcogenides (TMDCs) have been reported as potential candidates for next-generation electronic and optoelectronic devices, their native oxidation remains a major issue in achieving their long-term stability. Among TMDCs, ZrS2 shows superior electrical properties. However, ZrS2 has been reported to easily oxidize under ambient conditions. The oxidation processes deserve widespread attention. A detailed atomistic understanding into the oxidation mechanism is essential for designing better semiconductor devices. To accomplish this goal, we first optimized the reactive force field ReaxFF parameters for Zr/O/S using multi-objective genetic algorithm (MOGA) implemented in our EZFF force field training software, so that the new force field is able to give consistent charge values and bond-population dynamics with quantum-mechanically computed information. The detailed steps for force field xii optimization are provided along with the accuracy of the optimal force field. Later, we performed RMD simulations employing our optimized force field to study the oxidation processes of ZrS2. In our RMD simulations, (210) and (001) surfaces of ZrS2 are exposed to O2 atmosphere separately at a simulation temperature of 1200K. The oxidation rate at (210) surface was determined to be higher than that at (001) surface. The oxidation rate is highly dependent on the initial adsorption of oxygen molecules on the surface. Oxidation rate on ZrS2 (001) surface is found to be accelerated under higher pressure at a simulation temperature of 300 K. In order to investigate the pressure effect, the Deal-Grove model is discussed for its applicability in ZrS2 oxidation system, and the initial oxidation stage is examined using density functional theory. Our simulation also found layer-by-layer oxidation to amorphous-oxide-mediated continuous oxidation. Visualization of oxide surface plots are shown to elucidate this process. The atomistic understanding lays the foundation for describing oxidation of ZrS2. It is hoped that this understanding will be valuable for actual device processing. Second, the nonthermal switching of phase change materials for ultrafast memory devices is studied. The contrasting electrical and optical properties of the crystalline and amorphous phases of GeTe, a primary phase-change material in the (GeTe)x(Sb2Te3)1-x alloy family, are of great importance for non-volatile storage devices and advanced computing architectures. Here, by using nonadiabatic quantum molecular dynamics to investigate the evolution of GeTe excited states, we reveal a photoexcitation-induced picosecond nonthermal path for the loss of long-range order in GeTe. A valence electron excitation threshold of 4% is found to trigger local disorder by switching Ge atoms from octahedral to tetrahedral sites and promoting Ge-Ge bonding. The resulting loss of long-range order for a higher valence electron excitation fraction is achieved without fulfilling the Lindemann criterion for melting, therefore utilizing a nonthermal path. The photoexcitation- xiii induced structural disorder is accompanied by charge transfer from Te to Ge, Ge-Te bonding-toantibonding, and Ge-Ge antibonding-to-bonding change, triggering Ge-Te bond breaking and promoting the formation of Ge-Ge wrong bonds. These results provide an electronic-structure basis to understand the photoexcitation-induced ultrafast changes in the structure and properties of GeTe and other phase change materials. In the third part, we investigate the performance of MoO3 for hole doping of hydrogenated diamond for energy efficient high-power electronics. Hydrogen-terminated diamond field-effect transistors are promising candidates for next-generation high-power electronics. While hydrogenated diamond surface deposited with MoO3 was shown to exhibit high p-type surface conductivity, atomistic and electronic structures of the interface remain elusive, including the likely effects of oxygen vacancies. Here, we develop and optimize a reactive force field for the CH-Mo-O system, which describes charge transfer and bond formation processes with quantummechanical accuracy. Reactive molecular dynamics simulations are performed using this force field to elucidate the deposition mechanism of MoO3-x on hydrogenated diamond (111) surface. Electronic density of states alignment and charge transfer at the interface are studied using firstprinciples calculations based on density functional theory for selected thermalized structures taken from the reactive molecular dynamics trajectories. We observe shifting of the electronic density of states alignment and higher charge transfer for higher Mo oxidation state. Such atomistic and electronic details provide mechanistic understanding of the surface transfer doping process. 1 Chapter 1: Introduction In this thesis, I will use reactive and quantum molecular dynamics simulation mainly together with density functional theory calculations to investigate the behaviors of materials under different synthesis conditions. 1.1 ZrS2 Oxidation Following the seminal discovery of two-dimensional (2D) material, graphene—a single-atom thick, semi-metallic layer of carbon, transition metal dichalcogenides have gained significant interest. Graphene exhibits high carrier mobility and excellent molecular barrier properties (1, 2). However, pristine graphene is not suitable for field-effect transistors (FETs) (3) due to its zero band gap. Although band gap of graphene is tunable, bandgap engineering such as nanostructuring (4-6) and chemical functionalization (7) not only diminishes carrier mobility but also adds complexity to the design. In contrast, TMDCs with sizable band gaps around 1-2 eV (2, 8) are considered promising for FETs and optoelectronic devices. TMDCs are a class of materials with components of XY2, where X is a transition metal element from group IV, V and VI, and Y is a chalcogen element atom. Among TMDCs, ZrS2 is an indirect semiconductor material (9) and especially its 1T phase shows superior electrical properties for various applications (10, 11). Furthermore, Ultra-thin ZrS2 has been reported as potential material for photocatalysis and lithium-ion batteries (12), whereas monolayer ZrS2 has been predicted theoretically to be efficient water-splitting catalyst (13, 14). However, ZrS2 is known to oxidize under ambient conditions (15). Formation of the native oxide in TMDCs and their properties dictate device processing and their applicability. For example, 2 oxidation of TMDCs results in a reduction of on-state current in FET (16) and changes in work function (17). The long-term stability of TMDCs has been studied in different environments to slow down their degradation (18). Generally, most widely studied members of TMDCs such as MoS2, MoSe2, WS2 and WSe2 are air-stable, based on the oxidation study of their bulk crystals(19, 20). On the other hand, surfaces of bulk crystals HfSe2, MoTe2 and WTe2 are air-sensitive for rapid oxidation rate under ambient conditions (21-23). Monolayer MoS2 and WS2 are reported to undergo dramatic oxidation under ambient conditions (24). Qiang et al. systematically investigated the energetics of atomic-scale oxidation and degradation mechanisms of group VIB using density functional theory (DFT) calculations (25). However, oxidation mechanisms of ZrS2 remain less understood. Understanding oxidation mechanisms of layered semiconducting transition metal dichalcogenide is important not only for controlling native oxide formation but also for synthesis of oxide and oxysulfide products. 1.2 Photoexcitation-induced GeTe Local Disorder Phase-change materials (PCMs) exhibit remarkable differences in the resistivity of their crystalline and amorphous phases, which can be used to construct effective non-volatile memory devices (26). GeTe is a prototype binary PCM in the most promising (GeTe)x(Sb2Te3)1-x (GST) alloy family, capable of switching rapidly between its vacancy-free crystalline structure and its amorphous phase (27). GeTe is a narrow band-gap semiconductor and ferroelectric with the simplest conceivable structure containing just two atoms in a unit cell (28). Intense efforts have been made to understand the nature of the phase-change process, in particular changes in the bonding and local atomic arrangement. PCMs display unique characteristic features that make them ideal for memory applications, such as unusually high optical and electrical 3 contrast between their crystalline and amorphous phases, stability of both crystalline and amorphous phases, and fast switching between these states (29). The contrast in properties is reported to be resulted from resonance bonding in the crystalline phase (30-32), when elements with an average five valence electrons form six -bonds with the nearest neighbors using three porbitals; that is, bonds on each side of participating atoms are formed by sharing the same p-orbital. Long-range order is essential for resonance bonding to exist what leads to a drastic change in bonding when the material is amorphized. The resonance bonding in IV-VI materials is formed through the use of the backside lobes of the same p-orbitals that are used for the covalent amorphous structure. Since the p-orbitals have pronounced anisotropy, the bonds formed though their use in both crystalline and amorphous phases are strongly directional. Electronic structure calculations are instrumental in the description of bonding and the study of bonding changes in solids. Many electronic structure-based tools are available for this purpose including the crystal orbital Hamilton population (COHP), the electron localization function (ELF), and the charge density difference (CDD). Based on electronic structure calculations it was found that in amorphous GeTe, Ge atoms are sp3 hybridized while Te atoms feature nonbonding lonepair electrons. It was suggested that the formation and presence of Ge-Ge bonds stabilizes the amorphous state (33). It was suggested that relocating electrons from bonding to antibonding states by electronic excitation could potentially soften and destabilize the lattice structure at temperatures well below the melting temperature (34). Despite these theoretical developments, the dynamic sequence of atomistic events that lead to phase change and their time scales remain unknown. 1.3 Surface Transfer Doping of MoO3-x on Hydrogenated Diamond 4 Diamond is the next generation wide bandgap (5.5 eV) electronic device material. However, compared to other semiconductor materials, it is hard to achieve traditional substitutional doping in diamond lattice. Surface transfer doping (STD) has been reported to be potential solution. MoO3 is an effective electron accepting material, which improves the performance and stability of STD. MoO3 is an effective electron accepting material, which improves the performance and stability of STD. While hydrogenated diamond surface deposited with MoO3 was shown to exhibit high ptype surface conductivity, atomistic and electronic structures of the interface remain elusive, including the likely effects of oxygen vacancies. 5 Chapter 2: Methods In this thesis, I use classical molecular dynamics with reactive force field ReaxFF, ab initio molecular dynamics, and DFT to investigate oxidation dynamics in TMDCs. I use NAQMD and DFT to investigate laser induced amorphization in GeTe. MoO3-x/hydrogen/diamond structure is also investigated using reactive molecular dynamics, ab initio molecular dynamics, and DFT. 2.1 Molecular Dynamics (MD) Simulation Molecular dynamics (35-37) simulation play an important role in understanding the motions of atoms and molecules. It is done by solving Newton’s equation of motion to determine the positions and velocities under a given interatomic force field. The position vector could be expressed as and the velocity vector is described by . By numerically solving the Newton’s equations of motion, the trajectories of each atom are predicted. ⃗ () = ⃗ () (2.1) , where mi is the mass of the i th particle. Considering a trajectory, ⃗ (), and let ∆ be a small increment in time, the velocity of i th particle can be obtained as, 6 ⃗ () = ⃗ = lim ∆→0 ⃗ ( + ∆ 2 ) − ⃗ ( − ∆ 2 ) ∆ (2.2) The acceleration of the i th particle can be inferred from the change in velocity during that time, and can also be approximated by three consecutive positions with small increments, ⃗ () = lim ∆→0 ⃗⃗(+ ∆ 2 )−⃗⃗(− ∆ 2 ) ∆ (2.3) We can update the positions and velocities of a system of N particles with given interatomic force fields according to the above equations. Verlet algorithm and Velocity-Verlet algorithm are two most frequently used algorithms for updating positions and velocities. Verlet algorithm: ⃗ ( + ∆) = ⃗ () + ⃗ ()∆ + 1 2 ⃗ ()∆ 2 + 1 6 ⃗⃗ ⃛()∆ 3 + (∆ 4 ), (2.4) ⃗ ( − ∆) = ⃗ () − ⃗ ()∆ + 1 2 ⃗ ()∆ 2 − 1 6 ⃗⃗ ⃛()∆ 3 + (∆ 4 ), (2.5) By adding up and subtracting the above two equations, we have: ⃗ ( + ∆) = 2⃗ () − ⃗ ( − ∆) + ⃗ ()∆ 2 + (∆ 4 ), (2.6) ⃗ () = [⃗ (+∆)−⃗ (−∆)] 2∆ + (∆ 2 ) (2.7) Here is the Verlet algorithm for MD integration: Given ⃗ ( − ∆) and ⃗ (), 1. Compute ⃗ () as a function of {⃗ ()}, 2. ⃗ ( + ∆) ← 2⃗ () − ⃗ ( − ∆) + ⃗ ()∆ 2 , 3. ⃗ () ← [⃗ ( + ∆) − ⃗ ( − ∆)]/2∆ 7 Velocity-Verlet algorithm: ⃗ ( + ∆) = ⃗ () + ⃗ ()∆ + 1 2 ⃗ ()∆ 2 + (∆ 4 ), (2.8) ⃗ ( + ∆) = ⃗ () + 1 2 ∆(⃗ () + ⃗ ( + ∆)) + (∆ 3 ), (2.9) Given ⃗ () and ⃗ (), 1. Compute ⃗ () as a function of {⃗ ()}, 2. ⃗ ( + ∆ 2 ) ← ⃗ () + ∆ 2 ⃗ (), 3. ⃗ ( + ∆) ← ⃗ () + ⃗ ( + ∆ 2 ) ∆, 4. Compute ⃗ ( + ∆) as a function of {⃗ ( + ∆)}, 5. ⃗ () ← [⃗ ( + ∆) − ⃗ ( − ∆)]/2∆ Velocity-Verlet algorithm is used in our MD simulations for classical and quantum dynamics simulations. Statistical ensembles In MD simulations, ensemble represents the macroscopic conditions, which is a collection of the possible microstates with identical macroscopic properties of the thermodynamic system. A statistical ensemble can be characterized by a set of conserved thermodynamic observables such as volume , energy , temperature , and pressure . Three common ensembles are microcanonical (NVE) ensemble, canonical (NVT) ensemble, isothermal-isobaric (NPT) ensemble, and grand canonical (VT) ensemble. Though NVE ensemble is the most commonly used due to its simplicity, the choice of the ensemble for MD simulations is based on the specific experimental purpose. 2.2 Interatomic Force Field 8 In molecular dynamics, it is essential to compute the force on each particle which determine the position and velocities at the next time step. The atomic forces are usually considered as derived from the total potential energy. Let N denote the total number of atoms present in the system, then the potential energy can be written as = (1 ,… . . ) = (1 , 1 , 1 , … . . , , , … . . , , , ) (2.10) Force on the i th atom is ⃗⃗ = − (1 ,… . . ) ⃗ = − ( ⃗ , ⃗ , ⃗ ) (2.11) In general, the total interatomic potential energy could contain multiple terms corresponding to a specific physical effect. 2.2.1 Lennard-Jones potential For an interatomic potential used for inert gases, which is Lennard-Jones potential, shown in Figure 2.1. The relation between the potential energy and distance in dimensionless units could be expressed as () = ∑() = < ∑4 [( ) 12 − ( ) 6 ] < (2.12) 9 Figure 2. 1. Dimensionless Lennard-Jones potential. 2.2.2 Reactive force field (ReaxFF) A first-principle informed many-body potential, ReaxFF (38) is used for reactive molecular dynamics simulations. ReaxFF is an interatomic force field that is designed to accurately describe materials properties as well as chemical reactions and found many scientific and engineering applications. ReaxFF describes the bond breaking and formation based on the bond-order concept and the charge transfer between atoms using Charge Equilibrium (QEq) scheme. ReaxFF has been used to simulate diamond (39), metal oxides (40-49) and organic-inorganic interfaces (50-52). 10 Figure 2. 2. Energy terms included in ReaxFF where Ebond is the bond energy, Elp is the lone pair energy, Eover and Eunder are the over- and under- coordination energies, Eval is the valence angle energy, Epen is the penalty energy, Ecoa is the three-body conjugation energy, Etors is the torsion angle energy, Econj is the four-body conjugation energy, Ehbond is the hydrogen bond interaction energy, EvdWaals, ECoulomb and Elg are the van der Waals, Coulomb energies and dispersion energy, respectively. The energy terms for the potential could be divided into two classes, one being bond-orderdependent and the other being bond-order-independent, as shown in Figure 2.2. The bond-orderdependent terms require the bond order value, which is applied for the chemical bond formation and dissociation. For the bond-order-independent terms, ReaxFF considers shielded van der Waals (Morse Potential) and Coulomb interactions, EvdWaals and ECoulomb, with a finite range (53). EvdWaals is calculated between every atom. ECoulomb is also calculated between every atom, with polarizable charges updated every iteration. For atom i and atom j, let be the interatomic distance, then the bond order, , is calculated from the following equation: = + + = [1 ( ) 2] + [3 ( ) 4] + [1 ( ) 6], (2.13) ReaxFF energy: E = Elp+Eover+Eunder + Ebond + Eval+Epen+Ecoa Non-bonded + Etors+Econj + Ehbond + EvdWaals+ECoulomb Bonded 11 where terms are equilibrium bond lengths, and terms are empirical parameters. Bond orders are updated every iteration. Among the energy terms, Ebond, Eover, Eval, Etors, EvdWaals, ECoulomb are some general energy terms for the total energy of the system, Esystem. Different systems choose energy terms differently: (1) For noble gases, Esystem = EvdWaals; (2) for ionic materials: Esystem = EvdWaals + ECoulomb; (3) for metals: Esystem = Ebond + Eover + EvdWaals; and (4) for metal alloys: Esystem = Ebond + Eover + EvdWaals + ECoulomb. For special systems, other terms might be included for specific research interests: (1) EH-bond is for hydrogen bond correction, (2) Elp, lone pair energy, is penalty for breaking up lone pairs in O,N; (3) Eunder is for undercoordination energy which stabilizes undercoordinated atoms; (4) Epen is penalty for “allene”-type molecules (H2C=C=CH2); (5) Ecoa is for angle conjugation which stabilizes -NO2 groups; (6) Econj is torsion conjugation for general conjugation stability; (7) EC2 is C2 correction which destabilizes C=C. 2.3 Periodic Boundary Conditions (PBC) Periodic boundary conditions is used to remove surface effects in MD simulations. The motivation of using PBC is to accurately simulate bulk properties in a relatively small system. When an atom crosses an PBC boundary cell, an exact duplicate of that atom enters the cell through the opposite face, thus no physical walls exist at that boundary, preserving a constant number of atoms in the system. 2.4 Reactive Molecular Dynamics (RMD) 12 We use ReaxFF which can describe bond dynamics with near quantum accuracy for performing reactive molecular dynamics simulations. ReaxFF allows for RMD simulations on systems containing more than 1,000 atoms, which allows larger simulations scales than QM-methods. 2.5 Density Functional Theory (DFT) Based on Hohenberg-Kohn theorems, the density functional theory focuses on solving Kohn-Sham equations to obtain electron density and energy. In 1964, Hohenberg and Kohn proved that the electronic ground state is a functional of the electron density (54). They demonstrated that the ground-state properties of a many-electron system could be uniquely determined by an electron density reliant on merely three spatial coordinates. It laid the foundation for reducing the complex many-body problem of N electrons with 3N spatial coordinates to just 3 spatial coordinates, through utilizing functionals of the electron density. This theorem has been expanded to the time-dependent domain to develop time-dependent density functional theory (TDDFT), which enables the description of excited states. They also defined an energy functional for the system and proved that the ground-state electron density minimizes this energy functional. In Kohn-Sham equation, the Hamiltonian consists of the kinetic and potential terms, (− ℏ 2 2 Δ 2 + ()) (r) = (r), (2.14) where () is constructed by external potential (), the Hartree potential (), the exchange– correlation energy (). The exchange-correlation functional is to describe many-body effects beyond the Hartree approximation. Common exchange-correlation functional includes GGA (55), MetaGGA (56), Hybrid functional (57), DFT+U (58), DFT-D (59), etc. Pseudopotential formulation is used to describe the valence electrons (60). 13 Self-consistent field (SCF) iterations are performed to obtain Kohn-Sham wave functions and energies. 2.6 Quantum Molecular Dynamics (QMD) In QMD nuclei follow the classical molecular dynamics while interatomic forces are computed quantum mechanically within DFT scheme, and using the Hellmann-Feynman theorem (61). Apart from software like VASP, QXMD software is also used for performing quantum molecular dynamics simulations, and the nonadiabatic quantum molecular dynamics simulations for photoexcitation-induced GeTe amorphization project are performed using QXMD software (62, 63). QXMD use a NSC method to evaluate accurate forces at a moderate computational cost for an excited electronic state. In QXMD, the NSC Harris-Foulkes energy ENSC is defined as (64-66) ENSC = f s ys Hˆ KS ys s ys ys å - Edc , (2.15) where Edc is double-counting terms related to the Hartree and xc energies: Edc = - 1 2 ò drr(r)VH (r)+Exc [r] - ò drr(r)Vxc (r). (2.16) ENSC coincides with the energy term EKS at the time when the electron density reaches selfconsistency. The self-consistency will be judged whether the difference between in(r) and out(r) falls in a given tolerance. Even if the self-consistency is not reached, the energy term ENSC deviates from the self-consistent energy only by an amount of second order in the deviation of in(r) from the self-consistent density (r), while EKS deviates by an amount of first order. 14 The force Fk acting on the k-th atom is Fk = - ¶ ¶Rk ENSC + zl zm l<m Rl -Rm å æ è ç ç ö ø ÷ ÷ = Fk Hellmann-Feynman +Fk NSC , (2.17) where Rk and zk are the position and the valence, respectively, of the k-th atom. The NSC force Fk NSC is given by Fk NSC = drdr¢ dr(r¢) r - r¢ - ¶ ¶Rk rin (r) ì í î ü ý þ òò + drdr(r) - ¶ ¶Rk Vxc (r) ì í î ü ý þ ò , (2.18) where δ(r) = out(r) − in(r). If the electron density achieved self-consistently, since δ (r) = 0, there is no contribution from Fk NSC. To alleviate the computational cost of excited-state forces calculation, QXMD first uses a diagonal approximation Fk = - ¶ ¶Rk wI + zl zm l<m Rl - Rm å æ è ç ç ö ø ÷ ÷ ~ - ¶ ¶Rk CI,ais ias å 2 Fais Hˆ N Fais + zl zm l<m Rl - Rm å æ è ç ç ö ø ÷ ÷ , (2.19) QXMD evaluate the first term of Eq. (2.19) using the KS orbitals y¢ s and their occupation numbers f s ais in Φaiσ as - ¶ ¶Rk Fais Hˆ N Fais = - ¶ ¶Rk f s ais y¢ s Hˆ KS y¢ s y¢ s y¢ s - Edc s å ì í ï îï ü ý ï þï = Fk Hellmann-Feynman +Fk NSC,ais . (2.20) 15 QXMD used r¢ ais (r) = f s ais y¢ s (r) s å 2 . (2.21) as out(r), and the ground-state electron density as in(r) to evaluate Fk NSC,ais . The developers of QXMD propose a modified approach based on a mixing scheme for the electron density: They employ a mixing scheme for Fourier components (G) of (r) proposed by Kerker (67): r¢ ais (G) = rin (G)+ AG 2 G 2 + B r¢ ais (G)- rin [ (G)], (2.22) where A and B are parameters to be optimized. Based on a modified Harris-Foulkes approach, NSC excited-force is calculated using Eqs. (2.19) and (2.20), where r¢ ais (G) in Eq. (2.22) is used as out and the ground-state electron density as in(r) to evaluate the value of Fk NSC,ais in Eq. (2.18). The QXMD simulation code has been parallelized by a hybrid approach combining band decomposition and spatial decomposition. The iterative band-by-band minimization is carried out by the band decomposition. The electron density is also calculated by the band decomposition. In QXMD, the Gram-Schmidt orthonormalization of the KS orbitals is executed in the reciprocalspace decomposition scheme, while the Fourier components of the wave functions are distributed among multiple processors, thus all-to-all communications are essential to switch between the two schemes. Global communication is used to calculate the scalar products between the wave functions. Message passing interface (MPI) library is used for interprocessor communications. 16 2.7 Nonadiabatic Quantum Molecular Dynamics (NAQMD) To investigate photoexcitation effect, we use nonadiabatic quantum molecular dynamics simulation (63, 68-74). The method involves exciting an electron from ground state to an excited state and follow the motion of electrons and nuclei ae simultaneously. Nonadiabatic quantum molecular dynamics simulation is based on the fewest-switches surface-hopping (FSSH) approach (75, 76). Electrons occupy one of the excited eigenstates corresponding to the atomic configuration at the current time instance. Transitions between the excited states are guided by the transition probability determined by non-adiabatic coupling (NAC). The NAC between excited states is depicted by a density matrix, and its time evolution is calculated using time-dependent density functional theory (63, 77-82). 2.8 Visualization Molecular dynamics simulations and density functional theory calculations sometimes generate vast amounts of data, visualization and analysis tools for the simulations and calculations results play significant roles to understand the atomic and electronic mechanisms. There are several packages designed for this demand, such as OVITO pro (83), Vesta ((84)), VMD (85), Visit (86), Atomeye (87) and Paraview (88), etc. The choice of software tools will depend on the specific needs and preferences. 2.9 Machine Learning and High Performance Computing Machine learning has made significant development in the field of materials simulations, serving as a powerful tool and offering alternative approaches to accelerate the discovery and design of materials. The application includes property prediction, materials discovery, data generation and 17 augmentation, structure prediction, molecular design, high-throughput screening and materials informatics. Machine learning based interatomic potentials are used to describe the interactions between atoms in a material, which provide accurate and computationally efficient alternatives to traditional force fields. Moreover, active learning techniques can provide guidance for simulations or experiments to reduce uncertainty and improve the model's accuracy. High performance computing plays an important role in fully utilizing the computational capacity of computers. To enable larger and longer molecular dynamics simulations, it is necessary to apply high performance computing for a better performance of the programs. My collaborated work about using machine learning and about using high performance computing are published in the following papers: 1. Nazarova. A. L.; Yang, L.; Kuang, L.; Mishra A.; Kalia R. K.; Nomura K.; Nakano, A.; Nakano, A.; Vashishta, P., Dielectric Polymer Property Prediction Using Recurrent Neural Networks with Optimizations. Journal of Chemical Information and Modeling 2021, 61, 5, 2175-2186. 2. Ibayashi, H.; Razakh T. M.; Yang, L.; Linker, T.;Olguin, M.; Hattori, S.; Luo, Y.; Kalia R. K.; Nakano, A.;Nomura, K.; Vashishta, P., Allegro-Legato: Scalable, Fast, and Robust Neural-Network Quantum Molecular Dynamics via Sharpness-Aware Minimization. ISC High Performance 2023: High Performance Computing, 223–239. 18 Chapter 3: ZrS2 Oxidation In this section, most content are already published in the following papers: 1. Yang, L.; Jaramillo, R.; Kalia R. K.; Nakano, A.; Vashishta, P., Pressure-Controlled Layerby-Layer to Continuous Oxidation of ZrS2 (001) Surface. ACS Nano 2023, 17, 8, 7576– 7583 (89); 2. Yang, L.; Tiwari, S. C.; Jo, S. S.; Hong, S.; Mishra, A.; Krishnamoorthy, A.; Kalia, R. K.; Nakano, A.; Jaramillo, R.; Vashishta, P., Unveiling oxidation mechanism of bulk ZrS2. MRS Advances 2021, 6, 303-306 (90); 3. Jo, S. S.; Singh, A.; Yang, L.; Tiwari, S. C.; Hong, S.; Krishnamoorthy, A.; Sales, M. G.; Oliver, S. M.; Fox, J.; Cavalero, R. L.; Snyder, D. W.; Vora, P. M.; McDonnell, S. J.; Vashishta, P.; Kalia, R. K.; Nakano, A.; Jaramillo, R., Growth Kinetics and Atomistic Mechanisms of Native Oxidation of ZrSxSe2–x and MoS2 Crystals. Nano Letters 2020, 20 (12), 8592-8599 (91). 3.1 Introduction Transition metal dichalcogenides hold great promise for their prospective utilizations in future electronic and optoelectronic devices. Among TMDCs, semiconducting ZrS2 (9, 92) is especially suitable for microelectronics (10, 93, 94). In addition, ZrS2 has been studied for solar-energy applications due to its effective bandgap of 1.7 eV for absorbance of photons from visible-light region (400-760 nm) (93, 95-97). However, native oxidation poses a significant challenge for many TMDCs in ensuring their long-term stability, and ZrS2 is no exception from this concern (15, 89). To control the native oxidation, it is crucial to gain a comprehensive understanding of the atomistic 19 oxidation mechanisms of ZrS2. Moreover, oxidation products themselves hold great scientific and technological value. ZrO2 nanocrystals has shown chemical inertness, excellent thermal stability, high hardness (98, 99), and excellent optical properties, which make them excellent materials for fuel cells (100, 101), catalysis (102-104), sensors (105, 106), bioseparation (107), chromatography (108), and high-refractive-index nanocomposites (109). Tetragonal zirconium oxysulfide (ZrOS) nanopowder synthesized by sol-gel method was suggested for dielectric material and thin film photovoltaic applications (110). First-principles study by Zhang et al. (111) revealed that Janus monolayer ZrOS exhibits large dielectric permittivity, mechanical and dynamical stabilities, high and sharp absorption peaks in visible and ultraviolet (UV) light range, showing promising applications for dielectric semiconductor materials and visible and UV light sensor. The oxidation processes of ZrS2 have been studied recently. Li et al. (25) investigated the oxidation of monolayer ZrS2 and observed that surface vacancies and edges show high affinity towards the adsorption and activation of oxygen molecules. We developed new reactive force field parameters for Zr/O/S, and collaborated with Jo et al. (89) to study the oxidation mechanism. We observed a higher oxidation rate of ZrS2 than that of MoS2 and identified the key atomistic mechanisms: Rapid O2 adsorption and bond scission, followed by oxygen transport into the crystal via Zr-O bond switching and the collapse of van der Waals (vdW) gaps, leading to the formation of an intermediate amorphous oxy-sulfide. The further optimized force field parameters using multiobjective genetic algorithm (MOGA) (112, 113) are able to reproduce quantum-mechanically computed charge values and bond-population dynamics. Using RMD simulations with the optimized reactive force field parameters, we found that the oxidation is diffusion-controlled in the fast oxidization period (89, 90). It is also observed that (90) oxidation rate is higher on (210) 20 surface compared to that on (001) surface, and the oxidation in the fast oxidization period on both surfaces were identified as diffusion-controlled process. An outstanding and unsolved concern is how oxidation behaves under different environments (114), which is crucial for better device processing and future technology development. In humid environment, it has been proposed that H2O plays a role in the degradation of ZrS2 not only because an H2O molecule can adsorb on the metal to achieve continuous oxidation, but also because its robust polarization enhances the adsorption of nonpolar oxygen molecules (25, 115). To synthesize oxide-semiconductor heterostructures in emerging two-dimensional electronics, active control of oxidation is being vigorously explored, including plasma-, ultraviolet- and laser-assisted oxidation processing (116, 117). In such advanced growth chambers, active control of oxygen partial pressure can be used as another growth-control parameter. We performed RMD simulations to investigate how oxygen partial pressure affects oxidation of the ZrS2 (001) surface at atomistic level. The current work suggests that oxygen partial pressure is indeed a promising control parameter for such active oxidation for future 2D electronics. The oxidation rate is found to increase under higher O2 partial pressure, which is consistent with the kinetics. We observed that, as oxidation progresses, a crossover from layer-by-layer oxidation to amorphous-oxide-mediated continuous oxidation take place. During the first stage, layer-by-layer oxidation, the oxidation proceeds with reactive bond switching combined with rotation events, leading to the closure of vdW gaps. During the second stage, the kinetics of continuous oxidation and its pressure dependence are well-described by the conventional Deal-Grove model (118). Varied pressures are found to selectively expose distinct phases of oxide growth within a specified time frame, where the oxidation process transitions from layer-by-layer oxidation via reaction-limited linear to amorphous-oxide-mediated oxidation via diffusion-limited parabolic growth as pressure is raised. We also identified a trade-off between the 21 oxide growth rate and the oxide structure quality, along with the impact of pressure on the morphology of semiconductor/oxide interfaces. Overall, our findings investigated and showed the ability of adjusting pressure to effectively control the oxidation of layered vdW semiconductors at atomistic level. 3.2 Method The method for this section is put in relevant separate subsections. 3.3 Force Field Optimization 3.3.1 Training ReaxFF Parameters ReaxFF parameters are optimized in two steps. In the first stage, we optimized the reactive force field parameters to accurately describe the ground-state atomic charges of Zr, O and S atoms. In the second stage, we further optimized the reactive force field parameters for describing Zr-S and Zr-O bond-population dynamics. Training ReaxFF parameters for charges: To construct a training set for optimizing charges, we built 53 small clusters as shown in Figure 3.1. Mulliken charges were computed by performing DFT calculations using QChem software package (119). In DFT calculations, we used the B3LYP hybrid functional (120, 121), with dispersion correction from Grimme’s DFT-D3 method implemented (59). S and O atoms were treated with 6-31+G** basis set, while heavier Zr atoms were treated with LANL2DZ pseudopotential method. Mean square difference between ReaxFF charges and DFT-computed atomic charges was used as objective function in this training stage to optimize the ReaxFF parameters. 22 Figure 3. 1. Training set for optimizing charges. Training ReaxFF parameters for bond-population dynamics: To refine the reactive force field parameters to better approximate the ground-truth Zr-S and Zr-O bond-population dynamics in QMD simulations of ZrS2 oxidation, RMD simulations were performed with the same schedule for QMD. QMD simulations were performed using VASP (122), while RMD simulations were performed using the created reactive force field parameters in LAMMPS package (123). The same initial positions and velocities were provided in both QMD and RMD simulations for a Zr48S96 slab immersed in 72 O2 molecules. We performed reactive molecular dynamics simulations in the NVT ensemble at 1,200 K for 11.155 ps. The equations of motion were integrated numerically 23 using the time step of 1 fs, and the numbers of Zr-S, Zr-O, S-S, and S-O bonds were sampled at every time step. The numbers of bonds were recorded as a function of time for both RMD and QMD simulations. The root mean square of the deviation in the number of each bond was tabulated as a 4-column list, which was used as the multi-objective function to train the ReaxFF parameters. 3.3.2 Workflow for Force Field Training For multi-objective genetic algorithm training of ReaxFF parameters, we used the NSGAII algorithm (124) implemented in the EZFF force field training software (125). The mutation probability was set to be 0.2, and the population size was chosen to be 1,024 in this work. The initial population was generated using a random generator, where the parameter values were set within the specified range. The evolution process consists of four steps in each generation as shown in Figure 3.2: (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; and (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. This completes one iteration of the NSGAII algorithm. The stopping rule here was a maximum number of 200 iterations or no significant improvement in the fitness value for some fixed number of generations. 24 Figure 3. 2. Workflow for training bonds population. Since we need a large population size for the force field optimization, the workflow will encounter the file I/O bottleneck. To resolve the I/O issue, we utilized the client-server model employed by Mishra et al. (126). On a parallel computer, we let the head computing node, which acts as the server, run the main workflow script: Computing the charges or bond-populations and finally estimating the error corresponding to every client process. After the objective function is computed, the server writes all the error values to a file and executes the genetic algorithm procedure to produce a set of ReaxFF parameter offspring. The workflow spawns message passing interface (MPI) jobs (127), one per each force field to perform RMD simulation and then gives the information for error evaluation. 3.3.3 Force Field Accuracy Accuracy of the optimized ReaxFF parameters on reproducing charges: In the first stage of optimization, we train ReaxFF parameters for Zr, S, and O atomic charges. To evaluate the accuracy of the optimized force field parameters on reproducing atomic charges, we compared the ReaxFF charges and QM charges between the values obtained by ReaxFF description and the values by DFT calculation in QChem software. Figure 3.3 shows insignificant differences between 25 the two calculations, which confirms that the charge values obtained by optimized force field are consistent with DFT calculated charges. Figure 3. 3. Comparison of QM charges (black color) and ReaxFF charges (red color) for selected clusters. Green, yellow and red spheres represent Zr, S and O atoms, respectively. Accuracy of the optimized ReaxFF parameters on reproducing bond-population dynamics: In the second stage of optimization, we train Zr-S and Zr-O 2 two-body terms to correctly describe the bond breaking and bond formation during the oxidation reaction of ZrS2. Figure 3.4 shows the error (red curve for Zr-S, black curve for Zr-O) and uncertainty (shaded orange region for Zr-S, azure region for Zr-O) in the number of Zr-S bonds and Zr-O bonds during training in each epoch of genetic algorithm. The shaded area represents for the standard deviations in the number of bonds in the RMD simulations of the global Pareto optimals. After 130 iterations, we stopped the training process because error did not decrease further. 26 Figure 3. 4.Reduction of bond error as a function of the epoch in multi-objective genetic algorithm training. Figure 3.5 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 genetic-algorithm refinement of ReaxFF parameters. After the parameter tuning, Figures 3.5b, d exhibit much better agreement between QMD and RMD simulations compared with those before the training in Figures 3.5a, c. RMD simulations with the refined ReaxFF parameters thus provide quantitatively correct bondpopulation dynamics consistent with the ground-truth QMD. Figure 3. 5.Comparison of QMD (blue) and RMD (red) simulation results for time evolution of Zr-S (a, b) and Zr-O (c, d) bonds, both with initial (a, c) and optimized (b, d) RMD force fields. 27 As a comparison, we have also tested other genetic-algorithm variants like NSGAIII and epsNSGAII for the force field optimization. From Figure 3.6, we see that the force field trained by NSGAII algorithm is accurate enough to be used in RMD simulations to study the oxidation of ZrS2. Figure 3. 6. Comparison for the performance of force fields optimized by different algorithms. Time evolution of Zr-S and Zr-O bonds are plotted for QMD (blue) as well as for RMD (red) with force fields using (a) NSGAIII, (b) epsNSGAII, and (c) NSGAII algorithms. 3.4 Oxidation Mechanism 3.4.1 Comparison of Thermal Oxidation between MoS2 and ZrS2 The content in this subsection a collaborative work, and we are thankful to valuable experimental insights from experimentalists Prof. Rafael Jaramillo and Dr. Seong Soon Jo at MIT. 28 Method: To understand the rapid oxidation of Zr-based TMDs compared to MoS2, we performed RMD simulations of the oxidation of MoS2 and ZrS2 slabs. We simulated oxidation in an oxygen partial pressure of 1.16 kbar and at 800 K, to accelerate the process; as a result, the timescales of oxidation observed in simulation is not expected to match those observed in experiment. In Figures 3.7a, b we show the initial simulation setup (at time t = 0 ns) and the final snapshot (t = 2.5 ns) for MoS2 oxidation. The simulation box is periodic in the 〈100〉 and 〈010〉 directions (i.e., in-plane). We employed a similar setup for simulating ZrS2 oxidation, as shown in Figures 3.7c, d. Oxidation of MoS2 and ZrS2: Figures 3.7a-d show that ZrS2 oxidizes extensively compared to MoS2. To quantify the extent of oxidation, we show in Figures 3.7e, f the total numbers of M-S, M-O and O-O bonds (M = Mo or Zr). The increase in Mo-O with oxidation time is negligible compared to the rapid increase of Zr-O bonds. Concurrently, we observe a notable decrease in O-O bonds in ZrS2 slab, suggesting scission of O-O bonds followed by formation of Zr-O bonds. The bond counting data in Figure 3.7f suggest that ZrS2 oxidation progresses in two stages. During the first 0.1 ns, the number of Zr-O bonds rises rapidly, followed by much slower increase at longer times. The rapid initial reaction may be explained by energetically-favored O2 adsorption on the ZrS2 surface. The rapid increase in the number of Zr-O bonds is not accompanied by a similarly rapid decrease in the number of Zr-S bonds, as we discuss below. The subsequent, slower process of ZrO bond formation represents the growth of an amorphous oxy-sulfide that largely eliminates the vdW gaps. In contrast, we observe almost no decrease of O-O bonds in the MoS2 simulation. This suggests that MoS2 remains mostly inert and stable in the simulated environment. 29 Figure 3. 7. RMD simulations of MoS2 and ZrS2 oxidation. Snapshots of (a-b) MoS2 oxidation and (c-d) ZrS2 oxidation at time t = 0 and 2.5 ns, respectively. Spheres represent individual atoms: Blue (Mo), green (Zr), yellow (S) and red (O). Time evolution of corresponding bond count for (e) MoS2 and (f) ZrS2 oxidation. We analyze the process of ZrS2 oxidation in detail by considering atomic trajectories. We find that oxidation proceeds in six steps: (1) Adsorption of O2 molecules on the ZrS2 surface; (2) dissociation of O2 on ZrS2 surface, and formation of Zr-O bonds; (3) Zr pushed into vdW gaps between TMD layers due to formation of new Zr-O bonds; (4) formation of inter-layer Zr-S and Zr-O bonds, resulting in loss of the ZrS2 crystal structure (i.e. amorphization); (5) oxygen transport laterally within and vertically between layers assisted by Zr-O bond switching; and (6) formation and out-diffusion of SO2. In Figures 3.8a, b we show the cross-section and top-down view of the simulation at t = 2 ps, illustrating steps 1-2. We indicate O2 adsorption sites with black rectangles. 30 We find that the adsorption of O2 on the ZrS2 basal surface is energetically favorable, with adsorption energy Eads = -0.006 eV per O2 molecule, which we calculate by DFT simulations of a ZrS2 surface-O2 molecule system. In contrast, we find no stable, relaxed energetic minimum for a MoS2 surface-O2 molecule system, indicating that the adsorption energy is positive. Constrained MD simulations paint a similar picture, that O2 adsorption is comparatively favorable on ZrS2 but unfavorable on MoS2. After rapid adsorption of O2, we observe O-O bond breaking and the formation of Zr-O bonds. We also observe 2- and 3-coordinated oxygen atoms on the ZrS2 surface after O-O bond breakage, corresponding to 2-coordinated bridge oxygen, and 3-coordinated oxygen that are seeds for zirconia growth. A similar process was proposed for the oxidation of monolayer ZrS2 (25). Following the initial formation of Zr-O bonds, Zr is pushed out of the top layer and into the interlayer region. The protrusion of Zr atoms results in the collapse of the vdW gap between ZrS2 layers by the formation of new, interlayer Zr-S and Zr-O bonds, resulting in amorphous, nonlayered material. This is illustrated by a simulation view at t = 1 ns, presented in Figures 3.8c-d. The close-up in Figure 3.8d highlights the inter-layer Zr atoms and bonds. The expansion of these disordered regions into the vdW gaps provides pathways for oxygen transport further into the ZrS2 crystal. We find that oxygen transport proceeds by a process of Zr-O bond switching, illustrated step-by-step in Figure 3.8e, in which Zr and O atoms are numbered, and S atoms are omitted for clarity. The coordination of oxygen atom O1 first changes from 3 to 2 by the breakage of the Zr1- O1 bond. The subsequent creation of a new Zr4-O1 bond restores the 3-fold coordination. A similar process continues when the Zr2-O1 bond is broken and a new Zr5-O1 bond is formed. Oxygen thus traverses the system by successive breaking and forming Zr-O bonds. This bond-switching mechanism for oxygen transport is akin to Grotthuss mechanism which plays essential roles in 31 various processes in oxides (128-130). It is notable that the process of Zr-O bond switching and oxygen transport is not correlated with a decrease in the number of Zr-S bonds. This suggests that S atoms act as spectators during these initial steps of Zr-O bond formation, amorphization, and oxygen transport, and that charge neutrality is maintained by successive redox transitions at the S sites. These processes are also observed in first-principles QMD simulations, thus validating these atomistic mechanisms of ZrS2 oxidation. Figure 3. 8. Atomistic mechanisms of ZrS2 oxidation. Cross section (a) and top-down (b) views at t = 2 ps, showing O2 adsorption and Zr-O bond formation that involves both 2- and 3-fold coordinated oxygen atoms (enclosed by black rectangles). Cross section (c) and close-up (d) view at t = 2.5 ns, showing the formation of an amorphous matrix and the closing of the vdW gaps; S atoms are not shown in (d) for clarity. (e) Oxygen transport mechanism of Zr-O bond switching, including oxygen switching between 2- and 3-fold coordination: (i-ii) Breakage of Zr1-O1 bond and decrease in coordination number (3→2); (iiiii) formation of Zr4-O1 bond and increase in coordination number (2→3); (iii-iv) breakage of Zr2-O1 bond and reduced coordination number (3→2); (iv-v) formation of Zr5-O1 bond and coordination number recovery (2→3). 3.4.2 Comparison of Oxidation on ZrS2 (210) and (001) Surfaces Method: In order to compare the oxidation processes of ZrS2 slabs on (210) and (001) surfaces, we performed RMD simulations in NVT ensemble using the refined ReaxFF reactive force field 32 parameters. We applied periodic boundary conditions to both systems. To accelerate the oxidation process within accessible simulation time, the oxygen partial pressure was set to be 1.26 kbar and the reaction temperature was chosen to be 1,200 K. The initial configurations for both simulations are displayed in Figure 3.9a. Each system contained a six-layer ZrS2 slab with (210) or (001) surface exposed to O2 atmosphere. Oxidation on (210) and (001) surfaces: Figure 3.10b 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 on the right (210) surface, where zirconium atoms are directly exposed to O2 atmosphere and the adsorption of oxygen to Zr sites is thus possible. Consequently, the rate of oxidation is highly contingent on the initial adsorption, showing an anisotropic behavior. The subsequent oxidation pathway in both directions involves adsorption of oxygen on ZrS2 surface, followed by amorphization and the transport of oxygen into the bulk, resulting in a continuous oxidation. 33 Figure 3. 9. (a) Initial system setups of ZrS2 slabs with (210) and (001) surfaces. Green, yellow and red spheres represent Zr, S and O atoms, respectively. (b) Snapshots of oxidization process in both (210) and (001) simulations at times, t = 0.02, 0.2, 0.3 and 0.4 ns. Figure 3. 10. Oxide depth as a function of time on different surfaces: The right and left (210) surfaces (black) as well as top and bottom (001) surfaces (red). To provide a more detailed examination of the oxidation mechanism, we illustrate the evolution of the oxide growth depth at different surfaces in Figure 3.10. The general oxide growth trend is similar for both the right and left (210) surfaces, which is much faster than both top and bottom (001) surfaces. As previously discussed, the zirconium atoms are directly exposed to the O2 atmosphere on the right (210) surface, making them easier to be oxidized. According to Li et al (25), the adsorption energy of O2 on (210) edges tends to be more negative than on the pristine (001) surface for TMDCs. In their study, the pristine surface corresponds to the (001) surface mentioned here and the edge sites are situated on the (210) surface, though they considered Zredge with 50% S coverage, where the Zr edges are coordinated with five sulfur atoms. They calculated the adsorption energy of O2 on the pristine ZrS2 surface and the edge sites, which yielded values of -0.06 eV and -0.92 eV, respectively. Since Zr-edge coordinated with less sulfur atoms have more negative adsorption energy and more energy released by oxidation reaction, the (210) 34 surface oxidized much faster, where the zirconium atoms on the right panel of (210) surface are coordinated with 3 sulfur atoms. On the (001) surface, the oxide grows very slowly before a sudden increase at 0.2 ns. The increase is caused by the formation of a disordered ZrS2 region between ZrS2 layers, which facilitates extensive oxidation through the oxygen transport, as depicted in Figure 3.10. Both the top and bottom (001) surfaces show a similar oxidation trend later from 0.3 ns to 0.4 ns. While a similar fast oxide growth period is observed for both the (001) and (210) surfaces, the oxidation rate is slightly faster on (210) surface in that period. This difference can likely be attributed to a geometric effect, where spacing between ZrS2 layers provides channels for oxygen diffusion in the case of the (210) surface, facilitating faster oxidation. While this work compared oxidation behaviors on different ZrS2 surfaces and found that ZrS2 (001) surface oxidizes slower than ZrS2 (210) surface (90), the most compelling interest in ZrS2 oxidation (i.e., synthesizing oxide-semiconductor heterostructures for emerging 2D electronics) is predominantly concerned with (001) surfaces. Therefore, in the following subsections, I will be mainly focusing on investigating the oxidation behavior on ZrS2 (001) surfaces. 3.4.3 DFT Calculations of the Gibbs Free Energy Change of Adsorption Binding Energy: The Gibbs free energy change due to adsorption is ∆(, ) = 1 (⁄ − − ∆ − ), (3.1) where / and are the Gibbs free energies of the oxygen-adsorbed surface and clean surface, respectively. and are the chemical potentials of Zr and O atoms. is the area of the surface. The chemical potential of oxygen, , and its dependence on pressure, p, and temperature, T, are presented in Table 3.1. The Gibbs free energy change due to adsorption is simplified as: 35 ∆(∆) = 1 (NO/2 − ∆ − ∆), (3.2) where /2 is the binding energy of oxygen atoms on ZrS2(001) surface. We used DFT in the VASP software (122, 131) to calculate /2 . ZrS2 was modeled using a supercell, where a fivelayered ZrS2 slab is constructed with a vacuum region of 25 Å to prevent interaction between periodic images. The oxygen ad-layer structures were modeled using 2 × 2 surface unit cells. We let oxygen atoms be adsorbed on both sides of the slab for a higher accuracy to eliminate asymmetry. The oxygen atoms and all atoms in the outer two ZrS2 substrate layers were fully relaxed, whereas the central three layers were fixed in their bulk positions. We used the projector augmented wave (PAW) method (132) and Perdew-Burke-Ernzerhof (PBE) functional (55). The energy cut-off for plane-wave expansion was set to be 520 eV. The total energy and forces were converged to within 0.1 eV per atom and 1 meV/Å, respectively. Brillouin zone was sampled over 5 × 5 × 1 Monkhorst-Pack k-point meshes (133). The average binding energy per oxygen atom adsorbed on the surface, /2 , was calculated as /2 = − 1 (⁄ − − ∆ − 2 2 ), (3.3) where and ∆ are the number of oxygen atoms and the change in number of zirconium atoms, while ⁄ , , , and 2 are the energy of the adsorbate-substrate system, the clean surface, the zirconium atoms, and the free oxygen molecule. 36 Figure 3. 11. Calculated binding energies of oxygen on Zr sites of ZrS2 (001) surface for various coverages. The binding energy profile shows a thermodynamically favored initial oxygen adsorption on the Zr sites of ZrS2 (001) surface as shown in Figure 3.11. The binding energy first increases as an increase in oxygen coverage due to the dominance of the attractive interaction between oxygen and the ZrS2 (001) surface. However, at coverages of 0.75 and 1 ML, the binding energy profile shows a decreasing trend as the oxygen coverage increases, which suggests that the system reaches a certain level of saturation or equilibrium, where repulsive forces between oxygen might lead to a less favorable overall interaction. Chemical Potential of Oxygen: Assuming that N oxygen molecules at pressure and temperature have ideal gas behavior, the chemical potential is obtained as = ( ) , (3.4) through the Gibbs free energy, = ( ) , + ( ) , + ( ) , = − + + . (3.5) From the ideal gas equation of state, = , and thus 37 ( ) , = = . (3.6) When only varying pressure, we have (, ) − (, 0 ) = ∫ ( ) , 0 = ln 0 , (3.7) and thus (, ) = 1 2 2 (, ) = 1 2 [2 + ̃O2 ( 0 , ) + ln 2 0 ], (3.8) where 0 = 1 atm. Taking a reference state of (, ) to be the total energy of oxygen in an isolated molecule (134), 2 = ̃O2 ( = 0 K, ) = 0, (, ) = 1 2 [̃O2 ( 0 , ) + ln 2 0 ] = ̃O( 0 , ) + 1 2 ln 0 . (3.9) Since the surrounding O2 atmosphere forms an ideal-gas-like reservoir (135), ̃O( 0 , ) = ̃ −ℎ ( 0 , = 0) + 1 2 ∆(∆, 0 ,2 ) = 1 2 [(, 0 , O2 ) − (0 K, 0 , O2 )] − 1 2 [(, 0 , O2 ) − (0 K, 0 , O2 )]. (3.10) ̃O2 ( 0 , ) can be obtained using the enthalpy and entropy in thermochemical tables at 0 = 1 atm (see the ̃O( 0 , ) values in Table 3.1) (135, 136). Table 3. 1. Chemical potential of oxygen (K) 0 100 200 300 400 500 600 700 800 900 1000 ̃O( 0 , ) (eV) 0 -0.08 -0.17 -0.27 -0.38 -0.50 -0.61 -0.73 -0.85 -0.98 -1.10 (eV) / 0.0086 0.0172 0.0258 0.0345 0.0431 0.0517 0.0603 0.0689 0.0775 0.0861 ln( ∗ [atm]) / 18.568 19.727 20.888 22.049 23.209 23.596 24.204 24.659 25.272 25.530 38 To characterize the pressure dependence of the chemical potential, we introduce ∗ , for which ( ∗ , ) = 0. Table 3.1 shows the calculated ∗ values at various temperatures. The calculated values are consistent with the numbers given in Soon et al. (134). At temperatures 300, 600 and 900 K, they reported the pressure values to be 109 , 1010 and 1011 atm, respectively. Note from Table 3.1, ∗ = 20.888 =1.179× 109 ≅109 atm at 300 K, 23.596 =1.768× 1010 ≅1010 atm at 600 K, and 25.272 =0.945× 1011 ≅ 1011 atm at 900 K. From Figure 3. 12. Pressure dependence of the change of Gibbs free energy upon oxygen-atom adsorption. Calculated change of Gibbs free energy, ∆, upon oxygen adsorption on ZrS2(001) surface with varying oxygen chemical potential ∆. The corresponding pressure at various temperatures is shown as bars below the plot. Apart from the binding energy profile discussed previously, Figure 3.12 shows the change of the Gibbs free energy of adsorption, ∆ , on the change of the oxygen chemical potential, ∆ , at different pressures and temperatures. Gibbs free energy of adsorption decreases as the increase of the oxygen chemical potential. The oxygen adsorption is thermodynamically favored even for a pressure far below atmosphere. 39 3.4.4 Pressure-Controlled Layer-by-Layer to Continuous Oxidation My interests towards pressure effectstart when I was an undergraduate student and was performing chemical vapor deposition experiments, where pressure is one of the main controlling factors. Therefore, it is also exciting for me to think about using molecular dynamics simulations to understand pressure effect of ZrS2 oxidation. Method: We performed RMD simulations to study the oxidation of a ZrS2 (001) slab in oxygen (O2) atmosphere with varying pressure at temperature 300 K as shown in Figures 3.13 a and b. RMD simulations in the NVT ensemble were performed using LAMMPS software (137). While NVT simulation exhibits the same sequence of oxidation events as in the NPT simulation, the former corresponds more closely to advanced active-oxidation chambers that operate in a batchreaction mode (116, 117). The ReaxFF reactive force field was developed and optimized for ZrS2 oxidation as discussed in previous subsections (89, 90). An 8-layer ZrS2 (001) slab with 1,536 Zr atoms and 3,072 S atoms was placed in the middle of the simulation box, where the Cartesian zaxis is normal to the slab. O2 molecules were placed above and below the slab to form a gas atmosphere, where the number of O2 molecule was calculated to achieve the desired pressure. Both ZrS2 and O2 were heated and fully relaxed at a temperature of 300 K. The system temperature was controlled using Nosé-Hoover thermostat (138, 139). Periodic boundary conditions were applied in all directions. 40 Figure 3. 13. Initial and final configuration of RMD simulations of ZrS2 oxidation. Initial configuration in the parallel front view (a) and perspective side view (b). Final configurations of ZrS2 oxidation at time 0.6 ns under pressure p = 0.1 GPa (c), 0.08 GPa (d) and 0.06 GPa (e). Green, yellow and red spheres represent Zr, S and O atoms, respectively. Oxide growth dynamics: Figures 3.13, c-e, show snapshots of the system after 0.6 ns of oxidation at oxygen pressures of 0.1, 0.08 and 0.06 GPa. Oxidation proceeds the fastest at the highest pressure of 0.1 GPa, followed by that under 0.08 GPa, and slowest at 0.06 GPa. 41 Figure 3. 14. Snapshots of RMD simulations of ZrS2 oxidation at different pressures and times at temperature 300 K. (a-c) ZrS2 oxidation after time t = 0.2 ns at pressures 0.1 GPa (a), 0.08 GPa (b) and 0.06 GPa (c). (d-f) ZrS2 oxidation after time t = 0.6 ns at pressures 0.1 GPa (d), 0.08 GPa (e) and 0.06 GPa (f). Oxygen atoms are colored yellow. Zr atoms are colored according to the number of coordinated oxygen atoms, i.e., Zr atoms with no coordinated oxygen atoms are colored red, while the color of Zr atoms coordinated to more oxygen atoms is graded toward blue. To better visualize the oxide growth, surfaces of the oxidation front are shown in cyan color, and S atoms are omitted for clarity. To better represent the pressure dependence of the oxidation rate, Figure 3.14 shows snapshots of ZrS2 oxidation at various pressures and times, in which oxidation fronts are shown as surfaces and sulfur (S) atoms are omitted for clarity. Oxygen (O) atoms are colored yellow, and zirconium (Zr) atoms are colored according to the coordination number to O atoms. At time 0.2 ns, Figures 3.14, a-c, show a similar initial oxidation stage for all pressures. The pressure dependence becomes evident by time 0.6 ns. At the highest pressure of 0.1 GPa, oxidation of the second layer is completed by 0.6 ns, and the oxidation front has reached the third layer (Figure 3.14d). In contrast, the oxidation front has just reached the second layer at 0.08 GPa (Figure 3.14e) and is still in the first layer at 0.06 GPa (Figure 3.14f). 42 Figure 3. 15. Pressure dependence of ZrS2 oxidation dynamics. Time evolution of the oxide depth (a), number of Zr-O bonds (b), and number of Zr-S bonds (c). Blue, red and black curves represent oxidation at pressures 0.1, 0.08 and 0.06 GPa, respectively. To quantify the pressure dependence of the oxidation rate, Figure 3.15a shows the oxide depth as a function of time. We define the oxide depth as the average position of oxygen atoms in the solid, taking the origin as the ZrS2 surface. At 0.1 GPa (blue curve in Figure 3.15a), oxide growth begins slowly, accelerates around 0.16 ns, and then gradually slows down, thus exhibiting a two-stage oxidation behavior. Fast oxidation in the second stage is characterized by a parabolic shape, indicative of a diffusion-controlled process, as seen in our published research (90). At a lower 43 pressure of 0.08 GPa (red curve in Figure 3.15a), we observe similar behavior: Oxidation is slow at first, then becomes faster around 0.3 ns. At the lowest pressure of 0.06 GPa, oxide growth remains slow throughout. The slow oxidation process at 0.06 GPa is characterized by stepwise increase of oxide depth, indicative of layer-by-layer oxidization, as found in our previous work (90). Here, oxidation of each layer begins with a kink, followed by flattening of the oxide front with the completion of the layer’s oxidation, resulting in the observed stepwise increase in oxide depth. The oxidation process exhibits similar layer-by-layer growth kinetics at short times (i.e., before the onset of fast growth) for all pressures. Similar kink mechanism (i.e., nucleation of a spatially-localized growth front, followed by lateral growth to form a flat surface) is characterized by a low activation energy, thus observed in many processes such as domain-wall motion (140) and crack extension (141). It appears that, after a period of layer-by-layer oxidation, a transition occurs to amorphous-oxide-mediated continuous oxidation, with the transition occurring at earlier times with increasing pressure. 44 Figure 3. 16. Time evolution of Zr-O pair distribution function. Blue, red and black curves are for pressures 0.1, 0.08 and 0.06 GPa, respectively, at times t = 0.1 ns (a), 0.2 ns (b), 0.3 ns (c), 0.4 ns (d), 0.5 ns (e), and 0.6 ns (f). 45 To probe chemical reactions underlying the crossover of oxidation stages, Figures 3.15 b and c show time evolution of the numbers of Zr-O and Zr-S bonds, respectively. Bond order is calculated for all Zr-O and Zr-S bonds, and bond is counted only when bond order is larger than 0.3. Figure 3.15b shows increase in the number of Zr-O bonds as the oxidation progresses. The increase is slowest at the lowest pressure of 0.06 GPa. At the highest pressure of 0.1 GPa, there is an abrupt increase in the number of Zr-O bonds at 0.16 ns, which is consistent with the onset of the fast oxidation in Figure 3.15a. A notable observation is that the number of Zr-O bonds increases faster at 0.08 GPa than at 0.1 GPa. This could be due to disruptive and incomplete oxidation within short time at the highest pressure of 0.1 GPa, which in turn produces low-quality defective oxide. To quantify this effect, Figure 3.15 shows Zr-O partial pair distribution at various times for all pressures. At time t ≤ 0.2 ns, the first peak is located at 1.75 Å indicative of Zr-O chemical bond at all pressures. After 0.3 ns under the highest pressure of 0.1 GPa, however, we observe the lowering of the first peak at 1.75 Å accompanied by the development of large correlation at distance larger than 2 Å , signifying mechanical disordering of chemical structures. This explains the decrease in the Zr-O bond number after 0.3 ns (Figure 3.15b), while the oxide is still growing (Figure 3.15a); note that Zr-O pairs farther than 2 Å are not chemically bonded nor counted in Figure 3.15b. In fact, Figure 3.14d shows some surface Zr atoms to be undercoordinated (colored light blue). These surface Zr atoms are not fully coordinated with oxygen atoms, and thus oxygen atoms are able to move faster inward without being properly coordinated with Zr. Figure 3.15c shows a decrease in the number of Zr-S bonds as a function of time, where the decrease is fastest at the highest pressure of 0.1 GPa. Again, there is an abrupt change of the rate of chemical reactions at 0.16 ns for the 0.1 GPa case, which is consistent with the observations in Figures 3.15a and 3.14b. 46 Local oxide structures: Figure 3.17 shows the interfacial structures at 0.6 ns under 0.1 GPa -- locally disordered structures with Zr atoms bonding to different number of O atoms at different distances from the oxidation front. Figures 3.18 and 3.19 compare Zr-O bond-length and O-Zr-O bond-angle distributions, respectively, at 0.6 ns under various pressures. For calculating O-Zr-O bond-angle distribution, we have used the Zr-O bond-length cutoff of 2 Å to be consistent with ZrO bonds in ZrO2 (142, 143). We found that the Zr-O peak under 0.1 GPa is broader than those under lower pressures (0.08 and 0.06 GPa), which substantiate the structural disordering due to high pressure mentioned in the main text. This is likely due to the over-coordinated Zr and more Zr-O bonds with decreased bond order in the transient amorphous oxide. The broad bond-angle distributions in also indicate disordered oxide structures in the early stage of oxidation. Figure 3. 17. Snapshot of interface structure at time 0.6 ns under pressure 0.1 GPa. The color scheme is the same as that in Figure 3.1 in the main text. Zr-O bonds are drawn to highlight the progress of oxidation. 47 Figure 3. 18. Zr-O bond-length distribution at 0.6 ns under pressures (a) 0.1 GPa, (b) 0.08 GPa, and (c) 0.06 GPa. Figure 3. 19. O-Zr-O bond-angle distribution at 0.6 ns under pressures (a) 0.1 GPa, (b) 0.08 GPa, and (c) 0.06 GPa. Figure 3. 20. Snapshots of ZrS2 oxidation under 0.1 GPa at times t = 0.3 ns (a), 0.4 ns (b), and 0.5 ns (c). The color scheme is the same as that in Figure 3.14. 48 To elucidate the different stages of oxide growth, Figure 3.20 shows the oxidation front at different times for the 0.1 GPa simulation. At time t = 0.2 ns, oxidation is confined in the first layer shown in Figure 3.14a. By 0.3 ns (Figure 3.20a), the oxidation front has reached the second layer, and by 0.4 ns it has reached the third layer (Figure 3.20b). The oxide surface is not flat at this time, showing kinks in the middle. At time t = 0.5 ns (Figure 3.20c), the oxide surface has become flatter again, due to oxidation of ZrS2 in the middle of the frame. This suggests that a layer-by-layer like oxidation process continues even during the second stage of faster oxidation, but the oxidation front now encompasses multiple vdW layers. As a result, the average oxide depth increases continuously (Figure 3.14a) instead of stepwise despite the underlying kink-mediated layer-bylayer oxidation. A similar layer-by-layer oxidation process was reported in Cu2O nano-island growth (144) and twin-boundary-assisted oxidation in Ag nanocrystals (145). Diffusive Oxygen Transport: Figure 3. 21. Mean square displacement of oxygen atoms in the amorphous oxide after the onset of continuous oxidation at time 0.16 ns under pressure 0.1 GPa. 49 Figure 3. 22. Oxygen concentration profile along the slab-normal direction, z, at time 0.47 ns under pressure 0.1 GPa. Despite the high reactivity at the oxide front under highest pressure (0.1 GPa), oxygen is still transported by diffusion through continuous amorphous oxide and driven by oxygen-concentration gradient. According to Fick’s law, oxygen flux through the oxide layer is O = −/, where D is the diffusion coefficient, c is the oxygen concentration, and z is the coordinate along the slabnormal direction. To demonstrate the diffusive behavior, we have plotted the mean square displacement (MSD) of oxygen atoms in the amorphous oxide after the onset of continuous oxide growth (0.16 ns) as shown in Figure 3.21. The MSD exhibits a diffusive behavior, from which the diffusion coefficient is estimated as = 〈|∆⃗| 2 〉 6 = 4.1 × 10−4 (cm2 ⁄ /s), which is consistent with the fitted B parameter in the Deal-Grove model in the main text. The other ingredient of Fick’s law is oxygen-concentration gradient, /. Figure 3.22 shows the oxygen-concentration profile in the oxide layer between the O2 gas/oxide interface (left) and oxide/ZrS2 interface (right) at time 0.47 ns. The oxygen concentration was calculated by dividing the oxide layer into cells of side 4.5 Å. The figure indeed exhibits linear drop of the oxygen concentration in the oxide layer. Together, Figures 3.21 and 3.22 support the applicability of the Deal-Grove model to the continuous oxidegrowth stage as stated in the main text. Finally, it is worth noting that even the highest pressure of 0.1 GPa is orders-of-magnitude smaller than the elastic constant of ZrS2 in the surface normal 50 direction, 33~30 GPa (146). As a result, we have not observed mechanical narrowing of the van der Waals (vdW) gap. Number of oxygen atoms in the oxide is counted through the cell-based surface construction method. Table 3.2 shows the pseudo-code for this method. Table 3. 2. Pseudo-code for cell-based surface construction method % Lx, Ly, Lz are the box lengths in x, y, z axis. % lx, ly, lz are the cell lengths in x, y, z axis. % Px, Py, Pz are the number of cells in x, y, z axis. % NP is the total number of the cells, NP = Px×Py×Pz. % Cell stores the flag of cells, e.g. Cell[p] == 1 represents that cell with ID p has Zr atoms. % position.x, position.y, position.z are the position of specified atom. % (px, py, pz) is the cell vector ID of specified atom, and p is the cell sequential ID of the atom. % Total_count is the counting for the number of O atoms in oxide. % Initialize all the cells to be unoccupied — 0 is unoccupied Cell = np.zeros(NP) %Reset all the cells with Zr atoms to be occupied — 1 is occupied For all Zr atoms #Find cell id, p of this Zr atom: px = math.ceil(position.x/lx)-1; py = math.ceil(position.y/ly)-1; pz = math.ceil(position.z/lz)-1; p = px×Py×Pz + py×Pz + pz; Cell[p] = 1; % Scan all the O particles to construct the occupancy table Total_count = 0; For all O atoms #Find cell id, p of this O atom: px = math.ceil(position.x/lx)-1; py = math.ceil(position.y/ly)-1; pz = math.ceil(position.z/lz)-1; p = px×Py×Pz + py×Pz + pz; if Cell[p] == 1: Total_count += 1; Applicability of the Deal-Grove model: Oxide growth behavior after the sudden onset of fast oxidation is characterized by gradually decreasing oxidization rate similar to the Deal-Grove model of silicon thermal oxidation, i.e., linear growth controlled by interfacial chemical reaction followed by parabolic growth indicating a diffusion-controlled process (118), as shown in Figures 51 3.20 and 3.21. Oxidation of ZrS2 in the parabolic-growth stage has recently been proposed to be diffusion-controlled (89, 90). Figure 3.15a exhibits both parabolic (at 0.1 GPa) and linear (at 0.08 GPa) oxide growth in the second oxidation stage. Such linear-to-parabolic oxidation is conventionally described by the Deal-Grove model (118). The model assumes that oxidation proceeds by inward movement of oxidant species rather than by outward movement of reactants. The process applies after an initial transient period with the consequence that the fluxes of oxidant are identical at all time steps. Key assumptions for the Deal-Grove model are: (1) Oxidants are transported from the oxidizing gas to the outer surface, where they react or are adsorbed; (2) they are then transported across the oxide film towards the bulk; and (3) they finally react at the oxidation front to form a new layer of oxide. In applying these assumptions to the oxidation of ZrS2, we note differences that require examination. On one hand, it is important to clearly identify the reaction products, i.e., whether there are sulfur atoms remaining after oxidation, and when and how SO2 molecules are emitted if the final products include SO2. During our simulation timespan (< 0.6 ns), we did not observe SO2 formation. Accordingly, we do not consider the formation of gas products and their out-diffusion here. We can then assume that the oxidation process, A(solid) + B(gas) = C(solid), follows a first-order reaction during the simulation time window. On the other hand, there is a spatial gap between consecutive ZrS2 layers unlike bulk silicon, and thus the deviation from the Deal-Grove model resulting from the gap needs to be considered. In our simulation, the lattice constant for the surface-normal direction is 5.85 Å and the gap between vdW layers is 2.925 Å, thus an effective diffusion coefficient incorporating the effect brought by the gap should be considered. 52 Figure 3. 23. Fitting of the oxide growth depth vs. time for pressures (a) 0.1 GPa, (b) 0.08 GPa and (c) 0.06 GPa. Simulation results are shown blue, red and black, respectively, for 0.1, 0.08 and 0.06 GPa, while fits are shown as gray dashed lines. All fittings are with 95% confidence bound. The three stages of oxidation (I, IIa and IIb) are indicated by magenta, green and cyan backgrounds, respectively. In Figure 3.23, we present fits of the Deal-Grove model to our simulation results for pressures 0.08 and 0.1 GPa. The Deal-Grove equation describes the oxide growth behavior at both linear and parabolic growth stages: 0 2 + 0 = ( + ), (3.11) where 0 , and  represent oxide depth, oxidation time, and shift in the time to correct for the presence of the initial oxide layer. In our case,  can be identified as the onset time of the fast continuous oxidation in the second stage, after the slow layer-by-layer oxidation in the first stage shown in Figure 3.15a. We thus fit the data points using the Deal-Grove model only beyond the 53 oxide depth of 0 ( = 0.3 ns) = 0.76724 Å for 0.08 GPa and 0 ( = 0.16 ns) = 0.90787 Å for 0.1 GPa. At the lowest pressure of 0.06 GPa, the finite simulation time only places a lower bound of 0.6 ns for  (Figure 3.23c), before which Eq. (3.11) does not apply. Oxidation for > is wellmodeled by a linear law at 0.08 GPa (Figure 3.23b), while it exhibits both linear and parabolic behaviors at 0.1 GPa (Figure 3.23a). For long oxidation time ( ≫  and ≫ c = 2⁄4), Eq. (3.11) describes a parabolic growth law, 0 2 ≅ . The parabolic rate constant, , is proportional to the partial pressure of the oxidizing species in the gas phase according to Henry’s law (118). Accordingly, the crossover time, c = 2⁄4, which separates the linear ( < c ) and parabolic ( > c ) growth regimes (118), is inversely proportional to . At the highest pressure of 0.1 GPa, the estimated c is 0.25 ns and the fittings in Figure 3.23 exhibits a linear-to-parabolic crossover, which is expected from Eq. (3.11). On the other hand, only the first linear regime is observed within the simulated time at a lower pressure of 0.08 GPa, indicating a much longer crossover time, c > 0.6 ns. The pressure dependence of the crossover time is thus consistent with the Deal-Grove model. The short-time behavior of Eq. (3.11) for ≪  follows a linear law, 0 ≅ (⁄)( + ), which is dictated by interfacial reactions (118). From the linear fits in Figure 3.23, the estimated ⁄ = 87.3 and 361 cm/s, respectively at p = 0.08 and 0.1 GPa, which is again consistent with the pressure dependence expected from the Deal-Grove model. The parabolic regime is only probed for p = 0.1 GPa in our simulation data, which yields = 3.42 × 10−6 cm2 /s. While larger than typical diffusion constants in ambient thermal oxidation experiments, this value is an ordersof-magnitude smaller than those in high-pressure oxidation simulation of Al (147). Overall, the Deal-Grove model provides an excellent platform to mechanistically understand essential pressure control of the three stages of oxidation observed in our simulations: 54 • Stage I ( < ): Discrete layer-by-layer growth. • Stage IIa ( < < c ): Short-time continuous growth following linear law, which is dictated by interfacial reactions. • Stage IIb (c < ): Long-time continuous growth following parabolic law, which is dictated by oxidant diffusion. The two crossover times,  and c , are both pressure-dependent and thus are controllable by oxygen partial pressure. Our three simulations selectively expose these three stages of oxide growth within a given time window: Only stage I at 0.06 GPa; stages I + IIa at 0.08 GPa; and stages I + IIa + IIb at 0.1 GPa. In particular, layer-by-layer oxidation stage can be prolonged at lower pressures, which may be useful for device processing steps that require highly controlled and ultra-thin oxide formation. Initial oxidation stage: We study pressure-dependence of the initial oxidation process before time  , by computing the change of Gibbs free energy ∆ upon oxygen-atom adsorption on the ZrS2(001) surface. In Figure 3.12, we present ∆ for oxygen-atom adsorption on Zr sites as a function of oxygen chemical potential ∆ . The bars under the plot show the corresponding pressure at temperatures T = 300, 600 and 900 K. The calculated ∆ is negative, even below 1 atm, which suggests that initial oxidation is driven by strongly exothermic Zr-O binding. 55 Figure 3. 24. Snapshots of RMD simulation trajectory for ZrS2 oxidation under 0.1 GPa. Two examples in (a) and (b) show the transport of O atoms (labeled O1 in purple color) assisted by the Zr atoms (labeled Zr1 in blue color) bonded to inner layer. After the initial period of oxygen adsorption and layer-by-layer oxidation, we observe a process of more rapid oxygen diffusion and an amorphous oxide growth front that spans multiple vdW layers. We find that oxygen diffusion is intimately related to Zr-site disorder, as suggested by our previous RMD simulations (90). In Figure 3.24, we illustrate how disordered Zr atoms assist oxygen motion between vdW layers, using two instances drawn from the simulation at 0.1 GPa. In each case, a Zr atom (labeled Zr1) at the oxidation front has moved into a vdW gap, where it bonds with S atoms in adjacent layers. We observe that the O atom (labeled O1) moves across the vdW gap by transient bonding with Zr1, thus facilitating oxygen diffusion into the next ZrS2 layer. 56 We speculate that this mechanism is facilitated by the fact that sulfur is a soft ion that can adopt multiple oxidation states, and therefore can be readily oxidized and reduced by virtue of the motion of hard ions O2- and Zr4+. The formation of kinks in the oxide growth front (Figure 3.20) may result from the local acceleration of oxygen diffusion and subsequent oxide formation by Zr defects. This amorphous-oxide-mediated oxidation, in turn, explains the rapid increase in the oxide depth and the number of Zr-O bonds in Figure 3.15. Effects of Ensemble: To examine the effects of ensemble on the oxidation behavior, we have performed additional reactive molecular dynamics (RMD) simulations in the isothermal-isobaric (NPT) ensemble instead of those in the canonical (NVT) ensemble used in the main text. The NPT simulations exhibit the same sequence of oxidation events as in the original NVT simulations, including the crucial gap-closing mechanism. Figure 3.25 compares snapshots of NVT and NPT simulations at a pressure of 0.1 GPa, showing similar gap-closing events. Figure 3. 25. Snapshots of O transport through gap closing for ZrS2 oxidation under 0.1 GPa in (a) NVT and (b) NPT ensembles. 3.4.5 Critical Oxidation Moments 57 Method: In order to examine important oxidation events, we used the VASP software (122, 131) to run AIMD simulations for the ZrS2 oxidation. We performed DFT relaxation using the projector augmented wave (PAW) method (132) and Perdew-Burke-Ernzerhof (PBE) functional (55), with 500 eV plane wave cutoff. Brillouin zone was sampled over 1 × 1 × 1 Monkhorst-Pack k-point mesh (133). The system dimension is 12.644× 14.600 × 26.523 Å 3 , containing a Zr48S96 slab immersed in 72 oxygen molecules. We perform AIMD simulations at 1,200 K in NVT ensemble. The system temperature was controlled using Nosé-Hoover thermostat (138, 139). Periodic boundary conditions were applied in all directions. Oxidation events: 58 Figure 3. 26. Substitution of a S atom by an O atom. (a) Snapshots of an oxygen atom replacing a sulfur atom. (b) Schematic of the substitution. Figure 3.26 shows the substitution of a S atom by an O atom observed in the simulations. At 0.200 ps, one oxygen molecule adsorbs to the ZrS2 surface, the oxygen atom labeled ‘O1’ bonds to the sulfur atom labeled ‘S1’. With the subsequent lattice vibration, the zirconium atom labeled ‘Zr1’ 59 breaks one of its original Zr-S bonds and forms a new bond with O1 at 0.246 ps. The S1-Zr1, and S1-Zr3 bonds break at 0.300 ps. The S1-O1 bond and Zr1-O1 bond are rotating and at 0.321 ps O1-Zr3 bond forms. At 0.350 ps, S1-O2 bond forms. The S1-O1, Zr1-O1, and Zr3-O1 bonds keep rotating and at 0.420 ps, Zr2-O1 bond forms. At 0.450 ps, S1-Zr2 bond breaks, and the original Zr1-S1, Zr2-S1, and Zr3-S1 bonds become Zr1-O1, Zr2-O1, and Zr3-O1 bonds. In Figure 3.26b, we draw the schematic mechanism to illustrate the oxidation process. Figure 3. 27. Movement of a zirconium atom toward the inner layer. Figure 3.27 shows the process of a zirconium atom bonding to a sulfur atom of its inner layer. At 1.599 ps, the zirconium atom labeled ‘Zr1’ bonds to two oxygen atoms, labeled ‘O1’ and ‘O2’, together with three sulfur atoms, labeled ‘S5’, ‘S2’ and ‘S3’. At 1.804 ps, Zr1-O2 bond and Zr1- S5 bonds break, meanwhile Zr1 bonds to sulfur atom labeled ‘S1’. The movement of Zr1 could be 60 explained by the fact that the coordination number of Zr1 is decreased from 5 to 4, together with the bond formation with S1. At 1.849 ps, Zr1-S3 bond breaks and the coordination number of Zr1 is decreased to 3. At 1.879 ps, Zr1-S1 bond breaks and the distance between Zr1 to the next layer is short enough so that Zr1 bonds to the sulfur atom labeled ‘S4’ in the next layer. It should be noted that the sulfur atom S1 assists the movement of Zr1 from original layer toward the inner layer. With the formation of the Zr1-S4 bond, Zr1 is connecting two adjacent layers and closing the vdW gap, which further leads to the local disorder. This is in consistent with the previously reported vdW gap closing mechanism as shown in RMD simulations, which plays a key role for continuous oxidation at larger scale and longer time. 61 Figure 3. 28. Oxygen transfer in the ZrS2 inner layer through bond switching and structure rearrangement. Figure 3.28 shows how oxygen transfers in the ZrS2 inner layer, and the O atom labeled ‘O1’ is the target in this process. In the snapshot at 2.151 ps, one oxygen molecule bonds to the same Zr atom labeled ‘Zr1’. This oxygen molecule consists of one oxygen atom labeled ‘O1’ and another labeled ‘O2’, where both O1 and O2 bond to Zr1 while O2 bonds to another Zr atom labeled ‘Zr2’. At 2.238 ps, O1 atom also bonds to Zr2 atom, and immediately O1 bonds to a sulfur atom labeled ‘S1’. Later O1-O2 bond breaks at 2.284 ps. At 2.316 ps, O1 bonds to another Zr atom labeled ‘Zr3’ and O1-Zr2 bond breaks at 2.341 ps. After lattice vibration, S1 bonds to another Zr atom labeled 62 ‘Zr4’, thus forming a twisted O1-S1-Zr4-S2-Zr1 ring at 2.410 ps. It is noted that this five-member ring is in the ‘envelope’ shape, later this five-member ring rearranged its position to a flat shape as shown at 2.551 ps. With the bond-switching similar to Si-O system (128) and the rearranging for the shape of five-member ring, the oxygen atom O1 is transferred in the inner layer. Figure 3. 29. Charge redistribution and bonding characteristics during O1-S1-Zr4-S2-Zr1 ring rearrangement. Mulliken charges (a) and COHP (b) analysis for t = 2.341 ps. Mulliken charges (c) and COHP (d) analysis for t = 2.551 ps. The charge behaviors and bonding characteristics between the atoms in this O1-S1-Zr4-S2-Zr1 five-member ring are also studied. We computed the Mulliken charges and Crystal Orbital Hamilton Population (COHP) for time before the ring rearrangement at 2.341 ps and after the ring rearrangement at 2.551 ps. As shown in Figure 3.29, charge redistribution takes place during the atoms rearranging their location. Apart from the charge redistribution, we also observe a transition from anti-bonding state to non-bonding state for O1-S1 bond in this five-member ring during rearrangement. The COHP analysis results show strong anti-bonding state for O1-S1 bond before the ring rearrangement, which corresponds to the repulsive force between O1 and S1 atoms. The ring rearrangement brings the non-bonding state for O1-S1 bond, which is more stable than before rearrangements. 63 Figure 3. 30. Oxygen transfer in the ZrS2 inner layer through Zr-O bond rotation. In Figure 3.30, we found that when one Zr atom is moving toward the next layer, the spaces between it with the surrounding Zr and S atoms in the original layer increases. The oxygen atoms adsorbed on it will move through the space toward the next layer meanwhile keeping the Zr-O bond. The oxygen atoms nearby will also be affected thus move toward the next layer via the spaces though not directly bonded with the Zr atom. 3.5 Summary To study the ZrS2 oxidation, a first-principles based ReaxFF force field was first developed and optimized using multi-objective genetic algorithm in EZFF force field training software. A clientserver model was used for training this force field, where the main node performed MOGA training, the client nodes performed RMD simulations based on different force field. The optimized force field gives good agreements with ground truth calculation on atomic charges and bond-population dynamics. Using the optimized ReaxFF force field, I performed RMD simulations to study the oxidation mechanism of ZrS2. Our reactive molecular dynamics simulations supported by first- 64 principles calculations elucidate the atomistic oxidation mechanisms and our work suggests that oxygen partial pressure is indeed a promising control parameter for such active oxidation for future 2D electronics. We compared the oxidation rate between MoS2 and ZrS2, and found that ZrS2 has a much higher oxidation rate. We also observed an anisotropic oxidation rate, where ZrS2 (210) surfaces have a higher oxidation rate than on (001) surfaces. Furthermore, the oxidation process is found to be diffusion controlled. First-principles calculations show that the initial adsorption of oxygen atoms on the surface is driven by the large binding energy. We find that pressure can selectively expose three stages of oxide growth in the Deal-Grove model within a given time window: (1) initial incubation, (2) reaction-limited linear growth, and (3) diffusion-limited parabolic growth, with increased pressure. We also find a trade-off between the growth speed and the quality of the grown oxide structure, as well as pressure control of the morphology of semiconductor/oxide interfaces. After an initial stage of slow, layer-by-layer oxidation, local breakdown of van der Waals gaps accelerates the diffusion of oxygen, creating a kink-mediated amorphous oxide growth front, with kinetics well described by the conventional Deal-Grove model. The transition time from layer-by-layer to continuous oxidation, as well as that between reactionlimited linear and diffusion-limited parabolic oxidation within the Deal-Grove regime, is well controlled by pressure. Such atomistic elaboration in the well-established Deal-Grove framework is expected to apply to broad TMDC materials, thus rationally guiding experimental growth of vdW heterostructure devices. The ab initio molecular dynamics simulations reveal detailed picture for the oxygen substituting sulfur atoms and oxygen atoms transporting in the bulk. At even longer time scales (beyond our simulation window), we expect that the kinetics may deviate from the Deal-Grove model due to the out-diffusion of the oxidation byproducts, which is a topic of future study. Simulations performed at low temperatures does not extend for long enough 65 time to observe the final products, though simulations performed at a temperature of 1,500 K do observe the formation and out-diffusion of sulfur atoms coming out accompanied with oxygen atoms during the oxidation process. Much longer simulations, likely using accelerate molecular dynamics (AMD) will be required to encompass the entire oxidation stages for lower pressure cases. 66 Chapter 4: Photoexcitation-Induced GeTe Local Disorder This chapter is based on the published work: Yang, L.; Tiwari, S. C.; Fukushima, S.; Shimojo, F.; Kalia R. K.; Nakano, A.; Vashishta, P.; Branicio. P. S., Photoexcitation-Induced Nonthermal Ultrafast Loss of Long-Range Order in GeTe. J. Phys. Chem. Lett. 2022, 13, 43, 10230–10236 (148). 4.1 Introduction The (GeTe)x(Sb2Te3)1-x alloy family shows nanosecond fast crystalline-amorphous phase transition and is the most promising and studied phase change materials (26, 149-152). The amorphous phase of PCMs is conventionally obtained by melting the crystalline materials, by application of short laser or electric pulses, followed by quenching to room temperature. Alternatively, it has been reported that femtosecond laser photoexcitation can induce a nonthermal amorphization process (153-155). Using first-principles MD, Li et al. (155) showed that a 9% electronic excitation leads to amorphization of GST within several picoseconds at a temperature below 700 K. Tiwari et al. (155) demonstrated ultrafast nonthermal amorphization in Sb2Te3 induced by photoexcitation. Experimental reports also confirm ultrafast photoexcitation-driven phase-change in GST alloys (156). As a prototype binary PCM in the GST alloy family (27, 157), GeTe is known to adopt a bond hierarchy at low temperatures, where short and long Ge-Te bonds indicate strong and weak bonding in a Peierls distorted rock salt structure (158). Kolobov et al. (153) demonstrated that 67 small atomic misalignment could trigger the destruction of the subsystem of weaker bonds in GeTe and subsequent collapse of the ordered phase leading to nonthermal amorphization. In contrast with the crystalline phase, which is characterized by the presence of resonance bonding (30, 31, 159), the amorphous phase in PCMs displays primarily covalent bonding. In GeTe, while some Ge atoms preserve a local structure similar to that of the crystalline phase (160-162), which is referred to as “defective octahedral sites” or “pyramidal sites”, most Ge atoms acquire a tetrahedral bonding geometry, characterized by four equivalent Ge-Te bonds forming 109.5° bond angles in a typical sp3 hybridization (163). This raises the question about the formation mechanism of this tetrahedral geometry and the co-existence of octahedral- and tetrahedral-like sites in the amorphous phase (29, 34, 164, 165). The effect of laser pulses on a solid may not be described simply as a “heating effect” when the induced repulsive interatomic forces swiftly induce lattice disorder. That is the case when the laser pulse is intense, and the duration is as short as femtoseconds, which can generate effects in a time scale significantly below the characteristic phonon relaxation time (166). Nonetheless, the photoexcitation-induced nonthermal bond changes in GeTe leading to amorphization are not fully understood. Possible mechanisms, such as the phase flip of the wave function (166), and the formation of transient three-center bonds (167) have been proposed and need to be further validated. In addition, no direct investigation of the evolution of the GeTe structure with excited states leading to the loss of long-range order induced by photoexcitation has been performed so far to validate these views. Here, we apply nonadiabatic quantum molecular dynamics (NAQMD) simulations to investigate the evolution of GeTe structure with photoexcited states. NAQMD simulations show disorder in 68 GeTe via a nonthermal path within a few picoseconds when the photoexcitation level exceeds 4.0% of the valence electrons. The simulation results provide insights into the onset of the local disorder and the change in the GeTe bonding nature, preceding the nonthermal ultrafast loss of long-range crystalline order. 4.2 Method NAQMD simulations were performed on a supercell composed of 3×3×3 cubic GeTe crystalline unit cells with lattice constant a=b=c=6 Å, containing a total of 216 atoms. The system was first relaxed in a canonical ensemble at a temperature of T = 300 K using a Nosé–Hoover thermostat (168). The gamma point was used to sample the Brillouin zone. NAQMD simulations were then performed in a canonical ensemble at T = 700 K. A time step of 1.2 fs was used in the integration of the equations of motion. All simulations were performed using the plane-wave basis quantum molecular dynamics code QXMD (62, 63). In quantum molecular dynamics, the trajectories of all atoms are integrated based on interatomic forces calculated using the Hellmann-Feynman theorem in the density functional theory framework (54, 169-171). The GGA-PBE approximation is employed for the exchange-correlation energy, and dispersion forces are calculated with the DFT-D correction method (55, 172). The PAW method is used in the calculation of the electronic states (132). The plane-wave cutoff energy was set at 30 Ry for the electronic wave functions and 250 Ry for the electronic charge density. The dynamics of the excited states is modeled using the fewest-switches-surface-hopping method (63, 173, 174). In NAQMD, excited-state forces are computed using an extension of the Harris−Foulkes approach (63). Details of the NAQMD method are reported elsewhere (63). 4.3 Results 69 4.3.1 Structural Changes Mechanism Figure 4. 1. Changes in the structure of GeTe with photoexcitation. (a) GeTe structure after thermalization at 300 K. (b-e) GeTe structures after 5 ps of photoexcitation at different excitation levels n: n = 2.6% (b), 4.0% (c), 5.2% (d) and 7.5% (e). Green and yellow spheres represent Ge and Te atoms, respectively. In the GeTe cubic crystalline phase, each atom has a coordination number of six and displays an octahedral atomic arrangement typical of rocksalt structures with cation/anion sites occupied by Ge/Te atoms, respectively (175). Figure 4.1a shows the thermalized GeTe structure at 300 K, along with GeTe structures at 700 K after five ps of photoexcitation at different levels of excitation. NAQMD simulations are performed with n = 2.6%, 4.0%, 5.2%, and 7.5% valence electrons excited. The GeTe structure relaxed at 300 K shows a broader Ge-Te bond length peak, as shown in Figure 4.2, which can be explained by Peierls distortions (PD) (158, 164, 176-179). The distortion occurs along the [111] direction and was reported to lower the electronic energy (44 meV/atom) and increase the band gap (180) from 0.37 eV to 0.5 eV. In our simulations at 300 K, the thermal vibration destroys the well-defined arrangement of long and short bond lengths, and instead, the bond length distribution displays a broad distribution around ~3 Å (Figure 4.2). As shown in Figure 4.1b, excitation at n = 2.6% is insufficient to induce structural disorder. However, for photoexcitation at n = 4.0%, the results show local disorder, indicating that the threshold to trigger GeTe local disorder is achieved. A higher excitation level leads to extensive disorder and 70 widespread loss of the long-range crystalline order and the formation of Ge-Ge and Te-Te wrong bonds. Figure 4. 2. Change in Ge-Te bond length distribution after thermalization at 300 K. (a) Ge-Te bond length distribution in the perfect cubic structure at 0 K. (b) Ge-Te bond length distribution in the GeTe crystal thermalized at 300 K. As shown in Figure 4.2, the Ge-Te bond length for the cubic structure is 3.0 Å. After thermalization at 300 K, the bond length distribution single peak is split showing longer and shorter bonds in range between 2.6 Å and 3.5 Å. The split is caused by distortions of the structure occurring during the thermalization process and are consistent with the reported distorted rock-salt structure (158), where the Peierls distortion generates three short and three longer bonds, located at ~2.8 Å and ~3.2 Å. The thermal vibration at 300 K broadens the bond length distribution, which shows instead a single broad peak. 71 Figure 4. 3. Photoexcitation-induced structural change in GeTe. (a) Mean square displacement of Ge as a function of time; (b) Evolution of the fraction of wrong bonds. The black curve corresponds to adiabatic MD simulation, n = 0.0%. Cyan, red, blue, and green curves correspond to excitation of n = 2.6, 4.0, 5.2, and 7.5%, respectively. To characterize the structural change illustrated in Figure 4.1, the MSD of Ge atoms is calculated as a function of time (Figure 4.3a), while the evolution of the fraction of Ge-Ge and Te-Te “wrong” bonds is plotted in Figure 4.3b. The induced atomic diffusion leads to widespread bond breaking and loss of long-range order. It should be noted that such an amorphization takes place without any melting as the temperature of the system is kept at 700 K, well below the melting point (at ~1,000 K). That implies that by applying femtosecond laser photoexcitation, the crystalline to amorphous transition time could be reduced from hundreds to a few picoseconds using a nonthermal path. 72 The black curves in both plots represent the simulation with no excited electrons, i.e., n = 0.0%, where there exists neither atomic diffusion nor wrong bonds. The cyan line in Figure 4.3a, corresponding to n = 2.6%, overlaps the reference MSD black line indicating the preservation of the crystalline structure. The cyan line in Figure 4.3b shows negligible wrong bonds at n = 2.6%, in agreement with the MSD results. On the other hand, results for n = 4.0% (red curves) show a slightly higher MSD and a larger fraction of generated wrong bonds than that of n = 0.0% and 2.6%, which is consistent with the picosecond photoexcitation-induced localized disorder observed in Figure 4.1c. One can note a decrease in the MSD from time t = 2.4 ps to t = 4.8 ps, indicating the recovering of the crystalline structure, which is consistent with the decrease in the fraction of wrong bonds. For n = 5.2% (blue curves), the MSD increases to about 1.2 Å2 then remains constant, indicating that large atomic displacement can be induced by electronic excitation at n = 5.2% or higher. Correspondingly, the fraction of wrong bonds at n = 5.2% increases sharply in the first 2 ps and then remains constant at a relatively high value of ~0.35. For n = 7.5% (green curves), the MSD increases steadily during the simulation timespan. Correspondingly, the fraction of wrong bonds also increases during the whole simulation reaching almost 80%. 73 Figure 4. 4. Local disorder triggered by Ge diffusion at n = 4.0%. (a-c) Shift of Ge atom from octahedral to tetrahedral sites for t = 1.833 ps (a), 2.094 ps (b), and 3.942 ps (c). (d-f) Schematic of the corresponding local disorder mechanism. The results in Figures 4.1 and 4.3 indicate a threshold for the photoexcitation-induced disorder at n = 4%. It is instructive to investigate the associated disorder mechanism activated at this excitation level. Figure 4.4 shows the illustrations of the evolution of the structure at the local disorder site. The local structure change indicates the shift in a Ge atom position from an octahedral to a tetrahedral site in a nearest neighbor unit cell. The transition is identified as the trigger point for the local disorder, which is also found elsewhere (154, 166). From Figures 4.4a-c, at t = 1.833 ps, a Ge atom (labeled Ge1) is octahedrally bonded with six neighboring Te atoms. At t = 2.094 ps, as the bond lengths between Ge1-Te5 and Ge1-Te6 atoms increase, the bonds become weaker and finally rupture. At t = 3.942 ps, shown in Figure 4.4c, the Ge1-Te7 bond is also broken while a new Ge1-Te4 bond is formed, leading to a tetrahedrally coordinated Ge. While not shown in Figure 4.4, at longer times, the formation of an intermittent Ge1-Ge2 bond is also observed as the Ge1 atom moves along the [111] direction towards the Ge2 atom position. As shown in Figures 4.4d-f, 74 the schematic mechanism for this Ge octahedral to tetrahedral coordination transition is proposed: i) Bond breaking with two nearest neighbor Te atoms; ii) motion along the [111] direction; and iii) bond formation with diagonal Te atom with simultaneous bond breaking with another Te atom. The phase-flip mechanism proposed by Kolobov et al. (164) suggests that the wave-function phase could be changed at high-level excitation, resulting in the weakening of GeTe longer bonds, which would transition from bonding to antibonding states. The region with concentrated excited orbitals has a high possibility of experiencing such phase-flip explaining the repulsive interatomic forces leading to the observed disorder. This process is consistent with the observed increase in the MSD shown in Figure 4.3b. At the transition threshold, between 4% and 5.2% of valence electrons excited, the MSD values increase until the Lindemann criterion for melting is fulfilled, inducing localized nonthermal melting of the crystal structure. 4.3.2 Electronic Structure Calculations Figure 4. 5. GeTe density of states (DOS) at different simulation times for photoexcitation level n = 4.0%. (a) Total DOS at different times after photoexcitation. The black, red, blue, and green curves represent GeTe before photoexcitation (t = -0.1 ps) 75 and after photoexcitation at t = 0.01, 0.1, and 2.5 ps, respectively. (b) Partial DOS of cubic GeTe. (c) Partial DOS of the relaxed structure. (d-f) Partial DOS at excitation level n = 4.0% after t = 0.01 ps (d), 1 ps (e), and 2.5 ps (f). The blue, orange, yellow, purple, green, and cyan curves in Figures 4.5(b-f) represent Ge s, Ge p, Ge d, Te s, Te p, and Te d, respectively. Figure 4. 6. Atom and orbital projected band structure of cubic GeTe (̅). (a) Ge s, (b) Ge p, (c) Te s, (d) Te p. The colors indicate the weight of the projections. The Fermi level energy is set at 0 eV. (e) The Brillouin zone paths used. Figure 4.5 shows the total and partial DOS for GeTe structure during the simulation. One significant point in Figure 4.5a is the disappearance of the band gap after photoexcitation. By examining Figures 4.5b and 4.5c, the Te p is dominant in the valence band near the Fermi energy and Ge p is dominant in the conduction band near the Fermi energy in both the cubic and relaxed 76 GeTe structure. Figure 4.6 shows explicitly the projections of the atomic orbitals over the band structure, which is consistent with Figures 4.5b and 4.5c. After photoexcitation, as shown in Figures 4.5d-f, we observe a crossing of the Te p to the conduction band and the Ge p to the valence band. We also observe the crossing of Ge s from conduction band to valence band. The behavior of these crossings suggest that photoexcitation plays an important role on the behavior of the Te p, Ge p, and Ge s electrons. The purple and green curves in Figures 4.5b-f highlight the energy difference between Te 5s and Te 5p states, which turns Te s-p hybridization unlike. Figure 4. 7. Average Mulliken charge as a function of time for different photoexcitation levels. (a-d) photoexcitation for n = 2.6% (a), 4.0% (b), 5.2% (c) and 7.5% (d). The blue and red curves represent the average Mulliken charge for Ge and Te atoms, respectively. To better understand the evolution of the electronic state after photoexcitation we calculate the average Mulliken charge on Ge and Te atoms during the excited state dynamics as shown in Figure 4.7. Negative times indicate states before photoexcitation while positive times indicate the evolution after the instantaneous photoexcitation. The charge transfer from Te to Ge occurs in all cases at time t = 0 ps, when the photoexcitation is applied. For n = 5.2% and 7.5%, shown in 77 Figures 4.7c, d, the average charge reverse sign, indicating a strong charge transfer from Te to Ge. This is followed by a swift partial recovery of the Ge and Te ground state charge within 5ps. Figure 4. 8. Mulliken bond overlap as a function of time for Ge-Ge, Ge-Te, and Te-Te bonds. The black curve corresponds to time t = -0.1 ps (before photoexcitation). Red, blue, and green curves represent results at t = 0.01, 1.0, and 2.5 ps (after photoexcitation). The overall evolution of the excited electronic state is provided by the calculated Mulliken bond overlap as a function of time. Results for excitation at n = 2.6% and 7.5% are shown in Figure 4.8. Before photoexcitation at time t = -0.1 ps, Ge-Ge and Te-Te bond pairs are in antibonding states with peaks at -0.2 and -0.08, respectively. Meanwhile, the Ge-Te bond pair is in a bonding state as expected, with a peak at 0.25. For 2.6% excitation, at all times, the Ge-Ge, Ge-Te, and Te-Te states are preserved with only small changes in the peak positions and widths. In sharp contrast, at the large photoexcitation level of n = 7.5%, at time t = 1 ps and 2.5 ps, the Ge-Ge and Te-Te bond pairs show bonding characteristics, which is consistent with Ge-Ge and Te-Te wrong bond formation. That is particularly the case for Ge-Ge, which shows a well-defined peak at 0.5. In the case of TeTe, a positive tail indicates the presence of a non-negligible bonding state. In contrast, Ge-Te bond 78 pairs, while still displaying a mostly positive bonding state, feature now a negative antibonding state tail at t = 1 ps and 2.5 ps, consistent with Ge-Te bond scission. The combination of Ge-Te bond scission with Ge-Ge, and some limited Te-Te bond formation is consistent with the picture of amorphization of the GeTe structure. 79 Figure 4. 9. Crystal orbital Hamiltonian population bonding analysis of Ge-Te and Ge-Ge interaction. (a-c) Bond length (a), COHP (b) and (c) before Ge1-Ge2 bond formation for t = -0.001 ps. (d-f) Bond length (d), COHP (e and f) after Ge-Ge bond formation for t = 2.122 ps. (g-i) Bond length (g), COHP (h and i) after Ge-Ge bond formation for t = 3.942 ps. Since the formation of Ge-Ge wrong bonds is observed in the disordering process, it is worth considering the bond nature in the octahedral-like to the tetrahedral-like transition of Ge coordination. We examine the homopolar Ge-Ge bonding based on crystal orbital Hamiltonian population (COHP) (181) analysis together with the electron localization function (ELF) (182). A relaxed GeTe local structure at 300 K is shown in Figure 4.9a, displaying the same atoms in Figure 4.4. The COHP for the average Ge-Te bond in the structure is shown in Figure 4.9b, displaying a stable bonding state. The flat COHP curve with 0 value at the Fermi level for the Ge1-Ge2 pair show a nonbonding state in the relaxed structure in Figures 4.9c. During the disorder process, the homopolar Ge1-Ge2 bond formed, as shown in Figure 4.9d, at t = 2.122 ps. We find that the GeTe bonds are at a slightly antibonding state; meanwhile, the Ge1-Ge2 pair displays a weaker antibonding state, as shown in Figures 4.9e, f. To further understand the effect of Ge-Ge bonds in the presence of a tetrahedral structure, we consider the tetrahedral structure found at t = 3.942 ps in Figure 4.9g. We find that the COHP for all four Ge-Te pairs is at a slight antibonding state in Figure 4.9h. Concurrently, Ge1-Ge2 (separated by 2.69 Å), Ge1-Ge3 (separated by 2.48 Å), and Ge1-Ge4 (separated by 2.86 Å) pairs show an antibonding state, as shown in Figure 4.9i. The evolution of the excited state with dynamic generation of bonding and antibonding states for GeTe and Ge-Ge pairs is one of the reasons for the swift bond-breaking after formation, as observed in our simulation. 80 Figure 4. 10. COHP bonding analysis of Ge1-Te5 and Ge1-Te6 interaction. (a-c) COHP bonding analysis for t = -0.001 ps (a), 1.833 ps (b) and 2.094 ps (c). The red and blue curves in (a-c) represent the interactions of Ge1-Te5 and Ge1-Te6, respectively. To understand the changes in Ge-Te chemical bonding associated with the local structural change shown in Figure 4.4, we turn to the analysis of COHP, which is shown in Figure 4.10. The COHP for the Ge1-Te5 and Ge1-Te6 pairs before photoexcitation (Figure 4.10a) indicate a stable configuration state, i.e., bonding states are shown to the right and antibonding states to the left region as a function of energy and the value of COHP is close to zero at the Fermi level. In contrast, after photoexcitation at time t = 1.833 ps, the COHP for Ge1-Te5 and Ge1-Te6 bonds indicate antibonding states, which drive the bond breaking process at t = 1.833 ps (Figure 4.10b). At time t = 2.094 ps, the Ge1-Te5 and Ge1-Te6 bonds are broken and the corresponding nearly flat COHP curves (Figure 4.10c) indicate non-bonding states. The photoexcitation-induced amorphization process described here does not involve melting and quenching. It has been reported that the destruction of the weaker (long) Ge-Te bonds could be obtained by distorting the perfect GeTe structure leading to a subsequent collapse of the long-range order resulting in amorphization (153). The destruction of the longer bonds is induced by photoexciting electrons from high-energy states in the valence band to low-energy states in the 81 conduction band, enhancing the concentration of non-equilibrium charge carriers. These charge carriers populate antibonding states of the weaker metavalent bonds (183), as shown in the COHP curves in Figure 4.9 and 4.10, leading to structural relaxation and bond rupture, causing the collapse of the long-range order. This mechanism effectively opens a path for lowering the energy barrier for bond breaking and triggering structural disorder at a temperature lower than the nominal melting point, similar to photoexcitation preferential bond rupture in As2S3 (184) and GST (185). Figure 4. 11. ELF of Ge-Te and Ge-Ge interaction. (a) ELF before Ge1-Ge2 bond formation for t = -0.001 ps. (b) ELF after GeGe bond formation for t = 2.122 ps. (c,d) ELF after Ge-Ge bond formation for t = 3.942 ps. The scale is from 0 (blue) to 1 (red). ELF = 1 represents perfect covalent bonds and any values between 0.5 and 1 reveal covalent bonds of various bonding strength, although ELF = 0.5 gives a metallic system. The corresponding relaxed GeTe local structure at 300 K is shown in Figure 4.9a. From the zero COHP value in Figure 4.9c and near 0.5 ELF value in Figure 4.11a between the Ge1 and Ge2 atoms, we can safely conclude the nonbonding states of the Ge-Ge atom pair in the relaxed structure. During the disorder process at t = 2.122 ps, the homopolar Ge1-Ge2 bond is formed as shown in Figure 4.9d, and the Ge1-Ge2 pair display a higher probability of electron localization as shown in Figure 4.11b. The tetrahedral structure at t = 3.942 ps is shown in Figure 4.9g. The ELF analysis shows a clear Ge1-Ge2, Ge1-Ge3, and Ge1-Ge4 bond formation in Figures 4.11c, d. 82 Figures 4.8 and 4.9 show that Ge-Ge bonding does not lead to strong antibonding states. Otherwise, Ge-Ge bonding would result in a large repulsive force, which would be unfavorable for the approximation and bond formation between Ge1 and Ge2 atoms. Furthermore, with proper Ge-Ge distance, the bonding state of Ge-Ge pairs is reported in the disordered region, supporting its stability. These results are supported by bond-weighted distribution function calculations (186), showing that Ge-Ge bonds are chemically stronger in tetrahedral motifs than in octahedral motifs. It is also suggested that the unexpected Ge-Ge bonding interactions support the reshuffling of electrons from antibonding Ge-Te into bonding Ge-Ge contacts, lowering the energy and enhancing the stability of the structure (187). Figure 4. 12. Evolution of Kohn-Sham energy eigenvalues and the excited orbitals for photoexcitation level n = 4.0%. (a) KohnSham energy eigenvalues as a function of time. The red, yellow, and blue curves represent orbitals that are empty, singly occupied, and doubly occupied, respectively. (b-d) Excited orbitals for t = 1.833 ps (b), 2.094 ps (c) and 3.942 ps (d). 83 The results in Figure 4.12 show the change in the Kohn-Sham energy eigenvalues and provide visualizations of the excited orbitals as a function of simulation time. After excitation at n = 4.0%, the appearance of the singly occupied orbitals, as shown in Figure 4.12a, is a consequence of the evolution of the excited state. The decline of the blue and yellow lines in Figure 4.12a shows the partial recovery of the ground state within the 4.5 ps of simulation. Figures 4.12b-d show the concentration of excited high energy orbitals in the locally disordered region. Figure 4. 13. Charge density difference (CDD) of different GeTe structures. (a, d) Bond lengths and CDD analysis for cubic GeTe; (b, e) for distorted rock-salt (rhombohedral) GeTe; (c, f) for structure after relaxation at 300 K. The Ge-Te bonds in (a-c) are colored according to their bond length, the CDD maps in (d-f) show the difference in the electron charge distribution as a result of bonding in the different structures. The results presented so far incite further discussions of the electronic structure of GeTe. The electronic configurations for Ge and Te atoms are [Ar]3d104s2px 1py 1pz 0 and [Kr]4d105s25px 2py 1pz 1 , 84 respectively. Though Te has 6 valence electrons, the energy difference between the 5s and 5p orbitals is so large that s-p hybridization can be neglected, which is similar to the case of Sb in GST (188). While Ge has an empty p-orbital, Te displays a p-orbital lone pair. It is natural to consider the formation of a (dative) two-electron bond between Ge and Te atoms, which allows for the preservation of the p-orbitals. The tetrahedral coordination involves sp3 hybridization of Ge atoms and the lowering of the energy makes this a stable configuration in the amorphous phase of GeTe (189). The octahedral coordination involves only p-orbital electrons within a resonant bonding configuration. i.e., the six p-orbitals from Ge and Te atoms form 6 resonant bonds per atom, where the lobe of the p-orbital is shared on both sides (30, 158). In recent years, the delocalized bonding present in crystalline GeTe and other phase change materials have been discussed on the novel framework of “metavalent” bonding, which was proposed to address bonding configurations that feature both metallic and covalent bonding characteristics (190). To investigate the bonding in the theoretical cubic GeTe, we consider its structure and the corresponding charge density difference (CDD), as shown in Figures 4.13a, d. The CDD confirms that the bonding in cubic GeTe displays delocalized charges piling up between Ge and Te atoms, not fitting the traditional two center two electrons Lewis covalent bonding, and also not fitting the homogeneous charge distribution around ions of a typical metallic bonding. Crystalline GeTe, due to its electronic conductivity, coordination number, optical dielectric constant, Born effective charges, and mode specific Grüneisen parameters, is classified as an incipient metal fitting the metavalent bonding view, where the same electrons responsible for bonding are rooted at the displayed vibrational, optical, and electronic properties (190, 191). It has been claimed that the destabilization of the centrosymmetric three-center bonds along with the loss of long-range order is the key to fast switching speed in PCMs (192). The Te-Ge-Te bond 85 in GST is found to be three-center four-electron (3c-4e), with charge piling up on both sides of the center Ge atom. In contrast, the bonding in rhombohedral GeTe, which is the low temperature stable structure of GeTe, is found to be three-center two-electron (3c-2e), with charge piling up along the shorter Ge-Te bond (192). The optimized rhombohedral GeTe structure, as shown in Figure 4.13b, displays the perfect arrangement of short and long bonds. The charge density difference plot for this optimized structure is shown in Figure 4.13e, which agrees with the 3c-2e bonding view previously proposed. In our MD simulations, since the GeTe structure is thermalized at 300 K before phtoexcitation, the thermal vibration randomizes the distribution of the Ge-Te short and long bonds. However, the set of short and long bonds and the charge piling up along short bonds can be clearly identified for most Ge atoms, as shown in Figures 4.13c, f. The results shown in Figure 4.13f are also consistent with the 3c-2e bonds model found in the rhombohedral structure. 4.4 Summary In this chapter, to study the photoexcitation-induced GeTe amorphization, we performed NAQMD simulations and DFT calculations. The NAQMD simulations show that nonthermal structural disorder in GeTe occurs after photoexcitation at an excitation level of no less than 4.0%. Accompanying the photoexcitation, we observe the swift charge transfer from Te to Ge atoms. The Mulliken bond overlap show increased Ge-Ge and Te-Te bonding interaction together with weakened Ge-Te bonding interaction. The examination of local disorder at 4.0% excitation level shows the trigger to be a Ge atom diffusion from octahedral to tetrahedral sites. Meanwhile, the excited orbitals are found to be spatially concentrated in the locally disordered region. This work 86 sheds light on structural and electronic mechanisms for photoexcitation induced GeTe amorphization and explores the bonding nature during this process. 87 Chapter 5: Surface Transfer Doping of MoO3-x on Hydrogenated Diamond This chapter is under submission. 5.1 Introduction With a wide bandgap of 5.5 eV, high carrier mobility (4,500 cm2 /V s for electrons and 3,800 cm2 /V s for holes) (193), large breakdown electric field (exceeding 10MV/cm), high thermal conductivity (22 Wcm-1K-1 ), and resistance to extreme environmental conditions, diamond-based field-effect transistors are promising candidates for high-frequency and high-power electronics (194). To control the conductivity and transform the well-known insulator, diamond, to be semiconductor has been given rising attention (195). However, compared to other semiconductors, it is hard to achieve traditional substitutional doping in diamond tightly packed lattice. Though p-type doping by boron and n-type doping in phosphorous or nitrogen are reported, the charge carrier densities are still low at room temperature (196, 197). Surface transfer doping has been proposed to be potential solution, where doping is achieved by electron exchange at the interface between the diamond surface and the dopants (195, 198, 199). MoO3 and V2O5 are reported to be two of the most promising oxide materials for this STD process, immensely increasing their application for FETs (200). Diamond surface are reported to show negative electron affinity after hydrogen termination, which facilitates the electron transfer from diamond to the surface electron acceptors and creates a quasitwo-dimensional sub-surface hole gas (2DHG) (201, 202). To achieve hydrogenation of diamond, 88 methods including plasma-enhanced CVD (PECVD), thermal hydrogenation (203, 204) etc. are reported. After hydrogenation, exposing diamond with surface dopants is easy, and the effect of STD will not disappear because of low bulk conductivity (198). MoO3 is an effective electron accepting/ hole injection material, which improves the performance and stability of STD (200). While hydrogenated diamond surface deposited with MoO3 was shown to exhibit high p-type surface conductivity, atomistic and electronic structures of the interface remain elusive. There is also a lack of investigation on the effects of oxygen vacancies, which is easily shown in MoO3 when exposed to atmosphere (205). Here, we performed first-principles-informed RMD simulations (206) to study the interfacial structure of diamond-MoO3-x. Reactive molecular dynamics simulations are performed using this force field to elucidate the deposition mechanism of MoO3-x on hydrogenated diamond (111) surface. Electronic density of states alignment and charge transfer at the interface are studied using first-principles calculations based on density functional theory for selected thermalized structures taken from the reactive molecular dynamics trajectories. We observe shifting of electronic density of states alignment and higher charge transfer for higher Mo oxidation state. Such atomistic and electronic details provide mechanistic understanding of the surface transfer doping process. 5.2 Method The system dimension is 15.018×17.830× 30.000 Å 3 . The numbers of oxygen atoms are 74, 82, and 93 for MoO2.3, MoO2.6, and MoO2.9, respectively. The numbers of carbon, hydrogen, and molybdenum atoms are the same for different systems, and there are 240 carbon atoms, 120 hydrogen atoms, and 32 molybdenum atoms. A three-step simulation workflow is used to optimize interfacial structure and compute electronic structure. First, we performed RMD simulations to 89 generate thermalized interfacial structure. Second, we used DFT to further optimize thermalized interfacial structure. Finally, we conducted DFT calculations to compute electronic structures. RMD to generate fully thermalized interfacial structure: MoO3-x was gradually heated to 3,300K and then thoroughly melted at 3,300K using RXMD software (207). We placed two ‘momentum mirrors’ when melting the oxides. One mirror was placed at 1 Å above the H-diamond surface to avoid chemical bonds formation at the interface, while the other was placed on top of the oxides to avoid the oxide atoms flying away. Since the purpose of this step is to melt and fully thermalize the oxides, the carbon atoms and hydrogen atoms are fixed during this process. After this, amorphous MoO3-x is deposited on the hydrogen terminated diamond (111) surface and relaxed at 10K. DFT to further optimize thermalized interfacial structure: We used DFT through VASP software (122, 131) to optimize this thermalized interfacial structure. Bottom 3 layers of atoms in the Hdiamond are kept fixed. We performed DFT relaxation using the PAW method (132) and PBE functional (55), with 450 eV plane wave cutoff. Brillouin zone was sampled over 1 × 1 × 1 Monkhorst-Pack k-point mesh (133). DFT to compute electronic structures: We performed DFT calculations to compute electronic structure of interface to study the electronic density of states and charge transfer. For electronic structure calculations, a fine 3 × 3 × 1 Monkhorst-Pack k-point meshes (133) is used. 5.3 Results 90 5.3.1 Interfacial Structure Figure 5.1 shows the side views of the structures of all systems. MoO3 has two thermodynamic phases, one being stable and the other being metastable. In this study, we use molten oxides for this deposition process. We also note that the MoO3-x|H-diamond interface changes with x varying. Figure 5. 1. MoO3-x encapsulated H-diamond systems with x being 0.7 (a), 0.4 (b), and 0.1 (c). The grey, red, blue, green balls represent C, H, Mo, and O atoms, respectively. Figure 5. 2. Mo-O pair distribution for MoO3-x encapsulated H-diamond systems with x being 0.7 (a), 0.4 (b), and 0.1 (c). Figure 5.2 shows the Mo-O partial pair distribution for all systems. The Mo-O partial pair distribution left shifts as the vacancy level, x, increases, which suggests that as more vacancies are 91 introduced into the material, the atoms or particles that remain tend to get closer to each other on average. Atomic arrangement is becoming more compact or denser as vacancy level increases. This could be due to the atoms or particles rearranging themselves to minimize the overall energy of the system in response to the introduction of vacancies. Figure 5. 3. Mo-O bond length distribution for MoO3-x encapsulated H-diamond systems with x being 0.7 (a), 0.4 (b), and 0.1 (c). Figure 5.3 shows the Mo-O bond length distribution for all systems, respectively. Along with the atomistic arrangement getting more compact, the Mo-O bond length is decreasing with the increase of the vacancy level, indicating that the Mo-O bonds get stronger as the vacancy level increases. 5.3.2 Charge Transfer To investigate the charge transfer process quantitatively, we first computed the Bader charges (208, 209) of atoms to calculate the charge difference after deposition. Figure 5.4 shows the Bader charges before and after depositing oxides on top of hydrogenated diamond (111) surface. Upon deposition, the net computed electron charge for hydrogenated diamond is 1.56, 1.80, and 2.0 for 92 x = 0.7, 0.5, and 0.1, respectively, with MoO3-x having equal and opposite charge demonstrating MoO3-x as an electron acceptor for different levels of oxygen vacancy. This illustrates the net charge transfer shows a higher charge transfer for higher Mo oxidation state. Before deposition we consistently find that hydrogen is more electronegative (having negative Bader charge), which has been attributed to changing from loosely bound bonding to tightly bound C-H bonding upon hydrogenation of diamond (210). Furthermore, the Bader analysis indicates the hole transfer is primarily localized to hydrogen, with the hydrogen atoms becoming more positively charged compared to pre-deposition and the carbon atoms becoming more negatively charged. As oxygen is strongly electronegative the weak attractive interaction between oxygen and hydrogen at the interface drives this charge transfer. Table 5. 1. Charge transfer per elemental atoms due to deposition Elements MoO2.3/H-diamond MoO2.6/H-diamond MoO2.9/H-diamond C -0.012 -0.004 -0.003 H 0.037 0.023 0.023 Mo -0.021 0.005 -0.013 O -0.012 -0.025 -0.016 93 Figure 5. 4. Average Bader charges of atoms before and after deposition. Dimond-MoO3-x systems with x being 0.7 (a), 0.4 (b), and 0.1 (c). The grey, red, blue, green balls represent C, H, Mo, and O atoms, respectively. 94 Figure 5. 5. Charge density profiles for MoO2.9 deposited hydrogenated diamond. (a) Side view of charge density difference of MoO2.9 deposited hydrogenated diamond. Yellow and brown regions represent electron and hole accumulation (electron depletion), respectively. (b) Side view of iso-surface of charge density at the interface with iso-level of -0.08 signifies the hole distribution, shown in dark brown. The color scheme for atoms is the same as that in Figure 5.1. To visualize the charge transfer process due to deposition, the charge density difference, ∆, is computed using MoO2.9 as an example. In this work, charge density difference, ∆, is defined as: ∆ = 3−/− − 3− − − , (5.1) where 3−/− is charge density of oxides deposited hydrogenated diamond surface, while 3− , − represent charge density of oxides, hydrogenated diamond, respectively, before deposition. From Figure 5.5, we can observe that electrons are accumulated in oxide while holes are accumulated in the hydrogenated diamond, especially around the surface region. We can thus safely conclude that due to deposition, electrons are extracted from diamond surface to the oxide, leading to a positively charged hydrogenated diamond and hole accumulation as expected in the surface transfer doping model (201). The isosurface of the hole density in Figure 5(b) observed at the interface shows the hole accumulation in H-diamond. The isosurface also signifies the spatially extended nature of doped holes in H-diamond, which is consistent with excellent transport properties of MoO3−/H − diamond interfaces (211). 95 Figure 5. 6. Density of states for pristine hydrogenated diamond and MoO3-x deposited hydrogenated diamond systems with x being 0.7, 0.4, and 0.1. EF represents the Fermi level. Gray dashed and dash-dotted lines indicate the valence-band top and Fermi energies, respectively. We also computed the element-projected density of states distribution to examine the surface transfer doping process. Figure 5.6 shows the PDOS for hydrogenated diamond before and after deposition of MoO3-x. From Figure 5.6, Fermi level shifts across the carbon valance band, mostly containing carbon p orbitals, due to deposition, which suggests p-type doping effect of oxides and metallic characteristic of doped diamond after deposition, which is consistent with theoretical and experimental reports (201, 212). Meanwhile the energy range of holes in hydrogenated diamond is demonstrated by the dash-dotted and dashed lines which show the Fermi energy and valenceband maximum (VBM) energy, respectively. This energy difference between the VBM and the Fermi level represents the energy range of doped holes, and is a decreasing function of the vacancy level, x. This illustrates that the oxygen vacancy limits the hole-doping capability of MoO3-x. To further investigate the electronic structure, we show the electronic density of states for MoO2.9 deposited hydrogenated diamond alongside the density of states for the individual structures prior 96 to deposition in Figure 5.7. When the hydrogenated diamond surface contacts the oxides, the electrons transfer from diamond valence band to acceptor’s conduction band minimum, therefore the upward band bending occurs, and the surface conductivity is initiated. Figure 5. 7. (a) Density of states for pristine hydrogenated diamond and molten MoO2.9. (b) Density of states for MoO2.9 deposited hydrogenated diamond. 97 5.4 Summary In this work, we used reactive molecular dynamics simulations and density functional theory to study the deposition of MoO3-x on hydrogenated diamond (111) surface. We proposed a three-step workflow: (1) Performing RMD simulations to generate the thermalized interfacial structure, (2) employing DFT to further optimize the interfacial structure, (3) conducting DFT calculations to study the electronic structures. The atomistic arrangement is more compact and Mo-O bond is stronger as vacancy level increases. Difference in the Bader charges after deposition reveals the net charge transfer due to deposition. Our results show the molybdenum oxides as effective electron accepting materials, which are consistent with experiments. An increase in charge transfer for higher Mo oxidation state is observed, which results in resistance decrease and current increase in electrical device. This monotonic enhancement in charge transfer as a function of oxidation state provides guidance for engineering the STD process to maximize the charge transfer. 98 Chapter 6: Conclusions In this thesis, I studied the oxidation mechanism of ZrS2, photoexcitation-induced GeTe amorphization, and also examined the surface transfer doping behavior of MoO3-x deposited hydrogenated diamond. To study the oxidation mechanism of ZrS2, reactive force field was developed and validated, and we have further optimized the reactive force field parameters for Zr/O/S using MOGA. Using the optimized force field, we have performed RMD simulations to study oxidation of ZrS2. 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. Our RMD simulations supported by first-principles calculations elucidate the atomistic mechanisms and pressure-dependence of oxidation of pristine ZrS2 (001) surfaces. The initial adsorption of oxygen atoms on the surface is driven by the large binding energy. The work indicates that oxygen partial pressure is a promising control parameter for active oxidation for future 2D electronics. Within the specified time frame, we observed that pressure can selectively expose three stages of oxide growth according to the Deal-Grove model: (1) initial incubation, (2) reaction-limited linear growth, and (3) diffusion-limited parabolic growth, with increased pressure. A trade-off between the growth speed and the quality of the grown oxide structure is identified, and the morphology of semiconductor/oxide interfaces is found to be influenced through pressure control. Initial oxidation stage is slow layer-by-layer oxidation, followed by the local breakdown of van der Waals gaps which accelerates the diffusion of oxygen, creating a kink-mediated amorphous oxide growth front, with kinetics well described by the conventional Deal-Grove model. The transition time from layer-by-layer oxidation to continuous oxidation, as well as that between reaction-limited linear and diffusion-limited parabolic oxidation within the Deal-Grove regime, is effectively regulated 99 by adjusting the pressure. We expect that the kinetics may show deviation from the Deal-Grove model due to the out-diffusion of byproducts, which will be a topic of future study involving even longer time scales than our current simulation time window. Apart from oxygen partial pressure, active oxidation control subjects including plasma-, ultraviolet- and laser-assisted oxidation processing (116, 117) are worthy of intensive exploration. In the thesis, we also investigated the photoexcitation-induced nonthermal amorphization process in GeTe by performing NAQMD that clarify the time evolution of GeTe excited states, and the corresponding structural changes induced by different excitation levels. Simulations are performed by photo-exciting n = 2.6, 4.0, 5.2, and 7.5% of valence electrons and following the time evolution of GeTe electronic structure at a temperature of 700 K. We found that photoexcitation at n = 4.0% is able to trigger local disorder at picosecond time scale. The results show a swift charge transfer from Te to Ge after photoexcitation. The electronic density of states calculations show crossings of the Te p to the conduction band and the Ge p to the valence band, together with the crossing of Ge s from the conduction band to the valence band after photoexcitation. The behavior of these crossings combined with band structure analysis suggest that photoexcitation plays an important role on the behavior of the Te p, Ge p, and Ge s electrons. At the onset of the amorphization process, Ge atoms shift from octahedral to tetrahedral sites. Electronic excitation in GeTe leads to increased Ge-Ge and decreased Ge-Te bonding interaction. The amorphization process further features the formation of Ge-Ge ‘wrong bonds’, which seems to stabilize the GeTe amorphous Ge tetrahedral motif, facilitating the loss of long-range order. The presence of weak bonds whose rupture leads to collapse of the ordered phase is crucial for amorphization of GeTe. We also utilized RMD and DFT to study the surface transfer doping due to depositing MoO3-x on hydrogenated diamond (111) surface. As more oxygen vacancies are introduced into the material, 100 the local atomic arrangement is observed to become more compact, and the Mo-O bonds are found to get stronger. Further we employed first-principles calculations based on density functional theory to investigate the charge transfer and electronic structures. Bader charge calculations further confirm MoO3-x as effective surface electron acceptors. The change in Bader charge and notable shift of the electronic density of states alignment upon doping reveal that oxygen vacancy limits the hole-doping capability of MoO3-x. Spatial distribution of doped holes is characterized after deposition, which demonstrates a widespread distribution of hole density, so-called twodimensional hole gas, at the interface that can support exceptional transport properties. 101 References 1. A. D. Yoffe, Layer Compounds. Annual Review of Materials Science 3, 147-170 (1973). 2. A. D. Yoffe, Low-dimensional systems: quantum size effects and electronic properties of semiconductor microcrystallites (zero-dimensional systems) and some quasi-twodimensional systems. Advances in Physics 42, 173-262 (1993). 3. A. S. Mayorov et al., Micrometer-Scale Ballistic Transport in Encapsulated Graphene at Room Temperature. Nano Letters 11, 2396-2399 (2011). 4. M. W. Lin et al., Room-temperature high on/off ratio in suspended graphene nanoribbon field-effect transistors. Nanotechnology 22, 265201 (2011). 5. X. Li, X. Wang, L. Zhang, S. Lee, H. Dai, Chemically Derived, Ultrasmooth Graphene Nanoribbon Semiconductors. Science 319, 1229-1232 (2008). 6. M. Y. Han, B. Özyilmaz, Y. Zhang, P. Kim, Energy Band-Gap Engineering of Graphene Nanoribbons. Physical Review Letters 98, 206805 (2007). 7. R. Balog et al., Bandgap opening in graphene induced by patterned hydrogen adsorption. Nature Materials 9, 315-319 (2010). 8. J. A. Wilson, A. D. Yoffe, The transition metal dichalcogenides discussion and interpretation of the observed optical, electrical and structural properties. Advances in Physics 18, 193-335 (1969). 9. Y. Li, J. Kang, J. Li, Indirect-to-direct band gap transition of the ZrS2 monolayer by strain: first-principles calculations. Rsc Advances 4, 7396-7401 (2014). 10. G. Fiori et al., Electronics based on two-dimensional materials. Nat Nanotechnol 9, 768- 779 (2014). 11. W. Zhang, Z. Huang, W. Zhang, Y. Li, Two-dimensional semiconductors with possible high room temperature mobility. Nano Research 7, 1731-1737 (2014). 12. Y. Wen, Y. Zhu, S. Zhang, Low temperature synthesis of ZrS2 nanoflakes and their catalytic activity. RSC Advances 5, 66082-66085 (2015). 13. J. T. Jang et al., Ultrathin zirconium disulfide nanodiscs. J Am Chem Soc 133, 7636-7639 (2011). 14. R. J. Toh, Z. Sofer, M. Pumera, Catalytic properties of group 4 transition metal dichalcogenides (MX2; M = Ti, Zr, Hf; X = S, Se, Te). Journal of Materials Chemistry A 4, 18322-18334 (2016). 15. M. Zhang et al., Controlled Synthesis of ZrS2 Monolayer and Few Layers on Hexagonal Boron Nitride. Journal of the American Chemical Society 137, 7051-7054 (2015). 102 16. H. Qiu et al., Electrical characterization of back-gated bi-layer MoS2 field-effect transistors and the effect of ambient on their performances. Applied Physics Letters 100, 123104 (2012). 17. R. Addou et al., Impurities and Electronic Property Variations of Natural MoS2 Crystal Surfaces. ACS Nano 9, 9124-9133 (2015). 18. C. Su et al., Waterproof molecular monolayers stabilize 2D materials. Proceedings of the National Academy of Sciences 116, 20844-20849 (2019). 19. Y. Li, Z. Zhou, S. Zhang, Z. Chen, MoS2 Nanoribbons: High Stability and Unusual Electronic and Magnetic Properties. Journal of the American Chemical Society 130, 16739- 16744 (2008). 20. S. S. Grønborg et al., Synthesis of Epitaxial Single-Layer MoS2 on Au(111). Langmuir 31, 9700-9706 (2015). 21. G. Mirabelli et al., Air sensitivity of MoS2, MoSe2, MoTe2, HfS2, and HfSe2. Journal of Applied Physics 120, 125102 (2016). 22. A. Cruz, Z. Mutlu, M. Ozkan, C. S. Ozkan, Raman investigation of the air stability of 2H polytype HfSe2 thin films. MRS Communications 8, 1191-1196 (2018). 23. C. H. Lee et al., Tungsten Ditelluride: a layered semimetal. Sci Rep 5, 10013 (2015). 24. J. Gao et al., Aging of Transition Metal Dichalcogenide Monolayers. ACS Nano 10, 2628- 2635 (2016). 25. Q. Li et al., Unveiling chemical reactivity and oxidation of 1T-phased group VI disulfides. Physical Chemistry Chemical Physics 21, 17010-17017 (2019). 26. M. Wuttig, N. Yamada, Phase-change materials for rewriteable data storage. Nature materials 6, 824-832 (2007). 27. G. Bruns et al., Nanosecond switching in GeTe phase change memory cells. Applied Physics Letters 95, 043108 (2009). 28. P. Fons et al., Phase transition in crystalline GeTe: Pitfalls of averaging effects. Physical Review B 82, 155209 (2010). 29. A. V. Kolobov, P. Fons, J. Tominaga, M. Hase, Excitation-Assisted Disordering of GeTe and Related Solids with Resonant Bonding. The Journal of Physical Chemistry C 118, 10248-10253 (2014). 30. G. Lucovsky, R. M. White, Effects of Resonance Bonding on the Properties of Crystalline and Amorphous Semiconductors. Physical Review B 8, 660-667 (1973). 103 31. J. Robertson, K. Xiong, P. W. Peacock, Electronic and atomic structure of Ge2Sb2Te5 phase change memory material. Thin Solid Films 515, 7538-7541 (2007). 32. K. Shportko et al., Resonant bonding in crystalline phase-change materials. Nature Materials 7, 653-658 (2008). 33. A. V. Kolobov, P. Fons, J. Tominaga, Local instability of p-type bonding makes amorphous GeTe a lone-pair semiconductor. Physical Review B 87, 155204 (2013). 34. M. Wautelet, Cohesion of Solids under Very High Electronic Excitation Conditions. physica status solidi (b) 138, 447-456 (1986). 35. S. Yip, Handbook of materials modeling. (Springer, 2005). 36. H. Stanley et al., Unsolved mysteries of water in its liquid and glassy phases. Journal of Physics: Condensed Matter 12, A403 (2000). 37. A. Rahman, Correlations in the motion of atoms in liquid argon. Physical Review 136, A405 (1964). 38. A. C. T. van Duin, S. Dasgupta, F. Lorant, W. A. Goddard, ReaxFF: A Reactive Force Field for Hydrocarbons. The Journal of Physical Chemistry A 105, 9396-9409 (2001). 39. S. Yuan et al., Insights into the surface oxidation modification mechanism of nano-diamond: An atomistic understanding from ReaxFF simulations. Applied Surface Science 540, 148321 (2021). 40. A. C. T. van Duin et al., Development and Validation of a ReaxFF Reactive Force Field for Cu Cation/Water Interactions and Copper Metal/Metal Oxide/Metal Hydroxide Condensed Phases. The Journal of Physical Chemistry A 114, 9507-9514 (2010). 41. A. Włodarczyk, M. Uchroński, A. Podsiadły-Paszkowska, J. Irek, B. M. Szyja, Mixing ReaxFF parameters for transition metal oxides using force-matching method. Journal of Molecular Modeling 28, 8 (2021). 42. Y.-S. Zhang et al., Molecular dynamics simulations of the initial oxidation process on ferritic Fe–Cr alloy surfaces. RSC Advances 12, 9501-9511 (2022). 43. T. P. Senftle, R. J. Meyer, M. J. Janik, A. C. T. van Duin, Development of a ReaxFF potential for Pd/O and application to palladium oxide formation. The Journal of Chemical Physics 139, 044109 (2013). 44. D. Raymand, A. C. T. van Duin, M. Baudin, K. Hermansson, A reactive force field (ReaxFF) for zinc oxide. Surface Science 602, 1020-1031 (2008). 45. K. A. Roshan, M. K. Talkhoncheh, M. Y. Sengul, D. J. Miller, A. C. T. van Duin, Optimization of ReaxFF Reactive Force Field Parameters for Cu/Si/O Systems via Neural 104 Network Inversion with Application to Copper Oxide Interaction with Silicon. The Journal of Physical Chemistry C 127, 20445-20458 (2023). 46. D. Fantauzzi et al., Development of a ReaxFF potential for Pt–O systems describing the energetics and dynamics of Pt-oxide formation. Physical Chemistry Chemical Physics 16, 23118-23133 (2014). 47. H. Azizi Toupkanloo, M. Fathollahi, Molecular Dynamics Simulation of Al/NiO Thermite Reaction Using Reactive Force Field (ReaxFF). Physical Chemistry Research 5, 221-237 (2017). 48. F. Fiesinger et al., Development of a Mg/O ReaxFF Potential to describe the Passivation Processes in Magnesium-Ion Batteries**. ChemSusChem 16, e202201821 (2023). 49. W.-X. Song, S.-J. Zhao, Development of the ReaxFF reactive force field for aluminum– molybdenum alloy. Journal of Materials Research 28, 1155-1164 (2013). 50. T. P. Senftle et al., The ReaxFF reactive force-field: development, applications and future directions. npj Computational Materials 2, 15011 (2016). 51. M. Bhati, T. P. Senftle, Identifying Adhesion Properties at Si/Polymer Interfaces with ReaxFF. The Journal of Physical Chemistry C 123, 27036-27047 (2019). 52. F. A. Soria, W. Zhang, P. A. Paredes-Olivera, A. C. T. van Duin, E. M. Patrito, Si/C/H ReaxFF Reactive Potential for Silicon Surfaces Grafted with Organic Molecules. The Journal of Physical Chemistry C 122, 23515-23527 (2018). 53. A. C. T. van Duin, S. Dasgupta, F. Lorant, W. A. Goddard, ReaxFF: A reactive force field for hydrocarbons. J. Phys. Chem. A 105, 9396-9409 (2001). 54. P. Hohenberg, W. Kohn, Inhomogeneous Electron Gas. Physical Review 136, B864-B871 (1964). 55. J. P. Perdew, K. Burke, M. Ernzerhof, Generalized Gradient Approximation Made Simple. Physical Review Letters 77, 3865-3868 (1996). 56. J. Sun, A. Ruzsinszky, J. P. Perdew, Strongly Constrained and Appropriately Normed Semilocal Density Functional. Physical Review Letters 115, 036402 (2015). 57. J. Heyd, G. E. Scuseria, M. Ernzerhof, Hybrid functionals based on a screened Coulomb potential. The Journal of Chemical Physics 118, 8207-8215 (2003). 58. V. I. Anisimov, J. Zaanen, O. K. Andersen, Band theory and Mott insulators: Hubbard U instead of Stoner I. Physical Review B 44, 943-954 (1991). 59. S. Grimme, J. Antony, S. Ehrlich, H. Krieg, A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. The Journal of Chemical Physics 132, 154104 (2010). 105 60. N. Troullier, J. L. Martins, Efficient pseudopotentials for plane-wave calculations. Physical Review B 43, 1993-2006 (1991). 61. D. Marx, J. Hutter, Ab Initio Molecular Dynamics. (Cambridge University Press, Cambridge, UK, 2009). 62. F. Shimojo et al., QXMD: An open-source program for nonadiabatic quantum molecular dynamics. SoftwareX 10, 100307 (2019). 63. F. Shimojo et al., Large nonadiabatic quantum molecular dynamics simulations on parallel computers. Computer Physics Communications 184, 1-8 (2013). 64. J. Harris, Simplified method for calculating the energy of weakly interacting fragments. Physical Review B 31, 1770-1779 (1985). 65. W. M. C. Foulkes, R. Haydock, Tight-binding models and density-functional theory. Physical Review B 39, 12520-12536 (1989). 66. A. S. Torralba, D. R. Bowler, T. Miyazaki, M. J. Gillan, Non-self-consistent DensityFunctional Theory Exchange-Correlation Forces for GGA Functionals. Journal of Chemical Theory and Computation 5, 1499-1505 (2009). 67. G. P. Kerker, Efficient iteration scheme for self-consistent pseudopotential calculations. Physical Review B 23, 3082-3084 (1981). 68. N. L. Doltsinis, D. Marx, Nonadiabatic Car-Parrinello molecular dynamics. Physical Review Letters 88, 166402 (2002). 69. W. R. Duncan, C. F. Craig, O. V. Prezhdo, Time-domain ab initio study of charge relaxation and recombination in dye-sensitized TiO2. Journal of the American Chemical Society 129, 8528-8543 (2007). 70. C. P. Hu, H. Hirai, O. Sugino, Nonadiabatic couplings from time-dependent density functional theory: formulation in the Casida formalism and practical scheme within modified linear response. J Chem Phys 127, 064103 (2007). 71. E. Tapavicza, I. Tavernelli, U. Rothlisberger, Trajectory surface hopping within linear response time-dependent density-functional theory. Physical Review Letters 98, 023001 (2007). 72. S. Ohmura et al., Atomistic mechanisms of rapid energy transport in light-harvesting molecules. Applied Physics Letters 98, 113302 (2011). 73. W. Mou, S. Ohmura, F. Shimojo, A. Nakano, Molecular control of photoexcited charge transfer and recombination at a quaterthiophene/zinc oxide interface. Applied Physics Letters 100, 203306 (2012). 106 74. W. Mou et al., Enhanced charge transfer by phenyl groups at a rubrene/C60 interface. J Chem Phys 136, 184705 (2012). 75. J. C. Tully, Molecular dynamics with electronic transitions. J Chem Phys 93, 1061-1071 (1990). 76. J. R. Schmidt, P. V. Parandekar, J. C. Tully, Mixed quantum-classical equilibrium: Surface hopping. J Chem Phys 129, 044104 (2008). 77. E. K. U. Gross, W. Kohn, Local density-functional theory of frequency-dependent linear response. Physical Review Letters 55, 2850-2853 (1985). 78. I. Vasiliev, S. Ogut, J. R. Chelikowsky, Ab initio excitation spectra and collective electronic response in atoms and clusters. Physical Review Letters 82, 1919-1922 (1999). 79. C. Y. Yam, S. Yokojima, G. H. Chen, Linear-scaling time-dependent density-functional theory. Physical Review B 68, 153105 (2003). 80. S. Meng, E. Kaxiras, Real-time, local basis-set implementation of time-dependent density functional theory for excited state dynamics simulations. J Chem Phys 129, 054110 (2008). 81. V. Soto-Verdugo, H. Metiu, E. Gwinn, The properties of small Ag clusters bound to DNA bases. J Chem Phys 132, 195102 (2010). 82. M. A. L. Marques et al., Time-Dependent Density Functional Theory. Lecture Notes in Physics (Springer-Verlag, Berlin, Germany, 2006), vol. 706. 83. A. Stukowski, Visualization and analysis of atomistic simulation data with OVITO–the Open Visualization Tool. Modelling and Simulation in Materials Science and Engineering 18, 015012 (2010). 84. K. Momma, F. Izumi, VESTA: a three-dimensional visualization system for electronic and structural analysis. Journal of Applied Crystallography 41, 653-658 (2008). 85. W. Humphrey, A. Dalke, K. Schulten, VMD: Visual molecular dynamics. Journal of Molecular Graphics 14, 33-38 (1996). 86. H. Childs et al., VisIt: An End-User Tool for Visualizing and Analyzing Very Large Data. Proceed SciDAC, 1-16 (2011). 87. L. Ju, AtomEye: an efficient atomistic configuration viewer. Modelling and Simulation in Materials Science and Engineering 11, 173 (2003). 88. J. P. Ahrens, B. Geveci, C. C. Law, in The Visualization Handbook. (2005). 89. S. S. Jo et al., Growth Kinetics and Atomistic Mechanisms of Native Oxidation of ZrSxSe2–x and MoS2 Crystals. Nano Letters 20, 8592-8599 (2020). 90. L. Yang et al., Unveiling oxidation mechanism of bulk ZrS2. MRS Advances, (2021). 107 91. L. Yang, R. Jaramillo, R. K. Kalia, A. Nakano, P. Vashishta, Pressure-Controlled Layer-byLayer to Continuous Oxidation of ZrS2(001) Surface. ACS Nano 17, 7576-7583 (2023). 92. M. Moustafa, T. Zandt, C. Janowitz, R. Manzke, Growth and band gap determination of the ZrSxSe2 single crystal series. Physical Review B 80, 035206 (2009). 93. L. Li et al., Electrical Transport and High-Performance Photoconductivity in Individual ZrS2 Nanobelts. Advanced Materials 22, 4151-4156 (2010). 94. M. Hamada et al., ZrS2 symmetrical-ambipolar FETs with near-midgap TiN film for both top-gate electrode and Schottky-barrier contact. Japanese Journal of Applied Physics 60, SBBH05 (2021). 95. X. Wang et al., Large scale ZrS2 atomically thin layers. Journal of Materials Chemistry C 4, 3143-3148 (2016). 96. A. H. Reshak, S. Auluck, Theoretical investigation of the electronic and optical properties of ZrX2 (X=S, Se and Te). Physica B: Condensed Matter 353, 230-237 (2004). 97. L. Li et al., High-performance Schottky solar cells using ZrS2 nanobelt networks. Energy & Environmental Science 4, 2586-2590 (2011). 98. S. Zhou, L. Wu, Phase Separation and Properties of UV-Curable Polyurethane/Zirconia Nanocomposite Coatings. Macromolecular Chemistry and Physics 209, 1170-1181 (2008). 99. K. Luo, S. Zhou, L. Wu, G. Gu, Dispersion and Functionalization of Nonaqueous Synthesized Zirconia Nanocrystals via Attachment of Silane Coupling Agents. Langmuir 24, 11497-11505 (2008). 100. T. Koch, P. Ziemann, Zr-silicide formation during the epitaxial growth of Y-stabilized zirconia films on Si(100) and its avoidance by ion beam assisted deposition at a reduced temperature. Applied Surface Science 99, 51-57 (1996). 101. J. H. Shim, C.-C. Chao, H. Huang, F. B. Prinz, Atomic Layer Deposition of YttriaStabilized Zirconia for Solid Oxide Fuel Cells. Chemistry of Materials 19, 3850-3854 (2007). 102. T. M. Miller, V. H. Grassian, Environmental Catalysis: Adsorption and Decomposition of Nitrous Oxide on Zirconia. Journal of the American Chemical Society 117, 10969-10975 (1995). 103. Y. Li et al., Effect of calcium salts on isosynthesis over ZrO2 catalysts. Journal of Molecular Catalysis A: Chemical 175, 267-275 (2001). 104. D. Chen et al., Synthesis of Monodisperse Mesoporous Titania Beads with Controllable Diameter, High Surface Areas, and Variable Pore Diameters (14−23 nm). Journal of the American Chemical Society 132, 4438-4444 (2010). 108 105. D. Lu et al., A Novel Nanoparticle-Based Disposable Electrochemical Immunosensor for Diagnosis of Exposure to Toxic Organophosphorus Agents. Advanced Functional Materials 21, 4371-4378 (2011). 106. M. Zhou, A. Ahmad, Synthesis, processing and characterization of calcia-stabilized zirconia solid electrolytes for oxygen sensing applications. Materials Research Bulletin 41, 690-696 (2006). 107. A. Subramanian, P. W. Carr, C. V. McNeff, Use of spray-dried zirconia microspheres in the separation of immunoglobulins from cell culture supernatant. Journal of Chromatography A 890, 15-23 (2000). 108. B. Yan, C. V. McNeff, F. Chen, P. W. Carr, A. V. McCormick, Control of Synthesis Conditions to Improve Zirconia Microspheres for Ultrafast Chromatography. Journal of the American Ceramic Society 84, 1721-1727 (2001). 109. Y. Xia et al., Synthesis of Transparent Aqueous ZrO2 Nanodispersion with a Controllable Crystalline Phase without Modification for a High-Refractive-Index Nanocomposite Film. Langmuir 34, 6806-6813 (2018). 110. U. C. Uduh et al., Sol–gel synthesis, optical and structural characterization of ZrOS nanopowder. Journal of Sol-Gel Science and Technology 71, 79-85 (2014). 111. Y. Zhang, H.-X. Chen, L. Duan, J.-B. Fan, The structural, electronic, elastic, dielectric, dynamical, thermal and optical properties of Janus ZrOS monolayer: A first-principles investigation. Solid State Communications 327, 114207 (2021). 112. A. Jaramillo-Botero, S. Naserifar, W. A. Goddard, General Multiobjective Force Field Optimization Framework, with Application to Reactive Force Fields for Silicon Carbide. Journal of Chemical Theory and Computation 10, 1426-1439 (2014). 113. J. P. Larentzos, B. M. Rice, E. F. C. Byrd, N. S. Weingarten, J. V. Lill, Parameterizing Complex Reactive Force Fields Using Multiple Objective Evolutionary Strategies (MOES). Part 1: ReaxFF Models for Cyclotrimethylene Trinitramine (RDX) and 1,1- Diamino-2,2-dinitroethene (FOX-7). Journal of Chemical Theory and Computation 11, 381-391 (2015). 114. K. Reidy et al., Kinetic control for planar oxidation of MoS2. arXiv preprint arXiv:2211.16789, (2022). 115. Z. Hu et al., Water-Catalyzed Oxidation of Few-Layer Black Phosphorous in a Dark Environment. Angew Chem Int Ed Engl 56, 9131-9135 (2017). 116. S. Lai et al., HfO2/HfS2 hybrid heterostructure fabricated via controllable chemical conversion of two-dimensional HfS2. Nanoscale 10, 18758-18766 (2018). 117. Y. Y. Illarionov et al., Insulators for 2D nanoelectronics: the gap to bridge. Nature Communications 11, 3385 (2020). 109 118. B. E. Deal, A. S. Grove, General Relationship for the Thermal Oxidation of Silicon. Journal of Applied Physics 36, 3770-3778 (1965). 119. Y. Shao et al., Advances in molecular quantum chemistry contained in the Q-Chem 4 program package. Molecular Physics 113, 184-215 (2015). 120. C. Lee, W. Yang, R. G. Parr, Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys Rev B Condens Matter 37, 785-789 (1988). 121. A. D. Becke, A new mixing of Hartree–Fock and local density‐functional theories. The Journal of Chemical Physics 98, 1372-1377 (1993). 122. G. Kresse, J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Physical review B 54, 11169 (1996). 123. S. Plimpton, Fast parallel algorithms for short-range molecular dynamics. Journal of computational physics 117, 1-19 (1995). 124. K. Deb, A. Pratap, S. Agarwal, T. Meyarivan, A fast and elitist multiobjective genetic algorithm: NSGA-II. IEEE Transactions on Evolutionary Computation 6, 182-197 (2002). 125. A. Krishnamoorthy et al., Evolutionary multi-objective optimization and Pareto-frontal uncertainty quantification of interatomic forcefields for thermal conductivity simulations. Computer Physics Communications 254, 107337 (2020). 126. A. Mishra et al., Multiobjective genetic training and uncertainty quantification of reactive force fields. npj Computational Materials 4, 42 (2018). 127. W. Gropp, E. Lusk, A. Skjellum, Using MPI. (MIT Press, Cambridge, MA, ed. Third, 2014), pp. 336. 128. K.-i. Nomura, Y.-C. Chen, R. K. Kalia, A. Nakano, P. Vashishta, Defect migration and recombination in nanoindentation of silica glass. Applied Physics Letters 99, 111906 (2011). 129. F. Shimojo, S. Ohmura, R. K. Kalia, A. Nakano, P. Vashishta, Molecular Dynamics Simulations of Rapid Hydrogen Production from Water Using Aluminum Clusters as Catalyzers. Physical Review Letters 104, 126102 (2010). 130. N. Agmon, The Grotthuss mechanism. Chemical Physics Letters 244, 456-462 (1995). 131. G. Kresse, J. Hafner, Ab initio molecular dynamics for liquid metals. Physical Review B 47, 558-561 (1993). 132. P. E. Blöchl, Projector augmented-wave method. Physical Review B 50, 17953-17979 (1994). 110 133. H. J. Monkhorst, J. D. Pack, Special points for Brillouin-zone integrations. Physical Review B 13, 5188-5192 (1976). 134. A. Soon, M. Todorova, B. Delley, C. Stampfl, Oxygen adsorption and stability of surface oxides on Cu(111): A first-principles investigation. Physical Review B 73, 165424 (2006). 135. K. Reuter, M. Scheffler, Composition, structure, and stability of RuO2(110) as a function of oxygen pressure. Physical Review B 65, 035406 (2001). 136. JANAF thermochemical tables [electronic resource] / D.R. Stull and H. Prophet, project directors. D. R. Stull, H. Prophet, S. United States. National Bureau of, Eds., NSRDSNBS ; 37. (U.S. Dept. of Commerce, National Bureau of Standards, Washington, D.C, 1971). 137. A. P. Thompson et al., LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales. Computer Physics Communications 271, 108171 (2022). 138. S. Nosé, A unified formulation of the constant temperature molecular dynamics methods. The Journal of Chemical Physics 81, 511-519 (1984). 139. W. G. Hoover, Canonical dynamics: Equilibrium phase-space distributions. Physical Review A 31, 1695-1697 (1985). 140. F. J. Buijnsters, A. Fasolino, M. I. Katsnelson, Motion of Domain Walls and the Dynamics of Kinks in the Magnetic Peierls Potential. Physical Review Letters 113, 217202 (2014). 141. T. Zhu, J. Li, S. Yip, Atomistic Configurations and Energetics of Crack Extension in Silicon. Physical Review Letters 93, 205504 (2004). 142. P. Bouvier, E. Djurado, G. Lucazeau, T. Le Bihan, High-pressure structural evolution of undoped tetragonal nanocrystalline zirconia. Physical Review B 62, 8731-8737 (2000). 143. D. K. Smith, W. Newkirk, The crystal structure of baddeleyite (monoclinic ZrO2) and its relation to the polymorphism of ZrO2. Acta Crystallographica 18, 983-991 (1965). 144. M. Li et al., Unusual layer-by-layer growth of epitaxial oxide islands during Cu oxidation. Nature Communications 12, 2781 (2021). 145. Q. Zhu et al., Defect-driven selective metal oxidation at atomic scale. Nature Communications 12, 558 (2021). 146. Q. Zhao et al., Elastic, electronic, and dielectric properties of bulk and monolayer ZrS2, ZrSe2, HfS2, HfSe2 from van der Waals density-functional theory. physica status solidi (b) 254, 1700033 (2017). 111 147. T. Campbell et al., Dynamics of Oxidation of Aluminum Nanoclusters using Variable Charge Molecular-Dynamics Simulations on Parallel Computers. Physical Review Letters 82, 4866-4869 (1999). 148. L. Yang et al., Photoexcitation-Induced Nonthermal Ultrafast Loss of Long-Range Order in GeTe. The Journal of Physical Chemistry Letters 13, 10230-10236 (2022). 149. A. V. Kolobov, J. Tominaga, Chalcogenides: metastability and phase change phenomena. (Springer Science & Business Media, 2012), vol. 164. 150. V. L. Deringer, R. Dronskowski, M. Wuttig, Microscopic Complexity in Phase-Change Materials and its Role for Applications. Advanced Functional Materials 25, 6343-6359 (2015). 151. P. S. Branicio et al., Atomistic insights into the nanosecond long amorphization and crystallization cycle of nanoscale Ge2Sb2Te5: An ab initio molecular dynamics study. Physical Review Materials 2, 043401 (2018). 152. Z. Song, R. Wang, Y. Xue, S. Song, The “gene” of reversible phase transformation of phase change materials: Octahedral motif. Nano Research 15, 765-772 (2022). 153. A. V. Kolobov, M. Krbal, P. Fons, J. Tominaga, T. Uruga, Distortion-triggered loss of long-range order in solids with bonding energy hierarchy. Nature Chemistry 3, 311-316 (2011). 154. X.-B. Li et al., Role of Electronic Excitation in the Amorphization of Ge-Sb-Te Alloys. Physical Review Letters 107, 015501 (2011). 155. S. C. Tiwari et al., Photoexcitation Induced Ultrafast Nonthermal Amorphization in Sb2Te3. The Journal of Physical Chemistry Letters 11, 10242-10249 (2020). 156. E. Matsubara et al., Initial Atomic Motion Immediately Following Femtosecond-Laser Excitation in Phase-Change Materials. Physical Review Letters 117, 135501 (2016). 157. K. Singh et al., A review on GeTe thin film-based phase-change materials. Applied Nanoscience, (2021). 158. J. L. F. Da Silva, A. Walsh, H. Lee, Insights into the structure of the stable and metastable (GeTe)m(Sb2Te3)n compounds. Physical Review B 78, 224111 (2008). 159. D. Lencer, M. Salinga, M. Wuttig, Design Rules for Phase-Change Materials in Data Storage Applications. Advanced Materials 23, 2030-2058 (2011). 160. S. Caravati, M. Bernasconi, T. D. Kühne, M. Krack, M. Parrinello, Coexistence of tetrahedral- and octahedral-like sites in amorphous phase change materials. Applied Physics Letters 91, 171906 (2007). 112 161. J. Akola, R. O. Jones, Structural phase transitions on the nanoscale: The crucial pattern in the phase-change materials GeSbTe5 and GeTe. Physical Review B 76, 235201 (2007). 162. J. Hegedüs, S. R. Elliott, Microscopic origin of the fast crystallization ability of Ge–Sb– Te phase-change memory materials. Nature Materials 7, 399-405 (2008). 163. M. Krbal et al., Intrinsic complexity of the melt-quenched amorphous Ge2Sb2Te5 memory alloy. Physical Review B 83, 054203 (2011). 164. A. V. Kolobov et al., Understanding the phase-change mechanism of rewritable optical media. Nature materials 3, 703-708 (2004). 165. M. A. Paesler, D. A. Baker, G. Lucovsky, A. E. Edwards, P. C. Taylor, EXAFS study of local order in the amorphous chalcogenide semiconductor Ge2Sb2Te5. Journal of Physics and Chemistry of Solids 68, 873-877 (2007). 166. S. K. Sundaram, E. Mazur, Inducing and probing non-thermal transitions in semiconductors using femtosecond laser pulses. Nature Materials 1, 217-224 (2002). 167. A. V. Kolobov, P. Fons, J. Tominaga, Understanding Phase-Change Memory Alloys from a Chemical Perspective. Scientific Reports 5, 13698 (2015). 168. S. Nosé, A molecular dynamics method for simulations in the canonical ensemble. Molecular Physics 52, 255-268 (1984). 169. R. Car, M. Parrinello, Unified Approach for Molecular Dynamics and Density-Functional Theory. Physical Review Letters 55, 2471-2474 (1985). 170. W. Kohn, L. J. Sham, Self-Consistent Equations Including Exchange and Correlation Effects. Physical Review 140, A1133-A1138 (1965). 171. W. Kohn, P. Vashishta, in Theory of the Inhomogeneous Electron Gas, S. Lundqvist, N. H. March, Eds. (Springer US, Boston, MA, 1983), pp. 79-147. 172. S. Grimme, Semiempirical GGA-type density functional constructed with a long-range dispersion correction. Journal of Computational Chemistry 27, 1787-1799 (2006). 173. J. C. Tully, Molecular dynamics with electronic transitions. The Journal of Chemical Physics 93, 1061-1071 (1990). 174. M. E. Casida, D. Chong, Recent advances in density functional methods. Computational Chemistry: Reviews of Current Trends, (1995). 175. N. Yamada, T. Matsunaga, Structure of laser-crystallized Ge2Sb2+xTe5 sputtered thin films for use in optical memory. Journal of applied physics 88, 7020-7028 (2000). 176. S. Shamoto et al., Large displacement of germanium atoms in crystalline Ge2Sb2Te5. Applied Physics Letters 86, 081904 (2005). 113 177. R. Shayduk, W. Braun, Epitaxial films for Ge–Sb–Te phase change memory. Journal of Crystal Growth 311, 2215-2219 (2009). 178. W. Wełnic et al., Unravelling the interplay of local structure and physical properties in phase-change materials. Nature Materials 5, 56-62 (2006). 179. A. V. Kolobov et al., Why phase-change media are fast and stable: a new approach to an old problem. Japanese journal of applied physics 44, 3345 (2005). 180. J. F. Cornwell, Group theory and electronic energy bands in solids. (North-Holland Amsterdam, 1969). 181. V. L. Deringer, A. L. Tchougréeff, R. Dronskowski, Crystal Orbital Hamilton Population (COHP) Analysis As Projected from Plane-Wave Basis Sets. The Journal of Physical Chemistry A 115, 5461-5466 (2011). 182. A. Savin, R. Nesper, S. Wengert, T. F. Fässler, ELF: The Electron Localization Function. Angewandte Chemie International Edition in English 36, 1808-1832 (1997). 183. P. Kerres et al., Scaling and Confinement in Ultrathin Chalcogenide Films as Exemplified by GeTe. Small 18, 2201753 (2022). 184. S. N. Yannopoulos, M. L. Trunov, Photoplastic effects in chalcogenide glasses: A review. physica status solidi (b) 246, 1773-1785 (2009). 185. P. Fons et al., Photoassisted amorphization of the phase-change memory alloy Ge2Sb2Te5. Physical Review B 82, 041203 (2010). 186. V. L. Deringer et al., Bonding Nature of Local Structural Motifs in Amorphous GeTe. Angewandte Chemie International Edition 53, 10817-10820 (2014). 187. M. Küpers et al., Unexpected Ge–Ge Contacts in the Two-Dimensional Ge4Se3Te Phase and Analysis of Their Chemical Cause with the Density of Energy (DOE) Function. Angewandte Chemie International Edition 56, 10204-10208 (2017). 188. M. Wuttig, Phase change materials: Chalcogenides with remarkable properties due to an unconventional bonding mechanism. physica status solidi (b) 249, 1843-1850 (2012). 189. A. V. Kolobov, P. Fons, M. Krbal, J. Tominaga, Amorphous phase of GeTe-based phasechange memory alloys: Polyvalency of Ge-Te bonding and polyamorphism. physica status solidi (a) 209, 1031-1035 (2012). 190. B. J. Kooi, M. Wuttig, Chalcogenides by Design: Functionality through Metavalent Bonding and Confinement. Advanced Materials 32, 1908302 (2020). 191. J.-Y. Raty, M. Wuttig, The interplay between Peierls distortions and metavalent bonding in IV–VI compounds: comparing GeTe with related monochalcogenides. Journal of Physics D: Applied Physics 53, 234002 (2020). 114 192. A. V. Kolobov, P. Fons, J. Tominaga, S. R. Ovshinsky, Vacancy-mediated three-center four-electron bonds in GeTe-Sb2Te3 phase-change memory alloys. Physical Review B 87, 165206 (2013). 193. J. Isberg et al., High carrier mobility in single-crystal plasma-deposited diamond. (Reports). Science 297, 1670+ (2002). 194. M. Kasu, K. Ueda, Y. Yamauchi, A. Tallaire, T. Makimoto, Diamond-based RF power transistors: Fundamentals and applications. Diamond and Related Materials 16, 1010- 1015 (2007). 195. P. Strobel, M. Riedel, J. Ristein, L. Ley, Surface transfer doping of diamond. Nature 430, 439-441 (2004). 196. M. W. Geis et al., Progress Toward Diamond Power Field-Effect Transistors. physica status solidi (a) 215, 1800681 (2018). 197. K. Xing et al., MoO3 induces p-type surface conductivity by surface transfer doping in diamond. Applied Surface Science 509, 144890 (2020). 198. J. Ristein, Surface Transfer Doping of Semiconductors. Science 313, 1057-1058 (2006). 199. K. G. Crawford, I. Maini, D. A. Macdonald, D. A. J. Moran, Surface transfer doping of diamond: A review. Progress in Surface Science 96, 100613 (2021). 200. K. G. Crawford et al., Thermally Stable, High Performance Transfer Doping of Diamond using Transition Metal Oxides. Scientific Reports 8, 3342 (2018). 201. J. McGhee, V. P. Georgiev, Simulation Study of Surface Transfer Doping of Hydrogenated Diamond by MoO3 and V2O5 Metal Oxides. Micromachines 11, 433 (2020). 202. F. Maier, M. Riedel, B. Mantel, J. Ristein, L. Ley, Origin of Surface Conductivity in Diamond. Physical Review Letters 85, 3472-3475 (2000). 203. J. C. Angus et al., Chemical Vapour Deposition of Diamond. Philosophical Transactions: Physical Sciences and Engineering 342, 195-208 (1993). 204. H. Kawarada, Hydrogen-terminated diamond surfaces and interfaces. Surface Science Reports 26, 205-259 (1996). 205. J. Meyer, A. Shu, M. Kröger, A. Kahn, Effect of contamination on the electronic structure and hole-injection properties of MoO3/organic semiconductor interfaces. Applied Physics Letters 96, (2010). 206. S. B. Sinnott, D. W. Brenner, Three decades of many-body potentials in materials research. MRS Bulletin 37, 469-473 (2012). 115 207. K.-i. Nomura, R. K. Kalia, A. Nakano, P. Rajak, P. Vashishta, RXMD: A scalable reactive molecular dynamics simulator for optimized time-to-solution. SoftwareX 11, 100389 (2020). 208. G. Henkelman, A. Arnaldsson, H. Jónsson, A fast and robust algorithm for Bader decomposition of charge density. Computational Materials Science 36, 354-360 (2006). 209. R. F. W. Bader, R. F. Bader, Atoms in Molecules: A Quantum Theory. (Clarendon Press, 1990). 210. Y.-H. Kim, S. B. Zhang, Y. Yu, L. F. Xu, C. Z. Gu, Dihydrogen bonding, $p$-type conductivity, and origin of change in work function of hydrogenated diamond (001) surfaces. Physical Review B 74, 075329 (2006). 211. Y. Yang, F. A. Koeck, X. Wang, R. J. Nemanich, Surface transfer doping of MoO3 on hydrogen terminated diamond with an Al2O3 interfacial layer. Applied Physics Letters 120, 191602 (2022). 212. S. A. O. Russell et al., Surface transfer doping of diamond by MoO3: A combined spectroscopic and Hall measurement study. Applied Physics Letters 103, 202112 (2013).

Abstract (if available)

Abstract The end of Moore’s law poses enormous scientific and technological challenges. First-principles quantum molecular dynamics (QMD) and reactive molecular dynamics (RMD) simulations could revolutionize future semiconductor technology by providing fundamental mechanism at atomic and electronic scale underlining the synthesis and performance of these future semiconductor devices.

This thesis contains results from three projects: (1) Controlled oxidation of 2D transition metal dichalcogenide (TMD) to synthesize semiconductor-oxide heterostructure; (2) nonthermal switching of phase change material (PCM); and (3) performance of MoO3 for hole doping of hydrogenated diamond.

We investigate controlled oxidation of 2D transition metal dichalcogenide to synthesize semiconductor-oxide heterostructure: While transition metal dichalcogenides (TMDCs) have been reported as potential candidates for next-generation electronic and optoelectronic devices, their native oxidation remains a major issue in achieving their long-term stability. Among TMDCs, ZrS2 shows superior electrical properties. However, ZrS2 has been reported to easily oxidize under ambient conditions. The oxidation processes deserve widespread attention. A detailed atomistic understanding into the oxidation mechanism is essential for designing better semiconductor devices.

To accomplish this goal, we first optimized the reactive force field ReaxFF parameters for Zr/O/S using multi-objective genetic algorithm (MOGA) implemented in our EZFF force field training software, so that the new force field is able to give consistent charge values and bond-population dynamics with quantum-mechanically computed information. The detailed steps for force field optimization are provided along with the accuracy of the optimal force field. Later, we performed RMD simulations employing our optimized force field to study the oxidation processes of ZrS2. In our RMD simulations, (210) and (001) surfaces of ZrS2 are exposed to O2 atmosphere separately at a simulation temperature of 1200K. The oxidation rate at (210) surface was determined to be higher than that at (001) surface. The oxidation rate is highly dependent on the initial adsorption of oxygen molecules on the surface. Oxidation rate on ZrS2 (001) surface is found to be accelerated under higher pressure at a simulation temperature of 300 K. In order to investigate the pressure effect, the Deal-Grove model is discussed for its applicability in ZrS2 oxidation system, and the initial oxidation stage is examined using density functional theory. Our simulation also found layer-by-layer oxidation to amorphous-oxide-mediated continuous oxidation. Visualization of oxide surface plots are shown to elucidate this process. The atomistic understanding lays the foundation for describing oxidation of ZrS2. It is hoped that this understanding will be valuable for actual device processing.

Second, the nonthermal switching of phase change materials for ultrafast memory devices is studied. The contrasting electrical and optical properties of the crystalline and amorphous phases of GeTe, a primary phase-change material in the (GeTe)x(Sb2Te3)1-x alloy family, are of great importance for non-volatile storage devices and advanced computing architectures. Here, by using nonadiabatic quantum molecular dynamics to investigate the evolution of GeTe excited states, we reveal a photoexcitation-induced picosecond nonthermal path for the loss of long-range order in GeTe. A valence electron excitation threshold of 4% is found to trigger local disorder by switching Ge atoms from octahedral to tetrahedral sites and promoting Ge-Ge bonding. The resulting loss of long-range order for a higher valence electron excitation fraction is achieved without fulfilling the Lindemann criterion for melting, therefore utilizing a nonthermal path. The photoexcitation- induced structural disorder is accompanied by charge transfer from Te to Ge, Ge-Te bonding-to- antibonding, and Ge-Ge antibonding-to-bonding change, triggering Ge-Te bond breaking and promoting the formation of Ge-Ge wrong bonds. These results provide an electronic-structure basis to understand the photoexcitation-induced ultrafast changes in the structure and properties of GeTe and other phase change materials.

In the third part, we investigate the performance of MoO3 for hole doping of hydrogenated diamond for energy efficient high-power electronics. Hydrogen-terminated diamond field-effect transistors are promising candidates for next-generation high-power electronics. While hydrogenated diamond surface deposited with MoO3 was shown to exhibit high p-type surface conductivity, atomistic and electronic structures of the interface remain elusive, including the likely effects of oxygen vacancies. Here, we develop and optimize a reactive force field for the C-H-Mo-O system, which describes charge transfer and bond formation processes with quantum-mechanical accuracy. Reactive molecular dynamics simulations are performed using this force field to elucidate the deposition mechanism of MoO3-x on hydrogenated diamond (111) surface. Electronic density of states alignment and charge transfer at the interface are studied using first-principles calculations based on density functional theory for selected thermalized structures taken from the reactive molecular dynamics trajectories. We observe shifting of the electronic density of states alignment and higher charge transfer for higher Mo oxidation state. Such atomistic and electronic details provide mechanistic understanding of the surface transfer doping process.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Quantum molecular dynamics and machine learning for non-equilibrium processes in quantum materials

PDF

Molecular dynamics study of energetic materials

PDF

High throughput computational framework for synthesis and accelerated discovery of dielectric polymer materials using polarizable reactive molecular dynamics and graph neural networks

PDF

Nanomaterials under extreme environments: a study of structural and dynamic properties using reactive molecular dynamics simulations

PDF

In2O3 COVID-19 biosensors and two-dimensional materials: synthesis and applications

PDF

Neural network for molecular dynamics simulation and design of 2D materials

PDF

Controlling electronic properties of two-dimensional quantum materials: simulation at the nexus of the classical and quantum computing eras

PDF

Modeling of bicontinuous nanoporous materials and investigation of their mechanical and hydraulic properties

PDF

Molecular dynamics studies of protein aggregation in unbounded and confined media

PDF

Large-scale molecular dynamics simulations of nano-structured materials

PDF

Multi-scale quantum dynamics and machine learning for next generation energy applications

PDF

Molecular-scale studies of mechanical phenomena at the interface between two solid surfaces: from high performance friction to superlubricity and flash heating

PDF

Real time surface analysis of complex oxide thin films during pulsed laser deposition

PDF

Physics informed neural networks and electrostrictive cavitation in water

PDF

Thermal properties of silicon carbide and combustion mechanisms of aluminum nanoparticle

PDF

Molecular dynamics simulation study of initial protein unfolding induced by the photo-responsive surfactants, azoTAB

PDF

Shock-induced poration, cholesterol flip-flop and small interfering RNA transfection in a phospholipid membrane: multimillion atom, microsecond molecular dynamics simulations

PDF

Adiabatic and non-adiabatic molecular dynamics in nanoscale systems: theory and applications

PDF

Interactions of planetary surfaces with space environments and their effects on volatile formation and transport: atomic scale simulations

PDF

Phase change heterostructures for electronic and photonic applications

Asset Metadata

Creator Yang, Liqiu (author)

Core Title Reactive and quantum molecular dynamics study of materials: from oxidation and amorphization to interfacial charge transfer

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Chemical Engineering

Degree Conferral Date 2023-12

Publication Date 12/12/2024

Defense Date 12/12/2023

Publisher Los Angeles, California (original), University of Southern California (original), University of Southern California. Libraries (digital)

Tag density functional theory,molecular dynamics,OAI-PMH Harvest,oxidation,phase change material,surface transfer doping,TMDC

Format theses (aat)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Vashishta, Priya (committee chair)

Creator Email 1010520053@qq.com,liqiuy@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-oUC113788102

Unique identifier UC113788102

Identifier etd-YangLiqiu-12547.pdf (filename)

Legacy Identifier etd-YangLiqiu-12547

Document Type Thesis

Format theses (aat)

Rights Yang, Liqiu

Internet Media Type application/pdf

Type texts

Source 20231213-usctheses-batch-1114 (batch), University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the author, as the original true and official version of the work, but does not grant the reader permission to use the work if the desired use is covered by copyright. It is the author, as rights holder, who must provide use permission if such use is covered by copyright.

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA

Repository Email cisadmin@lib.usc.edu