University of Southern California Dissertations and Theses

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

(USC Thesis Other)

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

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content MULTI-SCALE QUANTUM DYNAMICS AND MACHINE LEARNING FOR NEXT GENERATION ENERGY APPLICATIONS By Thomas M. Linker 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 (PHYSICS) May 2023 Copyright 2023 Thomas M. Linker ii Acknowledgements I would like give thanks first and for most to the endless love and support of my parents Jon and Karen and my sister Kay. Their endless advice and the many scientific dinner conversations we have had over the years has shaped me into scientist I am today. Second, I would like to thank my advisors Profs. Aiichiro Nakano and Priya Vashishta. Working with two advisors has been very rewarding in being able to see multiple perspectives in all aspects of my scientific work. Besides the technical training they provided, their advice in forming my skills in scientific presentation and communication has been invaluable, and will most likely serve me in my career far more than any technical skills I have obtained in my PH.D. I would also like to thank Prof. Ken-ichi Nomura, the lifeblood of the CACS group. Without him no project in CACS would be successful. In a similar vain I am also extremely grateful for our collaborators from Kumamoto Profs. Fuyuki Shimojo and Kohei Shimamura. Without them many of the projects at CACS would not exist. I would also like to thank Dr. Subodh Tiwari for their advice and training at the beginning of my PH.D. Their help and guidance allowed me to quickly hit the ground running and make progress. Next, I would like to thank my undergraduate advisor Prof. Matt Beekman. Working with him I was able travel to two international conferences and publish two papers as an undergraduate, which greatly inspired me to pursue this PH.D. I am extremely grateful for the opportunities afforded to me by working with him and the advice and training given. Finally I would like to thank all members of CACS group and wide berth of students across the PH.D program here who I have made close friends with over these past 5 years. From hiking mountains in San Gabriel to taking tequila shots and bowling in Koreatown their friendship made this PH.D an enjoyable one. Without them I doubt I would have finished. iii Table of Contents Acknowledgements ......................................................................................................................... ii List of Tables ................................................................................................................................... v List of Figures ................................................................................................................................ vi Abstract .......................................................................................................................................... xi Chapter 1 : Introduction ................................................................................................................. 1 1.1 General Overview ................................................................................................................. 1 1.2 Development of polymer dielectrics with high dielectric break-down fields for high field capacitor applications. ................................................................................................................. 4 1.3 Optical, mechanical, and electrical control of emergent polarization topologies in in next generation ferroelectrics .............................................................................................................. 6 1.3.1 Optical Control ............................................................................................................... 7 1.3.2 Mechanical Control ........................................................................................................ 9 1.3.3 Electrical Control ......................................................................................................... 13 1.4 Modeling of neutron scattering experiments for development of theoretical frameworks for studying hydrogen storage candidates. ...................................................................................... 15 Chapter 2 : Methods ...................................................................................................................... 20 2.1 Molecular Dynamics Method. ............................................................................................. 20 2.2 Ab-initio Molecular Dynamics ............................................................................................ 21 2.2.1 Derivation from Quantum Mechanics: ......................................................................... 21 2.2.2 Born-Oppenheimer(BO) Molecular Dynamics ............................................................ 26 2.2.3 Ehrenfest Molecular Dynamics .................................................................................... 26 2.2.3 Surface-Hopping Molecular Dynamics ........................................................................ 27 2.3 Density-Functional Theory(DFT) ....................................................................................... 28 2.3.1 DFT Derivation ............................................................................................................ 28 2.3.2 Exchange Correlation Functional ................................................................................. 31 2.2.3 DFT and Periodic Boundary Conditions (PBCs) ......................................................... 32 2.2.4 Plane Wave Basis and Pseudopotentials ...................................................................... 32 2.2.5 Real-Time TDDFT with Maxwells equations. ............................................................. 34 2.3 Neural Network Quantum Molecular Dynamics ................................................................ 35 2.3.1 NNQMD for ground state PTO .................................................................................... 35 2.3.2 NNQMD for Excited state PTO ................................................................................... 40 2.3.3 MM and NN/MM Model Development ....................................................................... 42 2.3.4 NNQMD for Model for NH3 Vibrational Properties ................................................... 43 2.4 Path Integral Molecular Dynamics ...................................................................................... 45 Chapter 3 : NAQMD simulations for development of polymer dielectrics with high dielectric break-down fields for high field capacitor applications ................................................................ 46 iv 3.1 Studying field dependance of hot carrier dynamics in polyethylene (PE): ........................ 46 3.1.1 Simulation Results ........................................................................................................ 46 3.1.2 Discussion .................................................................................................................... 55 3.2 Mitigating Hot Carrier Damage with hexagonal Boron-Nitride(hBN) Nano-Coats. ......... 55 3.2.1 Simulation and Experimental results ............................................................................ 55 3.2.2 Discussion .................................................................................................................... 63 3.2.3 Further Simulation and Experimental Details ............................................................. 64 Chapter 4 : Multi-Scale NNQMD simulations for optical, mechanical, and electrical control of emergent polarization topologies in in next generation ferroelectrics. ......................................... 67 4.1 Optical Induction of topological defects in PTO. ............................................................... 67 4.1.1 Simulation Results ........................................................................................................ 67 4.1.2 Discussion .................................................................................................................... 86 4.1.3 Further Simulation Details .......................................................................................... 87 4.2 Symmetry Guided Mechanical control of Polar Skyrmions with hybrid NN/MM simulations ................................................................................................................................ 88 4.2.1 Simmulation Results .................................................................................................... 88 4.2.2 Discussion .................................................................................................................. 101 4.3 Induction and Ferroelectric Switching of Flux Closure Domains in Strained PbTiO3 with Neural Network Quantum Molecular Dynamics .................................................................... 102 Chapter 5 : PIMD simulations for modeling of neutron scattering experiments for development of theoretical frameworks for studying hydrogen storage candidates. ........................................ 111 5.1 Simulation and Experiment Results .............................................................................. 111 5.2 Discussion ..................................................................................................................... 124 5.2.3 Further Simulation and Experiment Details ............................................................... 126 Chapter 6 : Conclusions .............................................................................................................. 129 References: .............................................................................................................................. 130 v List of Tables Table 1. Summary of generated training data. First 6 columns describe box dimensions and temperature of QMD simulations. The 7 th column shows the length of trajectory generated and 8 th column shows the number of frames used from each trajectory and rate of frame skipped, i.e., only using every other or every 5 th frame from the MD trajectory. Frames were skipped in order to avoid oversampling in regions of similar atomic coordinates. .......................................................................................................................................................... 38 Table 2. Computed and experimental lattice constants for PTO. Elastic constants Cij and bulk modulus K are in GPa unit. .............................................................................................................................................................. 39 Table 3. Excited state NNQMD Model Training Performance and Correlation Coefficient. ...................................... 40 Table 4. Optimal Lennard-Jones parameters sab. 𝜖 = 7.2x10 -2 (eV). ........................................................................... 42 Table 5. Experimental and computed lattice parameters and elastic constants. ........................................................... 42 Table 6. Average displacements of each atom type after 60 fs and 500 fs along the c polarization axis. Oxygen types are labeled O1, O2, O3 in accordance with their symmetry positions in the PTO unit cell which is pictured in Fig. 2A. .............................................................................................................................................. 71 Table 7. Total Number of Atoms and System Dimensions of XS-NNQMD simulation ............................................. 84 Table 8. Fitted High Energy N-H(D) stretching modes to double exponential fits. ................................................. 115 vi List of Figures Figure 1. Pseudopotentials. (a) Shows local component of PAW pseudopotential for W compared to all-electron wave-function. (b) Shows comparisons of pseudo-wavefunction and all-electron wave-function. .................... 33 Figure 2.NNQMD model development for PbTiO3. (a) NNQMD framework for energy, force and stress prediction. (b)-(d) Root mean square error (RSME) for energy, force, and stress tensor as a function of training epochs. (e)-(j) Comparison of partial radial distribution functions between NNQMD and DFT using SCAN exchange-correlation functional. .................................................................................................... 37 Figure 3. Training (black) and test (red) performance for the excited-state NNQMD model. (a) RMSEs of the system energy (A) and aggregated atomic force (B). .......................................................................................... 40 Figure 4. Comparison of excited state NNQMD and QMD radial distribution functions. Radial distribution functions g(r) of excited-state QMD (black) and its NNQMD model (red) up to 6 (Å) cutoff distance. The systems are thermalized at 300 K using NVT ensemble. Overall, QMD and NNQMD show excellent agreements in their peak position, height and width of g(r). Due to the photo-induced structural change, the fine structures originated from the tetragonality of PTO crystal disappeared. .................................................... 41 Figure 5. Comparison of Excited State NNQMD and QMD bond angle distributions. Bond angle distribution functions of ground state QMD (black) and its NNQMD model (red). O-Pb-O, O-Ti-O and O-O-O angles distribution are shown with bond cutoff distance of 3Å. ..................................................................................... 41 Figure 6. Polyethylene (PE) under electric field. (a) (001) (or xy)-plane view of the PE slab, where PE chains are numbered from 1 to 36 with respect to their positions along the [100] direction. (b) Partial density of states (PDOS) projected onto the 36 PE chains under an electric field of 100 MV/m. Schematic of excitation is shown, where blue and red circles represent hole and electron, respectively. Blue and red arrows show the respective direction of travel of hole and electron. (c) The response of the local Kohn-Sham (KS) potential to the applied electric field. Blue and red curves show the KS potential profile along the [100] direction with and without electric field, respectively, whereas the black curve shows their difference. ........... 47 Figure 7. Hole relaxation and localization. (a)-(c): Time evolution of KS eigenvalues under electric fields 100 (a), 600 (b) and 900 (c) MV/m, respectively. The blue and red curves are the hole and electron eigenenergies, respectively. (d)-(f): Center-of-mass (COM) position of the hole, 𝑋GHole, along the [100] direction for 100 (d), 600 (e) and 900 (f) MV/m, respectively. (g)-(i) Chain participation number, CP, of the hole for 100 (g), 600 (h) and 900 (i) MV/m, respectively. ............................................................................ 49 Figure 8. Field Induced Localization. (a)-(c) show the expectation value of the position VBM wavefunction along the [100] or x direction during the first 50 fs for fields of 300,600, and 900 MV/m respectively. (d) shows the energy fluctuation of the VBM wavefunction versus the electric field. ............................................. 52 Figure 9. Bond stretching and breaking. (a) Three distinct bonds colored black, red, and blue present in PE chains. Dark gray and white spheres represent carbon (C) and hydrogen (H), respectively. (b)-(f): Time evolution of C-C bond lengths for surface chains 35(solid curve) and 36 (dashed curve). Black, red and blue curves represent same colored bond in (a). .................................................................................................. 54 Figure 10. (a) Surface potential decay and (b) trap energy, and (c) Weibull breakdown strength of Kapton and hBN coated Kapton. ............................................................................................................................................. 56 Figure 11 .Optimized structures and local potential profiles. (a) Optimized Al, hBN and PE structures used to compute the local potential and band alignment. Al, B, N, C and H atoms are colored in silver, blue, yellow, gray and white, respectively. (b)-(d) Laterally (xy) averaged local potential profiles along the surface- normal (z) direction. ............................................................................................................................................. 58 vii Figure 12. Band alignment and schematic of charg- injection experiment. (a) Computed band alignment with the PBE functional for Al, hBN and PE with VBM, CBM, Fermi-energy levels labeled. (b) Schematic of breakdown experiment, where hot carriers can either be: (I) directly injected into the polymer at electrode interface, or (II) caused by excitation of secondary electron-hole pairs in the bulk of the polymer. .................. 60 Figure 13. NAQMD simulation. (a) hBN-PE interface used in the NAQMD simulation. B, N, C and H atoms are colored mauve, lime, black and white, respectively. Initial excited electron and hole Kohn-Sham (KS) wave functions are colored red and blue, respectively. Electric field was applied along the interface-normal direction (z) using a sawtooth potential. (b) Laterally (xy-plane) averaged local potential with (red curve) and without (blue curve) the applied electric field. The difference displaying the linear voltage drop is colored in black. (c) Time evolution of the KS eigenvalues with the hole and electron energies colored in red and blue respectively. (d) Semi-log plot of the excitation energy vs. time (open circles) to estimate the relaxation rate by linear fit (dotted line). ............................................................................................................. 61 Figure 14. Carrier participation dynamics. (a-g) Angular-momentum (l)-dependent fractional contribution of each atomic orbital type to the hole wave function. The hole wave function hops between the C(p) and H(s) PE states and B(p) and N(p) hBN states before localizing to hBN. (h-n) l-dependent fractional contribution of each atomic orbital type to the excited electron wave function, which remains primarily localized in hBN for the entire simulation. ...................................................................................................................................... 63 Figure 15. Vector potential and absorbed energy. Time evolutions of the external vector potential due to a 200 nm laser pulse (A) and energy absorbed by the electronic system (B). ............................................................... 68 Figure 16. Electronic thermalization. Occupation numbers of KS orbitals as a function of KS energies at time t = 0.7 fs (A) and 2.4 fs (B). ................................................................................................................................... 69 Figure 17. Excitation as a function of fluence. Fraction of valence electrons that are excited as a function of laser fluence. ................................................................................................................................................................. 70 Figure 18. Photoinduced charge transfer. (A) Atom projected partial density of states (PDOS). States near the valence band maximum (VBM) are primarily composed of O and states near the conduction band minimum (CBM) are primarily composed of Ti, indicating expected O-Ti charge transfer upon excitation. (B) Isosurfaces of hole (blue) and excited electron (red) densities upon excitation. The hole density is primarily localized around O atoms, while the excited electron density is localized around Ti atoms. Spheres with black, silver and red color represent Pb, Ti and O atoms, respectively. (C) Instantaneous force vectors as a result of excitation. ............................................................................................................................................... 70 Figure 19. Bulk lattice and polarization dynamics. (A) Schematic of the initially activated TO1 phonon mode and the up-converted TO2 phonon mode. Oxygen atoms are labeled O1, O2, and O3 with respect to their symmetry position in the PTO unit cell. (B) Snapshot of oxygen rotations with magenta arrow labeling of the rotation orientation for 4 unit cells in ab plane 940fs after excitation. (C) ground and (D) excited state vibrational density of states (VDOS). VDOS was computed from Fourier transform the velocity auto- correlation function of the MD simulation at 300 K. (E) Ground-state and (F) excited-state phonon dispersion. Activated high-frequency optical branches are both softened and hardened to form a phonon band gap. (G) Excited-state polarization dynamics at 300 K. Average Ti polar displacements along the c axis are initially reversed, followed by a steady-state structure with near-zero average polarization. ................ 72 Figure 20. Ti polar displacement dynamics. (A) ground state and (B) excited state polarization dynamics at room temperature. Cyan line marks zero polarization. ................................................................................................. 73 Figure 21. Lattice evolution under excitation in NPT ensemble. (A) Evolution of the lattice constant ratios. (B) Evolution of the average unit-cell volume for the 4×4×4 supercell. Contraction of the lattice is consistent with the observed activation of the A1 TO1 and E TO2 modes. ......................................................................... 74 viii Figure 22. Mulliken bond-overlap population analysis. (A) and (B) show the Pb-O and Ti-O bond overlap populations, respectively, both in the ground and excited states. ........................................................................ 74 Figure 23. Oxygen charge dynamics. (A) O charge dynamics in the ground state within NVT ensemble. (B) and (C) O charge dynamics in the excited state within NPT and NVT ensembles, respectively. Shading represents standard deviation. .............................................................................................................................. 76 Figure 24. NAQMD simulation of PTO. (A) PTO unit cell with O types labeled. Spheres with black, silver and red colors represent Pb, Ti and O atoms, respectively. (B) Polarization dynamics for given unit cell. (C) Excited electron participation dynamics on Ti atom in given unit cell. (D) Hole participation dynamics on O atoms in given unit-cell with coloring given by O atom type as labeled in (a). The O1(green) and O2(magenta) curves are overlapping. .................................................................................................................. 77 Figure 25. Dual order parameters of the phase transition. (A) Ti-O partial radial distribution g(r) for both excited- state and ground-state trajectories. In the ground state, polarized structure Ti-O symmetry is broken, resulting in splitting of the first g(r) peak. The symmetry is restored upon optical excitation as the polarization is lost, resulting in a single peak. (B) O-O-O bond angle distribution for ground and excited states. Third peak splits in the excited state due to symmetry breaking from the optical induction of oxygen rotations. (C) and (D) diagram the phase change in terms of the ferroelectric (FE) order parameter 𝐷𝑧 and the tilt order parameter 𝜙𝑡𝑖𝑙𝑡. Optical excitation lowers the FE energy barrier and raises the tilt energy barrier ................................................................................................................................................................... 79 Figure 26. Excited-state topological defect dynamics. (A) XS-NNQMD workflow. (B) Initial configurations of bulk PTO crystal. Local polarization of each TiO6 structure is color-coded by their magnitude along c-axis. Panels (C)-(D) show the photoinduced nucleation of TiO6 tilt domains and the formation of polar stripes in the bulk PTO system at 24 picoseconds using XS-NNQMD simulation. The polar stripes are sandwiched by domains with different TiO6 cage tilt orientations. Panel (E) shows the initial configuration of PTO crystal with circular anti-polarization domains. The circular anti-polarization domains were first introduced by displacing Ti atom position and relaxed at low temperature using ground state NNQMD simulation. Emergent Block-type Skyrmion at each domain are shown in inset of panel (E) with the polarization vectors on a-b plane prior to optical excitation. Panels (F) and (G) show the nucleation, propagation, and formation of polar stripes in the system using XS-NNQMD simulation. Unlike in the bulk PTO crystal system, the polar stripes consist of multiple segments of remnant polarization of the original PTO crystal and the circular anti-polar domains. (H) Photoinduced domains and domain boundaries color-coded by the tilt orientation Ti-O6 cage as order parameter. (I) Atomic coordinates around the domain boundaries highlighting TiO6 cage tilt orientations in red. For clarity, only O-atom positions and lines connecting neighboring O atoms are shown. ......................................................................................................................... 81 Figure 27. Size-dependence of domain boundary dynamics. (A) and (B) show the domain dynamics for system 1 at 2.5ps and 7.5 ps respectively. Unit-cells are color-coded by the order parameters based on the tilt angle of TiO6 cage. The original PTO structure with zero tilt angle is shown blue. Red and green show positive and negative tilt angles due to the light-induced symmetry breaking described in the main text. (C) and (D) show the domain dynamics for system 2 and (E) and (F) show the dynamics system 3 again at 2.5ps and 7.5ps. (g) shows a zoom in on the domain structures of system 3. (H) shows the computed structure factor S(Q) at small Q for system 1 and system 3. System 1 is unable to capture the small Q region of S(Q) due to inadequate system size. .................................................................................................................................... 85 Figure 28. NN/MM model of STO/PTO super lattice ground state structure. (a) NNQMD framework used for NN engine. Atomic coordinates are used to construct symmetry functions which are fed into neural-network to predict energy, force, and stresses (b) STO-PTO heterostack of 8×16×8 STO-PTO-STO unit cells with periodic boundary conditions. Rc and RB represent the core and buffer radius to linearly mix the NNQMD and MM interactions at the boundary. (c) Represents the initial condition of circular oppositely polarized domain for STO/PTO super lattice. (d) Top perspective view of STO/PTO skyrmion structure after relaxation. (e) Front perspective view of the skyrmion after relaxation. (f) (001) Plane view of polarization ix vectors with background density colored by interpolated [001] polar displacements of Ti atoms Dz. (g) (001) Plane view of polarization vectors with background colored by the topological charge density q. .......... 89 Figure 29. Skyrmion evolution under uniform compression. (a) and (b) Top slice of skyrmion polarization at 0.0% and 1.0% compression along [001]. Blue double arrows are given as a guide to the eye. (c) and (d) (010) plane slice of Skyrmion structure at 0.0% and 2.7% strain. (e-f) Local c lattice constant as a function of position along the interface after compression and evolution of the compressed structure. (g)-(h) Local pressure profiles for the heterostack after compression and evolution of the compressed structure. (i) and (j) Evolution of polarization and local pressure under 3.0% compression. ......................................................... 92 Figure 30. Polarization dynamics under isotropic compression. (a) Snapshots of polarization structure before and after compression. (b) (001) plane perspective view of unit cells colored by [001] Ti polar displacements after compression. (c) and (d) (001) plane perspective slice of unit cells colored by in-plane ([010] and [100]) Ti polar displacements after compression. (e) Snapshots of polarization texture before and after unloading. (f) A schematic of ideal skyrmionium. (g) Skyrmionium formation visualized by Ti displacement vector on each PTO unit cells. ....................................................................................................... 95 Figure 31. Compressive skyrmionium-to-skyrmionium phase transition. (a) Compressive induction of anti- phased octahedral rotations under compression in the bulk of PTO structure. The rotation of Ti-O6 octahedra is observed through splitting O-O-O bond angle distribution which defines their relative orientation. (c) Polarization in the [001] direction after unloading. Concentric circular domain seeds creation of skyrmions with opposite topological charges that cancel out. (c) Topological charge density, which integrates to zero over the skyrmionium surface. ..................................................................................... 97 Figure 32. Strain wave dynamics after isotropic compression. (a) Pressure evolution under isotropic loading of the skyrmion system. Isotropic loading results in high pressure state with a pressure gradient within the PTO layers. This pressure gradient is slowly dissipated after compression resulting in uniform pressure for the PTO layer, and decline in pressure for the STO layer. (b) During the dissipation process the local pressure exhibits wave like dynamics within the middle of skyrmion in (001) plane. These wave like dynamics ultimately result in reversal of the polarization in some domains. (c) Schematic of strain wave induced polarization reversal. (d) Local z-strain profiles superimposed with polarization vectors. ................... 98 Figure 33. Skyrmionium formation at 200 K and 300 K. (a) Snapshots of polarization texture before loading at 200K. (b) Polarization texture under compression at 200 K. (c) Polarization map after unloading forming skyrmionium structure at 200 K. (d-f) Similar Skyrmionium formation seen at 300 K. ................................... 100 Figure 34. Switching Dynamics in 180 Degree Domain Wall(DW). (a) Initial domain wall structure at 300K. (b) and (c) show propagation DW switching 1ps and 5ps after turning on the electric field respectively. (d)-(e) show progression of switching of a given layer at the DW boundary. Switching starts with in an initial flip at the DW followed quickly by a 1D flip to the strong Ti-Ti repulsion of the facing polarized domains. This followed by 2D propagation in the layer. (i)-(k) show total polarization dynamics as a function of time for three different field strengths. ............................................................................................................................ 103 Figure 35. Evolution of polarization dynamics in 40x40x40 PTO unit-cell under applied electric field of 3MV/cm. Red color is polarization of unit-cells along the field direction and blue against. ............................ 105 Figure 36. Polarization rotation during switching. (a) Shows dynamics during switching process for a given unit- cell that resemble that of a Neele like rotation, while (b) shows more Bloch like rotation. The DW wall is along the x axis, and the polarization axis is the z axis in this figure. (c) Shows distribution of Ti polar displacements right at the point of switching (𝐷𝑧=0Å). ................................................................................. 106 Figure 37. Polarization dynamics in strained PTO. (a) shows ab plane lattice structure of PTO prior to application of the electric field and (b) shows plastically deformed ab plane after application of the electric field. (c) Shows the a displacement vector field of the Pb at atoms from there postions after application of the electric field. The atoms seems form Bloch like rotations as a result of the applied field. (d) Shows the polarization x profile after the application of the applied field. (e) and (f) shows the local strain profile after the applied field in the a/x and b/y direction respectively. ................................................................................................... 107 Figure 38. Switching of flux closure domain. (a) Potential energy profile during the switching process(with including the field energy). (b)-(d) Polarization profiles at the beginning, middle, and end of switching. At the height of the barrier a vortex state forms which is the intermediate state. .................................................. 109 Figure 39. Dynamic Structure Factor and Vibrational Spectrum for Solid Ammonia. (a) and (b) integrated dynamic structure factor S(Q,E) measured with SEQUOIA spectrometer for incident energies of 30 and 70 meV at various temperatures for ND3 (c) Integrated S(Q,E) measured with VISON spectrometer for NH3. (d)-(f) Full S(Q,E) spectrum measured for ND3 with incident energy 30 meV with increasing temperature from left to right. (g)-(i) Full S(Q,E) measured for ND3 with incident energy 70 meV with increasing temperature from left to right. ............................................................................................................................ 113 Figure 40. Intramolecular modes for ammonia. (a) and (b) illustrate molecular modes as a function of temperature measured for ND3. (c) Illustrate a double exponential fit high energy N-D stretching modes. (d) and (e) illustrate molecular modes as a function of temperature measured for NH3. (f) Double exponential fit for high energy N-H stretching modes, where shift in liquid phase becomes clearer. ............................................ 114 Figure 41. Computed vibrational spectrum for solid ammonia. (a) Computed INS spectrum for solid NH3 with various dispersion interactions in comparison to Vision experiment. Significant discrepancies are observed, especially for the librational band (as marked by the arrow with an illustration of a librational mode in solid ammonia). (b) Computed S(Q,E) with various exchange-correlation functionals and dispersion interactions in comparison to SEQUOIA experiment. Comparison of DFT simulations of ND3 spectrums is better than with NH3 spectrums indicating the role of nuclear quantum effects. ................................................................ 116 Figure 42. Evaluation of nuclear quantum effects with DFT. (a) INS spectra of solid ammonia simulated with various step sizes using the finite displacement methods (e.g., FDM-01 means a step size of 0.1 Å). The drifting librational band is an indication of its anharmonicity. (b) Potential energy profile of NH3 rotation in solid ammonia, solved by nudged elastic band method. The energy barrier is determined to be 167 meV and the associated NH3 quantum rotor has an excitation energy (n = 0→1) of 32 meV. (c) Illustration of PIMD simulation. While typical ab initio simulation treats atoms classically while electrons quantum- mechanically, PIMD simulation uses multiple replicas of each atom to mimic nuclear quantum effect. (d) INS spectra of solid ammonia from DFPT, ab initio PIMD (1-bead and 32-bead TRPMD) simulations, and VISION experiment. The PIMD simulation with 32 beads reproduces the experimentally observed peak positions very well, whereas the other two exhibit major discrepancies. .......................................................... 119 Figure 43. TRPMD simulations of ammonia using NNQMD. (a) Allegro model for NNQMD. (b) Comparison of computed vibrational spectrum using Allegro TRPMD and VISION spectrum which shows excellent agreement. DFT results at 0 K are shown for comparison. (c) Comparison of classical and TRPMD simulation of high-energy N-H stretching mode. Merging of Symmetric and Anti-symmetric splitting agrees with SEQUOIA experiment. (d) INS computation of N-H stretching mode with and without inclusion of mutli-phonon effects. (e) TRPMD ND3 VDOS for high-energy N-D stretch where peak splitting remains. (f) Computation of VDOS for NH3 liquid state, in which low energy peaks are all washed out. (g) Computed hardening of NH3 vibrational modes in liquid state consistent with SEQUOIA experiment. ........................................................................................................................................................ 123 xi Abstract First principles atomistic modeling of materials from quantum mechanics allows accurate prediction of material properties, but are seriously limited in their scope of applicability by their inherent computational complexity. In practice approximations of varying degrees are used to simulate any material of relevance. Through tackling three problems in the field of energy and materials science I will climb and connect multiple levels on the quantum hierarchal latter using multi-scale quantum dynamical simulations incorporating multiple levels of theory from real time time-dependent density functional theory (RT-TDDFT) for small unit-cells consisting of a few atoms to billion-atom molecular dynamics simulations for defects and disorder. During this climb I will demonstrate how incorporation of machine learning in conjunction with quantum simulations can help us. I will demonstrate how we can use quantum simulations to discover key features of materials electronic structure for inputs in machine-learning based material design frameworks as well as how we can use machine learning accelerate computation of quantum simulations. In particular, I will discuss the utilization of these methods for three problems: (1) Development of polymer dielectrics with high dielectric break-down fields, (2) Optical, mechanical and electrical control of polarization topologies in in next generation ferroelectrics, and (3) Modeling of neutron scattering experiments for studying hydrogen storage candidates. While seemingly quite different in nature, these problems require truly dynamical simulation methods that include multiple levels of physics. The multiscale quantum simulations presented here highlight how we can incorporate quantum mechanics and machine-learning for material discovery and design. 1 Chapter 1 : Introduction 1.1 General Overview The idea of modeling the matter around us from budling blocks stems from the ancient Greeks, with the word atom being derived from the Greek word atomos—uncuttable. While much has developed since then with the advent of the theories of quantum-mechanics, chemical bonding, etc., a large amount of scientific research being performed today stems from trying to understand matter around us from atoms and their fundamental interactions with each other. One of the most powerful tools in this regard, is the molecular dynamics (MD) method, with the first molecular dynamics simulation being performed by Aneesur Rahman in 1964 to study the dynamics of solid and liquid Argon 1 . Molecular dynamics is a method of atomistic simulation that treats atoms as classical point objects and follows the trajectory of all atoms within a given interaction potential energy surface by integrating Newton’s equations of motion 2 . Clearly this type of method is reliant on the accuracy of the interaction potential and that the nuclear quantum effects (NQEs) are minimal. The most accurate of these methods is quantum molecular dynamics 2 (QMD), which is an extension of MD that computes the potential energy surface from first principles by solving either the time-dependent or time-independent Schrödinger equation on the fly during the MD simulation. While one would like to be able to solve the Schrödinger equation directly, this is computationally intractable even with the most advanced super-computers today. Hence, most QMD simulations are based off Kohn-Sham density functional theory 2,3 (KS-DFT). DFT maps the intractable exponential scaling problem of solving the Schrödinger equation to solving the 2 O(N 3 ) scaling Kohn-Sham equations by taking advantage of the Hohenberg-Kohn 4 (HK) and Runge-Gross 5 (RG) theorems, which state that any observable can be written as an unique functional of the electron density. In addition, the majority of first principles methods assert the Born-Oppenheimer approximation which separates the nuclear and electronic motion, as the electronic time-scales are much shorter than that of the nuclei in chemical systems. One can see this trivially from the fact that the mass of the proton is about 1000 times heavier than that of the electron. While this works well for ground states, in excited state dynamics non-adiabatic coupling between the nuclear degrees and electronic degrees of freedom occurs and this approximation breaks down. In this case more advanced methods, and correspondingly more computationally expensive methods, that incorporate electronic dynamics and there coupling with nuclear degrees of freedom such as Ehrenfest and the surface hopping approaches based on formulations of time- dependent density functional theory or other time-dependent quantum-mechanical methods are required. Furthermore, these methods all assume classical nuclei which breaks down in the presence of hydrogen rich materials, high pressure, and low temperature dynamics. This too can be addressed with molecular dynamics through the so-called path-integral molecular dynamics (PIMD) approach, where one samples quantum partition function of the atomic system via multiple replica molecular dynamics simulations that are coupled together. This implementation can be on the order of 10-1000 times more expensive than underlying type of molecular dynamics methods depending on the method of implementation for PIMD. One can quickly see a zoo of quantum-mechanical approaches are forming, within which figuring how to navigate with one’s computational resources is quite daunting. For many 3 problems multiple scales of physics are required and one must take care when to use which approximations. In this work through addressing three problems in materials physics and energy science : (1) Development of polymer dielectrics with high dielectric break-down fields for high field capacitor applications, (2) Optical, mechanical and electrical control of emergent polarization topologies in in next generation ferroelectrics, and (3) Modeling of neutron scattering experiments for development of theoretical frameworks for studying hydrogen storage candidates, I will demonstrate how we can incorporate a powerful new tool, machine learning with quantum dynamical simulations to help us. In this regard in the first application I will demonstrate how we can use non-adiabatic quantum molecular dynamics (NAQMD) simulations to find key features in the electronic structure of polymers that are essential for understanding their here breakdown processes which can be used as proxies/features in machine-learning frameworks to predict better dielectric polymers. In the second application I will utilize multi- scale neural-network quantum molecular dynamics (NNQMD) simulations to understand optical, mechanical and electrical control of next generation ferroelectric materials. Finally in the third application I will discuss combining neutron-scattering experiments, Density-Functional Theory (DFT), and NNQMD based PIMD simulations to understand the role of nuclear quantum effects in ammonia vibrational dynamics. In the remainder of this chapter, I will provide a detailed introduction to these three major problems and how I have utilized multiple simulation frameworks to address them. In Chapter 2 I will give an overview of the simulation methods utilized in this thesis. In Chapters 3-5 I will discuss in detail the results for the 3 problems laid out in the remainder of this introduction. Finally in Chapter 6 I will give concluding remarks. 4 1.2 Development of polymer dielectrics with high dielectric break-down fields for high field capacitor applications. The use of polymers as dielectric materials in electronic circuits offer many advantages over typically used inorganic materials as they are highly flexible and are able to be cheaply produced 6–8 , but their application is severely limited by their electronic breakdown in the presence of high electric fields 9,10 . For inorganic materials, dielectric breakdown under electric field is well described by electron avalanche theory 11–13 . Formulated by the early works of Zener 14 , Frölich 15 and von Hippel 16 , avalanche theory describes the onset of electrical breakdown when the applied field accelerates electrons past the impact ionization energy, which in turn causes further excitation of hot carriers (electrons or holes). This avalanche process then results in structural damage of the material such as crack formation 8,9 . In this model, the critical field for breakdown can be determined by solving the energy balance equation between the energy gain by the applied field and the energy loss due to electron-phonon scattering 11–13 . Breakdown in organic materials like polymers is expected to follow an avalanche process similar to inorganic materials, but no such theory exists to predict the critical field. The difficulty partly arises from the strong coupling between charge carriers and atomic motions (or phonons) in polymers. In the case of inorganic solids, weak carrier-phonon coupling allows the use of perturbation theory to estimate the critical field 12,13,17,18 . In polymers, in contrast, strong carrier-phonon interaction leads to polaronic effects 19 , thus precluding the use of the well-established perturbation theory for critical-field estimation. Consequently, state-of-the-art computational screening for rational design of high-performance dielectric polymers heavily relies on a qualitative proxy (i.e. electronic band gap) of the breakdown field 7 . For quantitative estimation of the breakdown field in polymers, which is critically missing, it is essential to mechanistically understand field-induced hot carrier dynamics. In this 5 regard first-principles simulations based on nonadiabatic quantum molecular dynamics (NAQMD) that allow for modeling of strong carrier-phonon interaction can be used in order to develop a quantitative theory of breakdown of dielectric polymers. Previous NAQMD work had shown hot carriers can localize at polymer surfaces, thereby causing defect formation and carrier multiplication in a prototypical dielectric polymer, polyethene (PE) 20 . In chapter 3, I further show how near the breakdown field, that hot carriers can strongly couple to C-H vibrational modes in PE, resulting in rapid localization of hot carriers at the surface and damage to C-C bonds. This polaronic localization transition is analogous to the coupling of vibrational resonance with Stark shifts of electronic energy levels under electric field. 21 In addition, I use NAQMD to model how to use nano-coats to limit hot-carrier induced breakdown. To limit the effect of injected hot carriers, there has been a growing effort to utilize nanoscale two-dimensional (2D) materials with large bandgaps and favorable electrically insulating capability as coating materials. In the past decade, 2D materials with high breakdown fields, such as hexagonal boron nitride (hBN) (𝐸 ! ~1200 MV/m) 22 , have been widely used to perform nanocoating on polymer dielectrics for significantly enhanced breakdown strength 23–25 . The current hypothesis is that insulating 2D coatings prevent injection of hot charge carriers from the electrode, thus suppressing the leakage current. However, there is critical lack of mechanistic understanding about the behavior of the charge carriers at the electrode-polymer- coating interface nor a predictive framework to determine good potential nano-coating materials. I investigate the dynamics of hot carrier dynamics in hBN-polymer interfaces with hBN as a nano-coat along the charge injection interface. Collaborators at Uconn showed with surface- potential-decay and dielectric-breakdown measurements of the commercial polymer Kapton coated hBN could induce carrier trapping effect that leads to an enhanced break-down field. To 6 investigate the microscopic mechanisms behind this enhancement, NAQMD simulation was performed on a hBN-polymer electrode interface. Simulation reveals the formation of a deep potential well in the hBN coat, which can trap hot carriers, potentially mitigating their damage. The application of a hBN coat also results in a substantial decrease in the carrier relaxation rate in comparison to previous studies, which is indicative of weak carrier-phonon interaction. This weak carrier-phonon interaction suggests the potential of hBN to mitigate the detrimental polaronic effects seen to enhance the hot-carrier damage 26 . In addition, the signature deep well formed by the hBN layers can readily be calculated with a standard DFT formalism. This key feature can thus be used as a qualitative proxy for designing a rational search of effective coating materials. 3,8,20,21 1.3 Optical, mechanical, and electrical control of emergent polarization topologies in in next generation ferroelectrics The formation of topological defects originally seen in particle physics such as skyrmions and merons 30–32 in the spin fields of ferromagnetic materials has garnered attention over the past two decades in the field of spintronics for development next generation information transfer and memory devices such as race-track memory 33 . Recently, there has been growing effort to synthesize the electric analogue of these topological structures in ferroelectric materials with the most prominent example being the synthesis of polar vortex states and skyrmion bubbles in paraelectric/ferroelectric SrTiO3 (STO)/PbTiO3 (PTO) based nano-superlattices 34–39 . The emergence of such topological structures in these heterostacks relies on a subtle balance of the strain gradients and polarization fields within the materials, which make them ideal for exploring low power next generation devices, where optical, electrical and mechanical manipulation of these strain and polarization fields switches structures that are topologically protected from 7 thermal noise 40,41 . In chapter 4, I will address all three of these methods of control using multi- scale neural-network molecular dynamics (NNQMD) simulations. Introductions to these sub- topics and the methods used are described in the following three sub-sections. 1.3.1 Optical Control Optical control via ultrafast photoexcitation offers an exciting new avenue to control materials by drastically modifying their potential energy surface through the creation of electron- hole pairs, which can allow for order of magnitude faster phase changes in comparison to traditional methods as well as access to hidden non-equilibrium phases 38,42–46 . As mentioned above the perovskite oxide based nanostructures are suitable system of interest for study of via optical control due to their wide range of desirable electronic properties from traditional ferroelectricity to the formation of complex magnetic and polar topological structures such as skyrmions and merons 47–52 . These properties are usually dynamically controlled via traditional phase-change inducing methods such as heating and straining 50,51,53–58 , which has made perovskite oxides ideal materials for nanostructured devices 51,59 . Recent studies in non-equilibrium optical control of perovskite oxides have primarily focused on using weak THz pulses to couple with optical phonon modes to directly control polarization patterns 60–64 . Above-band gap laser pulse excitation with moderate fluence has also been shown to control nano-polarization topologies in SrTiO3 (STO)-PbTiO3 (PTO) superlattices 38,65 through introduction of electron-hole screening at the material interface. However, topological polarization control in light-induced far-from-equilibrium conditions has been less explored. This is largely due to that modeling the rich topological changes induced by non-equilibrium electron-ion dynamics under strong above band-gap excitation requires understanding of multiple spatiotemporal scales, which has remained elusive. 8 To investigate far-from-equilibrium optical control of complex polar topologies, I have implemented a multi-scale excited-state neural network quantum molecular dynamics framework (XS-NNQMD) to study strong optical excitation on the prototypical ferroelectric perovskite oxide PTO. In our multi-scale framework, excited electron-lattice dynamics is studied by quantum molecular dynamics (QMD) simulations that combine Fermi-occupation quantum molecular dynamics (FOQMD) to describe the effects of massive electronic excitations on atomic motions 66 , real-time time-dependent density functional theory (RT-TDDFT) along with Maxwell’s equations to describe light-matter interaction 67 and provide the initial excitation level for FOQMD simulation, and surface hopping based non-adiabatic quantum molecular dynamics (NAQMD) to include non-adiabatic coupling between electrons and phonons 68 , thereby assessing the relaxation time between the electronic and lattice temperatures until which FOQMD is justified. Ground- state QMD and excited FOQMD simulations were then used to train neural networks to perform large scale XS-NNQMD simulation to investigate atomistic mechanisms accountable for far-from- equilibrium control of large nano-polarization topologies in PTO. Despite remarkable success of machine-learning based MD simulations realizing atomistic simulations at scale while retaining quantum-mechanical accuracy 69 , NNQMD incorporating massive electronic excitations has been less explored. I find that strong photoexcitation lowers the ferroelectric energy barrier, which results in a transverse optical (TO) phonon driven transition to a non-polar symmetric structure at room temperature. Concomitant to the loss of ferroelectric order under far-from-equilibrium excitation, a new order emerges in the form of tilting of the oxygen octahedra. Such emergent order plays a key role in a wide variety of polar topological structures. Billion-atom XS-NNQMD simulations demonstrates how different regions in the crystal acquire different tilting orientations, which in 9 turn result in frustration of the tilting phase to create topological defects in the form of string-like domain boundaries. Along these frustrated regions, the ferroelectric order is preserved and forms polar strings. This type of phase transition is analogous to non-adiabatic Kibble-Zurek phase transitions, where rapid temperature quenching in a uniformly symmetric field introduces string- like domain boundaries in space 70,71 , but with two key distinctions. First, non-adiabatic transition here is induced through photoexcitation-driven change of energy landscape rather than conventional temperature quench. Second, the dynamics resembles more that of an inverse Kibble Zurek transition 72 , where the system is driven from an equilibrium anti-symmetric polarized state to a highly nonequilibrium non-polar symmetric state with topological defects emerging as result of the emergence of a hidden order parameter in the form octahedral tilting. Overall, our multi- scale quantum-simulation and machine-learning framework reveals pattern/defect formation during highly non-equilibrium phase changes (30-32) applied to perovskite-based ferroelectrics, but in the less explored regime of light-induced electronic excitation. The uncovered hidden parameter of light-induced octahedral rotations demonstrates the utility of the multi-scale framework for exploring ultrafast control of ferroelectric materials for applications within the emerging field of ferroelectric “topotronics” 40,41,73 . 1.3.2 Mechanical Control As stated prior, the emergence of such topological structures in these heterostacks relies on a subtle balance of the strain gradients and polarization fields within the materials, which make them ideal for exploring low power next generation devices, where optical, electrical and mechanical manipulation of these strain and polarization fields switches structures that are topologically protected from thermal noise 40,41 . The majority of efforts in topological manipulation of these vortex and skyrmion states have focused on direct manipulation of the 10 charge and polarization boundary conditions that give rise to skyrmion or vortex formation with static electric fields and ultrafast laser pulses, where topologically based switching has been experimentally and theoretically demonstrated 35,36,38,74–76 . While mechanical writing of ferroelectric polarization has been demonstrated 56 , direct mechanical manipulation of polar topological structures has been less explored, with primary efforts focused on strain engineering through modification of sample temperature, thickness and substrate lattice mismatch 51,77–79 . In particular, little is known about the stability and phase changes of these topological defects to externally applied mechanical stresses. As a wide variety of loading modes and symmetries can be adopted for mechanical stimuli that can directly control strain fields within these nanostructures, it is most compelling to investigate mechanical control of polar topology in order to establish a fundamental science of polar topological switching. Rational mechanical control of polar topologies requires fundamental understanding of the switching dynamics on picosecond-to-nanosecond time-scales on 10s-to-100s of nanometers spatial scales 40 . In particular, it is essential to describe the dynamical response of the strain gradients to mechanical manipulation and the corresponding response of the polarization through static and dynamical flexoelectric effects 80,81 . In this regard, simulation methods utilizing atomistic level dynamics offer an appropriate theoretical framework; however, one of the main difficulties in exploring such phase changes atomisticaly is the prohibiting computational cost in modeling the polarization dynamics within the ferroelectric layer with quantum-mechanical (QM) accuracy on the spatiotemporal scales that these topologies emerge while maintaining the appropriate boundary conditions applied by the paraelectric layer. To address this challenge, I employ a hybrid neural network 69 (NN) and Molecular Mechanics(MM) framework. My framework draws from the hugely successful field of physics-informed machine learning (ML) 82 11 and expands it to physics-ML hybridization in the spirit of the celebrated multiscale simulation approach called QM/MM which received the 2013 Nobel prize in chemistry 83,84 . In my NN/MM approach, a highly scalable neural network quantum molecular dynamics (NNQMD) 76,85–87 force field is trained using ab-intio molecular dynamics (MD) simulation data to accurately describe complex polarization dynamics within the ferroelectric layer, which in turn is embedded in a classical force field that correctly applies the essential strain and null-polarization boundary conditions imposed by the paraelectric layer. I applied the NN/MM framework to the study of mechanical manipulation of polar skyrmions in STO/PTO nanostructures under compressive strain. Individual MM and NN force fields were found to produce mechanical properties in good agreement for measured bulk STO and PTO, and the NN/MM model successfully describes the formation of polar skyrmions in STO/PTO nanostructures consistent with experiment and second principles and phase field simulations 34,35 . I examined the dynamics of these skyrmions in two strain modes with different symmetries: (1) uniaxial compression along the c polarization axis to mimic uniaxial loading along the polarization direction; and (2) isotropic compression to examine the effects of symmetry-preserving deformation of the skyrmion system. Under the symmetry-breaking uniaxial compression, I find a squishing-to-annihilation phase transition marked by interfacial depolarization and rotation of the TiO6 octahedra of PTO near the interface. The depolarization at the interface I determine to be due the compliance of the local c lattice constant to STO due to its stiffer C11 elastic constant, which eventually results in complete loss of polarization in the heterostack within picoseconds. Such mechanically- induced rapid depolarization of skyrmion structure could be of interest for pulsed-power technologies through the force-electric effect 88,89 . 12 Under the symmetry-preserving isotropic compression, on the other hand, polarization reduction occurs but annihilation is not seen. Skyrmion domains can become punctured, leading to the formation of concentric skyrmions with opposite topological charges upon unloading of the pressure. The combination of the two oppositely charged skyrmions form a topological composite called skyrmionium, 90 which has been investigated in magnetic materials for next generation spintronics based devices 32,90–94 but not in ferroelectric materials. To the best of my knowledge, this is the first demonstration of skyrmion-to-skyrmionium “reaction” in ferroelectric materials. Moreover, I explain such “topological mechano-chemistry” based on the well-known Landau-Lifshitz-Kittel scaling laws 95–97 for domain size versus thickness for ferroelectric and ferromagnetic domains. Namely, decreasing sample thickness size increases the number of domains within the sample, which is consistent with the new domain walls creating the skyrmionium topology in the thinner structure. The observed skyrmion-to-skyrmionium phase transition offers a mechanical pathway to introducing more regions in the superlattice with local polarization curl and divergence, which are anticipated to enhance the macroscopic dielectric properties of the superlattice 36 . In both strain modes, a phase transition occurs when the PTO lattice cannot maintain proper stress/strain gradients for the skyrmion due to the compression, but how this gradient is dissipated, and the corresponding dynamical response of the polarization, is fundamentally different depending on the symmetry of the applied strain modes. Overall, our machine learning based multiscale NN/MM simulations illuminate the complex interplay between the symmetry of the strain and polarization boundary conditions in the formation of topological defects, while further demonstrating the ability of ferroelectric materials to materialize topological defects seen in the sister field of ferromagnetism 40,41 . 13 1.3.3 Electrical Control After utilizing my developed machine-learning forcefields incorporated within multi- scale simulation frameworks to explore optical and mechanical manipulation of various topological structures in PbTiO3 based ferroelectric structures 76,98 , I wanted to incorporate electric field dynamics to study elecrrical control of polar topologies. To incorporate electric field dynamics within this machine-learning framework without having to retrain the model with dynamics in presence of an electric field I have incorporated a simple born effective charge extension. The Born effective charge tensor for an atom in a crystal represents the linear response of the electronic polarization of crystal with respect displacements from its equilibrium position 99 : 𝑍 "# =𝑄+𝛺 𝑑𝑃 " 𝑑𝑢 # 3 $%& (1.1) Where Q is the bare coulomb charge, 𝛺 is the unit cell volume, P is the electronic polarization, and u is the displacement of the atom from its equilibrium position. The same coefficient describes the linear response of the force on an atom from an applied electric field, because both can be connected to second-order derivative of the energy with respect to atomic displacements and a macroscopic electric field 99 : 𝑍 "# = 𝑑𝐹 " 𝑑𝐸 # (1.2) Thus, as long as one operates within a reasonably linear regime it is reasonable to approximate the change in the forces predicted from a neural network model 𝐹 '' " as 14 𝐹 '' " (𝐸)~𝐹 '' " (0)+8 𝑑𝐹 " 𝑑𝐸 # 𝐸 # +𝑂(𝐸 ( ) (1.3) Assuming 𝐹 '' " has learned the ground truth density functional theory forces (DFT) well used to compute the born effective charges. This will incorporate the change in forces due not only the electric-field ion interactions, but also screening and redistribution electron density due to 𝜖 ) . Here I incorporate this extension within my existing 76,98 neural network quantum molecular dynamics 85,100 (NNQMD) forcefield for PbTiO3 to examine its dynamics under electric field. We first utilize the method by examining the dynamics of 180° domain walls (DWs) in PbTiO3, which has been well studied 101,102 . Switching dynamics reveal a layer-by-layer nucleation limited switching mechanism from the domain wall that is consistent with conventional understandings of ferroelectric switching. During switching we see both Bloch and Neele like switching of the polarization at the domain wall. While there has much debate about role of Bloch and Neele rotations in FE switching, we find this unsurprising as the domain boundary is known to have both Neele and Bloch components of polarization. After validating the dynamics within the well-studied 180° DWs, we examine electric field dynamics of a biaxially expanded and c-axially compressed PbTiO3 film to mimic strain conditions often seen PbTiO3 based oxides superlattices and strained thin films. We find application of an electric field along the a axis of uniformly polarized strained PbTiO3 can form a1/a2 flux-closure domains (i.e formation of a flux closure in the ab plane). Flux closures are quadrants with closed head-tail dipole moments that become of vortex if continuous rotation around a vortex core is made. Flux closure domains have technological interest due there potential for high bit density in non-volatile ferroelectric random access memory applications, and are predicted to be useful for sensor and transducer applications 103 . The formation was due to a plastic deformation of lattice in the form of a slip mechanism of the atomic planes resulting in a 15 Bloch rotation of the atoms. We then examined the switching dynamics of the flux closure domain where an intermediate vortex state is found along the switching pathway. 1.4 Modeling of neutron scattering experiments for development of theoretical frameworks for studying hydrogen storage candidates. Properties of molecular crystals and liquids such as water and ammonia are dictated by complex of dynamics of hydrogen bonding, dispersion interactions, and many body polarization effects whose understanding remain a fundamental scientific problem 104–110 . Ammonia (NH3) has garnered much attention due not only to its rich “water like” nature but also to a wide range of applications such as room-temperature synthesis of metal chalcogenides and development of exotic liquid metals 111,112 . In particular, over the past few years there has been growing movement for development of green ammonia-based fuel technologies 113–115 . Ammonia has a higher energy density, at 12.7 MJ/L, than even liquid hydrogen, at 8.5 MJ/L. Liquid hydrogen must be stored at cryogenic conditions of –253 °C, whereas ammonia can be stored at a much less energy-intensive –33 °C. Furthermore, thanks to a century of ammonia use in agriculture, a vast ammonia infrastructure already exists. Worldwide, some 180 million metric tons of ammonia is produced annually, and 120 ports are equipped with ammonia terminals. Development of technologies based on ammonia will be reliant on our ability to understand and model the complex physical and chemical interactions that give rise to its unique properties. However, much is still not understood about the dynamics of ammonia. For example, whether ammonia hydrogen bonds in its solid and liquid phases is still matter of open debate 105–110 . Vibrational spectroscopy often offers a window into signatures of molecular crystals and liquid’s physical and chemical interactions 109,116–118 . The collective vibrational behavior of an ensemble of molecules (e.g., in clusters , solids and liquids) is critical for predicting their 16 thermodynamic behavior , with important implications on their applications in energy, biological, and pharmaceutical systems. Many relevant properties (such as heat capacity, melting, conformation, polymorph formation) are directly or indirectly driven by molecular vibrations (or phonons in solid state). The vibration also reflects the local environment of the molecule and the intermolecular interactions (e.g., van der Waals forces and hydrogen bonding). Therefore, it is also crucial for understanding physical and chemical processes such as adsorption and catalysis. Despite the importance, quantitative understanding of molecular vibrations is not always straightforward. Experimentally, widely used spectroscopy methods such as Raman and Infrared scattering are only able to measure a subset of the vibrational modes, due to the lack of momentum from photons, the selection rules, and the challenges to probe the low frequency region 117 (sub-THz). Inelastic neutron scattering (INS) is ideally suited to measure the full phonon density of states. Compared to Raman or infrared spectroscopy, INS has no selection rules, is sensitive to hydrogen in molecular systems, and has its unique strength in discerning low to intermediate frequency modes – these are the modes that are particularly sensitive to intermolecular interactions. In addition, calculation of INS is rigorous and straightforward if the dynamics of nuclei can be solved, and explicit treatment of electronic structure is not required. All these features make INS an appealing technique to study phonons in molecular solids and liquids. However, to the best of our knowledge there is no available INS data for solid and liquid ammonia of high enough quality to be rigorously compared to theoretical models. The currently only available data is from the work of Jack Carpenter et al. from 2004 119 , which does not have good enough resolution, especially in the energy range of inter-molecular interactions, to be 17 rigorously compared to physical models for ammonia in its liquid and solid phases. High quality data is necessary for development of accurate models of the dynamic behavior of molecular solids and liquids, as there are multiple complicating factors that require careful considerations such as van der Waals interactions, nuclear quantum effects (NQEs), and phonon anharmonicity. The role of van der Waals force is particularly relevant in molecular solids as it is a significant part, if not dominant part, of the intermolecular interactions 109 . The conventional density functional theory 120 (DFT) cannot describe van der Waals interactions and empirical corrections are often included, leaving additional uncertainties when modeling such systems. Moreover, most molecular solids contain light elements such as H, for which NQEs could be significant, especially at low temperatures (even though the impact can also be observed at room temperature). Conventional lattice dynamics or molecular dynamics treat nuclei as classical point particles with no spread and thus NQEs are not considered. Last but not least, molecular solids are usually “soft” and tend to exhibit phonon anharmonicity. Such anharmonicity could be coupled with NQEs, making the analysis even more demanding. In chapter 5 I report measured vibrational density of states and dynamic structure factor for deuterated and protonated ammonia along the solid to liquid phase transition with inelastic neutron scattering using SEQUOIA 121 and VISION spectrometers at Oak Ridge National Laboratory. I find strongly anharmonic behavior of the inter-molecular phonon dynamics in solid phase, but little change in the vibrational spectrum for the intra-molecular modes is observed as a function of temperature in solid phase. In liquid phase coherent excitations in the quasi-elastic scattering regime resembling acoustic phonons are observed in the dynamic structure factor, as well as hardening of the high energy N-H stretching modes that is indicative of a decrease in strength of inter-molecular interactions in the liquid phase. 18 These results are compared to DFT simulations. We find that the computed inelastic spectrums for ammonia depend strongly on the choice of dispersion correction for the exchange- correlation functional used, and all poorly agree with measured spectrum. Through frozen phonon and ab-initio path integral molecular dynamics (PIMD) simulations 122–125 we can qualitatively show these discrepancies arise from phonon anharmonicity coupled with nuclear quantum effects; however, quantitative agreement is lacking due to loss in resolution of the computed spectrum. This is due to long simulation run-time and system size are required to compute low energy vibrational spectrum from Fourier transform of correlation functions in PIMD simulations rather static matrix diagonalization used in the typical lattice dynamics approach. As the majority of the computational expense for ab-intio PIMD simulations comes from having to compute multiple replica DFT simulations, the computational cost can be greatly decreased if the underlying DFT simulations can be replaced by much cheaper computational models. In this regard, neural-network quantum molecular dynamics (NNQMD) simulations 69 based on machine learning offer a promising tool as they revolutionizing atomistic modeling of materials by following the trajectories of all atoms with quantum-mechanical accuracy at a drastically reduced computational cost. NNQMD can not only predict accurate interatomic forces but can capture quantum properties such as electronic polarization 100 and electronic excitation 76 , thus the ‘Q’ in NNQMD. A more recent breakthrough in NNQMD is drastically improved accuracy of force prediction over those previous models, which was achieved through rotationally equivariant neural networks based on a group theoretical formulation of tensor fields 126–128 . Here we show with this increased accuracy we can perform large scale NNQMD based PIMD simulations that are in good quantitative agreement with measured INS spectrums 19 for ammonia, cementing the role of NQE’s on the vibrational properties of ammonia. This type of approach is scalable to any type of material system, offering a robust method to quantify the role of nuclear quantum effects on material vibrational dynamics. 20 Chapter 2 : Methods 2.1 Molecular Dynamics Method. The Molecular dynamics methods follows the trajectory of all atoms at any given time from their initial positions {𝑟 * } and velocities {𝑣 * } using Newton’s equation of motion. 𝐹 * =−∇ + 𝑉({𝑟 * }) (2.1) 𝑎 * (𝑡) = 𝐹 * (𝑡) 𝑚 * (2.2) The trajectories of each atom are predicted by numerically solving the Newton’s equations of motion, where interatomic interactions between particles are defined by a parameterized inter-atomic force field that gives the total potential energy surface. Common example of such a forcefield is the Lennoard-Jones potential 𝑉(𝑟) = 8 𝑉G𝑟 *, H= *-, 8 4𝜖JK 𝜎 𝑟 *, M .( −K 𝜎 𝑟 *, M / N *-, (2.3) That depends on the atomic pair-wise distance 𝑟 *, , and free parameters 𝜎 and 𝜖 that determine the interaction length and strength. Numerical Integration of Newton’s equations of motion are typically done through the Verlet scheme: For each atom i 1. Compute 𝑣 * O𝑡+ 0 ( P=𝑣 * (𝑡)+ . ( 𝑎 * (𝑡)Δ 21 2. Compute 𝑟 * (𝑡+Δ) =𝑟 * (𝑡)+ 𝑣 * O𝑡+ 0 ( PΔ 3. Compute 𝑎 * (𝑡+Δ) from {𝑟 . (𝑡+Δ) ,…,𝑟 * (𝑡+Δ),… ,𝑟 ' (𝑡+Δ)} 4. Compute 𝑣 * (𝑡+Δ) =𝑣 * O𝑡+ 0 ( P+ . ( 𝑎 * (𝑡+Δ)Δ which is repeated for the desired number of time steps. 2.2 Ab-initio Molecular Dynamics Ab-initio molecular dynamics produces molecular forces from direct calculation of the underlying electronic structure rather than empirically or otherwise fitted classical potentials. In this framework the nucleons are treated as classical point particles whose forces are calculated quantum mechanically, while electrons are treated as fully quantomechanical objects and typically handled through the framework of density functional theory (DFT), but any quantum mechanical framework is valid. 2.2.1 Derivation from Quantum Mechanics: We will start from the non-relativistic Schrödinger equation for a general atomic system: 𝑖ℏ 𝜕 𝜕𝑡 Φ({𝒓 𝒊 },{𝑹 𝑰 }.𝑡)=𝐻Φ({𝒓 𝒊 },{𝑹 𝑰 }.𝑡) (2.4) where {𝒓 𝒊 } is the set of the electronic coordinates and {𝑹 𝑰 } the nuclear coordinates. The general multi-body Hamiltonian is: 𝐻 =−8 ℏ ( 2𝑀 3 ∇ 4 ( 3 −8 ℏ ( 2𝑚 5 ∇ + ( + * 8 𝑒 ( ]𝒓 𝒊 −𝒓 𝒋 ] − *-, 8 𝑒 ( 𝑍 3 |𝑹 𝑰 −𝒓 𝒊 | + 3,* 8 𝑒 ( 𝑍 3 𝑍 8 ]𝑹 𝑰 −𝑹 𝑱 ] 3-8 (2.5) 22 which can be more compactly written: 𝐻 =8 ℏ ( 2𝑀 3 ∇ 4 ( 3 +8 ℏ ( 2𝑚 5 ∇ + ( + * 𝑉 :;5 ({𝒓 𝒊 },{𝑹 𝑰 }.𝑡) =8 ℏ ( 2𝑀 3 ∇ 4 ( 3 +𝐻 5 ({𝒓 𝒊 },{𝑹 𝑰 }.𝑡) (2.6) For solving the Schrödinger equation we will make the ansatz that the wavefunction is separable in terms of the electronic and nuclear coordinates. This will lead to mean-field approximation used in many frameworks for solving the time-dependent Schrödinger Equation (TDSE). We will assume: Φ=Ψ({𝒓 𝒊 },𝑡)Υ({𝑹 𝑰 },𝑡) (2.7) Plugging our ansatz into the TDSE we arrive at: 𝑖ℏ⌊Υ({𝑹 𝑰 },𝑡)⟩ 𝜕 𝜕𝑡 ⌊𝛹({𝒓 𝒊 },𝑡)⟩+𝑖ℏ⌊Ψ({𝑹 𝑰 },𝑡)⟩ 𝜕 𝜕𝑡 ⌊Υ({𝒓 𝒊 },𝑡)⟩ = 8 ℏ ( 2𝑀 3 ∇ 4 ( 3 +8 ℏ ( 2𝑚 5 ∇ + ( + * 𝑉 :;5 ({𝒓 𝒊 },{𝑹 𝑰 }.𝑡) )|Ψ({𝒓 𝒊 },𝑡)Υ({𝑹 𝑰 },𝑡)⟩ (2.8) Multiplying through by ⟨𝛹({𝒓 𝒊 },𝑡) and asserting orthonormality we obtain: 𝑖ℏ 𝜕 𝜕𝑡 ⌊Υ({𝒓 𝒊 },𝑡)⟩ =⟨𝛹({𝒓 𝒊 },𝑡)|(−8 ℏ ( 2𝑀 3 ∇ 4 ( 3 −8 ℏ ( 2𝑚 5 ∇ + ( + * 𝑉 :;5 ({𝒓 𝒊 },{𝑹 𝑰 }.𝑡) ))|Ψ({𝒓 𝒊 },𝑡)Υ({𝑹 𝑰 },𝑡)⟩ =−8 ℏ ( 2𝑀 3 ∇ 4 ( 3 |Υ({𝑹 𝑰 },𝑡)⟩ +⟨𝛹({𝒓 𝒊 },𝑡)|8 ℏ ( 2𝑚 5 ∇ + ( + * 𝑉 :;5 ({𝒓 𝒊 },{𝑹 𝑰 }.𝑡) ))|Ψ({𝒓 𝒊 },𝑡)Υ({𝑹 𝑰 },𝑡)⟩ → 23 𝑖ℏ 𝜕 𝜕𝑡 ⌊Υ({𝒓 𝒊 },𝑡)⟩ =8 ℏ ( 2𝑀 3 ∇ 4 ( 3 |Υ({𝑹 𝑰 },𝑡)⟩+ ⟨Ψ({𝒓 𝒊 },𝑡)|𝐻 5 ({𝒓 𝒊 },{𝑹 𝑰 }.𝑡)|Ψ({𝒓 𝒊 },𝑡)⟩|Υ({𝑹 𝑰 },𝑡)⟩ (2.9) similarly, by multiplying through by ⟨Υ({𝑹 𝑰 },𝑡)| we obtain: 𝑖ℏ 𝜕 𝜕𝑡 ⌊Ψ({𝒓 𝒊 },𝑡)⟩ ==−8 ℏ ( 2𝑚 5 ∇ + ( 3 |Ψ({𝒓 𝒊 },𝑡)⟩+ ⟨Υ({𝑹 𝑰 },𝑡)|𝑉 :;5 ({𝒓 𝒊 },{𝑹 𝑰 }.𝑡)|Υ({𝑹 𝑰 },𝑡)⟩|Ψ({𝒓 𝒊 },𝑡)⟩ (2.10) These coupled equations form a time-dependent, self-consistent set of equations where the wavefunctions for the electrons and nucleons propagate in a mean field determined by an expectation value of the opposing wavefunction. To solve the equations and attempt to recover a form of classical-equations for the nucleons we need to be able to treat the nucleons as classical point particles. To do this we will first write the nucleon wave function in polar form: |Υ({𝑹 𝑰 },𝑡)⟩=𝐴({𝑹 𝑰 },𝑡)𝑒 *<({𝑹 𝑰 },A) ℏ (2.11) where A and S are real, A>0, and∫{𝑑𝑹 𝑰 } 𝐴 ( ({𝑹 𝑰 },𝑡)=1. Subbing into equation 2.9 we obtain 𝑖ℏ 𝜕|Υ({𝑹 𝑰 },𝑡)⟩ 𝜕𝑡 = m − 𝜕𝑆 𝜕𝑡 𝐴+𝑖ℏ 𝜕𝐴 𝜕𝑡 o𝑒 *<({𝑹 𝑰 },A) ℏ = 24 8 −ℏ ( 2𝑀 3 3 p𝐴∇ 4 ( 𝑒 *<({𝑹 𝑰 },A) ℏ +2∇ 3 𝐴∙∇ 3 𝑒 *<({𝑹 𝑰 },A) ℏ +𝑒 *<({𝑹 𝑰 },A) ℏ ∇ 3 ( 𝐴r +⟨Ψ({𝒓 𝒊 },𝑡)|𝐻 5 ({𝒓 𝒊 },{𝑹 𝑰 }.𝑡)|Ψ({𝒓 𝒊 },𝑡)⟩𝐴𝑒 *<({𝑹 𝑰 },A) ℏ = 𝑒 *<({𝑹 𝑰 },A) ℏ 8 (−𝑖ℏ(2∇ 3 𝐴∙∇ 3 𝑆+𝐴∇ 3 ( 𝑆)+𝐴(∇ 3 𝑆) ( +−ℏ ( ∇ 3 ( 𝐴) 2𝑀 3 3 + (⟨Ψ({𝒓 𝒊 },𝑡)|𝐻 5 ({𝒓 𝒊 },{𝑹 𝑰 }.𝑡)|Ψ({𝒓 𝒊 },𝑡)⟩𝐴)𝑒 *<({𝑹 𝑰 },A) ℏ (2.12) Equating real and imaginary parts we arrive at the following two equations: 𝜕𝑆 𝜕𝑡 + 8 (∇ 3 𝑆) ( 2𝑀 3 + 3 ⟨Ψ({𝒓 𝒊 },𝑡)|𝐻 5 ({𝒓 𝒊 },{𝑹 𝑰 }.𝑡)|Ψ({𝒓 𝒊 },𝑡)⟩ =ℏ ( 8 ∇ 3 ( 𝐴 2𝐴𝑀 3 3 (2.13) 𝜕𝐴 𝜕𝑡 +8 (2∇ 3 𝐴∙∇ 3 𝑆+𝐴∇ 3 ( 𝑆) 2𝑀 3 3 =0 (2.14) Equation 2.13 represents a continuity like equation for the wave amplitude of the nucleons, while equation 2.14 is much more interesting and actually captures the classical dynamics of our system we are interested in. If we make the semi-classical approximation and take the limit as ℏ →0 equation 2.13 becomes 𝜕𝑆 𝜕𝑡 +8 (∇ 3 𝑆) ( 2𝑀 3 + 3 ⟨Ψ({𝒓 𝒊 },𝑡)|𝐻 5 ({𝒓 𝒊 },{𝑹 𝑰 }.𝑡)|Ψ({𝒓 𝒊 },𝑡)⟩=0 (2.15) 𝜕𝑆 𝜕𝑡 +𝐻 EFGHH s{𝑹 𝑰 },{∇ 3 𝑆}H =0 (2.16) which we can now recognize as the classical Hamilton-Jacobi equations with 𝑆 being the action, the canonical momentum 𝑷 𝑰 =∇ 3 𝑆 =𝑀 3 𝑹 𝑰 ̇ , and classical potential 𝑉 = ⟨Ψ({𝒓 𝒊 },𝑡)|𝐻 5 ({𝒓 𝒊 },{𝑹 𝑰 }.𝑡)|Ψ({𝒓 𝒊 },𝑡)⟩. Using Hamilton’s equations of motion we can write: 𝑷 𝑰 ̇ =−∇ 4 𝐻 (2.17) 25 Implying 𝑀𝑹 𝑰 ̈ =−∇ 4 ⟨Ψ({𝒓 𝒊 },𝑡)|𝐻 5 ({𝒓 𝒊 },{𝑹 𝑰 }.𝑡)|Ψ({𝒓 𝒊 },𝑡)⟩ (2.18) As we can see the nucleons move in an effective classical field generated by the electrons. To find equations of motion that govern the electrons we will take the limit that nucleons become point particles, ie there probability amplitudes become delta functions: 𝑖ℏ 𝜕 𝜕𝑡 ⌊Ψ({𝒓 𝒊 },𝑡)⟩ =−8 ℏ ( 2𝑚 5 ∇ + ( 3 |Ψ({𝒓 𝒊 },𝑡)⟩ +⟨Υ({𝑹 𝑰 },𝑡)|𝑉 :;5 ({𝒓 𝒊 },{𝑹 𝑰 }.𝑡)|Υ({𝑹 𝑰 },𝑡)⟩|Ψ({𝒓 𝒊 },𝑡)⟩ (2.19) 𝑖ℏ 𝜕 𝜕𝑡 ⌊Ψ({𝒓 𝒊 },𝑡)⟩ =−8 ℏ ( 2𝑚 5 ∇ + ( 3 |Ψ({𝒓 𝒊 },𝑡)⟩ +wx 𝑑𝑹 𝑰 𝛿G𝑹 𝑰 −𝑹 𝑰 (𝑡)H 3 𝑉 :;5 ({𝒓 𝒊 },{𝑹 𝑰 }.𝑡)|Ψ({𝒓 𝒊 },𝑡)⟩ (2.20) 𝑖ℏ 𝜕 𝜕𝑡 ⌊Ψ({𝒓 𝒊 },𝑡)⟩ =−8 ℏ ( 2𝑚 5 ∇ + ( 3 |Ψ({𝒓 𝒊 },𝑡)⟩+𝑉 :;5 ({𝒓 𝒊 },{𝑹 𝑰 (𝑡)},𝑡)|Ψ({𝒓 𝒊 },𝑡)⟩ (2.21) 𝑖ℏ 𝜕|Ψ({𝒓 𝒊 },𝑡)⟩ 𝜕𝑡 =𝐻 5 (s{𝒓 𝒊 },{𝑹 𝑰 (𝑡)},𝑡H|Ψ({𝒓 𝒊 },𝑡)⟩ (2.22) 26 where {𝑹 𝑰 (𝑡)} represents the set of instantaneous positions of the nucleons governed by 2.18. These equations are the frame-work from which the forms the three main frameworks of first principle molecular dynamics: Born-Oppenheimer molecular dynamics, Ehrenfest molecular dynamics, and Surface Hopping Molecular Dynamics. 2.2.2 Born-Oppenheimer(BO) Molecular Dynamics In BO MD the electrons are confined to a fixed state configuration usually the ground state(however thermalized states can also be used). In this case the electrons are assumed to instantaneously respond to the nuclei motion and at each time step the ground state problem is solved. 𝐻 5 (s{𝒓 𝒊 },{𝑹 𝑰 (𝑡)},𝑡H|Ψ({𝒓 𝒊 },𝑡)⟩ =𝐸|Ψ({𝒓 𝒊 },𝑡)⟩ (2.23) 𝑀𝑹 𝑰 ̈ =−∇ 4 ⟨Ψ({𝒓 𝒊 },𝑡)|𝐻 5 ({𝒓 𝒊 },{𝑹 𝑰 }.𝑡)|Ψ({𝒓 𝒊 },𝑡)⟩ (2.24) After computing the force in equation 2.24 the atoms positions are updated and equation 2.23 is solved again to get the new ground state configuration with the updated atomic positions. In this case the timestep can be on the order of nuclear degrees of freedom (femto-second). BO MD is often referred to adiabatic MD. 2.2.3 Ehrenfest Molecular Dynamics Ehrenfest MD solves the TDSE and MD equations simultaneously on the fly during the MD simulation, and is a type of non-adiabatic MD. Typically one expands the electronic wave- function in an adiabatic basis. |Ψ(t)⟩ =8 𝑐 , (𝑡)]𝜙 , } , (2.25) Which one can plug into equations 2.18 and 2.22 to obtain 27 𝑀 3 𝐑 3 ̈ (𝑡) = 8 −]𝑐 , (𝑡)] ( , 𝛻 3 𝐸 , −8 𝑐 , ∗ (𝑡)𝑐 J (𝑡)G𝐸 , −𝐸 J H𝑑 ,J 3 , ,J (2.26) 𝑖𝑐 K ̇(𝑡) =𝑐 , (𝑡)𝐸 , −𝑖8 𝑐 J (𝑡)𝑹 𝑰 ̇ 𝑑 ,J 3 3,J (2.27) where 𝑑 ,J 3 = 𝜙 , ]𝛻 3 ]𝜙 J } (2.28) is the non-adiabatic coupling vector and 𝐸 , =𝜙 , ]𝐻]𝜙 , } (2.29) is the Adiabatic Energy. These equations are solved on the fly during the MD simulation which allows for changes in the occupation of the electronic states due to electron-phonon interaction. 2.2.3 Surface-Hopping Molecular Dynamics In the surface hopping approach instead of solving the TDSE and MD equations together on the fly, the probability of an electronic transition from one state to another is computed at each MD step from the non-adiabatic coupling vector. This can be seen from squaring and taylor expanding equation 2.27 𝑖𝑐 K ̇(𝑡) =𝑐 , (𝑡)𝐸 , −𝑖8 𝑐 J (𝑡)𝑹 𝑰 ̇ 𝑑 ,J 3 →|𝑐 K ̇(𝑡)| ( =−2𝑅𝑒J8 𝑐 , ∗ 𝑐 J 𝑹 𝑰 ̇ 𝑑 ,J 3 J,3 N 3,J (2.30) |𝑐 , (𝑡+𝑑𝑡)| ( − |𝑐 , (𝑡)| ( |𝑐 , (𝑡)| ( ≈− 2𝑑𝑡 |𝑐 , (𝑡)| ( Re8 𝑐 , ∗ 𝑐 J 𝑹 𝑰 ̇ 𝑑 ,J 3 J,3 (2.31) Where dt is much smaller than the MD time-step. Equation * represents the total population flux of state j. Thus we can write 𝑃 , →J (𝑡) ≈− 2𝑑𝑡 |𝑐 , (𝑡)| ( Re8 𝑐 , ∗ 𝑐 J 𝑹 𝑰 ̇ 𝑑 ,J 3 3 (2.32) 28 Where 𝑃 , →J (𝑡) is the flux of state j to state k. This can be used to compute probability of an inter-state transition between adiabatic states at each MD time-step. 2.3 Density-Functional Theory(DFT) DFT takes the unique approach of solving the electronic structure problem in terms of the charge density rather than the many body wavefunction, Similar to Hartree Fock approach , allows one to solve N one electron problems, which leads to an O(N 3 ) scaling for solving the problem, versus exponential scaling for directly solving the Schrödinger equation. All many- body effects are swept into EXC : Exchange-Correlation Functional. DFT is exact in theory, but the true EXC is unknown and requires approximation in practice. DFT is by far the most popular approach for solving the electronic structure problem in solid state systems. 2.3.1 DFT Derivation For any atomic system in the we can write the Hamiltonian 𝐻 =𝑇 +𝑈 +𝑉 (2.33) with: 𝑇 = 1 2 x 𝑑 M 𝑟 ∇𝜓 (𝒓) N ∇𝜓 (𝒓) (2.34) 𝑈 = 1 2 x x 𝑑 M 𝑟𝑑 M 𝑟 O 𝜓 (𝒓) N 𝜓 (𝒓 O ) N 𝜓 (𝒓)𝜓 (𝒓 O ) |𝒓−𝒓 O | (2.35) 𝑉 =x 𝑑 M 𝑟 𝑣(𝒓)𝜓 (𝒓) N 𝜓 (𝒓) (2.36) where we are using atomic units and 𝑣(𝒓) is the external potential. 𝜓(𝒓) N 𝜓(𝒓) are creation and annihilation operators. We know the external potential for a system fixes the Hamiltonian and hence fixes the ground state electron density 𝑛(𝒓) =Ψ]𝜓(𝒓) N 𝜓(𝒓)]Ψ}. We will first show that 29 the converse is true, a given electron density fixes the external potential 𝑣(𝒓),(up to an arbitrary constant), meaning the external potential is a unique functional of the electron density. Proof: Suppose there exists a 𝑣′(𝒓) such that it gives rise to a Ψ O (𝒓) that gives rise to the same density 𝑛(𝒓). We can write: 𝐸 O = ⟨Ψ O (𝒓)|𝐻 O |Ψ O (𝒓)⟩< ⟨Ψ(𝒓)|𝐻 O |Ψ(𝒓)⟩ =⟨Ψ(𝒓)|𝐻 O +𝑉−𝑉 O |Ψ(𝒓)⟩ (2.37) which implies: 𝐸 O <𝐸+x 𝑑 M 𝑛(𝒓)(𝑣(𝒓)−𝑣′(𝒓)) (2.38) by using the same logic for the primed variables, we can write: 𝐸 <𝐸 O −x 𝑑 M 𝑛(𝒓)(𝑣(𝒓)−𝑣′(𝒓)) (2.39) adding equation 2.38 and 2.39 yields the contradiction: 𝐸+𝐸 O <𝐸+𝐸 O (2.40) implying that 𝑣(𝒓) is a unique functional of 𝑛(𝒓). Thus the systems total energy is a unique functional of 𝑛(𝒓). Since the correct ground state wavefunction minimizes the energy functional E it follows trivially the functional is minimized by the correct electron density. Next, we will break up the energy functional into: 𝐸[𝑛] =𝑉[𝑛]+𝐹[𝑛] (2.41) where: 𝐹[𝑛] =𝑇 < [𝑛]+𝐻[𝑛]+𝐸 PE [𝑛] (2.42) 30 where 𝑇 < [𝑛] is the Kinect energy for non-interacting electrons, 𝐻[𝑛] is the Hartree mean field electron-electron coulomb interaction: 𝐻[𝑛]=𝑑 M 𝑟𝑑 M 𝑟 O 𝑛(𝒓 O )𝑛(𝒓) |𝒓−𝒓 O | (2.43) and 𝐸 PE [𝑛] is the Exchange-Correlation Functional. The Exchange-Correlation Functional is very difficult to calculate, and in practice is approximated by a multitude of different methods. See Exchange-Correlation Functional section. Assuming we know 𝐸 PE [𝑛] we minimize the energy functional using calculus of variations subject to the constraint ∫𝑑 M 𝑟 𝑛(𝒓)=𝑁 𝛿 𝛿𝑛(𝒓) [𝐸[𝑛(𝒓)]−𝜇x 𝑑 M 𝑟 𝑛(𝒓)]=0 (2.44) 𝛿𝑇 < [𝑛] 𝛿𝑛(𝒓) +x 𝑑 M 𝑟′ 𝑛(𝒓 O ) |𝒓−𝒓 O | +𝑣(𝒓)+ 𝛿𝐸 PE 𝛿𝑛(𝒓) = 𝜇 (2.45) we can rewrite this as: 𝛿𝑇 < [𝑛] 𝛿𝑛(𝒓) +𝑉 5QQ (𝒓)=𝜇 (2.46) We can now see that equation 2.46 would be the same Euler-Lagrange equation for a system of noninteracting electrons in an external effective potential. Meaning we can break the many electron system in N-equations of non-interacting electrons. We can write: 𝑇 < [𝑛] =− 1 2 8 𝜓 * ∗ ∇ ( 𝜓 * ' *%. (2.47) 𝑛(𝒓)= 8 𝜓 * ∗ 𝜓 * ' *%. 𝜃(𝜇−𝜖 * ) (2.48) Where 𝜇 is now recognized as the chemical potential and is fixed by the constraint ∫𝑑 M 𝑟 𝑛(𝒓)=𝑁. This yields the set of equations: 31 K− 1 2 ∇ ( +𝑉 5QQ (𝒓)M𝜓 * =𝜖 * 𝜓 * (2.49) 𝑉 5QQ (𝒓)=𝑣 5PA (𝒓)+𝑣 RGSAS55 G𝒓;𝑛(𝒓)H+𝑣 PE G𝒓;𝑛(𝒓)H (2.50) 𝑣 RGSAS55 G𝒓;𝑛(𝒓)H = x 𝑑𝐫′ 𝑛(𝐫′) |𝐫−𝐫 O | (2.51) 𝑣 PE G𝒓;𝑛(𝒓)H=𝜖 PE (𝒓)𝑛(𝒓)+ 𝑛(𝒓) 𝛿𝜖 PE (𝒓) 𝛿𝑛(𝒓) (2.52) These are known as the Kohn-Sham(KS) equations and have to be solved self-consistently : 1. Initialize Charge Density (Usually Sum of Atomic orbitals) 2. Solve KS Equations 3. Get New Charge Density 4. Repeat 2-3 Until Total Energy and Charge Density Stop Changing 5. Compute properties (Density of States, Atomic Forces,…) 2.3.2 Exchange Correlation Functional The development for exchange correlation functionals are based on the uniform electron gas whose correlation potential can be calculated by Monte Carlo Methods. Typical Flavors of the exchange correlation functional are • LDA : 𝐹G𝑛(𝑟)H. Local Functional of Density. Electrons locally behave like homogeneous electron gas. • GGA : 𝐹(𝑛(𝑟), 𝛻𝑛(𝑟)).Gradient Corrections applied to LDA. • Hybrid Functionals : Mix of Fock Exact Exchange energy and GGA/LDA exchange-correlation energy 32 2.2.3 DFT and Periodic Boundary Conditions (PBCs) For crystalline and other solid state systems that have some periodicity, we often use periodic boundary conditions to describe their bulk properties. In this case wavefunctions are subject to translational invariance 𝜓 :J (𝑟+𝑅)=𝑒𝑥𝑝(𝑖𝑘𝑟)𝑢 :J (𝑟) (2.53) 𝑢 :J (𝑟+𝑅) = 𝑢 :J (𝑟) (2.54) where k is Bloch vector that lies in the first Brillouin zone(BZ) of a periodic-crystal and R is a lattice vector. In this case desired quantities are expressed as integration over the first BZ 𝜌(𝑟) = 8 x 𝑑 M 𝑘𝑓 :J 𝜓 :J (𝑟) : (2.55) 𝑓 :J →𝑜𝑟𝑏𝑖𝑡𝑎𝑙 𝑜𝑐𝑐𝑢𝑝𝑎𝑡𝑖𝑜𝑛 In practice this integration is weighted sum over a “k-point” grid 𝜌(𝑟) =8 8 𝑤 :J 𝑓 :J 𝜓 :J (𝑟) J : (2.56) 𝑓 :J →𝑜𝑟𝑏𝑖𝑡𝑎𝑙 𝑜𝑐𝑐𝑢𝑝𝑎𝑡𝑖𝑜𝑛 𝑤 :J →𝑘𝑝𝑜𝑖𝑛𝑡 𝑤𝑒𝑖𝑔ℎ𝑡 𝑖𝑛 𝑖𝑛𝑡𝑒𝑔𝑟𝑎𝑡𝑖𝑜𝑛 𝑠𝑐ℎ𝑒𝑚𝑒 2.2.4 Plane Wave Basis and Pseudopotentials To solve the KS equations we must expand our individual electronic wave-functions in a basis set. For solid state systems the choice of plane-waves naturally arises as many band structures of elements can be interpreted in the nearly free electron picture, within which plane wave solutions naturally arise. In a plane wave expansion Bloch-wavefunctions, the charge density and the potential are expanded as 𝑢 :J (𝒓) = 1 (Ω E5FF ) . ( 8 𝐶 𝑮:J exp(𝑖𝑮∙𝒓) 𝑮 (2.58) 33 𝜌(𝒓) = 8 𝜌 𝑮 exp(𝑖𝑮∙𝒓) 𝑮 (2.59) 𝑉(𝒓)=8 𝑉 𝑮 exp(𝑖𝑮∙𝒓) U (2.60) One can quickly see this plane wave basis allows for fast evaluation of potentials in real and G space through FFTs. The plane waves are expanded up to wave-vectors within a cutoff energy, 𝟏 𝟐 (𝒌+𝑮) 𝟐 <𝐸 EXA . Plane waves offer the advantage of being “un-biased” : they are delocalized and do not depend on nuclear positions. Improvement of solution can be obtained by simply increasing the cutoff energy. However, due to the well-known “no free lunch” theorem, there must be a tradeoff when using plane-waves, and here is the need to use pseudopotentials. The singularity of the bare 1/𝑟 coulomb potential results in needing a large number of plane-waves in order to obtain convergence as the wave-function oscillates rapidly near the core of the potential, which is illustrated in figure 1b. Figure 1. Pseudopotentials. (a) Shows local component of PAW pseudopotential for W compared to all-electron wave-function. (b) Shows comparisons of pseudo-wavefunction and all-electron wave-function. 34 Pseudopotentials freeze the core-electrons which do not participate in bonding and mainly provide screening of the core for the chemically active valence electrons. Beyond the cutoff radius of pseudopotential rc , the pseudopotential and the all-electron(AE) potential agree, which is illustrated in figure 1a. Within the core the AE potential diverges while the pseudopotential is shielded. The most common types of pseudo-potentials used today are Norm Conserving Pseudo-Potentials(NCPPs) 129 , Ultrasoft Pseudo-Potentials (USPPs) 130 , and the Projector Augmented-Wave Method (PAW) 131 . 2.2.5 Real-Time TDDFT with Maxwells equations. To study the thermalization of electronic occupation numbers upon excitation by a laser pulse, we developed RT-TDDFT simulations coupled with Maxwells equations 67 into our quantum molecular dynamics (QMD) software, QXMD 132,133 . RT-TDDFT simulation follows the time evolution of Kohn-Sham (KS) wave functions {𝜓 HY (𝐫,𝑡)} at spatial position r and time t: 𝑖ℏ Z ZA 𝜓 HY (𝐫,𝑡) = § . ([ O ℏ * Z Z𝐫 + 5 E 𝐀(𝑡)P ( +𝑣 +]^ (𝐫)+𝑣 _`a [𝐫,𝑡;𝜌(𝐫,𝑡)]©𝜓 HY (𝐫,𝑡), (2.61) where ℏ is the Planck constant, c is the light speed, m and e are electron mass and charge, s and 𝜎 are KS-orbital and spin indices, A is the vector potential, 𝑣 +]^ is the ionic pseudopotential, and 𝑣 _`a is the Hartree and exchange-correlation potential which is a functional of the electron density, 𝜌(𝐫,𝑡)= ∑ 𝑓 HY |𝜓 HY (𝐫,𝑡)| ( HY . (2.62) In Eq. 2.62, 𝑓 HY are the initial occupation numbers of the KS orbitals, which we assume to be the ground-state occupation. The vector potential in Eq. (S1) is a sum of external and induced terms, 𝐀(𝑡)=𝐀 b`c (𝑡)+𝐀 +^d (𝑡). (2.63) Since experimental laser-spot sizes are typically on the order of micrometer or larger, we here consider an external laser pulse to be spatially uniform, with the form of 35 𝐀 b`c (𝑡) =− E e "#$ 𝐄 b`c sin(𝜔 b`c 𝑡)sin ( O fA g "#$ P (0<𝑡 <𝜏 b`c ) , (2.64) where 𝜔 b`c and 𝜏 b`c are the angular frequency and duration of the pulse, and 𝐄 b`c is the peak electric field. The induced vector potential, 𝐀 +^d (𝑡), satisfies the equation . E % Z % ZA % 𝐀 +^d (𝑡)= hf E 𝐣(𝑡), (2.65) where the spatially-averaged electric current is computed by an integration 𝐣(𝑡)=− 5 [i ∫ 𝑑𝐫²∑ 𝑓 HY Re³𝜓 HY ∗ (𝐫,𝑡) ℏ * ∇𝜓 HY (𝐫,𝑡)´ HY + 5 E 𝐀(𝑡)µ i (2.66) over the entire simulated volume Ω. In RT-TDDFT simulation, time-dependent KS equations, Eq. (S1), for electrons are solved concurrently with Eq. (S5) for electromagnetic field. 2.3 Neural Network Quantum Molecular Dynamics Neural-network quantum molecular dynamics (NNQMD) simulations 69 based on machine learning are revolutionizing atomistic modeling of materials by following the trajectories of all atoms with quantum-mechanical accuracy at a drastically reduced computational cost. NNQMD can not only predict accurate interatomic forces but can capture quantum properties such as electronic polarization 100 and electronic excitation 76 , thus the ‘Q’ in NNQMD. Here I will discuss the NNQMD model development for excited state and ground state PbTiO3 (PTO) dynamics, hybrid NN/MM forcefield development for skyrmions in STO/PTO superlattices, and forcefield development for ammonia for computing vibrational properties for neutron scattering experiments which are used in chapters 4 and 5. 2.3.1 NNQMD for ground state PTO 36 The neural-network quantum molecular dynamics (NNQMD) framework is outlined in Fig. 2a. The NNQMD framework employed utilizes a feed-forward neural network with two hidden layers and twenty nodes in each layer. I use a modified hyperbolic tangent as the activation function. To capture local atomic environments, we employ an energy model based on the symmetric functions proposed by Behler and Parrinello 134 : 𝐺 * SGj =∑ 𝑒 ;klm &' ;m ( n % ∙𝑓 E G𝑅 *, H , (2.67) 𝐺 * G:o =2 .;p ∑ G1+𝜆cos𝜗 *,J H p 𝑒 ;kqm &' % rm ') % rm &) % s ∙𝑓 E G𝑅 *, H∙𝑓 E G𝑅 ,J H∙𝑓 E (𝑅 *J ) , ,Jt* (2.68) 𝑓 E G𝑅 *, H= . ( (𝑐𝑜𝑠O fm &' m * P−1) if 𝑅 *, ≤𝑅 E ,else 0 (2.69) Equations 2.67 and 2.68 are the radial and angular symmetric functions for i-th atom as a function of the interatomic distance Rij between i-th and j-th atoms. Hyperparameters 𝜂, 𝑅 < , 𝜍, and 𝜆 are used to capture the local atomic environment. The cutoff function 𝑓 E (𝑅) in Eq. 2.69 smoothly truncates the symmetric functions at the cutoff distance Rc = 7.5 Å. We found that G rad with 𝜂 = {0.5, 1.0, 3.0} and Rs = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0} Å provides sufficient accuracy and computational efficiency for Pb, Ti and O atoms. Training was performed on a data set consisting of 1551 frames with various densities, tetragonalities and temperatures as describes in Table 1. Training data was generated in the canonical (NVT) ensemble using density functional theory (DFT) based quantum molecular dynamics (QMD) simulation on 4×4×4 unit-cells of PbTiO3 (PTO). 37 Figure 2.NNQMD model development for PbTiO3. (a) NNQMD framework for energy, force and stress prediction. (b)-(d) Root mean square error (RSME) for energy, force, and stress tensor as a function of training epochs. (e)-(j) Comparison of partial radial distribution functions between NNQMD and DFT using SCAN exchange-correlation functional. 38 Table 1. Summary of generated training data. First 6 columns describe box dimensions and temperature of QMD simulations. The 7 th column shows the length of trajectory generated and 8 th column shows the number of frames used from each trajectory and rate of frame skipped, i.e., only using every other or every 5 th frame from the MD trajectory. Frames were skipped in order to avoid oversampling in regions of similar atomic coordinates. 𝑎&𝑏 (Å) 𝑐 (Å) 𝛼 (°) 𝛽 (°) 𝛾 (°) Temp. (K) Total MD frames (dt =1.2fs) Used MD frames (skip rate) 3.90 4.15 90 90 90 300K 701 234(3) 3.90 4.15 90 90 90 1200K 468 234(2) 3.80 4.25 90 90 90 300K 750 150(5) 3.90 4.41 90 90 90 300K 589 120(5) 3.87 4.34 90 90 90 300K 659 132(5) 3.84 4.36 90 90 90 300K 745 248(3) 3.75 4.10 90 90 90 300K 750 150(5) 3.84 4.36 90 89.216 90 300K 750 150(5) 3.84 4.36 90 90 89.142 300K 662 133(5) Training data was generated using the Vienna Ab-intio Software Package (VASP) 135,136 . The SCAN meta-GGA exchange functional 137 was used, and electronic states were computed within the projector augmented wave method 138,139 . Pb 6s 2 6p 6 , Ti 3s 2 3p 6 3d 2 4s 2 , and O 2s 2 2p 6 electronic states were treated explicitly. The gamma point was used to sample the Brillouin zone. During training, energy, forces, and stress tensor were included in the loss function: 𝐿 = 𝑝 5 𝑁 QSG[5H 8 (𝐸 3 ''uvw −𝐸 3 uvw ) ( ' +,-./0 3%. + 𝑝 Q 3𝑁 QSG[5H 𝑁 GAx[H 8 8 8 G𝐹 3*J ''uvw −𝐹 3*J uvw H ( M J%. ' -12.0 *%. ' +,-./0 3%. + y 0 /' +,-./0 ∑ ∑ G𝜎 3* ''uvw −𝜎 3* uvw H ( / *%. ' +,-./0 3%. (2.70) Training was performed on 95% of the data set and 5 % was left out of training to be used as a test set for the model. The RSME for training (black) and test (red) sets are plotted in figure 2b-d for the respective components of the loss function in equation 2.70. Addition of stress tensor in 39 training was critical as accurate mechanical properties of the PTO structure were required for investigating the lattice dynamics under strain. We adopted an adaptive training scheduling scheme, where at first the forces and stresses were highly weighted then gradually turned down throughout training. To evaluate the model in structural properties, we first compared partial radial distribution functions at the experimental lattice constants at a temperature of 300 K, which are shown in Fig. 2, e-j. These figures demonstrate excellent agreement between the ground-truth DFT and trained NNQMD results. In addition, we examined bulk mechanical properties by computing the lattice constants, bulk modulus K, and elastic constant tensor 𝐶 *, ,FJ . For tetragonal PTO, there are six independent elastic constants C11, C12, C13, C33, C44 and C66 using Voigt notation 140 . Table 2 compares the NNQMD values with those of DFT calculations using SCAN and LDA exchange- correlation functionals (DFT-SCAN and DFT-LDA, respectively) as well as experimental values. 141,142 Elastic constants were computed using finite differences of stress tensor with respect to strain. DFT values were computed on a unit-cell with 16×16×16 k-point sampling with 550 eV energy cutoff. We see DFT-SCAN gives much better agreement with experiment in terms of mechanical properties, without suffering from large over-estimation in the tetragonality often seen in GGA-type functionals 143 . In-turn the NNQMD computed values show good agreement with experiment and the training data, indicating the model should be reliable in describing mechanical properties of PTO. a (Å) c/a C11 C33 C44 C66 C12 C13 K DFT- LDA 3.86 1.04 296.3 113.7 62.1 104.9 111.3 92.1 138 DFT- SCAN 3.84 1.125 231.7 69.4 59.3 105.3 99.9 71.8 113 NNQMD 3.87 1.095 276.7 55.5 58.4 95.1 117.0 85.2 128 Exp 8,9 3.90 1.064 237, 235 60, 105 69, 65.1 104,144 90,101 70, 98.8 104 Table 2. Computed and experimental lattice constants for PTO. Elastic constants Cab and bulk modulus K are in GPa unit. 40 2.3.2 NNQMD for Excited state PTO We performed Fermi-occupation quantum molecular dynamics (FOQMD) with electronic temperature informed from informed from RT-TDDFT to describe the effects of massive electronic excitations on atomic motions PTO. This generated the training set for NNQMD which consisted of of 3×3×3 and 4×4×4 unit cells and each system is thermalized at 300 K using NVT ensemble. A plane-wave basis was used with a cutoff energy of 410 eV for the wave functions and 3,400 eV cutoff for the charge density. Vanderbilt style ultra-soft pseudopotentials (USPPs) were used, and local density approximation (LDA) was used for the exchange-correlation. Only the gamma-point was sampled from the Brillouin zone. Atomic trajectories were obtained using the velocity-Verlet algorithm at time step of 50 a.u. (1.209 fs). In total, 500 frames of QMD trajectory were used. 95% of the total frames were used as training set and 5% frames were so for test. QMD simulations were carried out using the highly parallelized plane-wave based QXMD software 133 . Fig. 3 shows the training and test performance for the excited-state NNQMD model, and Table 3 summarizes RMSE and correlation coefficient between model prediction and the ground truth after 500 epochs. Here, the same set of feature vector and hyperparameters are used in the ground-state model. We have obtained well-converged RMSE and high correlations for the system energy and atomic forces for the excited state NNQMD model. Figure 3. Training (black) and test (red) performance for the excited-state NNQMD model. (a) RMSEs of the system energy (A) and aggregated atomic force (B). Table 3. Excited state NNQMD Model Training Performance and Correlation Coefficient. 41 RMSE Correlation Coefficient Energy (meV) Force (eV/Å) Energy (%) Force (%) Train Test Train Test Train Test Train Test 4.203 3.591 0.1178 0.1192 97.46 97.01 99.03 99.00 Figure 4. Comparison of excited state NNQMD and QMD radial distribution functions. Radial distribution functions g(r) of excited-state QMD (black) and its NNQMD model (red) up to 6 (Å) cutoff distance. The systems are thermalized at 300 K using NVT ensemble. Overall, QMD and NNQMD show excellent agreements in their peak position, height and width of g(r). Due to the photo-induced structural change, the fine structures originated from the tetragonality of PTO crystal disappeared. Figure 5. Comparison of Excited State NNQMD and QMD bond angle distributions. Bond angle distribution functions of ground state QMD (black) and its NNQMD model (red). O-Pb-O, O-Ti-O and O-O-O angles distribution are shown with bond cutoff distance of 3Å. 42 The excited-state NNQMD model is validated by comparing the pair distribution function g(r) and bond angle distribution against ground truth QMD result. We obtain excellent agreement in g(r) and bond angle between QMD and NNQMD (Figs. 4 and 5). Due to the photo-induced FE reduction in PTO crystal, the signatures of the off-center Ti atom position disappear in both g(r) and bond angle distributions. Emergent structural features due to the octahedral rotation, such as the split in the O-O-O peak around 120 degree, are successfully reproduced by NNQMD model. 2.3.3 MM and NN/MM Model Development To generate a molecular mechanics (MM) model for SrTiO3 (STO) for hybrid NN/MM simmulations which will be discussed in chapter 4, we used the Lennard-Jones (LJ) interaction, 𝑉G𝑟 *, H=4𝜖ÂO Y -3 S*, P .( −O Y -3 S*, P / Ã (𝑎,𝑏 ∈ {Sr,Ti,O}). (2.71) Parameters 𝜎 G! were fit to match the experimental lattice constant, and the interaction strength 𝜖 was fit to bulk modulus for STO. The model parameters are summarized in Table 4 and the computed lattice constants, bulk modulus K, and elastic constant tensor 𝐶 *, ,FJ for the model are summarized in Table 5 along with their comparison to experiment 144 . For cubic STO, there are 3 independent elastic constants C11, C12, and C44 in the Voigt notation 140 . sab (Å) Sr Sr 3.5532 Sr Ti 3.0771 Sr O 2.5125 Ti Ti 3.5532 Ti O 1.7766 O O 2.5125 Table 4. Optimal Lennard-Jones parameters sab. 𝝐 = 7.2x10 -2 (eV). 𝑎(Å) 𝐾(GPa) 𝐶 44 (GPa) 𝐶 45 (GPa) 𝐶 66 (GPa) Exp 144 3.905 174 317 103 124 MM 3.906 169 321 105 104 Table 5. Experimental and computed lattice parameters and elastic constants. We performed hybrid NN/MM simulations in a core+buffer approach analogous to divide- and-conquer approaches commonly used in linear-scaling QMD simulations. Here, a NN region 43 is defined with a core radius 𝑅 z and additional buffer radius 𝑅 { (>𝑅 z ). The forces 𝑭 on atoms from the two models are then determined by a linear support function 𝛼(𝑟) that forms a partition of unity over the entire simulation box: 𝐅 * (𝑟)=𝛼(𝑟)𝐅 * ''uvw (𝑟)+(1−𝛼(𝑟))𝐅 * vv (𝑟) (2.72) 𝛼(𝑟) =Í 1 (𝑟 <𝑅 z ) − . m 7 ;m 8 (𝑟−𝑅 { ) (𝑅 z ≤𝑟 <𝑅 { ) 0 (𝑟 ≥𝑅 | ) (2.73) Forces on the atoms are determined by the NNQMD model inside the core NN region within 𝑅 z , linearly mixed in the buffer region (𝑅 z ≤𝑟 <𝑅 { ) and purely MM outside. The NN region and corresponding support function 𝛼(𝑟) can be defined in 1, 2 or 3D depending on the geometry of the problem. For our simulations RC was set as the distance of the middle of PTO region to the start of the STO region. RB was set to RC plus the neural network cutoff distance. 2.3.4 NNQMD for Model for NH3 Vibrational Properties As discussed in chapter 1, I developed an NNQMD forcefield to perform large scale PIMD simulations of ammonia to understand the role of nuclear quantum effects on its vibrational properties. These NNQMD simulations were performed utilizing the recently developed group-theoretically equivariant neural-network forcefield model Allegro 128 . Allegro falls under the class of graph-based models where each atom i represents a node in the graph and edges are the interatomic distances 𝑟 *, , while retaining data locality of features to achieve high computational speed 126 . For each node i the model is designed to predict the edge/pair energies 𝐸 *, , and the total energy is the sum over all pair energies for each node. Forces are then computed through gradients of the total energy. 44 Based off the work of the state of the art neural equivariant interatomic potential (NEQUIP) 127 , Allegro utilizes equivariant features in the form linear embeddings of spherical harmonic projections of atomic pair wise atomic distances (𝑟 *, ). The model is “equivariant” in that only tensor operations on spherical harmonic projectors that respect symmetry of the Euclidean group E(3) are performed. Allegro also allows for scalable and faster computations in comparisons to previous works as the model is local with spherical harmonic embeddings only performed on neighbors within a cutoff distance. For further details the reader is directed to the methods and to Ref. 128. Training for the model was done using configurations from solid and liquid DFT-MD simulations allows as configurations from the PIMD simulations. Training for the Allegro neural network force-field was done using frames from DFT- MD simulations for the liquid (205 K) and solid phase(60 K), and PIMD simulations. These were performed on supercells containing 108 NH3 molecules (432 atoms) using the optPBE- vdW 145 functional with box lengths of 16.2Å and 15.39Å respectively. A distance cutoff of 5.0Å was used to train the model. Allegro implements equivariant features in the form of spherical harmonic projections of atomic pair wise atomic distances (𝑟 *, ). Allegro utilizes two latent spaces one for equivariant tensor features (𝑌 F[ ), and one for scalar invariant features 𝑟 *, ,𝑍 * ,𝑍 , , with 𝑍 * being the atomic number of species i. The equivariant latent space is updated through equivariant tensor products up to a maximum order l with a weighted embedding of spherical harmonics projections of atoms within the central atoms cutoff radius. The weights are determined by multi-layer perceptron acting on the scalar latent space features. The equivariant features of the same irreducible representation are linearly mixed, and then the scalar latent space is updated through tensor contraction of scalar and new linearly mixed equivariant features which is then passed through a multi-layer perceptron. This process is repeated for N layers and 45 then passed through a final multi-layer perceptron for the pair energy prediction. For more details on the model, we refer the reader to reference 128 and https://github.com/mir- group/allegro. 2.4 Path Integral Molecular Dynamics In all the above simulation techniques atomic nuclei were treated as classical particles. This approximation breaks down in the presence of materials with small mass elements (i.e hydrogen) under high-pressure and/or low temperature conditions. The incorporation of nuclear quantum effects for the nuclei can be treated semi-classically via path-integral molecular dynamics. If we assume the nuclei are Boltzmann particles, which allows us to neglect there exchange-interaction, we can write the quantum partition function for the system as the following limit 122–124 : 𝛧= lim 9→; DD(𝑀 < 𝑃/2𝜋𝛽) = 5 J𝑑𝑅 < > exp (−𝛽PP 1 2 𝑀 < 𝜔 9 5 (𝑅 ? > −𝑅 ? >@4 ) 5 +𝑉(𝑅 4 > , …𝑅 < > ,…,𝑅 A > ) A <B4 9 >B4 A <B4 9 >B4 ) (2.73) With 𝜔 } ( = } # % and 𝛽 =1/𝑘 | 𝑇. P here represents the number replica particles that are coupled together harmonically with frequency 𝜔 } with 𝑅 * } r. =𝑅 * . . Thus the path integral can be sampled from replica molecular dynamics simulations, with observable expectation values being computed as there average over the replica simulations. The more “quantum” the simulation is the more replicas are needed for convergence of computed observables. 46 Chapter 3 : NAQMD simulations for development of polymer dielectrics with high dielectric break-down fields for high field capacitor applications In this chapter I will discuss using NAQMD simulations to study hot carrier dynamics in polymers. In this part of the chapter, I will discuss a systematic study of the hot carrier dynamics in polyethylene under different field strengths where I find a polaronic localization transition near the breakdown field. Such polaronic localization transition has potential for critically missing prediction proxy for computationally screening dielectric polymers with high breakdown fields. In the second section of the chapter I will discuss using NAQMD simulations in conjunction with surface potential decay measurements performed at Uconn to illustrate the ability of boron nitride nano-coats to trap hot carriers and increase breakdown strength of polymers. Searching for materials with similar deep-well potential profiles we could lead to a computationally efficient way to determine good polymer coatings that can mitigate breakdown. 3.1 Studying field dependance of hot carrier dynamics in polyethylene (PE): 3.1.1 Simulation Results To systematically study the effects of electric fields on the hot carrier dynamics in PE precluding dielectric breakdown, we preformed NAQMD simulations on a PE slab consisting of 18´1´2 unit cells composed of 36 PE chains. Figure 6a shows (001) (or xy) plane view of the PE slab. PE chains are orientated along the [001] direction and are numbered from 1 to 36 with respect to their position along the [100] direction. An additional 15 Å of vacuum was applied along the [100] direction on both sides (total 30Å) as shown in Fig. 6a. Electric fields were applied in the [100] direction (left to right in Fig. 1a) using a sawtooth potential. Figure 6b shows the spatially projected partial density states (PDOS) on each polymer chain along the [100] direction for PE under a 100 MV/m electric field. In the PDOS we observe a superposition 47 bending due to Figure 6. Polyethylene (PE) under electric field. (a) (001) (or xy)-plane view of the PE slab, where PE chains are numbered from 1 to 36 with respect to their positions along the [100] direction. (b) Partial density of states (PDOS) projected onto the 36 PE chains under an electric field of 100 MV/m. Schematic of excitation is shown, where blue and red circles represent hole and electron, respectively. Blue and red arrows show the respective direction of travel of hole and electron. (c) The response of the local Kohn-Sham (KS) potential to the applied electric field. Blue and 48 red curves show the KS potential profile along the [100] direction with and without electric field, respectively, whereas the black curve shows their difference. the exposed surface 146,147 and tilting due to the linear voltage drop arising from applied electric field as shown in Fig. 6b. The resulting PDOS is tilted near the surface on the right edge of the supercell due to the applied positive electric field, while the left edge displays a flat band profile due to the mixture of band tilting and linear voltage drop. The linear voltage drop can be seen more clearly by examining the change of the local Kohn-Sham (KS) potential profile due to the applied electric field. Figure 6c compares the (100) plane-averaged local KS potential with an applied electric field of 100 MV/m (blue curve) with that without electric field (red curve). Their difference (black curve) represents change in potential due to the applied electric field. The local KS potential responds linearly to the applied electric field. To simulate carrier excitation, an electron is excited from the valence band in the bulk of the material to the conduction band, thereby creating an electron-hole pair. The hole behaves as a positively charged particle and thus moves towards x = Lx, indicated with the blue arrow in Fig. 1b, where Lx is the length of the supercell along the [100] direction. On the contrary, the excited electron moves towards x = 0 as schematically shown by the red arrow in Fig. 6b. Since the band profile near x = 0 is relatively flat as mentioned before, the electron experiences relatively less driving force by the electric field, and thus the hole dictates the breakdown processes. To investigate the effect of different electric fields, Fig. 7, a-c, shows time evolution of the KS eigenvalues for fields of 100, 600 and 900 MV/m, respectively. After excitation at time t 49 = 0, an electron-hole pair is generated. As a result, an unstable high-energy state is formed, Figure 7. Hole relaxation and localization. (a)-(c): Time evolution of KS eigenvalues under electric fields 100 (a), 600 (b) and 900 (c) MV/m, respectively. The blue and red curves are the hole and electron eigenenergies, respectively. (d)-(f): Center-of-mass (COM) position of the hole, 𝑋 C DEFG , along the [100] direction for 100 (d), 600 (e) and 900 (f) MV/m, respectively. (g)-(i) Chain participation number, CP, of the hole for 100 (g), 600 (h) and 900 (i) MV/m, respectively. causing the hole energy (blue curve) to relax to the valence band maximum (VBM) while the excited electron energy (red curve) to relax to the conduction band minimum (CBM) rapidly within 50-100 fs, as shown in Fig. 7, a-c. Under a 100 MV/m field the hole undergoes rapid and decoherent motion and does not localize to any particular electronic band or spatial location. The 50 rapid oscillation in the hole dynamics, shown in Fig. 7a, was found to result from high-frequency transitions between highly degenerate KS eigenstates. For E = 600 MV/m, this high frequency motion is still present, but after relaxation the hole localizes to the VBM from t = 125-155 fs and after t = 190 fs. During these periods the VBM wavefunction forms a gap state as it separates from the rest of the KS states. Under the application of a 900 MV/m field, a transition occurs from decoherent hole motion to coherent motion occurs, as the hole smoothly relaxes to the VBM with a clear dominant frequency. This is then followed by the immediate formation and sustainment of a gap state at t = 40 fs. Fourier transformation of the hole eigenenergy between 0 and 40 fs shows a peak near 200 THz, which is close to twice the C-H stretching frequency (92 THz). The second C-H harmonic reflects the deviation of the C-H bond length from its equilibrium values, which oscillates at twice the C-H stretching frequency as was found in a previous study of quaterthiophene on zinc-oxide surface 148 . The same harmonic is present in the hole time evolution under a field of 1200 MV/m. Similar C-H harmonic resonance may also be present in lower electric fields, but is washed out by the high frequency of electronic transitions due to KS state degeneracy. To understand the formation of these gap states, we further investigated the center-of- mass (COM) position of the hole along the [100] direction as a function of time. For the k th KS wavefunction Ψ J (𝒓,𝑡), COM position can be calculated as 𝑋 U J = ~ H (f Im ³LogO∫ 𝑑 M 𝑟 | i Ψ J (𝒓,𝑡)| ( expO (f* ~ H 𝑥PP´ , (3.1) where 𝐿 P is the length of the supercell along the [100] direction and Ω is the volume of the supercell. We also calculated the chain participation number (CP) of hole for each applied electric field: CP = ∑ 1/𝑝 * ( * , where 𝑝 * is the probability the hole is on the i th PE chain. Figure 7, d-f, shows the time evolution of the hole COM position, whereas Fig. 7, g-i, shows that of the CP 51 value. The 100 MV/m field does not induce any gap states, suggesting a highly delocalized hole in the system as shown in Fig. 7, d and h. The gap states formed under 600 and 900 MV/m fields correspond to spatial localization of the hole at the surface and its confinement to a single PE chain, indicated by the CP value near unity. In addition, the hole’s physical motion to the surface appears to undergo a similar coherent motion found in its eigenenergy evolution. We found a peak frequency near 200 THz in the carrier motion, which is consistent with earlier discussed C- H vibrational resonance in the hot carrier motion. This coherent frequency was also seen in PE under a 1200 MV/m field. The presence of this frequency in both the energy and spatial evolution of the hole indicates that C-H vibrational modes assist the carrier motion to the surface under these high electric fields. This transition from decoherent to coherent motion under increased electric fields indicates a critical transition occurring past 600 MV/m, which is near the reported breakdown field for PE. Experimental short-term breakdown of polymers is expected to follow a statistical distribution 149 , where the average experimental value for intrinsic breakdown is reported in the range 550-650 MV/m 10,150 . With increased electric field, we observed decoupling of the VBM eigenstate from the remaining KS states prior to hot-carrier relaxation, which displays an energy profile similar to the gap states formed when the hole localized on the surface of the slab. This indicates that high electric field induces localized states at surface of the PE slab before the hot- carrier relaxation. To quantify this localization transition, we examined the spatial position of the VBM wavefunction along the [100] direction for the first 50 fs of the simulation as shown in Fig. 8, a- c. As the field increases, the VBM wavefunction becomes increasingly localized to surface of 52 the slab. Figure 8. Field Induced Localization. (a)-(c) show the expectation value of the position VBM wavefunction along the [100] or x direction during the first 50 fs for fields of 300,600, and 900 MV/m respectively. (d) shows the energy fluctuation of the VBM wavefunction versus the electric field. By inducing a localized state at the surface of the slab, high electric fields increase the susceptibility for the hole to localize on the surface upon relaxation. For PE under an electric field of 1200 MV/m, the localization of the VBM wavefunction was delayed. This resulted in a delay of the emergence of the coherent carrier motion and in formation of a sustained gap state). A critical transition occurs past 600 MV/m, where the VBM wavefunction becomes extremely localized to the surface and decoupled from the other states, followed by coherent motion of the 53 hole to the surface. This transition is well quantified by the fluctuation (σ ) of the VBM energy before hole localization as a function of applied field (Fig. 8d). Figure 8d exhibits a phase- change like transition, where 900 and 1200 MV/m fields extremely dampen the VBM energy in comparison to the lower fields. This extreme dampening corresponds to a highly localized state at the surface decoupled from the remaining KS states. The degeneracy of KS states is broken by the extreme electric fields, resulting in coherent motion of the hole to surface followed by sustained localization. The extreme localization of hole on surface and formation of gaps state results in the chemical damage to the PE slab. Figure 9 shows the C-C bond lengths for PE chains 35 (solid curves) and 36 (dashed curves) as a function of time under all five applied electric fields. For a field of 100 MV/m, no localized gap state is formed and accordingly no bond stretching or breaking occurs at the surface. For 300 MV/m, two gap states were formed briefly at t = 90-115 fs and t = 195-210 fs resulting in a localization of the hole on polymer chain 35. These gap states resulted in weakening of C-C bonds henceforth increased C-C bond length on the 35 th chain. The 600 MV/m field induced the formation of localized gap states at t = 125-fs on polymer chain 35 leading to more extreme bond-stretching. Further at t = 200 fs, the hole localized on chain 36 for a longer time (50fs) leading to a C-C bond cleavage. Past the critical transition shown in Fig. 8d, we observed formation of a highly localized gap sate immediately after hole reaches the surface resulting in rapid destruction of C-C bonds at the surface for PE under 900 and 1200 MV/m fields as shown in Fig. 9, e and f, respectively. Electric field has a twofold effect on the chemical damage in PE by first inducing a localized occupied state at the surface that increases the susceptibility of the hole to localize at the surface after relaxation. After hole localization to the 54 surface, surface C-C bonds become depopulated making them more Figure 9. Bond stretching and breaking. (a) Three distinct bonds colored black, red, and blue present in PE chains. Dark gray and white spheres represent carbon (C) and hydrogen (H), respectively. (b)-(f): Time evolution of C-C bond lengths for surface chains 35(solid curve) and 36 (dashed curve). Black, red and blue curves represent same colored bond in (a). 55 susceptible to being broken by the electric field. Once broken, the resulting defects will change the charge dynamics inside the polymer triggering a positive feed-back breakdown process 151 . While weaker fields here caused limited chemical damage, long-term exposure to them could eventually result in breakdown due to the process known as polymer aging. Experimental evidence suggests chemical aging may occur as a result of energy release due to electron-hole recombination on time scales out of reach for NAQMD simulations 152,153 . 3.1.2 Discussion In summary, we have performed an exhaustive investigation of effects of electric fields in PE. We observed that high electric fields lift the degeneracy of the PE slab by inducing a highly localized state at the slab surface. This increases the susceptibility for the hot carrier to localize during carrier relaxation. Upon localization, a high energy gap state will form resulting in bond stretching and breaking. We found a critical transition in the response of PE to the applied electric field occurring after 600 MV/m, which is near the known intrinsic breakdown field for PE. This transition is marked by an extreme damping of the VBM energy fluctuation and the presence of C-H vibrational resonance in the hot carrier motion. The polaronic localization transition found here may provide a critically missing prediction method for computationally screening dielectric polymers with high breakdown fields. 3.2 Mitigating Hot Carrier Damage with hexagonal Boron-Nitride(hBN) Nano-Coats. 3.2.1 Simulation and Experimental results I next investigated utilizing NAQMD simulations to explain experimentally observed the carrier- trapping effect in an hBN-coated Kapton film by a surface potential decay (SPD) measurement. In the SPD measurement the polymer film was charged under a high voltage of 6 kV, and the surface potential was monitored by a Kelvin probe. Compared with uncoated Kapton, coated Kapton delivers a higher initial value of surface potential, indicating more carriers are trapped 56 due to the existence of hBN nanocoating, as shown in Figure 10a. While the surface potential decays with time in both samples, it could be found that the surface potential of the coated Kapton is always higher than that of the uncoated one, implying such carrier trapping effect could play a significant role in controlling the charge carrier dynamics. To further asses the effects of hBN to trap hot carriers, we calculated the trap energy based on the theory of isothermal surface potential decay (ISPD) 154 . As seen in Figure 10b, a larger trap density can be found after coating BN for both shallow and deep traps. This further demonstrate the charge trapping effect of BN. As surface potential decay rates can be highly sensitive to the test setup 155 , a detailed explanation of the test setup is provided in methods section. Nevertheless, the relative decay rates of uncoated and coated Kapton are expected to be similar regardless of test setup. The dielectric insulation of polymers is fundamentally dictated by charge carrier injection and transport. It is believed that the electrical conduction could be suppressed when charge carriers are trapped in deep wells, thus contributing to enhanced insulation properties. To investigate this, we performed dielectric breakdown measurements for coated and uncoated Kapton films. Figure 10c demonstrates that the breakdown strength of coated Kapton is higher than the uncoated polymer, which is due to the carrier trapping effect induced by the hBN Figure 10. (a) Surface potential decay and (b) trap energy, and (c) Weibull breakdown strength of Kapton and hBN coated Kapton. 57 nanocoating. In scanning electron microscope (SEM) image of the hBN coated Kapton is we could see the uniform distribution of BN coatings on the Kapton at a very large scale, and no obvious aggregation was seen. Similar to the surface potential decay, the exact value of the breakdown strength is highly associate with the setup of the experiment such as the ramping rate 156 . A detailed description of the experimental setup is provided in the methods section. To gain mechanistic understanding of how hBN nano-coats could enhance polymer breakdown, we performed static density functional theory (DFT) and NAQMD simulations to characterize the hot carrier dynamics in an electrode-hBN-polymer interface quantum mechanically. Polyethylene (PE) and aluminum (Al) were chosen as standard polymer and 58 electrode materials to assess the effectiveness of the hBN coating, as their interface is best characterized both experimentally and computationally. 19,20,26–28,157–160 As standard electrode/polymer databases use PE/Al as a baseline, with which other polymers and electrodes are assessed 26,27,29,158–160 , they can thus provide a practical computational-screening platform for determining coating materials in the future. While the band-alignment varies for different hBN- coated polymers, charge trapping by hBN is expected to remain operative. Thus, PE provides a computationally efficient benchmark for characterizing effective coating materials. Figure 11 .Optimized structures and local potential profiles. (a) Optimized Al, hBN and PE structures used to compute the local potential and band alignment. Al, B, N, C and H atoms are colored in silver, blue, yellow, gray and white, respectively. (b)-(d) Laterally (xy) averaged local potential profiles along the surface-normal (z) direction. 59 To understand the potential charge injection, we first computed the local potential profiles and band alignment of aluminum, hBN and PE using DFT 2,120,161 . To compute the potential profiles and band energies, each material was optimized separately with 30 Å of vacuum applied along the z axis (c crystallographic axis) with the Perdew–Burke–Ernzerhof (PBE) 162 styled generalized gradient approximation (GGA) within the projected augmented wave-vector method (PAW) 138 . Full details of the electronic structure calculation are described in the methods. The local potential energy, valence band maximum (VBM), conduction band minimum (CBM), and Fermi-level energies were then computed with respect to the mean vacuum energies for each system. Periodic boundary conditions (PBC) were applied along all 3 crystallographic directions. The optimized systems are illustrated in Figure 11a. For aluminum, structural optimization was performed for a supercell consisting of 2×3×5 symmetrized 4-atom cubic unit cells (instead of a 1-atom trigonal unit cell). For PE, a 2×5×4 supercell was optimized, where the polymer chain direction is along the y axis (b crystallographic axis) resulting in 8 total PE chains. For hBN, a 72-atom bilayer system (18 unit cells per layer) was optimized. Figure 11, b-d, shows the laterally (i.e., xy-plane) averaged local potential profiles along the surface-normal (z) direction for the optimized structures with the mean vacuum energy normalized to zero. A clear Stark deep potential well is formed in the hBN-bilayer in comparison to aluminum and PE, indicating the potential for hot carriers to be trapped in hBN layers in an Al-hBN-PE interface. Figure 12a shows the computed band alignment for the three structures. PE-hBN forms a straddling style band alignment and the aluminum Fermi level is between the hBN VBM-CBM gap. This indicates that high electric fields inject charge from aluminum to hBN in an Al-hBN- PE interface. While the band gap is underestimated with PBE functional for both hBN and PE, 60 their relative values both trend with experiment 163,164 . Figure 12b illustrates the typical breakdown experiment that we model with NAQMD simulation. Field-induced hot carriers can either be: (I) injected directly from the electrode interface into the polymer, or (II) introduced through secondary excitations in the bulk of the polymer. Here, we study the hot carrier dynamics described by process I at an hBN-PE interface, whereas process II has extensively been studied using NAQMD simulations 20,26 . To study injection of hot carriers from Al on a hBN-PE interface we preformed NAQMD simulation (see methods) on interfacial structure along the (110) plane with 2×5×2 unit-cells of PE sandwiched between two separate 72-atom bilayers of hBN, which is visualized in Figure 4a. A vacuum layer of thickness 15 Å was added to each end of the interface for a total of 30 Å separation to remove spurious interactions between periodic images of the slab. Aluminum was not included in the interface as it merely supplies the initial condition of injection of electrons and holes into the hBN layers, which is modeled in NAQMD as the generation of a high energy electron hole pair at the edges of the hBN-PE slab as illustrated in Figure 13a. Figure 12. Band alignment and schematic of charg- injection experiment. (a) Computed band alignment with the PBE functional for Al, hBN and PE with VBM, CBM, Fermi-energy levels labeled. (b) Schematic of breakdown experiment, where hot carriers can either be: (I) directly injected into the polymer at electrode interface, or (II) caused by excitation of secondary electron-hole pairs in the bulk of the polymer. 61 To simulate the field-induced hot carrier dynamics, a sawtooth potential was applied along the vacuum direction to simulate the presence of an electric field. An electric field strength of 600 MV/m was used, which is near the breakdown field for PE 10,150 . Figure 13b shows the local potential profiles of the slab with and without the applied electric field and the difference between them. There is a clear linear shift in local potential as a result of the applied sawtooth potential and the Stark deep potential wells of the hBN layers remain. After examining the effect of the field on the local potential, NAQMD simulation under electric field was run within the canonical ensemble at a temperature of 300 K with a time step of 0.2418fs. Figure 13. NAQMD simulation. (a) hBN-PE interface used in the NAQMD simulation. B, N, C and H atoms are colored mauve, lime, black and white, respectively. Initial excited electron and hole Kohn-Sham (KS) wave functions are colored red and blue, respectively. Electric field was applied along the interface-normal direction (z) using a sawtooth potential. (b) Laterally (xy-plane) averaged local potential with (red curve) and without (blue curve) the applied electric field. The difference displaying the linear voltage drop is colored in black. (c) Time evolution of the KS eigenvalues with the hole and electron energies colored in red and blue respectively. (d) Semi- log plot of the excitation energy vs. time (open circles) to estimate the relaxation rate by linear fit (dotted line). 62 Figure 13c shows time evolution of the Kohn-Sham (KS) eigenvalues along the non- adiabatic trajectory with the hole and electron energies colored in blue and red respectively. The excited hole and electron both slowly relax towards the VBM-CBM energy levels respectively (~200 fs). To better quantify the energy loss rate k, we calculated the average electronic excitation energy as a function time, Eexc(t)= ∑ Q & (A) & (A)( & (A)) & ;∑ (.;Q & (A)) & (A)(; & (A)) & ∑ Q & (A) & (A)( & (A)) & r∑ (.;Q & (A)) & (A)(; & (A)) & (3.2) where 𝑓 * (𝑡) is the occupation of the i th band at time t, 𝜖 * (𝑡) is the KS energy, and 𝜃 is the heavy- side function. The rate k was determined by a semi-log fit of the excitation energy as function of time over the exponential decay period of the energy (~200 fs) which is plotted in Figure 13d. The rate k = (2.26 ±.01) ×10 12 s -1 is orders-of-magnitude slower than reported in previous NAQMD study of sole PE 20 , which is indicative of weak carrier-phonon interaction. This indicates that introduction of the hBN layer could limit strong polaronic coupling seen to result in rapid damage to C-C bonds in previous study 26 . As the current study incorporated a much thinner PE layer than previous stuide, tunneling effects could play a role in the calcuated rate. However, as tunneling is anticpated to cause faster carrier relaxation and surface localization in this sudy of an hBN-PE interface, the calculated slower rate indicates that carrier trapping effect of hBN is far stronger than potential tunneling effects. 63 Figure 14. Carrier participation dynamics. (a-g) Angular-momentum (l)-dependent fractional contribution of each atomic orbital type to the hole wave function. The hole wave function hops between the C(p) and H(s) PE states and B(p) and N(p) hBN states before localizing to hBN. (h-n) l-dependent fractional contribution of each atomic orbital type to the excited electron wave function, which remains primarily localized in hBN for the entire simulation. We next examined the carrier dynamics by looking at the angular-momentum (l) dependent atomic orbital contributions to electron and hole wave functions as a function of time via Mulliken analysis 165 (Figure 14). Figure 14, a-g, illustrates the hole wave function alternates between being composed of C(p) and H(s) orbitals and being composed of B(p) and N(p) during the first 100 fs. This indicates that the hole wave function is hopping between the hBN and PE layers for the first 100fs. After the first 100 fs, the hole wave function is largely composed of N(p) orbitals, with small B(p) orbital composition, indicating that the hole localized and became trapped in the hBN layers for remainder of the simulation. The fractional excited electron composition is similarly plotted in Figure 14, h-n. The excited electron is almost solely composed of hBN states (primarily B(p) orbitals) throughout the entire simulation. This indicates that the electron remains trapped in hBN deep potential wells for the entire NAQMD trajectory and unable to damage the polymer. 3.2.2 Discussion 64 Surface potential measurement implies that hBN nanocoating could contribute to charge carrier trapping, and consequently improve the dielectric insulation of polymers, which is evidenced by the enhanced breakdown strength in hBN-coated Kapton. NAQMD simulations corroborate the ability of hBN layers to trap hot carriers at a hBN-polymer interface. The trapping of hot carriers is a result of deep potential wells formed in the hBN layers. Since local-potential profiles relative to the vacuum energy can be efficiently computed with standard DFT calculations, probing for this feature in other potential coating materials could provide efficient rational design principle for general screening of effective polymer coating materials. For assessing the effectiveness of coat for specific polymers it is also desirable to compute the compute the well potential profile of the interface and compare the relative depths, especially since charge transfer and polarization effects can happen at the interface with polymers containing nitrogen and oxygen groups. In addition as previously mentioned, more complex polymers are expected to have different band structures which could potentially change the charge carrier dynamics. In preliminary theoretical and experimental study of hBN coated polycarbonate, however, we have found a similar charge trapping effect of hBN despite the differing band- structures compared to the polymers discussed here. The study of polycarbonate goes beyond the scope of this work and will be published elsewhere. While the key point of this paper is the deep potential located in the BN coating layer, other experiments such as photocurrent injection measurement can also provide complementary information such as barrier height between electrode and polymer as well as the surface states. Recent study of hBN using photocurrent injection measurements has indeed showed the ability of hBN nano-coats to reduce surface states at polymer-metal interfaces 166 , further corroborating the carrier trapping reported here. 3.2.3 Further Simulation and Experimental Details 65 Experiments: To prepare coating precursor, hexagonal BN nanosheets (US Research Nanomaterials) with a diameter of 800 nm was first dispersed in ethanol by tip sonication for 30 min. Polyvinyl butyral was then dissolved in the above solution by stirring overnight, acting as the polymer binder in the nanocoating. The solution was tip-sonicated again just before the coating process. Kapton film with a thickness of 8 μm was purchased from DuPont for the experimental study. The film was cleaned by ethanol and dried at 80℃. Then, it was dipped into the coating precursor and stayed for 30 s before being lift up. The dip-coated film was dried at 80℃ overnight before conducting electric measurement. Dielectric breakdown tests were carried out by a resistor- capacitor (RC) circuit with a linear voltage ramp of 215 V/s. Ball-plate electrodes were adopted to apply high voltage and electrically ground the sample. More than 10 samples of coated and uncoated films were tested to obtain a statistic breakdown strength in Weibull plot. In the surface potential measurement, a pair of needle-to-plate electrodes were used to conduct corona charge. The distance between two electrodes is 6 mm. The film sample was placed on the plate electrode that is grounded. A DC voltage of 6 kV was applied on the needle electrode for 1 min. Then the sample was transferred under a Kelvin probe to observe the surface potential decay for 10 min. Density Functional Theory (DFT) Calculation of Band Alignment: Structure optimization for band alignment calculation was performed till the forces were within 5×10 -4 Hartree/Bohr. Energies and forces were computed using a plane wave basis within the projected augmented wave-vector (PAW) method 138 . Projector functions were generated for the 1s hydrogen, 2s and 2p carbon, 2s and 2p boron, 2s and 2p nitrogen, and 3s, 3p and 3d Al states. For all structures a plane wave cutoff of 35 Ry for the wavefunctions and 300Ry for the density was used. The gamma point was used to sample the Birullion zone. The 66 Perdew–Burke–Ernzerhof (PBE) 162 styled generalized gradient approximation (GGA) for the exchange-correlation functional was used, and Van der Walls interactions were included with the DFT-D method. 167 More accurate calculations of the bandgaps and VBM-CBM energies with respect to the vacuum energy can be obtained with the use hybrid functionals such as the Heyd- Scuseria-Ernzerhof (HSE06) hybrid functional, incorporation of Hubbard interactions in the local density component of the exchange correlation potential (LDA+U), or even many body perturbative Green’s function methods such as G0W0 or GWΓ . 164,168–171 . However, we were mostly interested in obtaining qualitative understanding of the charge injection at the metal-hBN- polymer interface to inform NAQMD simulation of the hot carrier dynamics, where these advanced methods become impractical to use. Calculations were performed in the highly parallelized plane wave based software QXMD 133 . NAQMD Simulation: NAQMD simulation was performed on a hBN-PE interface along the (110) plane with 2×5×2 unit cells of PE sandwiched between two separate 72 atom bilayers of hBN. Prior to NAQMD simulation, hybrid cell and atomic optimization was performed on the slab until the forces were within 5×10 -4 Hartree/Bohr and successive energy minimization was within 1×10 -7 Hartree/atom. NAQMD simulation 172–174 was performed with excited state forces on atoms computed in the framework of time-dependent DFT (TDDFT) 5 with phonon assisted excited transitions between electronic states modeled within the surface hopping approach 175,176 . 67 Chapter 4 : Multi-Scale NNQMD simulations for optical, mechanical, and electrical control of emergent polarization topologies in in next generation ferroelectrics. In this chapter we I will discuss using multi-scale NNQMD simulations for opto- mechano-electrical control of ferroelectric topologies in PbTiO3(PTO) based ferroelectric structures, as discussed in the introduction. This first part of the chapter will study the effect of strong optical excitation and subsequent topological defect formation in PTO. For this I incorporate a multi-scale framework that combine Fermi-occupation quantum molecular dynamics (FOQMD) ,real-time time-dependent density functional theory (RT-TDDFT) along with Maxwell’s equations to describe light-matter interaction 67 ,surface hopping based non- adiabatic quantum molecular dynamics (NAQMD) to include non-adiabatic coupling between electrons and phonons 68 , thereby assessing the relaxation time between the electronic and lattice temperatures until which FOQMD is justified. FOQMD simulations are then used to train neural networks to perform large scale excited state NNQMD (XS-NNQMD) simulations to study effect of optical exciation on large nano-polarization topologies in PTO. The second part of the chapter will discuss using hybrid neural network and molecular mechanics simulations to study mechanical control of polar skyrmions in SrTiO3(STO)/PbTiO3(PTO) superlattices, and the final part will discuss a simple born-effective charge extension to NNQMD to study electrical control of flux closure domains in strained PTO. 4.1 Optical Induction of topological defects in PTO. 4.1.1 Simulation Results To study the effect of strong optical excitation on PTO, we first performed FOQMD simulations, in which massive electronic excitation is represented by electronic occupation numbers that follow Fermi-Dirac distribution with an effective electronic temperature 66 . Following previous experimental and theoretical studies of non-thermal structural changes under 68 strong photoexcitation 177–179 , 5% of the valence electrons were promoted to the conduction band. We also performed RT-TDDFT simulations 67 to estimate the corresponding laser fluence. We considered an external laser pulse with a wave number of 200 nm to excite electrons above the band gap of 3.88 eV in PbTiO3 (PTO) 65 . It is worth noting that recent experimental studies have shown 400 nm laser pulses to be sufficient to create electron-hole pairs in PTO-based heterostructures 38,65 , but it is likely due to the band alignment at hetero-interfaces 180,181 . Fig. 15A and 15B show time evolutions of the external vector potential for a 200 nm laser pulse and energy absorbed by the electronic system, respectively, for the case of 𝐄 b`c = 0.2 a.u. and 𝜏 b`c = 6 fs. Figure 15. Vector potential and absorbed energy. Time evolutions of the external vector potential due to a 200 nm laser pulse (A) and energy absorbed by the electronic system (B). The absorbed energy in Fig. 15B generates non-equilibrium electron occupation numbers, which are subsequently thermalized toward the equilibrium Fermi-Dirac distribution due to electron-electron interaction. Fig. S2 shows the occupation number of each KS orbital for both spins, 𝑓 H (𝑡)=∑ | ∫𝑑𝐫𝜓 HY ∗ (𝐫,0)𝜓 H I Y (𝐫,𝑡)| ( 𝑓 HOY HOY , (4.1) 69 as a function of the KS energy at times t = 0.7 fs and 2.4 fs. At t = 0.7 fs, we see a sharp peak (shaded red) corresponding to photo-excited electrons, along with a dip (blue shade) corresponding to holes, which are superimposed on top of a step-function-like Fermi-Dirac distribution. At t = 2.4 fs, theses sharp features are absorbed into broadened Fermi-Dirac distribution. These results signify rapid thermalization of electronic occupation numbers toward thermalized Fermi-Dirac distribution within femtoseconds, thus justifying the use of Fermi-occupation quantum molecular dynamics (FOQMD) simulations. Figure 16. Electronic thermalization. Occupation numbers of KS orbitals as a function of KS energies at time t = 0.7 fs (A) and 2.4 fs (B). The number of excited electrons from the valence band to the conduction band can be estimated by summing 𝑓 H (𝑡) over initially unoccupied orbitals s. We have simulated 200 nm laser pulses with duration 2 fs for various values of electric field. Fig. 17 shows the fraction of valence electrons that are excited to the conduction band as a function of the corresponding laser fluence. With a fluence of 0.15 J/cm 2 , 5% of the valence electrons were promoted to the conduction band, which corresponds to the FOQMD simulations. We should note that this estimate is for a perfect 70 PTO crystal without defect or trap states. In actual experiments, much lower fluence would be required to achieve the same excitation level due to these factors. Figure 17. Excitation as a function of fluence. Fraction of valence electrons that are excited as a function of laser fluence. To study possible charge transfer under excitation, static ground-state and Fermi- occupation density functional theory (DFT) calculations were performed (see Methods for details). Fig. 18A shows the ground-state partial density of states (PDOS). Oxygen states are near the valence band maximum (VBM), whereas states near the conduction band minimum (CBM) are composed of Ti states, indicating that photoexcitation will cause charge transfer to from O to Ti. Figure 18. Photoinduced charge transfer. (A) Atom projected partial density of states (PDOS). States near the valence band maximum (VBM) are primarily composed of O and states near the conduction band minimum (CBM) are primarily composed of Ti, indicating expected O-Ti charge transfer upon excitation. (B) Isosurfaces of hole (blue) and excited electron (red) densities upon excitation. The hole density is primarily localized around O atoms, while the excited electron density is localized around Ti atoms. Spheres with black, silver and red color represent Pb, Ti and O atoms, respectively. (C) Instantaneous force vectors as a result of excitation. 71 Fig. 18B illustrates the expected charge transfer by plotting the difference in electron density between the excited and ground states. The hole density (blue) is primarily localized around O atoms, whereas the excited electron density (red) is primarily localized around Ti atoms. This charge transfer results in an attractive force between O and Ti atoms along the c polarization axis, directly opposing the ground-state polarization which is illustrated in Fig. 18C. The forces on the other atoms remained nearly zero due to the symmetry of the TiO6 octahedra. To investigate the lattice dynamics under excitation, we performed FOQMD simulations in the iso-thermal iso-baric (NPT) ensemble at a temperature of 300 K on a 4×4×4 PTO supercell as well computed the ground- and excited-state phonon dispersions. Table 6 lists the average O, Ti, and Pb displacements along the c polarization axis for the first 60 fs and 500 fs of the NPT-MD simulation. The O atoms are divided into 3 categories O1, O2, O3 in accordance to their positions in the 5-atom PTO unit cell which is labeled in Fig 19A. Time Pb Ti O1 O2 O3 60 fs -0.04 Å 0.09 Å -0.19 Å -0.20 Å -0.14 Å 500 fs -0.14 Å -0.26 Å -0.67 Å -0.68 Å -0.52 Å Table 6. Average displacements of each atom type after 60 fs and 500 fs along the c polarization axis. Oxygen types are labeled O1, O2, O3 in accordance with their symmetry positions in the PTO unit cell which is pictured in Fig. 2A. The displacements during the first 60 fs along with the instantaneous force at time t = 0 illustrated in Fig. 18C indicate an initial activation of A1 TO1 like optical phonons that could potentially reverse or erase the polarization of the PTO crystal which is diagrammed in Fig. 19A. Average displacements after 500 fs indicate up-conversion into the higher frequency E TO2 mode which is also illustrated in Fig. 19A. Collective thermal excitations of A1 and E TO modes have been previously deemed responsible for negative thermal expansion in PTO through c axis contraction, and their softening has been related to polarization and tetragonality reduction in PTO 182,183 . 72 Figure 19. Bulk lattice and polarization dynamics. (A) Schematic of the initially activated TO1 phonon mode and the up-converted TO2 phonon mode. Oxygen atoms are labeled O1, O2, and O3 with respect to their symmetry position in the PTO unit cell. (B) Snapshot of oxygen rotations with magenta arrow labeling of the rotation orientation for 4 unit cells in ab plane 940fs after excitation. (C) ground and (D) excited state vibrational density of states (VDOS). VDOS was computed from Fourier transform the velocity auto-correlation function of the MD simulation at 300 K. (E) Ground-state and (F) excited-state phonon dispersion. Activated high-frequency optical branches are both softened and hardened to form a phonon band gap. (G) Excited-state polarization dynamics at 300 K. Average Ti polar displacements along the c axis are initially reversed, followed by a steady-state structure with near-zero average polarization. Visual inspection of the trajectory also indicates the presence of oxygen rotations. Fig. 19B shows a snapshot of 4 unit cells in the ab plane of the trajectory after 940 fs with the tilt orientation highlighted with magenta arrows. We also examined the phonon dynamics by comparing the vibrational density states (VDOS) of the ground state and excited state simulations, which are plotted in Figs. 19C and 19D, respectively. Fig. 19D illustrates further phonon up-conversion beyond the initial TO1 and TO2 phonon activation. The VDOS also appears to form a gap between 10-15 THz upon excitation as well as loss of higher frequency modes. To further examine this, we plotted the phonon dispersion for the ground-state and excited state, which are shown in Figs. 2E and 2F. Hardening and softening of high-frequency optical branches leads to the formation of a 73 large phonon band gap. The presence of a large phonon gap has previously been attributed to quantum paraelectricity in STO 184 , further indicating the potential photoinduction of paraelectricity in PTO. Phonon instabilities also occurred at the R and M points (manifested as negative modes), which is consistent with octahedral rotations seen in the MD trajectory 185,186 . To examine the polarization dynamics, we plotted the average polar displacements 〈𝑫 〉 of the Ti atoms along the c polarization axis (Fig. 19G). 〈𝑫 〉 was computed as the average difference of the Ti atoms from the centroid of the TiO6 octahedra. In a centrosymmetric non-polar structure, Ti atoms will be located at the center of TiO6 octahedra, while in the polar structure they are displaced. During the first 200 fs, we observe an immediate reversal of the average polar displacements, followed by relaxation to a quasi-steady state structure with near-zero average polarization along the c axis. For comparison, the polar displacements for ground-state QMD trajectory are shown in Fig. 20, for which the crystal remains negatively polarized along the c axis. Figure 20. Ti polar displacement dynamics. (A) ground state and (B) excited state polarization dynamics at room temperature. Cyan line marks zero polarization. We also examined the time evolution of the ratio of the lattice constants and the average unit-cell volume, which are shown in Fig. 21. We found the lattice to contract and the lattice ratios 𝑐/𝑎 and 𝑎/𝑏 to equilibrate at 1 ps, further indicating the formation of a cubic non-polar phase. 74 Figure 21. Lattice evolution under excitation in NPT ensemble. (A) Evolution of the lattice constant ratios. (B) Evolution of the average unit-cell volume for the 4×4×4 supercell. Contraction of the lattice is consistent with the observed activation of the A1 TO1 and E TO2 modes. To better understand the formation of octahedral rotations, we performed Mulliken bond overlap analysis of the Ti-O and Pb-O bonded atoms for the excited state trajectory, which is plotted in Fig. 22. Figure 22. Mulliken bond-overlap population analysis. (A) and (B) show the Pb-O and Ti-O bond overlap populations, respectively, both in the ground and excited states. In comparison to a ground-state QMD trajectory, weakening of both Ti-O and Pb-O bonds is seen in the excited state; however, the relative Pb-O overlap in comparison to the Ti-O overlap increases indicating an increase in relative strength of the Pb-O interaction. This is consistent with the notion 75 that tilting in ABO3 perovskites is often thought to result from a necessity of stabilizing short A- O interactions 187 . Previous experimental study of the ABO3 perovskite EuTiO3 (ETO) has shown a similar photoinduction of octahedral tilting 188 , and the observed crossover between polar and octahedral tilting modes under photo-excitation is consistent with static DFT investigations of both PTO and BaTiO3 (BTO) cubic structures under photoexciation 186 . We also examined the effects of far-from equilibrium optical excitation on charge and bond dynamics in the crystal using Mulliken analysis 165 . Fig. 23 shows time evolution of the oxygen charges divided into oxygen atoms located perpendicular to the polarization axis 𝑂 and those parallel to polarization axis 𝑂 ∥ relative to Ti atom, with the shading representing the standard deviation. Fig. 23A for the ground state shows the existence of two distinct oxygen charges depending on the oxygen position with respect to the polarization axis. In the excited state within the iso-thermal iso-baric (NPT) ensemble, the charge is initially uniformly more positive across all oxygen atoms as charge is transferred from O to Ti upon optical excitation as shown in Fig. S5B. As the structure transforms to a cubic non-polar phase the two types of oxygen atoms become identically charged. If the lattice is not allowed to relax and dynamics are performed in the iso- thermal iso-volume (NVT) ensemble, the O charge still becomes more positive, but the two types of O atoms still remain distinctively charged. The charge dynamics are more frustrated with a large standard deviation as the lattice is being optically strained. 76 Figure 23. Oxygen charge dynamics. (A) O charge dynamics in the ground state within NVT ensemble. (B) and (C) O charge dynamics in the excited state within NPT and NVT ensembles, respectively. Shading represents standard deviation. We also quantified the effects of non-adiabatic coupling between electrons and phonons using NAQMD simulation 68,133 based on surface-hopping 175,176 . Fig. 24 shows time evolutions of polarization, fractional excited electron participation (fep), and fractional hole participation (fhp) on given unit-cell within a 4×4×4 PTO supercell with 5 carriers excited (neh = 0.48%). The fep and fhp are computed as the contribution of the atoms in a given unit cell to the excited state electron and hole wave-function respectively using Mulliken analysis 165 . Scattering of carriers with phonons will drive the localization of the excited electron and hole pairs from different lattice sites. When the electrons and holes are more concentrated on the given lattice site, the unit-cell polarization will trend towards to zero, but when the carriers scatter and de-localize from the unit- cell the polarization will rebound, resulting in large fluctuating unit-cell dipoles. Electron-electron scattering modeled by RT-TDDFT at high laser fluences suggests a more uniform scattering of carriers than that described by surface hopping. This results in more uniform charge transfer from O to Ti that can be well approximated with FOQMO. As the number of excited carriers is increased in the surface-hopping based approach, more lattice sites will be occupied with excited carriers 77 and the results will better match those of FOQMD as more lattice sites will become occupied with excited state carriers. To better model both electronic and phonon scattering at the experimental pump-probe time scales, a unified approach combining Ehrenfest dynamics based on RT-TDDFT and surface-hopping dynamics may shed light in this area, which will be focus of future works 189 . Figure 24. NAQMD simulation of PTO. (A) PTO unit cell with O types labeled. Spheres with black, silver and red colors represent Pb, Ti and O atoms, respectively. (B) Polarization dynamics for given unit cell. (C) Excited electron participation dynamics on Ti atom in given unit cell. (D) Hole participation dynamics on O atoms in given unit-cell with coloring given by O atom type as labeled in (a). The O1(green) and O2(magenta) curves are overlapping. The observed phase change can be best characterized by the order parameter 𝐷 , which describes the ferroelectric (FE) ordering of the system in terms of the Ti polar displacement, and 𝜙 A*FA , which describes the tilt ordering in terms of the tilt angle. A different 𝜙 A*FA can occur about each axis and thus in general constitutes three different order parameters; however certain tilt phases are forbidden by symmetry considerations 190 . We observe a loss of FE ordering, which is 78 manifested in merging of the split peaks in the ground-state Ti-O partial radial distribution function g(r) as shown in Fig. 25A. The loss of FE ordering is countered by emergence of tilt ordering, which can be seen as a splitting of the third peak in the O-O-O bond angle distribution as shown in Fig. 25B. The change in ordering can be described by photoexcitation-driven lowering of the FE potential well barrier (Fig. 25C) concurrently with raising of the tilt potential-well height (Fig. 25D). As the coupling of lattice motions is limited by the speed of sound within the material, this sudden onset of tilt ordering by photoexcitation in large nanostructures may result in topological defect formation as different regions acquire different tilt orientations with respect to the original tilt-free lattice. 79 Figure 25. Dual order parameters of the phase transition. (A) Ti-O partial radial distribution g(r) for both excited- state and ground-state trajectories. In the ground state, polarized structure Ti-O symmetry is broken, resulting in splitting of the first g(r) peak. The symmetry is restored upon optical excitation as the polarization is lost, resulting in a single peak. (B) O-O-O bond angle distribution for ground and excited states. Third peak splits in the excited state due to symmetry breaking from the optical induction of oxygen rotations. (C) and (D) diagram the phase change in terms of the ferroelectric (FE) order parameter 𝑫 𝒛 and the tilt order parameter 𝝓 𝒕𝒊𝒍𝒕 . Optical excitation lowers the FE energy barrier and raises the tilt energy barrier . To study photo-control of large nanostructures and potential formation of topological defects, we developed NNQMD models for both ground- and excited-state dynamics based on QMD and FOQMD training data, respectively, following workflow in Fig. 26A. Atomic positions from the FOQMD data are first used to compute highly scalable symmetry functions (see Methods Chapter) that are then fed into a neural network which is trained to the total free energy of the 80 quantum system including entropic contributions. Atomic forces are then computed by taking the derivative of the network and symmetry functions with respect to the atomic coordinate. The training data was generated in the canonical (NVT) ensemble at 300 K. While in the excited state the canonical ensemble will artificially strain the PTO lattice, experimental boundary conditions such as tensile strain can replicate this condition. There was no major structural change observed differences between NVT and NPT QMD simulations. Details on the NNQMD model training and validation are discussed in the methods Chapter. We applied the validated NNQMD model to study optically-induced phase change dynamics in both a large (749Å×749Å×12.46Å) uniformly polarized PTO crystal, Figs. 26B-D, and a PTO supercell (354Å×354Å×12.46Å) engrained with circular nano-domains of opposing polarization, Figs. 26E-G. The initial condition of opposing polarization inside the circular domains seeded the creation of Bloch type skyrmions in the ground-state NNQMD model as illustrated in the inlet of Fig. 26E. These type of Bloch skyrmions have been observed inside bulk PTO section of STO/PTO nanostructures 47 . While the initial condition of opposing circular polarized regions usually would require an interfacial material such as STO, recent experimental work has shown the ability to engrain arbitrary domains in ferroelectric crystals 191 . Thus, induction of Bloch-type skyrmions in a bulk PTO crystal is not out of the range of experimental possibility and provides a non-trivial ground-state topology whose dynamics under optical excitation can be 81 compared to the that of the optically excited dynamics in a single bulk crystal. The dynamics for both systems illustrate a wavelike nucleation of the nonpolar regions as optical excitation lowers the FE energy barrier; however, the polarization is partially preserved where the wave- propagations meet in the form of polar strings as illustrated in Figs. 26D and 26G. As optical Figure 26. Excited-state topological defect dynamics. (A) XS-NNQMD workflow. (B) Initial configurations of bulk PTO crystal. Local polarization of each TiO 6 structure is color-coded by their magnitude along c-axis. Panels (C)-(D) show the photoinduced nucleation of TiO 6 tilt domains and the formation of polar stripes in the bulk PTO system at 24 picoseconds using XS-NNQMD simulation. The polar stripes are sandwiched by domains with different TiO 6 cage tilt orientations. Panel (E) shows the initial configuration of PTO crystal with circular anti- polarization domains. The circular anti-polarization domains were first introduced by displacing Ti atom position and relaxed at low temperature using ground state NNQMD simulation. Emergent Block-type Skyrmion at each domain are shown in inset of panel (E) with the polarization vectors on a-b plane prior to optical excitation. Panels (F) and (G) show the nucleation, propagation, and formation of polar stripes in the system using XS-NNQMD simulation. Unlike in the bulk PTO crystal system, the polar stripes consist of multiple segments of remnant polarization of the original PTO crystal and the circular anti-polar domains. (H) Photoinduced domains and domain boundaries color-coded by the tilt orientation Ti-O 6 cage as order parameter. (I) Atomic coordinates around the domain boundaries highlighting TiO 6 cage tilt orientations in red. For clarity, only O-atom positions and lines connecting neighboring O atoms are shown. 82 excitation destroys the original FE topology, both the uniformly polarized and Bloch skyrmion systems gain similar polar string topologies, indicating this phase change is largely independent of the original PTO FE domain structure. In addition, we observed that the polarization-erasing waves correspond to the oxygen rotations. A snapshot of the rotation dynamics for the bulk crystal under optical excitation with each unit-cell colored by the parity of the tilt orientation about the c axis is shown in Fig. 26H. The red and green colors represent the two mirrored rotation orientations with respect to original lattice (which are illustrated in the x axis of Fig 25D), and the blue color represents the non-tilted unit-cells. While these two mirrored orientations belong to the same symmetry group, they can become frustrated if they meet at a boundary that is not commensurate to the periodicity of the tilt. Fig 26H clearly illustrates that polarization stripes occur at the meeting of two different tilt domains where the octahedral rotations are frustrated creating a tilt domain wall. A zoomed view on the local atomic structure of the TiO6 cages in the striped region is provided in Fig. 26I, which further illustrates the frustration of the TiO6 cage rotation at the domain boundaries. The frustration topologically preserves the c axis polarization in these regions as the FE and tilt phonon motions responsible for controlling these order parameters are non-linearly coupled. Similar frustration and symmetry breaking of different tilt domains has been seen at heterostructure interfaces with different tilt phases 192 and domain walls in improper ferroelectrics 193 . The wavelike propagation of octahedral rotations and formation of topological defects resembles Higgs condensation of a disordered field 70 . This type of transition has previously been observed at ferroelectric critical point of multiferroic hexagonal manganites 71 . In our case, instead of using rapid temperature quenching to induce the non-adiabatic phase-change, photoexcitation non-adiabatically changes the energy landscape by drastically lowering the FE energy barrier 83 while also creating a deep tilt potential well. In turn, topological defects emerge in the form of frustration of the tilt domains and polar strings. To investigate size dependence of the dynamics of photo-induced polar stripes and TiO6 tilt domains, we have performed a series of XS-NNQMD simulations using three different system sizes that consist of single PTO crystal and are thermalized at 150K. The system dimensions and the total number of atoms are summarized in Table 7. Fig. 27A shows the nucleation of photo- induced polar domains after 2.5 ps of XS-NNQMD simulation in system 1. Atoms are color-coded by their order parameter based on the TiO6 tilt angle. At 7.5 ps, the domains that have the same order parameter coalesce and grow into larger domains, at the same time, domains with different order parameters are well separated by the polar stripes; see Fig. S18b. Similarly, Fig. S28C and Fig. 27D shows the domain nucleation and growth for the larger system 2 at 2.5 ps and 7.5ps. By increasing the system size, we observe increasing heterogeneity in the domain nucleation and increase in domain size. For the largest system 3, we see a saturation in the domain size resulting in a large number nano-domains emerging on the multi-micrometer scale as illustrated in Fig. 27, E and Fig. 27F, which again shows the domain structures at 2.5 ps and 7.5 ps, respectively. A zoom in on domain structures for system 3 illustrated in Fig. 27G shows the substructures of polar stripes can extend over 100 nm spatial extent. This spatial extent was quantified through computation of the structure factor for the non-rotated polarized unit-cells: 𝑆(𝑸)=1+𝜌 j5Q5EA ∫𝑑 M 𝑟 𝑒 *𝑸∙𝒓 𝑔 * N/+/*1 ;* N/+/*1 (𝒓) (4.2) where 𝜌 j5Q5EA is the density of Ti atoms composing the polar stripes and 𝑔 * N/+/*1 ;* N/+/*1 (𝒓) is the radial distribution function for the Ti atoms composing the polar stripes, which is illustrated in 84 Fig. 27h. Due to the small system size, system 1 is unable to accurately capture small Q region of 𝑆(𝑸) corresponding to the large ~10-100 nm nano-domains that were only seen in system 3. For system 3 in this region a shoulder curve similar to that of polymers is seen in the 𝑆(𝑸) curve due to the string like nature of polar defects. Table 7. Total Number of Atoms and System Dimensions of XS-NNQMD simulation Number of Atoms System Dimensions (Å 3 ) System 1 1,737,280 998.9 × 749.2 × 12.5 System 2 2,949,120 1997.8 × 1498.4 × 12.5 System 3 1,045,954,560 37958.7 × 27969.5 × 12.5 85 Figure 27. Size-dependence of domain boundary dynamics. (A) and (B) show the domain dynamics for system 1 at 2.5ps and 7.5 ps respectively. Unit-cells are color-coded by the order parameters based on the tilt angle of TiO 6 cage. The original PTO structure with zero tilt angle is shown blue. Red and green show positive and negative tilt angles due to the light-induced symmetry breaking described in the main text. (C) and (D) show the domain dynamics for system 2 and (E) and (F) show the dynamics system 3 again at 2.5ps and 7.5ps. (g) shows a zoom in on the domain structures of system 3. (H) shows the computed structure factor S(Q) at small Q for system 1 and system 3. System 1 is unable to capture the small Q region of S(Q) due to inadequate system size. 86 4.1.2 Discussion In conclusion, our XS-NNQMD simulation incorporating massive light-induced electronic excited state allowed us to study ultrafast polarization control of PTO-based nanodomains. Simulation results reveal how intricate coupling between acoustic tilting and TO polarization controlling phonon modes under optical excitation can lead to unique topological structures in the form of polarization stripes. Photoinduction of octahedral tilts in PTO is similar to pervious experimental study of ETO 188 and the observed crossover of tilting and polar phonon dynamics under photo-excitation is consistent with previous static DFT study of the PTO cubic structure under excitation 186 . Coupling of octahedral tilting with polar phonon modes in perovskites and manganites has also been attributed to a similar form of spontaneous symmetry breaking in the form of Higgs and Goldstone like phonon excitations 194,195 that have potential analogues to the Kibble-Zurek like phase transition observed here. Here, we would like to emphasize two aspects. First, in contrast to conventional Kibble-Zurek transitions 70,71 , where the system is driven from disordered to ordered phases, a transformation analogous to an inverse Kibble-Zurek transition 72 is seen here, where polarization is driven from an equilibrium ordered phase to a disordered phase in highly non-equilibrium state, and string-like topological defects arise due to the emergence of a hidden order parameter in the form of octahedral tilting. The uncovered light-induced octahedral rotations demonstrate the power of this simulation framework for exploring ultrafast control of topological phase changes. Finally, the multi-scale XS-NNQMD framework employed here is not unique to PTO or ferroelectrics, and can be applied to broader light-induced phase changes for which adequate training data can be generated. By incorporating time-dependent carrier densities through calculation of electron-hole recombination rates, the framework will allow for investigation of 87 more extensive topological defects on the micro-second time scale. Overall, this work has introduced a powerful framework that incorporates the forefront of multi-scale quantum simulation and machine learning for exploring nonthermal phase changes, with a specific application to PTO, where a unique topologically protected polar string phase was found. Further investigation of such topological control in PTO and other ferroelectric materials could lead to development of novel topological ferroelectric devices on ultrafast timescales. 4.1.3 Further Simulation Details QMD Simulations: QMD simulation follows the trajectories of all atoms, while computing interatomic forces from first principles in the framework of DFT calculations 4,120 . In order to incorporate electronic excitations, FOQMD simulations 66 were carried out using the highly parallelized plane-wave based software, QXMD 132,133 . A 4×4×4 PTO tetragonal supercell was simulated using the gamma- point for Brillouin-zone sampling. A plane-wave basis was used with a cutoff energy of 30 Ry (410 eV) for the wave functions and 250 Ry (3,400 eV) cutoff for the charge density. Vanderbilt style ultra-soft pseudopotentials (USPPs) were used, and local density approximation (LDA) was used for the exchange-correlation functional 130 . MD simulation was carried out in the isothermal- isobaric (NPT) ensemble using a time step of 1.206 fs. Phonon Calculation: Phonon dispersion 4,1204,120was computed using DFT in the Vienna Ab-Initio Software Package (VASP) 135,136 . Calculations were performed on a 2×2×2 PTO supercell in the tetragonal structure at the experimental lattice constants, and a 8×8×8 k-point grid using Monkhorst-Pack sampling was used to sample the Brillouin zone 196 . LDA was used for the exchange-correlation functional 197 . A plane-wave basis with an energy cutoff of 800 eV 88 within the projected augmented wave-vector (PAW) method was used to calculate electronic states 138,139 . Ionic optimization was performed until forces were less than 5×10 -4 ev/Å. The dynamical matrix was calculated using finite differences, which was then used to calculate the phonon band structure and density states using the phonopy software 198 . 4.2 Symmetry Guided Mechanical control of Polar Skyrmions with hybrid NN/MM simulations 4.2.1 Simmulation Results To explore the effects of mechanical manipulation of STO/PTO superlattices, we first developed NN/MM forcefields to describe their interactions. A diagram of the workflow for the neural network quantum molecular dynamics (NNQMD) 76,85–87 framework used for the NN engine to describe the dynamics in the PTO nano-layer within the NN/MM model is illustrated in Fig. 28a. Atomic positions for a given lattice configuration from a density functional theory (DFT) based quantum molecular dynamics (QMD) trajectory are used to compute highly scalable symmetry functions 134,199 that are then fed into a neural network, which is trained to predict the total energy of the quantum system computed from DFT. The NNQMD model was evaluated in terms of structure prediction with computation of the radial distribution functions, and elastic properties with computation of bulk modulus and elastic constant tensor. The NNQMD model was found to have good agreement with the DFT data and experiment. Details on the DFT training-data generation as well as model training and validation are provided in the methods chapter. To apply the appropriate boundary conditions of paraelectric (zero) polarization field and strain applied by the cubic STO layer in an STO/PTO heterostack, an empirical interaction was fit to reproduce the bulk modulus and cubic lattice constant for STO so 89 as to describe its applied mechanical strain on the PTO layer correctly while maintaining a zero polarized field. The MM model was also found to have good agreement with experiment with elastic modulus tensor. Details of MM force field are provided in the methods chapter. Figure 28. NN/MM model of STO/PTO super lattice ground state structure. (a) NNQMD framework used for NN engine. Atomic coordinates are used to construct symmetry functions which are fed into neural-network to predict energy, force, and stresses (b) STO-PTO heterostack of 8×16×8 STO-PTO-STO unit cells with periodic boundary conditions. R c and R B represent the core and buffer radius to linearly mix the NNQMD and MM interactions at the boundary. (c) Represents the initial condition of circular oppositely polarized domain for STO/PTO super lattice. (d) Top perspective view of STO/PTO skyrmion structure after relaxation. (e) Front perspective view of the skyrmion after relaxation. (f) (001) Plane view of polarization vectors with background density colored by interpolated [001] polar displacements of Ti atoms Dz. (g) (001) Plane view of polarization vectors with background colored by the topological charge density q. 90 After validation of the prediction of desired bulk properties by the NNQMD model for PTO and MM model for STO, we examined if a combined NN/MM framework using the individually trained models correctly describe polar skyrmion formation in the heterostack. We created a heterostack such that along the [001] direction 16 unit cells of PTO were stacked in between 8 unit cells of STO on each side with periodic boundary conditions for a total of 32 unit cells (16STO/16PTO) as illustrated in Fig. 28b. The structure was expanded with 32×32 unit cells in the [100]×[010] directions. To handle the boundary between the STO/PTO region we employ a “core+buffer” approach commonly used in linear-scaling divide and conquer DFT simualtions 132 . From the center of PTO region, a core distance RC is defined for within which the forces and energies are determined by the NNQMD model. A second buffer distance RB ( > RC) is defined for which interactions of atoms within 𝑅 z <𝑧 <𝑅 { are combined between the MM model for STO and NNQMD model for PTO with linear support functions that form a partition of unity of the entire MD box. Beyond RB, the interaction is described purely by STO MM. This approach is similar to linear interpolation approaches used in second principles based methods 34 , and other buffer domain approaches in the rapidly emerging field of NN/MM force- field development 200 . Further details on the boundary condition are provided in the methods chapter. To seed the creation of the skyrmion a 3.2 nm circle was cut out in center of the PTO region superlattice and given opposite polarization, as illustrated in Fig. 28c. Low temperature simulation (10 K) in the canonical (NVT) ensemble was then performed to obtain the ground- state structure. Upon relaxation we found the formation of skyrmions in line with what has previously reported experimentally and with second-principles and phase-field simulation approaches 34,35 . The stability to thermal noise was tested up to 300 K with simulation in the 91 NVT ensemble, for which the skyrmion was found to be stable. A top perspective view of the skyrmion texture is provided in Fig. 28c, and a front perspective view highlighting the bubble texture is shown in Fig. 28d. Fig. 28e illustrates the smooth transition from negative to positive [001] Ti polar displacements in the (001) plane, which results in formation of polarization along the [100] and [010] axes. To quantify the topological character of the skyrmion, we computed the topological charge Q and density q in the (001) plane: 𝑄 = . hf ∬𝑞(𝑥,𝑦)𝑑𝑥𝑑𝑦 = . hf ∬𝒖∙G𝜕 P 𝒖×𝜕 𝒖H𝑑𝑥𝑑𝑦, (4.3) where u is the normalized dipole moment, and x and y correspond to the [100] and [010] directions. We found 𝑄 =1, which is consistent with formation of a skyrmion. Computation was performed on an interpolated grid using the method outlined in ref. 201 to ensure an integer topological charge. After validation of the NN/MM framework to describe the formation skyrmion domains in an STO/PTO super lattice, we examined their dynamics under compression. PTO possess a low-temperature tetragonal structure in the P4mm space group that is ferroelectric and high- temperature paraelectric cubic phase belonging to the Pm3 ã m space group. In hydrostatic compression experiments and DFT calculations, multiple pressure induced phase changes have been suggested for PTO, including the aforementioned Pm3 ã m paraelectric phase, as well as rhombohedral R3C,R3 ã C phases and pseudo-cubic/tetragonal I4/mcm and I4cm phases with latter 4 space groups adopting octahedral rotations 202–206 . Here, we manipulate the P4mm symmetry of PTO ground-state structure within the skyrmion superlattice through either breaking the symmetry via uniaxial compression along the polarization axis or preserving the symmetry through isotopically compressing the system. As STO has much higher rigidity than 92 PTO in terms of its elastic constants (both in our model and experimentally, see the methods chapter), the majority of the compression will occur with the PTO structure. We first examined the effect of uniaxial compression. The heterostack was loaded at a rate of 1.0% per picosecond (ps) and then relaxed for 12 ps at 10 K. The strain was applied by uniformly Figure 29. Skyrmion evolution under uniform compression. (a) and (b) Top slice of skyrmion polarization at 0.0% and 1.0% compression along [001]. Blue double arrows are given as a guide to the eye. (c) and (d) (010) plane slice of Skyrmion structure at 0.0% and 2.7% strain. (e-f) Local c lattice constant as a function of position along the interface after compression and evolution of the compressed structure. (g)-(h) Local pressure profiles for the heterostack after compression and evolution of the compressed structure. (i) and (j) Evolution of polarization and local pressure under 3.0% compression. 93 scaling the simulation box. Upon compression we saw an increase in the ab polar displacements resulting in an increased skyrmion radius. This is illustrated in Figs. 29a and 29b which show an (001) plane top slice of the Ti-polar displacements at 0% and 1.0% compression along the c ([001]). Beyond 2% compression the ab plane polarization of the top of the skyrmion begins to decline, and we found depolarization at the interface resulting in the top plane of the skyrmion being 2 unit cells lower in the PTO layer. In the depolarized unit cells, we also saw octahedral rotations, further indicating the formation of a high local pressure as well as of a I4/mcm phase. On the contrary, the bulk part of PTO remained un-rotated. The depolarization becomes more prominent as the degree of applied strain increases. This is highlighted in Fig. 29c and Fig. 29d, which show (010) plane view of the skyrmion texture at 0% and 2.7% compression. We observe a clear “squishing” of the skyrmion due to depolarization at the interface. Upon compression of up to 3.0% we saw annihilation of the skyrmion. To understand the interfacial depolarization and skyrmion annihilation mechanisms, we plotted the average local c-axis lattice constant as function of position along the interface direction for 0% strain for reference, and 2.5%, 2.7%, and 3.0% strain just after compression as shown in Fig. 29e. The same local c lattice constants were plotted after system was allowed to evolve and reach a steady state in Fig. 29f. As the elastic constant C11 of STO is nearly 6 times larger than the elastic constant C33 of PTO (both experimentally and within our model see supporting methods chapter for details), the change in the local c lattice constants is much more compliant to STO, resulting in the majority of the compression occurring in the PTO layer. Interestingly, just after compression by 2.5%, 2.7%, and 3.0% local strain profiles are very similar with a major decline in the strain gradient especially within the bulk of the PTO layer. At 94 2.5% and 2.7% strain an appreciable strain gradient can be maintained preserving the skyrmion structure but with depolarization at interface where the gradient results in a non-polar I4/mcm phase for PTO. At 3.0% compression, the PTO region lost its tetragonality, resulting in a uniform non-polar state and the annihilation of the skyrmion structure. To further understand the depolarization mechanism, we compare the local pressure after compression and relaxation for 2.7% and 3.0% compression. In the two systems, the high- pressure regions develop at the interface, in which the 3.0% compression system shows a slight increase in the interfacial pressure. While the 2.7% compression system maintains the interfacial pressure, we observe a rapid release of the accumulated interfacial pressure and resulting complete depolarization of PTO upon 3.0% compression, indicating the interfacial pressure exceeds a mechanical threshold. Fig. 29i and Fig. 29j, provide a side-by-side snapshots of the evolution of polarization and the local pressure. Interestingly the dissipation of local pressure gradient lags the depolarization front, indicating this is an inherently dynamic process where complex interactions between polarization and strain fields play a critical role. Our simulation shows a dynamic release of large pressure accumulated at the PTO/STO interface triggering the lattice deformation in the form of octahedral rotations throughout the entirety of PTO subsequently followed by the rapid depolarization of PTO. To examine the effect of temperature, we repeated the simulations at 300 K. The results were somewhat similar but the critical strain of annihilation of skyrmion was decreased to 2.5%. At compression of up to 2.0% significant depolarization was seen at the interface and at 2.5% the skyrmion was annihilated. We next examined the dynamics of skyrmions under isotropic compression, which is expected to preserve the tetragonality of the PTO layer under pressure. Again, the strain was 95 applied by uniformly scaling the simulation box. We first compressed each axis by 2.2% over 2.2 ps (i.e., volume compression of 6.5%) at low temperature to isolate the effects of compression from thermal noise. After loading the system was then relaxed over 15 ps. A snapshot of the skyrmion structure before and 15 ps after loading is provided in Fig. 30a. Unlike the uniaxial compression that results in the total depolarization and skyrmion annihilation, the isotropic compression causes reversal of polarization in some regions creating a concentric polar pattern out of the original skyrmion, as illustrated in Fig. 30a. Fig. 30b provides snapshots of the domain formation process using Ti displacement in the [001] direction. Immediately after compression the overall [001] polarization gradually reduced before concentric polarization pattern emerged. Figure 30. Polarization dynamics under isotropic compression. (a) Snapshots of polarization structure before and after compression. (b) (001) plane perspective view of unit cells colored by [001] Ti polar displacements after compression. (c) and (d) (001) plane perspective slice of unit cells colored by in-plane ([010] and [100]) Ti polar displacements after compression. (e) Snapshots of polarization texture before and after unloading. (f) A schematic of ideal skyrmionium. (g) Skyrmionium formation visualized by Ti displacement vector on each PTO unit cells. 96 Fig. 30c and Fig. 30d demonstrate the substantial loss of the in-plane polarization under the isotropic compression as well destroying the skyrmion topology. Upon isotropic loading of the skyrmion system, we found the induction of octahedral tilting throughout the PTO layers, indicative of an I4/mcm phase. Snapshots of 4 unit cells in the center of PTO lattice before and after loading are shown on the left panel of Fig. 31a. To further quantify the structural change, we computed the O-O-O bond angle distribution with a 3 Å cutoff, which defines the octahedra, as shown on the right panel of Fig. 31a. We see a clear splitting of the ~120° peak, which further confirms the global presence of octahedral rotations. We next examined the dynamics upon unloading, where we expanded the lattice back to its original lattice constants over 2.2 ps. Upon unloading the magnitude of [001] polarization increases and the [010] and [100] polarization are restored, creating multiple domains with in- plane polarization including the newly formed concentric pattern as shown in Fig. 30e. 97 Over the concentric polarization pattern the topological charge density goes from negative to positive resulting in total topological charge of 𝑄 =0, which is known as skyrmionium that consists of a pair of skyrmions with opposite topological charge 𝑄 =1 and 𝑄 =−1 90 and reside in a concentric manner. The polarization map of the skyrmionium obtained here (Fig. 30g) was found to be remarkably similar to x-ray map of isolated magnetic skyrmioniums seen in NiFe 93 . Figure 31. Compressive skyrmionium-to-skyrmionium phase transition. (a) Compressive induction of anti- phased octahedral rotations under compression in the bulk of PTO structure. The rotation of Ti-O 6 octahedra is observed through splitting O-O-O bond angle distribution which defines their relative orientation. (c) Polarization in the [001] direction after unloading. Concentric circular domain seeds creation of skyrmions with opposite topological charges that cancel out. (c) Topological charge density, which integrates to zero over the skyrmionium surface. 98 To understand the mechanism of the skyrmionium formation, we examined the evolution of the local pressure in the system. Upon loading we see the formation high pressure gradient within PTO as illustrated in Fig. 32a. Due to the higher rigidity of STO, the gradient is slowly dissipated in the compressed structure resulting in uniform pressure individually within the PTO and STO layers, with STO being lower pressure. During dissipation we observed the evolution of pressure waves in the center of skyrmion structure as illustrated in Fig. 32b. Dissipation of these waves to a uniform pressure state led to polarization reversal, indicating up-conversion of Figure 32. Strain wave dynamics after isotropic compression. (a) Pressure evolution under isotropic loading of the skyrmion system. Isotropic loading results in high pressure state with a pressure gradient within the PTO layers. This pressure gradient is slowly dissipated after compression resulting in uniform pressure for the PTO layer, and decline in pressure for the STO layer. (b) During the dissipation process the local pressure exhibits wave like dynamics within the middle of skyrmion in (001) plane. These wave like dynamics ultimately result in reversal of the polarization in some domains. (c) Schematic of strain wave induced polarization reversal. (d) Local z-strain profiles superimposed with polarization vectors. 99 acoustic waves to optical phonon modes. The polarization response was found to correlate most strongly with z component of the stress tensor. Two snapshots of the polarization profile overlayed with the local stress are shown in Fig. 32d, where in regions of high z-stress the polarization is reduced in comparison to lower stress regions. The flipping of the polarization via these dynamic waves can be understood through the diagram in Fig. 32c. If a region in the lattice is driven to high z stressed state by the dynamic strain propagation, the potential energy for polarization flipping is lowered allowing individual unit-cells to be driven over the barrier. After the wave dissipates and the strain becomes uniform, the flipped polarization of the unit-cell will become locked as the original energy barrier is restored. To examine the robustness of the skyrmion-to-skyrmionium phase transition pathway against strain rate and thermal noise, we thermalized the skyrmion system at 200 K and 300 K and repeated the isotropic loading and unloading simulation. At 200 K, a skyrmionium transition was observed with compression of 5.3% by total volume (Fig. 33, a-c) which is similar to the one observed at 10 K (Fig. 30). The dynamics under 300 K are illustrated in Fig. 33, d-f. Under loading of 4.5% we again observed the formation of similar concentric polarization domain indicative of a skyrmionium transition, indicating the existence of the skyrmion-to-skyrmionium transition at room temperature. The effect of strain rate was examined with the system at 200 K and a lower strain rate of 100 0.1%/ps. We found that the transition threshold was reduced to 4.9%, above which the original skyrmion transformed into fragments of multiple nano-domains upon unloading. Whether the skyrmion transforms into a nontrivial polarization pattern or fragmented nanodomains will likely depend on the domain size, strain rate, and percentage. From the well-known Landau-Lifshitz- Kittel scaling laws 95–97 , which dictate the number of domain wall size should grow proportionally to the square root of the film thickness, the formation of these new nanodomains on the creation of a thinner STO/PTO film is expected. The consistently seen formation of concentric circular domains can be attributed to the topologically-protected nature of skyrmion and that the concentric circular pattern is minimal switching pathway to forming more domains within the skyrmion structure. Figure 33. Skyrmionium formation at 200 K and 300 K. (a) Snapshots of polarization texture before loading at 200K. (b) Polarization texture under compression at 200 K. (c) Polarization map after unloading forming skyrmionium structure at 200 K. (d-f) Similar Skyrmionium formation seen at 300 K. 101 4.2.2 Discussion In conclusion we have demonstrated the ability to mechanically manipulate polar skyrmions through compressive strain. We found a squishing-to-annihilation transition pathway under uniaxial compression along the c polarization axis and skyrmion-to-skyrmionium transition pathway under isotropic compression. In both strain modes, depolarization occurs when the applied strain becomes too high to maintain the skyrmion structure. Under uniaxial compression dissipation of strain gradients results in the total depolarization and formation of a uniform non- polar I4/mcm phase in the PTO layers. In contrast under isotropic compression, the propagation of pressure waves through the system activates optical phonon modes that reverse the polarization inside the skyrmion forming a skyrmionium topological structure. Rapid depolarization of the superlattice under uniaxial compression could be of interest for pulse power generation technologies through the force-electric effect 88,89 , while the increased number domain walls formed with local polarization curl in the skyrmionium is anticipated to lead to more areas of local negative permittivity and resulting potential enhancement of dielectric constant of the superlattice 36 . It should be noted that the high strain rates employed here will be difficult to achieve in a laboratory setting. Since the effects we found are dictated by symmetry, we expect them to persist at lower strain rates, though the required strain requirement will most likely be decreased. Our results highlight how mechanically manipulating the subtle symmetry of strain and polarization boundary conditions can dynamically invoke different phase changes within these next generation ferroelectric materials, and further demonstrates the ability of topological defects formed in ferromagnetic materials to be realizable in ferroelectric materials 40,41 . While in principle, a well-developed traditional classical forcefield potentially could investigate the dynamics performed here, we did not find one in the literature that could be easily adapted for this 102 study. Due to the complex orbital hybridization that gives rise to ferroelectricity 50,207 in PbTiO3 it is very difficult to develop a traditional classical force field that describes both polarization dynamics and mechanical properties well. Since we had already developed a neural-network force field for PbTiO3 76 we decided to adapt it for this study rather trying to improve exsiting classical forcefields. Combined with previous studies of light 38,75,76 and electric field 35,36,74 manipulation, the NN/MM simulations performed here provide a promising outlook for broad optical- mechanical-electrical control of polar topologies within the emerging field of ferroelectric “topotronics” 40 . 4.3 Induction and Ferroelectric Switching of Flux Closure Domains in Strained PbTiO3 with Neural Network Quantum Molecular Dynamics We first report dynamics of PbTiO3 180-degree domain-walls under electric fields utilizing our existing NNQMD forcefield 76,98 , which has been discussed in prior sections of this chapter. As discussed the, NNQMD forcefield has been validated in terms of its elastic and structural properties, and in particular has examined the dynamics of skyrmion domains (which contain 180-degree domain-walls) in PbTiO3 consistent with state-of-the-art experimental imaging techniques and phase field simulations 34,208 . To examine the switching dynamics of PbTiO3 180-degree domain walls we performed molecular dynamics simulations on a 50×6×6 supercell with an electric field along the c polarization axis at 300K. Figure 34a shows the thermalized domain wall prior to application of the electric field. The domain wall structure is 103 diffuse with reduced polarization at the boundary which is consistent Figure 34. Switching Dynamics in 180 Degree Domain Wall(DW). (a) Initial domain wall structure at 300K. (b) and (c) show propagation DW switching 1ps and 5ps after turning on the electric field respectively. (d)-(e) show progression of switching of a given layer at the DW boundary. Switching starts with in an initial flip at the DW followed quickly by a 1D flip to the strong Ti-Ti repulsion of the facing polarized domains. This followed by 2D propagation in the layer. (i)-(k) show total polarization dynamics as a function of time for three different field strengths. with previous molecular dynamics study 101,102 . Evolution of the domain wall under a field of 3MV/cm is shown in figure 34b and 34c where the switching nucleates from the domain wall boundaries illustrating a linear progression of the switching dynamics. At the domain-wall 104 boundary initially a nucleation event occurs with the flip of the polarization within adjacent oppositely polarized layer as illustrated in figure 34d. This quickly causes the adjacent neighbors along the vertical c axis to quickly flip resulting in a 1d line flip (see figure 34e). This switching then propagates along the layer until the entire layer is switched as illustrated in figure 1f and 1h. The total polarization dynamics computed via born effective charges and local polar displacements is illustrated in figure 34i-k for three different field strengths. For the lowest field strength, the nucleation limitation of the switching becomes apparent as there are long periods for which the polarization is static before a nucleation event occurs. In addition to this rod like system which allowed us to focus on the dynamics at the domain wall, we also examined ferroelectric in a more bulk like system consisting of 40x40x40 PTO unit-cells that was uniformly polarized. As illustrated in Figure 35 there is initially a long delay between activation of the E-field and the first nucleation event (~17.5 ps). After the initial nucleation a column like domain form like the 1D line flip seen in the rod like system at the domain wall boundary. From this column like domain, we see a 2D propagation outward from the domain wall as well as the formation of new adjacent nucleation cites. This type of nucleation and growth shows many qualitative features similar to previous molecular dynamics 105 studies of ferroelectric switching of PbTiO3 101,102 . Figure 35. Evolution of polarization dynamics in 40x40x40 PTO unit-cell under applied electric field of 3MV/cm. Red color is polarization of unit-cells along the field direction and blue against. During the switching process we also examined the dynamics of the local Ti polar dipole- moments near the mid-point of switching (i.e local Pz=0). Figure 36a and 36b illustrate the time evolution of individual dipole-moment during the switching process under 3MV/m for the rod like 180 DW wall system, with each frame being separated by 50 fs. The x axis is along the domain-wall direction, and the z axis is the polarization axis. The dynamics in Figure 36a closely resemble that of a Neele like rotation around the polarization axis, while in Figure 36b the dynamics resemble those of a Bloch like rotation. This indicates that both Bloch and Neele rotations play a role in switching process. We further investigated this by categorizing the 106 individual dipole-moments were characterized as Bloch, Neele, or Ising like switching based on the local Ti dipole moment right mid-point of switching (i.e local Pz=0), which is illustrated in figure 2c. While our determination of Bloch, Neele, or Ising like switching is somewhat arbitrary, the large distribution of polarization components in figure 36c illustrates that in a thermalized system polarization rotation will be prominent during the switching dynamics. Figure 36. Polarization rotation during switching. (a) Shows dynamics during switching process for a given unit- cell that resemble that of a Neele like rotation, while (b) shows more Bloch like rotation. The DW wall is along the x axis, and the polarization axis is the z axis in this figure. (c) Shows distribution of Ti polar displacements right at the point of switching (𝑫 𝒛 =𝟎Å). We next examined the dynamics in a strained PTO lattice often seen in Ferroelectric-Dielectric superlattices. We compressed PTO by 2% along the c axis and expanded by 2% along the a and b axes. The system was then thermalized to 300K over 30 ps and then an electric field was applied. Figure 37a shows an ab profile of thermalized strained structure prior to application of an electric field of 3MV/cm and Figure 37b shows the structure several ps after the application. Figure 37b illustrates a slipping of the atoms in the ab plane resulting in plastic deformation of the 107 Figure 37. Polarization dynamics in strained PTO. (a) shows ab plane lattice structure of PTO prior to application of the electric field and (b) shows plastically deformed ab plane after application of the electric field. (c) Shows the a displacement vector field of the Pb at atoms from there postions after application of the electric field. The atoms seems form Bloch like rotations as a result of the applied field. (d) Shows the polarization profile after the application of the applied field. (e) and (f) shows the local strain profile after the applied field in the a/x and b/y direction respectively. lattice after application of the applied electric field. The slipping of the atomic planes was found to be a result of a Bloch like rotation of the due to the applied field, which is illustrated in vector field plot of the displacement of the Pb atoms after application of the applied field in figure 47c. 108 This resulted in the formation of flux closure domains in the ab plane with a net polarization along the applied field direction, which is illustrated in figure 37d. The flux closure domains form due to a trade off the applied electric field and strain in the system. Through compressing along c and expanding along the ab directions we have lowered the energy barrier to polarize in the system in the ab plane. However, applying field to uniformly polarized the system along either the a or b axis will drive the system to expand along the field direction. Since the system is still constrained, this causes the atomic planes to slip and form high regions of local strain, which is illustrated in figures 37e and 37f. The slipping seen here is somewhat similar to buckling dynamics see in ferroelectric and chirality switching of vortices in strained SrTiO3(STO)/PTO superlattices 74 . After removal of the electric field, the structure was locked and the net polarization remained along the applied field direction, indicating the system was truly plastically deformed and ferroelectrically switched. In dielectric/ferroelectric superlattices local polarization curl in the form of vortices and flux closure domains form due to both charge and strain boundary conditions applied on the ferroelectric layer 34 ; however topological structures have formed in strained PTO thin films from strain and surface interaction alone 51 . Our results demonstrate here how introduction of an external field can plastically deform the lattice introducing topological defects, that are maintained by strain boundary conditions alone after removal of the field. Such dynamics we expect to find relevant in electrical manipulations of strained PTO thin films, such as polar meron lattices in strained PTO 51 . 109 Figure 38. Switching of flux closure domain. (a) Potential energy profile during the switching process(with including the field energy). (b)-(d) Polarization profiles at the beginning, middle, and end of switching. At the height of the barrier a vortex state forms which is the intermediate state. We next examined the switching dynamics of the flux closure domain by applying a field against the polarization direction. Figure 38a shows the internal potential energy of the system as a function of time during the switching process and Figure 38b-c show the polarization profile at the beginning, middle, and end of the switching process. We can see that system switches from two different flux closure domains with different net polarization. In the center of the switching pathway we find a twin vortex state with zero net polarization, and its creation and annihilation is the rate determining step the switching process. In the twin vortex state all the strain for the system is concentrated at the vortex center, creating a high energy state that system must rotate through to switch the polarization. In conclusion we have utilized a simple born effective to explore electric field induced dynamics of PbTiO3 domain wall structures. Our simple extension effectively demonstrates a 110 nucleation limited switching of PTO 180 degree domain walls consistent previous studies and illustrate that polarization rotation plays a prominent role in the switching process. With this method we then examined the dynamics of a strained PTO lattice under electric fields. We demonstrate how an external field can plastically deform the lattice introducing topological defects, that are maintained by strain boundary conditions alone after removal of the field. Understanding such dynamics will be relevant in electrical manipulations in many applications involving strained PTO thin films. In addition, the simple Born effective charge extension we find will be relevant for implementing electric field dynamics into many existing so called 2 nd generation machine learning potentials 209 . While the “3 rd and 4 th” generation machine potentials incorporate some form of charge/coulomb interaction and dynamic charge transfer 209 , they need still need to incorporate the re-screening of the electron density from the applied field, which is taken into account by definition of the Born-effective charge to first order. This will either require retraining of the effective charges or incorporation of a polarizable model (such as a core-shell model), which we have yet to see be successfully incorporated into machine-learning potentials. The extension provided here offers a simple method that will still allow for reasonable dynamics. 111 Chapter 5 : PIMD simulations for modeling of neutron scattering experiments for development of theoretical frameworks for studying hydrogen storage candidates. In this chapter I will discuss performance inelastic neutron scattering(INS) experiments at Oak Ridge National Lab (ORNL) to understand the vibrational properties of Ammonia. I report the first inelastic neutron scattering measurement of vibrational properties of ammonia along the solid to liquid phase transition with high enough resolution for direct comparisons to ab-initio simulations. Theoretical analysis reveals the essential role of nuclear quantum effects (NQEs) for correctly describing the intermolecular spectrum as well as high energy intramolecular N-H stretching modes. This is achieved by training neural-network models using ab-initio path- integral molecular dynamics (PIMD) simulations, thereby encompassing large spatiotemporal trajectories required to resolve low energy dynamics while retaining NQEs. My results not only establish the role of NQEs in ammonia but also provide general computational frameworks to study complex molecular systems with NQEs. 5.1 Simulation and Experiment Results We performed INS experiments to obtain dynamic structure factor, S(Q,E), as a function of momentum Q and energy E for liquid and solid ammonia using SEQUOIA and VISION spectrometers at Oak Ridge National Laboratories. SEQOUIA is direct geometry spectrometer that has wide range of energy transfer from sub-meV range up to an eV, with a good energy resolution, thus covering all possible vibrations in solid and liquid ammonia. Neutron scattering angles at SEQUOIA are from 3 deg to 60 deg, thus providing relatively low neutron momentum transfer at high energy transfer, which is very important to measure the intramolecular modes of ammonia, despite large mean-squared displacement of its hydrogen atoms. VISION is an indirect 112 geometry neutron spectrometer optimized for chemical spectroscopy and molecular systems. It has a wide dynamic range and high resolution, especially in the range of inter-molecular modes for ammonia. Details on the experimental setup are provided in the further detail section at the end of the chapter. Fig. 39, a and b, shows measured integrated dynamic structure factor: 𝑆(𝐸) =x 𝑑𝑄𝑆(𝑄,𝐸) (5.1) measured with SEQUOIA spectrometer for deuterated ammonia (ND3) with incident energies in the range for intermolecular vibrations of ND3. The INS spectra for protonated ammonia (NH3) measured with VISION spectrometer at base temperature 5 K and 60 K are illustrated in Fig. 39c. Fig.39, a and b, demonstrates with increasing temperature past 100 K towards the melting point (195 K), peaks in the acoustic and optical regimes display a strong anharmonic softening. At 180 K, majority of features in the spectrum are lost due to anharmonic broadening and Debye-Waller effects. Past the melting transition all prominent features in the vibrational spectrum are washed out. For ND3 phonon dispersion can be seen in the dynamic structure due to the large coherent cross section of deuterium. Fig. 39d-f and Fig. 39g-i illustrate the dynamic structure factor measured using SEQUOIA spectrometer with increasing temperature from left to right for incident energies of 30 and 70 meV respectively. At 5 K, acoustic phonon branches can be seen propagating from Q ~ 2.15 and 3.75 Å -1 , with branch energy ~ 8.0 meV, and first optical branch is seen at ~10.5 meV. As seen with the integrated structure factor the branches become increasingly smeared and softened with increasing temperature due to strong phonon 113 anharmonicity and Debye Waller effects. Figure 39. Dynamic Structure Factor and Vibrational Spectrum for Solid Ammonia. (a) and (b) integrated dynamic structure factor S(Q,E) measured with SEQUOIA spectrometer for incident energies of 30 and 70 meV at various temperatures for ND 3 (c) Integrated S(Q,E) measured with VISON spectrometer for NH 3. (d)-(f) Full S(Q,E) spectrum measured for ND 3 with incident energy 30 meV with increasing temperature from left to right. (g)-(i) Full S(Q,E) measured for ND 3 with incident energy 70 meV with increasing temperature from left to right. We also examined the vibrational energies of the intra-molecular modes as a function of temperature, which are illustrated in Figs. 40a-b for ND3 and 40d-e for NH3. The low energy modes do not change as function of temperature, despite declines in intensity due to Debye Waller effects. Interestingly though for the high-energy N-H stretching modes we see a merging of the symmetric and anti-symmetric stretching peaks, but this peak split is resolved for ND3. The merging of the symmetric and anti-symmetric stretching peaks in the protonated sample was 114 determined to be a result of NQEs which is discussed later. In addition, an increase can be seen in the high energy N-H(D) stretching peak energy in liquid phase compared to solid phase, which indicates stronger inter-molecular interactions in the solid phase. Figure 40. Intramolecular modes for ammonia. (a) and (b) illustrate molecular modes as a function of temperature measured for ND 3. (c) Illustrate a double exponential fit high energy N-D stretching modes. (d) and (e) illustrate molecular modes as a function of temperature measured for NH 3. (f) Double exponential fit for high energy N-H stretching modes, where shift in liquid phase becomes clearer. To better understand the hardening of the N-H(D) stretching peak as a result of melting, we performed exponential fits of the neutron data with the addition of linear background term which are illustrated in Figs. 40d and 40f for ND3 and NH3 respectively. For ND3 we have compared the spectra at 5 K and just before and after the melting (at 180 and 200 K), when the phonon/vibrational populations are very similar (also the multiphonon neutron scattering should be similar). For ND3 we found an increase in 𝜔 . =293→297→301 meV and 𝜔 ( = 312→ 313→316 meV when transitioning from 𝑇 =5→180→200 K for the symmetric and 115 antisymmetric peaks. For NH3 little change was seen in the peak energy in solid phase upon further populating the phonons with increasing temperature but shift in solid to liquid phase was clear with 𝜔 . =417→418→421 meV and 𝜔 ( =445→445→453 meV for 𝑇 =5→ 180→200 K. In liquid phase at 200 K, we also saw coherent excitations resembling acoustic phonons propagating from first sharp diffraction peak (Q ~ 2.1Å -1 ) in ND3 due to quasi-elastic scattering near E = 0 meV. Similar Q dependence of the quasi-elastic scattering has been observed with inelastic x-ray scattering of water 210 . 5K 180K 205K NH3 293,312 297,313 301,316 ND3 417,445 418,445 421,453 Table 8. Fitted High Energy N-H(D) stretching modes to double exponential fits. To model the observed INS spectra for solid and liquid ammonia we performed DFT simulations. Fig 41a shows the computed INS intensities for NH3 based on lattice dynamics and density functional perturbation theory (DFPT) using the Perdew-Burke-Ernzerhof (PBE) implementation 162 of the generalized gradient approximation (GGA) for exchange-correlational functional with various methods of non-local dispersion corrections (see methods for details). These are plotted in comparison with the VISION experiment at 5 K. The OCLIMAX 211,212 software was used to compute the INS intensities from DFT simulations. The simulated spectra have a strong dependence on the dispersion correction model, and while the overall profile of the spectra seems to be similar, the exact peak positions vary significantly. The variation across different models is not surprising as the lower energy motions are more sensitive to intermolecular interactions. However, none of the simulated spectra agrees with experiment well. Comparisons for computed S(Q,E) dispersions for ND3 with SEQUOIA measurement are also 116 illustrated in Fig. 41b. For these comparisons we have also investigated the use of popular hybrid functional B3LYP 213 with inclusion of van der Waals interactions via the DFT-D 167 method, and the meta-GGA SCAN functional 137 , with rvv10 214 method for van der Waals interactions. The scaling parameter in the DFT-D method for the B3LYP-DFT-D simulations was adjusted to give near the experimental lattice constant for the NH3 crystal, similar to other scaling adjustment methods 215 . The hybrid functional performs similar to the optPBE 145 van der Waals implementation, which best matched the VISION experiment among the tested models, but the SCAN+rvv10 implementation performs much worse with too hardened acoustic branches, with much lower relative intensity. The different intensity indicates much different phonon modes predicted by SCAN in comparison to experiment as the intensity is proportional to the amplitude of the eigen-vectors 211 . Figure 41. Computed vibrational spectrum for solid ammonia. (a) Computed INS spectrum for solid NH 3 with various dispersion interactions in comparison to Vision experiment. Significant discrepancies are observed, especially for the librational band (as marked by the arrow with an illustration of a librational mode in solid ammonia). (b) Computed S(Q,E) with various exchange-correlation functionals and dispersion interactions in comparison to SEQUOIA experiment. Comparison of DFT simulations of ND 3 spectrums is better than with NH 3 spectrums indicating the role of nuclear quantum effects. One might find this surprising given SCAN’s success in describing dynamics of water 100,216,217 , and SCAN has also been reported to well describe the structure of liquid 117 ammonia in comparison to elastic neutron experiements 110 . SCAN+rvv10 also under-estimates the 0K lattice constant with 𝑎 zr.& =4.890Å. For comparison we found at 0 K 𝑎 ]c{ = 5.075Å, and we adjusted the scaling parameter in the DFT-D method to give 𝑎 {M r; = 5.033 Å at 0 K. The measured value at 2 K from elastic neutron scattering is 5.048 Å 218 . Similar type scaling adjustments have been made for B3LYP with DFT-D and given lattice parameters of 𝑎 {M r; =4.987 for ammonia crystal 215 . The above values do not include finite temperature and zero-point motion. Within the quasi-harmonic approximation which takes the harmonic zero-point vibration into consideration the optPBE predicted lattice constant is 5.098 at 2K. A similar small expansion is expected for the other predicted lattice constants. For SCAN it was also found that “harder” (i.e., lower cutoff radius) GGA pseudo- potentials than typically used were needed. For softer pseudo-potentials a much worse performance and an unphysical hyper-sensitivity of the phonon modes to the density was found. This is most likely due to issues of transferability of GGA pseudo-potentials to meta-GGA functionals, which is still an open issue in the DFT community 219,220 . All the above issues were seen in both SCAN and SCAN+rvv10; however, SCAN without Vander Waals correction did report a better agreement with the lattice constant (𝑎 z =4.956Å). While optPBE and the hybrid functional perform better in comparison to SEQUOIA measurement than SCAN for ND3, they still slightly overestimate optical modes, especially in the 25-55 meV range. We also found that optPBE computed radial distribution functions in liquid state illustrate good agreement with previous neutron experiments extracted from Ref. 107, especially the N-N distribution function which arises from inter-molecular interactions. Thus, the optPBE functional adequately describes structure of the liquid state and solid state 118 well, but fails to capture the inter-molecular vibrational spectrum, indicating something still missing that goes beyond the choice of exchange functional. By visualizing the vibrational modes for NH3, it is found that the peaks at ~40 meV in simulation (Fig. 41a) are due to NH3 libration, i.e., rotation of NH3 molecules around the “umbrella” axis. This presumably corresponds to the peak at ~30 meV in the VISION experiment. Such a major difference (~30%), regardless the van der Waals model used, is highly unusual for such a seemingly simple crystal. It points two possibilities: 1) the true van der Waals interaction cannot be described by any available model; or 2) this mode is highly anharmonic thus the harmonic approximation as used by DFPT failed. Since the discrepancy is seen with the data obtained at 5 K, the anharmonicity would essentially be caused by zero-point fluctuation of H atoms (rather than thermal motions). This hypothesis would also be consistent with the fact that while there is some discrepancy with the predicted S(Q,E) spectrum for ND3, it is smaller than that observed with NH3 measurements at VISION due to heavier mass of deuterium. To evaluate the contribution of anharmonicity, finite displacement method (FDM) was used to calculate the phonons. Different from DFPT which measures the second order derivative at the very bottom of the potential energy basin, FDM can probe the vicinity by adjusting the stride size. Interestingly, Fig. 42a shows that with larger displacement, the frequency of the librational band has a significant and systematic redshift, closing the gap between experiment and simulation. This is a clear indication that the librational modes are anharmonic, and the corresponding potential energy profile is non-parabolic. However, other parts of the spectra do not show a consistent improvement, because this method is artificial and do not solve the fundamental issue. The anharmonicity of a particular mode (in this case the NH3 libration) can also be evaluated by mapping out the potential energy profile corresponding to the mode, which 119 can be obtained by nudged elastic band calculations. The result in Fig. 42b shows that the three- fold potential energy profile has a barrier of about 167 meV. As a quantum rotor in this potential well, the excitation energies can be predicted using a quantum rotor model implemented in DAVE 221 . It is found that the excitation from n = 0 to 1 has an energy of 32 meV, consistent with the libration peak position observed in experiment. Figure 42. Evaluation of nuclear quantum effects with DFT. (a) INS spectra of solid ammonia simulated with various step sizes using the finite displacement methods (e.g., FDM-01 means a step size of 0.1 Å). The drifting librational band is an indication of its anharmonicity. (b) Potential energy profile of NH 3 rotation in solid ammonia, solved by nudged elastic band method. The energy barrier is determined to be 167 meV and the associated NH 3 quantum rotor has an excitation energy (n = 0→1) of 32 meV. (c) Illustration of PIMD simulation. While typical ab initio simulation treats atoms classically while electrons quantum-mechanically, PIMD simulation uses multiple replicas of each atom to mimic nuclear quantum effect. (d) INS spectra of solid ammonia from DFPT, ab initio 120 PIMD (1-bead and 32-bead TRPMD) simulations, and VISION experiment. The PIMD simulation with 32 beads reproduces the experimentally observed peak positions very well, whereas the other two exhibit major discrepancies. The above analysis highlights the problem but does not offer a general solution. Conventional DFT lattice dynamics or molecular dynamics cannot properly describe this NQE, and the quantum rotor model cannot be used to simulate the INS spectra or easily generalized to study other modes/systems. A promising solution is path integral molecular dynamics (PIMD), in which the quantum partition function is mapped to a classical analogue by using replicas (beads) connected by springs (ring polymers) to represent each atom 122 . The background of Fig. 42c shows a typical first principles-based simulation, where the atoms are treated classically and the electron charge density is treated quantum-mechanically to compute atomic forces, which is illustrated as gray iso-surfaces. In the foreground we have highlighted one NH3 molecule from a PIMD simulation of the same atomic configuration, where each atom has 32 replicas that are harmonically coupled together. The computation of the replica simulations is embarrassingly parallel, with only fixed nearest replica communication, and the major cost is computing the energy and forces for the atoms within each replica simulation, which is done from first principles. There are different flavors of path integral implementations, but not all of them are suitable for spectroscopic modeling through time correlation functions. For example, two implementations of PIMD, ring polymer molecular dynamics 124 (RPMD) and centroid molecular dynamics 123 (CMD), have been tested for simulation of infrared spectra. It is found that high frequencies modes tend to suffer from artifacts: the resonance problem in RPMD showing spurious peak splitting and the curvature problem in CMD showing red-shift and broadening of the stretching peak. A recently proposed method by adding a thermostat to the conventional RPMD 125 (TRPMD) provides a potential solution to address both issues. In Fig. 42f, four spectra 121 are compared, the experimentally measured INS spectrum, the simulated INS spectrum using DFPT, and two spectra simulated from TRPMD, with 1 bead (equivalent to conventional molecular dynamics without NQE) and 32 beads, respectively. The DFPT spectrum and the 1- bead TRPMD spectrum are similar; however, there is significant loss in resolution in the of simulated curve due to the small sample-size and limited run-time available with ab intio MD. The simulated spectrum from 32-bead TRPMD looks very different and is much closer to the experimentally measured spectrum. Specifically, the positions of the translational band, librational band, as well as the peaks at higher energy, are all in much better agreement with experiment. This lack of resolution leaves though much to be desired when trying to conclusively determine majority of spectrum discrepancies arise from NQEs. As most of the computational expense for ab intio PIMD simulations comes from having to compute multiple replica DFT simulations, the computational cost can be greatly decreased if the underlying DFT simulations can be replaced by much cheaper computational models. To circumvent this issue, we performed NNQMD based TRPMD simulations utilizing the recently developed group-theoretically equivariant neural-network forcefield model Allegro 128 , whose architecture is illustrated in Fig. 43a. Allegro falls under the class of graph-based models where each atom i represents a node in the graph and edges are the interatomic distances 𝑟 *, , while retaining data locality of features to achieve high computational speed 126 . For each node i the model is designed to predict the edge/pair energies 𝐸 *, , and the total energy is the sum over all pair energies for each node. Forces are then computed through gradients of the total energy. Based off the work of the state of the art neural equivariant interatomic potential (NEQUIP) 127 , Allegro utilizes equivariant features in the form linear embeddings of spherical harmonic projections of atomic pair wise atomic distances (𝑟 *, ). The model is “equivariant” in 122 that only tensor operations on spherical harmonic projectors that respect symmetry of the Euclidean group E(3) are performed. Allegro also allows for scalable and faster computations in comparisons to previous works as the model is local with spherical harmonic embeddings only performed on neighbors within a cutoff distance. For further details the reader is directed to the methods chapter and to Ref. 128. Training for the model was done using configurations from solid and liquid DFT-MD simulations as well as configurations from the PIMD simulations. Details on the model training are provided in the methods. With this framework we performed large-scale TRPMD simulations in crystalline phase in a ~30Å box consisting of 864 NH3 molecules (3456 atoms) at 60K. Fig. 43b shows the measured spectrum on vision spectrometer at 60 K and those computed with TRPMD utilizing the allegro model which shows excellent agreement. The 0 K DFT calculation which over- estimates the high energy optical modes is also shown for perspective. In addition, in Fig. 43c we computed the vibrational density of states for the high energy N-H stretching using classical MD and TRPMD from Fourier transform of velocity auto-correlation function, which demonstrates the merging of the symmetric and anti-symmetric seen in the SEQUOIA experiment is due to anharmonic broadening from zero-point motion. Further we can show the extended tail seen in SEQUOIA experiment is due to multi-phonon scattering through computation of INS spectrum within the incoherent approximation within OCLIMAX software which is illustrated in Fig. 43d. In addition, we found computation for ND3 vibrational spectrum did not show this peak merging and is illustrated in Fig. 43e. This is consistent with SEQUOIA measurement and further validates the role of zero-point motion in anharmonic broadening of this mode. 123 Figure 43. TRPMD simulations of ammonia using NNQMD. (a) Allegro model for NNQMD. (b) Comparison of computed vibrational spectrum using Allegro TRPMD and VISION spectrum which shows excellent agreement. DFT results at 0 K are shown for comparison. (c) Comparison of classical and TRPMD simulation of high-energy N-H stretching mode. Merging of Symmetric and Anti-symmetric splitting agrees with SEQUOIA experiment. (d) INS computation of N-H stretching mode with and without inclusion of mutli-phonon effects. (e) TRPMD ND 3 VDOS for high-energy N-D stretch where peak splitting remains. (f) Computation of VDOS for NH 3 liquid state, in which low energy peaks are all washed out. (g) Computed hardening of NH 3 vibrational modes in liquid state consistent with SEQUOIA experiment. We also examined the dynamics in liquid phase through PIMD simulations of 864 NH3 molecules in a cubic simulation box of side 32 Å at 205 K. We computed the vibrational density of state via Fourier transform of the velocity auto correlation function for liquid NH3 which are 124 illustrated in Fig. 43f. As seen with SEQUOIA experiment the low frequency intermolecular modes are all washed out and the computed values for the molecular peaks are in good agreement with the measured values. Fig. 43g shows the computed INS spectrum including multi-phonon scattering for solid and liquid NH3, where we also see a hardening in high energy N-H stretching modes in comparison to those computed in crystalline phase, which indicates weaker inter-molecular interactions in the liquid. As this frequency shift signifies stronger inter- molecular interaction in solid phase, which one typically attributes to hydrogen bonding, we examined the charge density overlaps in the two different phases. From DFT-MD simulations. the threshold for charge density overlap in crystalline phase was found to be on contours of 0.013 e/Å 3 , whereas in liquid phase molecules could be found sharing charge density iso- surfaces at almost twice the value seen in solid phase. Stronger charge density overlap in liquid would normally stronger inter-molecular interactions in solid phase. However, lifetime of the strong-intermolecular interactions with large density overlap in liquid phase (i.e., hydrogen-bond lifetime) is extremely short in liquid ammonia, with reported hydrogen-bond lifetimes in liquid ammonia have been on the order of ~100 fs 107,110 . Thus, in liquid phase there can instantaneous periods of strong inter-molecular interaction between a few molecules, but this interaction is transient. In contrast in solid phase, there is temporally constant weak interaction of the molecules. This on average constrains them more resulting in a softer stretching vibration, which is reflected in the measured and computed stretch peaks in the spectra 5.2 Discussion We have performed INS measurements on solid and liquid ammonia and compared the measurements to DFT simulations. We find NQE induced anharmonicity that fundamentally changes the predicted spectrum with conventional DFT simulations, which we illustrate through 125 neural network-based PIMD simulations using the TRPMD implementation of PIMD. In the liquid phase, PIMD simulations can reproduce the hardening of N-H stretching modes. The hardening was determined to be due to different spatial and temporal characters of the hydrogen bonds. In solid phase, the constant and percolated hydrogen bonding network makes the N-H stretching modes softer than in liquid phase, where brief periods of strong inter-molecular interaction is followed by periods of low/non-interaction. NQE-induced anharmonicity is expected to be quite common in molecular solids, especially in systems with “flexible” groups such as -OH, -CH3, -NH2, and for neutron spectroscopy which measures low to intermediate frequency range at low temperatures. In fact, significant discrepancies between simulated peak positions (harmonic approximation) and measured ones have been observed in many related systems, including ice, amino acids, drugs etc 222–226 . Interpretation of the discrepancies, when comparing theory and experiment, can be complicated and potentially misleading. For example, one can conceivably adjust the parameters in the dispersion correction or exchange-functional so that the simulated INS spectrum matches the experimental one, but if the simulation was done without considering NQE and/or anharmonicity, the agreement can be fortuitous and does not reflect the true intermolecular interactions. With the recent explosion of machine-learning based interatomic potentials, evaluation of NQEs in large spatiotemporal simulations at near DFT accuracy is possible, as demonstrated in a recent study of water/ice where NQEs are found to have a crucial contribution to the stability of ice Ih 227 . The framework presented here is anticipated to serve as a guide for proper analysis of INS experiments in the future. The major goal of spectral analysis is not to simply make a model that can replicate the spectrum but to have a model that can simulate the spectrum based on the correct physics, and thus also have extrapolation power by capturing the 126 fundamental interactions and nuclear quantum dynamics. This is extremely relevant for systems like ammonia where new technologies and applications such hydrogen storage will need large- scale simulation frameworks such as NNQMD simulations performed here. 5.2.3 Further Simulation and Experiment Details For SEQUOIA experiment the protonated ammonia was placed in a cylindrical annular container of 15 mm outer diameter, 0.05 mm thick and 50 mm high. The thickness for deuterated ammonia was 1 mm. Each sample was cooled to base temperature of 5 K and data was collected for several different incident energies. The samples were then heated to several subsequent temperatures for which the same incident energies were used. For VISION experiment, ammonia was condensed and solidified in an aluminum sample holder at 77 K in liquid nitrogen. It was then loaded to the instrument and further cooled to the base temperature 5 K. The INS spectra were then collected at a series temperatures including 5 K and 60 K. Density Functional Theory (DFT) calculations were performed using the Vienna Ab initio Simulation Package 135,136 (VASP). Calculations used the Projector Augmented Wave (PAW) method 138,139 to describe the effects of core electrons. Perdew-Burke-Ernzerhof 162 (PBE) implementation of the generalized gradient approximation (GGA), metaGGA SCAN 137 , and hybrid B3LYP 213 exchange-correlation functional functionals were used. Energy cutoff was 800 eV for the plane-wave basis of the valence electrons. The lattice parameters and atomic coordinates from literature 218 were used as the initial structure. For the SCAN functional it was found that hard versions of the GGA PAW potentials provided in the VASP package were required, and a 1000 eV cutoff was used in the plane wave basis. The unit-cell ground state structure was calculated on a 11×11×11 Γ-centered mesh for the unit cell (16 atoms). The total 127 energy tolerance for electronic energy minimization was 10 -8 eV, and for structure optimization it was 10 -7 eV. The maximum interatomic force after relaxation was below 0.001 eV/Å. Dispersion corrections with various models 145,228–230 (e.g., optPBE-vdW 145 functional) were applied. A 3×3×3 supercell (432 atoms) on Γ-point only was used for the calculation of phonons, using two different methods, DFPT and FDM. The vibrational eigen-frequencies and modes were then calculated by solving the force constants and dynamical matrix using Phonopy 198 . The OCLIMAX 211,212 software was used to convert the DFT-calculated phonon results to the simulated INS spectra. DFT based PIMD simulations (with TRPMD) were driven by i-Pi 231 , with the interatomic forces calculated by VASP (using optPBE-vdW 145 functional) on a 2×2×2 supercell (128 atoms). The simulations were performed with 1 bead (classical limit) and 32 beads, with a timestep of 0.25 fs and 40,000 total steps (10 ps trajectory). A 64-bead model was also tested (with shorter run time) for convergence. Temperature for the simulation was 60 K, and the reason for not running at base temperature (10 K) is because the number of beads needed for convergence is much higher at 10 K, making the simulation unfeasible. Training for the Allegro neural network force-field was done using frames from DFT- MD simulations for the liquid (205 K) and solid phase (60 K), and the above PIMD simulations. These were performed on supercells containing 108 NH3 molecules (432 atoms) using the optPBE-vdW 145 functional with box lengths of 16.2Å and 15.39Å respectively. A distance cutoff of 5.0Å was used to train the model. Allegro implements equivariant features in the form of spherical harmonic projections of atomic pair wise atomic distances (𝑟 *, ). Allegro utilizes two latent spaces: one for equivariant tensor features (𝑌 F[ ), and one for scalar invariant features 𝑟 *, ,𝑍 * ,𝑍 , , with 𝑍 * being the atomic number of species i. The equivariant latent space is updated 128 through equivariant tensor products up to a maximum order l with a weighted embedding of spherical harmonics projections of atoms within the central atom’s cutoff radius. The weights are determined by multi-layer perceptron acting on the scalar latent space features. The equivariant features of the same irreducible representation are linearly mixed, and then the scalar latent space is updated through tensor contraction of scalar and new linearly mixed equivariant features which is then passed through a multi-layer perceptron. This process is repeated for N layers and then passed through a final multi-layer perceptron for the pair energy prediction. For more details on the model, we refer the reader to reference 128 and https://github.com/mir- group/allegro. Large scale NNQMD-based MD and PIMD simulations were performed in supercells containing 864 molecules with box sizes of 30.24Å and 32.44Å for solid and liquid phase at 60K and 205 K respectively. A time-step of 0.25 fs was used and 20,000 MD steps (5 ps) of thermalization was performed before production runs of 40 ps (160,000 MD steps). Initialization of liquid structures was taken from replication of already melted simulations from DFT-MD. PIMD simulations were driven by iPi using the TRPMD implementation PIMD. Both the LAMMPS(https://www.lammps.org/ ) 232,233 and RXMD (https://magics.usc.edu/rxmd/ ) 234 software were used as the molecular force engines. The Liquid Lib Library 235 was used to compute the velocity auto-correlation functions. 129 Chapter 6 : Conclusions In this thesis we explored three problems in materials physics and energy sciences : (1) Development of polymer dielectrics with high dielectric break-down fields for high field capacitor applications, (2) Optical, mechanical and electrical control of emergent polarization topologies in in next generation ferroelectrics, and (3) Modeling of neutron scattering experiments for development of theoretical frameworks for studying hydrogen storage candidates. In the first problem we illustrated how we could use non-adiabatic quantum simulations to discovery key features in polymer and nano-coat electronic structure to understand and predict dielectric breakdown properties. Such key-features can be used as proxies in machine-learning frameworks to try and develop better dielectric polymers. In the second problem we developed multi-scale neural-network quantum molecular dynamics(NNQMD) simulation methods to understand and control topological phases in ferroelectric materials, and in the last problem we applied a similar NNQMD framework with addition of path-integral MD simulations in conjunction with neutron experiment’s to understand the role of nuclear quantum effects in the vibrational properties of hydrogen storage candidate ammonia. The thesis presented here highlights how we can use quantum simulation to not only inform machine learning models but also use machine-learning to accelerate quantum simulations. 130 References 1. Rahman, A. Correlations in the Motion of Atoms in Liquid Argon. Phys. Rev. 136, A405– A411 (1964). 2. Car, R. & Parrinello, M. Unified Approach for Molecular Dynamics and Density- Functional Theory. Phys. Rev. Lett. 55, 2471–2474 (1985). 3. Payne, M. C., Teter, M. P., Allan, D. C., Arias, T. A. & Joannopoulos, J. D. Iterative minimization techniques for ab initio total-energy calculations: molecular dynamics and conjugate gradients. Rev. Mod. Phys. 64, 1045–1097 (1992). 4. Hohenberg, P. & Kohn, W. Inhomogeneous Electron Gas. Phys. Rev. 136, B864–B871 (1964). 5. Runge, E. & Gross, E. K. U. Density-Functional Theory for Time-Dependent Systems. Phys. Rev. Lett. 52, 997–1000 (1984). 6. Baojin Chu, X. Z. et al. A Dielectric Polymer with High Electric Energy Density and Fast Discharge Speed. 313, 334 (2006). 7. Sharma, V. Rational Design of All Organic Polymer Dielectrics. Nat. Commun. 5, 4845 (2014). 8. Su, H., Strachan, A. & Goddard, W. A. Density Functional Theory and Molecular Dynamics Studies of the Energetics and Kinetics of Electroactive Polymers: Pvdf and P(Vdf-Trfe). Phys. Rev. B Condens. Matter Mater. Phys. 70, 64101 (2004). 9. Jones, J. P., Llewellyn, J. P. & Lewis, T. J. The Contribution of Field-Induced Morphological Change to the Electrical Aging and Breakdown of Polyethylene. 12, (2005). 10. Ieda, M. Dielectric Breakdown Process of Polymers. IEEE Trans. Electr. Insul. 15, 206 (1980). 11. Sun, Y., Bealing, C., Boggs, S. & Ramprasad, R. 50+ Years of Intrinsic Breakdown. IEEE Electr. Insul. Mag. 29, 8 (2013). 12. Sparks, M. et al. Theory of Electron-Avalanche Breakdown in Solids. Phys. Rev. B Condens. Matter Mater. Phys. 24, 3519 (1981). 13. Sun, Y., Boggs, S. A. & Ramprasad, R. The Intrinsic Electrical Breakdown Strength of Insulators from First Principles. Appl. Phys. Lett. 101, 132906 (2012). 14. Zener, C. A Theory of the Electrical Breakdown of Solid Dielectrics. Proc. R. Soc. 131 London, Ser. A 145, 523 (1934). 15. Frohlich, A. Theory of Electrical Breakdown in Ionic Crystals. Proc. R. Soc. A 160, 230 (1937). 16. Von Hippel, A. Electric Breakdown of Solid and Liquid Insulators. J. Appl. Phys. 8, 815 (1937). 17. Kim, C., Pilania, G. & Ramprasad, R. From Organized High-Throughput Data to Phenomenological Theory Using Machine Learning: The Example of Dielectric Breakdown. Chem. Mater. 28, 1304 (2016). 18. Kim, C. & Ramprasad, R. Dielectric Breakdown Field of Strained Silicon under Hydrostatic Pressure. Appl. Phys. Lett. 111, 112904 (2017). 19. Cubero, D. & Quirke, N. Computer Simulations of Localized Small Polarons in Amorphous Polyethylene. J. Chem. Phys. 120, 7772 (2004). 20. Kumazoe, H. et al. Hot-Carrier Dynamics and Chemistry in Dielectric Polymers. J. Phys. Chem. Lett. 10, 3937 (2019). 21. Stark, J. Observation of the Separation of Spectral Lines by an Electric Field. Nature 92, 401 (1913). 22. Hattori, Y., Taniguchi, T., Watanabe, K. & Nagashio, K. Layer-by-Layer Dielectric Breakdown of Hexagonal Boron Nitride. ACS Nano 9, 916–921 (2015). 23. Azizi, A. et al. High-Performance Polymers Sandwiched with Chemical Vapor Deposited Hexagonal Boron Nitrides as Scalable High-Temperature Dielectric Materials. Adv. Mater. 29, 1701864 (2017). 24. Cui, Z., Cao, Z., Ma, R., Dobrynin, A. V & Adamson, D. H. Boron Nitride Surface Activity as Route to Composite Dielectric Films. ACS Appl. Mater. Interfaces 7, 16913– 16916 (2015). 25. Cheng, S., Zhou, Y., Hu, J., He, J. & Li, Q. Polyimide films coated by magnetron sputtered boron nitride for high-temperature capacitor dielectrics. IEEE Trans. Dielectr. Electr. Insul. 27, 498–503 (2020). 26. Linker, T. M. et al. Field-Induced Carrier Localization Transition in Dielectric Polymers. J. Phys. Chem. Lett. 11, 352–358 (2020). 27. Huan, T. D. et al. A polymer dataset for accelerated property prediction and design. Sci. Data 3, 1–10 (2016). 28. Sharma, V. et al. Rational design of all organic polymer dielectrics. Nat. Commun. 5, 1–8 132 (2014). 29. Mannodi-Kanakkithodi, A., Huan, T. D. & Ramprasad, R. Mining Materials Design Rules from Data: The Example of Polymer Dielectrics. Chem. Mater. 29, 9001–9010 (2017). 30. Everschor-Sitte, K., Masell, J., Reeve, R. M. & Kläui, M. Perspective: Magnetic skyrmions—Overview of recent progress in an active research field. J. Appl. Phys. 124, 240901 (2018). 31. Tokura, Y. & Kanazawa, N. Magnetic Skyrmion Materials. Chem. Rev. 121, 2857–2897 (2021). 32. Göbel, B., Mertig, I. & Tretiakov, O. A. Beyond skyrmions: Review and perspectives of alternative magnetic quasiparticles. Phys. Rep. 895, 1–28 (2021). 33. P., P. S. S., Masamitsu, H. & Luc, T. Magnetic Domain-Wall Racetrack Memory. 320, 190–194 (2008). 34. Das, S. et al. Observation of room-temperature polar skyrmions. Nature 568, 368–372 (2019). 35. Han, L. et al. High-density switchable skyrmion-like polar nanodomains integrated on silicon. Nature 603, 63–67 (2022). 36. Das, S. et al. Local negative permittivity and topological phase transition in polar skyrmions. Nat. Mater. 20, 194–201 (2021). 37. Hsu, S.-L. et al. Emergence of the Vortex State in Confined Ferroelectric Heterostructures. Adv. Mater. 31, 1901014 (2019). 38. Stoica, V. A. et al. Optical creation of a supercrystal with three-dimensional nanoscale periodicity. Nat. Mater. 18, 377–383 (2019). 39. Du, K. et al. Manipulating topological transformations of polar structures through real- time observation of the dynamic polarization evolution. Nat. Commun. 10, 4864 (2019). 40. Tian, G., Yang, W. D., Gao, X. S. & Liu, J.-M. Emerging phenomena from exotic ferroelectric topological states. APL Mater. 9, 20907 (2021). 41. Seidel, J. Nanoelectronics based on topological structures. Nat. Mater. 18, 188–190 (2019). 42. Stojchevska, L. et al. Ultrafast Switching to a Stable Hidden Quantum State in an Electronic Crystal. Science (80-. ). 344, 177 LP – 180 (2014). 43. Lisowski, M. et al. Ultra-fast dynamics of electron thermalization, cooling and transport effects in Ru(001). Appl. Phys. A 78, 165–176 (2004). 44. Basov, D. N., Averitt, R. D. & Hsieh, D. Towards properties on demand in quantum 133 materials. Nat. Mater. 16, 1077–1088 (2017). 45. Krishnamoorthy, A. et al. Optical Control of Non-Equilibrium Phonon Dynamics. Nano Lett. 19, 4981–4989 (2019). 46. Linker, T. et al. Optically Induced Three-Stage Picosecond Amorphization in Low- Temperature SrTiO3. J. Phys. Chem. Lett. 11, 9605–9612 (2020). 47. Das, S. et al. Observation of room-temperature polar skyrmions. Nature 568, 368–372 (2019). 48. Yadav, A. K. et al. Spatially resolved steady-state negative capacitance. Nature 565, 468– 471 (2019). 49. Balachandran, P. V., Kowalski, B., Sehirlioglu, A. & Lookman, T. Experimental search for high-temperature ferroelectric perovskites guided by two-step machine learning. Nat. Commun. 9, 1668–1676 (2018). 50. Turik, A. V & Khasabov, A. G. On the origin of ferroelectricity in PbTiO3. Ferroelectrics 237, 65–71 (2000). 51. Wang, Y. J. et al. Polar meron lattice in strained oxide ferroelectrics. Nat. Mater. (2020). doi:10.1038/s41563-020-0694-8 52. Wang, L. et al. Ferroelectrically tunable magnetic skyrmions in ultrathin oxide heterostructures. Nat. Mater. 17, 1087–1094 (2018). 53. Haeni, J. H. et al. Room-temperature ferroelectricity in strained SrTiO3. Nature 430, 758– 761 (2004). 54. Smith, M. B. et al. Crystal Structure and the Paraelectric-to-Ferroelectric Phase Transition of Nanoscale BaTiO3. J. Am. Chem. Soc. 130, 6955–6963 (2008). 55. Xue, D. et al. Accelerated search for BaTiO 3 -based piezoelectrics with vertical morphotropic phase boundary using Bayesian learning. Proc. Natl. Acad. Sci. 113, 13301– 13306 (2016). 56. Lu, H. et al. Mechanical writing of ferroelectric polarization. Science. 335, 59–61 (2012). 57. Das, S. et al. Enhanced flexoelectricity at reduced dimensions revealed by mechanically tunable quantum tunnelling. Nat. Commun. 10, 1–7 (2019). 58. Wang, H. et al. Direct Observation of Huge Flexoelectric Polarization around Crack Tips. Nano Lett. 20, 88–94 (2020). 59. Rhim, S. M., Shin, M. C., Lee, S. G. & Ye, Z. G. Handbook of Advanced Dielectric, Piezoelectric and Ferroelectric Materials – Synthesis, Characterization and Applications. (2008). 134 60. Kozina, M. et al. Terahertz-driven phonon upconversion in SrTiO3. Nat. Phys. 15, 387– 392 (2019). 61. Li, X. et al. Terahertz field–induced ferroelectricity in quantum paraelectric SrTiO3. 364, 1079–1082 (2019). 62. Mankowsky, R., von Hoegen, A., Först, M. & Cavalleri, A. Ultrafast Reversal of the Ferroelectric Polarization. Phys. Rev. Lett. 118, 197601 (2017). 63. Qi, T., Shin, Y.-H., Yeh, K.-L., Nelson, K. A. & Rappe, A. M. Collective Coherent Control: Synchronization of Polarization in Ferroelectric PbTiO3 by Shaped THz Fields. Phys. Rev. Lett. 102, 247603 (2009). 64. Yang, T., Wang, B., Hu, J.-M. & Chen, L.-Q. Domain Dynamics under Ultrafast Electric- Field Pulses. Phys. Rev. Lett. 124, 107601 (2020). 65. Ahn, Y. et al. Photoinduced Domain Pattern Transformation in Ferroelectric-Dielectric Superlattices. Phys. Rev. Lett. 119, 1–6 (2017). 66. Weinert, M. & Davenport, J. W. Fractional occupations and density-functional energies and forces. Phys. Rev. B 45, 13709–13712 (1992). 67. Yabana, K., Sugiyama, T., Shinohara, Y., Otobe, T. & Bertsch, G. F. Time-dependent density functional theory for strong electromagnetic fields in crystalline solids. Phys. Rev. B 85, 45134 (2012). 68. Shimojo, F. et al. Large Nonadiabatic Quantum Molecular Dynamics Simulations on Parallel Computers. Comput. Phys. Commun. 184, 1 (2013). 69. Behler, J. Perspective: Machine learning potentials for atomistic simulations. J. Chem. Phys. 145, 170901 (2016). 70. Kibble, T. Phase-transition dynamics in the lab and the universe. Phys. Today 60, 47–52 (2007). 71. Lin, S.-Z. et al. Topological defects as relics of emergent continuous symmetry and Higgs condensation of disorder in ferroelectrics. Nat. Phys. 10, 970–977 (2014). 72. Yukalov, V. I., Novikov, A. N. & Bagnato, V. S. Realization of inverse Kibble–Zurek scenario with trapped Bose gases. Phys. Lett. A 379, 1366–1371 (2015). 73. Yang, W. et al. Quasi-one-dimensional metallic conduction channels in exotic ferroelectric topological defects. Nat. Commun. 12, 1306 (2021). 74. Piush, B. et al. Electric field control of chirality. Sci. Adv. 8, eabj8030 (2022). 135 75. Ahn, Y. et al. Photoinduced Domain Pattern Transformation in Ferroelectric-Dielectric Superlattices. Phys. Rev. Lett. 119, 57601 (2017). 76. Linker, T. et al. Exploring far-from-equilibrium ultrafast polarization control in ferroelectric oxides with excited-state neural network quantum molecular dynamics. Sci. Adv. 8, eabk262: 1-7 (2022). 77. Tang, Y. et al. Periodic Polarization Waves in a Strained, Highly Polar Ultrathin SrTiO3. Nano Lett. 21, 6274–6281 (2021). 78. Kim, K.-E. et al. Configurable topological textures in strain graded ferroelectric nanoplates. Nat. Commun. 9, 403 (2018). 79. Shao, Y.-T. et al. Emergent chirality in a polar meron to skyrmion transition revealed by 4D-STEM. Microsc. Microanal. 27, 348–350 (2021). 80. Zubko, P., Catalan, G. & Tagantsev, A. K. Flexoelectric Effect in Solids. Annu. Rev. Mater. Res. 43, 387–421 (2013). 81. Wang, B., Gu, Y., Zhang, S. & Chen, L.-Q. Flexoelectricity in solids: Progress, challenges, and perspectives. Prog. Mater. Sci. 106, 100570 (2019). 82. Karniadakis, G. E. et al. Physics-informed machine learning. Nat. Rev. Phys. 3, 422–440 (2021). 83. Warshel, A. & Levitt, M. Theoretical studies of enzymic reactions: Dielectric, electrostatic and steric stabilization of the carbonium ion in the reaction of lysozyme. J. Mol. Biol. 103, 227–249 (1976). 84. Warshel, A. Multiscale Modeling of Biological Functions: From Enzymes to Molecular Machines (Nobel Lecture). Angew. Chemie Int. Ed. 53, 10020–10031 (2014). 85. Shimamura, K. et al. Guidelines for creating artificial neural network empirical interatomic potential from first-principles molecular dynamics data under specific conditions and its application to α-Ag2Se. J. Chem. Phys. 151, 124303 (2019). 86. Krishnamoorthy, A. et al. Dielectric Constant of Liquid Water Determined with Neural Network Quantum Molecular Dynamics. Phys. Rev. Lett. 126, 216403 (2021). 87. Rajak, P. et al. Neural Network Quantum Molecular Dynamics, Intermediate Range Order in GeSe2, and Neutron Scattering Experiments. J. Phys. Chem. Lett. 12, 6020–6028 (2021). 88. Shkuratov, S. I. et al. Ultrahigh energy density harvested from domain-engineered relaxor ferroelectric single crystals under high strain rate loading. Sci. Rep. 7, 46758 (2017). 136 89. Gao, Z. et al. Giant power output in lead-free ferroelectrics by shock-induced phase transition. Phys. Rev. Mater. 3, 35401 (2019). 90. Zhang, X. et al. Control and manipulation of a magnetic skyrmionium in nanostructures. Phys. Rev. B 94, 94420 (2016). 91. Wang, J. et al. Magnetic skyrmionium diode with a magnetic anisotropy voltage gating. Appl. Phys. Lett. 117, 202401 (2020). 92. Kolesnikov, A. G., Stebliy, M. E., Samardak, A. S. & Ognev, A. V. Skyrmionium – high velocity without the skyrmion Hall effect. Sci. Rep. 8, 16966 (2018). 93. Zhang, S., Kronast, F., van der Laan, G. & Hesjedal, T. Real-Space Observation of Skyrmionium in a Ferromagnet-Magnetic Topological Insulator Heterostructure. Nano Lett. 18, 1057–1063 (2018). 94. Yang, J. et al. Robust formation of skyrmion and skyrmionium in magnetic hemispherical shells and their dynamic switching. Phys. Rev. B 104, 134427 (2021). 95. LANDAU, L. & LIFSHITZ, E. 3 - On the theory of the dispersion of magnetic permeability in ferromagnetic bodies Reprinted from Physikalische Zeitschrift der Sowjetunion 8, Part 2, 153, 1935. in (ed. PITAEVSKI, L. P. B. T.-P. in T. P.) 51–65 (Pergamon, 1992). doi:https://doi.org/10.1016/B978-0-08-036364-6.50008-9 96. Kittel, C. Theory of the Structure of Ferromagnetic Domains in Films and Small Particles. Phys. Rev. 70, 965–971 (1946). 97. Catalan, G. et al. Fractal Dimension and Size Scaling of Domains in Thin Films of Multiferroic BiFeO3. Phys. Rev. Lett. 100, 27602 (2008). 98. Linker, T. et al. Squishing Skyrmions: Symmetry-Guided Dynamic Transformation of Polar Topologies Under Compression. J. Phys. Chem. Lett. 13, 11335–11345 (2022). 99. Gonze, X. & Lee, C. Dynamical matrices, Born effective charges, dielectric permittivity tensors, and interatomic force constants from density-functional perturbation theory. Phys. Rev. B 55, 10355–10368 (1997). 100. Krishnamoorthy, A. et al. Dielectric Constant of Liquid Water Using Neural Network Quantum Molecular Dynamics. Phys Rev Lett 126, Just Accepted (2021). 101. Shin, Y.-H., Grinberg, I., Chen, I.-W. & Rappe, A. M. Nucleation and growth mechanism of ferroelectric domain-wall motion. Nature 449, 881–884 (2007). 102. Liu, S., Grinberg, I. & Rappe, A. M. Intrinsic ferroelectric switching from first principles. Nature 534, 360–363 (2016). 137 103. Tang, Y. L. et al. Observation of a periodic array of flux-closure quadrants in strained ferroelectric PbTiO3 films. Science. 348, 547–551 (2015). 104. Anderson, P. W. On the Anomalous Line-Shapes in the Ammonia Inversion Spectrum at High Pressures. Phys. Rev. 75, 1450 (1949). 105. Binbrek, O. S. & Anderson, A. Raman spectra of molecular crystals. Ammonia and 3- deutero-ammonia. Chem. Phys. Lett. 15, 421–427 (1972). 106. Nelson, D. D., Fraser, G. T. & Klemperer, W. Does Ammonia Hydrogen Bond? Science. 238, 1670–1674 (1987). 107. Boese, A. D., Chandra, A., Martin, J. M. L. & Marx, D. From ab initio quantum chemistry to molecular dynamics: The delicate case of hydrogen bonding in ammonia. J. Chem. Phys. 119, 5965–5980 (2003). 108. Fortes, A. D., Brodholt, J. P., Wood, I. G. & Vočadlo, L. Hydrogen bonding in solid ammonia from ab initio calculations. J. Chem. Phys. 118, 5987–5994 (2003). 109. Hoja, J., Reilly, A. M. & Tkatchenko, A. First-principles modeling of molecular crystals: structures and stabilities, temperature and pressure. WIREs Comput. Mol. Sci. 7, e1294 (2017). 110. Krishnamoorthy, A. et al. Hydrogen Bonding in Liquid Ammonia. J. Phys. Chem. Lett. 13, 7051–7057 (2022). 111. Buttersack, T. et al. Photoelectron spectra of alkali metal\&\#x2013;ammonia microjets: From blue electrolyte to bronze metal. Science (80-. ). 368, 1086–1091 (2020). 112. Henshaw, G., P. Parkin, I. & A. Shaw, G. Convenient, room-temperature liquid ammonia routes to metal chalcogenides. J. Chem. Soc., Dalt. Trans. 231–236 (1997). doi:10.1039/A605665B 113. Chehade, G. & Dincer, I. Progress in green ammonia production as potential carbon-free fuel. Fuel 299, 120845 (2021). 114. Hollevoet, L., De Ras, M., Roeffaers, M., Hofkens, J. & Martens, J. A. Energy-Efficient Ammonia Production from Air and Water Using Electrocatalysts with Limited Faradaic Efficiency. ACS Energy Lett. 5, 1124–1127 (2020). 115. Sutton, A. D. et al. Regeneration of Ammonia Borane Spent Fuel by Direct Reaction with Hydrazine and Liquid Ammonia. 331, 1426–1429 (2011). 116. Olschewski, M., Knop, S., Lindner, J. & Vöhringer, P. From Single Hydrogen Bonds to Extended Hydrogen-Bond Wires: Low-Dimensional Model Systems for Vibrational Spectroscopy of Associated Liquids. Angew. Chemie Int. Ed. 52, 9634–9654 (2013). 138 117. Parker, S. F., Lennon, D. & Albers, P. W. Vibrational Spectroscopy with Neutrons: A Review of New Directions. Appl. Spectrosc. 65, 1325–1341 (2011). 118. Sathyanarayana, D. N. Vibrational spectroscopy: theory and applications. (New Age International, 2015). 119. Carpenter, J., Micklich, B. & Zanotti, J. M. Neutron scattering measurements from cryogenic ammonia: a progress report. in ACoM-6 - 6th international workshop on advanced cold moderators Proceedings 236 (2004). 120. Kohn, W. & Sham, L. J. Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. 140, A1133–A1138 (1965). 121. Granroth, G. E. et al. SEQUOIA: A Newly Operating Chopper Spectrometer at the {SNS}. J. Phys. Conf. Ser. 251, 12058 (2010). 122. Feynman, R. P. & Hibbs, A. R. Quantum Mechanics and Path Integrals. McGraw-Hill, New-York. (1965). 123. Cao, J. & Voth, G. A. The formulation of quantum statistical mechanics based on the Feynman path centroid density. IV. Algorithms for centroid molecular dynamics. J. Chem. Phys. 101, 6168–6183 (1994). 124. Habershon, S., Manolopoulos, D. E., Markland, T. E. & Miller III, T. F. Ring-polymer molecular dynamics: Quantum effects in chemical dynamics from classical trajectories in an extended phase space. Annu. Rev. Phys. Chem 64, 387–413 (2013). 125. Rossi, M., Ceriotti, M. & Manolopoulos, D. E. How to remove the spurious resonances from ring polymer molecular dynamics. J. Chem. Phys. 140, 234116 (2014). 126. Thomas, N. et al. Tensor field networks: Rotation-and translation-equivariant neural networks for 3d point clouds. arXiv Prepr. arXiv1802.08219 (2018). 127. Batzner, S. et al. E(3)-equivariant graph neural networks for data-efficient and accurate interatomic potentials. Nat. Commun. 13, 2453 (2022). 128. Musaelian, A. et al. Learning Local Equivariant Representations for Large-Scale Atomistic Dynamics. arXiv Prepr. arXiv2204.05249 (2022). 129. Troullier, N. & Martins, J. L. Efficient pseudopotentials for plane-wave calculations. Phys. Rev. B 43, 1993–2006 (1991). 130. Laasonen, K., Pasquarello, A., Car, R., Lee, C. & Vanderbilt, D. Car-Parrinello molecular dynamics with Vanderbilt ultrasoft pseudopotentials. Phys. Rev. B 47, 10142–10153 (1993). 131. Blochl, P. E. Projector Augmented-Wave Method. Phys. Rev. B Condens. Matter Mater. 139 Phys. 50, 17953 (1994). 132. Shimojo, F. et al. A divide-conquer-recombine algorithmic paradigm for large spatiotemporal quantum molecular dynamics simulations. J. Chem. Phys. 140, (2014). 133. Shimojo, F. et al. QXMD: An open-source program for nonadiabatic quantum molecular dynamics. SoftwareX 10, 100307 (2019). 134. Behler, J. & Parrinello, M. Generalized Neural-Network Representation of High- Dimensional Potential-Energy Surfaces. Phys. Rev. Lett. 98, 146401 (2007). 135. Kresse, G. & Hafner, J. Ab initio molecular dynamics for liquid metals. Phys. Rev. B 47, 558–561 (1993). 136. Kresse, G. & Furthmüller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B 54, 11169–11186 (1996). 137. Sun, J., Ruzsinszky, A. & Perdew, J. P. Strongly Constrained and Appropriately Normed Semilocal Density Functional. Phys. Rev. Lett. 115, 36402 (2015). 138. Blöchl, P. E. Projector augmented-wave method. Phys. Rev. B 50, 17953–17979 (1994). 139. Kresse, G. & Joubert, D. From ultrasoft pseudopotentials to the projector augmented-wave method. Phys. Rev. B 59, 1758–1775 (1999). 140. Helnwein, P. Some remarks on the compressed matrix representation of symmetric second-order and fourth-order tensors. Comput. Methods Appl. Mech. Eng. 190, 2753– 2770 (2001). 141. Kalinichev, A. G., Bass, J. D., Sun, B. N. & Payne, D. A. Elastic properties of tetragonal PbTiO3 single crystals by Brillouin scattering. J. Mater. Res. 12, 2623–2627 (1997). 142. Li, Z., Grimsditch, M., Foster, C. M. & Chan, S.-K. Dielectric and elastic properties of ferroelectric materials at elevated temperature. J. Phys. Chem. Solids 57, 1433–1438 (1996). 143. Zhang, Y., Sun, J., Perdew, J. P. & Wu, X. Comparative first-principles studies of prototypical ferroelectric materials by LDA, GGA, and SCAN meta-GGA. Phys. Rev. B 96, 35143 (2017). 144. Beattie, A. G. & Samara, G. A. Pressure Dependence of the Elastic Constants of SrTiO3. J. Appl. Phys. 42, 2376–2381 (1971). 145. Klimeš, J., Bowler, D. R. & Michaelides, A. Chemical accuracy for the van der Waals density functional. J. Phys. Condens. Matter 22, 22201 (2009). 146. Zhang, Z. & Yates, J. T. Band Bending in Semiconductors: Chemical and Physical 140 Consequences at Surfaces and Interfaces. Chem. Rev. 112, 5520 (2012). 147. Kwon, S. K. et al. Surface Energy and Stress Release by Layer Relaxation. Phys. Rev. B Condens. Matter Mater. Phys. 72, 235423 (2005). 148. Mou, W., Ohmura, S., Shimojo, F. & Nakano, A. Molecular control of photoexcited charge transfer and recombination at a quaterthiophene/zinc oxide interface. Appl. Phys. Lett. 100, 203306 (2012). 149. Chauvet, C. & Laurent, C. Weibull Statistics in Short-Term Dielectric Breakdown of Thin Polyethylene Films. IEEE Trans. Electr. Insul. 28, 18 (1993). 150. Chen, G., Zhao, J., Li, S. & Zhong, L. Origin of Thickness Dependent DC Electrical Breakdown in Dielectrics. Appl. Phys. Lett. 100, 222904 (2012). 151. Sun, L., Marrocchelli, D. & Yildiz, B. Edge Dislocation Slows Down Oxide Ion Diffusion in Doped Ceo2 by Segregation of Charged Defects. Nat. Commun. 6, 6294 (2015). 152. Wang, C. C. et al. Computational Strategies for Polymer Dielectrics Design. Polymer (Guildf). 55, 979 (2014). 153. Huan, T. D. et al. Advanced Polymeric Dielectrics for High Energy Density Applications. Prog. Mater. Sci. 83, 236 (2016). 154. Du, B. X., Liang, H., Jin, L. & Zhang, C. Temperature dependent surface potential decay and flashover characteristics of epoxy/SiC composites. IEEE Trans. Dielectr. Electr. Insul. 25, 631–638 (2018). 155. Shen, W.-W., Mu, H.-B., Zhang, G.-J., Deng, J. & Tu, D.-M. Identification of electron and hole trap based on isothermal surface potential decay model. J. Appl. Phys. 113, (2013). 156. Zhou, C. & Chen, G. Space Charge and AC Electrical Breakdown Strength in Polyethylene. IEEE Trans. Dielectr. Electr. Insul. 24, 559 (2017). 157. Fukushima, S. et al. Effects of Chemical Defects on Anisotropic Dielectric Response of Polyethylene. AIP Adv. 9, 45022 (2019). 158. Chen, L., Huan, T. D., Quintero, Y. C. & Ramprasad, R. Charge injection barriers at metal/polyethylene interfaces. J. Mater. Sci. 51, 506–512 (2016). 159. Chen, L., Huan, T. D., Huzayyin, A., Quintero, Y. C. & Ramprasad, R. First-principles study of aluminum-polyethylene interfaces. in 2014 IEEE Conference on Electrical Insulation and Dielectric Phenomena (CEIDP) 887–890 (2014). doi:10.1109/CEIDP.2014.6995833 160. Chen, X. et al. First-principle investigation of the charge injection barriers of 141 polyethylene and polytetrafluoroethylene oligomers. J. Appl. Phys. 126, 35101 (2019). 161. Payne, M. C., Teter, M. P., Allan, D. C., Arias, T. A. & Joannopoulos, J. D. Iterative Minimization Techniques Forab Initiototal-Energy Calculations: Molecular Dynamics and Conjugate Gradients. Rev. Mod. Phys. 64, 1045 (1992). 162. Perdew, J. P., Burke, K. & Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 77, 3865 (1996). 163. Fujihira, M. & Inokuchi, H. Photoemission from polyethylene. Chem. Phys. Lett. 17, 554– 556 (1972). 164. Hinuma, Y., Grüneis, A., Kresse, G. & Oba, F. Band alignment of semiconductors from density-functional theory and many-body perturbation theory. Phys. Rev. B 90, 155405 (2014). 165. Mulliken, R. S. Electronic Population Analysis on LCAO–MO Molecular Wave Functions. I. J. Chem. Phys. 23, 1833–1840 (1955). 166. Wang, Y. et al. Tuning Surface States of Metal/Polymer Contacts Toward Highly Insulating Polymer-Based Dielectrics. ACS Appl. Mater. Interfaces 13, 46142–46150 (2021). 167. Grimme, S. Semiempirical Gga-Type Density Functional Constructed with a Long-Range Dispersion Correction. J. Comput. Chem. 27, 1787 (2006). 168. Krukau, A. V, Vydrov, O. A., Izmaylov, A. F. & Scuseria, G. E. Influence of the exchange screening parameter on the performance of screened hybrid functionals. J. Chem. Phys. 125, 224106 (2006). 169. Liechtenstein, A. I., Anisimov, V. I. & Zaanen, J. Density-functional theory and strong interactions: Orbital ordering in Mott-Hubbard insulators. Phys. Rev. B 52, R5467–R5470 (1995). 170. Shishkin, M. & Kresse, G. Self-consistent GW calculations for semiconductors and insulators. Phys. Rev. B 75, 235102 (2007). 171. Shishkin, M. & Kresse, G. Implementation and performance of the frequency-dependent GW method within the PAW framework. Phys. Rev. B 74, 35101 (2006). 172. Shimojo, F. et al. Large nonadiabatic quantum molecular dynamics simulations on parallel computers. Comput. Phys. Commun. 184, 1–8 (2013). 173. Meng, S. & Kaxiras, E. Real-time, local basis-set implementation of time-dependent density functional theory for excited state dynamics simulations. J. Chem. Phys. 129, 54110 (2008). 142 174. Shimojo, F. A Divide-Conquer-Recombine Algorithmic Paradigm for Large Spatiotemporal Quantum Molecular Dynamics Simulations. J. Chem. Phys. 140, 18A529 (2014). 175. Tully, J. C. Molecular dynamics with electronic transitions. J. Chem. Phys. 93, 1061–1071 (1990). 176. Craig, C. F., Duncan, W. R. & Prezhdo, O. V. Trajectory Surface Hopping in the Time- Dependent Kohn-Sham Approach for Electron-Nuclear Dynamics. Phys. Rev. Lett. 95, 163001 (2005). 177. Hase, M., Fons, P., Mitrofanov, K., Kolobov, A. V & Tominaga, J. Femtosecond structural transformation of phase-change materials far from equilibrium monitored by coherent phonons. Nat. Commun. 6, 8367 (2015). 178. Mannebach, E. M. et al. Ultrafast Electronic and Structural Response of Monolayer MoS2 under Intense Photoexcitation Conditions. ACS Nano 8, 10734–10742 (2014). 179. Tiwari, S. C. et al. Photoexcitation Induced Ultrafast Nonthermal Amorphization in Sb2Te3. J. Phys. Chem. Lett. 11, 10242–10249 (2020). 180. Mahjoub, R., Nagarajan, V. & Junquera, J. Structural and electronic properties of monodomain ultrathin PbTiO3/SrTiO3/PbTiO3/SrRuO3 heterostructures: A first-principles approach. J. Appl. Phys. 128, 244102 (2020). 181. Lebedev, A. I. Band offsets in heterojunctions formed by oxides with cubic perovskite structure. Phys. Solid State 56, 1039–1047 (2014). 182. Burns, G. & Scott, B. A. Lattice Modes in Ferroelectric Perovskites: PbTiO3. Phys. Rev. B 7, 3088–3101 (1973). 183. Jiao, Y.-C. et al. First-principles study of the negative thermal expansion of PbTiO3. Comput. Mater. Sci. 124, 92–97 (2016). 184. Choudhury, N., Walter, E. J., Kolesnikov, A. I. & Loong, C.-K. Large phonon band gap in SrTiO3 and the vibrational signatures of ferroelectricity in ATiO3 perovskites: First- principles lattice dynamics and inelastic neutron scattering. Phys. Rev. B 77, 134111 (2008). 185. Tomeno, I., Ishii, Y., Tsunoda, Y. & Oka, K. Lattice dynamics of tetragonal PbTiO3. Phys. Rev. B 73, 64116 (2006). 186. Paillard, C., Torun, E., Wirtz, L., Íñiguez, J. & Bellaiche, L. Photoinduced Phase Transitions in Ferroelectrics. Phys. Rev. Lett. 123, 87601 (2019). 143 187. Xiang, H. J., Guennou, M., Íñiguez, J., Kreisel, J. & Bellaiche, L. Rules and mechanisms governing octahedral tilts in perovskites under pressure. Phys. Rev. B 96, 54102 (2017). 188. Porer, M. et al. Ultrafast transient increase of oxygen octahedral rotations in a perovskite. Phys. Rev. Res. 1, 12005 (2019). 189. Tully, J. C. Perspective: Nonadiabatic dynamics theory. J. Chem. Phys. 137, 22A301 (2012). 190. Howard, C. J. & Stokes, H. T. Group-Theoretical Analysis of Octahedral Tilting in Perovskites. Acta Crystallogr. Sect. B 54, 782–789 (1998). 191. Langenberg, E. et al. Ferroelectric Domain Walls in PbTiO3 Are Effective Regulators of Heat Flow at Room Temperature. Nano Lett. 19, 7901–7907 (2019). 192. Qian, H. et al. Easy-to-use model to reveal the nature of octahedral rotation transformations in perovskites. Ceram. Int. 46, 4477–4483 (2020). 193. Catalan, G. & Scott, J. F. Physics and Applications of Bismuth Ferrite. Adv. Mater. 21, 2463–2485 (2009). 194. Marthinsen, A. et al. Goldstone-like phonon modes in a (111)-strained perovskite. Phys. Rev. Mater. 2, 14404 (2018). 195. Juraschek, D. M., Meier, Q. N. & Narang, P. Parametric Excitation of an Optically Silent Goldstone-Like Phonon Mode. Phys. Rev. Lett. 124, 117401 (2020). 196. Monkhorst, H. J. & Pack, J. D. Special points for Brillouin-zone integrations. Phys. Rev. B 13, 5188–5192 (1976). 197. Gunnarsson, O. & Lundqvist, B. I. Exchange and correlation in atoms, molecules, and solids by the spin-density-functional formalism. Phys. Rev. B 13, 4274–4298 (1976). 198. Togo, A. & Tanaka, I. First principles phonon calculations in materials science. Scr. Mater. 108, 1–5 (2015). 199. Behler, J., Lorenz, S. & Reuter, K. Representing molecule-surface interactions with symmetry-adapted neural networks. J. Chem. Phys. 127, 14705 (2007). 200. Lier, B., Poliak, P., Marquetand, P., Westermayr, J. & Oostenbrink, C. BuRNN: Buffer Region Neural Network Approach for Polarizable-Embedding Neural Network/Molecular Mechanics Simulations. J. Phys. Chem. Lett. 13, 3812–3818 (2022). 201. Berg, B. & Lüscher, M. Definition and statistical distributions of a topological number in the lattice O(3) σ-model. Nucl. Phys. B 190, 412–424 (1981). 202. Jabarov, S. G. et al. High-pressure effect on the ferroelectric-paraelectric transition in 144 PbTiO3. Phys. Solid State 53, 2300–2304 (2011). 203. Sani, A., Hanfland, M. & Levy, D. Pressure and Temperature Dependence of the Ferroelectric–Paraelectric Phase Transition in PbTiO3. J. Solid State Chem. 167, 446–452 (2002). 204. Frantti, J., Fujioka, Y. & Nieminen, R. M. Pressure-Induced Phase Transitions in PbTiO3: A Query for the Polarization Rotation Theory. J. Phys. Chem. B 111, 4287–4290 (2007). 205. Janolin, P.-E. et al. High-Pressure Effect on ${\mathrm{PbTiO}}_{3}$: An Investigation by Raman and X-Ray Scattering up to 63 GPa. Phys. Rev. Lett. 101, 237601 (2008). 206. Frantti, J. et al. Compression mechanisms of ferroelectric PbTiO3 via high pressure neutron scattering. J. Phys. Condens. Matter 30, 435702 (2018). 207. Cohen, R. E. & Krakauer, H. Electronic structure studies of the differences in ferroelectric behavior of batio3 and PbTiO3. Ferroelectrics 136, 65–83 (1992). 208. Pereira Gonçalves, M. A., Escorihuela-Sayalero, C., Garca-Fernández, P., Junquera, J. & Íñiguez, J. Theoretical guidelines to create and tune electric skyrmion bubbles. Sci. Adv. 5, eaau7023 (2019). 209. Behler, J. Four generations of high-dimensional neural network potentials. Chem. Rev. 121, 10037–10072 (2021). 210. Iwashita, T. et al. Seeing real-space dynamics of liquid water through inelastic x-ray scattering. Sci. Adv. 3, e1603079 (2017). 211. Cheng, Y. Q., Daemen, L. L., Kolesnikov, A. I. & Ramirez-Cuesta, A. J. Simulation of Inelastic Neutron Scattering Spectra Using OCLIMAX. J. Chem. Theory Comput. 15, 1974–1982 (2019). 212. Cheng, Y. Q., Kolesnikov, A. I. & Ramirez-Cuesta, A. J. Simulation of Inelastic Neutron Scattering Spectra Directly from Molecular Dynamics Trajectories. J. Chem. Theory Comput. 16, 7702–7708 (2020). 213. Stephens, P. J., Devlin, F. J., Chabalowski, C. F. & Frisch, M. J. Ab Initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields. J. Phys. Chem. 98, 11623–11627 (1994). 214. Peng, H., Yang, Z.-H., Perdew, J. P. & Sun, J. Versatile van der Waals Density Functional Based on a Meta-Generalized Gradient Approximation. Phys. Rev. X 6, 41005 (2016). 215. Civalleri, B., Zicovich-Wilson, C. M., Valenzano, L. & Ugliengo, P. B3LYP augmented with an empirical dispersion term (B3LYP-D*) as applied to molecular crystals. CrystEngComm 10, 405–410 (2008). 145 216. Chen, M. et al. Ab initio theory and modeling of water. Proc. Natl. Acad. Sci. 114, 10846 LP – 10851 (2017). 217. Zheng, L. et al. Structural, electronic, and dynamical properties of liquid water by ab initio molecular dynamics based on SCAN functional within the canonical ensemble. J. Chem. Phys. 148, 164505 (2018). 218. Leclercq, F., Damay, P. & Foukani, M. Structure of powder deuteroammonia between 2 and 180 K revisited: A refinement of the neutron diffraction pattern taking into account molecular reorientations: analysis of the diffuse intensity. J. Chem. Phys. 102, 4400–4408 (1995). 219. Yao, Y. & Kanai, Y. Plane-wave pseudopotential implementation and performance of SCAN meta-GGA exchange-correlation functional for extended systems. J. Chem. Phys. 146, 224105 (2017). 220. Furness, J. W., Kaplan, A. D., Ning, J., Perdew, J. P. & Sun, J. Accurate and Numerically Efficient r2SCAN Meta-Generalized Gradient Approximation. J. Phys. Chem. Lett. 11, 8208–8215 (2020). 221. Azuah, R. T. et al. DAVE: a comprehensive software suite for the reduction, visualization, and analysis of low energy neutron spectroscopic data. J. Res. Natl. Inst. Stand. Technol. 114, 341 (2009). 222. Li, J. Inelastic neutron scattering studies of hydrogen bonding in ices. J. Chem. Phys. 105, 6733–6755 (1996). 223. Moberg, D. R., Straight, S. C., Knight, C. & Paesani, F. Molecular Origin of the Vibrational Structure of Ice Ih. J. Phys. Chem. Lett. 8, 2579–2583 (2017). 224. Parker, S. F. & Haris, P. I. Inelastic neutron scattering spectroscopy of amino acids. Spectroscopy 22, 815790 (2008). 225. Cimas, A. & Gaigeot, M.-P. DFT-MD and vibrational anharmonicities of a phosphorylated amino acid. Success and failure. Phys. Chem. Chem. Phys. 12, 3501–3510 (2010). 226. Cherubini, M., Monacelli, L. & Mauri, F. The microscopic origin of the anomalous isotopic properties of ice relies on the strong quantum anharmonic regime of atomic vibration. J. Chem. Phys. 155, 184502 (2021). 227. Cheng, B., Engel, E. A., Behler, J., Dellago, C. & Ceriotti, M. Ab initio thermodynamics of liquid and solid water. Proc. Natl. Acad. Sci. 116, 1110–1115 (2019). 228. Dion, M., Rydberg, H., Schröder, E., Langreth, D. C. & Lundqvist, B. I. Van der Waals 146 Density Functional for General Geometries. Phys. Rev. Lett. 92, 246401 (2004). 229. Lee, K., Murray, É. D., Kong, L., Lundqvist, B. I. & Langreth, D. C. Higher-accuracy van der Waals density functional. Phys. Rev. B 82, 81101 (2010). 230. Klime Bowler, D. R. & Michaelides, A. Van der Waals density functionals applied to solids. Phys. Rev. B 83, 195131 (2011). 231. Kapil, V. et al. i-PI 2.0: A universal force engine for advanced molecular simulations. Comput. Phys. Commun. 236, 214–223 (2019). 232. Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 117, 1–19 (1995). 233. Thompson, A. P. et al. LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales. Comput. Phys. Commun. 271, 108171 (2022). 234. Nomura, K., Kalia, R. K., Nakano, A., Rajak, P. & Vashishta, P. RXMD: A scalable reactive molecular dynamics simulator for optimized time-to-solution. SoftwareX 11, 100389 (2020). 235. Walter, N. P., Jaiswal, A., Cai, Z. & Zhang, Y. LiquidLib: A comprehensive toolbox for analyzing classical and ab initio molecular dynamics simulations of liquids and liquid-like matter with applications to neutron scattering experiments. Comput. Phys. Commun. 228, 209–218 (2018).

Abstract (if available)

Abstract First principles atomistic modeling of materials from quantum mechanics allows accurate prediction of material properties, but are seriously limited in their scope of applicability by their inherent computational complexity. In practice approximations of varying degrees are used to simulate any material of relevance. Through tackling three problems in the field of energy and materials science I will climb and connect multiple levels on the quantum hierarchal latter using multi-scale quantum dynamical simulations incorporating multiple levels of theory from real time time-dependent density functional theory (RT-TDDFT) for small unit-cells consisting of a few atoms to billion-atom molecular dynamics simulations for defects and disorder. During this climb I will demonstrate how incorporation of machine learning in conjunction with quantum simulations can help us. I will demonstrate how we can use quantum simulations to discover key features of materials electronic structure for inputs in machine-learning based material design frameworks as well as how we can use machine learning accelerate computation of quantum simulations. In particular, I will discuss the utilization of these methods for three problems: (1) Development of polymer dielectrics with high dielectric break-down fields, (2) Optical, mechanical and electrical control of polarization topologies in in next generation ferroelectrics, and (3) Modeling of neutron scattering experiments for studying hydrogen storage candidates. While seemingly quite different in nature, these problems require truly dynamical simulation methods that include multiple levels of physics. The multiscale quantum simulations presented here highlight how we can incorporate quantum mechanics and machine-learning for material discovery and design.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

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

PDF

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

PDF

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

PDF

Simulation and machine learning at exascale

PDF

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

PDF

Large-scale molecular dynamics simulations of nano-structured materials

PDF

Atomistic modeling of the mechanical properties of metallic glasses, nanoglasses, and their composites

PDF

Designing data-effective machine learning pipeline in application to physics and material science

PDF

Theoretical modeling of nanoscale systems: applications and method development

PDF

Dynamic topology reconfiguration of Boltzmann machines on quantum annealers

PDF

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

PDF

Photoexcitation and nonradiative relaxation in molecular systems: methodology, optical properties, and dynamics

PDF

Molecular dynamics simulations of nanostructures

PDF

Out-of-equilibrium dynamics of inhomogeneous quantum systems

PDF

Physics informed neural networks and electrostrictive cavitation in water

PDF

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

$Multimillion-to-billion atom molecular dynamics simulations of deformation, damage, nanoindentation, and fracture in silica glass and energetic materials$

PDF

Multimillion-to-billion atom molecular dynamics simulations of deformation, damage, nanoindentation, and fracture in silica glass and energetic materials

PDF

Dynamics of wing cracks and nanoscale damage in silica glass

PDF

Sharpness analysis of neural networks for physics simulations

PDF

Heat-initiated oxidation of aluminum nanoparticles

Asset Metadata

Creator Linker, Thomas Matthew (author)

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

School College of Letters, Arts and Sciences

Degree Doctor of Philosophy

Degree Program Physics

Degree Conferral Date 2023-05

Publication Date 10/15/2023

Defense Date 03/28/2023

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

Tag machine learning,molecular dynamics,OAI-PMH Harvest,quantum simulations

Format theses (aat)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Nakano, Aiichiro (committee chair), Vashishta, Priya (committee chair), Branicio, Paulo (committee member), Haas, Stephan (committee member), Kalia, Rajiv (committee member)

Creator Email tlinker@usc.edu,tom.m.linker@gmail.com

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

Unique identifier UC113015582

Identifier etd-LinkerThom-11638.pdf (filename)

Legacy Identifier etd-LinkerThom-11638

Document Type Dissertation

Format theses (aat)

Rights Linker, Thomas Matthew

Internet Media Type application/pdf

Type texts

Source 20230417-usctheses-batch-1023 (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. The original signature page accompanying the original submission of the work to the USC Libraries is retained by the USC Libraries and a copy of it may be obtained by authorized requesters contacting the repository e-mail address given.

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