Effects of Gallium Doping in Garnet-Type Li7La3Zr2O12 Solid Electrolytes

Garnet-type Li7La3Zr2O12 (LLZrO) is a candidate solid electrolyte material that is now being intensively optimized for application in commercially competitive solid state Li ion batteries. In this study we investigate, by forcefield-based simulations, the effects of Ga doping in LLZrO. We confirm the stabilizing effect of Ga on the cubic phase. We also determine that Ga addition does not lead to any appreciable structural distortion. Li site connectivity is not significantly deteriorated by the Ga addition (>90% connectivity retained up to x = 0.30 in Li7−3xGaxLa3Zr2O12). Interestingly, two compositional regions are predicted for bulk Li ion conductivity in the cubic phase: (i) a decreasing trend for 0 ≤ x ≤ 0.10 and (ii) a relatively flat trend for 0.10 < x ≤ 0.30. This conductivity behavior is explained by combining analyses using percolation theory, van Hove space time correlation, the radial distribution function, and trajectory density.


■ INTRODUCTION
Solid electrolytes with high lithium ionic conductivity are now being actively investigated for application in commercially competitive all-solid-state rechargeable lithium-ion batteries. Among them, Li-based garnet oxides have shown promise for meeting the much needed safety and reliability requirements of today's commercial lithium ion batteries. 1−7 One of the garnet compositions being considered is the cubic Li 7 La 3 Zr 2 O 12 (LLZrO), this is due to its stability with elemental lithium and a total conductivity on the order of 10 −4 S/cm. 2 Its structure is usually defined in the Ia3̅ d space group with Li cations partiallly occupying 24d tetrahedral (Td) and 48g/96h octahedral (Oh) sites, La cations fully occupying 24c dodecahedral sites, Zr cations fully occupying 16a octahedral sites, and O anions fully occupying the 96h sites (see Figure 1). However, a more stable tetragonal symmetry (I4 1 /acd) has also been reported to exist at room temperature, with Li ordering in three crystallographic sites: one at the 8a site, which forms a subset of the cubic garnet 24d site, and the other two highly distorted 16f and 32g sites which, when combined, are equivalent to the cubic garnet 48g/96h sites. 8,9 This ordering is responsible for the low bulk Li + ion conductivity measured for tetragonal LLZrO, a value on the order of 1 × 10 −6 S/cm at room tempertaure. 8 Conductivity in LLZrO is known to vary depending on the synthesis route employed, with methods that are carried out at low 10−14 and high temperatures 15−24 yielding different values. Optimization efforts have focused on obtaining the cubic phase by promoting disorder across the Li sublattice; dopant incorporation is typically employed, targeting the framework cation sites (i.e., La and Zr). Some notable improvements were reported for Te-doped LLZrO (1.02 × 10 −3 S/cm), 15 Tadoped LLZrO (1.0 × 10 −3 S/cm), 16−19 Nb-doped LLZrO (8.0 × 10 −4 S/cm), 20 Sb-doped LLZrO (7.7 × 10 −4 S/cm), 21 and Sr-doped LLZrO (5.0 × 10 −4 S/cm). 22 Simultaneous substitution was also explored in the series Li 7+x−y (La 3−x A x )-(Zr 2−y Nb y )O 12 (where A is an alkali earth metal) and has shown that an optimum lattice parameter at a constant lithium content of 6.5 per formula unit (p.f.u.) can be obtained. 23 On the other hand, several simulation studies for LLZrO have been made almost hand in hand with experiments, focusing on understanding Li + ion diffusion mechanisms, phase transition, and stability. For example, in our previous density functional theory (DFT) molecular dynamics (MD) calculations, we found that the complex mechanism for self-diffusion of Li + ions in LLZrO proceeds in a cooperative manner and is governed by two crucial features: (i) the restriction imposed for occupied site-to-site interatomic separation and (ii) the apparent unstable residence of the Li + ion at the Td site due to local Li−Li repulsion effect. 25 Meier et al. further provided information to establish the difference between the Li + ion diffusion mechanism of the tetragonal and cubic phase through DFT-based MD and metadynamics simulations; the former has a fully collective nature or synchronous motion while the latter has an asynchronous motion governed by single-ion jumps and induced collective motion. 26 In another simulation using classical MD, with a BV-based Morse-type force field, Adams et al. predicted that for the garnet Li 7−x La 3 (Zr 2−x M x )O 12 system (x = 0, 0.25; M = Ta 5+ , Nb 5+ ): (i) the lithium distribution just above the cubic phase transition closely resembles that in the tetragonal phase and that (ii) pentavalent doping can enhance ionic conductivity by increasing the vacancy concentration and by reducing local Li ordering. 27 Wang et al. discovered through static energy minimization with Buckingham potentials, that the shape of energy probability distribution, can aid in unders t a n d i n g l i t h i u m d i s o r d e r / o r d e r e ff e c t s i n t h e Li 7−x La 3 Zr 2−x Ta x O 12 (x = 0−2) series. 28 Miara et al. investigated the effect of Rb and Ta doping on the ionic conductivity and stability of the garnet Li 7+2x−y (La 3−x Rb x )-(Zr 2−y Ta y )O 12 (0 ≤ x ≤ 0.375, 0 ≤ y ≤1) using ab initio-based calculations and concluded that (i) Rb or Ta doping does not change the topology of the migration pathways significantly, but instead acts primarily to change the lithium concentration, and that (ii) doping with larger cations will not provide a significant enhancement in performance. 29 Bernstein et al. has revealed by DFT and variable cell shape MD simulations that the strong dependence of the tetragonal phase stability on the simultaneous ordering of the Li + ions on the Li sublattice and a volume-preserving tetragonal distortion that relieves internal structural strain. 30 32 Overall, these efforts have significantly contributed to a better understanding of the measured properties of LLZrO-based materials.
Most doping strategies have been geared toward tuning conductivity without obstructing Li pathways, as highlighted above. However, substitution in LLZrO had also been carried out on the Li sublattice itself, directly impeding the conduction path of Li + ions. The most studied example of this is the intentional or unintentional chemical substitution of Li + by Al 3+ . 33−39 Two Li site vacancies are created for every addition of Al 3+ and this is made possible primarily by the small ionic radius of Al 3+ in the two Li coordination environments (r Al3+ = 0.39 Å vs r Li+ = 0.59 Å for tetrahedral coordination and r Al3+ = 0.535 Å vs r Li+ = 0.76 Å for octahedral coordination). However, to date, the site preference for Al 3+ has been ambiguous. Neutron diffraction data measurement made by Li et al. suggested an octahedral (48g)-site occupancy 35 but Wang et al. 36 and Buschmann et al. 37 concluded, based on their 27 Al MAS NMR spectroscopy results, that Al sits mostly in the 24d tetrahedral site. Geiger et al. 38 added that, according to their 27 Al MAS NMR spectral analysis, Al 3+ may also occupy a distorted site with 5-fold coordination which is presumably at the octahedral site; one of the two main resonances they found, which is assigned to an octahedral environment for Al 3+ , was attributed to the LaAlO 3 impurity phase. Similarly, Rettenwander et al. showed, by combining 27 Al NMR data and DFT calculation, that Al 3+ could sit, not only in a 24d site but also in a distorted 4-fold coordinated 96h site and even a 48g site with almost regular octahedral coordination; they also suggested that octahedral occupancy for Al 3+ may be stabilized depending on a number of factors such as dopant concentration, sintering temperature and time, heating rate, grain sizes, starting materials, etc. 39 Duvel et al. prepared Al-doped LLZrO samples by a mechanochemical route and discovered that with increasing Al content, Al ions can also substitute at non-Li cation sites, leading to the formation of La-and Zr-deficient LLZrO. 40 Huang et al. indicated a similar effect for Ge-doped LLZrO, according to a comparison of their experimental and simulated XANES spectra. 41 Recently, Ga-doped cubic LLZrO was successfully prepared by standard solid state synthesis 42−44 and a low-temperature sol gel approach. 45 Similar with Al 3+ , the ionic radius of Ga 3+ (0.47 and 0.62 Å for tetrahedral and octahedral coordination, respectively) is also comparable with that of Li + . This system has demonstrated comparable and even better structural and transport properties than the Al case. 42,43 In common with Aldoped LLZrO, the Ga distribution in the LLZrO framework is also difficult to ascertain. 71 Ga NMR spectroscopy conducted by Howard et al. for Li 5.5 Ga 0.5 La 3 Zr 2 O 12 (1 × 10 −4 S/cm) revealed a single broad peak with a chemical shift of ∼221 ppm which was assigned to Ga 3+ at tetrahedral sites. 42 Bernuy-Lopez et al. came to the same conclusion via their 71 Ga NMR spectrum for Li 5.5 Ga 0.15 La 3 Zr 2 O 12 which was sintered in a dry O 2 atm (1.3 × 10 −3 S/cm), from which a chemical shift was obtained with a value of 207 ± 10 ppm. 43 On the other hand, in the range 0.08 ≤ x ≤ 0.84 for Li 7−3x Ga x La 3 Zr 2 O 12 , Rettenwander et al. suggested, based on the relatively large NMR asymmetry parameter that they derived, η Q = 0.46 (3), that Ga 3+ mainly occupies a distorted 4-fold coordinated environment corresponding to the general 96h octahedral site. 44 Allen et al. 18 also proposed the 96h octahedral site occupancy, following the work of Geller et al. 46 who reported that Ga 3+ preferred 6-fold coordination in garnet.

Chemistry of Materials
Until now, except for the stabilizing effect on the cubic phase, there is no clear understanding of how foreign ion incorporation in the Li sublattice affects the structure, path topology and Li + ion dynamics of doped LLZrO. This is an important issue since the blocking effect by aliovalent dopants can progressively reduce path connectivity but conversely may also promote Li disordering for cubic phase stabilization due to the variation in the number of Li vacancies ( Figure 1 shows the linkage of Td and Oh sites). In this study, we investigated the effects of Ga doping in the solid solution series Li 7−3x Ga x La 3 Zr 2 O 12 (0 ≤ x ≤ 0.3) using atomistic simulations with Buckingham-type interatomic potentials. Our results offer practical insights beneficial to experimentalists, especially when formulating strategies for conductivity optimization in doped LLZrO materials.

■ COMPUTATIONAL DETAILS
Derivation of Interatomic Potentials. The crystal lattice was modeled using a classical Born description. 47 The total energy of the system was calculated by summing contributions from long-range electrostatic, short-range repulsive and van der Waals interactions. The long-range electrostatic interaction is given by Coulomb's law where q i and q j are the charges of ions i and j, respectively, ϵ 0 is the free space permittivity, and r ij is the distance between the ions. U ij was evaluated via the Ewald approach. 48 Nonformal charges were used, based on a partially ionic model; the ions' fractional effective charges (q eff ) were calculated from the product of their nominal charge (q nom ) and a defined ionicity parameter ζ, namely where ζ < 1. These effective charges are proportional to their nominal charges to ensure charge neutrality. For the short-range interactions, Φ ij , the Buckingham potential form was used 49 where A ij , ρ ij , and C ij are empirically derived parameters for interacting ions i and j. The cutoff distance, beyond which short-range interactions are considered negligible, was set to 10.5 Å. Using ζ = 0.70, the potential parameters (A, ρ, and C) were fitted against experimental lattice constants of several relevant oxides. Essentially, an error-based objective function, E fit , is minimized during the fitting procedure, as given by where N expt is the number of experimental parameters used to evaluate the fit, w i is the weighting factor (set to 1 for all parameters for equal weighting), and x i expt and x i calcd are the experimental and calculated parameter values, respectively. E fit was minimized using the Nelder− Mead simplex algorithm. 50 The constant-pressure energy minimization routine, in which the dimensions of the simulation box are relaxed together with atom positions, was performed with the GULP code using the Broyden−Fletcher−Goldfarb−Shanno algorithm. 51 Li−Ga Vacancy Configuration Sampling. In this study, a 3 × 3 × 3 supercell (containing 1944 available Li sites, 648 La atoms, 432 Zr atoms, and 2592 O atoms). For simplicity, the tetrahedral 24d and octahedral 48g/96h site cage were labeled as Td and Oh site, respectively. There exists an enormous number of conceivable arrangements for the Li, Li vacancy, and Ga species within the Li sublattice of a basic cubic (I) garnet structure, this increases even further with increasing model cell size. For practical reasons, a random sampling approach for the Li−Ga vacancy arrangements was used, the number of Li and Ga atoms were adjusted according to stoichiometry and with Ga content x varied in the range 0 ≤ x ≤ 0.30. At each x, a total of 16 000 random structures spanning from all-Td to all-Oh Ga occupancy were evaluated by a two-step screening procedure. In the first step, electrostatic/Coulomb energies were calculated. 48 In the second step, the lowest Coulomb energy structures for each discrete Ga occupancy ratio, denoted by Td/(Td+Oh), from the first step were optimized using the fitted Buckingham potentials.
Molecular Dynamics Simulation. All MD simulations were run with the DL_POLY code 52 with a time step of 1 fs. Three statistical mechanical ensembles were used, namely: (i) the microcanonical ensemble which maintains constant number of particles, volume, and energy (NVE), (ii) the canonical ensemble with constant number of particles, volume, and temperature (NVT), and (iii) the isothermal− isobaric ensemble with constant number of particles, pressure, and temperature (NPT). To ensure that the average temperature, pressure, and/or stress tensor are maintained close to the target conditions, the Nose-Hoover thermostat (for the NVT MD run) and barostat (for the NPT MD run) were used, with 0.8 and 1.2 ps as relaxation times, respectively. 53 A heating schedule was followed for MD production sampling: (i) 10 ps NVE MD run at 300 K with 5 ps equilibration, (ii) temperature ramping in the NPT ensemble from 300 K up to the target temperature incrementing at a rate of 100 K per 30 ps, and (iii) 500 ps NVT MD run at the target temperature. For the thermal expansion/contraction analysis, a similar heating schedule was followed but instead the run was performed in the NPT ensemble at a rate of 30 ps per 50 K step, with 100 ps sampling at the target temperature. Target pressure for constant pressure ensembles were set to 1 atm. Trajectory collection was carried out every 10 MD steps.

■ RESULTS AND DISCUSSION
Configuration Sampling. Incorporation of Ga onto the Li sublattice of the garnet LLZrO structure can be described using Kroger−Vink notation as follows where Ga Li ·· indicates Ga at the Li site with an effective defect charge of +2, O O X is O at the O site with a neutral charge, and V Li ′ is a vacancy at the Li site with an effective charge of −1. The Li−Li vacancy-Ga lattice configuration sampling space gives rise to a Coulomb energy distribution as displayed in Figure 2. As the Ga content increases (up to x = 0.30), the energy distribution related to the long-range interaction tends to widen and becomes more negative. This can be attributed to the decreased local repulsion in the Li sublattice due to the decreased cation occupancy and increasing mean separation of Li + -Li + and Li + -Ga 3+ . Meanwhile, the total (static) energy deviation (i.e., Coulomb energy + Buckingham potential term) is also found to be small with respect to Ga configuration for x   0.3 (less than 10 meV/atom with a variation of less than 15 meV/atom per x, see Figure S1 of the Electronic Supporting Information). This implies that Ga 3+ could sit in both Td and Oh sites without significant energy penalty, partly explaining the lack of consensus on the correct Ga 3+ site preference in the literature. 18,[42][43][44]46 Even with the inclusion of lattice vibration effects, this tendency is still preserved (see Figure S2 in the Supporting Information).

Chemistry of Materials
Unless specified otherwise, the lowest total energy structures at each value of x (Ga content) were used for the subsequent MD runs.
Effect of Ga 3+ toward Li Distribution, Structural Distortion, and Topology. As previously mentioned, it was suggested in experimental reports that Ga 3+ leads to the stabilization of cubic phase for LLZrO. To confirm this, we analyzed the Li distribution with respect to Ga 3+ addition. Figure 3a shows the Li−Li radial distribution function g(r) profiles of undoped (x = 0) and Ga-doped (x = 0.30) 3 × 3 × 3 supercell models for LLZrO (at 800 and 300 K). Vertical lines (dashed and solid) show the discrete peak locations for the tetragonal phase (note that due to cell edge mismatch, the first two peak lines from the left are shifted to the right with respect to the cubic phase). It should also be mentioned at this point that owing to the scale of the phase space involved (1944 Li sites for a 3 × 3 × 3 supercell), standard MD sampling and annealing would require a prohibitively long computation time to reach the exact Li ordering of the tetragonal phase (i.e., into the 8a, 16f, and 32g sites 25 ) when starting at the high-T Li disordered cubic phase. In other words, in the present computation, we would only observe the tendency for Li ordering in the cooling direction, that is, tetragonal signatures in the Li−Li g(r) profiles become more pronounced as we go down to low temperature range (starting from a high temperature phase) but the exact g(r) fingerprint of the fully ordered Li arrangement would not be reached. This is explained as follows: if the tetragonal phase is formed, the octahedral sites should be fully occupied and contributions for 16f Li−32g Li pairs should become prominent within the range 3.0 Å < r < 4.0 Å (first nearest neighbor, or first NN), 4.8 Å < r < 5.5 Å (2nd NN), and 6.3 Å < r < 7.0 Å (3rd NN); see the vertical solid lines in Figure 3a. Also, 1/3 of the Td sites should be occupied in a regular arrangement such that Td Li + do not have a first nearest Td site neighbors, that is, 8a Li−16e Li pair contribution (for first NN Td site pairs) should become weaker at r ≈ 4 Å. 25 In both models (x = 0 and x = 0.30), the first peak can be readily assigned to the characteristic Td-Oh Li intersite contribution. At 800 K, Li disordering is depicted in Figure 3a by the smoother peaks due to the ease of Li diffusivity ( Figure  3c, e). At 300 K, Li + becomes localized at specific sites, shown as peak narrowing ( Figure 3c) and peak splitting (Figure 3e) vis-a-vis the 800 K g(r) plots. Peaks for the 16f Li−32g Li pair contribution are slightly sharper (for first and second NN Oh Li + ) in the undoped model as compared to the doped case at 300 K; this stems from the increased occupancy at octahedral Li sites for the former. In addition, a decrease in g(r) becomes evident at r ≈ 4 Å, this can be related to a weakening 8a Li−16e Li pair contribution. Collectively, these observations clearly point to a tendency for Li ordering in the undoped LLZrO model at 300 K for both heating and cooling directions. Meanwhile, the low g(r) at r ≈ 4 Å for the doped model cannot be attributed to the Li ordering effect from the diminished first NN Td site Li pair contribution, but instead due to Ga occupancy in those sites (the selected lowest energy structure, such as for x = 0.30, mainly has Ga 3+ at Td sites). Because of Figure 3. (a) Calculated Li−Li radial distribution function g(r) profiles for undoped (x = 0; 300 and 800 K in the heating direction and 300 K in the cooling direction) and doped (x = 0.30; 800 K in the heating direction and 300 K in the cooling direction) 3 × 3 × 3 LLZrO cell models, (c) enlarged view of the g(r) plot in a highlighting the narrowing of the peak from 800 to 300 K in the range 1.5 Å ≤ r ≤ 3.0 Å (see arrows), and (e) enlarged view of the g(r) plot in a highlighting the peak splitting from 800 K (*) to 300 K (+) in the range 4.1 Å ≤ r ≤ 5.5 Å (see arrows). Vertical lines (dashed and solid) indicate discrete peak locations for the tetragonal phase (8a-16f-32g). Thermal expansion plots with respect to the monitored reaction coordinates (lattice constants c and a) for: (b) undoped, (d) 1-Ga atom, and (f) 2-Ga atom 1 × 1 × 1 (I cell) models (Ga 3+ added at the 24d site). Dashed lines show linear fitting for 600 K ≤ T ≤ 1000 K, in both heating (red) and cooling (blue) direction. this, Li + ions cannot assume the 8a-16f-32g arrangement for tetragonal phase formation in the Ga-doped model. Other evidence for the lack of ordering tendency in the Ga-doped model at 300 K are the broader peaks, the absence of a valley region at around r ≈ 3.0 Å, and shallower valleys at 4.5 Å < r < 5.0 Å and r > 6.0 Å. Thus, we verify from the g(r) plots the lack of Li ordering tendency and consequently, the suppression of the tetragonal phase in the Ga-doped model.

Chemistry of Materials
Another way of distinguishing the loss of tetragonal character, on Ga addition, is by monitoring the thermal expansion/contraction plots of undoped and doped LLZrO models. We started from a smaller 1 × 1 × 1 tetragonal (I4 1 / acd) cell (72 Li sites) with no Ga, 1 Ga, and 2 Ga atom(s). For the 1-Ga cell, we replaced 1 Li atom at the 8a site and made sure the first NN sites (16f or 32g) are empty, to minimize local repulsion (referred from the occupancy trend shown by lowest total energy structures). Similarly, for the 2-Ga cell, we replaced 2 Li atoms at the 8a site. The number of Li atoms was then adjusted accordingly for charge neutrality. We then monitored the cell edges, with the shortest one denoted as the lattice constant c and the average of the remaining two as lattice constant a, at every temperature step, and used them as key reaction coordinates for detecting the phase transition. We employed this condition because c can switch between the three cell edges because of the Li redistribution and cell symmetry, especially below the transition temperature T c . If the chosen coordinate only changes linearly with T, then only Li disordering exists and there is no phase transition. If a bend in the curve is observed, however, then it is identified as the tetragonal → cubic topotactic structural transition. Figure 3b, d, and f shows the evolution of lattice constants c and a with respect to temperature for an undoped cell, a cell with 1 Ga atom, and a cell with 2 Ga atoms, respectively. For the undoped cell, a bend in the thermal expansion plot is observed both in the heating (red) and the cooling (blue) direction. The onset of bending appears to start at ∼600 K; beyond this temperature the trend becomes linear. We assigned this bend as the phase transition point. This value is in good agreement with some experiments (>400 K) 27,38 and with a recent classical MD simulation (>450 K). 27 On the other hand, the bend in the thermal data in Figure 3b has disappeared in both the heating and cooling direction in Figure 3d, f (1 and 2 Ga atoms added, respectively), which signifies the suppression of the tetragonal phase. Thus, we again confirmed our prediction that the cubic phase can be effectively stabilized by Ga 3+ . The model cell we used with 2 Ga atoms (in an I cell) is equivalent to a doping level of x = 0.25, but stabilization has also been achieved experimentally at lower Ga-doping levels than this. 44 To determine the extent of cell expansion/contraction with increasing Ga 3+ , we plotted the average lattice constant with respect to x in Figure 4a. The plot shows good agreement between calculated and XRD-derived lattice constants (<1% difference). Also, Ga 3+ doping appears not to cause any severe lattice constant change; several experimental studies related to Li pathway blocking dopants for LLZrO have arrived at the same conclusion. 13,38,42,44,54 We also checked whether Ga can promote a more open local Li pathway by analyzing the void space when Li atoms are completely removed in the model cells. The Zeo++ software was used to perform geometry-based analysis of structure and topology of the void space inside the LLZrO framework; it can determine the bottleneck size by calculating the largest pore diameter along the Li pathway. 55 As a recap of the Li path topology, four Td faces related to the Td site (the pathway junction) form the path bottlenecks. Each of these faces are bounded by three O atoms as shown in Figure  4b. Results for the pathway bottleneck size variation are collected in Figure 4c. Clearly, no significant change is observed for the bottleneck size, the same with cell edges (Figure 4a)  Another important factor to consider is the extent of connectivity disruption in the Li sublattice by progressive blocking of Ga 3+ . For this analysis, we treated the system as a three-dimensional percolation problem. 56 Two cases need to be evaluated in order to reveal the site connectivity relationship: exclusive Td-site blocking and a combination of Td-and Ohsite blocking. A sufficiently large supercell (10 × 10 × 10) was used with 24 000 and 72 000 sites, respectively, in order to minimize periodic boundary condition effects. The results are plotted in Figure 4d. For the exclusive Td-site blocking case, the percolation threshold (p c ) for Ga 3+ is determined to be about ∼2.34 mol p.f.u. which in terms of the stoichiometry would lead to the full removal of Li + (∼7.02 mol p.f.u.). Up to x = 0.30 however, site connectivity is still >90%. For the simultaneous Td-and Oh-site blocking case, as expected, p c occurs at a much higher Ga content (∼4.68 mol p.f.u.), which is also way beyond the stoichiometric limit for Ga−Li substitution. On the basis of these results, we argue that siteblocking by Ga 3+ does not severely deteriorate connectivity and Li + can still form a contiguous path extending from one side of the cell to the other.
Lithium Ion Dynamics and Diffusion Mechanism. The motion of atoms were analyzed from their mean square displacements (MSDs). Li MSD profiles for undoped (x = 0) and doped (x = 0.30) LLZrO in Figure 5a, b, respectively, display increasing linear trends vs MD sampling time and diffusivity vs T. Framework atoms (La, Zr, and O) are noted to essentially just vibrate about their sites as confirmed by their zero-slope MSD trends (see Figures S2 and S3 in the Supporting Information). These results are consistent with other simulation studies. 25,26,29 In the case of Ga atoms (for x = 0.30), their MSD starts to increase at 1200 K (Figure 5c). This onset temperature for Ga ion diffusion does not vary with respect to Ga content (see Figure S4 in the Supporting Information). In general, the integrity of the garnet structure has been kept in all of the MD production runs for all doping levels considered, allowing us to reliably estimate the Li + ion diffusion coefficient and activation energy values.   To gain insights into the interplay between Li + ion diffusivity and Ga distribution, we compared the MSD slopes between exclusive Td-and exclusive Oh-blocking model cells. Note that these two extremes are predicted to be energetically comparable (see Figure S1 in Supporting Information). Results are displayed in Figure 5d for compositions x = 0.10 and x = 0.15 at 1000 K. Calculated diffusion coefficients are 2.53 × 10 −5 cm 2 /s and 2.80 × 10 −5 cm 2 /s for all-Td and all-Oh cases, respectively, for x = 0.10, and 2.33 × 10 −5 cm 2 /s and 2.63 × 10 −5 cm 2 /s for all-Td and all-Oh cases, respectively, for x = 0.15. In these two doping levels, the all-Oh Ga cases are systematically higher (by 11 and 13%, respectively) than the all-Td cases. This difference is consistent with (i) the spatial relationship of Td and Oh sites, that is, the more critical Td site forming a junction connecting four neighboring Oh sites and the Oh site only oppositely face-sharing two Td sites, and (ii) the generally preserved site connectivity of the pathway (>90%), independent of Ga site preference (see Figure 4d). Consequently, for a mix of Td-and Oh-site Ga occupancy, the increase in diffusivity relative to the all-Td Ga model is expected to have an upper bound defined by the all-Oh case.

Chemistry of Materials
The evolution of the Ga trajectory density for a Ga-doped LLZrO model (x = 0.30, 500 ps NVT MD production run with 1 fs step and 100 fs trajectory sampling) at increasing temperature, is shown in Figure 6. Below the onset temperature (T = 1000 K in Figure 6a), Ga atoms are strongly localized about their sites. At T = 1200 K in Figure 6b, the vibration signature becomes bigger and some hints of next-nearestneighbor-site hopping are becoming apparent, based on the increased Ga density. As the temperature further increases (Figure 6c, d), more Ga 3+ become delocalized from their initial sites, leading to limited Ga diffusion being possible with hopping up to third to fourth NN sites. This onset temperature may be exploited for conductivity optimization, either by investigating Ga redistribution in the bulk and/or Ga diffusion from surface to bulk at different heating rates and temperatures. Figure 7a displays the Arrhenius plots derived from MD simulations; experimental data are also included for comparison. Following the Nernst−Einstein relationship, the bulk Li + ion conductivity, σ Li , can be calculated from 57 where c Li is the carrier density for Li, z Li is Faraday's constant, R is gas constant, and T is temperature. The Li + ion activation energy (E a ) values are determined to be in the narrow range of 0.24−0.30 eV, indicating that Ga 3+ does not strongly affect the Li + diffusion barrier height. Meanwhile, room temperature bulk conductivity, as indicated in Figure 7b, is calculated to be in the range 1.42−6.08 × 10 −3 S/cm (the maximum at x = 0.02). We would like to emphasize here that the MD sampling procedure performed for x = 0 is in the high temperature region (see circle data points in Figure 7a), so the simulation box already has a cubic symmetry, and thus explaining the artificially high predicted bulk conductivity. On the other hand, in the heating direction (at 300 K), the conductivity for x = 0 with the tetragonal phase is at least 2 orders of magnitude lower in conductivity than the Ga-doped cases. It is also visible from the trend that a minimal Ga addition of x ≈ 0.05 is sufficient for the stabilization of the cubic phase. There appear to be two distinct regions in the high-T-derived conductivity plot (vs x). For 0.00 ≤ x ≤ 0.10, the conductivity tends to decrease in a linear fashion, whereas for 0.10 ≤ x ≤ 0.30, it tends to be flat. The existence of these two regions will be clarified later. At present, Bernuy-Lopez et al. have succeeded in achieving a conductivity of 1.3 × 10 −3 S/cm with their dense (94%) and air-/moisture-protected Li 6.55 Ga 0.15 La 3 Zr 2 O 12 ; 43 this value is very close to our prediction (1.6 × 10 −3 S/cm at x = 0.15). The strong agreement could mean two things: (i) that grain boundary resistance may be small at this doping level, and (ii) that higher conductivity values could still be achieved experimentally by pushing toward a much lower Ga 3+ doping concentration. However, conductivity in LLZrO-based solid electrolytes was shown to vary greatly in experiments due to differences in preparation conditions. Overall, important requirements were noted, such as ensuring high density for the pellet samples and the prevention of H 2 O/CO 2 uptake. A varying room temperature conductivity was measured by Howard et al. with 0.50 mol p.f.u. of Ga 3+ -doping and ∼80% dense samples (between 1 × 10 −3 and 1 × 10 −5 S/cm), 42 possibly with a contribution as well from H + /Li + exchange. 58 A further example of this variability are the results of Wolfenstine et al. in which they achieved ∼91% dense samples with 0.25 mol p.f.u of Ga 3+ , resulting in a total conductivity of 3.5 × 10 −4 S/cm, much lower than the 10 −3 S/cm range. 54 El Shinawi et al. 45 substituted up to 1.0 mol p.f.u. and achieved 92.5% density; they identified an increasing trend in the measured conductivity to about 5.4 × 10 −4 S/cm, although still <1 × 10 −3 S/cm. They determined that not all the Ga 3+ ions were incorporated into the garnet structure, and that Colored circle data points are MD simulation data above the predicted phase transition temperature (T c ≈ 600 K). (b) Predicted variation in bulk conductivity vs Ga content at room temperature as extrapolated from the high-temperature MD simulation data in a (colored circle data points) and as extracted from the 300 K MD run in the heating direction. a portion ended up reacting with excess Li + ions to form LiGaO 2 , which then resided at the grain boundaries; this secondary phase was argued to have aided in pellet densification. From these results, Ga 3+ can also be considered beneficial as a sintering aid. It is important to note that in all the experimental results quoted above, the total (bulk plus grain boundary) conductivity is reported, as it is difficult to resolve the individual bulk and grain boundary components from impedance spectroscopy, and thus values are extremely sensitive to the microstructure of the measured samples.

Chemistry of Materials
The two-region conductivity plot predicted in Figure 7b suggests a shift in the nature of the Li + ion diffusion within the garnet framework. To elucidate this, we checked the variation in Li−Li space-time correlation and analyzed it with other data such as Li trajectory density, number of Li vacancies, Li site blocking/connectivity/percolation and net mass transport change. For the Li−Li space-time correlation, the distinct part of the van Hove space-time correlation function was analyzed 59 where N is the number of Li+ ions, δ(•) is the threedimensional Dirac delta function, and r j and r i are displacements of particles j and i, respectively, at time t. G d (r,t) accounts for probability of finding particle j at location r after time t, in relation to the position of another particle i at the initial time t = 0. The main points relevant to the following discussion can be enumerated as follows: (i) when t = 0, G d (r,t) collapses to the more commonly known static g(r) plot with the first peak characterized by Td-Oh intersite distance, in collaborative hopping, r j (t) − r i (t = 0) progressively approaches zero with time and causes the G d (r,t) contribution near r = 0 to increase rapidly (ii) the local rearrangement of particles leads to a decreasing G d contribution over time. Figure 8 summarizes the G d (r,t) plots at different Ga content x (for 2−20 ps); corresponding initial (t = 0) g(r) plots are also presented for reference (shown as dashed curves). As can be seen, the probability build-up near r = 0 and concomitant decrease in the characteristic Td-Oh intersite distance (2 Å ≤ r ≤ 3 Å) are evident; these are characteristic features of the collaborative hopping of Li + ions. 27,60,61 On the other hand, the intermediate distance contribution to G d for second nearest site neighbors and beyond are shown to be insignificant (i.e., for r > 4 Å, there is no shift in the peaks and valleys). 27 However, the rate of increase for 1 Å ≤ r ≤ 2 Å is shown to be dependent on the doping level, that is, the G d (r,t) build-up rate generally decreases as Ga content increases.
To understand the dependence of the G d build-up rate for 1 Å ≤ r ≤ 2 Å on Ga 3+ content, we analyzed how Li + ions move within a given path segment formed by site linkages, such as Td 1 -Oh 1 -Oh 1 '-Td 2 -..., where Td1 is a 24d Td site, Oh1 is an Oh 96h splitting site linked to Td 1 , Oh 1 ′ is another 96h splitting site in the same Oh cage with Oh 1 , and Td 2 is another 24d Td site linked to Oh 1 ′ (see Figure 9a). In this path segment, the G d contribution from 1 Å ≤ r ≤ 2 Å can be accounted for by local hopping processes such as Td 1 (i)-Oh 1 (j) (Figure 9b), Oh 1 ′(i)-Td 2 (j) (Figure 9c), and cases for intermediate positions (eg., Figure 9d). For the G d build-up rate from 1 Å ≤ r ≤ 2 Å to decrease, two major effects can be considered: (i) the mean Li + −Li + separation has increased because of a reduced number of Li + charge carriers (ii) the degree of hopping by Li(j) into a vacant Li(i) site (as shown in Figure 9b−d) has been decreased by the increasing number of guest Ga 3+ cations, which trap the Li vacancies. Because there is no appreciable change in the Li−Li g(r) plot vs x (see Figure S6 in Supporting Information), the first region of the high-T-derived two-region conductivity plot in Figure 7b (0.00 ≤ x ≤ 0.10) may then be rationalized as being due to the second effect described above being imposed on the Li + -ion dynamics. As Li vacancies are expected to be strongly coupled to the substitutional Ga 3+ , these vacancies become inaccessible for a migrating Li + ion according to the following association reactions where {Ga Li ·· V Li ′ } · is a defect cluster formed by one Ga 3+ and one Li vacancy and {Ga Li ·· 2V Li ′ } X is a neutral cluster formed by one Ga 3+ and two Li vacancies. The Ga−Li g(r) plots confirm the presence of inaccessible Li sites around Ga 3+ (see Figure S5a in the Supporting Information); a radius enclosing up to the second NN sites tends to be occupied only with an average of ∼1. 20 Li + ions at any time (cutoff radius: 4.4 Å). For a classical ionic conductor, a conductivity maximum is usually observed, based on eq 4, when vacancies are progressively added. The peak in the conductivity trend for 0.00 ≤ x ≤ 0.10 (see Figure  7b) closely resembles this behavior. For the linear decrease in conductivity (for 0.05 ≤ x ≤ 0.10), higher order clusters may be involved as the Ga content is further increased, lowering the number of accessible vacancies where {2Ga Li ·· 4V Li ′ } X is a higher order neutral cluster formed by 2 Ga 3+ and 4 Li vacancies. However, the flat trend for 0.10 < x ≤ 0.30 cannot be explained in the same way, so other factors must act to prevent the conductivity decrease. A possible resolution to this, is that the Li motion becomes more directed and coordinated within the retained percolated pathways. Note that at higher x, the number of accessible or free Li vacancies (i.e., not strongly coupled with Ga 3+ ) and the average Ga 3+ -Li + repulsion are expected to increase; these are made possible by the fact that two new Li vacancies are created for every Ga 3+ ion incorporated.
The retention of percolated pathways with Ga 3+ doping (as confirmed in Figure 4d) is qualitatively shown in the Li trajectory comparison in Figure 10 between x = 0 and x = 0.30.
The 3D network pathway is clearly visible (Figure 10a, c) with a characteristic motif visible along the ⟨1 1 1⟩ direction. The trajectory of the x = 0 case is in strong agreement with a recent analysis of neutron powder diffraction data 62 as well as several simulation studies. [25][26][27]29 To make a clear comparison with the x = 0.30 case, we chose an appropriate slab section to expose the full 2D plane view of Li pathway connectivity, such as plane (3 6 4), which is shown in blue in Figure 10a and 10c. As can be seen, pathway features such as junctions formed by Td sites (LiO 4 ) and density lobes formed by Oh sites (LiO 6 ) can be observed in detail (Figure 10b, d). For the doped case (x = 0.30, Figure 10d), Li + ion transport pathways are evidently interrupted by Ga 3+ , as shown by some segments of the network illustrated by the trajectory densities being absent (for x = 0.30). However, there are still percolated pathways that remain, making long-range Li + -ion transport possible. Integration of results from the first peak of the Li−Li radial distribution function (RDF) plots g(r) have revealed no drop on the Li neighbor atom count with x (cutoff radius: 3.85 Å); ∼2.18 to ∼2. 25 Li atoms for x = 0 to x = 0.30 (see Figure S5b in the Supporting Information). This is consistent with having sustained or directed concerted Li + motion in the retained percolated pathways at higher x, as discussed in the previous paragraph.

■ CONCLUSION
The effects of Ga-doping on garnet-type LLZrO was successfully investigated through force-field-based simulations using Buckingham short-range potentials. We found that Ga 3+ can effectively stabilize the cubic phase. Ga 3+ incorporation does not change the lattice constant nor contribute to any significant structural distortion. The onset temperature for Ga diffusion is predicted to be about 1200 K. For the cubic lattice, two distinct regions in the conductivity behavior were observed for 0 ≤ x ≤ 0.30. For 0 ≤ x ≤ 0.10, a decreasing trend is noted, whereas for 0.10 ≤ x ≤ 0.30, a flat trend is observed. The former can be explained primarily by the increasing number of inaccessible Li vacancies (due to Ga 3+ −Li vacancy cluster formation) for use in initiating Li + ion migration and concerted Li + motion. The latter is likely governed by the subsequent increase in the number of accessible Li vacancies (at higher x) and the more directed motion of Li + ions in retained percolated pathways due to the increasing average Ga 3+ −Li + repulsion.

* S Supporting Information
Evaluation of fitted potential parameters, Tables S1−S3, Figures S1−S6. This material is available free of charge via the Internet at http://pubs.acs.org.

Notes
The authors declare no competing financial interest.

■ ACKNOWLEDGMENTS
R.J. is grateful for financial support from Nagoya Kogyokai Scholarship, Nagoya Institute of Technology. The present work was partially supported by JST, PRESTO-program, and MEXT program "Elements Strategy Initiative to Form Core Research Center" (Since 2012), MEXT, Ministry of Education Culture, Sports, Science and Technology, Japan. Crystal structures and isosurfaces were created using the VESTA software. 63