Ionic transport in solid electrolytes can often be approximated as ions performing a sequence of hops between distinct lattice sites. If these hops are uncorrelated, quantitative relationships can be derived that connect microscopic hopping rates to macroscopic transport coefficients; i.e. tracer diffusion coefficients and ionic conductivities. In real materials, hops are uncorrelated only in the dilute limit. At non-dilute concentrations, the relationships between hopping frequency, diffusion coefficient and ionic conductivity deviate from the random walk case, with this deviation quantified by single-particle and collective correlation factors, f and f_I, respectively. These factors vary between materials, and depend on the concentration of mobile particles, the nature of the interactions, and the host lattice geometry. Here, we study these correlation effects for the garnet lattice using lattice-gas Monte Carlo simulations. We find that, for non-interacting particles (volume exclusion only), single-particle correlation effects are more significant than for any previously studied three-dimensional lattice. This is attributed to the presence of two-coordinate lattice sites, which causes correlation effects intermediate between typical three-dimensional and one-dimensional lattices. Including nearest-neighbour repulsion and on-site energies produces more complex single-particle correlations and introduces collective correlations. We predict particularly strong correlation effects at x_Li=3 (from site energies) and x_Li=6 (from nearest-neighbour repulsion), where x_Li=9 corresponds to a fully occupied lithium sublattice. Both effects are consequences of ordering of the mobile particles. Using these simulation data, we consider tuning the mobile-ion stoichiometry to maximize the ionic conductivity, and show that the ‘optimal’ composition is highly sensitive to the precise nature and strength of the microscopic interactions. Finally, we discuss the practical implications of these results in the context of lithium garnets and other solid electrolytes.

1. Introduction

The ability of solid electrolytes to conduct electric charge through ion transport is central to their use in devices such as fuel cells and solid-state lithium-ion batteries [1–4]. In both cases, solid electrolytes with high ionic conductivities are desirable. In fuel cells high conductivities allow lower operating temperatures, reducing running costs and increasing operating lifetimes. In solid-state batteries, high conductivities allow faster charging rates and higher power outputs. Ionic conductivities depend on a number of factors, including the crystal structure, the chemical composition and the concentration of mobile ions [5]. Developing a quantitative understanding of how these factors interact is key to developing high-conductivity solid electrolytes for use in high-performance electrochemical devices.

Solid electrolytes can be considered to comprise two distinct sets of ions: ‘fixed’ ions that vibrate about their crystallographic sites, and ‘mobile’ ions that can move through the system. The fixed ion positions define a network of diffusion pathways through which the mobile ions move. Solid electrolytes with crystal structures in common have diffusion networks that are topologically equivalent, while electrolytes with different crystal structures have topologically distinct diffusion paths. While much research into solid electrolytes focuses on understanding differences in ionic conductivities within specific structural families, a complementary question considers how differences in crystal structure, and hence diffusion network topology, affect ionic transport.

Diffusion pathway geometries are defined by crystal structure, and therefore are a microscopic property of specific materials. The performance of solid electrolytes in devices, however, is characterized by macroscopic transport coefficients: diffusion coefficients and ionic conductivities; which describe ensemble averages over all microscopic diffusion processes. Understanding the differences in ionic conductivity between solid electrolytes depends on resolving the quantitative relationships that link these two perspectives; in doing so, connecting the microscopic picture of specific ion-diffusion mechanisms to the macroscopic properties of long-ranged mass and charge transport.

In many solid electrolytes, the microscopic transport of ions can be approximated as a sequence of discrete ‘hops’ between distinct lattice sites.¹ If these hops are independent, every ion follows a random walk. The tracer diffusion coefficient, D*, and ionic conductivity, σ, can then be expressed in terms of the average hop rate per atom, ² via [8,9]

1.1

and

1.2

where a is the characteristic hop distance, C is the mobile-ion concentration and q is the charge of the mobile ions. Equations (1.1) and (1.2) can be combined to give the Nernst–Einstein relation, which relates D* and σ:

1.3

These three equations provide quantitative relationships between the hop rate, , tracer diffusion coefficient, D* and ionic conductivity, σ. Their derivation, however, depends on the assumption of independent hops, which holds only in the limit of very low carrier concentrations, or for fully non-interacting mobile ions [10].

Practical solid electrolytes typically have high carrier concentrations, and interparticle interactions can be significant. Under these conditions, individual hopping probabilities depend on the positions of nearby ions, and hops are no longer statistically independent. Instead, ion trajectories are correlated, and the system dynamics deviates from random walk behaviour [8,11–13]. Correlations between hops made by any single ion modify the relationship between average hop rate per atom, , and tracer diffusion coefficient, D*, which becomes

1.4

where f is a single-particle correlation factor that accounts for the deviations from random walk behaviour. Correlations between hops made by different ions modify the relationship between and σ, which becomes

1.5

where f_I is a collective or ‘physical’ correlation factor [7,10,14]. The relationship between σ and D* now differs from Nernst–Einstein behaviour (equation (1.3)) by the ratio of these correlation factors:

1.6

The inverse ratio f/f_I is commonly referred to as the Haven ratio, H_R [10,15].

Empirical relationships between microscopic hopping rates and macroscopic transport coefficients can be obtained, in principle, by combining experimental data for , D* and σ. Ion hopping rates can be measured in NMR or muon spin-relaxation experiments [16–21], diffusion coefficients obtained from tracer diffusion experiments [22], and ionic conductivities extracted from impedance spectroscopy [23,24]. Computational methods provide an increasingly useful complement to experimental studies of solid electrolytes. First principles calculations of vibrational partition functions and barrier heights along diffusion pathways can be used to obtain hopping rates ab initio [25,26]. Molecular dynamics simulations can be used to directly calculate diffusion coefficients and ionic conductivities [27]. Often, however, one or more of are unknown, and it is necessary to derive these from the other, known, properties. In principle, quantitative conversions between are possible via equations (1.4)–(1.6), provided the correlation factors f, f_I (and hence also H_R) are known.

For many simple crystal lattices, the correlation parameters have been calculated [10,28]. For more complex crystal structures, however, these parameters are often still unknown. A common approximation, therefore, is to assume that correlation effects can be neglected, which allows the simpler equations (1.1)–(1.3) to be used. This approximation is equivalent to assuming dilute-limit non-interacting behaviour. In solid electrolytes, where ionic motion exhibits strong correlation effects, however, this can introduce quantitative errors when processing data.

In this study, we report lattice-gas Monte Carlo simulation of ionic transport on the garnet lattice, performed to quantify correlation effects for this lattice geometry. The garnet lattice provides a model for diffusion pathways in the ‘lithium-garnets’, [29,30]. This family of solid lithium-ion electrolytes has attracted significant attention as candidate electrolytes for all-solid-state lithium-ion batteries [1,31–33]. The garnet crystal structure has an unusual three-dimensional network of lithium diffusion pathways, consisting of interlocking rings [34]. Each ring comprises twelve alternating tetrahedral and octahedral sites. Each tetrahedral site is coordinated to four octahedral sites, and each octahedral site is coordinated to two tetrahedral sites, with the tetrahedral sites acting as nodal points connecting adjacent rings (figure 1).

Aliovalent substitution of the M and M^′ cations allows the lithium stoichiometry to be tuned across a broad range. A theoretical lithium stoichiometry of x_Li=9 corresponds to a fully occupied lithium-site lattice. Ionic conductivities vary enormously as a function of x_Li, with σ increasing by approximately 10⁹ between Li₃Tb₃Te₂O₁₂ and Li_6.55La₃Zr₂Ga_0.15O₁₂ [1,30]. It remains an open question how the lithium diffusion coefficient and ionic conductivity vary with lithium stoichiometry. It is also not known to what extent the unusual diffusion pathway topology affects ionic transport. Resolving these questions is critical for the optimization of ionic conductivity for this family of materials.

Structural considerations and published data both suggest lithium garnets exhibit significant correlation effects. The low connectivity of the two-coordinate octahedral sites means blocking effects are expected to be considerable [34]. Short distances between neighbouring lattice sites of approximately 2.4 Å suggest strong lithium–lithium repulsion, expected to be particularly significant at high lithium stoichiometries [35–38]. The presence of two non-equivalent sets of lattice sites is also a factor. Non-interacting lithium ions would be expected to occupy octahedral and tetrahedral sites in a 2 : 1 ratio, reflecting the relative site populations. Neutron data, however, show that at low lithium content (x_Li=3) only tetrahedral sites are occupied [39], while at higher lithium content (x_Li=5→7) octahedral sites become preferentially occupied [30,37]. Experimental conductivities depend nonlinearly on x_Li [40], and deviate from ideal values predicted from muon-spin–spectroscopy hopping rates (via equation (1.2)) [20]. Further evidence for correlated transport in lithium garnets comes from computational studies. A variety of correlated diffusion processes have been observed in molecular dynamics simulations [41–44], and calculated diffusion coefficients and ionic conductivities show non-Nernst–Einstein behaviour (H_R<1) [45]. These results collectively indicate the existence of significant interactions, either between lithium ions or between these ions and the host lattice. The quantitative effects of correlation in lithium garnets, however, are not known, and consequently studies often assume uncorrelated behaviour when extrapolating between hop rates, diffusion coefficients and ionic conductivities [20,21,23,41,46–54].

Here, we present a computational study of these correlation effects, using lattice-gas kinetic Monte Carlo simulations of diffusion on a garnet lattice, across a range of model Hamiltonians. We calculate f and f_I as functions of lithium stoichiometry, first for a non-interacting volume-exclusion model,³ and then for models that include on-site single-particle energies and/or repulsive nearest-neighbour interactions. In addition to self- and collective-correlation factors, we present site occupation populations, diffusion coefficients and reduced ionic conductivities for this range of simulation models. Our results illustrate how different interactions contribute to non-ideal behaviour, and modify the relationships between particle hopping rate, diffusion coefficient and ionic conductivity.

We find that for non-interacting particles (volume exclusion only) single-particle correlation effects are more significant than for any previously studied three-dimensional lattice. This is attributed to the presence of two-coordinate octahedral sites, which produce correlation effects intermediate between typical three- and one-dimensional lattices. Including nearest-neighbour repulsion or on-site energy differences gives more complex single-particle correlation behaviour and introduces collective correlations. In particular, we find strong correlation effects at x_Li=3 (due to site energy differences) and x_Li=6 (due to nearest-neighbour repulsion). Both effects correspond to mobile particles ordering over the lattice, with associated sharp decreases in diffusion coefficients and ionic conductivities. By analysing our simulation data, we consider the question of tuning the mobile-ion stoichiometry to maximize the ionic conductivity. We show this does not have a straightforward answer, and the optimal stoichiometry is highly sensitive to the choice of interaction parameters. Finally, we discuss the practical implications of these results in the context of garnet-structured and other solid electrolytes.

2. Methods

Lattice-gas Monte Carlo simulations describe the diffusion of a set of mobile ions populating a host lattice, expressed as a graph of interconnected sites [56]. Every lattice site is either occupied or vacant, and during a simulation the mobile ions hop from site to site. These hops are randomly selected, with relative probabilities that satisfy the principle of detailed balance and represent the underlying model Hamiltonian. The simplest model considered here is a non-interacting volume-exclusion model [55]. Double occupancy of sites is forbidden, and allowed hops are all equally likely. The non-interacting model allows the pure geometric effect of the lattice to be evaluated, but neglects other interactions that may be important in experimental systems. We therefore also consider the effects of nearest-neighbour interactions between mobile ions, described by a nearest-neighbour repulsion energy, E_nn, and of interactions between single ions and the lattice, described by site-occupation energy differences between tetrahedral and octahedral sites, E_tet, E_oct. The energy of any configuration of occupied sites, , is given by

2.1

where is the number of occupied nearest neighbour sites for (occupied) site j. For interacting systems, the relative probability of hop i depends on the change in total energy if this hop was selected, ΔE_i, according to the scheme of Metropolis et al. [57]:

2.2

For our interacting systems, the change in energy for each candidate hop can depend on the change in the number of nearest-neighbour interactions and on the change in site-occupation energy when moving from a tetrahedral to octahedral site (or vice versa):

2.3

where ΔE_site=E_oct−E_tet. At each simulation step, one hop is randomly selected according to the set of relative probabilities. The corresponding ion is moved, and a new set of relative hop probabilities is generated for the following simulation step.

In the limit of a large number of hops, the tracer- and collective-correlation factors can be evaluated as

2.4

where 〈R²〉 is the mean-squared displacement of the mobile ions, N is the total number of hops during the simulation [25] and

2.5

where is the net displacement of all mobile particles. In both cases, the denominators correspond to the limiting behaviour for uncorrelated diffusion.

To allow time-dependent properties to be evaluated, such as average site occupations and transport coefficients, we perform our simulations within a rejection-free kinetic Monte Carlo scheme [58]. At each simulation step, k, the set of relative hop probabilities, , are converted to rates, , by scaling by a common prefactor ν^′=10¹³ s⁻¹. After selecting a hop, the simulation time is updated by , where Q_k is the ‘total rate’; and u is a uniform random number u∈(0,1].

Our lattice-gas kinetic Monte Carlo simulations were performed using the lattice_mc code [59]. Simulations were performed for an ideal cubic 2×2×2 garnet lattice, with 384 octahedral sites and 192 tetrahedral sites. The lattice-site coordinates were generated from the cubic high-temperature Li₇La₃Zr₂O₁₂(LLZO) structure,⁴ using the centres of the octahedra and tetrahedra defined by the oxide sublattice. In cubic LLZO, each octahedron available to lithium contains a ‘split’ pair of distorted 96h sites, separated by 0.81 Å . The construction used here considers each octahedron as a single ideal 48g site. The graph of diffusion pathways includes connections between nearest-neighbour sites only, i.e. all connections are between neighbouring tetrahedra–octahedra pairs. For each simulation, n_Li mobile ions are initially randomly distributed across the lattice sites. We perform 1000 simulation steps for equilibration, followed by 10 000 production steps.

For each set of model parameters, , simulations were performed across the full range of possible lithium stoichiometry. For a 2×2×2 garnet supercell, the maximum lithium content of x_Li=9 corresponds to n_Li=576. For each set of interaction parameters, data were collected as an average over 5000 independent trajectories.

3. Results

3.1. Non-interacting particles and geometric effects

We first examine the geometric effect of the garnet lattice by considering non-interacting particles, where any deviations from random walk behaviour are solely due to blocking effects. Figure 2 shows, as a function of x_Li, (panel a) the calculated self- and collective-correlation factors, f and f_I, (panel b) average tetrahedral and octahedral site occupations, n_tet and n_oct, (panel c) tracer and ‘jump’ diffusion coefficients, D* and D_J, respectively, and (panel d) a reduced ionic conductivity, σ^′ (equation (3.3)).

In the single-particle limit, x_Li→0, both correlation factors equal 1. There are no blocking effects, and particles follow a random walk. With increasing concentration, however, the single-particle diffusion deviates from random walk behaviour. The tracer correlation factor, f, decreases from f=1 in the single-particle limit to f=0.25 in the single-vacancy limit x_Li→9, showing approximately linear dependence on x_Li.⁵

The magnitude of the tracer correlation effect for different lattice geometries can be quantified by considering f in the limit of a single vacancy, f_v. Table 1presents f_v values previously calculated for simple three-dimensional lattices [60], and for a one-dimensional chain [7], alongside our result for the garnet lattice. The garnet lattice value of f_v=0.25 is less than for all previously studied three-dimensional lattices, and is a factor of two smaller than the next lowest (the diamond lattice). This indicates particularly strong site-blocking effects. For a general set of three-dimensional lattices, as the number of nearest neighbours of each lattice site, z, decreases, f_v also decreases, as site-blocking effects become more significant [28]. The garnet lattice has both four-coordinate (tetrahedral) and two-coordinate (octahedral) sites, and ion hopping follows an alternating tet→oct→tet→oct sequence. The calculated value of f_v=0.25 is halfway between the values for a four-coordinate three-dimensional diamond lattice (f_v=0.5) and for the two-coordinate one-dimensional chain (f_v=0) [7]. This suggests that the low value of f_v for the garnet lattice is a consequence of the low coordination of the lattice sites, in particular the local one-dimensional coordination at the octahedral sites, which act as bottlenecks for long-ranged diffusion.

For any non-interacting system, the hops made by different particles are uncorrelated, and f_I=1 for all x_Li; hence H_R=f. There are also no correlations between site occupations, and the mobile particles are randomly distributed over the available octahedral and tetrahedral sites, with a 2 : 1 occupation ratio that reflects the underlying lattice geometry.⁶

We also calculate three explicit measures of ionic transport in this system.⁷Figure 2c shows the tracer diffusion coefficient, D* (equation (1.4)) and the ‘jump diffusion coefficient’, D_J [5], calculated as

3.1

At a fixed temperature, D_J is proportional to the mobility and measures the ease with which the mobile particles collectively migrate. Both D* and D_J decrease monotonically from x_Li=0 to x_Li=9 (x=0→1), as progressively fewer vacancies are available to accommodate hopping ions. For the non-interacting system there are no correlations between hops made by different particles, and the jump diffusion coefficient is proportional to (1−x) (in the garnet lattice, x=1 corresponds to a stoichiometry of x_Li=9) [5,55]. The tracer diffusion coefficient, however, is affected by correlations between hops made by individual particles, and varies as D*∝(1−x)f.⁸ The ionic conductivity of a system depends on both the charge-carrier concentration, and the ionic mobility, which is proportional to D_J. We quantify the relative effect of carrier concentration on ionic conductivity by considering a reduced conductivity, σ^′,⁹ given by

3.3

For any non-interacting system, σ^′∝x(1−x), giving a maximum at x=0.5, corresponding to x_Li=4.5 in the garnet lattice (figure 2d).

3.2. Interacting particles

The conceptual simplicity of the non-interacting system makes it a useful starting point for understanding the factors affecting ionic transport in different lattices. In particular, purely geometric effects can be resolved. In real lithium-garnet materials, however, interactions exist between lithium ions, and between lithium ions and the host lattice, and these can significantly affect ion transport. Lithium ions are positively charged, and can be expected to experience mutual electrostatic repulsion. The different oxygen-coordination environments of the octahedral and tetrahedral sites can be expected to produce a preference for occupation by lithium at one site versus the other [63]. Within the lattice-gas Monte Carlo scheme, we consider these two factors by introducing, first, nearest-neighbour repulsion, and second, an octahedral versus tetrahedral site preference.

3.2.1. Nearest-neighbour repulsion

For lithium–lithium repulsion, we consider a simplified model with only nearest-neighbour repulsion. The energy of lithium at each site is now proportional to the number of occupied neighbouring sites, and individual hop probabilities depend on whether they increase or decrease the total number of nearest-neighbour pairs. Figure 3 presents results from simulations performed for E_nn=0.0–3.0 kT. Repulsive nearest-neighbour interactions disfavour simultaneous occupation of adjacent pairs of sites, which promotes ordering of particles into alternating occupied–vacant–occupied–vacant configurations. This ordering causes the single-particle correlation behaviour to deviate from that of the non-interacting system, and also introduces collective correlations between the mobile ions [10]. f and f_I both have their non-interacting values in the empty and fully occupied lattice limits: x→0 and x→1. In a lattice with only one crystallographic site, complete ordering would occur at half site occupancy, corresponding to x_Li=4.5 for the garnet lattice. We see that f and f_I approximately follow this trend (figure 3a,b), both decreasing at intermediate x_Li values as E_nn increases. Superimposed on this general shape, for larger E_nn values, both correlation factors sharply decrease at x_Li=6, i.e. two-thirds occupancy. Because f and f_I do not change uniformly as E_nn is increased, the Haven ratio H_R develops structure. Above x_Li=6, corresponding to stoichiometries of typical lithium-stuffed garnets, nearest-neighbour repulsion reduces H_R even further from the already low value for the non-interacting system.

The garnet lattice contains octahedral and tetrahedral sites in a 2:1 ratio. In the non-interacting system, the average site occupancies follow this ratio at all values of x_Li (figure 2b). Introducing repulsive nearest-neighbour interactions increases the probability that octahedra are occupied relative to tetrahedra. Because octahedral sites are two-coordinate, compared to the four-coordinate tetrahedral sites, occupying octahedral sites minimizes the number of unfavourable nearest-neighbour interactions. This effect is strongest at two-thirds site occupation (x_Li=6) where a sufficiently large E_nn drives the system into a fully ordered arrangement with all the octahedral sites filled and all the tetrahedral sites empty. In this fully ordered system, correlation effects are maximized: a single ion hopping from octahedron to tetrahedron is blocked from further forward motion, and must return to its starting position unless the blocking ion moves first, disrupting the local ordering.¹⁰ Diffusion is only possible for groups of particles undergoing highly correlated collective motion [44,64]. Both tracer diffusion and ionic conductivity are strongly reduced compared to their values in the non-interacting system. The collective correlation effects (f_I<1) are visible in the reduced conductivity, σ^′, which decreases relative to the non-interacting system across the full x_Li range, with a particularly strong decrease at x_Li=6.

3.2.2. Asymmetric site-occupation energies

In the non-interacting model, not only do mobile ions not interact with each other (excepting volume exclusion), but also there are no interactions between the mobile ions and the host lattice. Identifying a site as octahedral or tetrahedral only has relevance for defining the connectivity of the lattice graph. Mobile ions show an equal preference for octahedral and tetrahedral sites, with average occupations following a simple 2:1 ratio. This behaviour contrasts with experimental observations. Neutron data for lithium-garnets such as Li₃Y₃Te₂O₁₂reveal that at x_Li=3 the lithium ions exclusively occupy the tetrahedral sites [39].¹¹ This suggests that at relatively low lithium concentrations, there is an energetic penalty for occupying octahedral rather than tetrahedral sites.¹² We model this difference in site-occupation energies by including an on-site term ΔE_site=E_oct−E_tet. To investigate the effect of this ion–lattice interaction on ion dynamics and site occupations we performed a series of simulations for otherwise non-interacting particles, with ΔE_site=0–5 kT.

The effect of ion–lattice interactions qualitatively mirrors the effect of nearest-neighbour interactions (figure 4). Both single-particle and collective correlation factors are lower than their non-interacting values, average site occupancies deviate from those in the ideal system, and the reduced ionic conductivity decreases. Here, however, the strongest correlations emerge at x_Li≈3. As ΔE_siteincreases, tetrahedral sites are preferentially occupied with respect to octahedral sites, contrasting with the opposite behaviour observed for increasing E_nn. In the limit T→0, this again results in a fully ordered arrangement of ions, now with all the tetrahedral sites filled and all the octahedral sites empty. The Haven ratio, H_R, shows less variation compared to the non-interacting result, with only a small decrease for x_Li<3.

3.2.3. Combined site inequality and nearest-neighbour repulsion

In real lithium-garnet electrolytes, lithium ions can be expected to interact with both the host lattice and with each other. To explore the behaviour when both nearest-neighbour and site-occupation interactions are present, we performed simulations to map the {x_Li,ΔE_site,E_nn} parameter space. The data from these calculations are presented in figures 9–12. With both interactions present, the ion dynamics and site occupation statistics are more complex, with specific details that depend on the precise values of both interaction terms. The general features, however, are illustrated by considering the subset E_nn=ΔE_site (figure 5). The correlation factors, f and f_I, both show sharp decreases at x_Li=3 and at x_Li=6, in both cases corresponding to ordered arrangements of lithium ions throughout the lattice. As in the previous single-interaction models, the ordering at x_Li=3 corresponds to filled tetrahedra and empty octahedra (due to ΔE_site), and the ordering at x_Li=6 corresponds to filled octahedra and empty tetrahedra (due to E_nn). The average site occupation switches sharply from pure tetrahedral occupation to pure octahedral occupation in the range x_Li=3→6. The reduced ionic conductivity, σ^′, is depressed most strongly at lithium stoichiometries corresponding to the ordered arrangements of ions, again, mirroring the results for single interactions.

3.2.4. Tuning lithium stoichiometry to maximize ionic conductivity

One challenge regarding lithium-garnet solid electrolytes is the question of identifying specific compositions with high ionic conductivities. For garnets with stoichiometries Li_xA₃B₂O₁₂, the lithium content can be continuously varied by choosing appropriate A and B cations, or by substituting Li⁺ with small hypervalent cations such as Al³⁺ or Ga³⁺. Different lithium stoichiometries can exhibit very different ionic conductivities. For example, at x_Li=3 (e.g. Li₃Y₃Te₂O₁₂) room temperature conductivities are typically too low to measure [1,30,39] while, for x_Li≈6.5 (e.g. Li_6.55Ga_0.15La₃Zr₂O₁₂), conductivities as high as 1.3×10⁻³ S cm⁻¹ have been reported [66,67]. One strategy for identifying lithium garnets with high ionic conductivities is to consider whether there is an ‘optimal’ lithium stoichiometry that maximizes the ionic conductivity [35,48,65,68–73]. A conceptually related question concerns how the ionic conductivity depends on the distribution of lithium ions over tetrahedral and octahedral sites [29,35,71,74]. The lithium distribution is itself a function of the lithium stoichiometry, modulated by the interactions experienced by the lithium ions, as seen above for the model Hamiltonians including nearest-neighbour repulsion and site-occupation energies.

For a non-interacting lattice gas, the ionic conductivity varies with the mole fraction of mobile particles, x, as

3.4

The (1−x) term is a ‘blocking’ factor, due to volume exclusion [55]. The conductivity varies parabolically, as seen in the non-interacting system results presented above (figure 2d). For the lithium garnets this would give a maximum ionic conductivity at x_Li=4.5. In real systems, the mobile ions are subject to additional interactions that introduce collective correlations, and the variation in ionic conductivity with mole fraction of mobile particles becomes

3.5

Because f_I is itself a function of x, this gives non-trivial overall concentration dependence that cannot be described analytically. The concentration dependence of f_I is an emergent property of the specific interactions the lithium ions are subject to, which indicates that the mobile-ion concentration that maximizes the ionic conductivity in turn depends on the details of the lithium-ion interactions.

To explore this relationship in the model systems considered here, we can identify the maximum reduced ionic conductivity as a function of lithium stoichiometry; ; for each interaction parameter set . The resulting surface in parameter space is plotted in figure 6.

As suggested by equation (3.5), the value of x_Li that maximizes the ionic conductivity is strongly dependent on the interaction parameter values, due to their effect on the f_I(x_Li). For the non-interacting system (E_nn=0, ΔE_site=0.0), . For non-zero interaction parameters, however, ranges from less than 1.5 to greater than 7.5. Interestingly, the site-occupation energy difference has little effect in the limit of zero nearest-neighbour interactions. As E_nn increases, deviates from the non-interacting value. At low values of ΔE_site, increasing nearest-neighbour repulsion causes the optimal x_Li to decrease. This is associated with the strong suppression of collective ion transport close to x_Li=6 (cf. figures 3 and 12). At high values of ΔE_site, however, increasing nearest-neighbour repulsion causes the optimal x_Li to increase, reaching a maximum value of approximately 7.5. Under these conditions, the preference to occupy tetrahedral sites dominates, and ion transport rates are most strongly decreased close to x_Li=3 (cf. figures 4and 12).

4. Summary and discussion

By considering ionic transport in solid electrolytes as effected by particles undergoing random hops, simple analytical relationships can be derived that quantitatively connect microscopic hop rates to macroscopic transport coefficients (cf. equations (1.1)–(1.3)). In real solid electrolytes, these equations are exact only in the dilute limit. At moderate mobile-ion concentrations, ion hops are not independent; instead, they are correlated. The probability of a specific hop occurring depends on the particular arrangement of other nearby ions. Correlations between hops modify the quantitative relationships between hop rates and transport coefficients, with deviations from random walk behaviour expressed via the single-particle and collective correlation factors, fand f_I, respectively. Quantifying these correlation factors allows accurate conversions between microscopic (hopping rates) and macroscopic (tracer diffusion coefficients and ionic conductivities) transport data. These factors also provide information about the ionic transport process: the single-particle correlation factor quantifies the efficiency with which individual ions move through the electrolyte structure; the collective correlation factor provides an equivalent measure for the efficiency of mass or charge transport.

The simplest cause of correlation effects is volume exclusion, where occupied sites are unavailable to adjacent ions. This effect causes sequential hops of single particles to become correlated. The precise value of f depends on the concentration of mobile ions and the geometry of the host lattice. For this reason, lattice geometry can be key to understanding different behaviours between structural families of solid electrolytes. Explicit interactions between the mobile ions, or between the ions and the lattice, produce additional single-particle correlation effects. These interactions can also promote ordering of the mobile ions, which causes hops by different particles to become correlated. In real solid electrolytes, therefore, microscopic ionic transport depends on both lattice geometry and the nature of interactions acting on the mobile ions.

In this study, we have explored this behaviour for the garnet lattice, which provides a model for the lithium diffusion network in lithium-garnet solid electrolytes. From a theoretical perspective, this lattice possesses intriguing topological features. Previous theoretical and computational analyses of correlated ionic transport in crystalline lattices have considered only lattices where all sites are geometrically equivalent. The garnet lattice, however, contains both four-coordinate tetrahedral sites and two-coordinate octahedral sites, arranged in an open three-dimensional network of interconnected rings (cf. figure 1).

To study correlation effects in the garnet lattice, we have performed lattice-gas kinetic Monte Carlo simulations [59]. These consider the host structure as an idealized lattice, and describe ion interactions through simple model Hamiltonians. Rather than seeking an explicit description of a single material, as one might by using, e.g. first-principles or classical molecular dynamics [41,42,44], here we focus on understanding general behaviour as a function of lattice geometry and mobile-ion stoichiometry, and how this changes in response to conceptually simple, but physically motivated, microscopic interactions.

We find that, for the non-interacting (volume exclusion only) system, the single-particle correlation effects due to the lattice geometry are more significant than for any previously studied three-dimensional lattice (table 1). We propose that this is a consequence of the lattice containing two-coordinate octahedral sites, which act as bottlenecks to diffusion, producing correlation effects intermediate between those of simple three-dimensional and one-dimensional lattices.

Explicit interactions acting on the mobile ions (here, nearest-neighbour repulsion and site-occupation energy differences) produce stronger single-particle correlation effects with a complex variation with x_Li. These explicit interactions also promote ordering at set mobile-ion concentrations, which manifests as large collective correlation effects, with f_I→0 as T→0 [10]. The precise mobile-ion stoichiometry where ordering occurs depends on the lattice geometry and the explicit form of the interaction energy term. Ordering occurs for mobile-ion stoichiometries that are commensurate with the stoichiometry and symmetry of lattice sites, where ordering minimizes the ion-interaction energy. In the cases considered here, nearest-neighbour repulsion promotes ordering at x_Li=6, with all octahedral sites occupied and all tetrahedral sites vacant. A site-occupation energy that favours tetrahedral site occupation promotes ordering at x_Li=3, with all tetrahedral sites occupied and all octahedral sites vacant.

The ionic conductivity depends on the mobile-ion concentration directly, through the mole fractions of mobile ions and vacant sites, and indirectly, through the collective correlation factor (cf. equation (3.5)). Because the form of f_I(x_Li) depends on the interaction energy term, for interacting systems there is no simple expression that gives the mobile-ion concentration that maximizes the ionic conductivity. The simulations presented here show that is in fact very sensitive to the type and strength of mobile-ion interactions, with the ‘optimal’ lithium stoichiometry varying in the range x_Li=1.5–7.5 within the range of parameters we have considered. Even within the simplified models studied here, therefore, ionic transport on the garnet lattice exhibits correlation effects that are both more significant than predicted for simple three-dimensional lattices, and that show a complex dependence on mobile-ion stoichiometry.

The prediction that lithium-garnet solid electrolytes exhibit strong correlation effects is consistent with the observation of highly cooperative diffusion processes in first-principles simulations [41,42]. Because the quantitative correlation behaviour is sensitive to the mobile-ion interactions, this raises the question of how the interactions in real lithium-garnet electrolytes might map to the effective interactions considered here. This sensitivity also raises the possibility of tuning ionic conductivities through isovalent substitution within the host lattice. As an example, host lattices containing different cations will have different lattice parameters, and different distances between neighbouring tetrahedral and octahedral sites. This will modify both ΔE_site and E_nn, with consequential non-trivial effects on f_I, and hence on σ.

Although f and f_I are sensitive to the interaction parameters, their ratio, f/f_I=H_R, is less so. The calculated Haven ratios can therefore be used to improve the quantitative nature of conversions between tracer diffusion coefficients and ionic conductivities, via the modified Nernst–Einstein relation (equation (1.6)). Figure 7shows the calculated Haven ratios for all parameter sets considered in our study. Also plotted is the average Haven ratio across parameter sets as a function of lithium stoichiometry.

The sensitivity of f and f_I to the interaction details means that estimating these correlation factors for specific materials purely from this study is challenging. It is clear, however, that assuming that f=f_I=1 is likely to introduce quantitative errors when extrapolating between hop rates and tracer diffusion coefficients or ionic conductivities. One qualitative observation is that where ion interactions are present, the resulting correlation effects will increase as the temperature decreases. One consequence of this temperature-dependent correlation is that Arrhenius plots of tracer diffusion coefficients or ionic conductivities may not give straight lines. Instead, as temperatures tend to zero, and correlation effects become more significant, they are expected to curve downwards. Figure 8shows Arrhenius plots of relative ionic conductivities at x_Li=6, calculated for non-interacting particles, using equation (1.2), and for interacting particles subject to nearest-neighbour repulsion. In both cases, a microscopic activation energy of 0.3 eV is used. For the interacting case, f_I is interpolated from our simulation results. The nearest-neighbour repulsion energy is chosen as the Coulomb energy for two point charges occupying neighbouring sites, at a separation of 2.4 Å , using a typical garnet relative permittivity of ϵ_r=50 [75]. The Arrhenius plot for the data calculated assuming uncorrelated hopping gives a straight line, and a linear fit to obtain an ‘observed’ activation energy recovers the microscopic activation energy of 0.3 eV. The data calculated including the temperature-dependent collective correlation factor, however, fall below the first dataset—the collective correlation decreases the ionic conductivity relative to the ideal value—and this effect becomes more significant as the temperature decreases and the correlation effects strengthen. Fitting to the low-temperature regime (less than 1000 K) gives an ‘observed’ activation energy of 0.42 eV. The additional temperature dependence in the collective correlation means that the observed activation energy cannot be directly equated with the microscopic activation energy.¹³

One of the limitations of this study is that it uses a fixed predetermined lattice geometry. The ordering predicted at x_Li=3 and x_Li=6 occurs at lithium stoichiometries that are commensurate with the lattice symmetry and site ratios. In both cases, the ordered lithium configuration, with either tetrahedral or octahedral sites fully occupied and the alternate site fully vacant, preserves the lattice symmetry. In real materials lattice distortions are possible, and ordering of mobile ions can occur in concert with lattice symmetry breaking. In the lithium garnets, the prototypical example is the low-temperature tetragonal phase of Li₇La₃Zr2O₁₂ (LLZO) [34,78]. This material is cubic at high temperature (T>600 K), but at lower temperatures, it undergoes a tetragonal distortion, associated with the lithium ions ordering to occupy all the octahedral sites and one-third of the tetrahedral sites, accompanied by a decrease in ionic conductivity of two orders of magnitude. This is another example of ions ordering at low temperature, with low ionic conductivities as a consequence of the resulting strong correlation effects [44]. Again, this ordering occurs at a stoichiometry commensurate with the lattice symmetry. In the case of LLZO, ordering is promoted at x_Li=7 because the accompanying tetragonal distortion lowers the crystal symmetry; in each lattice ring, the six tetrahedral sites, equivalent by symmetry in the cubic lattice, become an inequivalent set of (2+4) paired sites.

Lithium ordering coupled to symmetry breaking has also been predicted at other lithium stoichiometries by Kozinsky et al. [79], who performed a group-theoretical analysis combined with first-principles energy calculations. Interestingly, this study predicted an ordered phase at x_Li=6 with a lower symmetry than the parent cubic phase, with octahedra and tetrahedra occupied in a 3 : 1 ratio, as well as ordered phases at other lithium stoichiometries, again accompanied by spontaneous symmetry breaking and lattice distortion. Because we have restricted our study to the ideal cubic garnet lattice, our results provide no information about possible alternate ordered phases that might be energetically favoured in distorted garnet lattices (e.g. at x_Li=6) or that might appear at stoichiometries where we do not predict ordering. A complete description of the order–disorder phase behaviour in lithium garnets would need to include not only the ideal cubic lattice, but also symmetry-broken lattices. Studying this behaviour within a lattice-gas Monte Carlo simulation scheme would require a more sophisticated approach than used here.

A second limitation is the use of the Metropolis scheme for calculating hopping probabilities (equation (2.2)). This approach considers hops to be barrierless, with hopping probabilities that depend only on the energy differences between initial and final configurations. In real materials, ion hopping is an activated process, and ions move across potential energy barriers. Expressions for hopping probabilities that take these barrier heights into account are expected to give more accurate kinetics, but require parameters that describe typical barrier heights, and how these are affected by the instantaneous local arrangements of mobile atoms [61,62]. For specific materials, these parameters can be derived from first-principles calculations [80–83]. For this study, we have focused on a broad description of the geometry effects in the garnet lattice. Including hopping barriers in a general scheme would significantly increase the dimensionality of the available parameter space, making a full analysis of lattice geometry effects impractical. It is apposite, however, to consider to what extent the results presented here, using the simpler Metropolis scheme, might differ from equivalent calculations that do account for hopping barrier effects. Addressing this question in the context of specific garnet materials will be the subject of a future study.

Interestingly, the behaviour described here, that in interacting systems, ordering is predicted at particular stoichiometries commensurate with the lattice symmetry, which manifests as strong correlation effects and greatly reduced transport coefficients, is qualitatively similar to results obtained for other lattice geometries using lattice-gas Monte Carlo models that do include barrier terms. For example, Murch and Thorn have modelled the effects of site-energy differences and nearest-neighbour repulsion in the two-dimensional honeycomb lattice, using a fixed transition barrier when deriving their hopping probabilities, and observed ordering and strong correlation effects at half occupancy [84–86]. The observation of qualitatively similar behaviour using models that do account for hopping barriers, albeit for different lattice geometries, suggests that these effects are not strongly dependent on the precise scheme. It should also be noted that any ordering of the mobile particles, which is the physical origin of these correlation effects, is independent of any transition barriers. The equilibrium distribution of particles depends only on the relative energies of different configurations, and is therefore exactly described (for a given Hamiltonian) by the Metropolis scheme.

A third consideration is that ion transport is assumed to be effected by a sequence of discrete hops made by individual ions. Although this is a good model for ionic transport in a large number of solid electrolytes, this is not always the case. In particular, so-called ‘superionic’ solid electrolytes exhibit diffusion mechanisms where ions move through highly concerted ‘liquid-like’ processes [87,88]. For solid lithium-ion electrolytes with particularly high ionic conductivities, such as those typically of interest for all-solid-state lithium-ion batteries, it is not known to what extent ion transport proceeds by concerted rather than single-ion diffusion mechanisms.¹⁴ In the case of the lithium garnets, data are limited and confined to cubic LLZO. Meier et al. performed a first-principles metadynamics study of cubic LLZO and identified a concerted diffusion process in their simulation trajectory [42], and a recent first-principles study by He et al. showed that concerted diffusion processes in this material can have lower potential energy barriers than single-ion hopping processes [89]. Support for single-ion hopping, however, comes from a study by Chen et al., who performed classical molecular dynamics simulations of LLZO [73]. By decomposing their simulation trajectories into sequences of single-ion hops, these authors showed that diffusion can be modelled as a Poisson process, which is a characteristic signature of an independent hopping process [64].

The question of contributions from concerted diffusion processes is not only pertinent to high-conductivity systems, but can also be important in ordered phases with low ionic conductivities. Under strong ordering of mobile ions, correlation effects may sufficiently impede ion transport by single-particle hopping that alternate concerted mechanisms become the dominant ion transport process [64]. This is believed to be the case for the low-temperature tetragonal phase of LLZO, with lithium transport effected by highly concerted motion of groups of ions moving around the lattice rings [44]. In the context of developing a theoretical framework that can quantitatively connect microscopic diffusion processes in solid electrolytes to macroscopic transport coefficients, a general treatment of concerted diffusion mechanisms remains an intriguing problem.

Data accessibility

Electronic supplementary material for this study is available as a GitHub repository, published under the CC-BY-SA-4.0 licence [76]. This repository contains (i) the complete dataset used to support the findings of this study, (ii) example scripts for running lattice_mc simulations on a garnet lattice and collating output data and (iii) a Jupyter notebook containing the code used to generate figures 2–12. The lattice_mc code is available under the MIT licence [59]. All data used to support the findings of this study are available as part of the electronic supplementary material dataset available at http://dx.doi.org/10.5281/zenodo.821870.

Competing interests

B.J.M. has no competing interests.

Funding statement

B.J.M. acknowledges support from the Royal Society (UF130329).

Footnotes

This article has been edited by the Royal Society of Chemistry, including the commissioning, peer review process and editorial aspects up to the point of acceptance.