First-principle DFT calculations were conducted using the Quantum-ESPRESSO software (gpu-enabled version 7.2), based on plane-wave DFT (Giannozzi et al., 2009; 2020). Ultrasoft pseudopotentials were utilized in all DFT calculations (Vanderbilt, 1990). For thermodynamic calculations, the PBE exchange-correlation functional was applied (Perdew et al., 1996), while structural parameters were calculated using the PBEsol exchange-correlation functional (Perdew et al., 2008). The expected improvement in describing structural parameters with PBEsol functional, as expected by its design (Perdew et al., 2008), was demonstrated in many of our previous studies (e.g., Blanca Romero et al., 2014; Connor et al., 2021; Kowalski et al., 2021; Ting and Kowalski, 2023). The supercell unit was set as Li₅₆La₂₄Zr₁₆O₉₆ for pure LLZO. Calculations were performed with a plane-wave cutoff energy of 50 Ry and a 4 × 4 × 4 Monkhorst-Pack k-point mesh (Monkhorst and Pack, 1976) to ensure convergence. Structural optimizations were carried out with convergence thresholds of 10^–5 for energy and 10^–4/a ₀ (where a ₀ is the Bohr radius) for forces.

The investigated systems were based on cubic-LLZO (c-LLZO) with space group Ia-3d and tetragonal-LLZO (t-LLZO) with space group I4₁/acd. In pure LLZO, both structures feature 8-fold coordinated La sites and 6-fold coordinated Zr sites. The Li sites differ between these two crystal structures. In c-LLZO, Li occupies tetrahedral 24d (Li1 site) and octahedral 96h (Li2 site), as depicted in Figure 2A. In t-LLZO, there are three Li sites: tetrahedral 8a (Li1), and two groups of octahedral sites, 16f (Li2) and 32g (Li3), shown in Figure 2B. In the initial DFT structure of c-LLZO, Li distribution is considered at positions 24d and 48g, with post-relaxation displacement to 96h (Figure 2C).

In our calculations, various dopants were selected to investigate their entropic effects, with compositions detailed in Table 1. Initially, partial Zr⁴⁺ was replaced with a A ⁵⁺ element, reducing the Li number from seven to six to induce more Li-vacancies and consequently higher configurational entropy. The number of Li atoms was kept constant in subsequent compositions. This enabled us to focus on the investigation of entropic contribution from the dopant concentration. In HEG2, alongside A ⁵⁺, partial Zr⁴⁺ was substituted with B ⁵⁺ and C ⁴⁺. HEG3 involved further replacement of partial La³⁺ with D ³⁺, gradually increasing configurational entropy.

Element distribution in multi-doped LLZO was categorized into Li-vacancy and metal-ion groups. Two methodologies were utilized: the Special Quasirandom Structure and a Coulomb-energy-based approach. SQS, executed via the ATAT package (van de Walle et al., 2013), distributes elements in a minimal-size supercell to resemble a fully random solid solution (Zunger et al., 1990). The Coulomb-energy-based (CE-based) method arranges atoms to minimize the Ewald Summation for Coulomb interactions, optimized through Monte Carlo sampling (Binninger et al., 2020). When species share the same oxidation number, we preferentially selected structures with higher symmetry. This strategic choice is underpinned by insights suggesting that structural symmetry plays a role in enhancing stability. Specifically, symmetry in a structure can lead to a reduction in the system’s degrees of freedom and overall energy, making the structure inherently more stable (Borisov et al., 2019). We calculated the 50 structures with the lowest Coulomb energies using DFT and chose the one with the minimum computed energy for our analysis.

With two groups of distribution and two methodologies, we systematically explored element distribution in four distinct ways: I. both groups using SQS; II. both using the CE-based method; III. one group via SQS and the other through CE-based method; IV. the inverse of arrangement III. These configurations are referred to as Li_SQS/M_SQS, Li_CE/M_CE, Li_SQS/M_CE, and Li_CE/M_SQS, respectively, where ‘Li’ refers to the Li-vacancy group and ‘M’ to the metal-ion group.

In calculating configurational entropy, Li-site vacancies are treated as distinct species. Notably, in c-LLZO’s cubic phase, despite Li occupying 96h positions, the placement of two Li atoms in adjacent sites is energetically unfavorable, leading to their treatment as a single site. This reduces the site count from 96h to 48g (Figure 2C). Thus, for configurational entropy calculations in the cubic phase, the total number of Li sites is considered to be 72, with 24 at Li1 and 48 at Li2 sites. The Li concentration in tetragonal and cubic phases is calculated based on the distribution of Li atoms over available sites.

The high-entropy garnets were synthesized by conventional solid-state reaction. The starting materials Li₂CO₃, La(OH)₃, ZrO₂, TiO₂, Ta₂O₅, Sb₂O₅, and Nd₂O₃ were mixed in stoichiometric amounts with 5 wt% Li₂CO₃ in excess by ball-milling in hexane. After dried, the powder was placed on a gold sheet inside an alumina crucible for calcination at 850°C for 12 h in air. After cooled down, the calcined powder was ball-milled again in hexane. The fine powder was pressed to form pellets at 13.7 MPa by using uniaxial hydraulic press, and the pellets were subsequently pressed at 150 MPa in an isostatic press. The obtained pellets were sandwiched between two gold sheets and sintered in an alumina crucible, and sintered in air at 1,050°C for 12 h. The sintered pellets were collected and stored in the Argone-filled glovebox.

The samples were characterized for their phase purities by X-ray diffraction (XRD) on a Bruker D4 Endeavor device (Bruker, Germany) using Cu Kα radiation equipped with a 1D detector LynxEye. The qualitative phase analyses were done using Panalytical’s HighScore software, and LeBail profile fitting for lattice parameters of the investigated samples were carried out using the software Jana2006 based on the structure model in the cubic space group of Ia-3d (ICSD 158372).

The Li-ion conductivities of HEGs were measured by electrochemical impedance spectroscopy (EIS). Prior to the EIS measurements, the HEGs pellets were polished on 400, 800, 1,500 and 4,000 grid sandpapers, and sandwiched by two lithium metal foils with a diameter of 6 mm. Such Li|HEG|Li symmetric assembly was then sealed in a pouch bag with two Ni bars as current collectors for the EIS measurements. EIS was measured in the frequency range from 10 MHz to 1 Hz with an electrical field perturbation of 50 mV using an impedance analyzer by Novocontrol Technologies, or in the frequency range from 1 MHz to 0.1 Hz with an electrical field perturbation of 10 mV using a frequency response analyzer by Solartron Analytical. A fitting of the impedance spectrum was conducted in the software ZView (Scribner) to obtain bulk conductivity.

Table 2 lists the studied dopants, their ionic radii (Shannon, 1976) and reported doping positions from prior studies. Dopant selection was based on the similarity of their ionic radii to either Zr⁴⁺ or La³⁺ and their doping positions. The formation Enthalpies from oxides (ΔH_f,ox) of the considered compounds were calculated using total DFT energies of reactants and products. We then selected those with the lowest ΔH_f,ox for further evaluation. Atomistic structures of products in the cubic phase were generated using SQS. For HEG1, HEG2, and HEG3, ΔH_f,ox was determined by the energy difference between reactants and products, as described in the following reaction: 3 L i 2 O + 3 - x 4 La 2 O 3 + x 2 D 2 O 3 + 2 - y - z - w ZrO 2 + y 2 A 2 O 5 + z 2 B 2 O 5 + w C O 2 → Li 6 La 3 − x D x Zr 2 − y − z − w A y B z C w O 12 The results in Figure 3 show that Sb has the lowest formation enthalpy among the elements considered for element A in Li₆La₃ZrAO₁₂, leading to the identification of HEG1 as Li₆La₃ZrSbO₁₂. For HEG2, Ta/Ti combination was by 5 kJ/mol slightly more favorable compared to the Nb/Ti combination. Therefore, HEG2 was considered as Li₆La₃Zr_0.5Sb_0.5Ta_0.5Ti_0.5O₁₂. Nd was chosen as the dopant D in HEG3 based on its thermodynamic stability. Notably, Sb-doped LLZO (HEG1) demonstrated the lowest ΔH_f,ox relative to HEG2 and HEG3. A subtle change in ΔH_f,ox was observed when transitioning from three-element-doped HEG2 to four-element-doped HEG3. This trend suggests that careful dopant selection leads to minor changes in ΔH _f,ox , but enhancing the configurational entropy (ΔS _conf ) with the addition of more dopant types.

Following the computational predictions, we proceeded with the synthesis of the compounds. HEG1, with composition Li₆La₃ZrSbO₁₂, was successfully synthesized in a cubic phase. However, challenges emerged during the electrochemical performance assessment of HEG2. A potential issue identified was the high titanium content in HEG2, which might cause a reduction of Ti⁴⁺ to Ti³⁺, leading to a transformation of the garnet from an electrical insulator to a conductor. To address this, we modified the composition of HEG2, reducing the Ti content from Li₆La₃Zr_0.5Sb_0.5Ta_0.5Ti_0.5O₁₂ to Li₆La₃Zr_0.75Sb_0.5Ta_0.5Ti_0.25O₁₂ and correspondingly adjusting the Zr content. For HEG3, we maintained the mole fraction in the Zr-sublattice and introduced Nd doping to the La-sublattice. This resulted in the Li₆La_1.5Nd_1.5Zr_0.75Sb_0.5Ta_0.5Ti_0.25O₁₂ model compound.

The distribution of elements in the cubic phase was modelled in four distinct ways, as described in Section 2.1: Li_SQS/M _SQS, Li_CE/M _CE, Li_SQS/M _CE, and Li_CE/M _SQS. The computed entropies illustrated in Figure 4 provide insights into the relative stability and phase preference across these different distribution methods. Contrary to our initial expectations, the Li_CE/M _CE atomic distribution resulted in the highest computed energies (enthalpies) for all three compounds, indicating the least stability of such arrangement. As opposite, the Li_SQS/M _SQS method resulted in the lowest energy, suggesting the fully disordered structure as the most stable one for each of the considered compounds. We note that this is in contrast to a common expectation of disordered compounds having higher formation enthalpies than ordered compounds with identical stoichiometry (e.g., Kowalski, 2020). After structural relaxation, we assessed lattice parameters and the position of Li to determine if the structures were in cubic or tetragonal phases. We emphasize the importance of examine the Li position, since the lattice parameters for relaxed structures can vary even in the cubic phase. We note that Li-vacancy distribution could influence phase stability, particularly when the number of dopant elements is low. For example, the relaxation of HEG1 structure that was determined from electrostatic calculations (Li_CE/M _CE) transitioned from cubic to tetragonal phase, while HEG2 and HEG3 maintained their cubic phase.

Analysis of Li distribution before and after structural relaxation shows that the site occupancies and displacement degrees depend on the distribution methods, as shown in Table 3. The CE-based method tended to fill the Li1 site and distribute remaining Li atoms to the Li2 site, contrasting with the SQS method where occupancies of Li1 and Li2 sites were identical. After relaxation, some Li atoms migrated from tetrahedrally coordinated Li1 to octahedrally coordinated Li2 sites. The displacement, particularly at Li1 sites, appear to be influenced by the occupancy after relaxation. For each compound, the higher the occupancy at the Li1 site, the greater the displacement degree observed, as reported in Table 3. This relationship was not found at Li2 sites. These findings suggest that the occupation of the Li1 site may influence structural displacement, affecting other properties such as ionic conductivity, activation energy, and diffusion paths. This appears to contrast with previous reports, which indicate that variations in Li + ion arrangements do not significantly impact DFT results (Rettenwander et al., 2014). Our results, however, underscore the importance of Li distribution in not only influencing phase stability but also other intrinsic material properties. The choice of atom distribution method, therefore, is not trivial and demands careful consideration.

We further investigate the influence of Li- and M-group atomic distribution on the thermodynamic stability of the phases. We uncover a nuanced relationship between changes in the atomic distributions of Li- and M-groups and their thermodynamic stability impacts. Figure 5 illustrates the enthalpy variation from different sublattice distribution methods for HEG1, HEG2, and HEG3 compounds. Specifically, we find that variation of energy is more sensitive to the M-group’s atomic distribution than the Li-group’s arrangement. This phenomenon is evident in the energy difference, where shifts in the M-group distribution yielded larger variances (Figure 5). For instance, the energy difference is more pronounced in HEG1 when the M-group atoms change from an SQS to a CE distribution, compared to the same change in the Li-group. The result suggests that the structural stability and electronic characteristics of the compounds are markedly more sensitive to the configuration of the M-group atoms than of the Li-group. This sensitivity could be attributed to the higher oxidation state of M ions (e.g., Zr⁴⁺, Ta⁵⁺) compared to Li⁺, making the electrostatic effects more impactful on enthalpy. These findings underscore the importance of the atomic configuration of dopant in modulating the compound’s energy profile and emphasize the importance of targeted distribution method selection to optimize material properties.

To investigate the influence of Coulomb interactions on the structural stability of high entropy garnets, the Ewald summation was performed for both initial and relaxed states of HEG1, HEG2, and HEG3. Figure 6 shows the obtained significant variations in Coulomb energy. These variations were more pronounced in structures with increased number of doping elements, particularly for those configured with Li_SQS/M _SQS. The complexity introduced by increasing number of dopants contributed to this trend, highlighted by the increased displacement within the tetrahedral Li1 site (details in Table 3), with displacement degrees of 0.133, 0.199, and 0.214 for the respective structures. An example that underscores this observation is the comparison between the relaxed structures of HEG1 and HEG3. Despite identical Li occupancies, the larger displacement degree in HEG3 — a consequence of the greater number of doping elements—especially at the Li1 site, leads to extended interatomic distances. This, in turn, results in a more substantial impact on electrostatic energy. Hence, it is evident that increased displacement of Li atoms significantly influencing the electrostatic energy.

A direct correlation between Coulomb and total DFT energies for the unrelaxed structures of high entropy garnet compounds is found, which is in line with the findings presented in Binninger et al. (2020). However, this correlation diminishes after structural geometry relaxation. For the relaxed structures across the three high entropy garnet compounds, the SQS structure exhibited the lowest computed enthalpies, whereas the CE-based model resulted in the highest computed enthalpies. In contrast, results from Ewald summation revealed an opposing trend, showing that both before and after structural relaxation, the CE-based model consistently yielded the lowest Coulombic energies. The inconsistency between DFT and Ewald summation results is likely attributable to changes in interatomic distances and electronic structures that occur during the relaxation process. These changes facilitate charge redistribution and electron delocalization, which in turn diminish the applicability of the localized point-charge model for the Coulomb energy in the relaxed state. Thus, our findings suggest that while CE-based atomic distributions effectively identify atomic configurations with the lowest energy dominated by Coulomb interactions, their predictive accuracy for determining the lowest-energy configurations decreases in materials with higher covalency and structural complexity, such as high-entropy garnet LLZO. In these materials, the CE-based methods are less reliable for determining the most stable, low-energy configurations.

In order to gain a deeper understanding of the thermodynamic properties of high-entropy garnet compounds HEG1, HEG2, and HEG3, we performed calculations of their configurational entropy and compared these findings with their enthalpy values. Figure 7 shows the enthalpy differences between the tetragonal and cubic phases of these compounds. For pure LLZO, the enthalpy of the tetragonal phase is lower by ∼ 6.5 kJ/mol compared to the cubic phase. This observation is in line with the computational findings reported by Bernstein et al. (2012). Intriguingly, upon doping LLZO with additional elements, a change in phase stability is found; the cubic phase becomes progressively more stable than the tetragonal phase as the number of doping elements increases, as shown from HEG1 to HEG3. While previous studies have primarily emphasized the role of lithium vacancies in stabilizing the cubic phase of LLZO (Thompson et al., 2014; Shin et al., 2015), our findings contribute a novel perspective by focusing on the number of doping elements. We have established, for the first time through DFT calculations, that increasing the variety of doping elements, rather than altering lithium concentrations, significantly enhances the cubic phase’s stability. This finding underscores the impact of doping on the thermodynamic stability of these garnet compounds, indicating a preference for the cubic phase structure with increased number of doping elements.

Afterwards, we compare the configurational entropies of the garnet compounds, as detailed in Table 4. For the fully ordered LLZO in its tetragonal phase, the configurational entropy ( Δ S c o n f t e t ) is zero. In contrast, the disordered cubic phase exhibits a configurational entropy of ∼ 30 J/mol/K per formula unit. This significant difference in configurational entropy allows us to estimate the phase transition temperature from the tetragonal to cubic phase. This is done by assuming that ΔG = 0, so T = ΔH/ΔS (Li et al., 2015). Assuming that the total entropy change (ΔS) is solely due to the configurational entropy (ΔS_conf) and neglecting the other entropy contributions (e.g., vibrational and electronic), the calculated temperature required to overcome the enthalpy difference (ΔH_cubic−tet) of 6.5 kJ/mol is ∼ 500 °C. However, we note that the experimentally reported phase transition temperature is around 650 °C (Larraz et al., 2013; Matsui et al., 2013). This discrepancy between the calculated and experimental values may be attributed to the neglected entropy contributions and/or not accurately accounting for the exchange-correlation electronic effects in DFT. Li et al. (2015) demonstrated that applying standard DFT method results in underestimation of disordering temperatures for a series of pyrochlore compounds. Additionally, the theoretical value of ΔS_conf represents the maximum possible value of configurational entropy; in reality, this value is expected to be lower due to the partial ordering preferences within the structure (Navrotsky, 2010; Finkeldei et al., 2017; Kowalski, 2020).

The analysis further reveals that introducing Li vacancies and increasing the number of dopants leads to an increase in the configurational entropy. For the compounds HEG1, HEG2, and HEG3, regardless of the variety and concentration of dopants, a consistent difference in configurational entropy—approximately 13 J/mol/K per formula unit—is observed between the cubic and tetragonal phases. This difference is mainly due to the greater number of possible Li sites in the cubic phase compared to the tetragonal phase, as reflected in Eq. 2. However, it is important to note that this entropy difference is relatively minor when compared to the enthalpy difference between the two phases. These observations suggest that while configurational entropy does play a role in stabilizing the cubic phase, enthalpy is the main driving force in the formation of this phase. We note, however, that entropy and enthalpy of a partially or fully disordered materials are correlated (e.g. Kowalski, 2020). This finding emphasizes the significant role of enthalpy in dictating the phase stability in high-entropy garnet compounds, particularly under varying conditions of dopant concentrations and structural configurations.

It has been suggested that the ideal Li content to form a conductive garnet with cubic structure is 5 < Li < 6.6 per formula unit (Bernstein et al., 2012). This can be seen, for example, in Li_{7– 3x–y}M _x La₃Zr_2−y B _y O₁₂ (M = Al, Ga, B = Ta) (Allen et al., 2012; Thompson et al., 2014; Baklanova et al., 2018). At Li contents larger than 6.6, the cubic structure undergoes a reduction of symmetry to a tetragonal polymorph and the ionic conductivity decreases. Our result is aligned with this observation. For Li content of 7, the stable structure appears to be tetragonal and as Li content reduces to 6, cubic phase is favorable. Jung et al. (2022) discusses that the cubic phase stabilization in the garnet with multiple dopants is presumably caused by an entropy effect rather than the enthalpy effect for a constant lithium content of 7.0 in Li₇La₃M₂O₁₂ (M = Zr, Hf, Sn, Sc, Ta). In our calculation, we observe a different scenario. The main stability driver is enthalpy, and configurational entropy provides extra stabilization of cubic phase.

The high-entropy garnet compounds HEG1, HEG2, and HEG3 were synthesized utilizing conventional solid-state reaction methodologies. The X-ray diffraction (XRD) patterns of HEG1, HEG2 and HEG3 are shown in Supplementary Figure S1. Additionally, Thermogravimetric Analysis/Differential Thermal Analysis (TGA/DTA) was performed to examine the thermal stability of HEG2, confirming the thermal stability of the garnet material (see Supplementary Figure S2). The experimentally measured lattice parameters were in close agreement with theoretical predictions, as reported in Table 5, underscoring the effectiveness of the PBEsol functional for accurate prediction of structural parameters. This alignment between experimental and computational results not only validates the synthesis process but also reinforces the reliability of the computational methods employed in this study.

Subsequent to synthesis, we conducted a comprehensive analysis of the electrochemical properties of these compounds, focusing on bulk Li-ion conductivity. The Electrochemical Impedance Spectroscopy (EIS) spectra of HEG1, HEG2, and HEG3 are presented in Supplementary Figure S3. The results, illustrated in Figure 8, revealed that HEG2 exhibited superior ionic conductivity compared to its counterparts. In contrast, HEG3, despite exhibiting the most significant displacement indicative of higher local disorder (see Table 3), displayed the lowest conductivity. This outcome challenges the previously suggested hypothesis by Zeng et al. (2022) regarding the enhancement of ionic conductivity via local disorder and site percolation in high-entropy materials. The minimal influence of entropy and the amount of doping elements on ionic conductivity suggests that these parameters may not be the primary drivers of conductivity in high-entropy garnet compounds. Detailed experimental data encompassing the influence of cation doping on the microstructure and additional electrochemical characterizations are presented in Ye et al. (2022).

Our findings indicate a more intricate interplay between local structural disorder, site percolation, and ionic conductivity than previously anticipated. This complexity warrants a deeper investigation. The insights from this analysis not only contribute to our understanding of the electrochemical behavior of high-entropy garnets but also open avenues for future research, particularly in optimizing these materials for advanced energy storage applications.