STUDY OF DIFFUSION AND CONDUCTION IN LITHIUM GARNET OXIDES LixLa3Zrx-5Ta7-xO12 (x=5-7) WITH INTERATOMIC POTENTIALS By Jin Dai A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Materials Science and Engineering—Doctor of Philosophy 2023 ABSTRACT Lithium-ion batteries, based on the pioneering work of three Nobel Laureates, are everywhere in our lives from portable electronics, electric vehicles, to grid storage. However, they currently employ liquid electrolytes containing flammable organic solvents that could lead to a fire if batteries are overheated. Solid electrolytes, also called fast-ion conductors or superionic conductors, are alternatives with the uttermost safety. Among various solid electrolytes, lithium garnet oxides are a promising family of materials due to their high ionic conductivity and electrochemical stability. This work discusses the study of diffusion and conduction in LixLa3Zrx- 5Ta7-xO12 (x=5-7) garnet oxides using computational methods. We developed two new generations of interatomic potentials, induced dipole, and machine learning, for this composition series. We compared them with existing interatomic potentials in terms of force/virial error against density- functional theory, prediction of phase transition, self-diffusivity, and ionic conductivity, and found machine learning interatomic potentials have the best accuracy. We then applied machine learning interatomic potentials to investigate the temperature and composition dependence of diffusion and conduction in bulk materials and the influence of grain boundary structure on ionic conductivity. We believe that the atomic insight obtained from this work could be worthwhile in understanding the bottleneck of materials performance and could provide guidance on further improvements. ACKNOWLEDGEMENTS First of all, I would like to express my sincere gratitude to my advisor, Dr. Wei Lai, for providing this precious opportunity to purse this work at Michigan State University. He has set the best role model as a researcher and has been a wonderful mentor during my Ph.D. study. All those weekly meetings and discussions with him not only showed me his erudition, intelligence, and patience but also guided me through the difficulties in the process. I would never be able to complete this work without his help, and I could not imagine having a better Ph.D. advisor. Thanks again with all my heart for everything he did these years. I also want to thank other members of my thesis committee, Dr. Alexandra Zevalkink, Dr. Xianglin Ke, and Dr. Hui-Chia Yu for their guidance and advice on this thesis. I’d like to thank my former and current colleagues, Dr. Mathew Klenk, Dr. Junchao Li, Dr. Qian Chen, Yue Jiang, and Yining He. I am so lucky to join this group and it was such a pleasure working with them. They are very nice and patient and always willing to help. Matt helped me a lot when I first joined the group. Junchao helped me a lot when I was a teaching assistant for MSE381. Yue’s help in the experimental part has led to some work presented here. Many thanks to Qian, we joined this group at the same time and spent so much time together studying and working. Your help is priceless and thanks for all the memories we shared. Yining, even though we have not spent so much time working together, thank you for the help in need, and wish you all the success at MSU in the future. I’d like to thank all my friends and families for all of their support over the years. In the end, special thanks to my cute daughter, Shutong, it is such a blessing to have you in my life. You are the sunshine to me and your lovely smile helped me go through all those tough times. iii TABLE OF CONTENTS CHAPTER 1 Introduction ......................................................................................................... 1 CHAPTER 2 Background ....................................................................................................... 11 CHAPTER 3 Model comparison ............................................................................................ 36 CHAPTER 4 Study of diffusion and conduction in lithium garnet oxides LixLa3Zrx-5Ta7-xO12 (x=5-7) by machine learning interatomic potential ...................................................................... 49 CHAPTER 5 Grain boundary contributions to Li-ion transport in Li7La3Zr2O12 ................... 68 CHAPTER 6 Dissertation conclusions and future work ......................................................... 88 BIBLIOGRAPHY ......................................................................................................................... 91 iv CHAPTER 1 Introduction 1.1 The world’s energy issues and energy storage systems Energy plays an important role in human lives, and most of them are generated from the fossil fuels such as coal, natural gas, petroleum, etc. However, the ever-increasing demand for energy and the limitation of fossil fuel due to the accessibility and climate change require the shift of electricity generation from fossil fuels to renewable energy sources such as wind, solar, hydropower, etc. According to the world energy outlook 2022, the global energy share of electricity generation from solar photovoltaic (PV) and wind will increase from 10% in 2021 to 40% by 2030, and 70% by 2050[1]. Therefore, to efficiently use intermittent clean energy, there is a motivation to develop energy storage (ES) systems that store the electrical energy when the demand is low and release it when the demand is high. Energy storage (ES) technologies can be classified into five scientific categories, mechanical, thermal, electrical, chemical, and electrochemical[2, 3]. Mechanical ES includes four types, pumped-storage hydroelectricity (also known as pumped hydro), compressed air energy storage (CAES), liquid air energy storage (LAES), and flywheels. Thermal ES includes thermal chemicals, sensible thermal, and latent thermal. Chemical ES stores the energy/power in the form of gas, called power to gas or PtG, like hydrogen and synthetic natural gas (SNG). Electrochemical ES such as rechargeable batteries convert electrical energy from a power grid stored in its electrodes into chemical energy and then convert it back to electric energy by an electrochemical reduction-oxidation (redox) reaction. Electrical ES like supercapacitors store energy electrostatically without chemical reactions, utilizing high surface area electrode materials and thin dielectrics to achieve higher capacitance. Figure 1a summarized the range of discharged time at a rated power and energy capacity 1 of ES systems. It can be seen that different technologies have different characteristics and can be mapped to the most suited applications. For example, PtG has the highest energy capacity and the longer discharge time at rated capacity and therefore is suitable for storing a large amount of energy that is discharged over long periods, however, they are still in the research and demonstration stage. Pumped hydro has the highest energy capacity and a daily based energy to power shift ability and therefore is most suited to support the national grid energy storage, in addition, it is the most commercially mature technology and currently occupies over 97% of the total storage capacity in operation[4]. Supercapacitors and batteries have lower energy capacity and shorter discharge time at rated power and therefore are suited for short-term and higher-power ES systems. Even though supercapacitors and batteries are similar as they both store and release electrical energy, their applications are different. Supercapacitors can release a large amount of energy in a very short time, however, their energy density is still very low as shown in Figure 1b, and therefore are suited for a small burst of power. In comparison, batteries have high energy density and can convert energy on an hourly basis and therefore are commonly used in solar and wind energy storage to compensate for the day-night load imbalance. Figure 1 (a) Comparison of the range of the discharge time at a rated power and energy capacity of different energy storage (ES) systems[4]. (b) Ragone plot of performance ranges of various energy storage devices[5]. 2 1.2 Li-ion battery Li-ion batteries (LiBs) are the most common ES system and have been widely used in our daily life from portable electronics to electric vehicles. They outperform all the other rechargeable batteries such as nickel–metal hydride (Ni-MH), Ni-cadmium (Ni-Cd), and lead acid in terms of the gravimetric and volumetric energy densities (Figure 2)[6, 7]. In 2019, the Nobel Prize in Chemistry was awarded jointly to M. Stanley Whittingham, John B. Goodenough, and Akira Yoshino for the development of rechargeable LiBs. Figure 3 shows the three-generation batteries designed by them. Figure 2 Energy density for different rechargeable batteries[6]. In 1972, Prof. Whittingham developed the first generation of rechargeable Li metal battery using titanium disulphide (TiS2) as the positive electrode, Li metal as the negative electrode, and lithium perchlorate in dioxolane as the electrolyte with the potential of 2 V. TiS2 is a light material and has a high energy density that meets the requirement of the battery electrode. In addition, TiS2 is a layered structure that easy for ions to intercalate/deinterclate. Using Li metal as a negative electrode arising from the fact that Li is the lightest metal (0.534 g cm-3) and the most electropositive (–3.04 V vs standard hydrogen electrode), which enables the high theoretical 3 specific capacity (3860 mA h g-1) and high-output voltage and therefore high-energy density, respectively. Unfortunately, this trailblazing work soon encountered a safety issue, during each subsequent charge-discharge recycle, Li dendrite formation at the interface of Li-metal/liquid electrolyte propagates and penetrates the barrier and contacts the positive electrode, causing a short circuit and leading to an explosion. In 1980, Prof. Goodenough expanded Whittingham’s work and discovered that by replacing the metal sulphide (TiS2) in the positive electrode with a metal oxide (LiCoO2), which also has a layered structure, lightweight but a higher potential of 4V, the capacity of battery could be doubled. This discovery is the decisive step in our “mobile living” life. LiCoO2 is still the dominant material for the positive electrode nowadays. In the meantime, Prof. Yoshino made another significant breakthrough by using petroleum coke as the negative electrode. Petroleum coke is also a layered structure and enables Li-ions to intercalate/deinterclate. Since the presence of Li is ionic instead of metallic state, this type of battery is called a Li-ion battery. This rocking-chair technology solved the dendrite issue in Li-metal batteries. Specifically, during the discharging (spontaneous) process, Li-ions deintercalated from the layered graphite carbon (anode) and migrate across the electrolyte to intercalate into the layered structure of metal oxide (cathode), in the meanwhile, the electrons flow from the anode to the cathode through the external circuit to compensate for electroneutrality. The reactions are: 𝑑𝑖𝑠𝑐ℎ𝑎𝑟𝑔𝑒 Anode reaction: 𝐿𝑖𝑥 𝐶 → 𝑥𝐿𝑖 + +𝑥𝑒 − + 𝐶 𝑑𝑖𝑠𝑐ℎ𝑎𝑟𝑔𝑒 Cathode reaction: 𝐿𝑖1−𝑥 𝐶𝑜𝑂2 + 𝑥𝐿𝑖 + +𝑥𝑒 − → 𝐿𝑖𝐶𝑜𝑂2 𝑑𝑖𝑠𝑐ℎ𝑎𝑟𝑔𝑒 Overall reaction: 𝐿𝑖𝑥 𝐶 + 𝐿𝑖1−𝑥 𝐶𝑜𝑂2 → 𝐶 + 𝐿𝑖𝐶𝑜𝑂2 This process is in reverse for the charging process. LiBs apply organic liquid as electrolytes, most of which are solutions of lithium salt in mixtures of two or more solvents, e.g., 1M solutions of LiPF6 dissolved in 1:1 volume EC (ethylene 4 carbonate) and PC (propylene carbonate). These kinds of organic liquid electrolytes are flammable and have poor electrochemical stability. Therefore, to solve this safety issue, there is a need for the development of non-flammable solid-state electrolytes (SSEs) for LiBs[8-10]. Figure 3 The most powerful batteries developed by the Nobel Prize winners. All the batteries undergo a spontaneous discharging state. The Nobel Prize in Chemistry 2019 - Popular information. 1.3 Motivation for the study 1.3.1 Seeking a promising solid electrolyte To be used as a solid-state electrolyte (SSE), the material needs to possess high ionic conductivity at room temperature, negligible electronic conductivity, excellent electrochemical stability with both electrodes, especially with Li metal, environmentally benign, low cost, and 5 easiness of preparation[11]. Figure 4 compared the ionic conductivity of some well-known SSEs at different temperatures along with the organic liquid electrolytes (1M solutions of LiPF6 dissolved in 1:1 volume EC: PC), which is the goal of SSEs[12-19]. It can be seen that none of these Li-ion conductors can provide both high ionic conductivity and a wide electrochemical window. Even though the LGPS exhibits remarkably high ionic conductivity, it has a very limited electrochemical stability window from 1.7 to 2.1 V, as shown in Figure 4b[20]. Overall, Li garnet oxides are believed to be a promising candidate for solid electrolytes because of their relatively high ionic conductivity and good electrochemical stability. Figure 4 (a) Ionic conductivity as a function of the temperature of a few important solid-state electrolytes from literatures[12-19, 21]. The organic liquid electrolyte is shown for comparison. SC: single crystal. (b) the electrochemical window of a few solid electrolytes, the oxidation potential to fully delithiate the material is marked by the dashed lines. Edited from [20]. 1.3.2 Experimental versus computational studies First of all, experimental samples often contain impurities, for example, they could be Al- contaminated when alumina crucibles are used in the preparation process[22]. Secondly, lithium garnet oxide Li7La3Zr2O12 exhibits a tetragonal-cubic phase transition after being exposed to H2O/CO2 during material synthesis and characterization[23]. In addition, normally experimental 6 samples are added with an arbitrary amount of extra Li to compensate for the Li loss in the preparation process, i.e. usually inaccurate Li composition. On the other hand, computational studies have advantages over experiments as simulation materials are pure phases and without the complication of structure defects like grain boundaries. 1.3.3 Accuracy versus efficiency of different computational methods A lot of computational methods have been used to study SSE. These methods generally fall into three categories, density-functional theory (DFT), semi-empirical (SE) and density- functional tight-binding (DFTB) potential, and empirical interatomic potential (IP). The relationship between these methods is shown in Figure 5[24]. Density-functional theory (DFT) is one of the most widely used computational tools for studying the properties of solid-state materials. It generates the energy, force, and stress of atoms and cells based on fundamental principles of quantum mechanics without any empirical parameters and therefore is the most accurate computational method. However, as seen in Figure 5, DFT is very slow at a time scale of ps and therefore is restricted to model systems containing typically hundreds of atoms. For example, running a molecular dynamics (MD) trajectory of 1-ps (e.g., 1000 calculations for a timestep of 1 fs) for a unit lithium garnet oxide cell with about 200 atoms takes about 15 hours, which means running a 1-ns (i.e., a common trajectory length using the IP method) would take 625 days[24]. Therefore, it is unrealistic to run DFT-based MD for big structures and at low temperatures. For example, lithium garnet oxide structures with grain boundaries (GBs), which will be discussed in Chapter 5, could contain more than ten thousand atoms. Since it usually needs to run much longer trajectories to capture the property information at low temperatures due to the slow mobility of atoms, DFT has been widely used in studying the self-diffusivity and ionic conductivity at high temperatures and then extrapolating the room temperature values according 7 to the Arrhenius relationship. However, this methodology could be inaccurate. For example, Miara et al.[25] performed a DFT-MD simulation for Li6.75La3Zr1.75Ta0.25O12 at the temperature range of 600-1500 K and extrapolated the room temperature ionic conductivity of 1.2×10-2 S/cm, which is one order magnitude higher than the hot-pressing sample (0.87×10-3 S/cm)[26]. The second category includes the SE and DFTB. DFTB has several versions, e.g., self- consistent charge (SCC) DFTB (also referred to as DFTB2), which is an approximation method of the Kohn-Sham DFT method. It has less accuracy but 2 orders of magnitude faster compared to the standard DFT. xTB[27] and PM6[28], the variants of DFTB and SE, respectively, have similar or even worse force/virial accuracy than the Core model[24]. Besides, the computational cost of this group is still very high[24]. The third group is IP models, including physically motivated models such as Core, core- shell (CS), induced dipole (ID) and fluctuating charge, and pure mathematic machine-learning interatomic potential (MLIP) models. Compared to DFT and DFTB, IP methods are faster but less accurate. In recent years, a lot of efforts have been devoted to developing accurate IP models. In particular, MLIP has become an increasingly popular tool in the study of materials science. 8 Figure 5 Schematic of different computational methods based on their length and time scales[24]. 1.4 Thesis overview With the introduction and motivation discussed above, this thesis is organized into the following five chapters. Chapter 2 provides the background of several inorganic solid electrolytes, the structure of lithium garnet oxides, computational methods applied in this work, and methods for property calculations of lithium garnet oxides. In Chapter 3, we compared MD simulation results of Li7La3Zr2O12 from three IP models, i.e., Core, CS, and ID, including the lattice parameters, self- diffusivity, jump length, residence time, ionic conductivity, dielectric constant, and elastic modulus. In Chapter 4, we constructed new potentials of DFTB and MLIP models and compared them with the IP models through force error, lattice parameter, and radial distribution function by 9 performing MD simulations of Li6La3TaZrO12. It turns out that the MLIP model approaches the DFT accuracy and matches the ionic conductivity of literature-reported experimental values. Then we used the MLIP model to study the transport of lithium garnet oxide series LixLa3Zrx-5Ta7-xO12 (x=5-7) and solve the problem of what composition gives the highest ionic conductivity. In Chapter 5, we used the MLIP model to study the GB contribution to the ionic conductivity of lithium garnet oxide Li7La3Zr2O12. In Chapter 6, we summarized the present work and discussed the future direction. 10 CHAPTER 2 Background 2.1 Inorganic solid-state electrolytes During the last few decades, a wide variety of inorganic materials have been explored and studied to be used as possible solid-state electrolytes (SSEs) in all-solid-state Li-ion batteries (ASSLiB)[11, 29-33]. Among them, Li-β-alumina, Li3N, lithium superionic conductor (LISICON), amorphous lithium phosphorus oxynitride (LiPON), perovskite-type oxides, Li garnet-type oxides, and Li10GeP2S12 (LGPS) have attracted a great deal of attention. Some of the crystal structures are shown in Figure 6. Figure 6 Historical developments of Li-ion solid electrolytes. 2.1.1 Li-β-alumina In 1967, Yao and Kummer first discovered β-alumina, which has the empirical formula Na2O•11Al2O3. It is a hexagonal structure with the lattice parameter of a=5.58 and c=22.45 Å. The 11 self-diffusivity of Na+ at 25 °C and 300 ℃ are 4.0×10-7 and 1.0×10-5 cm2/s, respectively[34]. They also synthesized Li-β-alumina by exchanging cation Na+ with Li+, which shows an expansion in the lattice (a=5.593 and c=22.642 Å) despite the smaller size of Li+ with the ionic conductivity of one to two orders of magnitude lower[34]. However, it is difficult to prepare a pure phase. 2.1.2 Li3N Li3N crystal is a hexagonal structure with space group P6/mmm (a=3.648 and c=3.875 Å)[35]. It consists of Li2N layers perpendicular to the c axis and connected by Li atoms occupying the sites between the respective N atoms. Li3N attracted high interest owing to its high Li content and its layered structure building up of sequential Li2N and Li layers. The ionic conductivity of a single crystal at 300 K is 1.2×10-3 S/cm perpendicular to the c-axis and 1.0×10-5 S/cm along the c-axis, with an activation energy of 0.29 and 0.49 eV, respectively[18]. However, its low thermodynamic decomposition potential (0.44 V) and sensitivity to moisture limit its use as a solid electrolyte[36]. 2.1.3 LISICON and thio-LISICON LISICON was first proposed by Hong in 1978 for a new Li-ion superionic conductor Li14Zn(GeO4)4 [14], which is a member of Li16-2xZnx(GeO4)4 (0.5 ≤ x≤ 3.5) series. It is an orthorhombic structure with space group Pnma (a=10.828, b=6.251, c=5.140 Å). The structure has a rigid 3D network of Li11Zn(GeO4)4, in which Li, Ge, and Zn are located in tetrahedral sites, and the remaining Li+ occupy octahedra sites. The Li+ diffusion pathway is formed through the tetrahedral and interstitial octahedral sites. It has an ionic conductivity of 0.13 S/cm at 300 ℃. Then, in 2000, Kanno et al. synthesized a family of thio-LISICON by replacing O2- in LISICON 12 with S2-[37], among which Li4+x+δ(Ge1−δ′−xGax)S4 (x=0.25) exhibits the highest ionic conductivity that reaches up to 6.5×10-5 S/cm at ambient temperature. 2.1.4 LiPON The amorphous lithium phosphorus oxynitride (LiPON) thin film was first studied by Bates et al. in 1992 and successfully incorporated into rechargeable thin-film lithium batteries[38]. LiPON was fabricated by sputtering Li3PO4 in N2. Its composition can be represented by xLi2O• yP2O5•zPON, where PON is phosphorous oxynitride. LiPON has a conductivity of ~2.3×10-6 S/cm at 25 ℃[16]. 2.1.5 Perovskite-type In 1993, Inaguma et al. first discovered cubic perovskite Li0.34(1)La0.51(1)TiO2.94(2) (a=3.871 Å) with the bulk and total ionic conductivity of 1×10-3 and 2×10-5 S/cm, respectively at room temperature[15]. The ideal perovskite has a formula of ABO3, where A and B refer to twelve- coordinated and six-coordinated cations, respectively. It is a cubic structure with space group Pm- 3m. The octahedral BO6 occupies the corner of the cube, leaving the center for cation A and vacancies. The introduction of Li modifies the structure and vacancies and the ionic conductivity changes with Li concentration. However, the major disadvantage of perovskite is the facile reduction of Ti4+ into Ti3+ via a chemical reaction with low-potential anodes such as Li metal. 2.1.6 Garnet-type In 2003, Thangadurai et al. first reported Li garnet oxides Li5La3M2O12 (M=Nb, Ta), both of which exhibit the same order of magnitude of bulk ionic conductivity of ~10-6 S/cm at 25 ℃[39]. Since then, lithium garnet oxides have attracted extensive interest, and tremendous work has been devoted to studying the structure and properties of lithium garnet oxides. The details of the 13 structure will be introduced in section 2.2. The properties of the lithium garnet family LixLa3Zrx- 5Ta7-xO12 (x=5-7) will be discussed in the following chapters using theoretical calculations. 2.1.7 LGPS In 2011, Kamaya et al. synthesized a new lithium superionic conductor, Li10GeP2S12 with a structure different from thio-LISICON[12]. It is a tetragonal phase with a lattice parameter of a=8.71771(5) and c=12.63452(10) Å. The Li ions move in one-dimensional z-axis and exhibit ionic conductivity of over 10-2 S/cm at 27 ℃ (as shown in Figure 4a), which is comparable to those of the commercial organic liquid electrolytes being used in LiBs. However, it is highly sensitive to moisture and is unstable against Li metal anodes. 2.2 Lithium garnet oxides Li garnet oxides are based on garnets with the general structural formula, A3B2(SiO4)3 or Si3A3B2O12 (space group #230, Ia-3d), where A and B refer to eight-coordinated and six- coordinated cation sites, respectively (as shown in AO8 dodecahedra and BO6 octahedra in Figure 7a). Li garnet oxides replace Si in the garnets by Li. A large number of compounds have been found with the Li garnet structure, in which element A could be replaced by La, Nd, Ba, Sr, etc., and element B could be Ta, Zr, Nb, Te, etc.[24, 40]. In LixA3B2O12, there are two kinds of polyhedra to accommodate Li atoms, i.e., LiO6 octahedron (Oh) and LiO4 tetrahedron (Td). Each LiO4 Td is surrounded by four LiO6 Oh and each LiO6 Oh is linked by two LiO4 Td. The Td and Oh are connected by a triangular oxygen face called the bottleneck for Li ion transport. There are 24 Td and 48 Oh in each cell and provides 72 sites for Li atoms in total. However, due to the exclusion principle, a cluster of Td-Oh-Td with three Li neighboring atoms does not exist, therefore the maximum concentration of Li per unit formula is 7.5 instead of 9 as there are 8 formulas in each unit cell. 14 Li garnet oxides have two types of phases, cubic phase (space group #230, Ia-30) and tetragonal phase (space group #142, I41/acd). Specifically, there are two special model materials for the cubic phase, endmember Li3Nd3Te2O12 (Li3-phase) and intermediate member Li5La3Ta2O12 (Li5-phase). In the Li3-phase, Li atoms locate exclusively and fully in Td-24d sites, as shown in Figure 7e, in which the blue rectangle represents LiO6 Oh, the pink squares represent LiO4 Td, and the orange dots represent Li atoms. The absence of Li atoms in the neighboring Oh sites reduces the repulsion between Li atoms and thus leads to very low ionic conductivity of ~10-5 S/cm at 600 °C with high activation energy (> 1 eV)[41, 42]. As Li content increase, extra Li are incorporated into the Oh sites, and due to the electrostatic repulsion, not all Td/Oh sites can be occupied at the same time, resulting in Li atoms occupying both Td and Oh sites randomly (as shown in Figure 7f). For example, in the Li5-phase, the occupancies of Td and Oh are 0.802(4) and 0.43(2), respectively[43]. Even though the average distance between Td-24d and Oh-48g sites is about 2 Å, the majority of the Li atoms are found in off-center Oh-96h sites due to strong repulsion between neighboring atoms, suggesting that Li atoms jump ~2.5 Å to the next site[44], as shown in Figure 7b. For the tetragonal phase, a typical model is the practical upper endmember Li7La3Zr2O12 (Li7-phase) at room temperature (as shown in Figure 7c, d, g), which contains 7 Li per unit formula by replacing pentavalent Ta in Li5-phase with tetravalent Zr. In the tetragonal phase, Td-24 sites are separated into Td-8a (dark pink) and Td-16e (light pink) sites and Oh-48g sites are separated into Oh-16f and Oh-32g sites. The tetragonal phase is an ordered structure, in which Td-8a, Oh-16f, and Oh-32g are fully occupied by Li atoms while Td-16e sites are fully empty. The complete ordering of Li-vacancies in a tetragonal structure leads to low ionic conductivity, which is about 2 orders of magnetic lower than that of the cubic phase[45]. However, when Li7La3Zr2O12 is exposed to the air, the Li-ordering can be disrupted by protons due to Li+/H+ 15 exchange and induce the tetragonal-to-cubic phase transition[23]. Furthermore, Li7La3Zr2O12 undergoes an intrinsic phase transition when the temperature increase. It has been reported that the tetragonal-to-cubic phase transition occurs between 373 to 923 K depending on the purity of the sample[24]. Figure 7 (a) 3D schematics of Li garnet oxide (LixA3B2O12) structures. (b) and (c) 2D schematic of a polyhedral ring of cubic and tetragonal phase, respectively. (d) Characteristic of 8-member helix in a unit cell of tetragonal phase along the crystallographic a or b direction. (e-f) 2D schematic of two endmembers of Li garnet oxides Li3Nd3Te2O12, LixA3B2O12 (3 𝑅𝑐 where 𝑅𝑖𝑗 is the distance between atoms 𝑖 and 𝑗. The ACSFs in this work include the radial and angular functions with forms: 2 𝐺𝑖𝑟𝑎𝑑𝑖𝑎𝑙 = ∑ 𝑒 −𝜂(𝑅𝑖𝑗−𝑅𝑠 ) ∙ 𝑓𝑐 (𝑅𝑖𝑗 ) 𝑗 𝑎𝑙𝑙 2 2 2 𝑎𝑛𝑔𝑢𝑙𝑎𝑟 𝐺𝑖 =2 1−𝜉 ∑ (1 + 𝜆𝑐𝑜𝑠𝜃𝑖𝑗𝑘 )𝜉 ∙ 𝑒 −𝜂(𝑅𝑖𝑗+𝑅𝑖𝑘+𝑅𝑗𝑘) ∙ 𝑓𝑐 (𝑅𝑖𝑗 ) ∙ 𝑓𝑐 (𝑅𝑖𝑘 ) ∙ 𝑓𝑐 (𝑅𝑗𝑘 ) 𝑗,𝑘≠𝑖 The radial function is a two-body potential, as a summation of Gaussian multiplied by the cutoff function. The parameters 𝜂 and 𝑅𝑠 determine the width and center shift of the Gaussian functions. In this work, 𝜂 has eight values (listed in Table 1) and 𝑅𝑠 was set to 0, hence there are eight radial SFs for each pair. Table 1 Radial function parameters[84]. GNo. η (Å-2) 1 0.003214 2 0.035711 3 0.071421 4 0.124987 5 0.214264 6 0.357106 7 0.714213 8 1.428426 The angular function is a three-body potential, a summation of the cosine function of the angle 𝜃𝑖𝑗𝑘 centered at atom 𝑖 multiplied by three cutoff functions. The parameter 𝜂 represents the distance from the central atom, the parameter 𝜉 describes the angular resolution, and the parameter 24 𝜆 shifts the maxima of the cosine function to 𝜃𝑖𝑗𝑘 = 0° and 𝜃𝑖𝑗𝑘 = 180°. The parameter values are listed in Table 2, the combination of the three parameters yields 18 angular SFs for each triplet. Table 2 Angular function parameters[84]. η (Å-2) ξ λ 0.000357 1 +1 0.028569 2 -1 0.089277 4 Figure 8b is an architecture of a 2-3-3-1 neural network [83], which shows how NN describes the potential energy surface (PES) of each atom. There are four layers, which are the left input layer, the right output layer, and in between two hidden layers. The three nodes in each hidden layer have (𝑘) no physical meaning. The value of node 𝑖 in layer k is noted as 𝑦𝑖 ( 𝑖 starts from 1 and 𝑘 starts (0) (0) 𝑇 from 0). Therefore, the input 𝑮 = (𝐺1 , 𝐺2 )𝑇 corresponds to 𝒚(0) = (𝑦1 , 𝑦2 ) , which are used (1) (1) (1) 𝑇 to compute the first hidden layer values 𝒚(1) = (𝑦1 , 𝑦2 , 𝑦3 ) and thereafter the second hidden (2) (2) (2) 𝑇 layer values 𝒚(2) = (𝑦1 , 𝑦2 , 𝑦3 ) , in the end, output the energy, 𝐸. The vector of 𝒚(k) (𝑘 ≥1) are calculated through the intermediate vector 𝒙(𝑘) , 𝑇 𝑇 𝒙(𝑘) = (𝒂(𝑘−1) ) 𝒚(𝑘) + (𝒃(𝑘−1) ) where 𝒂(𝑘−1) the weight matrix and 𝒃(𝑘−1) the bias weight matrix, the elements of which need to (𝑘) be optimized during the training process. And then the value 𝑦𝑖 is computed through (𝑘) (𝑘) 𝑦𝑖 = 𝑓 (𝑘) (𝑥𝑖 ) where 𝑓 is the activation function (also called transfer function) for the hidden layers, the effect of different activation functions will be discussed later in Chapter 4. For the output layer, a linear function is used. 25 2.3.5 Molecular dynamics simulation Most of the work presented in this dissertation is the result of the molecular dynamics (MD) simulations. In classical MD, the atoms are treated as classical particles and the motions of the particles follow Newton’s second law as in classical mechanics. The force is the derivative of the potential energy 𝑑2 𝒓𝑖 (𝑡) 𝑭𝑖 = −∇𝑈 = 𝑚𝑖 𝑑𝑡 2 where 𝑈 is the potential energy, 𝑖 represents the 𝑖th atom, and 𝑚𝑖 , 𝑭𝑖 , and 𝒓𝑖 are the mass, force, and coordinate of the 𝑖th atom at time 𝑡. Since 𝑭𝑖 of each step can be calculated based on the coordinates and potential energy, the acceleration of each atom, 𝒂𝑖 , is 𝑭𝑖 divided by 𝑚𝑖 . With the initial assignment of coordinates and velocities, the coordinates of the first step with a short time ∆𝑡 can be predicted. Combining the Verlet integration, the coordinates of the next step can always be predicted, and then the velocity can be calculated. 1 𝑟(𝑡 + ∆𝑡) = 𝑟(𝑡) + 𝑣(𝑡)∆𝑡 + 𝑎(𝑡)∆𝑡 2 2 𝑎(𝑡) + 𝑎(𝑡 + ∆𝑡) 𝑣(𝑡 + ∆𝑡) = 𝑣(𝑡) + ∆𝑡 2 This process of calculating force, acceleration, coordinate, and velocity is repeated until the required steps are performed. Forces in classical MD simulations normally come from interatomic potentials. If the force comes from the electronic structure calculations based on quantum mechanics, the MD simulation is called ab initial MD (AIMD) or first principles MD (FPMD). 2.4 Training interatomic potentials Potential parameters of IP models are obtained by fitting them to either experimental crystal structures or quantum mechanical training data. In this work, the potential parameters of Core and CS models were trained by fitting to lattice parameters and atomic positions of the 26 experimental structures in our previous work through the GULP package [85, 86]. Potential parameters for the ID model were derived from the force-matching method[87, 88] which fits interatomic potential forces to the first-principles forces of each atom. DFT-MD trajectories of binary oxides Li2O, ZrO2, La2O3, and Ta2O5 together with the lithium garnet oxides were used with the lattice parameters, atomic coordinates, and atomic forces included in the training, as shown in Figure 9a. DFTB and MLIP models were trained from the DFT-MD trajectories of lithium garnet oxides with energy, force, and cell-virial, see Figure 9b. The reliability of the potential parameters was assessed by root mean squared error (RMSE) of atomic force and cell- virial values from the trained potential against DFT. Figure 9 Structures of (a) binary oxides and (b) one member in lithium garnet oxide family Li6La3ZrTaO12 used for performing DFT-MD to train ID model, DFTB, and MLIP models. For the atomic charges, the average charge for each atom species of Core and ID models was obtained using the Density Derived Electronic and Chemical (DDEC) approach [89-91], the CS model used 80% of formal charge, the DFTB model used Mulliken atomic charge, and the 27 MLIP model is a mathematic model without charge. The details and values of potential parameters and atomic charges are listed in Chapters 3 and 4. 2.5 Methods for property calculations 2.5.1 Self-diffusivity As MD simulations record the time evolution of atoms concerning coordinates and velocities, time correlation functions are important methods to extract the property information. There are two types of correlation functions. The first one is the incoherent (also called self-part) van Hove correlation function (VHCF), which is the probability density of finding a particle in the vicinity of 𝒓 after traveling time 𝑡: 𝑁𝛼 𝛼 1 𝐺𝑖𝑛𝑐 (𝒓, 𝑡) = 〈∑ 𝛿[𝒓 + 𝒓𝛼𝑖 (0) − 𝒓𝛼𝑖 (𝑡)]〉 𝑁𝛼 𝑖=1 where 𝑁𝛼 the number of atoms of species 𝛼, and 𝒓𝛼𝑖 is the position of 𝑖th atom of species 𝛼. The angle brackets represent the ensemble average. 𝛼 The atomic diffusion was investigated by the spatial-Fourier transform of 𝐺𝑖𝑛𝑐 (𝒓, 𝑡), called intermediate scattering function or incoherent density correlation function[92], which was calculated according to 𝑁𝛼 𝛼 1 𝛼 𝛼 𝐹𝑖𝑛𝑐 (𝒌, 𝑡) = 〈 ∑ 𝑒 −𝑖𝒌∙𝒓𝑖 (0) 𝑒 𝑖𝒌∙𝒓𝑖 (𝑡) 〉 𝑁𝛼 𝑖=1 where 𝒌 is the wave vector of the simulation cell as 𝒌𝑛𝑥𝑛𝑦𝑛𝑧 = 𝑛𝑥 𝒂∗ + 𝑛𝑦 𝒃∗ + 𝑛𝑧 𝒄∗ where 𝒂∗ , 𝒃∗ , 𝒄∗ are reciprocal unit vectors and 𝑛𝑥 , 𝑛𝑦 , and 𝑛𝑧 are all integers. Figure 10a gives an example of the correlation time dependence of the function of a specific 𝒌. By fitting the curve with a (Kohlrausch-Williams-Watts) function[93] 28 𝐾𝑊𝑊 (𝑘)]𝛽(𝑘) 𝐼 𝐾𝑊𝑊 (𝑘, 𝑡) = 𝑒 −[Γ where Γ 𝐾𝑊𝑊 is the relaxation rate and 𝛽 is the stretching parameter conforming to 0 < 𝛽 < 1 [94]. 1 − 𝛽 is called coupling parameter by Ngai[95], which depends on the intermolecular interaction. The mean relaxation rate, 〈Γ〉, was obtained from Γ 𝐾𝑊𝑊 (𝑘)𝛽 〈Γ〉 = 𝐺(𝛽 −1 ) where 𝐺 is the Gamma function. The k dependence of relaxation rate 〈Γ〉 is shown in Figure 10a. Two diffusion mechanisms can be used to analyze the Li-ion self-diffusivity. If Γ is linear dependence of 𝑘 2 , it means that Li-ions diffuse continuously in the material and the Fickian diffusion model can be used to extract the self-diffusivity, 𝛤(𝑘) = 𝐷𝑘 2 where 𝐷 is the self-diffusivity. Otherwise, the jump-diffusion models such as Chudley-Elliot (CE), Hall-Ross (HR), and Singwi-Sjolander (SS) models can be used to extrapolate the relation time, jump length, and self-diffusivity[96]. As shown in Figure 10a, the vertical arrows indicate vibrate and the horizontal arrows indicate diffuse/jump to the next site, the atoms (orange balls) vibrate for a residence time 𝜏 and then jump to the next site with jump length 𝑙. In this work, the SS jump- diffusion model[96, 97] was applied, 1 𝑘 2 〈𝑙2 〉/6 Γ 𝑆𝑆 (𝑘) = 𝜏 1 + 𝑘 2 〈𝑙2 〉/6 where 𝜏 is the residence time, 〈𝑙2 〉 the mean square jump length square and the self-diffusivity can be calculated through 𝐷 = 〈𝑙2 〉/(6𝜏) 29 Figure 10 (a) Correlation time dependence of incoherent density correlation function. Relaxation rate as a function of 𝑘 2 and illustration of diffusion models. Correlation time dependence of (b) mean squared displacement (MSD) and (c) velocity autocorrelation function (VACF)[24]. The second group of the time correlation function is based on Einstein-Helfand (EH) formula[24, 92] 1 𝐾𝐸𝐻 = lim 〈|𝑨(𝑡) − 𝑨(0)|2 〉 𝑡→∞ 2𝑡 where 𝑨(𝑡) could be incoherent or coherent versions of position or charge-weighted position. For example, the incoherent mean squared displacement (MSD) is the version of the position 𝑁𝛼 1 〈|𝛿𝒓𝛼𝑖 (𝑡)|2 〉 = ∑〈|𝒓𝛼𝑖 (𝑡) − 𝒓𝛼𝑖 (0)|2 〉 𝑁𝛼 𝑖=1 where 𝑁𝛼 the number of atoms of species 𝛼, and 𝒓𝛼𝑖 is the position of 𝑖th atom of species 𝛼. And therefore the self-diffusivity can be calculated from the slope of MSD vs. time plot as shown in Figure 10b 30 1 𝐷 = lim 〈|𝛿𝒓𝛼𝑖 (𝑡)|2 〉 𝑡→∞ 6𝑡 2.5.2 Ionic conductivity Ionic conductivity is one of the most important properties of solid-state electrolytes. Here, we are introducing two methods, one is the coherent MSD of the charge-weighted EH equation[98] 2 𝑁 𝑁 1 1 𝜎= 𝑙𝑖𝑚 〈(∑ 𝑧𝑖 𝑒𝑟𝑖 (𝑡) − ∑ 𝑧𝑖 𝑒𝑟𝑖 (0)) 〉 3𝑘𝐵 𝑇𝑉 𝑡→∞ 2𝑡 𝑖=1 𝑖=1 where 𝑘𝐵 the Boltzmann constant, 𝑇 the temperature, 𝑉 the volume of the simulation cell, 𝑁 the total number of atoms in the simulation cell, 𝑧𝑖 and 𝒓𝑖 the charge and coordinate of the 𝑖th atom, 𝑒 the elementary electron charge, and 𝑡 the time. This coherent ionic conductivity is directly related to experimental measurement. The other method is based on the Green-Kubo (GK) formula ∞ 𝐾𝐺𝐾 = ∫0 〈𝐴̇(0) ∙ 𝐴̇(𝑡)〉dt which is mathematically equal to the 𝐾𝐸𝐻 formula[24]. Figure 10c illustrates the incoherent velocity autocorrelation function (VACF) as a GK formula, which yields the self-diffusivity of species 𝛼 as 𝑁𝛼 ∞ 1 𝐷= ∫ ∑〈𝒗𝛼𝑖 (𝑡) ∙ 𝒗𝛼𝑖 (0)〉𝑑𝑡 3𝑁𝛼 0 𝑖=1 where 𝑁𝛼 the number of atoms of species 𝛼 and 𝒗𝛼𝑖 is the velocity of the 𝑖th atom of species 𝛼. Similarly, the coherent charge current correlation function ∞ 𝑁 𝑁 1 𝜎= ∫ 〈∑ 𝑧𝑖 𝑒𝒗𝑖 (𝑡) ∙ ∑ 𝑧𝑖 𝑒𝒗𝑖 (0)〉 𝑑𝑡 3𝑘𝐵 𝑇𝑉 0 𝑖=1 𝑖=1 31 where 𝑁 the number of total atoms in the simulation cell and 𝒗𝑖 is the velocity of the 𝑖th atom. It can be separated into two parts: the longitudinal component 𝐿,𝑐 (𝒌, 𝑒2 1 𝐶 𝑡) = 〈(𝒌 ∙ 𝑱(𝒌, 0))(𝒌 ∙ 𝑱(−𝒌, 𝑡))〉 𝑉𝑘𝐵 𝑇 𝑘 2 and the transverse component 𝑇,𝑐 (𝒌, 𝑒2 1 𝐶 𝑡) = 〈(𝒌 × 𝑱(𝒌, 0)) ∙ (𝒌 × 𝑱(−𝒌, 𝑡))〉 𝑉𝑘𝐵 𝑇 2𝑘 2 where 𝑁 𝑱(𝒌, 𝑡) = ∑ 𝑧𝑖 𝒗𝑖 (𝑡)𝑒 −𝑖𝒌∙𝒓𝑖 (𝑡) 𝑖=1 is the collective charge current, 𝑒 the elementary electron charge, 𝑉 the volume of the simulation cell, 𝑘𝐵 the Boltzmann constant, 𝑇 the temperature, 𝑞𝑛 and 𝒗𝑛 are the charge and velocity of the 𝑛 th atom, respectively. The coherent charge current correlation can be transformed to the frequency domain through the Fourier Transform (FT) ∞ 𝑆𝐿,𝑇 (𝒌, 𝜔) = ∫ 𝐶𝐿,𝑇 (𝒌, 𝑡)𝑒 𝑖𝜔𝑡 𝑑𝑡 0 where 𝜔 is the frequency. The ionic conductivity can be extrapolated at the long wavelength and zero frequency 𝜎 = lim 𝑆𝑇 (𝒌, 𝜔) 𝜔→0 Even though the ionic conductivity calculated through the GK formula is mathematically equivalent to that from the EH formula, it can be seen that GK starts with a strong oscillation and drop to 0 in a few picoseconds, which suggests that the MD trajectory should be saved more frequently to capture the oscillation features while eliminating or reducing the dependence of the saving frequency of integral. 32 In order to illustrate that the ionic motion is positively cooperated instead of uncorrelated, we also compared the coherent ionic conductivity with the Nernst-Einstein (NE) ionic conductivity that calculated from self-diffusivity 𝐷𝐿𝑖 𝜎𝑁𝐸 = (𝑧𝐿𝑖 𝑒)2 (𝑐 ) 𝑘𝐵 𝑇 𝐿𝑖 where 𝑐𝐿𝑖 the lithium concentration. The ratio of 𝜎𝑁𝐸 to 𝜎 is known as the Haven ratio: 𝜎𝑁𝐸 𝐻𝑅 = 𝜎 The Haven ratio is 1 for ideal solutions where the motion of ions is independent and uncorrelated. 2.5.3 Dielectric constant Since the conductivity is related to dielectric function as 𝜔𝜀0 [𝜀𝐿,𝑇 (𝒌, 𝜔) − 1] = 𝑖𝜎𝐿,𝑇 (𝒌, 𝜔) where 𝜀0 is the vacuum permittivity, 𝜔 is the frequency, and 𝜀𝐿,𝑇 (𝒌, 𝜔) is the dielectric constant. Therefore, the real part of the dielectric function can be obtained through the imaginary part of conductivity as 𝑅𝑒[𝜀𝑇 (𝒌, 𝜔) − 1] = −𝐼𝑚[𝜎𝑇 (𝒌, 𝜔)]/(𝜔𝜀0 ) 2.5.4 Elastic modulus The elastic moduli were investigated by the momentum correlations that can be separated into the longitudinal and transverse components as 1 𝐶𝐿,𝑚 (𝒌, 𝑡) = 〈(𝒌 ∙ 𝑱(𝒌, 0))(𝒌 ∙ 𝑱(−𝒌, 𝑡))〉 𝑘2 1 𝐶 𝑇,𝑚 (𝒌, 𝑡) = 〈(𝒌 × 𝑱(𝒌, 0)) ∙ (𝒌 × 𝑱(−𝒌, 𝑡))〉 2𝑘 2 where 33 𝑁 𝑱(𝒌, 𝑡) = ∑ 𝑚𝑛 𝒗𝑛 𝑒 −𝑖𝒌∙𝒓𝑛 (𝑡) 𝑛=1 is the collective momentum, 𝑚𝑛 the atomic mass of the 𝑛th atom, 𝒗𝑛 the velocity of the nth atom. The Fourier Transform (FT) of longitudinal and transverse forms are ∞ 𝑆𝐿,𝑇 (𝒌, 𝜔) = ∫ 𝐶𝐿,𝑇 (𝒌, 𝑡)𝑒 𝑖𝜔𝑡 𝑑𝑡 0 where 𝜔 is the frequency. Figure 11 gives an example of the momentum correlations of both longitudinal and transverse components of Li7La3Zr2O12 using the induced dipole (ID) model at 1100 K. Figure 11 (a) Time domain and (b) frequency domain of momentum correlation function of the induced dipole (ID) model at 1100 K. Black lines are Lorentz fit to locate the peak position (dash lines). We located the peaks (𝜔) of 𝑆𝐿,𝑇 (𝒌, 𝜔) and assumed linear dispersion of 𝜔 on 𝑘 to obtain longitudinal and transverse sound velocities 𝑣𝐿 and 𝑣𝑇 . For the cubic phase, the elastic constants 𝐶11 , 𝐶44 and 𝐶12 were calculated by the following formulas[99]: (𝑣𝐿2 𝜌)[100] = 𝐶11 (𝑣𝑇2 𝜌)[001] = 𝐶44 34 1 (𝑣𝑇2 𝜌)[110] = (𝐶11 − 𝐶12 ) 2 Note that the (𝑣𝑇 )[110] was averaged from the two peaks. For the tetragonal phase, 𝐶11 , 𝐶44 and 𝐶12 were calculated in the same way as in the cubic phase, while 𝐶33 and 𝐶66 were calculated by the following formulas[99]: (𝑣𝐿2 𝜌)[001] = 𝐶33 (𝑣𝑇2 𝜌)[100] = 𝐶66 We assumed that 𝐶12 = 𝐶23 as they differ by only 2% from Deng et al[100]. The bulk and shear modulus were calculated according to [100] 1 𝐵 = [(𝐶11 + 𝐶22 + 𝐶33 ) + 2(𝐶12 + 𝐶23 + 𝐶31 )] 9 1 𝐺= [(𝐶 + 𝐶22 + 𝐶33 ) − (𝐶12 + 𝐶23 + 𝐶31 ) + 3(𝐶44 + 𝐶55 + 𝐶66 )] 15 11 In the cubic phase, 𝐶11 = 𝐶22 = 𝐶33 , 𝐶12 = 𝐶23 = 𝐶31 , and 𝐶44 = 𝐶55 = 𝐶66 . In the tetragonal phase, 𝐶11 = 𝐶22 , 𝐶44 = 𝐶55 , and 𝐶23 = 𝐶31 . 2.5.5 Grain boundary formation energy Grain boundary (GB) formation energy, 𝛾𝐺𝐵 , is the energy needed to form the grain boundary from the bulk material per unit area according to 1 𝛾𝐺𝐵 = (𝐸 − 𝑛𝐸𝑏𝑢𝑙𝑘 ) 2𝐴 𝐺𝐵 where 𝐴 the area of GB plane, 𝐸𝐺𝐵 the total energy of the GB structure, 𝑛 the number of a unit cell of the bulk material in the GB structure, and 𝐸𝑏𝑢𝑙𝑘 is the energy of the unit cell of the bulk material. 35 CHAPTER 3 Model comparison 3.1 Motivation As introduced in Chapter 2, Li7La3Zr2O12 (Li7) is a special member of the lithium garnet family LixLa3Zrx-5Ta7-xO12 (x=5-7) due to the phase transition at high temperatures. Hence, it has been widely studied using interatomic potential (IP) models. For example, Adams et al.[101] performed MD simulation on Li7 with the core model and Morse potential. Chen et al.[102] used the core model to study the diffusion characteristics of Li7 based on the Buckingham potential. In the previous work of our group[85, 103], we used the core model to study the size effect on the MD simulation of Li7 and used a polarizable core-shell (CS) model to study the local structure and dynamics of Li7. Burbano et al.[104] used the polarizable induced dipole (ID) model to perform an MD simulation of Li7 to study the mechanism of its low ionic conductivity. However, no work has been done to compare these models. In this chapter, we compared the accuracy of three models, i.e., Core, CS, and ID models in reproducing the phase transition, residence time, jump length, self-diffusivity, ionic conductivity, dielectric constant, and elastic moduli of fast ion conductor Li7. 3.2 Computational details 3.2.1 Potential parameters In the Core model, ionic interactions include the long-range Columbic potential and short- range Buckingham potential. The charges and potential parameters of the Core model were taken from our group’s previous study[85]. CS model interactions include the long-range Columbic potential, short-range Buckingham and Morse potential, as well as the Dick-Overhauser core-shell polarizable potential for the O atom[105]. The harmonic spring constant between the oxygen core and shell is 19.64 eV Å-2. The O shell (OS) mass was set to 0.2 a.u. The charges of each species were set as 80% of their formal charges. 36 In the ID model, we assigned the electronic polarizability to the O, La, and Zr species. Potential parameters were derived from the force-matching method. The force-matching was carried out using the FIST code as implemented in CP2K[106, 107]. The training set was from the first-principles MD simulations based on DFT, which were performed using the generalized gradient approximation (GGA-PBEsol)[108] and the projector augmented wave (PAW)[109] method as implemented in the Vienna Ab Initio Simulation Package (VASP)[110]. The unit cell was simulated using the NPT ensemble at 1200 K. The energy cutoff of 500 eV was used for the plane wave and the Brillouin zone was sampled at the Γ-point only. The production time was 3 ps with a time step of 1 fs. The trajectories after the first 0.15 ps were used to perform the force- matching. We extracted 190 configurations, i.e. one configuration every 15 fs. The obtained average structure was used for the subsequent static calculation. The average charge for each atom species was obtained using the Density Derived Electronic and Chemical (DDEC) approach [89- 91]. The charge of each species and the potential parameters of the Core, CS, and ID model are listed in Table 3, Table 4 and Table 5, respectively. Table 3 Interatomic potential parameters for the Core model. Buckingham pair interaction with O Species q (e) A (eV) ρ (Å) C (eV Å6) La 2.50 2075.26 0.326 23.25 Li 2.65 1087.29 0.260 0.00 O 1.00 4870.00 0.267 77.00 Zr -1.65 1650.32 0.311 5.10 Table 4 Interatomic potential parameters for the CS model. Buckingham potential Morse potential Species q (e) Interaction pair A (eV) ρ (Å) C (eV Å6) E0 (eV) k (Å-1) r0 (Å) La 2.4 La-OS 1828 0.328 0 0.066 1.471 2.423 Zr 3.2 Zr-OS 2135 0.302 0 0.241 2.742 2.016 Li 0.8 Li-OS 370 0.291 0 0.081 0.870 1.937 O -2.116 Li-Li 1470 0.239 0 OS 0.516 37 Table 5 Interatomic potential parameters for the ID model. Polarization potential Buckingham pair interaction with O Species q (e) α (Å3) Ion-pair bij (Å-1) cij (Å-1) A (eV) ρ (Å) C (eV Å6) La 2.54 1.71 La-O 3.344 2.315 2512.92 0.314 0 Zr 2.93 0.76 Zr-O 3.155 1.233 3034.17 0.282 0 Li 1 - Li-O 3.772 1.729 507.79 0.298 0 O -1.71 1.74 3.2.2 Molecular Dynamics simulations MD simulations were performed with the DL_POLY Classic package[111] for both Core and CS models, while for the ID model, it was carried out using the FIST code as implemented in CP2K[106, 107]. A 2×2×2 supercell of Li7La3Zr2O12 was used for all three models. There are 1536 atoms in the Core and ID model and 2304 atoms (including O shell) in the CS model. For Core and CS modes, the velocity-Verlet integration was employed with a time step of 1 and 0.25 fs, respectively. The constant number, stress, and temperature (NσT) ensemble (Berendsen thermostat and barostat relaxation times were set to be 0.1 and 0.5 ps, respectively) simulations were performed for 75 ps and the last 35 ps were averaged to obtain the lattice parameters. The equilibration time of the constant number, volume, and energy (NVEe) ensemble was set to 100 ps with the atomic velocities rescaled every 0.1 ps. Subsequent NVE ensemble simulations were performed from 500 ps at high temperatures and to 70 ns in the case of low mobility of Li+ in the tetragonal phase at 800 K. For the ID model, the constant number, pressure, and temperature (NPT) ensemble was performed for 50 ps and the average lattice parameters were extracted from the last 40 ps. The Nose thermostat and barostat relaxation times were set to 0.05 and 0.25 ps, respectively. The NVE equilibrium (NVEe) process was 50 ps and the maximum accepted deviation of the temperature 38 from the desired target temperature was 50 K. The NVE process was set from 200 ps to 10 ns for different temperatures. The time step was 1fs for NPT, NVEe, and NVE process. 3.3 Results and discussion 3.3.1 Lattice parameters First, we checked if all three IP models could reproduce the tetragonal to the cubic phase transition of Li7 observed experimentally. A wide variety of transition temperatures were reported for Li7 experimentally ranging from 373 to 923 K[23, 53, 112-114]. The difference could be resulting from the impurity/dopant (i.e., Al3+) incorporation and/or H2O or CO2 exposure during material synthesis and characterization. For example, Matsuda et al.[53] found that the phase transition temperature of Li7 decreases as the Al content increases and the reversible phase transition between the tetragonal phase and cubic phase occurs at 640 °C (913 K) for pure Li7 without Al. Larraz et al.[115] reported that the presence of H2O could induce the tetragonal to cubic phase transition due to the H+/Li+ exchange mechanism. On the other hand, undoped Li7 not exposed to humidity at any moment undergoes a reversible tetragonal to cubic phase transition at about 645 °C (918 K). Given these observations, we only plotted the experimental values from Matsuda et al.[53] that applies to pure Li7 in Figure 12, along with lattice parameters obtained from the Core, CS, and ID models. Lattice parameters of all three models increase with increasing temperature as expected in both tetragonal and cubic phases, with the phase transition occurring at about 900 to 1000 K. Both Core and ID models show thermal expansion and phase transition behaviors consistent with experiments while the CS model shows lower lattice parameters in 𝑎. 39 Figure 12 Lattice parameters from Core, CS, and ID models as a function of temperature. XRD experimental values from Matsuda et al. [33] are shown for comparison. 3.3.2 Self-diffusivity, jump length, and residence time Next, we compared the long-range diffusion parameter, i.e., self-diffusivity, and short- range parameters such as jump length and residence time. Figure 13 gives an example of the incoherent density correlation function of the ID model at 1100 K for selected 𝑘. It can be seen that for species La, Zr, and O, there is a rapid decrease in the ps range followed by a stable non- zero plateau, which suggests that these species only vibrate in their sites and provide a framework for Li atoms. In contrast, species Li exhibits an exponential decay from 1 to 0, especially faster at high k values, suggesting Li atoms diffuse in the structure. Ignoring data for the first 5 ps that 𝛼 undergo short-time dynamics, we fit the 𝐼𝑖𝑛𝑐 (𝑘, 𝑡) curve of Li species using stretched exponential Kohlrausch-Williams-Watts function. Figure 14 shows an example of the stretching parameters 𝛽 and mean relaxation rate as a function of 𝑘 2 of ID model at 1100 K. The behavior of 𝛽 reveals that the Li transport in the material is a complex cooperative dynamic. As shown in Figure 14a, 𝛽 decays from ~1 at small 𝑘 values and converges at high 𝑘 values. Similar behavior can also be found at all the other temperatures we studied. At small 𝑘 40 (longer distance), the interaction between any pair of atoms is weak so their motion is uncooperative or not coupled. At higher 𝑘 (shorter distance), the interaction is cooperative or coupled. The plateau starts at 𝑘 2 ~2.0 Å−2, corresponding to a cooperative length of ~4.4 Å, which is the approximate distance between octahedral sites across the tetrahedral site[44]. Figure 14b exhibits the fit of the Fickian diffusion model and Singwi-Sjolander (SS) jump-diffusion model to the MD data. Figure 13 Incoherent density correlation function of species La, Zr, O, and Li for selected 𝑘 values of the ID model at 1100 K. Figure 14 (a) The stretching parameter 𝛽 and (b) the mean relaxation rate 〈Γ〉 as a function of 𝑘 2 of ID model at 1100 K. 41 The residence time (𝜏) and jump length (𝑙 ) obtained from the SS model at different temperatures are shown in Figure 15. All Core, CS, and ID models show a sharp increase in residence time (Figure 15a) at the phase transition temperature. CS and ID models have close values in the cubic phase. Activation energies obtained from an Arrhenius fit for Core, CS, and ID models are 0.26, 0.26, and 0.31 eV, respectively. In the tetragonal phase, all three models show high activation energy values above 1.7 eV. The mean jump length (Figure 15b) obtained from Core, CS, and ID models at higher temperatures (cubic phase) are very close, i.e., 1.4~1.7 Å, which is roughly the distance between 24d (tetrahedral) and 48g (octahedral) sites in the cubic phase[103]. The jump length exhibits a significant increase at the phase transition temperature. In the tetragonal phase (800 K), the jump lengths for CS and ID models are about 2~2.5 Å, which is roughly the distance between 8a (tetrahedral) and 32g (octahedral) sites[45, 103]. Figure 15 (a) Mean residence time (𝜏) of Li diffusion in LLZO obtained from Core, CS, and ID models in both cubic (c-Li7) and tetragonal (t-Li7) phases. Solid lines are an Arrhenius fit. (b) Mean jump length of Li diffusion obtained from Core, CS, and ID models. Figure 16 compares Li-ion self-diffusivity in Li7 of the Core, CS, and ID models along with other MD simulation results from the literatures. Our calculated values from the Core model are close to those from Jalem et al. using first-principles molecular dynamics (FPMD) simulation 42 with 1×1×1 cell size[116], while values from the CS and ID models are a few times higher. The activation energy of the Core model for the cubic phase is 0.24 eV, which is also close to that from Jalem et al. at the temperature range of 873 to 1473 K, 0.26 eV. The activation energy of CS and ID models for the cubic phase are 0.18 and 0.20 eV, respectively. The self-diffusivity obtained from the CS model at high temperatures is close to that from Burbano et al. also using an ID model but with a smaller activation energy of 0.14 eV[104]. The difference is likely to result from the lower phase transition temperature of 623 K, as opposed to 900-1000 K, in their work. Experimentally, the Li self-diffusivity at ambient temperature from Kuhn et al. using the spin- lattice relation Nuclear Magnetic Resonance (NMR) measurement is on the order of 10-14 cm2/s[117]. Although our simulation was unable to calculate Li self-diffusivity at room temperature, extrapolation of our data to the room temperature yields values on the order of 10-25 ~10-20 cm2/s, which are much lower than the experimental result. Figure 16 Self-diffusivity as a function of temperature for Core, CS, and ID models in both cubic (c-Li7) and tetragonal (t-Li7) phases. Solid lines are Arrhenius fit. The FPMD simulation data by Jalem et al.[116] and classic MD simulation data by Burbano et al.[104] are shown for comparison. 43 3.3.3 Ionic conductivity and dielectric constant Figure 17a gives an example of the transverse component of coherent charge current density correlation function 𝐶 𝑇,𝑐 (𝒌, 𝑡) of the ID model at 1000 K for the three smallest 𝑘. It was found that the Fourier transform of the longitudinal component is a peak that goes to zero at zero frequency (not shown) but the transverse component, 𝑆 𝑇 (𝒌, 𝜔) shown in Figure 17b, has a non- zero intercept suggesting that Li7 is an ionic conducting material. The ionic conductivity was extrapolated from the transverse component. To manifest that 𝑆 𝑇 (𝒌, 𝜔) is 𝑘 independent, it is important to separate out all transverse outputs as shown in Figure 17c. The first number in the legend denotes the 𝑘 value in Å−1 while the second and third series of three digits denote the 𝑘 direction and its transverse direction. For example, the transverse polarization of vector [100] has two directions, i.e., [010] and [001]. It is apparent that the outputs of the two 𝑘 vectors are similar to each other, which indicates that 𝑆 𝑇 (𝒌, 𝜔) is independent of 𝑘 vector. Thus, we averaged all six values of the first 𝑘 vector and took the result of the smallest frequency as the ionic conductivity. Figure 17 (a) Time and (b) frequency domain of the transverse coherent charge current density correlation function for the first three 𝑘 values of the ID model at 1000 K. (c) All transverse outputs for first two 𝑘 vectors at low frequencies. Figure 18 compares the ionic conductivities of Core, CS, and ID models at different temperatures along with the experimental results from the literatures. It is obvious that all three models have a kink at about 900 K corresponding to the phase transition temperature. The 44 experimental measured ionic conductivities from Wang et al.[23] and Matsui et al.[118] also exhibit a sharp decrease at phase transition temperature between 600 °C (873 K) and 650 °C (925 K). Activation energies of the cubic phase for Core, CS, and ID models, i.e., Core-c, CS-c, ID-c, are 0.28, 0.13, and 0.17 eV, respectively, while values from Wang et al.[23] and Matsui et al.[118] are 0.13 and 0.12 eV, respectively. For the tetragonal phase, our simulations are unable to reliably calculate the ionic conductivities at lower temperatures due to the low Li+ mobility. A simple extrapolation of our data will yield lower values than those from Wang et al.[23], Matsui et al.[118], and J. Awaka et al.[45]. Figure 18 Ionic conductivity as a function of temperature for Core, CS, and ID models. Literature data from experimental measurements by Wang et al. [4], Matsui et al. [37], and Awaka et al. [6] are shown for comparison. Figure 19a gives an example of the imaginary part of transverse conductivity 𝜎𝑇 (𝒌, 𝜔) over 𝜀0 for the first 𝑘 value (0.24 Å-1) as a function of frequency for the ID model at 800 K. The linear fit (inset of Figure 19a) yields a slope of 𝜀𝑟 (0) − 1 where 𝜀𝑟 (0) is the static (zero frequency) dielectric constant. The results of 𝜀𝑟 (0) of all three models at different temperatures are summarized in Figure 19b. In the tetragonal phase, the dielectric constants of the three models 45 decrease with increasing temperatures, as expected for typical materials due to thermal agitation. On the other hand, the dielectric constants of the three models increase with increasing temperatures and we believe this is due to mobile Li ions. The dielectric constant of tetragonal and cubic Li7 is not known but S. Narayanan et al.[119] reported that Li5+2xLa3Nb2-xYxO12 (x=0.25, 0.5, and 0.75) have dielectric constants around 60 below room temperature. Figure 19 (a) The imaginary part of transverse conductivity 𝜎𝑇 (𝒌, 𝜔) over 𝜀0 as a function of frequency for the ID model at 800 K. The inset shows details at low frequencies, and the solid line is a linear fit. (b) The dielectric constant at zero frequency 𝜀𝑟 (0) as a function of temperature for Core, CS, and ID models. 3.3.4 Bulk modulus and shear modulus Last but not the least, the mechanical properties of solid-state electrolytes (SSEs) also play an important role in the fabrication and operation process of all-solid-sate LiBs. According to Monroe and Newman, the shear modulus of SSEs that are two times higher than that of Li metal would suppress Li dendrite growth in the SE/Li metal interface[120]. For this reason, we also calculated and compared the bulk modulus (𝐵) and the shear modulus (𝐺) of the three models, as shown in Figure 20. It can be seen that Li7 exhibits a high shear modulus, which is way higher than that of Li metal (~8.5 GPa). Both 𝐵 and 𝐺 decreases with increasing temperature, suggesting 46 the material softens upon heating. The Core model has the highest values of both 𝐵 and 𝐺, while the CS and ID models have similar elastic moduli. Experimental values at 300 K from Ni. et al.[121] are presented in Figure 20 and they are similar to our computational values. Figure 20 Elastic modulus as a function of temperature for Core, CS, and ID models. Experimental values from Ni et al. [39] are shown for comparison. 3.4 Conclusions In this chapter, we compared one non-polarizable (core/Core) and two polarizable (core- shell/CS and induced-dipole/ID) models for the molecular dynamics simulation of a fast-ion conductor, Li7La3Zr2O12 (Li7). Overall, we think it is difficult to explicitly point out which model is the best, as all three models can reproduce various structural and dynamic properties (lattice parameters, diffusivity, conductivity, dielectric constants, and elastic moduli) of Li7 reasonably well. Due to various limitations of existing codes, we used different methods to obtain potential parameters, e.g., arbitrary vs first-principles charge, Buckingham vs Buckingham+Morse short- range interaction, empirical fit vs force-matching, etc. We believe this also complicates the direct comparison of different potential models. On the other hand, we think the ID model combined with the force-matching approach is the most advanced one as it includes the polarization through 47 atomic polarizability and takes from first-principles results as the input. The CS model includes the polarization through a spring-connected shell where parameters of shell charge and spring constant are almost entirely empirical. Looking forward, we think that the ID model can be improved by incorporating the anisotropy in the atomic polarizability and fluctuating charges, where the charge will vary with the local chemical environment. Furthermore, only pair potentials were considered in these three models while many-body potentials will be more accurate to describe the mixed metallic/covalent nature of the bonding. 48 CHAPTER 4 Study of diffusion and conduction in lithium garnet oxides LixLa3Zrx-5Ta7- xO12 (x=5-7) by machine learning interatomic potential 4.1 Motivation As ionic conductivity is one of the most important properties of solid electrolytes, many experimental results have been reported for various compositions of LixLa3Zrx-5Ta7-xO12 (x=5-7), including both Al contaminated (through alumina crucibles), Al-doped, and nominally Al-free phases. Figure 21 summarized some work of Al-free and Al-contaminated experiments. For example, Wang et al.[22] and Li et al.[52] used alumina crucibles and found that ionic conductivity reaches a maximum at x=6.7 (0.96×10-3 S/cm) and x=6.4 (1.0×10-3 S/cm), respectively. The difference between these two is that the former is for the bulk and the latter for the total. On the other hand, Buschmann et al.[46], Inada et al. [48], Matsuda et al. [54] and Yi et al.[51] synthesized Al-free compositions of LixLa3Zrx-5Ta7-xO12 and found that ionic conductivity reaches a maximum at x=6 (2.6×10-4 S/cm), x=6.5 (6.1×10-4 S/cm), x=6.6 (4.7×10-4 S/cm), and x=6.7 (1.03×10-3 S/cm), respectively. Furthermore, Kataoka et al.[19] synthesized single crystals and found x=6.6 has the highest ionic conductivity of 1.1×10-3 S/cm. Thus, it appears that it is controversial as to what is the composition for the maximum conductivity in LixLa3Zrx-5Ta7-xO12 even in Al-free samples. In recent years, the machine learning interatomic potentials (MLIP) based on the artificial neural network (NN) model has received considerable attention as they show a combination of DFT accuracy and IP efficiency[79, 122-124]. In this chapter, we applied MLIP-based MD simulations to investigate lithium garnet series LixLa3Zrx-5Ta7-xO12 (x=5, 6, 6.5, 6.75, and 7), named Li5, Li6, Li6.5, Li6.75, and Li7 afterward. First, we compare the force/virial errors of MLIP, along with other commonly applied methods such as Core, ID, and DFTB, and check how well they can reproduce the self-diffusivity and ionic conductivity of Li6. Second, we present the MLIP 49 error against DFT for all compositions of Li5, Li6, Li6.5, Li6.75, and Li7 and show how this MLIP is approaching DFT accuracy. Third, we examine the temperature dependence of diffusion and conduction of various compositions including the defect structure. Finally, we turn to composition dependence and focus on the question of what composition gives the highest ionic conductivity, by considering both experimental literature and present computational results. Figure 21 Summary of composition dependence of ionic conductivity of LixLa3Zrx-5Ta7-xO12 (x=5- 7) from the literature [19, 22, 46, 48, 51, 52, 54]. SC: single crystal. 4.2 Construction of MLIP, Core, ID, and DFTB potentials 4.2.1 Construction of machine learning interatomic potential The MLIP potential of lithium garnet series LixLa3Zrx-5Ta7-xO12 was trained with the DFT- MD trajectories of Li5, Li6, and Li7 with energy, force, and cell-virial included in the training. DFT-MD simulations for all compositions were performed using the Vienna ab initio simulation package (VASP)[110] within the projector augmented-wave (PAW)[109] approach using the PBEsol[108] exchange-correlation functional. The plane-wave energy cutoff was 500 eV and a 50 single Γ point was sampled. First, we performed MD simulations, using an NPT ensemble, at 1200 K for 3 ps with a 1fs time step. For each composition, configurations were sampled every 4 fs after skipping the first 100 fs. 90% of this DFT dataset was used for training and 10% was used for validation. In the meanwhile, we used the average structure of the DFT-MD trajectory to obtain the atomic charges in the DDEC6 scheme[91], which are listed in Table 6, to be used in the calculation of the ionic conductivity for the MLIP model. Table 6 DDEC6 charge of each composition. Li5 Li6 Li6.5 Li6.75 Li7 La 1.87 1.88 1.88 1.84 1.89 Ta 2.32 2.33 2.33 2.23 — Zr — 2.23 2.23 2.18 2.25 O -1.20 -1.26 -1.29 -1.28 -1.32 Li 0.83 0.82 0.82 0.81 0.81 MLIP potential was trained with the SNU Interatomic Machine-learning PotentiaL packagE-version Neural Network (SIMPLE-NN) package[122], with atom-centered symmetry functions (ACSF)[79] as a descriptor and neural network (NN) (two hidden layers with 30 nodes in each layer) as the regressor. To develop the ideal potential for the Li garnet series LixLa3Zrx- 5Ta7-xO12, we started by studying the potential parameter effect on Li5, weighting factors, learning rate, and transfer functions. The criterion of the fitting is based on the loss function or the root mean square error for energy, force, and stress, which can be defined as (angle bracket represents the average over all frames): 𝐸𝐷𝐹𝑇 −𝐸𝐼𝑃 2 (𝑓𝐷𝐹𝑇 −𝑓𝐼𝑃 )2 (𝑆𝐷𝐹𝑇 −𝑆𝐼𝑃 )2 𝐹𝐿𝑜𝑠𝑠 = 𝑤𝐸 〈( ) 〉 + 𝑤𝑓 〈 〉 + 𝑤𝑆 〈 〉, 𝑁 3𝑁 6 where 𝑤𝐸 , 𝑤𝑓 , and 𝑤𝑆 are the weighting factors of energy, force, and stress. In simple-NN, the weighting factors of energy and force were set as default, i.e., 𝑤𝐸 = 1 and 𝑤𝑆 = 0.00001. In other words, the goal of the fitting is to minimize the force term, 51 𝑁 3 (𝑓𝐷𝐹𝑇 − 𝑓𝐼𝑃 )2 1 = ∑ ∑[𝑓𝐷𝐹𝑇 (𝑖, 𝑘) − 𝑓𝐼𝑃 (𝑖, 𝑘)]2 3𝑁 3𝑁 𝑖=1 𝑘=1 where 𝑁 is the number of atoms, and 3 means x, y, and z in three directions. First, we compared different values of 𝑤𝑓 (i.e., 0.05, 0.1, and 0.5) and finally decided to set 𝑤𝑓 to 0.1. We also compared two different activation functions used in the hidden layers, i.e., the 1 hyperbolic tangent 𝑓(𝑥) = 𝑡𝑎𝑛ℎ(𝑥) and the sigmoid function 𝑓(𝑥) = . Figure 22 1+𝑒𝑥𝑝(−𝑥) shows the root mean squared error (RMSE) of atomic forces and cell virials. It is found that the sigmoid function has a smaller RMSE compared to tanh function. Figure 22 The comparison of atomic forces and cell virials between MLIP (scatter points) against DFT (black lines) for Li5La3Ta2O12 using tanh (top) and sigmoid (bottom) activation functions with the same learning rate and 𝑤𝑓 . The NN potential parameters were optimized through the loss function for the training set i.e., configurations of DFT-MD). However, these structures contain various deformations, surfaces, and defects, meaning the descriptor is highly inhomogeneous and biased. This unbalanced training could significantly undermine the accuracy and reliability of NN potential, and therefore the Gaussian density function (gdf) has been added to cure the training bias and increase the 52 transferability of the potential[122, 125]. Figure 23 compared the lattice parameters of Li5 from MLIP-MD using sigmoid and sigmoid+gdf activation functions in the training, it can be seen that the lattice parameters calculated from sigmoid+gdf exhibit a linear dependence of temperature. Figure 23 Lattice parameter of Li5La3Ta2O12 using different activation functions, i.e., sigmoid and sigmoid+gdf. 4.2.2 Construction of Core, ID, and DFTB potentials The Core model potentials were taken from the work of Chen et al.[126], which was fitting from the experimental data. ID model (80 % of formal charge) of Li6 was trained with CP2K[106, 107] on the DFT-MD trajectory including only forces but all configurations, as we did previously in Chapter 3 for Li7. DFTB parameterization for the Li-La-Zr-Ta-O chemical space was obtained with the Tight-binding Approximation-enhanced Global Optimization (TANGO) method[127] with CP2K as the DFT package including both energy and force. 4.3 Computational and experimental details 4.3.1 Molecular Dynamic Simulations As mentioned in the Motivation, there was Li loss during the high-temperature synthesis and arbitrary amounts of lithium precursor were often added to compensate for this. It is thus worthwhile comparing the normal (nominal composition) to defect structures with Li loss. In this 53 work, in addition to Li5, Li6, Li6.5, Li6.75, and Li7, we also added the Li5.88La3Zr1.5Ta0.5O11.69 (10 mol% Li2O off of Li6.5, named Li6.5defect afterward). MD simulations of the Core and MLIP models were performed using Large Scale Atomic/Molecular Massively Parallel Simulator (LAMMPS)[128] with 2×2×2 supercells. MD simulations of the ID and DFTB models were performed using CP2K with 2×2×2 and 1×1×1 supercells, respectively. Details of MD simulations including the ensemble, thermostat, barostat, etc. are listed in Table 7, along with the approximate simulation cost. DFT-MD simulation of 2×2×2 supercells is expected to be 5 or 6 orders of magnitude slower than MLIP. Table 7 MD simulation details. The approximate simulation time of four models for 1 ns of trajectory is shown (Core, ID, and MLIP models were using 64 CPUs while the DFTB model was using 28 CPUs). It is worth noting that the DFTB model is using 1×1×1 while other models are using 2×2×2 supercells. Ensemble Parameters Core ID DFTB MLIP thermostat (ps) 0.05 0.1 0.1 0.1 barostat (ps) 0.25 0.5 0.5 0.5 NPT timestep (fs) 1 1 1 1 lattice parameter average range (ps) 30 – 50 30 – 50 30 – 75 30 – 50 thermostat (ps) — — 1 — NVT timestep (fs) — — 1 — trajectory length (ns) — — 0.5 – 3 — NVE thermostat (ps) 0.1 0.1 — 0.1 equilibration trajectory length (ps) 50 50 — 50 timestep (fs) 1 1 — 1 NVE trajectory length (ns) 0.5 – 300 0.5 – 8 — 0.5 – 40 approximate simulation time per ns NVT/NVE 0.3 40 170 20 trajectory (h) 4.3.2 Experiments Powders of Li6 and Li6.5 were prepared with a solid-state reaction where stoichiometric quantities of LiOH·H2O (Alfa Aesar, 98%), La2O3 (Alfa Aesar, 99.9%), ZrO2 (Alfa Aesar, 99.7%) and Ta2O5 (Alfa Aesar, 99.85%) were used as raw materials. La2O3 powders were heated at 900 °C for 12 h in a MgO crucible (all crucibles and lids used in this work were MgO to avoid Al 54 contamination). An extra 10 wt% LiOH·H2O was added to compensate for lithium loss during calcination. The powders were wet-milled in a roller mixer for 12 h in polyethylene jars filled with 2-propanol and then the slurry was dried by infrared heating. The mixed powders were calcined at 1000 °C for 12 h with a heating and cooling rate of 2 °C/min in MgO crucibles covered by MgO lids. Afterward, another extra 10 wt% or 5 wt% LiOH·H2O were added into the calcined powders of Li6 and Li6.5, respectively and ball-milled again, and finally pressed into 13mm diameter pellets and sintered at 1150 °C for 12 h. The silver paste was applied to both sides of the pellets as blocking electrodes and heated at 700 °C. Impedance tests at a temperature range of 22 – 700 °C were performed using the same setup as in our previous work[23] with airflow. A sinusoidal voltage with an amplitude of 10 mV was applied for the frequency range of 3 MHz to 1 Hz. 4.4 Results and discussion 4.4.1 Comparing Different Models for Li6 Before studying the whole series of Li5, Li6, Li6.5, Li6.75, Li7, and Li6.5defect, we want to compare the results of MLIP vs other commonly utilized models, using a representative composition of Li6. We will compare the model force/virial against DFT values, and lattice parameters, self-diffusivity, and ionic conductivity against experimental values. 4.4.1.1 Force and virial Atomic force (3 for each atom) and cell-virial (6 for each structure) values from the trained MLIP, ID, and DFTB models against DFT values out of 145 structures are shown in Figure 24. Deviation from DFT results is assessed by the root mean squared error (RMSE). RMSE of force and virial for the MLIP model are 0.14 eV/Å and 0.20 GPa, respectively, which are lower than that of Core, DFTB, and ID models. 55 Figure 24 The comparison of (a) atomic forces and (b) cell virials of Core, DFTB, ID, and MLIP (scatter points) against DFT (black lines) for Li6La3ZrTaO12. 4.4.1.2 Lattice parameters of Li6La3ZrTaO12 using different models We first checked the temperature dependence of lattice parameters, shown in Figure 25 to see if these four models could reasonably predict the thermal expansion of Li6La3ZrTaO12 (Li6). We include diffraction results from various works of literature [46, 54, 86, 129] in the figure for comparison. Most experimental measurements of lattice parameters were performed at room temperature. At room temperature, the lattice parameters of MLIP, Core, and ID models have an error between 0.1 - 0.3 %, while the DFTB model’s error is ~0.6 %. If we use the 10 K and 300 K measurement results from Wang et al.[86] as a reference, we see that all four models predict a nearly constant thermal expansion and similar expansion coefficients. Specifically, the thermal expansion coefficients for the MLIP, Core, ID, and DFTB models are 1.74×10-5, 1.51×10-5, 1.57 ×10-5, and 1.91×10-5 K-1, respectively. 56 Figure 25 Lattice parameter as a function of the temperature of Li6La3ZrTaO12 from different models. Literature values of Al-free samples at room temperature[46, 54, 86, 129] are shown for comparison. XRD: X-ray diffraction; ND: neutron diffraction. Solid lines are linear fit. 4.4.1.3 Diffusion and conduction properties comparison among different models Self-diffusivity and ionic conductivity from the four models are plotted in Figure 26, along with experimental results from single crystal measurements. All four models predict similar self- diffusivity and conductivity at high temperatures. However, if we examine the whole temperature, the MLIP model agrees with experiments the best, while the ID and Core models significantly over- and under-estimates, respectively, the low-temperature diffusion and conduction. While it is understandable that the Core model did not perform well due to the very large force error (Figure 24), it is surprising to see the prediction from the ID model, as its force error is comparable to that of MLIP. We believe that this indicates that the ID model lacks transferability. The DFTB model predicts diffusion and conduction similar to those of MLIP, but its force error is much higher. Thus, it is possible some self-cancellation leads to its reasonable performance. In addition, the DFTB model is much more computationally expensive compared with other three models. As shown in Table 7, the simulation time for a 1 ns trajectory is approximately 170 hours for a 1×1×1 supercell, while the other models take less than 40 hours for supercells that are 8 times bigger. 57 Figure 26 (a)Self-diffusivity and (b) ionic conductivity of Li6La3ZrTaO12 as a function of inverse temperature from different models. Single crystal (SC) data from Stanje et al.[129] are shown for comparison. Self-diffusivity from MLIP is slightly higher than experimental values at low temperatures. It is worth noting that the work of Stanje et al.[129] only measured the relaxation time from nuclear magnetic resonance experiments and assumed a jump length of 2 Å. As will be shown later in Figure 28b, we found the jump length is 2.3 Å at 400 K. If 2.3 instead of 2 Å was used, the agreement of MLIP to experiment will be even better. 4.4.2 Force/virial Error and Lattice Parameter Force and virial error against DFT (i.e., PBEsol) for the whole series studied in the present work are shown in Table 8. To understand the meaning of these numbers, we also compared the error between two different exchange-correlation functional, i.e., PBE vs PBEsol, and found the force error was about 0.1 eV/Å. In other words, we can argue the DFT error might be around 0.1 eV/Å and our MLIP model is approaching DFT accuracy but much faster (DFT will be 5 to 6 orders of magnitude slower based on our estimate). Furthermore, it is transferable as it shows similar errors for different compositions, i.e., Li5 to Li7. 58 Table 8 The root mean squared error (RMSE) values of atomic forces, and cell virials between MLIP and DFT for all the compositions studied. RMSE between PBE and PBEsol are shown for comparison. Force (eV/Å) Virial (GPa) Li5 PBE vs Li5 PBEsol 0.10 2.57 MLIP vs Li5 PBEsol 0.15 0.23 MLIP vs Li6 PBEsol 0.14 0.20 MLIP vs Li6.5 PBEsol 0.12 0.21 MLIP vs Li6.75 PBEsol 0.12 0.22 MLIP vs Li7 PBEsol 0.11 0.22 Before applying this MLIP model to study the diffusion and conduction properties, we want to check if it can reproduce the structure of different compositions in terms of lattice parameters. The temperature dependence of lattice parameters is shown in Figure 27a. For cubic phases such as Li5, Li6, Li6.5, and Li6.5defect, lattice parameters increase linearly with increasing temperatures with similar thermal expansion coefficients. The lattice parameter of Li6.5 and Li6.5defect are almost the same. For Li7, the tetragonal to cubic phase transition takes place at a temperature range between 1000 and 1100 K, which is slightly higher than the experimental value of 913 K[53]. For Li6.75, the tetragonal to the cubic phase transition is around 700-800 K. Figure 27b shows calculated lattice parameters at room temperature as a function of x (Li content) in LixLa3Zrx-5Ta7-xO12. Experimental data of single crystals from Kataoka et al.[19], Al- free powder samples from Buschmann et al.[45], Matsuda et al.[46], Awaka et al.[54], and Thompson et al.[130] are shown in the plot for comparison. For cubic phases, the lattice parameter increases linearly with increasing Li content. Our simulation data agree very well with experimental values reported in the literatures. Both simulation and experimental results suggest that the cubic to tetragonal transition takes place at around x=6.6 at room temperature, except for Kataoka et al.[19] where it was cubic even at x=6.8. 59 Figure 27 (a) Lattice Parameters as a function of the temperature of Li5La3Ta2O12 (Li5), Li6La3ZrTaO12 (Li6), Li6.5La3Zr1.5Ta0.5O12 (Li6.5), Li6.75La3Zr1.75Ta0.25O12 (Li6.75), Li7La3Zr2O12 (Li7), and Li5.88La3Zr1.5Ta0.5O11.69 (Li6.5defect). Solid lines are a guide to the eye. (b) Lattice parameters from MLIP-MD at room temperature versus x (Li content). Literature values of single crystals (SC)[19] and Al-free powders[45, 46, 54, 130] are shown for comparison. 4.4.3 Temperature Dependence of Self-diffusion Li diffusion was investigated by the incoherent density correlation function as introduced in Chapter 2.2, from which the residence time, jump length, and self-diffusivity (𝐷𝐿𝑖 ) could be extracted. Figure 28a-c summarizes the residence time, jump length, and self-diffusivity of Li5, Li6, Li6.5, Li6.75, and Li7 at the temperature range of 400 – 1200 K. In terms of self-diffusivity (Figure 28c), temperature dependence is roughly Arrhenius with activation energies (0.26 - 0.33 eV), except for Li7, where the cubic (0.27 eV) and tetragonal (1.33 eV) phase have very different activation energies. Li5 has the highest diffusivity at high temperatures while Li6.5 has the highest diffusivity at low temperatures. Self-diffusivity is related to residence time (𝑡) and jump length (𝑙) as 𝑙2 /(6𝑡). Figure 28b indicates that jump lengths have a narrow range of 1-3.5 Å, with shorter jumps at higher temperatures. The nearest neighbor tetrahedral and octahedral sites in lithium garnets are 1.5-2.5 Å[44]. The plot of residence time, the average time a Li spends on a site before 60 jumping to the next site, largely follows Figure 28c, as the jump length is relatively temperature insensitive. Figure 28 Inverse temperature dependence of (a) residence time, (b) jump length, and (c) self- diffusivity of Li of Li5La3Ta2O12 (Li5), Li6La3ZrTaO12 (Li6), Li6.5La3Zr1.5Ta0.5O12 (Li6.5), Li6.75La3Zr1.75Ta0.25O12 (Li6.75), and Li7La3Zr2O12 (Li7). Solid lines are a guide to the eye and shades in (b) are the range of error. (d-e) Self-diffusivity as a function of inverse temperature for Li5 and Li6.5. Experimental data from pulse-field gradient NMR measurements of single crystal (SC)[21, 131, 132] and QENS measurements [44] are shown for comparison. To show a comparison of self-diffusivities between the results calculated from our MD simulations and that from the experiments, we plotted the results of Li5 and Li6.5 separately in Figure 28(d-e). A comparison for Li6 was shown previously in Figure 26a. Computed self- diffusivity values are slightly lower than those from quasi-elastic neutron scattering (QENS) measurements[44] for Li5, while they are slightly higher than those from pulse-field NMR measurements[21, 131, 132] for single crystals of Li6.5. In Figure 28e, we also compared the normal and defect structure of Li6.5. Self-diffusivities of these two structures are almost the same 61 at high temperatures, e.g., 900 -1200 K, while the defect structure has noticeably lower diffusivity at low temperatures. 4.4.4 Temperature Dependence of Ionic Conductivity Similar to the self-diffusivity plot, we present in Figure 29a the ionic conductivity of Li5, Li6, Li6.5, Li6.75 and Li7, and in Figure 29d-f the comparison to experiments. It is usually assumed an Arrhenius relation between ionic conductivity and reverse temperature from a limited temperature range, however, examination of our wide-temperature experimental measurement of ionic conductivities of Li6 and Li6.5 from this work, along with that from Stanje et al.[129] and Jin et al.[133] as in Figure 29d/e, clearly suggests a non-Arrhenius behavior. Given this, we chose a common non-Arrhenius model, Vogel-Tammann-Fulcher (VTF) equation[134-136] with 𝜎 = 𝐴 ∙ 𝐵 𝑒𝑥𝑝 [− ], where 𝐴 is a pre-exponential factor, 𝐵 the artificial activation energy, 𝑇0 the Vogel 𝑇−𝑇0 temperature below which ions are not mobile. We applied the VTF fit to our MD data, except for Li7 where we applied the Arrhenius fit for the cubic (0.13 eV) and tetragonal (1.15 eV) phase separately. The transition temperature 𝑇0 is plotted in Figure 29b. Its values are below 200 K and decrease with increasing Li content, which suggests that higher Li content might lead to higher conductivity at low temperatures. If we extrapolate MD data through the VTF fit to lower temperatures, we can see they agree well with single crystal data from experiments, except for Li6.75. 62 Figure 29 (a) Ionic conductivity as a function of inverse temperature for Li5La3Ta2O12 (Li5), Li6La3ZrTaO12 (Li6), Li6.5La3Zr1.5Ta0.5O12 (Li6.5), Li6.75La3Zr1.75Ta0.25O12 (Li6.75), and Li7La3Zr2O12 (Li7). Solid lines are fitted to the VTF equation (except for Li7 where two Arrhenius fits are applied to the high and low temperatures). Error bars of Li6.75 and Li7 represent the ionic conductivity in the ab and c directions. (b) Vogel temperature (T0) of different Li content. (c) 𝐻𝑅 for the composition series. Solid lines are a guide to the eye. (d-f) Comparison of results from (a) and experimental values for Li6, Li6.5/Li6.5defect (Li5.88La3Zr1.5Ta0.5O11.69), Li6.75. Experimental values include single crystal (SC)[19, 21, 129, 131], Al-free[51], and Al contaminated powders (through alumina crucibles)[133]. In Figure 29d, we also compared the Nernst-Einstein conductivity (𝜎𝑁𝐸 ) and (coherent) ionic conductivity of Li6. It shows that the coherent ionic conductivity is higher than the Nernst- Einstein conductivity, which indicates the ionic motion is positively cooperated, i.e., ions moving together is faster than moving alone. This positive cooperation was called concerted in the studies of Li7[116, 137]. The Haven ratio is plotted in Figure 29c for all compositions, with values between 0.1 and 0.4. A Haven ratio of ~0.3 was reported by He et al.[137] (900 K) and Mottet et al.[138]. Temperature dependence of Haven ratio indicates a higher value (less cooperation) at higher temperatures. This is intuitive as it is expected that the thermal agitation tends to break the 63 atomic cooperation, which makes the atomic motion more independent. In Figure 29e, we compared the ionic conductivity of Li6.5 and Li6.5defect. Similar to the self-diffusivity comparison in Figure 28e, the ionic conductivity of normal structure is higher. 4.4.5 Finite-size effects on Li-ion transport We further investigated the finite-size effect on self-diffusivity and ionic conductivity of Li6.5La3Zr1.5Ta0.5O12 at a temperature range of 800 -1200 K to examine if the results converged with respect to the supercell size. The initial cells of 1×1×1 (188 atoms) and 3×3×3 (5076 atoms) were adopted to compare with the current 2×2×2 (1504 atoms) structures. The comparison of the finite-size effect is displayed in Figure 30. It is found that 1×1×1 cell has higher values of ionic conductivity at 900 and 1200 K. The difference between 2×2×2 and 3×3×3 supercells at high temperatures is ignorable. Note that, the simulation time of 3×3×3 supercells are more than three times of 2×2×2 supercells, hence, considering the computational cost, we believe that the 2×2×2 supercell is the best option. Figure 30 Finite-size effect on (a) self-diffusivity and (b) ionic conductivity for 1×1×1, 2×2×2, and 3×3×3 cells of Li6.5La3Zr1.5Ta0.5O12. 64 4.4.6 Composition Dependence To examine the composition dependence of diffusion, we plotted in Figure 31a the self- diffusivity at 1200 K and 400 K to represent high and low temperatures, respectively. At 1200 K, the self-diffusivity decreases with increasing Li content, suggesting a vacancy mechanism. At high temperatures, each Li-ion has enough thermal energy to move around. However, as Li ions must move through a network of interconnected tetrahedral and octahedral sites, they will be competing for the vacant sites. At low temperatures, e.g., 400 K, Li ions do not have sufficient thermal energy so only some of them are moving around. Adding more Li will increase the overall Li-Li repulsion and initiate more Li to move, leading to overall increasing diffusivity with increasing composition. The composition dependence of ionic conductivity is shown in Figure 31b. It is worth noting that we could not detect Li conduction at 400 K for Li7 so MD simulation suggests the maximum conductivity at 400 K occurs at Li6.75. If we compare Figure 31b with Figure 21 presented in the motivation for room temperature experimental values from various works of literature, we can see the overall composition dependence looks similar. The difference is the peak composition. For example, results from Buschmann et al.[46] indicate the peak composition is Li6. Yi et al.[51] showed that Li6.7 has the highest ionic conductivity. Inada et al.[48] found Li6.75 to be cubic and its ionic conductivity is lower than that of Li6.5. Matsuda et al.[54] reported the peak at Li6.6 and Li content 6.625 to 7 to be tetragonal. While these reports are all on polycrystalline samples, Kataoka et al.[19] synthesized single crystals and found Li6.2 to Li6.8 to be cubic with a peak conductivity at Li6.6. From the perspective of experiments, single crystals are probably preferred over polycrystals as they are less susceptible to the influence of density and grain boundary segregation (e.g., space charge layer). However, an arbitrary amount of Li excess had to be added to all these samples, which makes the composition used for conductivity measurement a 65 little ambiguous after an arbitrary amount of Li loss during the high-temperature processing. Our MD results suggest that Li loss, e.g., Li6.5defect, will decrease the ionic conductivity. Finally, experimental samples are also subjected to H2O and CO2 contamination to form LiOH and Li2CO3 and it is reasonable to expect such contamination will be more severe when Li is more mobile, i.e., high Li content. Nonetheless, it is probably reasonable to argue the maximum conductivity occurs somewhere between 6.6 and 6.8 based on combined experimental and computational results. Figure 31 Composition dependence of (a) self-diffusivity and (b) ionic conductivity at 1200 and 400 K. 4.5 Conclusions In this chapter, we studied the diffusion and conduction of lithium garnet oxides LixLa3Zrx- 5Ta7-xO12 (x = 5, 6, 6.5, 6.75, and 7) with MD simulations based on a machine learning interatomic potential (MLIP). We found that this MLIP is approaching DFT accuracy in force error. It is more accurate than other commonly applied models such as Core, induced dipole (ID), and DFTB and predicts diffusion and conduction properties that agree with single crystal experimental data much better than the other three models, when using Li6La3ZrTaO12 as a model material. This MLIP also produces self-diffusivity and ionic conductivity values agreeing with experimental data across the 66 x composition range. Examination of computational and experimental ionic conductivities suggest that the temperature dependence is non-Arrhenius which can be fitted to a Vogel-Tammann- Fulcher equation with Vogel temperature below room temperature. Examination of computational and experimental values together also suggests that the maximum room-temperature conductivity occurs between x=6.6 to 6.8 with the precise composition depending on Al/H2O/CO2 contamination, Li loss, grain boundary, density, etc. 67 CHAPTER 5 Grain boundary contributions to Li-ion transport in Li7La3Zr2O12 5.1 Motivation Like other oxide-based solid electrolytes, most of the synthesized lithium garnet samples are generally polycrystalline consisting of a large number of grains separated by grain boundaries (GBs). Although Li garnet oxides have been considered as one of the most promising solid electrolytes due to their chemical stability and high ionic conductivity, their practical application has been hindered by the GBs. After the first report of lithium garnet oxide Li5La3M2O12 (M=Nb, Ta)[39], there has been extensive research to improve its ionic conductivity. Experimentally, the ionic conductivity of Li garnet oxides is commonly measured using electrochemical impedance spectroscopy (EIS). Figure 32 shows a typical Nyquist plot of Li7La3Zr2O12 (named Li7 afterward) at 18 °C[139]. The left depressed semicircle in high frequency represents the bulk resistance (𝑅𝑏 ), the right one represents the GB resistance (𝑅𝑔𝑏 ), and the slope line represents the blocking effect of Au paste on the Li+ motion. The solid line represents the fitted data based on an equivalent circuit model of (𝑅𝑏 𝑄𝑏 )(𝑅𝑔𝑏 𝑄𝑔𝑏 )(𝑄𝑒𝑡 ), where 𝑄𝑏 , 𝑄𝑔𝑏 , and 𝑄𝑒𝑡 are the constant phase element (CPE) contributions of 𝑅𝑏 , 𝑅𝑔𝑏 , and electrode. It is observed that the contribution of GB to the ionic conductivity is about 44 %, which highly reduces the total ionic conductivity of the material. This characteristic has also been found in some other oxide SEs[140-143]. Even though only a few works have been done to study the effect of the GBs on the transport in Li7, there are some conflicting conclusions. Yu et.al. first performed MD simulations of three symmetric tilt GB structures using a Core model with IP based on Morse-type interactions derived from softBV bond valence[144]. Shiiba et al. also performed Core model MD simulations using Buckingham potential with partial charge[145]. In addition to the symmetric GBs, they also 68 studied GBs with mixed GB planes. Both of these two research groups found that GBs retard the Li+ diffusion. The difference is Yu et al. et al believe that higher energy GBs penalize the ionic conductivity to a greater degree whereas Shiiba et al. proposed the opposite conclusion. Recently, Gao et al. studied the GB effects via DFT-MD[146]. In contrast to both Yu and Shiiba et al., Gao et al. found that GB can also increase ionic conductivity. Based on the discrepancy, it is important and meaningful to employ our newly developed MLIP to study the GB contribution in the transport of lithium garnet oxides. Figure 32 Impedance plot of Li7La3Zr2O12 measured in air at 18 °C. The open circles are experimental values and the solid line represents the fitted data based on the equivalent circuit of (RbQb)(RgbQgb)(Qel)[139]. 5.2 Computational details 5.2.1 MLIP potential After comparing our MLIP-MD results with the experiments from the works of literature in Chapter 4, we have found that our newly developed MLIP provides reliable predictions for the properties of lithium garnet oxides, making it a valuable tool for further research, e.g., GB effect in lithium garnet oxides. It should be mentioned that the MLIP used in this chapter was trained from Li7 to study the GBs of Li7, named MLIP-Li7. The force error of MLIP-Li7 for each species 69 and all the atoms and virial error for all the atoms against DFT are shown in Figure 33. The force error of La, Zr, O, and Li is 0.11, 0.17, 0.10, and 0.05 eV/Å, respectively, with an overall error of 0.10 eV/Å. The virial error is 0.19 GPa. The force and virial error are slightly smaller than the values of MLIP (0.11 eV/Å and 0.22 GPa, see Table 8) obtained from the training set of Li5, Li6, and Li7, suggesting the potential is reliable. Figure 33 The comparison of atomic forces and cell virials between MLIP-Li7 (scatter points) against DFT (black lines) for Li7La3Zr2O12. To test the performance of MLIP-Li7, we used it to run MD simulations of Li7 within NPT ensembles and compared the lattice parameters with the results obtained from MLIP presented in Chapter 4, as shown in Figure 34a. There is a negligible difference in the values and the same tetragonal-cubic phase transition temperature, which is around 1000 K. Note that, in Chapter 4, we performed NVE instead of NVT ensemble, therefore we conduct NVT ensemble using both potentials and extracted the ionic conductivity of Li7, as shown in Figure 34b. It can be seen that the values are almost identical, the ionic conductivity decreases dramatically as a tetragonal phase, 70 and the activation energy of MLIP-Li7 is 1.24 and 0.15 for the tetragonal and cubic phases, respectively, which are very close to the values of MLIP. This confirms that MLIP-Li7 is reliable and appropriate to be applied to study the GB effect on the transport of Li7. Figure 34 (a) Lattice parameter and (b) ionic conductivity of Li7La3Zr2O12 extracted from MLIP and MLIP-Li7 based MD trajectories. 5.2.2 GB structure construction The initial GB structures were constructed using the software Materials Studio. Three symmetric tilt grain boundary geometries with low-order coincidence site lattices (CSL), i.e., Σ 1(110), Σ3(112), and Σ5(310), have been studied for Li7. The selection of symmetric tilt GBs studied in this work was based on the low-energy GB orientation in BCC metals [147, 148]. As symmetric GB is a special case in the plane defect, for comparison, we have also studied the mixed GB structure formed by planes (110) and plane (211) (named mixGB(110)(211) afterward). The structural parameters of all the GB structures studied in this chapter are listed in Table 9. Figure 35 illustrates the initial structures of the GB models using relatively smaller structures. As Zr occupies the BCC-type sublattice, each model was provided with a structure that only contains Zr atoms to illustrate the coincident-site nature of the GB model. It can be seen that there are two GB 71 planes in each GB structure, one is in the center of the structure along the a-axis, and the other one is at the a-axis’s boundary. Table 9 Lattice parameters, number of atoms, and number of formula units of each GB structure. Structure model a (Å) b (Å) c (Å) No. of atoms No. of unit cells of Li7 Σ1(110)-4nm 73.95 26.21 18.54 3072 16 Σ3(112)-3nm 64.13 37.07 22.70 4608 24 Σ5(310)-4nm 84.94 26.21 41.44 7680 40 mixGB(110)(211)-3nm 69.70 24.46 37.07 5376 28 Σ3(112)-5nm 96.23 37.07 22.70 6912 36 Σ3(112)-10nm 192.54 37.07 22.70 13824 72 Σ1(110)-9nm 185.16 26.21 37.07 15360 80 Σ5(310)-8nm 165.77 26.21 41.45 15360 80 mixGB(110)(211)-9nm 191.00 24.46 37.07 14592 76 72 Figure 35 Initial structures of three symmetric tilt GBs, i.e., (a, c) Σ1(110), (b, d) Σ3(112), and (e, g) Σ5(310), and one mixed GB structure, i.e., (f, h) mixGB(110)(211). As Zr occupies the BCC- type sublattice in Li7La3Zr2O12, (a, b, e, and f) show the GBs with only Zr atoms. The GB interfaces are guided by the dotted lines. There are two GB planes in each GB structure, one is in the center of the structure along the a-axis and one at the a-axis’s boundary. 73 5.2.3 MD simulations The training data is collected from DFT-MD trajectories of Li7 with energy, force, and cell- virial, therefore the potential will be named MLIP-Li7 afterward to differentiate from the MLIP trained from the reference data set of Li5, Li6, and Li7 in Chapter 3. Details of training can be found in Chapter 4. Geometry optimization and MD simulations were both performed using LAMMPS. MD simulations were conducted in the temperature range of 700 - 1200 K. NPT ensemble simulation was performed to obtain the lattice parameters. The timestep was 1 fs, thermostat and barostat relaxation time were set to 0.1 and 0.5 ps, respectively. The NVT ensemble was performed to extract the transport behavior. The trajectory length is from 100 ps – 20 ns depending on the temperature. 5.3 Results and discussions 5.3.1 Grain boundary energy After the geometry optimization, we performed NPT of GB structures at a temperature range of 700 – 1200 K and calculate the GB formation energy. The evolution of the total energy of each GB structure at different temperatures is shown in Figure 36. The average of the equilibrate states was used to calculate the GB formation energy at this temperature. Figure 37 summarized the GB formation energies as a function of temperature for both symmetric and asymmetric GB structures. It can be seen that the fluctuation of each GB structure at different temperatures is very small. The mixed GB structures exhibit 2 – 3 times higher GB formation energy than the symmetric GB structures. For symmetric GBs, Σ3(112) has the lowest GB energy while Σ5(310) has the highest GB energy, which is consistent with the results from the classic IP model simulations[144, 145]. 74 Figure 36 Total energy evolution of GB structures Σ1(110), Σ3(112), Σ5(310) and mixGB(110)(211) at a temperature range of 700 - 1200 K. Figure 37 GB energy as a function of the temperature for both symmetric and mixed GB structures. The literature values from Yu[144], Shiiba[145], and Gao[146] are shown for comparison. 75 5.3.2 Grain boundary composition In order to compare the properties between the GB and the bulk regions in the GB structures, we performed MD simulations using the larger cell for each GB structure shown in Table 9. Figure 38 compares the nominal composition and the average composition from NVT 800 K for each GB structure. We divided the supercells into several bins along the a-axis (normal to the GB plane) and calculated the number of each species in each bin. For symmetric GBs, the average number of Zr in each bin is still the same, whereas the average composition of La is different in different GBs, i.e., uniform in Σ3(112) but reduced in the GB region of Σ1(110) and Σ5(310). In addition, the number of Li and O reduces in the vicinity of the GB compared to the bulk region for all the GBs, forming an ‘M’, suggesting that there is a reduction of Li and O at the GBs. The reduction of Li population at the GBs suggests that Li deficient sites formed at the GB. This is consistent with the observation from Shiiba et al. [145] who also calculated the population of Li atoms at the GB and believed that the trapping Li vacancies formed at the GB[145]. However, it is different from the calculation from Yu et al.[144] as they found there is Li and O accumulation at the GBs and suggested that the cosegregation of Li and O forms a new phase of Li2O[144]. Gao et al. calculated the Li interstitial formation energy using DFT-MD and also found there is Li accumulation at the GBs[146]. On the other hand, for mixed GB structure mixGB(110)(211), Li and O atoms tend to accumulate in the bulk of grain (211). 76 Figure 38 Atom numbers of Li, La, Zr, and O as a function of position normal to the GB plane (a- axis). The light solid gray lines represent the nominal composition and the colorful dashed lines with symbols represent average structures at 800 K. 5.3.3 Charge With the obtained composition in Figure 38 and the DDEC6 charge in Table 6, we calculated the average charge in each bin for GB structures at 800 K. Both Σ3(112) and Σ5(310) exhibit charge accumulation at the GBs, the width of this space charge area is around 3 nm. However, this characteristic has not been found in Σ1(110) and mixGB(110)(211), where the space charge region spans the whole cell. 77 Figure 39 Average charge as a function of position normal to the GB plane (x-axis) of MD 800 K for GB-containing structures. 5.3.4 Ionic conductivity For ionic conductivities, we first studied the grain size effect by comparing three different grain sizes of the Σ3(112) model. The grain length along the b and c axis (along the GB plane) are the same, whereas it is different along the a-axis (normal to the GB plane), which are 3.21, 4.81, and 9.63 nm named as 3, 5, and 10 nm, respectively (see Table 9). Figure 40 compares the ionic conductivity of these three models with the bulk Li7 in two perpendicular directions separately. Again, crystal Li7 has a tetragonal-cubic phase transition at 1000 K. Above this temperature, the ionic conductivities of the GB structures are lower than that of the bulk, while under this temperature, it is vice versa. Based on the width of GB area (as shown in the density map later) and the ionic conductivity of bulk Li7 from Figure 34, we calculated the local ionic conductivity 78 of Σ3(112) of different grain size at 1000 K, as shown in Figure 40c. It indicates that the grain size has a significant effect on the behavior of GBs with respect to Li-ion transport. The GB structure with the biggest grain size behaves more like the bulk material in both directions, confirming that GB models with small grain sizes overestimate the GB effect of the Li-ion transport in the material. Figure 40 Ionic conductivity as a function of the temperature of different grain sizes of Σ3(112). (a) normal to and (b) along the GB plane. (c) ionic conductivity as a function of grain size. Figure 41a-d shows the ionic conductivity of different GB structures along and normal to the GB plane. The smallest cell size for each GB structure in Table 9 was used to clearly see the GB effect. It can be seen that Σ3(112) exhibits obvious anisotropic characteristics especially when temperature decreases. The ionic conductivity normal to the GB plane is much lower than that along the GB plane, and it is almost four times lower at 700 K. This observation is consistent with the analysis from Yu et al.[144] and Shiiba et al.[145] while completely opposite to the conclusion from Gao et al. where Σ3(112) exhibits a high intergranular diffusion[144]. Σ1(110) exhibits similar behavior with a small difference between these two perpendicular directions. However, Σ5(310) shows an almost isotropic characteristic at high temperatures, and interestingly, the ionic conductivity normal to the GB plane is even slightly lower. For the mixGB(110)(211), it shows an anisotropic behavior with lower ionic conductivity normal to the GB plane. Figure 41e-f compared the ionic conductivity of each GB structure normal to the GB plane and along the GB plane, 79 respectively. First, the mixGB(110)(211) exhibits an obviously lower ionic conductivity at high temperatures in both directions, and it is vice versa at low temperatures. For symmetric GBs, the conductivity normal to the GB plane increases with increasing GB energy, i.e., Σ3(112) to Σ1(110) and to Σ5(310). There is no direct correlation between conductivity along the GB plane with the GB energy, but Σ5(310) shows the smallest value. Figure 41g compared the average ionic conductivity of different GB structures with the bulk. For symmetric GBs, there is a kink between 900 and 1000 K. At the temperature range of 1000 – 1200 K, the ionic conductivities of all the GB structures are smaller than the bulk, suggesting that the ionic conductivity degrades at the GBs in the cubic phase. At the temperature range of 700 – 900 K, the ionic conductivities of all the GB structures are higher than the bulk, suggesting that the existence of GB changes the transport properties in the tetragonal phase. The activation energies of Σ3(112), Σ1(110), and Σ5(310) at high temperatures are 0.26, 0.17, and 0.20 eV, respectively, which are all higher than the bulk (0.15 eV), while activation energies at low temperatures are 0.63, 0.70, and 0.51 eV, respectively, which are lower than the bulk (1.24 eV). The activation energy difference between high and low temperatures is not observed in the literatures, and therefore only one value of activation energy was presented for each GB structure, as listed in Table 10. Interestingly, the ionic conductivity of mixGB(110)(211) shows an Arrhenius behavior with an activation energy of 0.29 eV. Overall, for symmetric GBs in the cubic phase, the higher the GB formation energy, the higher the ionic conductivity and smaller activation energy, and vice versa. This is consistent with the conclusion from Shiiba et al.[145] but opposite to that from Yu et al.[144, 145]. 80 Figure 41 Comparison of ionic conductivity between normal to and along the GB plane for GB structures (a-d) Σ1(110), Σ3(112), Σ5(310), and mixGB(110)(211). (e-f) Comparison of ionic conductivity among different GB structures normal to and along the GB plane, respectively. (g)Comparison of average ionic conductivity of each GB structure with bulk material. Solid lines are Arrhenius fit. Error bars of the bulk structure represent the ionic conductivity in the ab and c directions in the tetragonal phase. Table 10 Activation energy, Ea (eV), of current work and literatures [144-146]. Activation energies of high T are calculated from 1000 – 1200 K for all the structures, and of low T from 700 – 900 K for the GB structures and 800 – 1000 K for the bulk. The GB energy was obtained by averaging over 700 – 1200K. Current work Structures GB energy (J/m2) Yu (Core) Shiiba (Core) Gao (DFT) High T Low T bulk — 0.15 1.24 0.52 0.18 0.25 Σ1(110) 0.74 0.17 0.70 — 0.32 0.49 Σ3(112) 0.63 0.26 0.63 0.52 0.41 0.20 Σ5(310) 1.05 0.20 0.51 0.64 0.31 — mixGB(110)(211) 2.23 0.29 — — — 5.3.5 Nuclear density map Figure 42 and Figure 43 show nuclear density maps of Li-ions with an isosurface level of 0.2 Å-3 viewed from along the GB plane (i.e., along the b-axis) of different GB structures at temperatures of 1000 and 800 K, respectively. These maps provide a visualization of the pathway 81 of atoms by dividing the supercell into voxels with a length of roughly 0.1 Å and calculating the probability density of each type of atom in each voxel. The figures show that the GB structures exhibit different characteristics in the vicinity of the GB compared to the bulk area. The widths of the GB area of Σ1(110), Σ3(112), Σ5(310), and mixGB(110)(211) at 1000 K are approximately 1.61, 1.22, 1.86, and 1.89 nm, respectively, as indicated by the black dash frame. These widths are consistent with the results from Tenheff et al.[149], who synthesized Li7 pellets using the hot- pressed technique and observed a grain boundary width of ~2 nm in pellets with a relative density of 98% and a grain size of 5 μm. According to Figure 43, all the GB structures exhibit a more discontinuous network for Li in both GB and the bulk areas due to the lower thermal energy at 800 K compared to 1000 K. Moreover, the width of the GB area at 800 K is slightly longer than that at 1000 K, which is around 2 nm. The first two GBs, Σ1(110) and Σ3(112), show a clear discontinuity along the a-axis (normal to the GB plane), implying that the ionic conductivity in this direction is low, which could be attributed to the fact that Li-ions require high energy to diffuse across the GB plane. Consequently, the overall ionic conductivity in the GB region is reduced. This behavior is illustrated in Figure 41a-b and Figure 41d. On the other hand, Σ5(310) exhibits a discontinuity along both the a- and c-axes, which explains the relatively small difference in ionic conductivity along and normal to the GB plane, as shown in Figure 41c. Finally, mixGB(110)(211) also displays a discontinuity along the a-axis, which suggests that the ionic conductivity in this direction is low as well, leading to reduced overall ionic conductivity in the GB region, as shown in Figure 41f. 82 Figure 42 Li nuclear density map of (a-d) Σ1(110), Σ3(112), Σ5(310), and mixGB(110)(211) derived from MD trajectories at 1000 K with an isosurface level of 0.2 Å-3. 83 Figure 43 Li nuclear density map of (a-d) Σ1(110), Σ3(112), Σ5(310), and mixGB(110)(211) derived from MD trajectories at 800 K with an isosurface level of 0.2 Å-3. 84 5.3.6 Self-diffusivity Figure 44 compares the self-diffusivity of Li-ions in bulk Li7 and GB-containing structures at 1000 K. The figure uses boxplots to illustrate the range, median, and quartiles of diffusivity values, with the boxes representing the interquartile range (IQR) and the whiskers extending to the minimum and maximum values that fall within 1.5 times the IQR. The black boxes represent the average diffusivity in the a, b, and c directions, while the green and orange boxes represent the diffusivity along and normal to the GB plane, respectively. Outliers beyond the range of 1-99% values are marked with open circles. Figure 44 Boxplots of Li-ion self-diffusivity of bulk Li7La3Zr2O12 and GB-containing structures at 1000 K. The figure shows that the bulk Li7 has the highest average self-diffusivity, which is expected since there are no GBs to hinder the movement of ions. In contrast, the mixGB(110)(211) structure has the lowest average self-diffusivity and a larger range of box and whisker compared to the symmetric GB structures. Among the symmetric GB structures, Σ3(112) shows a longer whisker than the other two, suggesting a wider range of self-diffusivity values. The diffusivity 85 normal to the GB plane has a larger range and longer whisker compared to that along the GB plane, and the average diffusivity normal to the GB plane is generally smaller than that along the GB plane, except for Σ5(310). This observation is consistent with the conclusion drawn from Figure 41, which showed that Σ1(110) and Σ3(112) exhibit anisotropic diffusion, while Σ5(310) exhibits isotropic diffusion. Figure 45 shows the scatter plots of structures extracted from MD trajectories over a time period of 2 ns at a temperature of 1000 K for the fastest and slowest 0.5% Li ions in terms of average diffusivities shown in Figure 44. In the case of symmetric grain boundaries, the slowest- moving Li-ions appear to be confined to the boundary itself, while the fastest-moving ions are located in the bulk area and move around freely. However, in the case of mixed GBs, the behavior of the Li-ions appears to be dramatically different. Specifically, the fastest-moving 0.5% of ions are found exclusively in one type of grain (211), while the slowest-moving 0.5% of ions are found exclusively in another type of grain (110). This is likely to be related to reduced Li content in the (110) grain shown in Figure 38. Figure 45 Scatter plots of structures extracted from the MD trajectories over 2 ns at 1000 K for the fastest and slowest 0.5% Li-ions. 86 5.4 Conclusions In this chapter, the effects of GB on lithium garnet oxide Li7La3Zr2O12 were comprehensively studied. It was found that the existence of GBs significantly affects the composition and Li+ transport behavior in the vicinity of GBs regardless of the grain orientation and GB types. Firstly, the GB effect is strongly associated with the grain size, the GB-contained structures act more like bulk materials as grain size increases. Secondly, it shows that the GB effect depends on the crystal phase of the lithium garnet oxide. At high temperatures the bulk of Li7 is a cubic phase, the structure with the high GB formation energy delivers lower ionic conductivity. Especially, the ionic conductivity is significantly reduced at randomly mixed GBs compared to symmetric GBs. However, at low temperatures the bulk of Li7 is a tetragonal phase, the existence of the GB changed the ordered structure in the bulk area and increased the total ionic conductivity of the structure. 87 CHAPTER 6 Dissertation conclusions and future work 6.1 Dissertation conclusions Lithium garnet oxides are a family of promising candidates for all-solid-state Li-ion batteries. In this thesis, we studied Li-ion self-diffusion and coherent ionic conduction of LixLa3Zrx-5Ta7-xO12 (x=5-7) using MD simulations. We applied several IP models including Core, CS, ID, and MLIP models as well as a DFTB model. The DFTB model has high accuracy but is still too slow to study complex systems with thousands of atoms. Overall, we think the ID model combined with the force-matching approach is the most advanced physical IP model as it includes the polarization through atomic polarizability and takes from first-principles results as the input. On the other hand, MLIP is a pure mathematic model approaching DFT accuracy in force error. Both of these two models can study the lithium garnet oxides and even complex GB systems, in the meanwhile, they are 2 to 3 orders of magnitude faster than DFT. In terms of which composition has the maximum ionic conductivity, there are some different conclusions in experimental works, which could be induced by the contamination of the sample from the alumina crucible, low relative density, grain boundary defects in the polycrystalline, and uncertainty of the composition from arbitrarily amount of extra Li added to compensate for the Li loss in the preparation process. However, there is no such issue in computational simulations, which enables the accuracy of the results. Our MLIP-MD simulations reveal that the maximum room-temperature conductivity occurs between x=6.6 to 6.8. Besides, temperature-dependent ionic conductivity shows the Vogel-Tammann-Fulcher (VTF) equation should be used to fit the data and extract the ionic conductivity at room temperature instead of the Arrhenius fit. 88 As most of the samples are polycrystalline, the effect of grain boundary is significantly important, and in-depth research is imperative. According to the GB models we built, including both symmetric and mixed ones, the Li-ion migration is hindered at the GBs despite the type of GBs at high temperatures for Li7 where the bulk material is a cubic phase. Specifically, the mixed GB with higher GB formation energy penalizes Li-ion conduction more than the symmetric GBs with lower GB formation energies. Our model also reveals that the contribution of GB in Li-ion transport along and normal to the GBs is dependent on the type of grain orientation. Besides, the grain size plays an important role in the conductivity too, the bigger the grain size the higher the ionic conductivity. 6.2 Future work Due to the cubic-tetragonal phase transition of Li7, it is virtually impossible to study the GB influence in Li7 at room temperatures through MD simulations directly. Based on our MD simulations, Li6.5 is a cubic phase and has relatively higher ionic conductivity in lithium garnet oxide families, thus it is essential to study the GB influence on the ionic conductivity of Li6.5 under the guidance of our preliminary work of Li7. In this thesis, we studied the diffusion and conduction in lithium garnet oxides at the temperature range of 400 – 1200 K through MD simulations and extrapolate the ionic conductivity at low temperatures using VTF equation. However, there is another method can be tried to investigate the dynamics of Li-ions at low temperatures directly, i.e., kinetic Monte Carlo, which approximates the transport of Li-ion as a sequence of discrete ‘hops’ between distinct lattice sites[150]. As lithium garnet oxides are electrochemical stable with Li metal which enables the Li metal to be used as a negative electrode, the computational study of the interface between the 89 garnet/Li metal has attracted much attention[151-153]. However, as far as we know, there is no report of the study using MLIP, which makes it worth using our MLIP model to study the interface properties. 90 BIBLIOGRAPHY [1] I. E. Agency. Available: https://www.iea.org/reports/world-energy-outlook-2022 [2] B. Dunn, H. Kamath, and J. M. Tarascon "Electrical Energy Storage for the Grid: A Battery of Choices." Science 2011, 334, 928-935. http://dx.doi.org/10.1126/science.1212741 [3] G. L. Soloveichik, "Battery Technologies for Large-Scale Stationary Energy Storage," in Annual Review of Chemical and Biomolecular Engineering, Vol 2. vol. 2, J. M. Prausnitz, Ed., ed Palo Alto: Annual Reviews, 2011, pp. 503-527. [4] W. E. Council. Available: www.worldenergy.org [5] A. Afir, S. M. H. Rahman, A. T. Azad, J. Zaini, M. A. Islan, and A. Azad "Advanced materials and technologies for hybrid supercapacitors for energy storage - A review." Journal of Energy Storage 2019, 25, 24. http://dx.doi.org/10.1016/j.est.2019.100852 [6] J. M. Tarascon and M. Armand "Issues and challenges facing rechargeable lithium batteries." Nature 2001, 414, 359-367. http://dx.doi.org/10.1038/35104644 [7] M. Armand and J. M. Tarascon "Building better batteries." Nature 2008, 451, 652-657. http://dx.doi.org/10.1038/451652a [8] J. B. Goodenough and Y. Kim "Challenges for Rechargeable Li Batteries." Chemistry of Materials 2010, 22, 587-603. http://dx.doi.org/10.1021/cm901452z [9] B. Scrosati and J. Garche "Lithium batteries: Status, prospects and future." Journal of Power Sources 2010, 195, 2419-2430. http://dx.doi.org/10.1016/j.jpowsour.2009.11.048 [10] V. Etacheri, R. Marom, R. Elazari, G. Salitra, and D. Aurbach "Challenges in the development of advanced Li-ion batteries: a review." Energy & Environmental Science 2011, 4, 3243-3262. http://dx.doi.org/10.1039/c1ee01598b [11] S. Ramakumar, C. Deviannapoorani, L. Dhivya, L. S. Shankar, and R. Murugan "Lithium garnets: Synthesis, structure, Li+ conductivity, Li+ dynamics and applications." Progress in Materials Science 2017, 88, 325-411. http://dx.doi.org/10.1016/j.pmatsci.2017.04.007 [12] N. Kamaya, K. Homma, Y. Yamakawa, M. Hirayama, R. Kanno, M. Yonemura, et al. "A lithium superionic conductor." Nature Materials 2011, 10, 682-686. http://dx.doi.org/10.1038/nmat3066 [13] J. T. Kummer "Beta-alumina Electrolytes." Journal of the Electrochemical Society 1972, 119, 141-175. [14] H. Y. P. Hong "Crystal Structure and Ionic Coductivity of Li14Zn(GeO4)4 and Other New Li+ Superionic Conductors." Materials Research Bulletin 1978, 13, 117-124. http://dx.doi.org/10.1016/0025-5408(78)90075-2 91 [15] Y. Inaguma, L. Q. Chen, M. Itoh, T. Nakamura, T. Uchida, H. Ikuta, et al. "High Ionic Conductivity in LIthium Lanthanum Titanate." Solid State Communications 1993, 86, 689- 693. http://dx.doi.org/10.1016/0038-1098(93)90841-a [16] X. H. Yu, J. B. Bates, G. E. Jellison, and F. X. Hart "A stable thin-film lithium electrolyte: Lithium phosphorus oxynitride." Journal of the Electrochemical Society 1997, 144, 524- 532. http://dx.doi.org/10.1149/1.1837443 [17] P. E. Stallworth, J. J. Fontanella, M. C. Wintersgill, C. D. Scheidler, J. J. Immel, S. G. Greenbaum, et al. "NMR, DSC and High Pressure Electrical Conductivity Studies on Liquid and Hybrid Electrolytes." Journal of Power Sources 1999, 81, 739-747. http://dx.doi.org/10.1016/s0378-7753(99)00144-5 [18] U. V. Alpen, A. Rabenau, and G. H. Talat "Ionic Conductivity in Li3N Single Crystals." Applied Physics Letters 1977, 30, 621-623. http://dx.doi.org/10.1063/1.89283 [19] K. Kataoka and J. Akimoto "Lithium-ion conductivity and crystal structure of garnet-type solid electrolyte Li7-xLa3Zr2-xTaxO12 using single-crystal." Journal of the Ceramic Society of Japan 2019, 127, 521-526. http://dx.doi.org/10.2109/jcersj2.19022 [20] Y. Z. Zhu, X. F. He, and Y. F. Mo "Origin of Outstanding Stability in the Lithium Solid Electrolyte Materials: Insights from Thermodynamic Analyses Based on First-Principles Calculations." Acs Applied Materials & Interfaces 2015, 7, 23685-23693. http://dx.doi.org/10.1021/acsami.5b07517 [21] A. Dorai, N. Kuwata, R. Takekawa, J. Kawamura, K. Kataoka, and J. Akimoto "Diffusion coefficient of lithium ions in garnet-type Li6.5La3Zr1.5Ta0.5O12 single crystal probed by 7Li pulsed field gradient-NMR spectroscopy." Solid State Ionics 2018, 327, 18-26. http://dx.doi.org/10.1016/j.ssi.2018.10.016 [22] Y. X. Wang and W. Lai "High Ionic Conductivity Lithium Garnet Oxides of Li7-xLa3Zr2- xTaxO12 Compositions." Electrochemical and Solid State Letters 2012, 15, A68-A71. http://dx.doi.org/10.1149/2.024205esl [23] Y. X. Wang and W. Lai "Phase transition in lithium garnet oxide ionic conductors Li7La3Zr2O12: The role of Ta substitution and H2O/CO2 exposure." Journal of Power Sources 2015, 275, 612-620. http://dx.doi.org/10.1016/j.jpowsour.2014.11.062 [24] W. Lai "Transport in Lithium Garnet Oxides as Revealed by Atomistic Simulations." Annual Review of Materials Research 2022, 52, 305-330. http://dx.doi.org/10.1146/annurev-matsci-081720-115334 [25] L. J. Miara, S. P. Ong, Y. F. Mo, W. D. Richards, Y. Park, J. M. Lee, et al. "Effect of Rb and Ta Doping on the Ionic Conductivity and Stability of the Garnet Li7+2x-y(La3-xRbx)(Zr2- yTay)O12 (0 <= x <= 0.375, 0 <= y <= 1) Superionic Conductor: A First Principles Investigation." Chemistry of Materials 2013, 25, 3048-3055. http://dx.doi.org/10.1021/cm401232r 92 [26] J. L. Allen, J. Wolfenstine, E. Rangasamy, and J. Sakamoto "Effect of substitution (Ta, Al, Ga) on the conductivity of Li7La3Zr2O12." Journal of Power Sources 2012, 206, 315-319. http://dx.doi.org/10.1016/j.jpowsour.2012.01.131 [27] J. Dai, Y. Jiang, and W. Lai "Study of diffusion and conduction in lithium garnet oxides LixLa3Zrx-5Ta7-xO12 by machine learning interatomic potentials." Phys Chem Chem Phys 2022, 24, 15025-15033. http://dx.doi.org/10.1039/d2cp00591c [28] S. Grimme, C. Bannwarth, and P. Shushkov "A Robust and Accurate Tight-Binding Quantum Chemical Method for Structures, Vibrational Frequencies, and Noncovalent Interactions of Large Molecular Systems Parametrized for All spd-Block Elements (Z=1- 86)." Journal of Chemical Theory and Computation 2017, 13, 1989-2009. http://dx.doi.org/10.1021/acs.jctc.7b00118 [29] C. W. Wang, K. Fu, S. P. Kammampata, D. W. McOwen, A. J. Samson, L. Zhang, et al. "Garnet-Type Solid-State Electrolytes: Materials, Interfaces, and Batteries." Chemical Reviews 2020, 120, 4257-4300. http://dx.doi.org/10.1021/acs.chemrev.9b00427 [30] B. R. Tao, C. J. Ren, H. D. Li, B. S. Liu, X. B. Jia, X. W. Dong, et al. "Thio-/LISICON and LGPS-Type Solid Electrolytes for All-Solid-State Lithium-Ion Batteries." Advanced Functional Materials 2022, 32, 25. http://dx.doi.org/10.1002/adfm.202203551 [31] K. Takada "Progress and prospective of solid-state lithium batteries." Acta Materialia 2013, 61, 759-770. http://dx.doi.org/10.1016/j.actamat.2012.10.034 [32] J. C. Bachman, S. Muy, A. Grimaud, H. H. Chang, N. Pour, S. F. Lux, et al. "Inorganic Solid-State Electrolytes for Lithium Batteries: Mechanisms and Properties Governing Ion Conduction." Chemical Reviews 2016, 116, 140-162. http://dx.doi.org/10.1021/acs.chemrev.5b00563 [33] P. Knauth "Inorganic solid Li ion conductors: An overview." Solid State Ionics 2009, 180, 911-916. http://dx.doi.org/10.1016/j.ssi.2009.03.022 [34] Y. F. Y. Yao and J. T. Kummer "ION EXCHANGE PROPERTIES OF AND RATES OF IONIC DIFFUSION IN BETA-ALUMINA." Journal of Inorganic & Nuclear Chemistry 1967, 29, 2453-&. [35] A. Rabenau and H. Schulz "Re-evaluation of lithium nitrides structure." Journal of the Less-Common Metals 1976, 50, 155-159. http://dx.doi.org/10.1016/0022-5088(76)90263- 0 [36] A. Rabenau "Lithium nitride and related materials - case sudy of the use of modern solid- state research techniques." Solid State Ionics 1982, 6, 277-293. http://dx.doi.org/10.1016/0167-2738(82)90012-1 [37] R. Kanno, T. Hata, Y. Kawamoto, and M. Irie "Synthesis of a new lithium ionic conductor, thio-LISICON-lithium germanium sulfide system." Solid State Ionics 2000, 130, 97-104. http://dx.doi.org/10.1016/s0167-2738(00)00277-0 93 [38] J. B. Bates, N. J. Dudney, G. R. Gruzalski, R. A. Zuhr, A. Choudhury, C. F. Luck, et al. "Electrical properties of amorphous lithium electrolyte thin-films." Solid State Ionics 1992, 53, 647-654. http://dx.doi.org/10.1016/0167-2738(92)90442-r [39] V. Thangadurai, H. Kaack, and W. J. F. Weppner "Novel fast lithium ion conduction in garnet-type Li5La3M2O12 (M = Nb, Ta)." Journal of the American Ceramic Society 2003, 86, 437-440. http://dx.doi.org/10.1111/j.1151-2916.2003.tb03318.x [40] V. Thangadurai and W. Weppner "Li(6)ALa(2)Ta(2)O(12) (A=Sr, Ba): Novel garnet-like oxides for fast lithium ion conduction." Advanced Functional Materials 2005, 15, 107-112. http://dx.doi.org/10.1002/adfm.200400044 [41] V. Thangadurai, S. Narayanan, and D. Pinzaru "Garnet-type solid-state fast Li ion conductors for Li batteries: critical review." Chemical Society Reviews 2014, 43, 4714- 4727. http://dx.doi.org/10.1039/c4cs00020j [42] M. Xu, M. S. Park, J. M. Lee, T. Y. Kim, Y. S. Park, and E. Ma "Mechanisms of Li+ transport in garnet-type cubic Li3+xLa3M2O12 (M = Te, Nb, Zr)." Physical Review B 2012, 85, 5. http://dx.doi.org/10.1103/PhysRevB.85.052301 [43] E. J. Cussen "The structure of lithium garnets: cation disorder and clustering in a new family of fast Li+ conductors." Chemical Communications 2006, 412-413. http://dx.doi.org/10.1039/b514640b [44] M. J. Klenk, S. E. Boeberitz, J. Dai, N. H. Jalarvo, V. K. Peterson, and W. Lai "Lithium self-diffusion in a model lithium garnet oxide Li5La3Ta2O12: A combined quasi-elastic neutron scattering and molecular dynamics study." Solid State Ionics 2017, 312, 1-7. http://dx.doi.org/10.1016/j.ssi.2017.09.022 [45] J. Awaka, N. Kijima, H. Hayakawa, and J. Akimoto "Synthesis and structure analysis of tetragonal Li7La3Zr2O12 with the garnet-related type structure." Journal of Solid State Chemistry 2009, 182, 2046-2052. http://dx.doi.org/10.1016/j.jssc.2009.05.020 [46] H. Buschmann, S. Berendts, B. Mogwitz, and J. Janek "Lithium metal electrode kinetics and ionic conductivity of the solid lithium ion conductors "Li7La3Zr2O12" and Li7- xLa3Zr2-xTaxO12 with garnet-type structure." Journal of Power Sources 2012, 206, 236- 244. http://dx.doi.org/10.1016/j.jpowsour.2012.01.094 [47] A. Logeat, T. Koohler, U. Eisele, B. Stiaszny, A. Harzer, M. Tovar, et al. "From order to disorder: The structure of lithium-conducting garnets Li7-xLa3TaxZr2-xO12 (x=0-2)." Solid State Ionics 2012, 206, 33-38. http://dx.doi.org/10.1016/j.ssi.2011.10.023 [48] R. Inada, K. Kusakabe, T. Tanaka, S. Kudo, and Y. Sakurai "Synthesis and properties of Al-free Li7-xLa3Zr2-xTaxO12 garnet related oxides." Solid State Ionics 2014, 262, 568- 572. http://dx.doi.org/10.1016/j.ssi.2013.09.008 [49] N. Hamao, K. Kataoka, N. Kijima, and J. Akimoto "Synthesis, crystal structure and conductive properties of garnet-type lithium ion conductor Al-free Li7-xLa3Zr2-xTaxO12 94 (0 <= x <= 0.6)." Journal of the Ceramic Society of Japan 2016, 124, 678-683. http://dx.doi.org/10.2109/jcersj2.16019 [50] Y. Gong, Z. G. Liu, Y. J. Jin, J. H. Ouyang, L. Chen, and Y. J. Wang "Effect of sintering process on the microstructure and ionic conductivity of Li7-xLa3Zr2-xTaxO12 ceramics." Ceramics International 2019, 45, 18439-18444. http://dx.doi.org/10.1016/j.ceramint.2019.06.061 [51] M. Y. Yi, T. Liu, X. N. Wang, J. Y. Li, C. Wang, and Y. C. Mo "High densification and Li-ion conductivity of Al-free Li7-xLa3Zr2-xTaxO12 garnet solid electrolyte prepared by using ultrafine powders." Ceramics International 2019, 45, 786-792. http://dx.doi.org/10.1016/j.ceramint.2018.09.245 [52] Y. T. Li, J. T. Han, C. A. Wang, H. Xie, and J. B. Goodenough "Optimizing Li+ conductivity in a garnet framework." Journal of Materials Chemistry 2012, 22, 15357- 15361. http://dx.doi.org/10.1039/c2jm31413d [53] Y. Matsuda, K. Sakamoto, M. Matsui, O. Yamamoto, Y. Takeda, and N. Imanishi "Phase formation of a garnet-type lithium-ion conductor Li7-3xAlxLa3Zr2O12." Solid State Ionics 2015, 277, 23-29. http://dx.doi.org/10.1016/j.ssi.2015.04.011 [54] Y. Matsuda, Y. Itami, K. Hayamizu, T. Ishigaki, M. Matsui, Y. Takeda, et al. "Phase relation, structure and ionic conductivity of Li7-x-3yAlyLa3Zr2-xTaxO12." Rsc Advances 2016, 6, 78210-78218. http://dx.doi.org/10.1039/c6ra13317g [55] M. Botros, R. Djenadic, O. Clemens, M. Moller, and H. Hahn "Field assisted sintering of fine-grained Li7-3xLa3Zr2AlxO12 solid electrolyte and the influence of the microstructure on the electrochemical performance." Journal of Power Sources 2016, 309, 108-115. http://dx.doi.org/10.1016/j.jpowsour.2016.01.086 [56] C. Bernuy-Lopez, W. Manalastas, J. M. L. del Amo, A. Aguadero, F. Aguesse, and J. A. Kilner "Atmosphere Controlled Processing of Ga-Substituted Garnets for High Li-Ion Conductivity Ceramics." Chemistry of Materials 2014, 26, 3610-3617. http://dx.doi.org/10.1021/cm5008069 [57] R. Wagner, G. J. Redhammer, D. Rettenwander, A. Senyshyn, W. Schmidt, M. Wilkening, et al. "Crystal Structure of Garnet-Related Li-Ion Conductor Li7-3xGaxLa3Zr2O12: Fast Li-Ion Conduction Caused by a Different Cubic Modification?", Chemistry of Materials 2016, 28, 1861-1871. http://dx.doi.org/10.1021/acs.chemmater.6b00038 [58] J. F. Wu, W. K. Pang, V. K. Peterson, L. Wei, and X. Guo "Garnet-Type Fast Li-Ion Conductors with High Ionic Conductivities for All-Solid-State Batteries." Acs Applied Materials & Interfaces 2017, 9, 12461-12468. http://dx.doi.org/10.1021/acsami.7b00614 [59] S. Ohta, T. Kobayashi, and T. Asaoka "High lithium ionic conductivity in the garnet-type oxide Li7-XLa3(Zr2-X, NbX)O12 (X=0-2)." Journal of Power Sources 2011, 196, 3342-3345. http://dx.doi.org/10.1016/j.jpowsour.2010.11.089 95 [60] D. W. Wang, G. M. Zhong, W. K. Pang, Z. P. Guo, Y. X. Li, M. J. McDonald, et al. "Toward Understanding the Lithium Transport Mechanism in Garnet-type Solid Electrolytes: Li+ Ion Exchanges and Their Mobility at Octahedral/Tetrahedral Sites." Chemistry of Materials 2015, 27, 6650-6659. http://dx.doi.org/10.1021/acs.chemmater.5b02429 [61] L. Dhivya, N. Janani, B. Palanivel, and R. Murugan "Li+ transport properties of W substituted Li7La3Zr2O12 cubic lithium garnets." Aip Advances 2013, 3, 21. http://dx.doi.org/10.1063/1.4818971 [62] X. T. Liu, Y. Li, T. T. Yang, Z. Z. Cao, W. Y. He, Y. F. Gao, et al. "High lithium ionic conductivity in the garnet-type oxide Li7-2xLa3Zr2-xMoxO12 (x=0-0.3) ceramics by sol- gel method." Journal of the American Ceramic Society 2017, 100, 1527-1533. http://dx.doi.org/10.1111/jace.14736 [63] P. Bottke, D. Rettenwander, W. Schmidt, G. Amthauer, and M. Wilkening "Ion Dynamics in Solid Electrolytes: NMR Reveals the Elementary Steps of Li+ Hopping in the Garnet Li6.5La3Zr1.75Mo0.25O12." Chemistry of Materials 2015, 27, 6571-6582. http://dx.doi.org/10.1021/acs.chemmater.5b02231 [64] P. Hohenberg and W. Kohn "Inhomogeneous Electron Gas." Physical Review B 1964, 136, B864-B871. http://dx.doi.org/10.1103/PhysRev.136.B864 [65] W. Kohn and L. J. Sham "SELF-CONSISTENT EQUATIONS INCLUDING EXCHANGE AND CORRELATION EFFECTS." Physical Review 1965, 140, 1133-&. http://dx.doi.org/10.1103/PhysRev.140.A1133 [66] J. P. Perdew, K. Burke, and M. Ernzerhof "Generalized gradient approximation made simple." Physical Review Letters 1996, 77, 3865-3868. http://dx.doi.org/10.1103/PhysRevLett.77.3865 [67] U. von Barth "Basic density-functional theory - an overview." Physica Scripta 2004, T109, 9-39. http://dx.doi.org/10.1238/Physica.Topical.109a00009 [68] J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, et al. "Restoring the density-gradient expansion for exchange in solids and surfaces." Physical Review Letters 2008, 100, 4. http://dx.doi.org/10.1103/PhysRevLett.100.136406 [69] M. J. Klenk and W. Lai "Effect of exchange-correlation functionals on the density functional theory simulation of phase transformation of fast-ion conductors: A case study in the Li garnet oxide Li7La3Zr2O12." Computational Materials Science 2017, 134, 132- 136. http://dx.doi.org/10.1016/j.commatsci.2017.03.039 [70] M. Elstner, D. Porezag, G. Jungnickel, J. Elsner, M. Haugk, T. Frauenheim, et al. "Self- consistent-charge density-functional tight-binding method for simulations of complex materials properties." Physical Review B 1998, 58, 7260-7268. http://dx.doi.org/10.1103/PhysRevB.58.7260 96 [71] T. Frauenheim, G. Seifert, M. Elstner, T. Niehaus, C. Kohler, M. Amkreutz, et al. "Atomistic simulations of complex materials: ground-state and excited-state properties." Journal of Physics-Condensed Matter 2002, 14, 3015-3047. http://dx.doi.org/10.1088/0953-8984/14/11/313 [72] P. Koskinen and V. Makinen "Density-functional tight-binding for beginners." Computational Materials Science 2009, 47, 237-253. http://dx.doi.org/10.1016/j.commatsci.2009.07.013 [73] D. Porezag, T. Frauenheim, T. Kohler, G. Seifert, and R. Kaschner "Construction of tight- binding-like potentials on the basis of density-functional theory: Application to carbon." Physical Review B 1995, 51, 12947-12957. http://dx.doi.org/10.1103/PhysRevB.51.12947 [74] P. Cieplak, F. Y. Dupradeau, Y. Duan, and J. M. Wang "Polarization effects in molecular mechanical force fields." Journal of Physics-Condensed Matter 2009, 21, 21. http://dx.doi.org/10.1088/0953-8984/21/33/333102 [75] S. W. Rick, S. J. Stuart, and B. J. Berne "Dynamical fluctuating charge force fields: Application to liquid water." Journal of Chemical Physics 1994, 101, 6141-6156. http://dx.doi.org/10.1063/1.468398 [76] J. Dai, Q. Chen, T. Glossmann, and W. Lai "Comparison of interatomic potential models on the molecular dynamics simulation of fast-ion conductors: A case study of a Li garnet oxide Li7La3Zr2O12." Computational Materials Science 2019, 162, 333-339. http://dx.doi.org/10.1016/j.commatsci.2019.02.044 [77] S. Tazi, J. J. Molina, B. Rotenberg, P. Turq, R. Vuilleumier, and M. Salanne "A transferable ab initio based force field for aqueous ions." Journal of Chemical Physics 2012, 136, 12. http://dx.doi.org/10.1063/1.3692965 [78] K. T. Tang and J. P. Toennies "An improved simple model for the van der Waals potential based on universal damping functions for the dispersion coefficients." Journal of Chemical Physics 1984, 80, 3726-3741. http://dx.doi.org/10.1063/1.447150 [79] J. Behler "Atom-centered symmetry functions for constructing high-dimensional neural network potentials." Journal of Chemical Physics 2011, 134, 13. http://dx.doi.org/10.1063/1.3553717 [80] M. Rupp, A. Tkatchenko, K. R. Muller, and O. A. von Lilienfeld "Fast and Accurate Modeling of Molecular Atomization Energies with Machine Learning." Physical Review Letters 2012, 108, 5. http://dx.doi.org/10.1103/PhysRevLett.108.058301 [81] F. Faber, A. Lindmaa, O. A. von Lilienfeld, and R. Armiento "Crystal structure representations for machine learning models of formation energies." International Journal of Quantum Chemistry 2015, 115, 1094-1101. http://dx.doi.org/10.1002/qua.24917 [82] A. P. Bartok, R. Kondor, and G. Csanyi "On representing chemical environments." Physical Review B 2013, 87, 16. http://dx.doi.org/10.1103/PhysRevB.87.184115 97 [83] M. Hellström and J. Behler, "Neural Network Potentials in Materials Modeling," in Handbook of Materials Modeling: Methods: Theory and Modeling, W. Andreoni and S. Yip, Eds., ed Cham: Springer International Publishing, 2020, pp. 661-680. [84] N. Artrith and A. Urban "An implementation of artificial neural-network potentials for atomistic materials simulations: Performance for TiO2." Computational Materials Science 2016, 114, 135-150. http://dx.doi.org/10.1016/j.commatsci.2015.11.047 [85] M. J. Klenk and W. Lai "Finite-size effects on the molecular dynamics simulation of fast- ion conductors: A case study of lithium garnet oxide Li7La3Zr2O12." Solid State Ionics 2016, 289, 143-149. http://dx.doi.org/10.1016/j.ssi.2016.03.002 [86] Y. X. Wang, A. Huq, and W. Lai "Insight into lithium distribution in lithium-stuffed garnet oxides through neutron diffraction and atomistic simulation: Li7-xLa3Zr2-xTaxO12 (x=0- 2) series." Solid State Ionics 2014, 255, 39-49. http://dx.doi.org/10.1016/j.ssi.2013.11.017 [87] F. Ercolessi and J. B. Adams "Interatomic Potentials from First-Principles Calculations: The Force-Matching Method." Europhysics Letters 1994, 26, 583-588. http://dx.doi.org/10.1209/0295-5075/26/8/005 [88] A. Carre, S. Ispas, J. Horbach, and W. Kob "Developing empirical potentials from ab initio simulations: The case of amorphous silica." Computational Materials Science 2016, 124, 323-334. http://dx.doi.org/10.1016/j.commatsci.2016.07.041 [89] N. G. Limas and T. A. Manz "Introducing DDEC6 atomic population analysis: part 2. Computed results for a wide range of periodic and nonperiodic materials." Rsc Advances 2016, 6, 45727-45747. http://dx.doi.org/10.1039/c6ra05507a [90] T. A. Manz "Introducing DDEC6 atomic population analysis: part 3. Comprehensive method to compute bond orders." Rsc Advances 2017, 7, 45552-45581. http://dx.doi.org/10.1039/c7ra07400j [91] T. A. Manz and N. G. Limas "Introducing DDEC6 atomic population analysis: part 1. Charge partitioning theory and methodology." Rsc Advances 2016, 6, 47771-47801. http://dx.doi.org/10.1039/c6ra04656h [92] J.-P. Hansen and I. R. McDonald, Theory of simple liquids: with applications to soft matter, 4 ed.: Academic Press, 2013. [93] G. Williams and D. C. Watts "NON-SYMMETRICAL DIELECTRIC RELAXATION BEHAVIOUR ARISING FROM A SIMPLE EMPIRICAL DECAY FUNCTION." Transactions of the Faraday Society 1970, 66, 80-+. http://dx.doi.org/10.1039/tf9706600080 [94] K. Trachenko and V. V. Brazhkin "Collective modes and thermodynamics of the liquid state." Reports on Progress in Physics 2016, 79, 36. http://dx.doi.org/10.1088/0034- 4885/79/1/016502 98 [95] K. L. Ngai "Short-time and long-time relaxation dynamics of glass-forming substances: a coupling model perspective." Journal of Physics-Condensed Matter 2000, 12, 6437-6451. http://dx.doi.org/10.1088/0953-8984/12/29/316 [96] H. Jobic and D. N. Theodorou "Quasi-elastic neutron scattering and molecular dynamics simulation as complementary techniques for studying diffusion in zeolites." Microporous and Mesoporous Materials 2007, 102, 21-50. http://dx.doi.org/10.1016/j.micromeso.2006.12.034 [97] K. S. Singwi and A. Sjolander "DIFFUSIVE MOTIONS IN WATER AND COLD NEUTRON SCATTERING." Physical Review 1960, 119, 863-871. http://dx.doi.org/10.1103/PhysRev.119.863 [98] E. Helfand "TRANSPORT COEFFICIENTS FROM DISSIPATION IN A CANONICAL ENSEMBLE." Physical Review 1960, 119, 1-9. http://dx.doi.org/10.1103/PhysRev.119.1 [99] A. S. Nowick and B. S. Berry, Anelastic Relaxation in Crystalline Solids: Academic Press, 1972. [100] Z. Deng, Z. B. Wang, I. H. Chu, J. Luo, and S. P. Ong "Elastic Properties of Alkali Superionic Conductor Electrolytes from First Principles Calculations." Journal of the Electrochemical Society 2016, 163, A67-A74. http://dx.doi.org/10.1149/2.0061602jes [101] S. Adams and R. P. Rao "Ion transport and phase transition in Li7-xLa3(Zr2-xMx)O12 (M = Ta5+, Nb5+, x=0, 0.25)." Journal of Materials Chemistry 2012, 22, 1426-1434. http://dx.doi.org/10.1039/c1jm14588f [102] C. Chen, Z. H. Lu, and F. Ciucci "Data mining of molecular dynamics data reveals Li diffusion characteristics in garnet Li7La3Zr2O12." Scientific Reports 2017, 7, 8. http://dx.doi.org/10.1038/srep40769 [103] M. Klenk and W. Lai "Local structure and dynamics of lithium garnet ionic conductors: tetragonal and cubic Li7La3Zr2O7." Physical Chemistry Chemical Physics 2015, 17, 8758-8768. http://dx.doi.org/10.1039/c4cp05690f [104] M. Burbano, D. Carlier, F. Boucher, B. J. Morgan, and M. Salanne "Sparse Cyclic Excitations Explain the Low Ionic Conductivity of Stoichiometric Li7La3Zr2O12." Physical Review Letters 2016, 116, 6. http://dx.doi.org/10.1103/PhysRevLett.116.135901 [105] P. J. Mitchell and D. Fincham "Shell model simulations by adiabatic dynamics." Journal of Physics-Condensed Matter 1993, 5, 1031-1038. http://dx.doi.org/10.1088/0953- 8984/5/8/006 [106] J. VandeVondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing, and J. Hutter "QUICKSTEP: Fast and accurate density functional calculations using a mixed Gaussian and plane waves approach." Computer Physics Communications 2005, 167, 103-128. http://dx.doi.org/10.1016/j.cpc.2004.12.014 99 [107] J. Hutter, M. Iannuzzi, F. Schiffmann, and J. VandeVondele "CP2K: atomistic simulations of condensed matter systems." Wiley Interdisciplinary Reviews-Computational Molecular Science 2014, 4, 15-25. http://dx.doi.org/10.1002/wcms.1159 [108] G. I. Csonka, J. P. Perdew, A. Ruzsinszky, P. H. T. Philipsen, S. Lebegue, J. Paier, et al. "Assessing the performance of recent density functionals for bulk solids." Physical Review B 2009, 79, 14. http://dx.doi.org/10.1103/PhysRevB.79.155107 [109] G. Kresse and D. Joubert "From ultrasoft pseudopotentials to the projector augmented- wave method." Physical Review B 1999, 59, 1758-1775. http://dx.doi.org/10.1103/PhysRevB.59.1758 [110] G. Kresse and J. Furthmuller "Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set." Physical Review B 1996, 54, 11169-11186. http://dx.doi.org/10.1103/PhysRevB.54.11169 [111] W. Smith and T. R. Forester "DL_POLY_2.0: A general-purpose parallel molecular dynamics simulation package." Journal of Molecular Graphics 1996, 14, 136-141. http://dx.doi.org/10.1016/s0263-7855(96)00043-4 [112] C. A. Geiger, E. Alekseev, B. Lazic, M. Fisch, T. Armbruster, R. Langner, et al. "Crystal Chemistry and Stability of "Li7La3Zr2O12" Garnet: A Fast Lithium-Ion Conductor." Inorganic Chemistry 2011, 50, 1089-1097. http://dx.doi.org/10.1021/ic101914e [113] M. Kotobuki, K. Kanamura, Y. Sato, and T. Yoshida "Fabrication of all-solid-state lithium battery with lithium metal anode using Al2O3-added Li7La3Zr2O12 solid electrolyte." Journal of Power Sources 2011, 196, 7750-7754. http://dx.doi.org/10.1016/j.jpowsour.2011.04.047 [114] H. Buschmann, J. Dolle, S. Berendts, A. Kuhn, P. Bottke, M. Wilkening, et al. "Structure and dynamics of the fast lithium ion conductor "Li7La3Zr2O12"." Physical Chemistry Chemical Physics 2011, 13, 19378-19392. http://dx.doi.org/10.1039/c1cp22108f [115] G. Larraz, A. Orera, and M. L. Sanjuan "Cubic phases of garnet-type Li7La3Zr2O12: the role of hydration." Journal of Materials Chemistry A 2013, 1, 11419-11428. http://dx.doi.org/10.1039/c3ta11996c [116] R. Jalem, Y. Yamamoto, H. Shiiba, M. Nakayama, H. Munakata, T. Kasuga, et al. "Concerted Migration Mechanism in the Li Ion Dynamics of Garnet-Type Li7La3Zr2O12." Chemistry of Materials 2013, 25, 425-430. http://dx.doi.org/10.1021/cm303542x [117] A. Kuhn, S. Narayanan, L. Spencer, G. Goward, V. Thangadurai, and M. Wilkening "Li self-diffusion in garnet-type Li7La3Zr2O12 as probed directly by diffusion-induced 7Li spin- lattice relaxation NMR spectroscopy." Physical Review B 2011, 83, 11. http://dx.doi.org/10.1103/PhysRevB.83.094302 100 [118] M. Matsui, K. Takahashi, K. Sakamoto, A. Hirano, Y. Takeda, O. Yamamoto, et al. "Phase stability of a garnet-type lithium ion conductor Li7La3Zr2O12." Dalton Transactions 2014, 43, 1019-1024. http://dx.doi.org/10.1039/c3dt52024b [119] S. Narayanan, A. K. Baral, and V. Thangadurai "Dielectric characteristics of fast Li ion conducting garnet-type Li5+2xLa3Nb2-xYxO12 (x=0.25, 0.5 and 0.75)." Physical Chemistry Chemical Physics 2016, 18, 15418-15426. http://dx.doi.org/10.1039/c6cp02287a [120] C. Monroe and J. Newman "The impact of elastic deformation on deposition kinetics at lithium/polymer interfaces." Journal of the Electrochemical Society 2005, 152, A396- A404. http://dx.doi.org/10.1149/1.1850854 [121] J. E. Ni, E. D. Case, J. S. Sakamoto, E. Rangasamy, and J. B. Wolfenstine "Room temperature elastic moduli and Vickers hardness of hot-pressed LLZO cubic garnet." Journal of Materials Science 2012, 47, 7978-7985. http://dx.doi.org/10.1007/s10853-012- 6687-5 [122] K. Lee, D. Yoo, W. Jeong, and S. Han "SIMPLE-NN: An efficient package for training and executing neural-network interatomic potentials." Computer Physics Communications 2019, 242, 95-103. http://dx.doi.org/10.1016/j.cpc.2019.04.014 [123] L. Himanen, M. O. J. Jager, E. V. Morooka, F. F. Canova, Y. S. Ranawat, D. Z. Gao, et al. "DScribe: Library of descriptors for machine learning in materials science." Computer Physics Communications 2020, 247, 12. http://dx.doi.org/10.1016/j.cpc.2019.106949 [124] W. W. Li, Y. Ando, E. Minamitani, and S. Watanabe "Study of Li atom diffusion in amorphous Li3PO4 with neural network potential." Journal of Chemical Physics 2017, 147, 10. http://dx.doi.org/10.1063/1.4997242 [125] W. Jeong, K. Lee, D. Yoo, D. Lee, and S. Han "Toward Reliable and Transferable Machine Learning Potentials: Uniform Training by Overcoming Sampling Bias." Journal of Physical Chemistry C 2018, 122, 22790-22795. http://dx.doi.org/10.1021/acs.jpcc.8b08063 [126] F. Aguesse, W. Manalastas, L. Buannic, J. M. L. del Amo, G. Singh, A. Llordes, et al. "Investigating the Dendritic Growth during Full Cell Cycling of Garnet Electrolyte in Direct Contact with Li Metal." Acs Applied Materials & Interfaces 2017, 9, 3808-3816. http://dx.doi.org/10.1021/acsami.6b13925 [127] M. Van den Bossche, H. Gronbeck, and B. Hammer "Tight-Binding Approximation- Enhanced Global Optimization." Journal of Chemical Theory and Computation 2018, 14, 2797-2807. http://dx.doi.org/10.1021/acs.jctc.8b00039 [128] S. Plimpton "Fast Parallel Algorithms for Short-Range Molecular Dynamics." J. Computa. Phys. 1995, 117, 1-19. http://dx.doi.org/10.1006/jcph.1995.1039 101 [129] B. Stanje, D. Rettenwander, S. Breuer, M. Uitz, S. Berendts, M. Lerch, et al. "Solid Electrolytes: Extremely Fast Charge Carriers in Garnet-Type Li6La3ZrTaO12 Single Crystals." Annalen Der Physik 2017, 529, 9. http://dx.doi.org/10.1002/andp.201700140 [130] T. Thompson, A. Sharafi, M. D. Johannes, A. Huq, J. L. Allen, J. Wolfenstine, et al. "A Tale of Two Sites: On Defining the Carrier Concentration in Garnet-Based Ionic Conductors for Advanced Li Batteries." Advanced Energy Materials 2015, 5, 9. http://dx.doi.org/10.1002/aenm.201500096 [131] K. Kataoka and J. Akimoto "High Ionic Conductor Member of Garnet-Type Oxide Li6.5La3Zr1.5Ta0.5O12." Chemelectrochem 2018, 5, 2551-2557. http://dx.doi.org/10.1002/celc.201800679 [132] K. Hayamizu, Y. Terada, K. Kataoka, and J. Akimoto "Toward understanding the anomalous Li diffusion in inorganic solid electrolytes by studying a single-crystal garnet of LLZO-Ta by pulsed-gradient spin-echo nuclear magnetic resonance spectroscopy." Journal of Chemical Physics 2019, 150, 9. http://dx.doi.org/10.1063/1.5089576 [133] Y. Jin, K. Liu, J. L. Lang, D. Zhuo, Z. Y. Huang, C. A. Wang, et al. "An intermediate temperature garnet-type solid electrolyte-based molten lithium battery for grid energy storage." Nature Energy 2018, 3, 732-738. http://dx.doi.org/10.1038/s41560-018-0198-9 [134] K. Xu "Nonaqueous liquid electrolytes for lithium-based rechargeable batteries." Chemical Reviews 2004, 104, 4303-4417. http://dx.doi.org/10.1021/cr030203g [135] K. M. Diederichsen, H. G. Buss, and B. D. McCloskey "The Compensation Effect in the Vogel-Tammann-Fulcher (VTF) Equation for Polymer-Based Electrolytes." Macromolecules 2017, 50, 3832-3841. http://dx.doi.org/10.1021/acs.macromol.7b00423 [136] G. Y. Gu, S. Bouvier, C. Wu, R. Laura, M. Rzeznik, and K. M. Abraham "2-Methoxyethyl (methyl) carbonate-based electrolytes for Li-ion batteries." Electrochimica Acta 2000, 45, 3127-3139. http://dx.doi.org/10.1016/s0013-4686(00)00394-7 [137] X. F. He, Y. Z. Zhu, and Y. F. Mo "Origin of fast ion diffusion in super-ionic conductors." Nature Communications 2017, 8. http://dx.doi.org/10.1038/ncomms15893 [138] M. Mottet, A. Marcolongo, T. Laino, and I. Tavernelli "Doping in garnet-type electrolytes: Kinetic and thermodynamic effects from molecular dynamics simulations." Physical Review Materials 2019, 3. http://dx.doi.org/10.1103/PhysRevMaterials.3.035403 [139] R. Murugan, V. Thangadurai, and W. Weppner "Fast lithium ion conduction in garnet-type Li7La3Zr2O12." Angewandte Chemie-International Edition 2007, 46, 7778-7781. http://dx.doi.org/10.1002/anie.200701144 [140] S. Sasano, R. Ishikawa, K. Kawahara, T. Kimura, Y. H. Ikuhara, N. Shibata, et al. "Grain boundary Li-ion conductivity in (Li0.33La0.56)TiO3 polycrystal." Applied Physics Letters 2020, 116, 4. http://dx.doi.org/10.1063/1.5141396 102 [141] B. B. Chen, C. Q. Xu, and J. Q. Zhou "Insights into Grain Boundary in Lithium-Rich Anti- Perovskite as Solid Electrolytes." Journal of the Electrochemical Society 2018, 165, A3946-A3951. http://dx.doi.org/10.1149/2.0831816jes [142] J. A. Dawson, P. Canepa, T. Famprikis, C. Masquelier, and M. S. Islam "Atomic-Scale Influence of Grain Boundaries on Li-Ion Conduction in Solid Electrolytes for All-Solid- State Batteries." Journal of the American Chemical Society 2018, 140, 362-368. http://dx.doi.org/10.1021/jacs.7b10593 [143] J. A. Dawson, P. Canepa, M. J. Clarke, T. Famprikis, D. Ghosh, and M. S. Islam "Toward Understanding the Different Influences of Grain Boundaries on Ion Transport in Sulfide and Oxide Solid Electrolytes." Chemistry of Materials 2019, 31, 5296-5304. http://dx.doi.org/10.1021/acs.chemmater.9b01794 [144] S. Yu and D. J. Siegel "Grain Boundary Contributions to Li-Ion Transport in the Solid Electrolyte Li7La3Zr2O12 (LLZO)." Chemistry of Materials 2017, 29, 9639-9647. http://dx.doi.org/10.1021/acs.chemmater.7b02805 [145] H. Shiiba, N. Zettsu, M. Yamashita, H. Onodera, R. Jalem, M. Nakayama, et al. "Molecular Dynamics Studies on the Lithium Ion Conduction Behaviors Depending on Tilted Grain Boundaries with Various Symmetries in Garnet-Type Li7La3Zr2O12." Journal of Physical Chemistry C 2018, 122, 21755-21762. http://dx.doi.org/10.1021/acs.jpcc.8b06275 [146] B. Gao, R. Jalem, H. K. Tian, and Y. Tateyama "Revealing Atomic-Scale Ionic Stability and Transport around Grain Boundaries of Garnet Li7La3Zr2O12 Solid Electrolyte." Advanced Energy Materials 2021, 12, 202102151. http://dx.doi.org/10.1002/aenm.202102151 [147] M. A. Tschopp, K. N. Solanki, F. Gao, X. Sun, M. A. Khaleel, and M. F. Horstemeyer "Probing grain boundary sink strength at the nanoscale: Energetics and length scales of vacancy and interstitial absorption by grain boundaries in alpha-Fe." Physical Review B 2012, 85, 21. http://dx.doi.org/10.1103/PhysRevB.85.064108 [148] M. Rajagopalan, M. A. Tschopp, and K. N. Solanki "Grain Boundary Segregation of Interstitial and Substitutional Impurity Atoms in Alpha-Iron." Jom 2014, 66, 129-138. http://dx.doi.org/10.1007/s11837-013-0807-9 [149] W. E. Tenhaeff, E. Rangasamy, Y. Y. Wang, A. P. Sokolov, J. Wolfenstine, J. Sakamoto, et al. "Resolving the Grain Boundary and Lattice Impedance of Hot-Pressed Li7La3Zr2O12 Garnet Electrolytes." Chemelectrochem 2014, 1, 375-378. http://dx.doi.org/10.1002/celc.201300022 [150] B. J. Morgan "Lattice-geometry effects in garnet solid electrolytes: a lattice-gas Monte Carlo simulation study." Royal Society Open Science 2017, 4, 21. http://dx.doi.org/10.1098/rsos.170824 103 [151] J. Goo, X. Y. Guo, Y. T. Li, Z. H. Ma, X. X. Guo, H. Li, et al. "The Ab Initio Calculations on the Areal Specific Resistance of Li-Metal/Li7La3Zr2O12 Interphase." Advanced Theory and Simulations 2019, 2, 7. http://dx.doi.org/10.1002/adts.201900028 [152] B. Gao, R. Jalem, and Y. Tateyama "Surface-Dependent Stability of the Interface between Garnet Li7La3Zr2O12 and the Li Metal in the All-Solid-State Battery from First-Principles Calculations." Acs Applied Materials & Interfaces 2020, 12, 16350-16358. http://dx.doi.org/10.1021/acsami.9b23019 [153] X. G. Han, Y. H. Gong, K. Fu, X. F. He, G. T. Hitz, J. Q. Dai, et al. "Negating interfacial impedance in garnet-based solid-state Li metal batteries." Nature Materials 2017, 16, 572- +. http://dx.doi.org/10.1038/nmat4821 104