Dennis Wiedemann, Mazharul M. Islam, Thomas Bredow and Martin Lerch

Diffusion Pathways and Activation Energies in Crystalline Lithium-Ion Conductors

De Gruyter | 2017

Abstract

Geometric information about ion migration (diffusion pathways) and knowledge about the associated energy landscape (migration activation barriers) are essential cornerstones for a comprehensive understanding of lithium transport in solids. Although many lithium-ion conductors are discussed, developed, and already used as energy-storage materials, fundamental knowledge is often still lacking. In this microreview, we give an introduction to the experimental and computational methods used in our subproject within the research unit FOR 1277, “Mobility of Lithium Ions in Solids (molife)”. These comprise, amongst others, neutron diffraction, topological analyses (procrystal-void analysis and Voronoi–Dirichlet partitioning), examination of scattering-length density maps reconstructed via maximum-entropy methods (MEM), analysis of probability-density functions (PDFs) and one-particle potentials (OPPs), as well as climbing-image nudged-elastic-band (cNEB) computations at density-functional theory (DFT) level. The results of our studies using these approaches on ternary lithium oxides and sulfides with different conduction characteristics (fast/slow) and dimensionalities (one-/two-/three-dimensional) are summarized, focusing on the close orbit of the research unit. Not only did the investigations elucidate the lithium-diffusion pathways and migration activation energies in the studied compounds, but we also established a versatile set of methods for the evaluation of data of differing quality.

1 Introduction

A variety of experimental methods is available to characterize the mobility of lithium ions in solids: e.g. electrical conduction phenomena and activation energies may be deduced from impedance spectroscopy; nuclear magnetic resonance (NMR) spectroscopy shows the dimensionalities of migration pathways and differences between crystallographically unique ion positions; tracer-diffusion experiments allow for the investigation of self-diffusion on a macroscopic scale. Yet, none of these provide direct information on the geometry of diffusion pathways.

In light of this shortcoming, variable-temperature neutron diffraction holds a unique position amongst experimental approaches. Because scattering takes place at the atomic nuclei, it is sensitive towards the otherwise (e.g. for X-ray techniques) notoriously elusive lithium ions. Furthermore, it provides a precise space- and time-averaged picture of the crystal structure, which is not influenced by electron-density shifts due to bonding, polarization, or other interactions. Diffusion-induced displacement is primarily coded in the intensity of the reflections at high diffraction angles. Data may be collected either from powder or single-crystal specimen to derive a crystal-structure model using the Rietveld method or refining on structure-factor magnitudes. Subsequently, several methodologies allow for the mapping of diffusion pathways, the assessment of their relative probability, and for the evaluation of associated effective energy barriers. In our recent work, we have focused on three approaches (named in order of increasing complexity): topological analyses, mapping of the scattering-length density (SLD) as reconstructed via maximum-entropy methods (MEM), and the evaluation of the probability-density function (PDF) as well as the effective one-particle potential (OPP).

These experimental techniques are ideally complemented by solid-state quantum-chemical methods. The present state of the art of theoretically modeling ion mobility in solids is based on density-functional theory (DFT) combined with periodic supercells. These methods allow for a pointwise calculation of the potential-energy hypersurface (PEHS) along a local lithium migration pathway, which is defined as the minimum-energy pathway (MEP). The MEP can be calculated within various approximations. The widely used nudged-elastic-band (NEB) method starts from a linear interpolation of the atomic positions between initial and final ion configuration [1], [2], so that both crystal structures (within the selected unit cell) must be known in advance. The MEP is calculated for a small number of images (usually three to seven) between these endpoints. For each image, the atomic positions of all or a selected set of atoms are relaxed under the restriction that the structures are equidistant with respect to the intrinsic reaction coordinate (IRC) connecting the initial and final configuration. This is achieved by projecting out the MEP tangential parts of the atomic forces in the energy minimization.

In the climbing-image NEB (cNEB) variant, this constraint is relaxed for the image with the highest energy to obtain a better approximation of the saddle point [3]. The PEHS around a known structure can be explored with the dimer method [4]. A search for a saddle point is initiated either in a random direction or along the normal vector with the smallest eigenvalue obtained for the initial minimum structure. This approach allows the detection of unexpected MEPs but is computationally much more demanding than the cNEB method, which is therefore preferred, in case the local migration process is known. Traditional transition-state (TS) optimization procedures based on the Newton–Raphson method tend to fail due to numerical instabilities if the initial structure is not close enough to that of the TS. Thus, a reasonable compromise is to combine NEB and Newton–Raphson methods by taking the image with the highest energy from a cNEB calculation as a starting structure for a TS optimization.

A completely different approach to simulate migration processes is the molecular dynamics (MD) technique at DFT level (aiMD). The calculated atomic forces are used to solve Newton’s equations of motion. The velocities and updated atomic positions are calculated numerically with discrete time steps, e.g. by the velocity Verlet algorithm [5]. In principle, the aiMD approach can find all possible migration processes for a given structure. However, since the activation barriers for lithium migration in solids are usually large, a jump process is a rare event and an extremely large number of time steps in the order of 105–106 has to be calculated, which limits the applicability for solids with large unit cells. It should be noted that such and even larger time scales are accessible with classical MD simulations (cMD), in which forces are calculated from empirical pair potentials. With cMD, the results are thus strongly dependent on the quality of the empirical parameters. On the other hand, due to advances in numerical-method development and computer hardware, aiMD simulations of solids have become an important tool in theoretical solid-state chemistry [6]. With the metadynamics technique [7], it is possible to accelerate MD simulations of rare events by several orders of magnitude. The method, however, requires knowledge about the final structure in order to define a generalized reaction coordinate that distinguishes the initial and final points. Metadynamics simulations provide the free energy of activation, but not the exact TS structure. The latter is only approximated by the snapshot with the highest energy and must be further refined as mentioned before.

In this microreview, we present a general introduction to our experimental and computational methods, as well as a selection of results acquired in the framework of FOR 1277: “Mobility of Lithium Ions in Solids (molife)” funded by Deutsche Forschungsgemeinschaft (DFG).

2 Methods of evaluation

2.1 Topological analyses

2.1.1 Common aspects

Topological methods divide ionic conductors into a static framework and mobile interstitial species via partitioning of the direct crystallographic space. They rely only on accurately determined framework-ion positions, while the interstitial containing the mobile ions is treated as void and classified by its ability to host or conduct them [8]. Additional prior knowledge (e.g. electron-density distributions from ab-initio calculations or ionic radii from statistical analyses/data mining) is inferred to warrant this classification. Topological methods are especially appealing because of their high speed, low computing cost, and presentiveness – be it as a heuristic tool to scan for possible migration pathways or as a standalone method for lower-quality data.

To understand the underlying problem, let us assume we want to topologically analyze a crystal of known structure. Chemical intuition tells us that this crystal is somehow “made up of” its components – be it molecules or small ions. But how exactly do we partition the crystallographic space, so that any point of the asymmetric unit is assigned to a component or a void? In other words, where exactly does one component end and another one or a void begin? The answer to this, at first glance, rather trivial question is that there is no “correct” or even single best way for partitioning. But two approaches have become particularly well-established for the investigation of fast ion conductors during the last 10 years: the Voronoi–Dirichlet partitioning (VDP) and the procrystal-void analysis (PVA). These methods are of great value

  • as supportive approaches, which are more powerful then crystal-chemical rules of thumb,

  • for unruly datasets, which do not allow for the use of high-end methods, and

  • as heuristic means to search for and categorize migration pathways.

When dealing with disorder, some heed has to be paid: most crystal structures are derived from diffraction data and thus represent a time and space average (becoming, e.g. manifest in sites with mixed or partial occupation). Topological analyses, on the other hand, deal with ordered structures and their local features. Because of this, ordered supercells with localized defects must be set up to reflect probable situations. Not only can this be tedious (high number of possibilities) or biased (arbitrary choice of models), but the mere substitution/deletion of an atom also introduces a systematic error. The configuration around the modeled defect will be energetically unfavorable because it is not locally “relaxed”, but a representation of the predominant situation in the whole crystal. Therefore, it must be kept in mind that topological analyses are somewhat more approximate when dealing with disordered systems.

2.1.2 Procrystal-void analysis

Although procrystal analysis was originally developed for the examination of Hirshfeld surfaces in molecular crystals [9], [10], it also performed well in the characterization of voids in metal-organic/covalent organic frameworks [11] and migration pathways in ionic solids [12]. In its course, the so-called “procrystal” is constructed from the known crystal structure by centering spherically averaged electron density ρ at ( r ) at the atomic positions. These densities only depend on the atom/ion present and are derived from ab-initio calculations. To investigate potential diffusion pathways, the set 𝕄 of mobile species is not included in the construction of the procrystal, thus leaving only the static framework 𝔽 and void space:

ρ pro ( r ) = i F ρ i at ( r ) = i ρ i at ( r ) i M ρ i at ( r )

The procrystal density ρ pro ( r ) is defined at any position r and conveniently given in Hartree atomic units (a0: Bohr radius) [13]:

[ ρ pro ( r ) ] = 1 a .u . = 1 a 0 3 = 6 .748334504(5) Å 3

Integrated facilities to calculate and visualize procrystal-density isosurfaces from crystal-structure data are available in the fairly state-of-the-art program crystalexplorer [14]. An isosurface for a given value (isovalue) will be a set of continuous closed surfaces, which may, however, seem to be capped at the boundaries for reasons of representation. The higher the isovalue is, the closer the surface is to the atomic centers. Below a certain threshold, there will be no associated surfaces, because the procrystal density is finite everywhere. So, when successively raising the isovalue from zero, a first isosurface will appear and show the regions with the lowest procrystal density: the largest voids (see Figure 1, left). The enclosed volume will then grow to show more features and eventually form continuous channels, layers, or networks: the most probable migration pathways (especially if they comprise the positions of the excluded mobile species; see Figure 1, center). More features – less probable pathways – may ensue before the isosurface closely covers the immobile framework or even everts to close around isolated molecular species (see Figure 1, right).

Fig. 1: Procrystal-void surfaces of the defective Ti3O72− framework in ramsdellite-like Li2Ti3O7 (orthorhombic, Pnma) at 422°C for ρpro=0.004 a.u. (left), 0.006 a.u. (center), and 0.02 a.u. (right) [15]. Gray: titanium, red: oxide ions (arbitrary radii), blue: procrystal-density isosurface (caps not drawn for clarity, inner side in gray), orange circles: bottleneck positions. View approximately along [010].

Fig. 1:

Procrystal-void surfaces of the defective Ti3O72− framework in ramsdellite-like Li2Ti3O7 (orthorhombic, Pnma) at 422°C for ρpro=0.004 a.u. (left), 0.006 a.u. (center), and 0.02 a.u. (right) [15]. Gray: titanium, red: oxide ions (arbitrary radii), blue: procrystal-density isosurface (caps not drawn for clarity, inner side in gray), orange circles: bottleneck positions. View approximately along [010].

In this way, it is possible to discriminate probable from less probable pathways. But the choice of the isovalue, at which an indicated pathway can be considered viable, may seem somewhat arbitrary at first glance. However, it has been shown that a value of 0.002 a.u. compares well to smoothed van-der-Waals surfaces of small to medium-sized molecules and their electron densities from ab-initio Hartree–Fock calculations (6-31G* level) [16]. An isovalue of 0.0003 a.u., on the other hand, is fitting for the investigation of porous materials [11] and lithium diffusion pathways are observed at 0.004–0.006 a.u. and 0.002–0.004 a.u. in layered and other structures, respectively [12]. (In our experience, poor to moderate lithium-ion conductors fall into the higher range, even if not layered).

According to the simplified picture of the PVA, the most probable pathway of migration between two points is the one of the lowest procrystal density. The position of the density maximum along it signifies the bottleneck of diffusion and the associated density (or its difference to the path’s minimum) should be related to the activation energy. Nevertheless, to our knowledge, a sound correlation has not yet been found [12].

2.1.3 Voronoi–Dirichlet partitioning analysis

The VDP analysis takes a completely different approach to the partitioning of space and the concept of voids. It is based on purely mathematical concepts first introduced into crystallography by Niggli in 1927 [15] and does not rely on prior knowledge (except for the crystal structure) in its initial steps. The idea behind it is simple: any point is assigned to the atomic center that is nearer than all the others. Also, the algorithm for construction is straightforward:

  • connect any center to the centers in its first coordination sphere with straight line segments,

  • construct a perpendicular bisecting plane for every line segment, and

  • cut the planes at intersections with other planes.

The resulting atomic Voronoi–Dirichlet polyhedra (see Figure 2, left) are centered at the atomic positions and tessellate crystallographic space. The centers of elementary voids are the points the furthest away from the at least four adjacent atomic centers – the vertices of the Voronoi–Dirichlet polyhedra. They are connected by straight elementary channels, which run as far away as possible from at least three adjacent atomic centers. Thus, they correspond to polyhedral edges (see Figure 2, center) [17]. But being “far away” is certainly not enough to classify as a crystallographic void structure, which can host or conduct an ion. Elementary voids/channels can be

Fig. 2: Details from the VDP analysis of the AlO2− framework in γ-LiAlO2 (tetragonal, P41212) at 72°C. Voronoi–Dirichlet polyhedra of the aluminum and oxide ions (left), initial void structure (center), conduction pattern after deletion of improbable voids and channels (right). Blue: lithium, gray: aluminum, red: oxide ions (Slater radii), black: elementary voids and channels (arbitrary radii). View approximately along [001].

Fig. 2:

Details from the VDP analysis of the AlO2 framework in γ-LiAlO2 (tetragonal, P41212) at 72°C. Voronoi–Dirichlet polyhedra of the aluminum and oxide ions (left), initial void structure (center), conduction pattern after deletion of improbable voids and channels (right). Blue: lithium, gray: aluminum, red: oxide ions (Slater radii), black: elementary voids and channels (arbitrary radii). View approximately along [001].

  • too small,

  • ill-shaped,

  • surrounded by like-charged ions,

  • too narrow,

  • isolated, or a dead-end.

At this point, additional prior knowledge must be inferred. The size and shape of an elementary void is assessed using its spherical-domain radius Rsd (radius of an isochoric sphere) and its second moment of inertia G3. Typical values required by known ions are tabulated [18]. The influence of like-charged ions is typically estimated from the solid angle Ω of interaction with them. Finally, the width of a channel is characterized via the adjacency radius Rad that should not exceed the typical anion–cation distance Rca by large. In this way, VDP analysis takes the nature of framework and migrating ions into account to transparently classify diffusion pathways according to their probability (cf. Figure 2, right). The necessary tools for VDP calculations are integrated into the well-maintained and powerful program ToposPro [19].

2.2 Mapping of scattering-length densities

Maps of the scatterer density are a valuable tool to (visually) identify diffusion pathways without the need to actually model them. At high temperatures, thermally activated migration of scatterers (i.e. the conducting ions) sets in. In the space- and time-averaged picture of diffraction, this movement manifests as a “smearing” of previously localized scatterer density along the diffusion paths. For neutron experiments, the scatterers are the atomic nuclei and their neutron-scattering lengths are the physical observables. Maps of the SLD are constructed from observed structure factors including their phases so that a reliable structure solution comprising the mobile ions is a prerequisite.

Classical maps acquired by Fourier synthesis often suffer from noise and – especially in the case of truncation in the high-angle regime – termination artifacts (see Figure 3, center). The reconstruction via maximum-entropy methods is a means to overcome this (see Figure 3, right). MEM is a versatile approach to the estimation of a model from a limited amount of information by maximizing information entropy under constraints consistent with observed physical quantities [21]. According to the principle of maximum entropy, the resulting probability distribution is “the least biased estimate possible on the given information” [22] and thus best represents the current state of knowledge. Being part of Bayesian statistics, the underlying theory stems from statistical mechanics in the 19th century and was further developed in the contexts of information theory and thermodynamics.

Fig. 3: SLD maps of ramsdellite-like Li2Ti3O7 (orthorhombic, Pnma) at 24°C. Located densities from calculated structure factors (left), Fourier transform of observed structure factors (center), and MEM-reconstructed densities (right). Dark blue: disordered lithium, gray: titanium, red: oxide ions (arbitrary radii); yellow and light blue: isosurfaces for +0.3 fm Å−3 and –0.3 fm Å−3, respectively. View approximately along [010]. Image generated using Vesta [20].

Fig. 3:

SLD maps of ramsdellite-like Li2Ti3O7 (orthorhombic, Pnma) at 24°C. Located densities from calculated structure factors (left), Fourier transform of observed structure factors (center), and MEM-reconstructed densities (right). Dark blue: disordered lithium, gray: titanium, red: oxide ions (arbitrary radii); yellow and light blue: isosurfaces for +0.3 fm Å−3 and –0.3 fm Å−3, respectively. View approximately along [010]. Image generated using Vesta [20].

In the crystallographic case, MEM yields the maximum variance of calculated structure factors within standard deviations of observed structure factors, thus producing feature-deprived scatterer-density maps that are less prone to misinterpretation [23], [24]. For their construction, the unit cell is first divided into N voxels, in which the scatterer density is regarded as uniform. With the normalized (to a sum of unity) scatterer density ρ′, the dimensionless information (or Shannon) entropy H is then formulated. This quantity differs from the thermodynamic (or Gibbs) entropy S only by a constant divisor of kB, the Boltzmann constant. If a normalized prior scatterer density τ′ is to be inferred, the relative entropy (or Kullback–Leibler divergence [25]) D is used instead:

H = S k B = i = 1 N ρ i ln ρ i

D ( ρ | | τ ) = i = 1 N ρ i ln ( ρ i τ i )

Using the method of Lagrange multipliers, H or D is then maximized with the sum over all normalized scatterer densities constrained to unity. Further constraints (mostly featuring a criterion of convergence instead of no tolerance) enter the maximization at this point. Amongst these, the so-called F constraint is the most important, as it introduces the observed structure factors and an uncertainty-based (and possibly lattice-plane dependent) weighting scheme. When dealing with powder profile data, a G constraint handling reflection overlap may be used in addition. Both of these depend on the calculated structure factors, which are the Fourier transforms of the variated densities ρ′. The hence necessary iterative process is implemented in the program Dysnomia [26], successor to Prima [27]. Classically, the data are treated with separate F and G constraints, which are quadratic with respect to the structure-factor or intensity differences (i.e. a least-squares condition). Dysnomia, in contrast, uses F and G constraints of multiple orders n=2, 4, 6, …, 16 and invariably merges those of the same order into generalized constraints Cn. These are then linearly combined using relative weights λn, which have to be adjusted by the user to warrant a Gaussian distribution of the final structure-factor and intensity differences. The challenges employing this method are

  • to assess the validity and adjust the standard deviations provided with the input data (by putting in a correction factor E);

  • to test and adjust the relative weights λn, so that the central moments of the distributions approach unity;

  • to check the physical reasonability of the final scatter density.

If these are overcome, the resulting map is in best possible agreement with the measured data and contains a minimum of artifacts (the smoother appearance has led to the term “MEM filtering”; cf. Figure 3, right). Its interpretation may, however, be tedious, as chunks of SLD can often not be assigned to one and only one species. Fortunately, lithium is one of the rare elements with negative scattering lengths in its natural mixture due to the large proportion of lithium-7. This increases the contrast to the SLD associated with the framework, being mostly positive. [1]

2.3 Probability-density functions and effective one-particle potentials

For well-behaved high-quality datasets, mapping of the SLD is not the end of the line. Instead, it is possible to model the “smearing” of scatterer density. This ultimately leads not only to the visualization of diffusion pathways but also to the calculation of effective migration barriers. To this end, the expression of the Debye–Waller factor (DWF) T ( q ) , describing intensity changes due to atomic displacement, has to be extended. In its common harmonic (i.e. Gaussian) formulation, it is fit to model thermal displacement of scatterers by vibration ( q : scattering vector, U: tensor of atomic mean-square displacement):

T harm ( q ) = e x p ( 2 π 2 q T U q )

Referring to a Cartesian basis, this leads to the representation of atoms as ellipsoids known so well from standard crystal-structure analysis. Now, if the DWF is to describe more complex SLD distributions caused by static (e.g. from incomplete ordering during cooling) or dynamic disorder (e.g. from thermal diffusion), it has to be extended with anharmonic terms. For this purpose, the Gram–Charlier or quasi-moment expansion is the most often used approach (qi: components of the scattering vector) [28]:

T anharm ( q ) = T h a r m ( q ) [ 1 + ( 2 π i ) 3 3 ! C j k l q j q k q l + ( 2 π i ) 4 4 ! D j k l m q j q k q l q m + ]

In the general case (site symmetry: 1), it introduces ten cubic, 15 quartic, dimensionless anharmonic tensorial coefficients Cjkl, Djklm, … per atom, which enter into the refinement as parameters. By far not all of them are significant so that they require careful evaluation. This also makes clear that high-quality data are necessary, especially in the high scattering-angle regime, where the effects of the DWF are most prominent.

The same way scatterer density and structure factors constitute a Fourier transform pair, an atom’s DWF and its probability-density function (PDF) [2] p ( u ) do. The latter describes the differential probability of finding the atom displaced by a vector u from its reference position ( [ p ( u ) ] = 1  Å 3 ) . It can be analytically calculated if the DWF takes the form of a Gram–Charlier expansion [ H ( u ) : Hermite polynomials] [29]:

p anharm ( u ) = p harm ( u ) [ 1 + 1 3 ! C j k l H j k l ( u ) + 1 4 ! D j k l m H j k l m ( u ) + ]

As PDF maps are derived from a model, it is possible to draw them for an ensemble of chosen atoms (e.g. all lithium ions in a unit cell). The resulting joint PDF is a very useful tool to visualize diffusion pathways and find the bottlenecks of migration (see Figure 4, left). Facilities to calculate DWFs including anharmonic terms and to derive PDFs from them are, e.g. provided by the program Jana2006 [31].

Fig. 4: Optimal diffusion path between two adjacent lithium positions in γ-LiAlO2 and position of the OPP maximum at 1475 K. Left: PDF isosurface for lithium probability density of 0.0125 Å−3 (pink) with section through Li1, Li1′, and Max [gray; ca. (0 10 47)]. Right: 3D OPP map for the same section (e1, e2: orthogonal coordinates with origin at Li1; dark gray: areas of infinitely high potential). Reprinted with permission from [30]. © 2016 American Chemical Society.

Fig. 4:

Optimal diffusion path between two adjacent lithium positions in γ-LiAlO2 and position of the OPP maximum at 1475 K. Left: PDF isosurface for lithium probability density of 0.0125 Å−3 (pink) with section through Li1, Li1′, and Max [gray; ca. (0 10 47)]. Right: 3D OPP map for the same section (e1, e2: orthogonal coordinates with origin at Li1; dark gray: areas of infinitely high potential). Reprinted with permission from [30]. © 2016 American Chemical Society.

A further step enables the evaluation of activation barriers for migration along certain paths. From the joint PDF, the effective OPP V ( u ) , as experienced by ions on average, can be calculated straightforwardly (T: temperature, w: occupation):

V ( u ) = k B T ln [ i w i p i ( u ) max u ( i w i p i ( u ) ) ]

In this approach, every ion is treated as an individual Einstein oscillator subject to the Boltzmann distribution (true at the classical high-temperature limit). The probe for the OPP is the migrating ions themselves. If they are to scan the potential surface, a sufficient number of ions must be able to populate regions near the potential barriers along the pathway. If this is not the case, barriers will inadvertently be overestimated, as hardly any ion even had enough kinetic energy to overcome them. Cum grano salis, the maximum of the OPP along a certain path is the associated migration activation energy Em, which does not comprise the defect formation energy Ef [32].

Jana2006 provides the opportunity to calculate the OPP along a polygonal chain with up to ten defined coplanar points. For two-dimensional data, we have developed the open-source program CalcOPP [33], which is written in FORTRAN 2003 and available as binary for Windows and Linux systems under an MIT license. With its help, it is possible to calculate OPPs from the PDFs within planes defined in Jana2006 and put them out in a format that is easily read into common plotting software for visualization (cf. Figure 4, right).

Very recently, efforts have been made to calculate OPPs directly from the MEM-reconstructed SLD [34], [35]. The appeal of this method lies in bypassing the cumbersome and potentially error-prone modeling of the PDF. However, it requires the assignment of parts of the SLD to the mobile species under investigation. For materials, in which migrating lithium ions are the only scatterers with a negative scattering length, this is trivial: all negative SLD is attributed to them. In other cases, one would have to discriminate SLD regions dominated by the mobile species from those that are not (during PDF modeling, this happens automatically). Unfortunately, software tools to conduct such studies are not yet available to the public.

2.4 cNEB computation

The cNEB calculations were performed with the plane-wave program package Vasp [36]. The atomic forces entering the cNEB procedure are calculated analytically at DFT level. In order to fulfill the Hellmann–Feynman theorem, which is the basis of the gradient calculation, the electronic wave functions have to be variationally optimized. In Vasp, the Bloch functions are expanded as Fourier series of plane waves. The expansion order is determined by the largest kinetic energy of a plane wave, the energy cutoff Ecut. This parameter therefore uniquely determines the variational quality of the wave function. Another parameter that determines the number of variational degrees of freedom in the wave function is the number of special k points in the irreducible part of the Brillouin zone in reciprocal space. It is determined by the procedure of Monkhorst and Pack [37], which employs the shrinking factors si as parameters. The reciprocal lattice vectors are subdivided into si intersections. The total number of k points is, therefore, s1×s2×s3. Both parameters, Ecut and si, were converged in preliminary studies for each system. In most plane-wave methods, the core electrons are replaced by effective core potentials. Vasp employs the projector-augmented wave (PAW) method [38] to project the local pseudo wave-functions from the plane-wave Bloch functions. This approach has been shown to provide higher numerical stability and improved accuracy for calculations of structural properties compared to norm-conserving and ultra-soft pseudopotentials. For all our calculations, the solid-state standard density-functional PBE [39], [40] has been used.

Figure 5 shows the principle of a NEB calculation. It starts with linear interpolation of the initial and final structure. This is simply achieved using the Cartesian or fractional coordinates of all atoms in the selected unit cell. If the linear interpolation leads to unphysically short interatomic distances or if a nonlinear migration pathway is expected, it is possible to introduce an intermediate structure (as was done in the example shown in Figure 6). During pathway optimization, modified atomic forces are used. The negative energy gradient V i of each image is decomposed into the tangential force F i and the radial contribution V i with respect to the arc connecting to the previous and the next image:

Fig. 5: Schematic contour plot of a potential-energy hypersurface (left; black points: initial and final structure) and construction of modified atomic forces (right; white points: interpolating images, black arrow: negative energy gradient, red arrow: tangential force, blue arrow: radial contribution, gray: connection to previous [i–1] and next [i+1] image).

Fig. 5:

Schematic contour plot of a potential-energy hypersurface (left; black points: initial and final structure) and construction of modified atomic forces (right; white points: interpolating images, black arrow: negative energy gradient, red arrow: tangential force, blue arrow: radial contribution, gray: connection to previous [i–1] and next [i+1] image).

Fig. 6: Example of a lithium migration path in 3R-LixTiS2. In the initial structure (left), one octahedral lithium site is vacant (white sphere). In the intermediate state (center), the lithium ion (red sphere) occupies a tetrahedral interstitial site. In the final structure (right), the ion has migrated to the formerly vacant octahedral site. Adapted with permission from [41]. © 2015 American Chemical Society.

Fig. 6:

Example of a lithium migration path in 3R-LixTiS2. In the initial structure (left), one octahedral lithium site is vacant (white sphere). In the intermediate state (center), the lithium ion (red sphere) occupies a tetrahedral interstitial site. In the final structure (right), the ion has migrated to the formerly vacant octahedral site. Adapted with permission from [41]. © 2015 American Chemical Society.

V i = V i + F i

During convergence of the images to the MEP, only the radial component is applied.

The calculations of migration pathways are usually performed with supercell models to reduce artifacts from periodic boundary conditions, such as direct and indirect defect–defect interaction, restriction of lattice relaxation, and enforced concerted movement of ions. Depending on the system under investigation, all possible defect models have to be explicitly considered via supercells with atomic vacancies or displacements, e.g. to model Frenkel defects. A possible situation with an empty lithium site in 3R-LixTiS2 as well as an associated migration event is illustrated by Figure 6. The principle shape of the potential curve along the corresponding reaction coordinate is shown in Figure 7. Herein, the occupation of the tetrahedral site is assumed to lead to a local energy minimum. Since initial and final structures are equivalent in this case, the MEP is symmetric and contains two maxima, which are represented by the cNEB images (black points) with highest energy. The distances between the images of highest energy are no longer forced to be identical. The cNEB procedure removes the spring forces between the image with the highest energy and its two neighbor images. In this way, the geometry of the image with highest energy is relaxed with fewer restrictions and can come closer to the transition structure.

Fig. 7: Schematic drawing of the minimum-energy path connecting the initial and final structures of Fig. 6. Three images are assumed between initial and intermediate structure for the cNEB calculation. The energy difference between the second image (signifying the TS) and the initial structure is an approximation to the migration activation energy Em. The MEP connecting intermediate and final state is assumed to mirror the first part.

Fig. 7:

Schematic drawing of the minimum-energy path connecting the initial and final structures of Fig. 6. Three images are assumed between initial and intermediate structure for the cNEB calculation. The energy difference between the second image (signifying the TS) and the initial structure is an approximation to the migration activation energy Em. The MEP connecting intermediate and final state is assumed to mirror the first part.

In principle, the cNEB procedure should find the MEP between any two given structures. However, since it is a gradient-based method, it is possible that the “MEP” obtained with NEB or cNEB corresponds to a local energy minimum or even a path over a higher-order saddle point. For this reason, it is necessary to perform full force-constant–matrix calculations for the cNEB maximum-energy images in order to verify their true TS nature. The lattice-vibration calculation can then also be used to calculate activation enthalpies and free activation energies by employing statistical thermodynamics. Also, the frequency prefactor is calculated in this way. Altogether, this allows for a (static) estimate of lithium-ion diffusivity.

3 Examples of application

3.1 Ultra-slow and 3D: γ-LiAlO2

We successfully applied the methods described above to different crystalline ternary lithium chalcogenides set apart by their conduction characteristics (fast/slow) and the connectivity of migration paths (1D/2D/3D). γ-LiAlO2 (tetragonal, P41212/P43212; see Figure 8) is an example of an ultra-slow three-dimensional lithium conductor – densely packed and lacking structural cation defects. It was investigated in different forms with a variety of techniques, such as NMR and impedance/conductivity spectroscopy on micro- to nanocrystalline [42] or amorphous powders [43], as well as on differently oriented single crystals [44], [45], [46]. Even dynamic mechanical analysis (DMA) has been performed [47].

Fig. 8: Crystal structure of γ-LiAlO2 with view approximately along [100] (left) and [001] (right); ions with arbitrary radii, unit cell in black. Reprinted with permission from [30]. © 2016 American Chemical Society.

Fig. 8:

Crystal structure of γ-LiAlO2 with view approximately along [100] (left) and [001] (right); ions with arbitrary radii, unit cell in black. Reprinted with permission from [30]. © 2016 American Chemical Society.

Our own work on this system started with a computational study of interstitial lithium-diffusion pathways [48]. The cNEB method implemented in the Vasp program was applied. The PEHS was calculated at PBE level using a plane-wave basis set defined by an energy cutoff of 400 eV. The supercell models contained one empty lithium lattice site, either by simply removing a litihium atom or by generating a Lii+VLi Frenkel defect pair. The calculated migration pathways correspond to a lithium ion moving from a regular lattice site to the empty lattice site (VLi) via an interstitial position. The calculated migration barriers for a jump process to the nearest empty lithium site were 0.65 eV (assuming a lithium point-defect) and 0.64–0.93 eV (assuming a Frenkel defect with decreasing LiiVLi distance). Migration to second- and third-nearest VLi sites have much larger barriers and are therefore less likely.

Our experimental investigations began with the evaluation of two single-crystal neutron diffraction datasets previously acquired at 25 and 770°C [49]. Although the high-temperature data were of low quality, maps of the MEM-reconstructed SLD hinted at strongly curved lithium pathways connecting adjacent positions. In an ensuing multi-method study [30], we assessed the state of research on lithium diffusion in different crystalline forms of γ-LiAlO2 and added tracer-diffusion experiments, conductivity spectroscopy, and high-temperature powder neutron diffraction. PVA and VDP (see Figure 2) identified the aforementioned pathway as the most probable one; PDF analysis confirmed this assumption (see Figure 4). Evaluation of the OPP revealed a migration barrier as high as 0.72(5) eV. [3] The difference between this and activation energies determined by other methods (1.0−1.2 eV) is too small for the formation of, e.g. Frenkel defects from an ideal domain, meaning that defects (twin boundaries, impurity dopants, etc.) must be abundant.

3.2 Fast and 2D: LixTiS2

In dependence on temperature, lithium content, and synthesis conditions, lithium titanium disulfides exist in two layered polytypes – differing in the anion stacking sequence (1T and 3R, see Figure 9) – and a cubic modification (c-LixTiS2). With x<1, the former are fast two-dimensional ion conductors within their structurally defective lithium layers. Ion diffusion in 1T-LixTiS2 (trigonal, Pm1), a paradigmatic battery material, was extensively investigated using conductivity and, in particular, NMR spectroscopy [50], [51], [52], [53], [54] supported by computation [55], [56].

Fig. 9: Crystal structures of 1T- (left) and 3R-LiTiS2 (right); ions with arbitrary radii, unit cell in black. Reprinted with permission from [41]. © 2015 American Chemical Society.

Fig. 9:

Crystal structures of 1T- (left) and 3R-LiTiS2 (right); ions with arbitrary radii, unit cell in black. Reprinted with permission from [41]. © 2015 American Chemical Society.

We succeeded in the reliable synthesis and structural characterization of fully lithiated (albeit not pure) 3R-LiTiS2 (trigonal, Rm), which computation showed to be metastable with respect to the 1T polytype by 4 kJ mol−1 at 0 K [57]. Pure phases with lower lithium content (x=0.7, 0.9) were shown to reversibly transform between 1T and 3R in the range of 400–600°C [58]. This transformation, also recently investigated via high-temperature NMR spectroscopy [59], is governed by a complex interplay of polarization and electron correlation that stabilize the 1T polytype at low temperature with respect to the electrostatically preferred 3R form.

Because our powder neutron-diffraction study on Li0.7TiS2 and Li0.9TiS2 required high temperature to thermally activate lithium-ion diffusion, it dealt with the 3R phases exclusively [41]. We used PDF/OPP analysis to test three probable paths: two along different directions strictly in plane and one including connecting tetrahedral voids above/below the lithium-containing layer (see Figure 6). A path only slightly bent out of the lithium plane showed to be preferred with migration barriers of 0.5 and 1.5 eV for x=0.7 and 0.9, respectively (see Figure 10; MEM-reconstructed SLD maps verified the soundness of the PDF/OPP analysis). Computation corroborated the preference found, but the calculated energy barrier was significantly lower than the experimental one (0.46 eV at 0 K) at higher lithium content. We attributed this discrepancy to the different nature of the methods (space-time average vs. individual-ion perspective).

Fig. 10: Section of the lithium OPP for Li0.7TiS2 (left) and Li0.9TiS2 (right) at 600°C in view of (010). Absolute Cartesian coordinates e1 along line between next neighboring lithium positions and e2 along c axis. Red, gray, and white areas: low, high, and infinitely high OPP, respectively; green line: preferred diffusion pathway. Adapted with permission from [41]. © 2015 American Chemical Society.

Fig. 10:

Section of the lithium OPP for Li0.7TiS2 (left) and Li0.9TiS2 (right) at 600°C in view of (010). Absolute Cartesian coordinates e1 along line between next neighboring lithium positions and e2 along c axis. Red, gray, and white areas: low, high, and infinitely high OPP, respectively; green line: preferred diffusion pathway. Adapted with permission from [41]. © 2015 American Chemical Society.

It should be noted that the erstwhile assumption of a pathway preference as a function of the lithium content was based on preliminary results and proved untenable [60]. For our research, 3R-LixTiS2 was a case of serendipity: similar materials in scope yielded unevaluable diffractograms (2H-LixNbS2) – caused by severe preferred orientation and concomitant selective line broadening due to stacking faults – or proved to be unstable (1T-LixCrS2).

3.3 Very fast and quasi-1D: ramsdellite-like Li2Ti3O7

The structure of ramsdellite-like Li2Ti3O7 (or, with Z=1, Li2.286Ti3.429O8; orthorhombic, Pnma) can be described as a stuffed derivative of the VO2 type or a defect-homeotype of Li0.5TiO2. As such, it consists of a framework mainly made up of TiO6 octahedra and spacious channels along b (in Pnma standard setting), in which interstitial disordered lithium ions reside (see Figure 11). Because of the presence of numerous voids, very fast conduction along b occurs. However, along the other crystallographic directions, ion movement is still considerable, making ramsdellite-like Li2Ti3O7 a quasi–one-dimensional (or strongly anisotropic) conductor. This has been shown and rationalized by a variety of diffraction, conductivity spectroscopy [62], and NMR studies [63]. Nonetheless, the lithium distribution between interstitial channels and framework has been a matter of debate for a long time.

Fig. 11: Crystal structure of ramsdellite-like Li2Ti3O7 in general view (left) and channel detail (right). View approximately along b, ions with arbitrary radii, unit cell in black. Reprinted from [61], © 2016, with permission from Elsevier.

Fig. 11:

Crystal structure of ramsdellite-like Li2Ti3O7 in general view (left) and channel detail (right). View approximately along b, ions with arbitrary radii, unit cell in black. Reprinted from [61], © 2016, with permission from Elsevier.

We acquired neutron powder diffractograms of the material at high temperature, which confirmed the recently favored formulation as [Li25]i[(Ti30.5)O7]f (□: vacancy, i: interstitial, f: framework) [61]. Because of the severe disorder of the interstitial lithium ions and the onset of thermally activated decomposition into the stable phases Li4Ti5O12 and rutile, we were unable to conduct a rigorous PDF/OPP analysis. Investigation of the MEM-reconstructed SLD (cf. Figure 3), however, showed the channels with their energetically very similar cation positions to be the major pathway for lithium diffusion. PVA and VDP (see Figure 1) identified an additional minor path through framework cation defects as accountable for the less efficient conduction along the other crystallographic directions.

Our computational study on ramsdellite-like Li2Ti3O7 considered a large number of possible lithium/titanium distributions over interstitial and framework sites in a Li16Ti24O56 supercell in order to mimic the complex structure [64]. The two most stable models were selected for cNEB calculations of lithium migration. Various pathways to empty nearest and next-nearest neighbor sites both along the channels and in directions perpendicular to the channels were considered. The calculated migration barriers showed large variations with respect to the lithium hopping-distance. The ranges, 0.20–0.92 eV along the channels and 0.30–0.85 eV in perpendicular directions, would indicate three-dimensional diffusion in contrast to the experimental results discussed above. However, the energy variation of initial and final structures of a migration process is much smaller along the channels than in the other directions. We therefore concluded that, although short-range lithium migration is possible between the channels, the lithium ions are trapped in energetically favored sites in this case. This effect is less pronounced along the channels, which leads to quasi-1D diffusion.

3.4 Strongly disordered: LiBi3S5

With LiBi3S5 (monoclinic, C2/m), we were able to selectively synthesize and characterize an announced, but hitherto elusive compound [65]. The material had been deemed a potentially fast one-dimensional lithium-ion conductor; to our dismay, NMR studies showed it to be moderate at best. In this case, we used neutron powder diffraction as a tool for structure elucidation: X-ray diffractograms had already revealed the structure to be a bismuth-disordered variant of the monoclinic AgBi3S5 type, but only neutron radiation could show the concomitant disorder in the lithium sublattice. Employing PVA and VDP as standalone methods, we identified viable pathways along the channels in b direction including tetrahedral voids for idealized configurations (see Figure 12). This revealed the reason for the worse-than-expected diffusivity: disordered bismuth ions clog the pathways and effectively prohibit long-range transport. For an experimental substantiation, neutron diffraction at high temperatures would be the method of choice – if one could overcome the challenges posed by cation disorder and low symmetry (reflection overlap).

Fig. 12: Detail of ordered model of LiBi3S5 with migration path. Blue: bismuth, green: lithium, yellow: sulfide ions, red: elementary voids and channels; ions with arbitrary radii, unit cell in black. Reprinted from [65], © 2016, with permission from Elsevier.

Fig. 12:

Detail of ordered model of LiBi3S5 with migration path. Blue: bismuth, green: lithium, yellow: sulfide ions, red: elementary voids and channels; ions with arbitrary radii, unit cell in black. Reprinted from [65], © 2016, with permission from Elsevier.

4 Conclusion

Diffusion pathways and activation energies in crystalline lithium-ion conductors can be formidably investigated using neutron diffraction as experimental and cNEB calculations as computational methods. We have established a versatile set of evaluation approaches that can be adapted to data quality and necessary information output: topological analyses (PVA and VDP), examination of MEM-reconstructed SLD maps, and PDF/OPP modeling. Ideally, these are complemented by additional techniques, such as NMR relaxometry, conductivity/impedance spectroscopy, tracer-diffusion experiments, or MD computations to verify scientific soundness.

We used these methods systematically to characterize lithium diffusion in paradigmatic ternary lithium oxides and sulfides with different conduction properties (fast/slow) and dimensionalities (one-/two-/three-dimensional), like γ-LiAlO2, 1T- and 3R-LixTiS2, ramsdellite-like Li2Ti3O7, and LiBi3S5. We feel that, with our research, we have only just scratched the surface and that more substance classes and evaluation methods (e.g. the inspection of the bond-valence sum mismatch) are waiting for scientists to dig deeper.

Acknowledgements

Financial support by the Deutsche Forschungsgemeinschaft (FOR 1277: “Mobilität von Lithiumionen in Festkörpern [molife]”, TP 2, LE 781/16-2, BR 1768/5-2) is gratefully acknowledged.

References

1. G. Henkelman, H. Jónsson, J. Chem. Phys. 113 (2000) 9978. Search in Google Scholar

2. D. Sheppard, P. Xiao, W. Chemelewski, D. D. Johnson, G. Henkelman, J. Chem. Phys. 136 (2012) 074103. Search in Google Scholar

3. G. Henkelman, B. P. Uberuaga, H. Jónsson, J. Chem. Phys. 113 (2000) 9901. Search in Google Scholar

4. P. Xiao, D. Sheppard, J. Rogal, G. Henkelman, J. Chem. Phys. 140 (2014) 174104. Search in Google Scholar

5. D. Frenkel, B. Smit, Understanding Molecular Simulation – From Algorithms to Applications, Academic Press, San Diego, USA (2002). Search in Google Scholar

6. R. Car, M. Parrinello, Phys. Rev. Lett. 55 (1985) 2471. Search in Google Scholar

7. M. Iannuzzi, A. Laio, M. Parrinello, Phys. Rev. Lett. 90 (2003) 238302. Search in Google Scholar

8. M. Avdeev, V. B. Nalbandyan, I. L. Shukaev, in V. V. Kharton (Ed.), Solid State Electrochemistry I: Fundamentals, Materials and their Applications, Wiley-VCH, Weinheim, Germany (2009), p. 227. Search in Google Scholar

9. F. L. Hirshfeld, Theor. Chim. Acta 44 (1977) 129 Search in Google Scholar

10. M. A. Spackman, P. G. Byrom, Chem. Phys. Lett. 267 (1997) 215. Search in Google Scholar

11. M. J. Turner, J. J. McKinnon, D. Jayatilaka, M. A. Spackman, CrystEngComm 13 (2011) 1804. Search in Google Scholar

12. M. Ø. Filsø, M. J. Turner, G. V. Gibbs, S. Adams, M. A. Spackman, B. B. Iversen, Chem. Eur. J. 19 (2013) 15535. Search in Google Scholar

13. D. B. Newell, B. N. Taylor, P. J. Mohr, CODATA Recommended Values of the Fundamental Physical Constants: 2014 (2015), available from world wide web: https://dx.doi.org/10.5281/zenodo.22826 [cited July 29, 2016]. Search in Google Scholar

14. S. K. Wolff, D. J. Grimwood, J. J. McKinnon, M. J. Turner, D. Jayatilaka, M. A. Spackman, CrystalExplorer – Crystal Structure Analysis with Hirshfeld Surfaces, v. 3.1 (2012), available from world wide web: http://www.hirshfeldsurface.net. Search in Google Scholar

15. P. Niggli, Z. Kristallogr. 65 (1927) 391. Search in Google Scholar

16. A. S. Mitchell, M. A. Spackman, J. Comput. Chem. 21 (2000) 933. Search in Google Scholar

17. V. A. Blatov, G. D. Ilyushin, O. A. Blatova, N. A. Anurova, A. K. Ivanov-Schits, L. N. Dem’yanets, Acta Crystallogr., Sect. B Struct. Sci. 62 (2006) 1010. Search in Google Scholar

18. V. A. Blatov, Crystallogr. Rev. 10 (2004) 249. Search in Google Scholar

19. V. A. Blatov, A. P. Shevchenko, D. M. Proserpio, Cryst. Growth Des. 14 (2014) 3576. Search in Google Scholar

20. K. Momma, F. Izumi, J. Appl. Crystallogr. 44 (2011) 1272. Search in Google Scholar

21. F. Izumi, Solid State Ionics 172 (2004) 1. Search in Google Scholar

22. E. T. Jaynes, Phys. Rev. 106 (1957) 620. Search in Google Scholar

23. A. Senyshyn, H. Boysen, R. Niewa, J. Banys, M. Kinka, Y. Burak, V. Adamiv, F. Izumi, I. Chumak, H. Fuess, J. Phys. D Appl. Phys. 45 (2012) 175305. Search in Google Scholar

24. M. Yashima, J. Ceram. Soc. Jpn. 117 (2009) 1055. Search in Google Scholar

25. S. Kullback, R. A. Leibler, Ann. Math. Stat. 22 (1951) 79. Search in Google Scholar

26. K. Momma, T. Ikeda, A. A. Belik, F. Izumi, Powder Diffr. 28 (2013) 184. Search in Google Scholar

27. F. Izumi, R. A. Dilanian, in S. G. Pandalai (Ed.), Recent Research Developments in Physics Volume 3, Transworld Research Network, Trivandrum, India (2002), p. 699. Search in Google Scholar

28. K. N. Trueblood, H.-B. Bürgi, H. Burzlaff, J. D. Dunitz, C. M. Gramaccioli, H. H. Schulz, U. Shmueli, S. C. Abrahams, Acta Crystallogr. Sect. A Found. Crystallogr. 52 (1996) 770. Search in Google Scholar

29. W. Kuhs, Acta Crystallogr. Sect. A Found. Crystallogr. 48 (1992) 80. Search in Google Scholar

30. D. Wiedemann, S. Nakhal, J. Rahn, E. Witt, M. M. Islam, S. Zander, P. Heitjans, H. Schmidt, T. Bredow, M. Wilkening, M. Lerch, Chem. Mater. 28 (2016) 915. Search in Google Scholar

31. V. Petříček, M. Dušek, L. Palatinus, Z. Kristallogr. – Cryst. Mater. 229 (2014) 345. Search in Google Scholar

32. H. Boysen, Z. Kristallogr. 218 (2003) 123. Search in Google Scholar

33. D. Wiedemann, CalcOPP – Calculation of 2D OPP from PDF Data, v. 1.6.1 (2015), available from world wide web: https://dx.doi.org/10.6084/m9.figshare.1461721. Search in Google Scholar

34. M. Monchak, T. Hupfer, A. Senyshyn, H. Boysen, D. Chernyshov, T. Hansen, K. G. Schell, E. C. Bucharsky, M. J. Hoffmann, H. Ehrenberg, Inorg. Chem. 55 (2016) 2941. Search in Google Scholar

35. D. A. Weber, A. Senyshyn, K. S. Weldert, S. Wenzel, W. Zhang, R. Kaiser, S. Berendts, J. Janek, W. G. Zeier, Chem. Mater. 28 (2016) 5905. Search in Google Scholar

36. G. Kresse, J. Furthmüller, Phys. Rev. B Condens. Matter Mater. Phys. 54 (1996) 11169. Search in Google Scholar

37. H. J. Monkhorst, J. D. Pack, Phys. Rev. B Solid State 13 (1976) 5188. Search in Google Scholar

38. G. Kresse, D. Joubert, Phys. Rev. B Condens. Matter Mater. Phys. 59 (1999) 1758. Search in Google Scholar

39. J. P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 78 (1997) 1396. Search in Google Scholar

40. J. P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865. Search in Google Scholar

41. D. Wiedemann, M. M. Islam, S. Nakhal, A. Senyshyn, T. Bredow, M. Lerch, J. Phys. Chem. C 119 (2015) 11370. Search in Google Scholar

42. E. Witt, S. Nakhal, C. V. Chandran, M. Lerch, P. Heitjans, Z. Phys. Chem. 229 (2015) 1327. Search in Google Scholar

43. J. Rahn, E. Witt, P. Heitjans, H. Schmidt, Z. Phys. Chem. 229 (2015) 1341. Search in Google Scholar

44. S. Indris, R. Uecker, P. Heitjans, Diffus. Fundam. 2 (2005) 50.1. Search in Google Scholar

45. S. Indris, P. Heitjans, R. Uecker, B. Roling, J. Phys. Chem. C 116 (2012) 14243. Search in Google Scholar

46. S. Indris, P. Heitjans, R. Uecker, T. Bredow, Phys. Rev. B Condens. Matter Mater. Phys. 74 (2006) 245120. Search in Google Scholar

47. J. Langer, D. Wohlmuth, A. Kovalcik, V. Epp, F. Stelzer, M. Wilkening, Ann. Phys. (Berlin, Ger.) 527 (2015) 523. Search in Google Scholar

48. M. M. Islam, T. Bredow, J. Phys. Chem. Lett. 6 (2015) 4622. Search in Google Scholar

49. D. Wiedemann, S. Indris, M. Meven, B. Pedersen, H. Boysen, R. Uecker, P. Heitjans, M. Lerch, Z. Kristallogr. Cryst. Mater. 231 (2016) 189. Search in Google Scholar

50. P. Heitjans, A. Schirmer, S. Indris, in P. Heitjans, J. Kärger (Eds.), Diffusion in Condensed Matter, Springer, Berlin, Germany (2005), p. 367. Search in Google Scholar

51. M. Wilkening, P. Heitjans, Defect Diffus. Forum 237–240 (2005) 1182. Search in Google Scholar

52. M. Wilkening, W. Küchler, P. Heitjans, Phys. Rev. Lett. 97 (2006) 065901. Search in Google Scholar

53. M. Wilkening, P. Heitjans, Diffus. Fundam. 6 (2007) 40.1. Search in Google Scholar

54. M. Wilkening, P. Heitjans, Phys. Rev. B Condens. Matter Mater. Phys. 77 (2008) 024311. Search in Google Scholar

55. T. Bredow, P. Heitjans, M. Wilkening, Phys. Rev. B Condens. Matter Mater. Phys. 70 (2004) 115111. Search in Google Scholar

56. M. M. Islam, T. Bredow, Z. Phys. Chem. 226 (2012) 449. Search in Google Scholar

57. S. Nakhal, M. Lerch, J. Koopman, M. M. Islam, T. Bredow, Z. Anorg. Allg. Chem. 639 (2013) 2822. Search in Google Scholar

58. D. Wiedemann, S. Nakhal, A. Senyshyn, T. Bredow, M. Lerch, Z. Phys. Chem. 229 (2015) 1275. Search in Google Scholar

59. M. Wilkening, Philos. Mag. 95 (2015) 861. Search in Google Scholar

60. D. Wiedemann, S. Nakhal, A. Senyshyn, M. Lerch, Z. Anorg. Allg. Chem. 640 (2014) 2342. Search in Google Scholar

61. D. Wiedemann, S. Nakhal, A. Franz, M. Lerch, Solid State Ionics 293 (2016) 37. Search in Google Scholar

62. J. B. Boyce, J. C. Mikkelsen Jr., Solid State Commun. 31 (1979) 741. Search in Google Scholar

63. J. Heine, M. Wilkening, P. Heitjans, Diffus. Fundam. 11 (2009) 47.1. Search in Google Scholar

64. M. M. Islam, P. Heitjans, T. Bredow, J. Phys. Chem. C 120 (2016) 5. Search in Google Scholar

65. S. Nakhal, D. Wiedemann, B. Stanje, O. Dolotko, M. Wilkening, M. Lerch, J. Solid State Chem. 238 (2016) 60. Search in Google Scholar