Molecular Dynamics Simulations of Ionic Motions: Dynamic Heterogeneity as a Basis of Studies ofNanostructured Materials
Molecular Dynamics Simulations of Ionic Motions: Dynamic Heterogeneity as a Basis of Studies of Nanostructured Materials
In this chapter, the mechanism of ionic motions in silicate systems derived from MD simulations will be discussed. The time development of mean-squared displacement (MSD) of ions [1, 2] (corresponding to frequency dependence of ionic motion [1-7]) shows several characteristic time regions [8, 9]. Caged ion dynamics in early time (nearly constant loss region (NCL)) region take important roles in the enhancement of dynamics found in porous lithium disilicate in the NVE condition, and details of short- time dynamics will be discussed in [8, 9] and Chapters 6 and 7. Behaviors of characteristic times found in the MSD of ions can be a basis to understand how dynamics are affected by further modifications such as introducing pores. That is, they are useful to distinguish possible mechanisms.
Molecular Dynamics of Nanostructures and Nanoionics: Simulations in Complex Systems Junko Habasaki
Copyright © 2021 Jenny Stanford Publishing Pte. Ltd.
ISBN 978-981-4800-77-8 (Hardcover), 978-1-003-04490-1 (eBook) www.jennystanford.com
Characteristics of Ionic Motions in Ionically Conducting Melts and Glasses
Dynamics of ions are the complex ones having many questions to be solved and are continued to be discussed in papers [1-5, 10], reviews [11-14], and books [6, 7]. In recent works, several points have been addressed and made clearer, where molecular dynamics simulations play roles to gain a nanoscopic understand of these complex behaviors.
One of the challenging problems is the relationships between AC conductivity and DC conductivity, where many functional forms are assumed empirically. Another interesting problem in ionics is the mechanism of the "mixed alkali effect" [6, 7, 15-20], which is extensively examined by molecular dynamics simulations [21-25] as well as experiments and theoretical methods.
In this chapter, we focus on the former phenomena observed in single alkali silicates and related systems. As already mentioned, these results are important as a basis to consider the effect of modification of Li ions by introducing pores, as will be discussed in Chapters 6 and 7.
In experiments, it is well known that conductivity spectra show complex behaviors and let's start a discussion from this behavior.
An example of this behavior is shown in Fig. 5.1 [4, 5]. The ac conductivity o'(со) of ionic systems shows a constant value at low frequencies. The value is called dc conductivity. In higher frequency, several regions are observed. The higher frequency region than vxl showing up in the conductivity spectrum and/or susceptibility is called nearly constant loss (NCL) region. Details of the assignment of frequencies characterizing regions will be explained later.
From this kind of common observations, the real part of the conductivity is often expressed as follows [26-28]:
![]()
where odc is the dc conductivity, A is a temperature-dependent parameter and n is a fractional exponent which usually lies between 0.6 and 1. In this expression, ac and dc seem to be independent and overlap.
Figure 5.1 Log10 a'{у) versus log (frequency/Hz) plot of data at T = 293 К of 0.80LiF-0.20Al(P03)3 obtained by Kulkami et al. (Kulkami, A. R., Myer, R, Schnider, U., Lunkenheimer, R, Loidl, A. (1999), Ionic Conductivity and Relaxations in Li and (Li and K) Fluorophosphate Glasses, Ceramic Transactions, 92, p. 115) to show the nearly linear frequency dependence of the NCL by the line with slope equal to 0.99. The line with slope 0.44 is the high-frequency dependence of logio cr'(v) corresponding to the Kohlrausch fit to the electric modulus data (shown in the inset) with n = 0.44. The dashed line in the main figure is the fully cooperative ion conductivity relaxation. The crossover frequencies and vx2 as well as the independent hopping frequency Vo are determined from the data in the same manner as explained in Fig. 5.1 (in the following reference). Reprinted from Ngai, K. L., Leon, C. (2002), Cage decay, near constant loss, and crossover to cooperative ion motion in ionic conductors: Insight from experimental data, Phys. Rev. B, 66, 064308, pp. 1-11.
Combination of this expression with Kramers-Kronig relations leads to the frequency dependence of the complex conductivity, which can be written as
![]()
This expression is explained as the sum of high-frequency permittivity term, jcoe0£x, due to the contributions from all the polarization at higher frequencies, and the power-law frequency dependence (at w > a>c), where coc representing a crossover frequency from the dc conductivity plateau (at w < ojc).
Many experimental data exhibiting such behaviors can be found in books [6, 7]. Origin of such universal behavior of the frequency dependence of the conductivity has been a target of debates for a long time [28-31].
Mean-Squared Displacement and Characteristic Time Regions
In molecular dynamics simulations, mean-squared displacement (MSD) of ions or atoms can be obtained using the following relations:
where, г,- is the positional vector of the /-th particle.
To learn the origin of these behaviors in the frequency domain, it is useful to examine the time development of mean-squared displacement of ions.
MSD of ions as a function of time shows several distinctly different regions with different functional forms. Figure 5.2 shows the typical behaviors of ionic motion, that is, MSD of Li ions for lithium metasilicate at 700 К obtained by MD simulation. In this figure, characteristic times, txl and tx2, are defined for representing changes in the functional forms in MSD. The value r0 obtained based on the following relation of the coupling model [4, 5] is also shown for comparison.
![]()
This is a representative example of time (or frequency)- dependent dynamics, where the method can cover the wide time ranges from ~fs to several ~tens ns.
The complex conductivity o*(oj) is in the frequency domain can be connected to the MSD in the time domain by the following equation [32, 33]:

where Np is a number density of mobile ions, qc (= Zi e) is a charge of the ion, к is the Boltzmann constant, HR is the Haven ratio [34], and Г is temperature.
Figure 5.2 Mean square displacement of Li ions in lithium metasilicate at 700 K. Characteristic times to separate different time regions are shown. The value r0 of the Coupling model [7] is also shown. Reprinted with permission from Habasaki, J., Ngai, K. L., Hiwatari, Y. (2002), Molecular dynamics study of cage decay, near constant loss, and crossover to cooperative ion hopping in lithium metasilicate, Phys. Rev. E, 66, 021205, pp. 1-11. Copyright (2002) by the American Physical Society.
Based on this relation, characteristic times of MSD can explain why the different regions are found in conductivity spectra. The events in each region can be examined further in MD as discussed later.
The Haven ratio in Eq. (5.5) can be experimentally obtained from the ratio of tracer diffusion and conductivity. It can also be represented by the ratio
as shown in the linear-response theory by Kubo [28]. The Haven ratio is regarded as an indicator not only of the geometrical correlations but also of the cooperativity in the ionic motions because it includes the cross term of ions i and j.
Needless to say, this expression can be used to calculate the ratio in MD in principle. However, one should note that it is not easy to obtain the Haven ratio in good statistics from MD data especially at low temperatures, because of the cooperative character of the denominator and existence of rare events. The same is true for conductivity itself, although several attempts have been done so far.
Even a MSD for the self-motion tends to show large fluctuations due to the contribution of rare events with cooperative nature. If the calculation time is not long enough, rare events will not be assessed. From the long run, many time series can be taken for analyses, while it may not work well to observe rare events by many independent short runs, because of an extremely large number of runs needed for this purpose [6].
Figure 5.3 Characteristic times found in MSD of Li ions to separate different functional forms. This example was taken for lithium disilicate glass at 800 K. Dashed and dotted lines (fitted to power laws) for Li ions at each temperature are for NCL and power law regions, respectively.
In Fig. 5.3, a schematic description of the MSD with characteristic times of Li ions in glassy Li2Si205 systems is shown for the data at 800 K. It is separated by characteristic times, txl, tx2, and tdif. Comparable time regions as in lithium metasilicate are found in lithium disilicate, Li2Si205, in the silica-rich region, and therefore the same assignments are applicable. In this figure, a concise explanation of each region and the relation between time scale and length scales are also depicted. Here the distance Rdlf corresponds to the mean distance among Li ion sites.
In silica-rich systems, the Haven ratio is known to be only slightly smaller than 1, and therefore changes in diffusive motion are representative of those in the conductivity.
By using MD, an examination of a wide region can be done relatively inexpensively. If one would like to examine the relaxation phenomena experimentally, usually combinations of several expensive apparatuses will be required.
Events in Each Characteristic Time Region
Each time-region in these MSDs can be characterized as follows by using van Hove functions, fractal dimension analyses [20], intermediate scattering functions and cage decay functions [1, 2, 10] for lithium silicate systems with several compositions [6]. It is worth to mention that quite similar behaviors are found in many systems, including non-ionic systems [6, 7].
Four distinct time regions for MSD of Li ions are found and time developing of it corresponds to the processes of ionic motion, from the localized motion within a cage (I) to final diffusive regimes (IV) through the regions with first jump motion of some ions (II), and power law regions (III), where bifurcated motions of fast and slow ions are mixing.
Region I; NCL region: After the initial overshoot of the ballistic motion corresponding to the boson peak, a nearly flat region assigned to the NCL is found, where the MSD increases with time as t° with a * 0.07. In this time region, ions are still caged by surrounding oxygen atoms and further shells of ions and atoms.
This behavior corresponds to the frequency dependence of the real part of the conductivity o’(or)
Thus, the time scale of txl is related to the waiting time of jump motions.
Region II; Primitive jump region: Initial jumps of several ions are found in this region, where the MSD rises more rapidly than f° of the regime (I). In this intermediate time regime of about txl < t < tx2, MSD is nearly proportional to time and hence regarded as a quasi-diffusive regime. This statement can be easily confirmed by a linear scale plot.
It is interesting to note that the jump motion behaves as independent jumps here because the observation of squared displacement is usually done along with the direction of the first jump. After that, it is modified by angles among successive jumps. The jump can be a part of the successive jumps. In this sense, the difference between the first jump and the following jumps depends on when one starts the observation [19].
Region III; Power law region: Power law behaviors of MSD are found in the region tx2< t < tdif. In this region, MSD depends on t6 with 9 < 1 and the small value of 9 means the existence of slower dynamics compared with normal diffusion with the Gaussian type. Both fast and slow ions are found here and the slope 9 is affected by the back-correlated motion of jumps considerably in the glass. The exponent 9 is a function of the fractal dimension of the random walk [35], dw which is an index to represent the complexity of trajectories. For the mean behavior of ions, 9 is expressed by Eq. (5.7) [36].

Region IV: Diffusive region: Diffusion coefficient of ions can be calculated after tdif, where MSD reaches the square of the typical neighboring distance of ion sites and the MSD shows a linear t dependence in that region. Thus, the time scale of tdif determines the elementary process of diffusive motion modified by correlations among successive motions. One should note that the activation energy calculated from the diffusion coefficients is not for elemental jumps but the mean behaviors of both fast and slow ions. With decreasing temperature, each time- region becomes longer and longer. Namely, the diffusive regime is found after 2 ns at 800 К for Li ions, while it is found after ~15 ns at 600 K. The validity of these assignments can be confirmed by the functional forms and van Hove functions as well as observed values of
In Chapter 6, the characteristics of each time-region in lithium disilicate will be compared with the motion in the modified system by introducing pores.
In the following sections, details of ionic motions will be examined from caged ion dynamics to diffusive dynamics step by step.
Many examples are taken from works for both lithium disilicate and lithium metasilicate systems and related ones. Of course, there are some differences among systems, but the following discussions hold for these cases with some shifts in time and temperature regions.