Analysis of Diffusion in Solid-State Electrolytes through MD Simulations, Improvement of the Li-Ion Conductivity in β-Li3PS4 as an Example.

Molecular dynamics simulations are a powerful tool to study diffusion processes in battery electrolyte and electrode materials. From molecular dynamics simulations, many properties relevant to diffusion can be obtained, including the diffusion path, amplitude of vibrations, jump rates, radial distribution functions, and collective diffusion processes. Here it is shown how the activation energies of different jumps and the attempt frequency can be obtained from a single molecular dynamics simulation. These detailed diffusion properties provide a thorough understanding of diffusion in solid electrolytes, and provide direction for the design of improved solid electrolyte materials. The presently developed analysis methodology is applied to DFT MD simulations of Li-ion diffusion in β-Li3PS4. The methodology presented is generally applicable to diffusion in crystalline materials and facilitates the analysis of molecular dynamics simulations. The code used for the analysis is freely available at: https://bitbucket.org/niekdeklerk/md-analysis-with-matlab. The results on β-Li3PS4 demonstrate that jumps between bc planes limit the conductivity of this important class of solid electrolyte materials. The simulations indicate that the rate-limiting jump process can be accelerated significantly by adding Li interstitials or Li vacancies, promoting three-dimensional diffusion, which results in increased macroscopic Li-ion diffusivity. Li vacancies can be introduced through Br doping, which is predicted to result in an order of magnitude larger Li-ion conductivity in β-Li3PS4. Furthermore, the present simulations rationalize the improved Li-ion diffusivity upon O doping through the change in Li distribution in the crystal. Thus, it is demonstrated how a thorough understanding of diffusion, based on thorough analysis of MD simulations, helps to gain insight and develop strategies to improve the ionic conductivity of solid electrolytes.

Introduction

To prevent further global warming by greenhouse gas emissions, it is necessary to move from fossil fuels toward renewable energy sources. For transport applications, other energy carriers such as hydrogen and batteries are considered. Of the current technologies that can replace fossil fuels in vehicles, batteries result in the lowest greenhouse gas emissions,[1] especially if renewable sources are used for the energy production.

However, safety concerns and the limited range of current battery electric vehicles are slowing down their implementation. Solid-state batteries are a promising technology[2,3] based on the much lower flammability risks, the higher energy density on the cell level, and lower self-discharge rate.

One of the prerequisites toward the realization of solid state batteries is the development of highly conductive solid electrolytes. In recent years several materials have been discovered which show conductivities comparable to liquid electrolytes. Room temperature ionic conductivities in the order of 10–3 S/cm have been reported in a range of lithium containing compounds such as[2] LLTO (Li3La2/3–TiO3)), argyrodites, LGPS (Li10GeP2S12), and LATP (Li1+AlTi2–(PO)4). Fewer sodium containing compounds with such high conductivities are known, most likely the result of less intensive research in this area, but several have been established, including β-alumina[4] and Na3PS4.[5] Recently, it has been shown that there is a relation between phonons, high ionic mobility, and low electrochemical stability in lithium-ion conductors,[6] which shows that the combination of high ionic conductivity and large electrochemical stability is challenging,[2,7] although electrochemical stability does not have very strict requirements. A solid electrolyte can be successful if its decomposition products are stable and have a reasonable ionic conductivity and low electronic conductivity,[7] similar to the functioning of solid-electrolyte interface (SEI) layers at electrodes in liquid electrolytes. The complex demands on solid electrolytes necessitates fundamental research toward solid electrolyte properties and new solid electrolyte materials.

Computer simulations are playing an important role in understanding and directing materials design toward improved battery performance. For example, calculations have shown that the electrochemical stability of solid electrolytes is enhanced by passivating decomposition products,[7] that strain can enhance diffusion,[8] how Li-ion diffusion can be increased in antiperovskites[9] and Na3PS4,[5,10] and why bond frustration is beneficial for Li-ion diffusion.[11]

In solid state electrolytes the high concentration of diffusing atoms, 31 mol/L in β–Li3PS4, can lead to complex interactions and diffusion behavior. Diffusion can involve collective jumps[12] and lattice vibrations,[10,13] which are all included in molecular dynamics (MD) simulations by taking into account all possible motions of ions and their interactions. Furthermore, MD simulations can show unanticipated diffusion behavior,[14] whereas static calculations (e.g., nudged elastic band) are limited by the imagination of the researcher. To understand diffusion in solid-state electrolytes MD simulations are thus a powerful tool, allowing the dynamic diffusion processes to be studied in detail.

Although MD simulations have been shown to provide understanding of complicated diffusion processes, however, typically the only property that is extracted is the tracer diffusivity, from which the activation energy is calculated by assuming Arrhenius behavior. A thorough analysis of MD simulations is able to give much more detailed results, potentially providing more understanding and concrete direction toward the design of improved conductivities.[10,15,16] To make such thorough analysis of MD simulations more easily available we present an approach, here demonstrated for β-Li3PS4, that allows to extract the detailed diffusional properties based on a MD simulation and the crystalline structure of the studied material. The approach determines jump rates, activation energies of different jumps, attempt frequency, vibrational amplitude, radial distribution functions, possible collective motions, site occupancies, tracer diffusivity, and the correlation factor.

The first part of this paper describes the approach that is followed to obtain the diffusional properties from a single MD simulation. In the second part MD simulations on β-Li3PS4 are analyzed, exemplifying how the developed approach helps in understanding diffusion in solid state electrolytes, and how this provides direction to design new and improved solid electrolyte materials. The Matlab code used for the analysis of MD simulations is freely available online.[17]

Information from MD Simulations

After performing a MD simulation the position of all the atoms at every time step is known. Typically, this result is used to determine the tracer diffusivity (D*) via the mean squared displacement:[18−21]where r(t + t0) is the displacement of a single atom with respect to the starting position (r(t0)), t the simulated time, N the number of diffusing atoms, and d the number of diffusion dimensions. Provided that the atomic displacement is significantly larger than the vibration amplitude this gives reliable values for the tracer diffusion, although the tracer diffusivity is only an approximation of the ionic diffusion, and taking ion correlations and the displacement of the center of mass into account leads to more precise results.[22]

Using the diffusivity and the Nernst–Einstein relation, assuming that the Haven ratio is equal to one,[22] the conductivity (σ) can be approximated[18,22] bywhere n is the diffusing particle density, e the elementary electron charge, z the ionic charge, kB Boltzmann’s constant, and T the temperature in Kelvin. The ionic conductivity determined in this way gives a good indication whether the material has an ionic conductivity which is high enough to make it suitable as a solid electrolyte.

However, to get a thorough understanding of the ionic diffusion in the material much more properties related to the diffusion process can be obtained from a single MD simulation, including:

amplitude of vibrations

attempt frequency

site occupations

jump rates

correlation factor

activation energies

collective jumps

radial distribution functions

Amplitude of Vibrations

Atomic vibrations in a crystal are the “back and forth” movement of an atom around a (meta)stable position. From a MD simulation the position of all atoms is known at any time, hence the absolute displacement can be obtained per time step. The atomic vibrations can be determined by monitoring the derivative of the absolute displacement, a change in the direction of movement gives a change in sign of this derivative, and atomic vibrations can thus be monitored. The vibrational amplitudes (A) are obtained by determining the change in absolute displacement (r) of each atom (i) while the derivative of the displacement keeps the same sign (between time steps ta and tb):Doing this for all the atoms of interest throughout the MD simulation gives a distribution of vibrational amplitudes, an example of which is shown in Figure . By fitting a Gaussian function to the obtained distribution the standard deviation in vibrational displacement can be obtained, providing an estimate of the average amplitude of vibrations in the crystal. Because the 3D distribution is known, anisotropic vibrational amplitudes can also be determined from a MD simulation in this way.

Figure 1

Histogram showing the vibrational amplitude of Li ions in β-Li3PS4 at 600 K, with the fitted Gaussian (solid red line) and the standard deviation (±0.495 Å, dotted green line).

As shown in Figure , there are some vibrations with amplitudes above 1 Å, which are probably caused by ionic jumps between different crystallographic sites. However, these large amplitude vibrations are only a small percentage of the total number of vibrations and will not significantly influence the Gaussian fit.

Attempt Frequency

The derivative of the absolute displacement can also be used to determine the vibration time of an atom or, via a Fourier transformation, the vibration frequency. The vibrational spectrum is obtained by determining the derivative of the absolute displacement (Δr) per atom (i) at every time step (t):which is transformed to a frequency spectrum via a Fourier transformation. The frequency spectrum of each (diffusing) atom is then combined to obtain the vibrational spectrum, as shown for Li ions in β-Li3PS4 in Figure . From the frequency spectrum the attempt frequency can be obtained if we consider every vibration of a diffusing atom as an attempt. Thus, we define the average vibration frequency as the attempt frequency (ν*), which is necessary to determine the relation between jump rates and activation energies.

Figure 2

Vibration frequency spectrum of Li-ions in β-Li3PS4 at 600 K, the average frequency of 8.29 (±0.46) × 1012 Hz is shown by the solid (±dotted) red line.

A high attempt frequency does not necessarily lead to to a high diffusivity, as high attempt frequencies are often related to a high activation energy,[23] a phenomenon which is known as the Meyer–Neldel rule.[24] Note that the presented approach determines a frequency spectrum based on the vibrations of single atoms. The obtained spectrum thus differs from a phonon spectrum, which depends on the collective vibrations of the atoms. To test the robustness of the presented approach, we have also determined the attempt frequency with a high and a low cutoff frequency; as shown in the Supporting Information, both have a negligible effect on the average attempt frequency.

The approach of obtaining the attempt frequency presented here is very different from the usual approaches. Several definitions of the attempt frequency exist,[25] all of which require the determination of the transition state and calculation of the phonon spectrum for the stable and the transition state.[26] However, often these calculations are not performed, and a “standard value” of 1 × 1013 Hz is used.[25,27]

In comparison to other methods the method presented here is straightforward, by using only the information that is already present in a MD simulation. Furthermore, because the attempt frequency is obtained from a single MD simulation, the influence of temperature, structural parameters, etc., on the attempt frequency are included and can be investigated.

Site Occupancy

In crystalline ionic conductors diffusion occurs through transitions between relatively stable sites. Typically these crystallographic sites are known from diffraction experiments, but if these are not known the sites can be extracted from a MD simulation through data mining.[28] The condition used for site occupancy is that the distance of the ion to the center of the crystallographic site is smaller than the site radius. At present the site radius is defined as twice the vibrational amplitude.

A defined site radius can lead to atoms that are not at a defined site, but by using the system and temperature dependent vibrational amplitude as a measure of site radius this effect is minimized. Instead, this can be used to determine what happens while an atom jumps between sites, which can be used to obtain information about the transition state. Furthermore, an “undefined” region between sites prevents the counting of small movements at site borders as jumps between sites, thus giving more reliable jump rates.

The site radius can also lead to overlapping sites, but this is accounted for in the analysis code. If two sites show overlap the site radius (of all sites) is reduced to half the distance between the two sites, thus preventing overlapping sites.

Jump Rates

When the crystallographic sites are known, detecting the transitions between sites that occur in a MD simulation is straightforward. Counting the number of jumps (J) between (types of) sites provides the mean jump rate (Γ) usingwhere N is the number of diffusing atoms, the subscript i defines a type of jump, and t is the simulation time. To get an estimate of the uncertainty in the jump rate the MD simulation is divided into ten different parts. By assuming that the correlation between consecutive parts is negligible, the differing jump rates in each part allow for the calculation of the standard deviation.

As demonstrated recently[15] determination of the different jump rates in a crystal provides direct insight in which jump process is rate-limiting for diffusion. This information can be used to design crystal structures with larger atomic diffusivity. Because NMR relaxation experiments can directly probe the jump rates, comparison with the jump rates from MD simulations can be used to validate the MD simulations,[29] or to better understand the complex results from NMR experiments.[30]

Using the Einstein–Smulochowski relation the jump rates are related to the jump rate diffusivity (DJ):where i are the different types of jumps, a is the jump distance of jump type i, and d the number of diffusion dimensions. The sum is over all types of jumps, since in most solid electrolytes several different types of jumps contribute to macroscopic diffusion.

In comparison to tracer diffusivity (D*) the jump rate diffusivity is usually an overestimate. The overestimation is caused by back-and-forth jumps, which cancel each other in the tracer diffusivity but are both counted in the jump diffusivity, and by the angles between consecutive jumps,[31] which causes the displacement to be lower as the total jump distance. To get an estimate of how effectively jumps contribute to macroscopic diffusion the correlation factor (f) can be calculated:[27,31]The correlation factor can be used to determine the diffusion mechanism (under certain conditions).[31]

Activation Energies

A simple way to describe the temperature dependence of diffusion is the activation energy, which is usually obtained by fitting an Arrhenius equation to diffusion data at various temperatures. This assumes that Arrhenius behavior is obeyed over the fitted temperature range, which assumes that there is no change in the material properties which determine diffusion over the studied temperature range. However, this is often incorrect,[32] especially when the studied temperature range is large.

Non-Arrhenius behavior usually leads to an underestimate of the activation energies at room temperature if extrapolated values from high temperatures are used. It is hard to say how large the errors are, since this will depend on the material[32] and the temperature range of the extrapolation.

To overcome this problem, we can calculate the activation energy at a given temperature, on the basis of which the temperature dependence of the activation energy can be determined via MD simulations at various temperatures. The activation energy (ΔE) for a type of jump (i) can be determined by the jump rates in a MD simulation via:[33]where kB is Boltzmann’s constant, T the temperature in Kelvin, Γeff the effective jump frequency, and ν* the attempt frequency, which is obtained using the approach described in the Attempt frequency section and assumed to be isotropic.

Entropy effects are naturally included when extracting the activation energy from a MD simulation because the temperature is above zero Kelvin, the obtained activation energy thus includes the activation enthalpy and entropy of the ionic jump, in contrast to activation energies from NEB or bond valence calculations. Since the activation energy is determined by the ratio between the effective jump frequency and attempt frequency, the activation energy can be seen as a measure for the jump probability: if the activation energy is low, the jump probability is high, and vice versa.

As shown in Figure , a jump from site A to B can have a different energy barrier as the reverse jump due to the difference in site energy, even though the number of A–B and B–A jumps will be the same in equilibrium. Because the number of jumps is equal, simply using eq to determine the jump frequency would give the same activation energy for A–B and B–A jumps, which is clearly incorrect.

Figure 3

Energy landscape and the corresponding ion density, with sites A, B, and C.

The exponential term in eq means that the activation energy is determined by the amount of successful jumps divided by the number of attempted jumps. Thus, the amount of time an atom spends at a certain site should be taken into account to correctly determine the activation energy. The effective jump rate, Γeff, thus differs from the jump rate in eq by taking into account the fraction of time that the diffusing atoms occupy a type of site (oj):where t is the total simulated time, and N the number of diffusing atoms. By incorporating site occupancy the sites with lower occupancy will have lower activation energies, correctly representing the difference in activation energy between A–B and B–A jumps in Figure . Furthermore, the difference in activation energy between back and forth jumps provides the energy difference between two sites, which can be used to predict changing site occupancies with temperature. Since jump and attempt frequencies are temperature dependent the activation energy may also be a function of temperature. Such non-Arrhenius behavior can be investigated by performing MD simulations at different temperatures.

Collective Jumps

The large concentration of diffusing atoms in solid electrolytes, 31 mol/L in β-Li3PS4, is likely to result in interactions between the diffusing Li atoms. This potentially causes collective jump processes,[12] which may have a severe impact on macroscopic diffusion.[34,35]

Due to the high ion concentration and rapid movement of ions in a solid electrolyte the energy landscape experienced by an ion changes continuously,[36,37] due to its own movements and the movements of neighboring ions. As shown schematically in Figure , the Coulombic repulsion of neighboring atoms can make energy barriers disappear, or new local energy minima can be created. In this way, the movement of one ion can cause the movement of a second ion, thus causing collective jumps.

Figure 4

Schematic pictures of the energy landscape due to the crystal lattice (top) and the Coulomb repulsion of another Li ion (middle) during a collective jump, giving the energy landscape experienced by the gray Li ion (bottom). Step 1: situation before a collective jump. Step 2: the other ion moves to the right and makes the barrier between the gray ion and the next site disappear. Step 3: situation after the collective jump.

MD simulations are a powerful tool for investigating complicated collective jump processes (the process that we call “collective” jumps is also referred to as “correlated” or “concerted” jumps in the literature). Knowing the position and time of jumps allows to determine if jumps are correlated in time and space, although determining the atomic process behind a correlated jump will require further study.

At present, the spatial condition for correlated motion is assumed to be slightly larger as the largest jump distance between Li sites in the crystal. For β-Li3PS4 this is 4.5 Å, in which 4 Å between the bc planes is the largest jump distance.

For collective motions to occur, atoms should move in the same direction within a certain time. Because atoms change their direction after a vibration, it is unlikely that atoms still show collective behavior after a large number of vibrations. The time scale for collective motions should thus be on the order of an average atomic vibration, which is equal to the period of the attempt frequency ( seconds).

A minimal condition for collective motion thus is that two jumps occur within the time required for a single vibration, and therefore we use a time of seconds. This time condition thus gives a lower limit of the number of collective jumps, but this is enough for the present purpose of determining whether collective motions are of importance for diffusion. Clearly, the conditions that define transitions as collective are debatable and should be chosen carefully for each material.

Radial Distribution Functions

The atomic environment determines the forces and energy barriers that govern the behavior of diffusing atoms. The Radial Distribution Function (RDF) can be used to reveal the density (g) of an element versus distance (r) with respect to another element during the MD simulation:where A is the number of atoms of the element at the center of the RDF, B the number of atoms of the other element, ntot the total amount of simulation steps, and r(n) the distance between atom x and y at time step n.

Because the position of all atoms is available at all time steps of a MD simulation, RDFs can be readily obtained for any site, atom or element. For example, this has shown to be useful for Na3PS4, where the RDFs from MD simulations suggested that Na vacancies are essential for increased Na diffusion.[10]

Summary

Summarizing, if the crystallographic sites are known, detailed diffusion properties can be extracted from MD simulations. In principle a single MD simulation, in which each type of jump occurs a significant number of times, already provides detailed insight into diffusion, including the diffusion path, attempt frequency, jump rates, activation energies, collective motions, atomic environments, and the correlation factor.

For a thorough understanding MD simulations at several temperatures might be necessary, for instance in the case of non-Arrhenius behavior, or to investigate the reliability of the results over a range of temperatures. Extracting the described information is beneficial for understanding of the diffusion process, allowing for a targeted approach to design and prepare new materials with enhanced properties, as will be demonstrated in this paper for the Li-ion conductor β-Li3PS4.

Example: β-Li3PS4

Li3PS4 has been a well-known Li-ion conductor since the 1980s,[38] but interest grew after experiments with nanosized crystals showed a Li-ion conductivity of 1.6 × 10–4 S/cm,[39] approaching the value that is required for solid state Li-ion batteries. Three polymorphs of Li3PS4 have been reported,[40] the low-temperature γ-phase, the β-phase at intermediate temperatures, and the high temperature α-phase. The β-phase shows the highest room temperature conductivity of the three polymorphs,[39] and is thus most interesting for application as a solid electrolyte. A beneficial property of β-Li3PS4 is its apparent stability against Li metal,[39] although DFT calculations report otherwise.[7,41] Li3PS4 can be prepared via a solvent route,[42,43] resulting in a conductivity of 3.3 × 10–4 S/cm,[42] enabling coating of cathode materials. In this way, no additional solid electrolyte material needs to be added in the cathodic mixture,[43] resulting in a larger effective energy density in combination with a small interface resistance.

The β-phase of Li3PS4 is reported to crystallize in the Pnma space group (no. 62)[38,40,44] in which Li ions occupy 4b, 4c, and 8d positions. Studies[38,40,44] investigating the structure of β-Li3PS4 report significantly different Li-ion positions and occupancies. Neutron diffraction[44] indicates that the coordinates of the Li-ion 4c position strongly depend on temperature, potentially explaining the differences between X-ray diffraction studies.[38,40]

Based on the larger sensitivity to Li ions of neutrons compared to X-rays, the Li positions determined from neutron diffraction[44] at 413 K are used for the analysis of the present MD simulations on β-Li3PS4.

Effect of Li Vacancies and Li Interstitials

The introduction of Li vacancies has been suggested to be beneficial for Li-ion conductivity in β-Li3PS4,[13] whereas the high ionic conductivity of the isostructural[45] compound Li10GeP2S12 (= Li3.33Ge0.33P0.67S4) suggests that introducing extra Li ions in β-Li3PS4 can also lead to an increased Li-ion conductivity. To study the effect of both Li vacancies and Li interstitials on the diffusion mechanism DFT MD simulations were performed for β-Li3PS4, β-Li2.75PS4, and β-Li3.25PS4 at 450, 600, and 750 K.

The compositions of β-Li2.75PS4 and β-Li3.25PS4 were chosen because this leads to two Li vacancies/interstitials in the super cell. The two Li vacancies/interstitials were divided over the two bc planes present in the used super cell, in which the Li ions show fast diffusion.

The introduction of Li interstitials leads to a significant increase in the occupancy of 4c sites, as shown in Figure S2 and reported by Lepley et al.[46]

Li-Ion Diffusion

The diffusion paths from simulations at 600 K, shown in Figure , demonstrate that diffusion takes place along the b axis via 4b–4c jumps, via intraplane jumps along the c axis through 4b–8d and 4c–8d jumps, and 8d–8d jumps in the a direction are responsible for interplane jumps. As shown in Figure b in stoichiometric β-Li3PS4 the most jumps occur along the b axis, followed by jumps in the bc planes, and relatively few transitions occur between the bc planes. This indicates that Li-ion diffusion occurs primarily within the bc planes. The MSD shown in Figure S1 shows a different picture, with the b direction having the largest MSD, followed by the c and a direction. This indicates that the correlation factor (eq ) is lower for faster diffusion directions, similar to the results of Marcolongo et al.[22] in LGPS.

Figure 5

Jump diffusion paths at 600 K for (a) β-Li2.75PS4, (b) β-Li3PS4, and (c) β-Li3.25PS4. Li-ion sites are shown by 4b = blue, 4c = green, and 8d = black. The jump types are shown by 4b–4c = red, intralayer = pink, and interlayer = cyan, thicker lines correspond to larger jump rates.

When Li vacancies or Li interstitials are introduced the Li-ion diffusion within the bc plane remains similar to the stoichiometric composition β-Li3PS4, whereas the amount of interplane jumps increases significantly, resulting in three-dimensional diffusion.

The beneficial effect of the three-dimensional diffusion is reflected in the tracer diffusivity, shown in Figure . The Li-ion diffusivity in β-Li2.75PS4 is almost an order of magnitude larger than in β-Li3PS4. Introducing Li interstitials by creating β-Li3.25PS4 also results in a larger diffusivity, especially at the lowest simulated temperature.

Figure 6

Tracer diffusivity from the current MD simulations, Phani et al.[13] and Yang et al.[49]

On the basis of the tracer diffusivity the conductivity of β-Li3PS4 is 1 × 10–2 S/cm at 450 K, comparable to impedance experiments[42] at the same temperature. Extrapolating the Li-ion diffusivity of β-Li3PS4 to 110 °C results in a Li-ion diffusivity of 1 × 10–8 cm2/s, close to the values reported by NMR experiments: 3.0 × 10–8 cm2/s at 100 °C[47] and between 1 × 10–6 and 1 × 10–8 cm2/s at 120 °C.[48]

The results from the current MD simulations on β-Li3PS4 are comparable to the values reported previously,[13,49] except at 450 K. This anomaly is most likely caused by the shorter simulation times of the previous studies, which can lead to an overestimation of the tracer diffusion,[21] especially at low temperatures. At 750 K all the MD simulations show a similar value for the diffusivity, which can be explained by the melted lithium sublattice[13] at this temperature. After melting the lithium ordering over the different crystallographic sites disappears, which seems to have a larger impact as the deviating stoichiometries investigated here. The detrimental effect of Li ordering has also been reported for the solid electrolyte Li7La3Zr2O12, in which Li ordering can reduce the Li conductivity by several orders of magnitude.[50]

Jump Rates

The differences in tracer diffusivities between the three compositions can be explained by the rate-limiting jump mechanism. The most frequent jump process is the 4b–4c transition, the rate of which is comparable between the three compositions, as shown in Figure . However, to obtain three-dimensional diffusion paths in β-Li3PS4 interplane jumps are necessary, the rate of which is significantly different for the three compositions, also shown in Figure . With lower temperature, these differences increase, and in β-Li3PS4 at 450 K, no interplane jumps occurred.

Figure 7

Jump rates for the 4b–4c and interplane jumps.

Because two-dimensional diffusion processes have a smaller correlation factor compared to three-dimensional processes[18,31,51] the tracer diffusivity in β-Li3PS4 is significantly lower, even though the jump rate of the fastest diffusion process is similar.

Activation Energies

The activation barriers for 4b–4c and intraplane jumps obtained from the MD simulations are shown in Figure . At 600 K the interplane 8d–8d jumps show activation energies of 0.40 ± 0.02 eV for β-Li3PS4, 0.31 ± 0.06 eV for β-Li3.25PS4, and 0.29 ± 0.02 eV for β-Li2.75PS4. It should be noted that other jump processes also occur in the simulations, however, their significantly larger activation energy indicates that these will not contribute significantly to Li-ion diffusion and are therefore left out of the current analysis.

Figure 8

Activation energy at 600 K for (a) 4b–4c and (b) intraplane jumps.

Over the simulated temperature range some activation energies remain constant, whereas others show significant changes, as shown in Figures S4 and S5. The changes in activation energies could be caused by non-Arrhenius behavior of some types of jumps, but the loss of Li ordering at 750 K might also be a reason for the changing activation energies. Furthermore, the different compositions sometimes show differing behavior for the same type of jump. To say something about the (non-)Arrhenius behavior of β-Li3PS4 further study is thus required.

To validate the activation energies from MD simulations comparison with experimental values is important, however, a wide distribution in values is reported based on electrochemical experiments: 0.16,[40] 0.32,[42] 0.36,[39] and 0.47[43] eV. NMR experiments resulted in activation energies of 0.40 eV for macroscopic diffusion and 0.09 eV for local jumps.[48] Given this wide distribution of values a comparison of experimental activation energies with the present simulations seems unreasonable.

Simulations also report a wide range of activation energies. NEB calculations on β-Li3PS4 report activation energies of 0.3 eV along the a axis and 0.2 eV along the b and c axis,[46,52] whereas other NEB calculations[53] report 0.26 eV along the a and b axes and 0.08 eV for collective Li-ion jumps in the b direction, bond-valence calculations report values of 1.0 eV along the a axis and 0.8 eV in the bc plane,[54] and an Arrhenius fit to MD simulations[13] gives an activation energy of 0.35 eV. The results from NEB calculations and MD simulations are comparable, whereas bond-valence calculations appear to overestimate the activation energy. The activation energies from the present MD simulations indicate that diffusion along the b axis is most facile, followed by diffusion along the c axis, and along the a axis diffusion is most difficult, in agreement with results from neutron diffraction.[44]

Collective Jump Processes

Given the large lithium concentration of 31 mol/L in β-Li3PS4, Li ions can be expected to interact strongly with each other. Yang et al.[53] reported the presence of collective jumps in β-Li3PS4 by showing that the activation energy for diffusion along the b axis is just 0.08 eV for two Li ions moving collectively, whereas it is 0.26 eV for a single Li ion.

To determine the importance of collective jumps in the present MD simulations, we have determined the amount of collective jumps as described in the Collective Jumps section. Comparison of the collective jumps with the total amount of jumps reveals that the percentage of collective jumps strongly depends on the temperature and Li concentration, as shown in Figure . The percentage of collective jumps displays a strong increase with temperature, with 65–80% of the jumps occurring collectively at 750 K. Although at 450 K the simulations show fewer collective jumps, still 24% of the jumps are collective in Li2.75PS4.

Figure 9

Percentage of collective jumps in the MD simulations

Surprisingly, the lower Li concentration in Li2.75PS4 leads to a higher percentage of collective jumps. As shown in Figure , the jump rate is similar in all three compositions, indicating that collective jumps occur more frequently when there are more Li vacancies, but the mechanism behind this effect is unclear.

To study the effect of the time condition on the determination of collective jumps, we varied the time condition between 0.25 and 5 times the attempt frequency, the results of which are shown in Figure .

Figure 10

Number of collective jumps versus the time condition (in units of the attempt frequency).

This shows that the amount of collective jumps determined using the presented method strongly depends on the specified time condition. If the collective jumps would purely be due to coincidences a linear relationship between the number of collective jumps and the time condition is expected, but this is not observed in Figure . When using a shorter time condition, the number of collective jumps increases more rapidly as when using longer time conditions, indicating that collective jumps are occurring. However, the strong dependence of this analysis on the time condition specified for collective jumps shows that quantifying the amount of collective jumps will require a more thorough understanding of collective behavior before conclusions can be drawn. However, the large percentage of collective jumps indicates that collective jump processes may have a significant effect on the Li-ion diffusion in β-Li3PS4, especially at elevated temperatures. Although the larger amount of collective jumps is also caused by the higher amount of jumps at elevated temperatures, the movements of a Li ion will be influenced by jumps of other Li ions nearby, and increasing the temperature thus increases the amount of collective motion.

In Figure S6, the number of collective jumps are shown per combination of jump type during the MD simulation at 600 K. This shows that in β-Li3PS4 the collective jumps are primarily simultaneous 4c–4b jumps and simultaneous 4b–4c jumps. In Li3.25PS4, collective 4b–4c and 4c–4b jumps also occur most frequently, but additionally, 4b–8d jumps collective with 4b–4c jumps occur and collective interplane jumps are predicted. In the simulations of Li2.75PS4 different collective behavior is observed, with the combination of 4b–8d jumps and 4b–4c jumps occurring most frequently. Collective 4b–4c jumps also occur frequently, but significantly less compared to the other compositions.

Collective jumps involving more than two Li ions also occur in the current MD simulations, in some cases involving up to 5 atoms. The collective movement of multiple atoms is complex and difficult to analyze. However, it should be anticipated that collective motion of several ions induces large ionic conductivities, as observed in Li10GeP2S12,[35] making this an interesting and relevant research area.

Attempt Frequency

The attempt frequencies obtained from the MD simulations are shown in Figure a. Using the presented method of obtaining attempt frequencies results in values between 7.6 × 1012 and 9.3 × 1012 Hz, comparable to reported attempt frequencies, between 1 × 1012 and 1 × 1013 Hz, for other materials by experiments[23,55] and DFT calculations.[25,26] This consistency indicates that the presented method of determining the attempt frequency from MD simulations, by a straightforward Fourier transformation of the ionic velocity, can be used to determine the attempt frequency.

Figure 11

(a) Attempt frequencies and (b) vibration amplitudes from the MD simulations.

Figure a demonstrates that the attempt frequency decreases with increasing temperature, this decrease with increasing temperature can be understood by the change in the vibrational amplitude, shown in Figure b. For Li3.25PS4 and Li3PS4 the vibrational amplitude increases by approximately 35% between 450 and 750 K, whereas the speed of the atoms increases by just 29% (based on: ). The average vibration time therefore increases, leading to a decreasing attempt frequency. Furthermore, as shown in Figure the attempt frequency and vibration amplitude in Li2.75PS4 differ by approximately 10% from the other two compositions. This effect can be explained by the lower Li concentration, which leads to less repulsive interactions between Li ions, thus allowing for bigger vibrations and a lower attempt frequency.

This example demonstrates that there can be a significant temperature dependence on the attempt frequency and that relatively small changes in the crystal structure can have a significant effect on the attempt frequency. Thereby the presented analysis indicates that consideration of the attempt frequency and its dependence on structure and temperature is of significant importance in quantifying and understanding ionic diffusion.

Effects of Doping

The MD simulations on β-Li3.25PS4 show that creating Li interstitials in β-Li3PS4 significantly increases three-dimensional diffusion, raising the macroscopic Li-ion conductivity, in line with experimental work.[45] Li vacancies also increase three-dimensional diffusion, and the MD simulations indicate significantly larger Li-ion conductivity compared to the introduction of Li interstitials. However, we are unaware of work which has explored the impact of Li vacancies on Li-ion diffusion in β-Li3PS4. To determine the impact of introducing Li vacancies by doping, we performed MD simulations on β-Li2.75PS3.75Br0.25.

Additionally, the impact of oxygen doping is investigated, because this is also reported as a strategy to improve the Li-ion diffusivity.[47,52,54] However, O doping does not change the Li content. It has been proposed that the smaller radius of O2– in comparison to S2– opens up a diffusion path along the a axis, but the mechanism of the oxygen induced larger Li-ion diffusivity remains unclear. To gain further understanding of how oxygen doping increases the Li-ion diffusivity, we performed MD simulations on β-Li3PS3.75O0.25.

The structures with dopants were created by replacing two S atoms by Br and O atoms, respectively. The dopants were placed the maximum possible distance apart from each other in the supercell, to minimize the interaction between dopant atoms. The introduction of the dopant atoms causes only minor changes in the lattice parameters, less than 2% in both cases. The dopants only cause local distortions in the crystal, as shown by the large similarity in the RDF’s of the S atoms in the doped and nondoped structures in Figure S3, which is consistent with previous calculations on O-doped β-Li3PS4.[52]

Predicting the stability of solid electrolytes is hard, as shown by the varying decomposition voltages and decomposition products between papers which calculated the stabilities of various solid state electrolytes,[7,41] and it becomes even more complicated when including dopants. In this paper we focus on which diffusion properties can be determined from MD simulations on solid electrolytes, and determination of the stability of the doped structures is thus beyond the scope of this paper.

However, for the O-doped β-Li3PS4 it has been shown that the addition of oxygen is beneficial for the stability.[52] And the Na analogue (Na3PS4) of Li3PS4 with halogen doping has been made experimentally[5] and is thus shown to be (meta)stable, and we expect that β-Li4PS4 with Br doping is also (meta)stable.

Br Doping

The jump diffusion path from the MD simulation of β-Li2.75PS3.75Br0.25 is shown in Figure a. As should be anticipated from the simulations of β-Li2.75PS4, Figure a shows that Br doping significantly increases the amount of three-dimensional diffusion. In β-Li2.75PS3.75Br0.25, the tracer diffusivity results in 1.56 × 10–6, 1.01 × 10–5, and 3.71 × 10–5 cm2/s at 450, 600, and 750 K, respectively, comparable to the diffusivities of β-Li2.75PS4. The activation energies for diffusion along the bc plane in the Br-doped composition at 600 K, shown in Figure , differ by just 0.02 eV from the MD simulation with Li vacancies. The activation energy for interplane jumps is 0.29 ± 0.03 eV in β-Li2.75PS3.75Br0.25 and 0.29 ± 0.02 eV in β-Li2.75PS4.

Figure 12

Jump diffusion paths at 600 K for (a) β-Li2.75PS3.75Br0.25 and (b) β-Li3PS3.75O0.25. Li-ion sites are shown by 4b = blue, 4c = green, and 8d = black. The jump types are shown by 4b–4c = red, intralayer = pink, and interlayer = cyan, thicker lines correspond to larger jump rates.

Figure 13

Activation energies at 600 K in Br- and O-doped β-Li3PS4 for (a) 4b–4c and (b) intraplane jumps.

The similar activation energies in β-Li2.75PS4 and β-Li2.75PS3.75Br0.25 indicates that the primary cause of the high Li-ion diffusivity in Br-doped β-Li3PS4 are the Li vacancies. For other dopants that introduce Li vacancies similar results are thus expected, suggesting that there are various ways to increase the Li-ion diffusivity in β-Li3PS4.

O Doping

The jump diffusion paths from the present MD simulations on β-Li3PS3.75O0.25, shown in Figure b, demonstrate that O doping also increases three-dimensional diffusion, as was predicted.[52,54] Compared to β-Li3PS4 doping with oxygen leads to a larger Li-ion diffusivity, and tracer diffusivities of 6.96 × 10–7, 7.59 × 10–6, and 2.52 × 10–5 cm2/sec at 450, 600, and 750 K, respectively. The introduction of oxygen results in smaller activation energies compared to β-Li3PS4, as shown in Figure . The biggest impact is observed for the interplane jumps, which have an activation energy of just 0.31 ± 0.03 eV in β-Li3PS3.75O0.25. NEB calculations on β-Li3PS4 with O doping[52] show activation energies of approximately 0.2 eV in the b and c direction, which is comparable to the results presented in Figure . In the a direction NEB calculations show that the activation energy is 0.2 eV close to the O dopant, and 0.3 eV far away from the O dopant. From the MD simulations no distinction is made between jumps in the a direction close or far from the O dopant, but the 0.31 ± 0.03 eV is close to the value reported by the NEB calculations.

It is surprising that the activation energies for the O-doped structure are comparable to the Li-rich β-Li3.25PS4, even though the introduction of oxygen does not affect the Li concentration. To investigate this, we calculated the radial distribution functions (RDFs) for oxygen and sulfur in β-Li3PS3.75O0.25 using eq and shown in Figure . The smaller ionic radius of oxygen[56] results in a smaller O–Li distance compared to the S–Li distance. However, the Li density of the first coordination shell in the RDF is significantly lower around the O atoms. Integrating the Li density up to 3.5 Å shows that (on average) there are 2.9 Li atoms in the first coordination shell of O atoms, and 3.5 Li atoms in the first coordination shell of S atoms.

Figure 14

O–Li and S–Li distribution in β-Li3PS3.75O0.25 at 600 K, calculated using eq .

This implies that it is unfavorable for Li ions to be near the O atoms in β-Li3PS3.75O0.25, and these Li ions must be accommodated elsewhere within the crystal structure. Oxygen doping thus changes the distribution of Li ions in the crystal, which is shown to be beneficial for the Li-ion diffusivity in β–Li3PS4.

It is usually assumed that a higher polarizability leads to higher diffusivity[57] via lower activation energies caused by lattice softening.[23] In the case of β-Li3PS3.75O0.25 the higher diffusivity caused by O atoms, which have lower polarizability compared to S atoms, demonstrates that the opposite can also occur. Integrating the RDFs in Figure shows that there are fewer Li ions near the O atoms as near the S atoms, which indicates that the site energy near the O atoms is higher. If the transitions state energy stays the same this lowers the activation energy. In this case, the less-polarizable O atoms are thus beneficial for Li-ion diffusion.

Conclusions

The present approach demonstrates that from a single MD simulation key properties for ionic diffusion can be obtained, through which a thorough understanding of diffusion can be developed. The thorough analysis of MD simulations presented is a general approach that can be applied to all crystalline ionic conductors, which can help to build understanding of diffusional processes in solid state electrolytes, and provide direction to the design of new and improved solid electrolyte materials. The Matlab code developed for the analysis of MD simulations is freely available online.[17] The example of DFT MD simulations on β-Li3PS4 indicates that Li-ion jumps between bc layers limit the macroscopic conductivity. Adding Li interstitials or Li vacancies significantly promotes these transitions, increasing macroscopic Li-ion diffusion. Li vacancies can be introduced through Br doping at the sulfur sites, which is predicted to result in an order-of-magnitude larger Li-ion conductivity in β-Li3PS4. Furthermore, it is shown that oxygen doping at the sulfur site changes the Li distribution in the crystal, rationalizing the increased Li-ion diffusivity that has been reported.

DFT Calculations

The DFT simulations were performed using VASP,[58] using the GGA approximation[59] and the PAW–PBE basis set.[60] A cutoff energy of 400 eV was used for simulations containing oxygen, and 280 eV for the other simulations. The β-Li3PS4-phase crystallizes in the orthorhombic space group Pnma (no. 62), with lattice parameters of a = 12.82, b = 8.22, and c = 6.12 Å at 637 K.[40] The crystal structure as measured by Homma et al.[40] was used as a starting point for the structure minimizations. A 1 × 1 × 2 super cell was used in the calculations, testing of finite size effects was not done due to computational limitations. The fractional site occupancies reported experimentally were approximated as much as possible in the super cell, while maximizing the Li–Li distances. After minimization without symmetry restrictions, the lattice angles were close to 90° in all cases, and the lattice parameters changed by less as 2%, with the a and c parameters showing a small increase, whereas the b parameter decreased slightly. During the minimizations, a k-point mesh of 4 × 6 × 4 was used, which was reduced to a k-point mesh of 1 × 2 × 1 for the MD simulations. The total simulation time of the MD simulations was 500 ps, with 2 fs time steps. The first 2.5 ps were used as equilibration time and were thus not used for the analysis. Simulations were performed in the NVT ensemble, with temperature scaling after every 50 time steps.