Equilibration of ion distribution at polymer/ceramic interfaces

Abstract

Composite solid electrolytes (CSEs), involving ceramic and polymer electrolytes, are a promising candidate for use in solid-state batteries. These non-homogeneous electrolyte systems have altered ion distributions compared to their respective bulk phases. Due to high energy barriers for ion exchange between phases, achieving interfacial ion equilibration can be an issue when classical molecular dynamics (MD) techniques are employed at temperatures relevant for battery applications (e.g. 300–360 K), as the accessible timescales are limited to <μs. In this work, we present MD simulations at elevated temperatures (400–700 K) to reach equilibrium ion distributions in CSEs and evaluate whether these high-temperature simulations can also provide insight into the structural and dynamic properties at lower temperatures. Specifically, we selected Li₇La₃Zr₂O₁₂ (LLZO) as the ceramic and lithium bis(trifluoromethanesulfonyl)imide (LiTFSI) dissolved in polyethylene oxide (PEO) as the polymer electrolyte. We demonstrate that, at temperatures equal to or above 600 K, the interfacial structures reach a state of equilibrium that is essentially independent of temperature, involving a significant Li⁺ transfer to the ceramics relative to the bulk phase. Below 600 K, the Li⁺ distribution is not equilibrated. Furthermore, even below 600 K, the TFSI⁻ distribution within the polymer phase follows the Li⁺ distribution in the two phases. This results in apparent equilibrium ion distributions at temperatures below 600 K. From the equilibrium Li⁺ free energy curve at 700 K, together with the transition rates at 600 K and 700 K, we can estimate the free energy barrier for Li⁺ entering the LLZO phase at lower temperatures, and thus predict the expected hopping rates at experimentally relevant temperatures. The results suggest that simulations at 400 K or lower require initial high-temperature simulations to achieve an appropriate lithium partitioning. The resulting structure can then be used for simulations at lower temperatures.

---

1 Introduction

Composite polymer electrolytes (CPEs) consist of blends of two components: ceramic particles and salt solutions in polymers, respectively. These materials are promising for utilization in different solid-state battery configurations.¹ While having been explored for decades, and often showing improved mechanical and ion conducting properties as compared to the corresponding neat (ceramic-free) solid polymer electrolytes (SPEs),² the area has been rejuvenated in recent years when ceramic particles with intrinsic ionic conductivity – i.e., active ceramic electrolytes – have been used for CPEs.³ Numerous such ceramic ionic conductors have been explored, encompassing oxides, sulfides, glasses, etc.⁴ The goal is then to combine the most useful properties of ceramic electrolytes – primarily high bulk ionic conductivity – with those of SPEs, i.e. mechanical flexibility and surface adhesion. However, the overall outcome of these studies is still inconclusive: large uncertainties exist in whether or not the inclusion of ceramic particles actually contributes to improvements, and also what the ionic conduction mechanisms actually are. It has been questioned if the ceramic particles are at all active in these CPE systems, i.e. through contribution to the overall ion transport by their bulk conductivity.⁵ The variety of systems studied is also large, and many system properties can be varied for the particles (type, size, loading), the polymer (type, structure, molecular weight) and the salt (type, concentration), just to mention a few variables. This intrinsic diversity, together with the general difficulty of studying ionic transport in composite solid-state systems using most experimental techniques, makes it difficult to draw affirmative conclusions. These are strong reasons for employing modeling techniques to help resolving these issues and obtain a deeper microscopic understanding.

Considering the comparatively long time-scales involved in mass transport processes in both polymers and ceramics, and that the atomic-level structure-dynamic properties dictate the involved mechanisms, classical molecular dynamics (MD) methodology appears as the method of choice for studying these materials.⁶ While not capturing redox activity, but which is of less concern in the involved systems, the technique can supply production runs of hundreds of nanoseconds, thereby capturing the diffusive regime for the ionic species. MD has been employed in numerous studies of both polymer electrolyte^7–10 and ceramic electrolyte^11–13 systems, and emerging work is being performed on CPEs.^14–19 In this latter case, however, shortcomings appear, partly due to the MD technique itself, but also due to the complexity of the system. First, atomic-scale MD, using classical force fields, is adopted for systems of some hundreds of thousand atoms. This is far from being able to capture the bulk properties of these materials, which consist of many nano- or micron-sized particles embedded in a polymer matrix. Instead, simulations are limited to capturing effects of the ceramic and polymer in the very interfacial region between these domains, generally <10 nm thick layers of the different materials. Second, as is well known, the surface chemistry of ceramic particles is often non-uniform, and the materials are often partially covered by adsorbed species or thin interfacial layers. The surface symmetry can also vary substantially, with several ceramic surfaces being exposed to the polymer bulk. Furthermore, the polymers used in these MD simulations are generally short-chain analogues of the high-molecular-weight counterparts used in experimental studies. This will influence chain ordering and dynamics. As has been shown in MD work,¹⁷ these uncertainties give rise to fundamentally different results when studying the structural properties at the ceramic–polymer interfaces. Third, the force fields employed are themselves generally developed for more simplistic systems, and merging them in studies of more complex systems should require sophisticated tailoring and proper bench-marking.⁹ This is, however, often not done.

In this study, mass transport and ionic distribution in CPEs are in focus. Whereas both of these phases are cation conductive, anions can only exist in the polymer phase or at the ceramic surfaces. Both cations and anions can be characterized by a certain equilibrium distribution in their respective bulk systems. Initially, when the polymer and ceramic are combined, anions and cations reorganize their positions to reach an equilibrium distribution in the new two-phase system, taking the initial difference in chemical potential between these phases into account. Anticipating that the inter-phase equilibration is much slower than the intra-phase equilibrium, one can define three equilibrium processes; (i) Li⁺ and TFSI⁻ distribution equilibration in the polymer phase, (ii) Li⁺ distribution in the ceramic phase, and (iii) the general Li⁺ distribution between the polymer–ceramic phases. In other words, this gives rise to the polymer and ceramic intra-phase ion distributions, respectively, as well as an inter-phase Li⁺ distribution in the system. The intra-phase and inter-phase equilibrium processes can be independent from each other, or coupled.

From the general challenges in modelling polymer–ceramic system, this work sets out to address two of them: the limit of accessible MD simulation time and the possibility to achieve faster equilibration. Nowadays, simulations usually reach a few hundred nanoseconds, which for some systems may be too short to reach the equilibrium state due to a large maximum relevant free energy barrier, for example between the polymer and ceramic phase. One way to circumvent this limitation is to perform MD simulations at elevated temperatures. Thereby, higher energy states in the free energy landscape can be explored, and the dynamics can be sped up without compromising the physical and structural properties.

The most relevant physical insights about ion transport in CPEs should be acquired at temperatures ranging from 300 to 360 K, equivalent to operating conditions for solid-state batteries. In practice, due to the limitations of MD simulations, ion exchange between phases in non-homogeneous systems at room temperature is rare, and the ions remain within the phase where they were initially inserted. Alternatively, phase exchange can occur sporadically, together with a gradual accumulation and depletion of ions in certain regions of the system. Equilibrium is not reached in any of these scenarios. This raises questions if it is possible to study ion transport in the CPE at, e.g., 400 K? And if not, can higher temperature simulations, e.g. at 700 K, be used as a strategy to gain physical insight about the system at lower temperatures?

Guided by these questions, we here aim to explore the possibilities of achieving robust information on structure-dynamic properties in CPEs as a prototype of a non-homogeneous system. Specifically, we focus on the LLZO–PEO_x LiTFSI system, i.e. particles of Li₇La₃Zr₂O₁₂ (LLZO) inserted into a polymer electrolyte composed of lithium bis(trifluoromethanesulfonyl)imide (LiTFSI) dissolved in polyethylene oxide (PEO). This is, arguably, the most well-studied CPE system in scientific literature. The study thereby targets making insights on how a robust equilibration can be carried out for MD simulations of CPEs, and which parameters and phenomena that can be appropriately determined from the extracted MD data. For the specific analysis, the free energy landscape turns out to be a viable observable. This, in turn, lays a solid foundation for further computational exploration of CPE materials.

2 Methods

2.1 Interface model details

The general features of the atomistic model used for PEO:LiTFSI/LLZO has been described in previous work,¹⁷ but with some additional adjustments as indicated below. The model consists of a LLZO slab, on both sides covered by LiTFSI in PEO. The system was created by first relaxing layer by layer. The general MD settings, including the system size, are largely taken from our previous work on this system.¹⁷ A snapshot, including the system size, is shown in Fig. 1.

Fig. 1 Snapshot of the simulation box of the system at 700 K at the beginning of the NVT run, including the dimensions of the box. Li⁺ are shown in orange, TFSI⁻ in purple, the polymer in green and the atoms from the LLZO excluding the Li⁺ in grey. The snapshot was produced using the software VMD.³³

Ceramic electrolyte. The LLZO slab was based on the cubic LLZO structure, which was taken from the ICSD crystallography database.²⁰ Using the python module Pymatgen “surface”, a slab with the dimensions of 3 × 3 × 3 unit cells, a (100) orientation and a c-shift of 0.9 was created.^21,22 Due to the asymmetry of the crystal structure, the generated slab possesses two different surfaces. Based on previous work¹⁷ the side that is initially terminated by oxygen atoms is defined as the LEFT side. The side terminated by lithium ions is referred to as the RIGHT side. Based on neutron-diffraction experiments,²³ in the slab lithium atoms were distributed randomly among the two different coordination sites which are present in LLZO: the 24d and the 96h sites using the respective occupation numbers, 0.542 and 0.448.

Polymer electrolyte. The LiTFSI-doped PEO polymer structure was built using an Enhanced Monte Carlo (EMC) code.²⁴ The z-axis was chosen to be the axis perpendicular to the interface of the ceramic material. The simulation box contained 45 symmetric polymer chains, all terminated with methyl groups and each being 20 monomers long. One single salt concentration was investigated, 1 : 10 (Li⁺ : EO). The EMC code allowed simultaneous creation of the mixture consisting of 90 LiTFSI ion pairs and the PEO chains.

Force fields. Interactions between the LLZO atoms were described by Buckingham and Coulomb potentials with parameters determined according to Jalem et al.¹² The PEO polymer was described by bonded and non-bonded (Lennard-Jones (LJ) and Coulomb potentials) interactions with the OPLS-AA²⁵ force field as defined in the EMC²⁴ software. The bonded and LJ parameters for the TFSI⁻ anions were adopted from Lopes et al.²⁶ It is worth noting that Li⁺ in LLZO force field possesses a partial charge equal to 0.7.¹² For that reason all Li⁺ and TFSI⁻ charges were scaled by 0.7. In fact, the charge scaling approach to simulate such polymer systems has been validated elsewhere,⁹ and is beneficial because it accounts for the polarization effect. The bond, angle and improper interactions were described harmonically, while the dihedral interactions were set as multi/harmonic.

To describe cross-interactions between the polymer and ceramic phases, the same method as in previous work on the same system was employed.¹⁷ Thereby, the Universal Force Field (UFF) was used to assign additional LJ parameters for the LLZO atoms. Since all atoms had associated LJ parameters, Lorenz–Bertheloth rules were used for the poly-nuclear interactions.

Electrostatic interactions were modeled by a standard Ewald summation as implemented in LAMMPS, and fixed partial charges were assigned to all atoms in the system. The cutoff for all non-bonded interactions was set to 12 Å.

2.2 Simulation details

All MD simulations were performed using the LAMMPS package.²⁷ The full simulation details can be found in the Fig. S1, including a work flow chart and examples of compression and decompression schemes for system equilibration.

General MD settings. Periodic boundary conditions were applied in all three dimensions. The timestep was set to 2 fs and the mass of the hydrogen atoms was changed to the mass of a deuterium atom for faster equilibration. A Nosé–Hoover thermostat and barostat, as defined in LAMMPS, were used to control temperature and pressure. The LLZO slab and the PEO phase were equilibrated separately because of the much longer equilibration time of the polymer phase compared to the ceramic phase.

1. LLZO slab equilibration. First, a minimization was performed to avoid initial Li⁺–Li⁺ overlap in the crystal. Next, a NVT simulation with a gradual increase of temperature from 100 K to 800 K was performed. The final equilibration run was an NpT simulation at 800 K and 1 atm.

2. PEO wall equilibration. The polymer was inserted above one flat wall, made of equally spaced dummy atoms interacting through LJ potentials with parameters (ε = 0.05 kcal mol⁻¹, σ = 2.193 Å). This was performed to restrict the polymer chains from crossing the boundary along the z-axis, to facilitate the following assembly with the LLZO slab. In order to equilibrate the polymer system confined by the wall, a procedure inspired by a 21 step equilibration process according to Larsen et al.,²⁸ was followed. Here, a series of simulations including compression and decompression were performed (multiple NpT and NVT ensemble simulations at varying temperatures and pressures). The snapshot at 10 ns of the final NpT simulation was extracted and the wall dummy atoms were removed.

The equilibrated ceramic slab and the salt-doped polymer phases were combined using the Moltemplate package.²⁹ To equilibrate the assembled polymer/ceramic system, two steps were performed:

1. PEO phase equilibration with fixed LLZO phase. This step was performed to allow PEO and TFSI⁻ adjustment to the LLZO slab without affecting the surfaces themselves, thereby avoiding excessive strain at the interface. This step involved a series of compression and decompression steps (multiple NpT and NVT ensemble simulations) employing increasing pressure and temperature conditions; see more details in the Fig. S1. The final snapshot in each step was extracted to start the following step, as described below.

2. All-atom equilibration. This step involved compression and decompression (using multiple NpT and NVT ensemble simulations) with temperature and pressure that were adjusted to the final NVT production run. For example, a system that was intended to be simulated at 400 K, most of the runs were performed at 400 K (the same holds for 500, 600 and 700 K). The 1000 K NVT simulation, however, constitutes an exception and was run for about 100 ps. The equilibration at the final temperature was performed for 20.2 ns. The last snapshot was extracted to start the final NVT production run, as described below.

Finally, a NVT production run was performed. 400 ns of trajectory were obtained for the simulations at 400 and 500 K and 580 ns for the simulation at 600 and 700 K.

2.3 Atom density and free energies

Atom density profiles were generated, describing the local distribution of atoms along the z-axis of the simulation box. As stated above, the z-axis was defined as being perpendicular to the interface and z = 0 is defined as being at the center of the LLZO. For the calculation samples were taken every 25 ps starting at 200 ns of each simulation of each simulation. The simulation box was divided into 520 bins of equal width along the z-axis. The mean number of atoms per species i in each bin was calculated and divided by the volume of the bin V. This is described in eqn (1) with T_tot being the total number of time steps that were part of the calculation.

The free energy profiles of Li⁺ and TFSI⁻ ions at the LLZO–SPE interfaces were calculated from the mean atom density profiles ρ using the relation

F = −k_B·T·ln(ρ) (2)

Due to the low density of Li⁺ in the interfacial region, sampling was increased to every 2 ps and the size of the bins was set to equal 0.2 Å for all simulations.

The free energy was normalized to equal zero inside the LLZO crystal. This was achieved by subtracting the value of the free energy in the bin at z = 15.3 Å.

For each region described in Section 3.1, an average value of atoms/molecules found in the layer was calculated,

where n(j,l,t) is a population of a specie j in a corresponding layer l at time step t, and t_max is the total number of time steps included in the calculation.

2.4 Displacement variance

The variance of displacement (MSD_Var) was calculated according to

Var_k(Δr_i) = 〈Δr_i²〉_k − (〈Δr_i〉_k)² (4)

where i is an ion type (Li⁺ or TFSI⁻), and Δr_i is a particle displacement in x, y or z direction. The variables in brackets, 〈..〉, indicate an average over multiple starting times (t₀) and all i particles in a k-th layer at t₀. The layers along the z-axis were selected following the interfacial structure of the ions, as indicated in Section 3.1. The specific positions of the k-th layer boundaries, used to obtain the z-resolved Var_k(Δr_i), are indicated in the Table S1.

Following ref. 9, we have not calculated the common mean-square displacement (MSD) function. Due to the presence of the ceramic surfaces in the xy-plane, the particles near the slab encounter a restriction in movement along the z-axis, which results in a net drift towards the polymer bulk. This, in turn, artificially increases their MSD. This shift is taken into account by the squared linear term in eqn (4).

2.5 Equilibration diagnosis

In the context of MD simulations, a physical property can be considered equilibrated when its average value does no longer change with time.³⁰ This applies, e.g., to the potential energy or density of a system measured as a function of time. The focus of this paper is the equilibration of the ion distribution in the non-homogeneous CPE system.

At equilibrium between two states A and B (e.g., reflecting different regions of the system) the number of transitions (e.g. hopping processes of ions) from A und B (S_AB) in a given time interval is basically the same as the number of reverse transitions (S_BA), reflecting detailed balance principle. When starting a simulation, the initial condition is typically a non-equilibrium condition. Before equilibration is reached, one expects S_AB ≠ S_BA. Furthermore, the diagnosis of equilibration requires that some transitions have occurred at all. Here we suggest an empirical equilibration coefficient 0 ≤ EC ≤ 1 such that larger values of EC correspond to a higher degree of equilibration. Two conditions for equilibration are implemented. First, the number of transitions in both directions needs to be similar, i.e. S_AB ≈ S_BA since otherwise the system is still in the process of equilibration. Second, the total number of transitions must not be too small. Specifically, we choose

with λ = 10. EC can only approach the upper limit if the total number of transitions is at least of the order of λ. The choice of this parameter is rationalized in the SI although there is no strict argument why its value should be exactly 10. Anyway, none of the conclusions in this work are changed if we had chosen, e.g., λ = 5. Furthermore, the choice of the min-function guarantees that significant differences between S_AB and S_BA reduce the value of EC.

For a complex energy landscape with significant substructures, as shown in Fig. 2, equilibration may first occur only on a local scale so that the EC value can differ for the different transitions. At the long-time limit, the total system should naturally have reached equilibrium. As will become clear in the following analysis, the energy landscape from Fig. 2 reflects some basic properties of the current system.

Fig. 2 Sketch of a model potential energy landscape where a high barrier separates the high-energy region on the left side (phase 1) and the low-energy region on the right side (phase 2). It is supposed to reflect expected key features of CPEs, namely relative fast inter-phase and much slower intra-phase equilibration. Left: Initial population (blue bars). Middle: Population at a time interval at intermediate times, including a sketch of the number of transitions during that time interval. Only a few inter-phase transitions to the lower-energy right region have occurred. Right: Population at the long-time limit. On these longer time scales several forward–backward inter-phase transitions are seen. The quantity EC (equilibrium coefficient) is defined between 0 and 1 (see text for definition) and expresses whether the transitions across a free energy barrier are in equilibrium (EC approaching 1) or not (EC close to 0).

2.6 Charge integral

For a charge distribution ρ_q(z) one can calculate the accumulated charge distribution

The notation is guided by the observation that according to the Poisson equation E_eff(z) is proportional to the local electric field E(z) in z-axis if z_b is located in a neutral bulk region, where the electric potential is zero. The change in sign of E_eff(z) can be attributed to charge reversal.

To calculate E_eff(z) first the charge distribution ρ_q(z) was determined. The simulation box was divided into bins with a width of 0.326 Å along the z-axis. Taking samples every 50 ps starting from 200 ns of the NVT trajectory the charges of each individual atom in the bin were summed up and the average over time was calculated.

Then the accumulated charge distribution was determined following eqn (6)z_b was chosen to be the point farthest from both interfaces, which is equivalent to the end of the simulation box on the LEFT and RIGHT side and leads to a symmetric set up.

3 Results

3.1 Ion distribution

Structural regimes. To help understanding the general system characteristics, the atomic distribution within the simulation box is first analyzed. This serves as a basis for dividing the simulation box into specific regions that are of importance for further analysis. A couple of discrete regions are defined based on the average number density of Li⁺ and TFSI⁻ that reside in these regions. This leads to regions that are partly overlapping, as can be seen in the colored bars at the top and bottom of Fig. 3. The exact boundaries are listed in Table 1. Structural regimes are defined based on the system run at 700 K, since the ion distribution has reached equilibrium at this temperature, as will be discussed in Section 3.2. Data for the systems between 400 K and 600 K can be found in the Fig. S2 and S3.

Fig. 3 Mean atom density profile of Li⁺, TFSI⁻ and O_PEO at 700 K including definition of regions in the top and bottom bars.

Table 1 Discretization of the simulation box into layers perpendicular to the surface defined by ranges z ∈ {..} Å

Atom	PEO bulk LEFT	1st layer LEFT	LLZO region	1st layer RIGHT	PEO bulk RIGHT
Li^+	−55, −32	−32, −23	−23, 23	23, 32	32, 55
TFSI^−	−55, −26	−26, −19	—	19, 26	26, 55

---

If focusing on the TFSI⁻ distribution, there are three distinct regions: the bulk LLZO (free from TFSI⁻), the bulk SPE, and a TFSI⁻-enriched layer at the very surface region on LLZO. Here, the TFSI⁻ density shows a maximum between 19 < |z| < 26, which is named TFSI^1stlayer. At the position of TFSI^1stlayer, the O_PEO density is very small. It increases slightly, close to the LLZO surface, indicating some small direct interaction between LLZO and PEO.

If focusing on Li⁺, on the other hand, three somewhat different regions appear, depending on that a Li-free region appears, next to the interface at 23 < |z| < 25, in the TFSI⁻-enriched layer. This defines an LLZO region, which also possesses a stoichiometric excess of Li as compared to the starting configuration. The periodic structure of Li⁺ ions in the LLZO crystal is clearly visible between −23 < z < 23. In the following, this region is called Li^LLZO,excess.

Outside of the Li⁺-depleted region is then a region of the SPE where Li⁺ displays a higher concentration. A Li⁺ peak is formed around |z| = 30, i.e. in the SPE phase. This region is defined as Li^1stlayer (23 < |z| < 32). A noticeable maximum of O_PEO atom density is also observed in that region, likely due to the coordination strength of Li⁺ to the polymer ether oxygens. Finally, there is a bulk SPE region where Li⁺ has an atomic density corresponding to a stable Li⁺ concentration throughout.

Average population in the structural regions. In order to analyze specific numbers of ions residing in the structural regions defined above, it is useful to study the average ion population. This is calculated according to Eqn 3, and thereby shows the average ion populations N(j) with j ∈ Li^excess, Li^1stlayer and TFSI^1stlayer. For Li^LLZO,excess we concentrate on the excess lithium ions, i.e. subtracting the number of lithium ions, which were originally in the LLZO slab, from the total number N(Li^LLZO,excess). This number is denoted by N(Li^excess). Fig. 4 summarizes the data, where each data point represents an average number of the respective ions in the previously defined regions, Table 1, at its specific simulation temperature.

Fig. 4 Li^excess, Li^1stlayer and TFSI^1stlayer average population analysis. The analysis was performed on the trajectory starting from 200 ns of the NVT run. Snapshots were extracted every 1 ns.

At all temperatures studied, Li⁺ are observed to accumulate in the Li^LLZO,excess region. In other words, the number of Li⁺ exceeds the stoichiometric number initially present in the crystal structure of LLZO. This implies that the LLZO region becomes positively charged, giving rise to a finite number of excess Li⁺. In turn, the PEO phase becomes depleted from Li⁺. Similarly, TFSI⁻ accumulate at the LLZO surfaces in TFSI^1stlayer.

A closer look shows that for both N(Li^excess) and N(TFSI^1stlayer) there is a drastic reduction at 400 K and a smaller reduction at 500 K as compared to the expectation from extrapolation of the high-temperature behavior. It will be argued below that this is due to the non-equilibrium effects.

Charge compensation at the interface. Interestingly, the number of TFSI⁻ at the surfaces exceeds the number of Li⁺ in Li^LLZO,excess. For a more general analysis of the charge distribution, we show in Fig. 5 the integral of the charge density as a function of the upper boundary at 700 K. The lower boundary has been chosen to lie in the bulk region (see Section 2.6). The z-axis, in Fig. 5, for the RIGHT interface has been shifted so that the positions of the respective peaks, occurring for both sides, agree as well as possible. The charge distribution close to both interfaces is very similar, except for the region at the Li-surface layer. This indicates a general mechanism of self-assembly in the polymer phase, independent of the specific details and structure of the LLZO surface.

Fig. 5 Integrated total charge around the LEFT and RIGHT LLZO–SPE interface for the simulation at 700 K.

We observe that E_eff(z) shows oscillations which even vary between positive and negative values. This is a clear signature of overcharging.³¹ As discussed in³¹ overcharging reflects a complex interplay of steric and electrostatic interactions, typically involving the presence of charged ions/molecules with quite different sizes, albeit typically involving multivalent ions. Note that in the present study only monovalent ions are present. The overcharging seen in Fig. 5 is consistent with the observation that N(TFSI^1stlayer) is larger than N(Li^excess), i.e. the negative compensation of the accumulated positive charge in LLZO region is larger than required to achieve charge neutrality. Indeed, on the RIGHT LLZO–SPE interface there are six locations where the electric field changes its sign, highlighting a relatively complex structural arrangement. Finally, it appears that the high N(TFSI^1stlayer) can exist despite the strong electrostatic anion–anion repulsion, possibly because it is partially compensated by the positive charge in the LLZO region and the formation of Li^1stlayer.

3.2 Equilibration analysis

Fig. 6(a) and (b) show the evolution of the N(Li^excess) and N(TFSI^1stlayer) as a function of time and temperature, allowing the tracking of the ionic transitions in and out of these regions. This is key for being able to assess the equilibrium properties of the system. The N(Li^excess) and N(TFSI^1stlayer) values are much larger than zero already from the very start, i.e. t = 0 ns. That is due to the equilibration procedure that preceded the production run, as described in Section 2.2.

Fig. 6 Li^LLZO,excess and TFSI^1stlayer time-dependent population analysis. The analysis was performed on the full NVT run. Snapshots were extracted every 1 ns.

As mentioned in Section 2.5, there are two criteria for evaluating whether the ion distribution has reached equilibrium: there should be a similar number for forward–backward jumps in or out of a specific region, and also a sufficient number of such transitions must occur. Starting with the Li^LLZO,excess region, no forward–backward jumps occur for 400 K and 500 K during 400 ns (see Fig. 6(a)), implying that the equilibrium condition is not reached. Instead, the equilibrium criterion seems to be fulfilled at the highest temperature, and to some degree at 600 K. In contrast, in TFSI^1stlayer equilibrium appears to be reached for all temperatures since in Fig. 6(b), the population of TFSI^1stlayer shows frequent transitions in and out of that region for the timescale of 1 ns or less at all temperatures. Thus, it is only at higher temperatures (in particular 700 K) that both TFSI^1stlayer and Li^LLZO,excess appear to fulfill the equilibrium criteria.

These observations can be quantified by the equilibration coefficient EC which can vary between 0 and 1 and where values close to 1 indicate equilibration. The results for EC are shown in Fig. 7. A few important observations can be made: (i) TFSI^1stlayer is equilibrated for all temperatures because EC(TFSI^1stlayer) ≈ 1; (ii) Li^LLZO,excess is well equilibrated at T = 700 K since EC(Li^LLZO,excess) ≈ 0.6; (iii) at T = 600 K, Li^LLZO,excess is not fully equilibrated considering EC(Li^LLZO,excess) ≈ 0.3; and (iv) Li^LLZO,excess at 500 K or below is not equilibrated at all because EC(Li^LLZO,excess) = 0. Thereby, it can be stated that at all temperatures TFSI^1stlayer reached an apparent equilibrium distribution. The apparent equilibrium implies that it is adjusted to the number of Li⁺ in Li^LLZO,excess. Moreover, there is a notable asymmetry between the behavior of Li⁺ and TFSI⁻ ions: the energy barrier for changes in the N(Li^LLZO,excess) is significantly higher than that for N(TFSI^1stlayer). As a consequence the EC(Li^LLZO,excess) value strongly decreases upon decreasing temperature as outlined above. In Section 3.3 this dependence is related to the properties of the free energy landscape.

Fig. 7 Equilibration coefficient (top) calculated for Li^LLZO,excess and TFSI^1stlayer. The analysis was performed on the trajectory starting from 200 ns of the NVT run, snapshots were extracted every 1 ns. S_t (bottom) is a sum of forward–backward transitions counted for a specific region. S_t at 600 and 700 K are used in the following transition rate calculations, hence they are indicated explicitly with the arrows.

A visual inspection of the time-dependence of N(Li^excess) and N(TFSI^1stlayer) at 600 and 700 K in Fig. 6 indicates that both populations are positively correlated; i.e. TFSI⁻ accumulation at the LLZO surfaces follow Li⁺ accumulation in Li^LLZO,excess. This correlation is quantified in Fig. 8. Here, it is seen – to a good approximation – that an additional Li⁺ in Li^LLZO,excess gives rise to an additional TFSI⁻ in TFSI^1stlayer. Together with the observation of much less forward–backward transitions for Li^LLZO,excess as compared to TFSI^1stlayer, this shows that the accumulation of ions in Li^LLZO,excess is the driving force for the increasing number of ions in TFSI^1stlayer. This aligns well with the observation of significant charge overcompensation, which also indicates a strong driving force of TFSI⁻ accumulation as a consequence of Li⁺ accumulation in Li^LLZO,excess.

Fig. 8 Correlation of TFSI^1stlayer (sum of both LLZO sides) and Li^LLZO,excess at 600 and 700 K. The analysis was performed on the trajectory starting from 200 ns of the NVT run. Snapshots were extracted every 1 ns.

The above analysis supports the following scenario which is compatible with the sketch in Fig. 2: the accumulation of Li⁺ in Li^LLZO,excess is the slowest process in the LLZO/PEO:LiTFSI system, governing the complete equilibration. This is only reached at 700 K. Already on shorter time scales, N(TFSI^1stlayer) manages to acquire its equilibrium value within the polymer phase which, however, is conditioned on the instantaneous number of Li⁺ in Li^LLZO,excess. This is the consequence of the strong coupling, seen in Fig. 8.

These results may also explain the strong reduction of N(Li^excess) and N(TFSI^1stlayer) in Fig. 4 at 400 K and (to some degree) at 500 K. Due to the slow equilibration process, N(Li^excess) did not have the time to equilibrate during the simulation time scale of 400 ns. Although intra-phase equilibration of TFSI⁻ is fast, its conditioning on N(Li^excess) also explains the strong reduction of N(TFSI^1stlayer) at the lowest temperatures.

3.3 Free energy landscape analysis

The equilibration properties of the Li and TFSI distributions can be rationalized via a free energy analysis. Transitions from PEO to LLZO occur on a very short time scale but are rare events. This implies that the interfacial region becomes poorly sampled. Therefore, the sampling interval needs to be carefully adjusted. Since basically no transitions have been observed at 400 K and 500 K, see Fig. 6, this analysis is restricted to 600 K and 700 K. The free energy for T = 600 K and T = 700 K is shown in Fig. 9(a)).

Fig. 9 Li⁺ free energy barriers at the LEFT and RIGHT LLZO–SPE interface (a) at 700 K and 600 K based on the NVT trajectory starting at 200 ns; (b) same data as in (a) but the barrier in the PEO phase for 600 K is shifted (shift values: LEFT (−0.18 Å, 360 K), RIGHT (+0.18 Å, +360 K). The average slope of the broken lines is 1850 K Å⁻¹.

To characterize the barrier height ΔF for the transition from the polymer phase to the Li^LLZO,excess, the free energy difference between Li^1stlayer, corresponding to the minimum of the free energy data around z ≈ ±28 Å and the free energy barrier around z ≈ ±23 Å is determined. For the simulations at 700 K, it turns out that ΔF slightly varies when comparing the LEFT and the RIGHT side of the LLZO slab (ΔF(LEFT) = 8000 K vs. ΔF(RIGHT) = 7500 K). In what follows we use the Boltzmann averaged value of ΔF = 7700 K. Furthermore, for the prediction of the final rate, information about the shape of the minimum of the 1st Li peak is required. Starting from a minimum at z₀, it is fitted by a parabola: F(x) = k·(z − z₀)². Thereby, we obtain k(LEFT) = 4.9 × 10⁻²¹ J Å⁻² with k(RIGHT) = 5.2 × 10⁻²¹ J Å⁻² and the average k = 5.1 × 10⁻²¹ J Å⁻². By employing transition state theory, the hopping rate from the 1st Li peak to the LLZO surface can be predicted via the relation

Γ = κν exp(−βΔF) (7)

where ΔF is the free energy barrier, ν is the vibrational frequency and κ ≤ 1 is the transmission coefficient which serves as a correction factor. With the relation k = (1/2)mω² = 2mπ²ν² the frequency can be expressed as

Using the lithium mass of 1.15 × 10⁻²⁶ kg results in ν = 1.5 × 10¹² s⁻¹. Without the correction factor κ this yields Γ = 2.5 × 10⁷ s⁻¹ at T = 700 K. With this rate, it becomes possible to estimate the expected number of transitions. The Li^1stlayer (sum of both LLZO interfaces) contains on average 15 ions, see Fig. 4. Thus, the expectation is to observe 2.5 × 10⁷·15·380 × 10⁻⁹ ≈ 140 jumps per 380 ns. However, in the actual MD simulations, only 7 transitions onto the LLZO surface are observed during the 380 ns of simulation time. Note that the 14 transitions, mentioned in Fig. 7, also take into account the transitions back to the polymer region. The results could be reproduced with a transmission factor of κ = 0.05. Formally, this could be rationalized by choosing an effective lithium mass that is approximately a factor of (1/0.05)² ≈ 400 higher than the elementary mass. Such an increase could be explained by a strong binding of the lithium ions to the polymers. Indeed, as shown in Fig. 3, there is a direct interaction between lithium ions and PEO in the first lithium layer.

Remarkably, the number of excess ions hardly changes, if at all, when comparing 600 K and 700 K, see Fig. 9(a). This implies that there is a well-defined limit of excess lithium ions, beyond which the positive charging starts to become prohibitive. It is likely that the energy gain of lithium ions to enter the LLZO due to the van der Waals and Coulomb interaction is counteracted by the overall positive charge. The work to identify possible mechanisms is in progress. Anyway, in the equilibrium situation one can thus expect that also at the barrier, the presence of several lithium ions is prohibitive as well.

When analyzing the free energy landscape at 600 K, slight differences are observed for the polymer region relative to the LLZO region as compared to the data at 700 K, see also Fig. 9(a). A careful analysis shows that the free energy landscape can be superimposed in the PEO region after two transformations: (i) shift of free energy by 360 K and (ii) shift of the z-axis by Δz = 0.18 Å. To understand the shift of free energy, one needs to take into account that the number of lithium ions in the first peak is basically temperature-independent. This is consistent with the observation that the number of excess ions is temperature independent as well. The reasons were already discussed above. In contrast, for a constant free energy a significant decrease of the number of lithium ions would have been expected. This directly explains the shift in the free energy. Anyway, this shift does not influence the barrier for the transitions from the polymer to the LLZO. The second transformation reflects that the simulation box in the NpT ensemble is slightly shorter for lower temperatures, see Fig. S4. The nearly perfect agreement of the free energy data after these transitions, as highlighted in Fig. 9(b), clearly shows that there is basically a temperature independent free energy landscape in the polymer regime. However, it is expected that the reduced distance between the LLZO regime and the first lithium peak reduces the barrier height. Indeed, close to the free energy barrier the free energy curve in the polymer regime can be approximated by a straight line with slope 1850 K Å⁻¹. Thus, a value of Δz = 0.18 Å reduces the barrier by approximately 0.18·1850 K = 330 K.

In general, the free energy difference ΔF = 7700 K between the barrier and the minimum has an entropic and an enthalpic contribution, i.e. −k_BTΔS + ΔE. Qualitatively, ΔS expresses possible differences in the phase space of possible lithium sites when comparing the barrier region with the first lithium peak. Since sites in the crystal are quite localized as compared to the z-region of the first lithium peak and since the barrier largely reflects also properties of the crystal, one may expect that ΔS < 0, reflecting a reduction of available phase space for lithium ions close to the barrier. An estimation of this decomposition of ΔF can be performed by comparing the number of lithium hopping processes from the polymer region to the LLZO region for both temperatures 600 K and 700 K. Since we are mainly interested in the dynamics close to the equilibrium value of N(Li^1stlayer), we consider the number of transitions to the LLZO region if N(Li^excess) ≥ 42. For this purpose, the transitions at 600 K are only recorded in the final 300 ns. There we observe 2 transitions from the polymer to the LLZO. Together with the 7 transitions during 380 ns at 700 K we estimate Γ_est(700 K)/Γ_est(600 K) ≈ (7/2)·(300/380) ≈ 3. Obviously, due to the small number of transitions this value has quite some uncertainty. Therefore, we define γ = Γ_est(700 K)/Γ_est(600 K) and perform the subsequent analysis for general γ and finally check the outcome for γ = 3.

Taking into account the additional barrier reduction, discussed above, one has

and

−700 K·ΔS + ΔE = 7700 K (10)

The first equation yields ΔE = (4200 log(γ) + 2310) K and the second equation ΔS ≈ 6log(γ) − 7.7 which for γ = 3 gives rise to ΔS = −1.1. This would imply that there is indeed a reduction of the phase space, available for lithium ions around the barrier. However, one should take into account the uncertainty of the value of γ which is mainly determined by the small number of transitions at T = 600 K. Considering also the options of just a single or three transitions, we approximate the uncertainty of γ as 2 ≤ γ ≤6.

As stated previously, it can be seen that TFSI⁻ crosses from the bulk regime to the TFSI^1stlayer much more frequently than Li⁺. To rationalize this observation, the z-dependent free energy for TFSI⁻ is also determined; see Fig. 10. It turns out that the free energy barrier height to reach TFSI^1stlayer from the bulk region of the polymer electrolyte is approx. a factor of 3 smaller than for the Li⁺ transitions (in the range between 2500 K and 3000 K). Thus, in agreement with the numerical observations, the low free energy barrier enables fast TFSI⁻ transitions. In particular, this allows N(TFSI^1stlayer) to reach the apparent equilibrium and to follow N(Li^LLZO,excess), as discussed above.

Fig. 10 TFSI⁻ free energy barriers for 600 K and 700 K based on the NVT trajectory starting at 200 ns.

3.4 Temperature dependence

LLZO–SPE composite electrolytes are generally considered to be operated at elevated temperatures as compared to Li-ion batteries based on liquid electrolytes, but temperatures as high as 400 K are likely too high for most applications. Therefore, MD simulations targeting realistic temperatures need to be performed at lower temperatures than those performed here, while higher temperatures might be necessary to achieve equilibration and sufficient statistical sampling. In order to determine whether the high temperature simulations can be used to provide a reliable insight about the low temperature processes, the temperature dependence of relevant properties needs to be analyzed.

Charge distribution. Fig. 11 shows the integrated charge around the LEFT and RIGHT interface of the LLZO slab for simulations at 400 K, 500 K, 600 K and 700 K. The charge distribution is similar at all temperatures, despite the reduced N(Li^excess) at low temperatures. This is promising since structural details, characteristic at low temperatures, can be found already at 700 K. To achieve a good overlap of the positions of the individual peaks, the z-axis has been shifted as compared to the simulations at 700 K (600K: 0.2 Å, 500 K: 0.35 Å, 400 K: 0.45 Å). This reflects the different box lengths in the NpT ensemble, as already discussed in the context of Fig. 9. The impact of the shifting procedure is shown in the SI.

Fig. 11 Integrated total charge around LEFT and RIGHT interface and for simulations at 400 K, 500 K, 600 K and 700 K with z-axes shifted individually for each temperature (see text for explanation).

Dynamics parallel to the LLZO surface. Fig. 12 and 13 shows the spatially resolved Li⁺- and TFSI⁻-displacement variance for individual bins along the z-axis. The MSD_Var data are analyzed at 1 ns. Thereby, the Li⁺ and TFSI⁻ mobility characterized by displacement parallel to the LLZO surface, along the x-axis, is seen. Mobilities of the ions, perpendicular to the LLZO surfaces, are shown in Fig. S6 and S7. In general, Fig. 12 shows that at all temperatures studied, the lowest mobility of Li⁺ is observed directly at the PEO and LLZO interface. This holds specifically for the regions z ∈ {±15, ±23} that cover the surface of LLZO and z ∈ {±23, ±32} that correspond to Li^1stlayer. Moreover, the Li⁺ mobility at 700 K is similar in the bulk LLZO and the bulk PEO phase. This is in contrast to the simulations at 400 K. Naturally, the MSD_Var is significantly smaller at lower temperatures. However, Li⁺ moves 3 times faster in the ceramic phase than in the polymer phase, indicating a much higher activation energy for the lithium diffusivity in the polymer phase. Due to the fast intra-phase equilibration and the mild sensitivity of diffusivity on salt concentration, these observations are expected to remain if equilibrium simulations were possible also at 400 K.

Fig. 12 Displacement variance of Li⁺ at 1 ns in x-axis (parallel to LLZO surface) as a function of ion position perpendicular to the LLZO surface. The vertical gray dashed lines indicate the bins used for the calculation of specially resolved Var(Δr). Last 300 ns of NVT trajectories at 400 and 700 K were taken for the analysis.

Fig. 13 Displacement variance of TFSI⁻ at 1 ns in x-axis (parallel to LLZO surface) as a function of ion position perpendicular to the LLZO surface. The vertical gray dashed lines indicate the bins used for the calculation of specially resolved Var(Δr). 300–400 ns of NVT trajectories at 400 and 700 K were taken for the analysis.

In Fig. 13, spatially resolved MSD_Var at 1 ns are plotted for TFSI⁻. Qualitatively, a similar behavior is observed as compared to Li⁺ − MSD_Var in Fig. 12, namely a strong reduction of mobility of the anion is seen at the ceramic–polymer interface, i.e. in the TFSI^1stlayer region, as compared to the PEO bulk region. This reduction, by one order of magnitude, is already visible at 700 K whereas for Li⁺ such a strong reduction is only observed at 400 K.

Rate prediction from free energy landscape. The discussion of the free energy landscape has been formulated such that it is also applicable to lower temperatures. To estimate the hopping rate at 400 K we use the general relation

Γ(400 K) = κν exp(ΔS)exp(−β(ΔE − Δz·1850 K Å⁻¹)). (11)

Here we use ν = 1.5 × 10¹² s⁻¹ and κ = 0.05. Furthermore, from the data in Fig. S5 we take Δz = 0.45 Å. Finally, we take the general expressions of ΔE and ΔS as derived above. Straightforward calculation yields

Γ(400 K) = 0.8 × 10⁶γ^−4.5 s⁻¹ (12)

This yields Γ(400 K, γ = 2) ≈ 4 × 10⁴ s⁻¹, Γ(400 K, γ = 3) ≈ 6 × 10³ s⁻¹, Γ(400 K, γ = 6) ≈ 3 × 10² s⁻¹. This implies that for the case γ = 3 the probability for a single transition is approx. 3% during the final 380 ns of the simulation. This rationalizes that in the trajectories in Fig. 6 no transition has been observed at 400 K. Naturally, this range of possible rates between 4 × 10⁴ s⁻¹ and 3 × 10² s⁻¹ should only be considered as a first reasonable estimate, based on the assumption that our insight about the free energy landscape can be quantitatively extrapolated to lower temperatures.

4 Discussion

Taking into account all arguments presented in the results section, we argue that ion transport in the non-homogeneous systems should not be studied for systems that are equilibrated at the usual 400 K (and even 500 K). This is due to the low transition rate associated with the high energy barrier (see Section 3.3) which prevents ions from reaching a sufficient number of phase exchanges to reach an equilibrium distribution, and may therefore lead to inadequate interpretations of the results.

Several MD studies of the LLZO/PEO interface have been performed at temperatures lower than 700 K and simulation times not exceeding 400 ns.^14–17 This allows for a direct assessment of the interpretation of these results in light of the importance of equilibrated ion distribution between the phases. For example, umbrella sampling has been used before¹⁵ to calculate the free energy profiles of Li⁺ when leaving the LLZO surface (Fig. 3(a) and (b)). However, the density profiles, as shown in Fig. 2(a),¹⁵ indicate that only about 3 Li⁺ entered Li^LLZO,excess from the polymer phase and that Li^1stlayer is not formed. In Fig. 4, it is shown that about 28 and 38 Li⁺ entered LLZO from the PEO phase, in the simulations at 400 K and 500 K, respectively. Based on the behavior at 600 K and 700 K, these values should be only considered as a lower boundary. Thus, the free energy profiles in¹⁵ are likely not providing the actual equilibrium behavior.

Moreover, one may still question whether the observed N(Li^excess) values at 400 K or 500 K (Fig. 4) were artificially influenced by the elevated temperatures used during the initial equilibration procedure, as discussed in Section 2.2. We argue that at 400 K and 500 K, the equilibrium value of N(Li^excess) should be even higher. The estimated transition rate at 400 K suggests that about 1 ms simulation is needed at 400 K to observe a jump of another Li⁺ to Li^LLZO,excess. Hence, it is not possible to verify our claim directly, but we can provide two indirect arguments that support the statement: (1) a single transition of Li⁺ to Li^LLZO,excess is observed in the NVT run at 500 K (see Fig. 6), which implies that saturation of LLZO with additional Li⁺ has not been achieved and it is more likely for Li⁺ to enter the LLZO crystal than leave it. The direction of phase exchange when Li⁺ enters the LLZO would be very unlikely if the equilibrium value of N(Li^excess) were much smaller than observed, in extreme case close to zero. (2) At 600 K and 700 K the N(Li^excess) values are basically identical. Hence we envision that the values will not be drastically different also for the lower temperatures.

The z-dependence of the MSD suggests that there is no increase of the MSD close to the interface as compared to the bulk regime. Due to the finite size of our simulation box, we cannot strictly exclude that the MSD would decrease again further in the bulk. Since, however, in our simulations the oscillations of E_eff(z) have basically decayed far away from the interface, we believe that we already see the bulk MSD. As a consequence, the regions close to the interfaces cannot be responsible for an increased overall conductivity upon adding LLZO to the polymer electrolyte. This somehow contrasts to the conclusions drawn from.¹⁶

The impact of the polymer chain length on ion dynamics in a bulk polymer electrolyte has been studied in previous work.³² It has been shown that, concerning the binding properties of lithium, polymers chains as long as 10 monomers behave similarly to much longer polymer chains. This implies that binding properties of 20-monomer-long polymer chains represents the Li⁺ binding sufficiently well. Moreover, regarding simulations of Li⁺ in a composite system Kozdra et al.¹⁷ have compared a 900-monomer-long polymer chain with 45 chains 20-monomer-long. In the system with the 900-monomer-long polymer chain the anion distribution is strikingly asymmetric as compared to the 20-monomer-long polymer chains which suggests that in the systems with longer polymer chains not only the Li⁺ distribution among phases is not equilibrated but also the equilibrium anion distribution may be challenging to achieve at 400 K. This aligns with the fact that longer chains require longer equilibration times. When it comes to the ceramic–polymer interface simulations, in¹⁷ several simulation setups are tested and the most pronounced difference of residence time and the density profiles are observed at the polymer–ceramic junction. However, the trends are the same, excess Li⁺ arrive and remain the longest at the LLZO surfaces. This indicates some generality of this observation, although very likely not all details of the interface are accurately captured.

We could quantify the underlying reason of this slow equilibration at low temperatures by the presence of a high free-energy barrier between LLZO and the polymer electrolyte. Since the number of excess ions is identical for 600 K and 700 K one may hypothesize that also the equilibrium structure is largely temperature independent. Furthermore, at 400 K the Li⁺ diffusivity is significantly higher than in the polymer phase which is an important aspect when using LLZO as active material in hybrid electrolytes.

In the present case the equilibrium ion distribution could finally be reached for 600 K but not below. Note that in principle this temperature may somewhat depend on system properties such as salt concentration and, e.g., choice of the force field. However, the impact of the free energy barrier for Li⁺ transport between both phases and the subsequent implications for an appropriate simulation setup, relying on long simulations at very high temperatures, is likely to be generalizable.

As a consequence, the important message of this study is that simulations at very high temperatures (e.g. 700 K for the present case) should be used as a strategy to gain physical insight about the non-homogeneous solid-state electrolyte systems at lower temperatures. The two main arguments are: (1) simulations at high temperature are necessary to reach the equilibrium ion distribution and structural features at the interface in such systems, (2) the charge distribution in the polymer phase at 700 K resembles the charge distribution at 400 K, indicating sufficient level of similarity when comparing this property. From a practical perspective one might use the equilibrated structure at 700 K and continue the simulation at the low temperature of interest.

5 Conclusion

Ion distribution in a non-homogeneous system was studied using classical molecular dynamics simulations at elevated temperatures. An interface built between a LLZO ceramic Li⁺ ion conductor and a LiTFSI-doped PEO polymer electrolyte was used as a model system to investigate the impact of ion distribution in the simulation box on selected static and dynamic properties.

Population analysis of Li⁺ and TFSI⁻ in the selected regions of the simulation box show that the equilibrium ion distributions are not reached in simulations at T = 400 K and 500 K. This statement is confirmed with the analysis of the free energy landscape. Hence, high temperature simulations (or other strategies) are necessary to determine the equilibrium distribution of ions in the non-homogeneous systems such as LLZO/PEO:LiTFSI. A suitable strategy to reach the adequate interfacial structure could be to equilibrate the system at 700 K and then gradually quench it to 400 K. Effectively, this means imposing the 700 K structure to 400 K simulation. This is considered to be a better strategy than starting a 400 K simulation with randomly distributed LiTFSI molecules in the polymer phase.

On the basis of the statements above, care needs to be taken to equilibrate the ion distribution in such interfacial systems. The reason for this is that the ion distribution has a significant impact on the structure and ion mobility at the interface. Here, the structure of the interface is described by the positions of the peaks in the density profiles and the number of species in the distinct regions in the simulation box. In addition, equilibrium ion distribution at the interface is crucial for the accurate determination of the energy barrier, and thus for Li⁺ entering the ceramic phase from the polymer phase. Therefore, it is expected that the ion dynamics should only be studied once the ion distribution has reached equilibrium in the system. Otherwise, computational artifacts may arise in the results due to the non-equilibrium effects.

Author contributions

Melania Kozdra: conceptualization, software, validation, formal analysis, investigation, writing – original draft, writing – review & editing, visualization, supervision; Laura Hölzer: software, validation, formal analysis, investigation, writing – original draft, writing – review & editing, visualization; Daniel Brandell: resources, supervision, writing – review & editing, funding acquisition; Andreas Heuer: conceptualization, validation, formal analysis, resources, writing – original draft, writing – review & editing, supervision, funding acquisition.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article, namely examples of necessary input scripts for the conducted LAMMPS simulations, have been uploaded to github at the following link: https://github.com/laurahoelzer/equilibration_of_ion_distribution_at_polymer_ceramic_interface.

Supplementary information (SI) is available. See DOI: https://doi.org/10.1039/d5cp01988e.

Acknowledgements

MK and DB gratefully acknowledge generous support by the Swedish Energy Agency (project 50098-1) and STandUP for Energy, as well as computational infrastructure provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS), and AH and LH by the German Federal Ministry of Education and Research (BMBF) within the “FB2-TheoDat” (grant: 03XP0435A/E) as well as computer time on the Palma cluster at the University Muenster. Finally, we would like to acknowledge helpful discussions with Diddo Diddens and Gourav Shukla.

Notes and references

1. S. C. Sand, J. L. M. Rupp and B. Yildiz, A critical review on Li-ion transport, chemistry and structure of ceramic–polymer composite electrolytes for solid state batteries, Chem. Soc. Rev., 2024, 54, 178–200.
2. F. Croce, G. B. Appetecchi, L. Persi and B. Scrosati, (Al₂O₃, TiO₂) nanocomposite polymer electrolytes for lithium batteries, Nature, 1998, 394, 456–458.
3. A. G. Nguyen and C. J. Park, Insights into tailoring composite solid polymer electrolytes for solid-state lithium batteries, J. Membr. Sci., 2023, 675, 121552.
4. Z. Cheng, T. Liu, B. Zhao, F. Shen, H. Jin and X. Han, Recent advances in organic-inorganic composite solid electrolytes for all-solid-state lithium batteries, Energy Storage Mater., 2021, 34, 388–416.
5. I. Cuevas, K. Elbouazzaoui, M. Valvo, J. Mindemark, D. Brandell and K. Edström, Boron surface treatment of Li₇La₃Zr₂O₁₂ enabling solid composite electrolytes for Li-metal battery applications, ChemSusChem, 2024, 18(3), e202401304.
6. H. Zhang, F. Chen and J. Carrasco, Nanoscale modelling of polymer electrolytes for rechargeable batteries, Energy Storage Mater., 2021, 36, 77–90.
7. O. Borodin and G. D. Smith, Mechanism of ion transport in amorphous poly(ethylene oxide)/LiTFSI from molecular dynamics simulations, Macromolecules, 2006, 39(4), 1620–1629.
8. D. J. Brooks, B. V. Merinov, W. A. Goddard, B. Kozinsky and J. Mailoa, Atomistic description of ionic diffusion in PEO-LiTFSI: effect of temperature, molecular weight, and ionic concentration, Macromolecules, 2018, 51(21), 8987–8995.
9. A. Thum, D. Diddens and A. Heuer, Impact of charged surfaces on the structure and dynamics of polymer electrolytes: insights from atomistic simulations, J. Phys. Chem. C, 2021, 125(46), 25392–25403.
10. H. Gudla, Y. Shao, S. Phunnarungsi, D. Brandell and C. Zhang, Importance of the ion-pair lifetime in polymer electrolytes, J. Phys. Chem. Lett., 2021, 12(35), 8460–8464.
11. C. H. Chen and J. Du, Lithium ion diffusion mechanism in lithium lanthanum titanate solid-state electrolytes from atomistic simulations, J. Am. Ceram. Soc., 2014, 98(2), 534–542.
12. R. Jalem, M. J. D. Rushton, W. Manalastas, M. Nakayama, T. Kasuga, J. A. Kilner and R. W. Grimes, Effects of gallium doping in garnet-type Li₇La₃Zr₂O₁₂ solid electrolytes, Chem. Mater., 2015, 27(8), 2821–2831.
13. B. Akkinepally, I. N. Reddy, T. J. Ko, K. Yoo and J. Shim, Dopant effect on Li⁺ ion transport in NASICON-type solid electrolyte: insights from molecular dynamics simulations and experiments, Ceram. Int., 2022, 48(9), 12142–12151.
14. M. R. Bonilla, F. A. García Daza, P. Ranque, F. Aguesse, J. Carrasco and E. Akhmatskaya, Unveiling interfacial Li-ion dynamics in Li₇La₃Zr₂O₁₂/PEO(LiTFSI) composite polymer–ceramic solid electrolytes for all-solid-state lithium batteries, ACS Appl. Mater. Interfaces, 2021, 13(26), 30653–30667.
15. M. R. Bonilla, F. A. García Daza, H. A. Cortés and E. Akhmatskaya, On the interfacial lithium dynamics in Li₇La₃Zr₂O₁₂:poly(ethylene oxide)(LiTFSI) composite polymer–ceramic solid electrolytes under strong polymer phase confinement, J. Colloid Interface Sci., 2022, 623, 870–882.
16. H. A. Cortés, M. R. Bonilla, E. E. Marinero, J. Carrasco and E. Akhmatskaya, Anion trapping and ionic conductivity enhancement in PEO-based composite polymer–Li₇La₃Zr₂O₁₂ electrolytes: the role of the garnet Li molar content, Macromolecules, 2023, 56(11), 4256–4266.
17. M. Kozdra, D. Brandell, C. M. Araujo and A. Mace, The sensitive aspects of modelling polymer–ceramic composite solidstate electrolytes using molecular dynamics simulations, Phys. Chem. Chem. Phys., 2024, 26(7), 6216–6222.
18. Y. Chen, D. Liang, E. M. Y. Lee, S. Muy, M. Guillaume, M. D. Braida, A. A. Emery, N. Marzari and J. J. de Pablo, Ion transport at polymer–argyrodite interfaces, ACS Appl. Mater. Interfaces, 2024, 16(36), 48223–48234.
19. F. Scharf, A. Krude, P. Lennartz, M. Clausnitzer, G. Shukla, A. Buchheit, F. Kempe, D. Diddens, P. Glomb, M. M. Mitchell, T. Danner, A. Heuer, A. Latz, M. Winter and G. Brunklaus, Synergistic enhancement of mechanical and electrochemical properties in grafted polymer/oxide hybrid electrolytes, Small, 2024, 20, 2404537.
20. G. Bergerhoff, I. D. Brown and F. Allen, Crystallographic databases, IUCr, 1987, 360, 77–95.
21. W. Sun and G. Ceder, Efficient creation and convergence of surface slabs, Surf. Sci., 2013, 617, 53–59.
22. R. Tran, Z. Xu, B. Radhakrishnan, D. Winston, W. Sun, K. A. Persson and S. P. Ong, Surface energies of elemental crystals, Sci. Data, 2016, 3, 160080.
23. H. Xie, J. A. Alonso, Y. Li, M. T. Fernández-Díaz and J. B. Goodenough, Lithium distribution in aluminum-free cubic Li₇La₃Zr₂O₁₂, Chem. Mater., 2011, 23(16), 3587–3589.
24. P. J. in’t Veld and G. C. Rutledge, Temperature-dependent elasticity of a semicrystalline interphase composed of freely rotating chains, Macromolecules, 2003, 36(19), 7358–7365.
25. W. L. Jorgensen, D. S. Maxwell and J. Tirado-Rives, Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids, J. Am. Chem. Soc., 1996, 118(45), 11225–11236.
26. J. N. Canongia Lopes and A. A. H. Pádua, Molecular force field for ionic liquids composed of triflate or bistriflylimide anions, J. Phys. Chem. B, 2004, 108(43), 16893–16898.
27. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott and S. J. Plimpton, LAMMPS – a flexible simulation tool for particlebased materials modeling at the atomic, meso, and continuum scales, Comput. Phys. Commun., 2022, 271, 108171.
28. G. S. Larsen, P. Lin, K. E. Hart and C. M. Colina, Molecular simulations of PIM-1-like polymers of intrinsic microporosity, Macromolecules, 2011, 44(17), 6944–6951.
29. A. I. Jewett, D. Stelter, J. Lambert, S. M. Saladi, O. M. Roscioni, M. Ricci, L. Autin, M. Maritan, S. M. Bashusqeh, T. Keyes, R. T. Dame, J. E. Shea, G. J. Jensen and D. S. Goodsell, Moltemplate:a tool for coarse-grained modeling of complex biological matter and soft condensed matter physics, J. Mol. Biol., 2021, 433(11), 166841.
30. D. Frenkel and B. Smit, Understanding molecular simulation, Academic Press, San Diego, 2002.
31. W. M. de Vos and S. Lindhoud, Overcharging and charge inversion: finding the correct explanation(s), Adv. Colloid Interface Sci., 2019, 274, 102040.
32. J. Chattoraj, M. Knappe and A. Heuer, Dependence of ion dynamics on the polymer chain length in poly(ethylene oxide)-based polymer electrolytes, J. Phys. Chem. B, 2015, 119(22), 6786–6791.
33. W. Humphrey, A. Dalke and K. Schulten, VMD: Visual Molecular Dynamics, J. Mol. Graphics, 1996, 14, 33–38.

---