(NCL) in Lithium Metasilicate Glass at Low Temperatures—Anisotropic

Article pubs.acs.org/JPCC

Nearly Constant Loss (NCL) in Lithium Metasilicate Glass at Low TemperaturesAnisotropic and Dynamical Caging from Molecular Dynamics Simulations Junko Habasaki*,† and K. L. Ngai‡ †

School of Materials and Chemical Technology, Tokyo Institute of Technology, Nagatsuta 4259, Yokohama 226-8502, Japan IPCF-CNR, Università di Pisa, Largo Bruno Pontecorvo 3, I-56127 Pisa, Italy

‡

ABSTRACT: Ions in ionically conducting materials show complex but universal trends in the time development of the mean squared displacements (MSD) of ions, which correspond to the complex frequency dependence conductivity spectra or susceptibility. Molecular dynamics simulations enable us to characterize the details of such motions including the initial nearly logarithmic or weak power law time dependence of the MSD of caged ions, showing up in susceptibility as the nearly constant loss (NCL) and the important connections of the caged ion dynamics to ion hopping dynamics in ionically conducting materials. In previous works, the motion of ions in the NCL region for lithium metasilicate was characterized at 700 and 500 K. Since each time region of MSD becomes longer and longer with decreasing temperature, the NCL region will cover the entire observation time during the MD run at lower temperatures, and thus enabling a more detailed study of the properties of the NCL. This paper reports results from new simulations at temperatures lower than ever done before down to 100 K of long time scale. Temperature dependence of the time regime of the MSD corresponding to the NCL and the cage decay that terminates it are obtained from the self-part of the Van Hove function. The results at 100 K clearly show the behavior expected as well as new features. These include anisotropy of the shapes of cages, and some characteristics of the ionic motion within fluctuating cages, which are visualized by the simulations. The potential field in the NCL region is explained by the multipole expansions. The dynamics of ions in the caged regime are characterized by anharmonic motions within cages, and from which the NCL originates. Moreover, the dynamics shows strong intermittency. The relation to similar behaviors in glass forming liquids is discussed.

■

INTRODUCTION

molecules (or basic structural units) in other nonionic materials. Recently, the importance of the caged ion dynamics as well as caging of atoms in glass forming metallic alloys is pointed out in several works on the relationship between nanostructures and dynamics, which are examined for both fundamental and application purposes.4,13−16 For example, enhancement of ion dynamics was predicted in molecular dynamics simulations under NVE conditions of nanoporous lithium disilicate,13 where the enhancement is found to be related to the structural changes of cages (coordination polyhedra), accompanying the shortening of the time scale of caging. Mean squared displacements (MSD) of ions as a function of time show several distinctly different regions roughly separated by characteristic times, tx1, tx2, and tdif as defined in a previous work.5 The MSD related to the NCL ends at the first characteristic time tx1, which corresponds to the onset of the jump motion of ions or particles to neighboring sites, and is followed by the start of another power law dependence of MSD

There are a number of questions on the complex dynamics and mechanism of ionics which remain to be answered and is called into attention in reviews and books.1−4 Recently several points have been addressed and made clearer.3,4 However, some questions will remain unresolved until there is consensus among the researchers of the answers given. Meanwhile, molecular dynamics simulations continue to serve as a research tool to gain microscopic understanding of these complex behaviors. In particular, it is the caged ion dynamics5,6 seen by simulations in terms of the complex time development of the mean square displacement (MSD) of the ions, or the corresponding frequency dependent conductivity and susceptibility spectra. At short times, the MSD is weakly time dependent and describable by either a logarithmic or weak power law function, and as a function of frequency corresponds in susceptibility to the nearly constant loss (NCL). The origin and physical mechanism of the short time MSD and the high frequency NCL is still a matter of debate.7−12 History of the NCL as well as results from recent works are summarized in ref 4. As shown therein, existence of the NCL is known to be rather general, and that is NCL of ionics also exists as NCL of © XXXX American Chemical Society

Received: April 3, 2017 Revised: May 27, 2017

A

DOI: 10.1021/acs.jpcc.7b03132 J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C

■

after tx2, and finally the long-term diffusive region starting after tdif.4,5,17,18 With decreasing temperature or increasing pressure, each time region is expected to be longer and longer due to thermal activation governing their characteristic times.6,7 Of course, the rigorous dynamics at lower temperatures is not a simple time-scaled version of the dynamics at higher temperatures,19 and new features can appear. If the NCL region becomes long enough to cover the entire region of observation time, then the MSD becomes flatter and the intermediate scattering function and/or other related correlation functions decay little. In this time region, essentially no ion can jump out of the cage, although in the case of ionically conducting glasses near room temperature, there is still a low probability for jumps. In previous MD simulations, we analyzed the motion of ions in lithium metasilicate glass at temperatures equal to or higher than 500 K5,6 In the present work, MSD corresponding to the NCL are obtained over a more extended time region at lower temperatures. To characterize the ion dynamics in the NCL regime under these conditions, temperature dependences of the characteristics times as well as diffusive motions are analyzed to give a detailed description of the dynamics of the ions while confined within the cages. The objective of the present paper is to bring out more microscopic properties of the dynamics of ions and its temperature dependence in the short time/high frequency region by performing molecular dynamics simulations at low temperatures.

Article

RESULTS AND DISCUSSIONS Mean squared displacements of ions are useful to characterize the time development of dynamics from caged region to diffusive motions.3−7 The typical behavior is summarized below. In Figure 1(a), MSDs of Li ions in molten and glassy

■

MOLECULAR DYNAMICS SIMULATIONS Simulations of xLi2O-(1-x)SiO2 with x = 0.5 in the molten and glassy states were performed as in previous works5,6,20−24 for the system with 3456 particles. The potential function, which consists of charge interactions, repulsive potential of the Gilibert-Ida type25,26 and an r−6 term, was used. Potential parameters based on ab initio molecular orbital calculations of SCF level developed by one of us24 are used. In MD simulations, Coulombic force was calculated by the Ewald method, and the usual periodic boundary (PBC) conditions were imposed. Validity of the parameters of the potential used to simulate the dynamics is important to obtain the characteristic times in the present work. So far, several authors have successfully compared dynamics found in experiments and MD simulations. For example, diffusion coefficients and viscosities obtained by Gonçalves and Rino27 using our parameters for lithium disilicate systems are found to be comparable to the experimental values shown in ref 28. It means that the time scales of dynamics derived in the present model are reasonable. Further details of comparison including the results of different cooling rate can be found in ref 4. Starting from the equilibrated structure at 1000 K obtained in our previous work, target temperatures from 100 to 900 K are attained by NPT ensemble using a cooling rate of ∼0.1 K/ ps using the Nose-Hoover methods. In this paper, our attention is mainly focused on the results in the glassy state within the range of 100−500 K. After (quasi) equilibration had been reached at the target temperature, successive runs were performed at each temperature during several ns to several tens of ns by NVE ensemble. In the NVE runs, system volume and axis lengths were adjusted to those obtained by NPT runs under atmospheric pressure. The role of different extended conditions such as pressure in the glass transition problem is argued in a separate paper29 for the case of a soft-spheres system.

Figure 1. (a) Mean squared displacement of Li ions in Li2SiO3 at 100, 300, 500, 600, 700, 800, and 900 K (from bottom to top). (b) Schematic description of several regions found in MSD (The data are taken at 800 K). Dashed lines are for showing power laws. The first power law region corresponds to NCL.

Li2SiO3 systems at 100, 300, 500, 600, 700, 800, and 900 K (from bottom to top) are shown. After the initial overshoot of the ballistic motion corresponding to the boson peak, the time dependence of MSD gradually changes. Ionic motion of Li ions found in MSD can be broken up to consist of several time regions as represented for the data at 800 K in Figure 1(b). Regions are separated by the characteristic times, tx1, tx2, and tdif (defined in ref 5), as indicated. Region I is found in an early time regime between approximately 0.2 and 2 ps at 800 K, where the MSD has a power law dependence, ⟨r2⟩ ∝ ctα, with small positive exponent α ≈ 0.2. It corresponds to the frequency dependence of the real part of the conductivity σ ′(ω) ∝ ω1‑α and the imaginary part of dielectric susceptibility, ε″(ω) ∝ ω−α. Because of the small α value, region I corresponds to a nearly constant loss (NCL). In this region, practically all ions are still confined within the cages formed by B

DOI: 10.1021/acs.jpcc.7b03132 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C matrix oxygen atoms and further away shells. In region II, the MSD is nearly proportional to time, where the initial jumps to the neighboring sites are observed for several ions. This linearity may not be clear in a log−log scale plot but it is easy to see it in a linear scale plot as shown in a figure in ref 30. Region II ends at tx2. In region III, MSD shows a power law with the exponent less than 1. After the end of region III at tdif, a diffusive region is attained. As shown in Figure 1(a), at 300 and 100 K, MSDs are nearly flat after the overshoot at these lower temperatures. This region corresponds to the NCL, where almost all ions are located within the initial cages over a longer time region. The caging of ions by the surrounding atoms will be proven later by the selfpart of the Van Hove functions. Naturally, the contribution of jump motions becomes smaller with decreasing temperature, and almost all ions are located within some cage at low temperatures at the initial times of observations. Nevertheless, we note that the initial time can be arbitrarily chosen, and it is possible to start the observation from the beginning of the jump or at the saddle point. This probability is low but contribution to MSD is non-negligible at 300 K. If one used a short observation time for analysis and a short time run, then such contribution of rare events starting at early times might be completely missing. In this case, MSD in the short time is systematically smaller than that obtained for long time scale MSD, namely a gap appeared in the MSD. One needs a large time window to average the MSD or Van Hove function of short time scale to include such contributions. Thereby, it is noteworthy that the time regions are for the mean dynamics of all ions and not for individual ions. Naturally, the results from glass-forming materials such as that in ref 31 on thermodynamic scaling of α-relaxation time stems from the secondary relaxation is for the mean dynamics and not for individual dynamics. As clarified by the comparison with the self-part of the Van Hove functions of the corresponding time region,3−6 a small increase of MSD after this time (the end of NCL) is the beginning of the (successful) jump motion region. Non-decaying Part Obtained from the Van Hove Function Found in the NCL Region at Low Temperatures. The Van Hove function is used for further characterization of the MSD. Characteristics of the MSD in the NCL region at lower temperature regions are examined by the selfpart of the Van Hove function defined as follows:33

Figure 2. (a) Self-part of the Van Hove function of Li ions in Li2SiO3 at 300 K plotted in a linear scale. (b) Corresponding curves plotted in a log−log scale. Curves from inner to outer are for 0.8, 8, 39.2, 80.0, 560, and 1520 ps, respectively. The curves in (a) overlap each other, and the development of the second peak is negligibly small while the time development is clearer in the double logarithmic plot in (b). Contribution at around 3 Å by the first jumps is found at the longest time in the plot.

Cage Decay Function Obtained from the Van Hove Function. The cage decay function Ncage(t) is defined as follows:5 Ncage(t ) ≡

i=1

rc

4πr 2Gs(r , t ) dr

(2)

and is shown in Figure 3. Here, the summation is taken within a distance rc from the initial cage (less than 1.7 Å), which corresponds to the first minimum of the self-part of the Van Hove function observed between 500 and 800 K. The same cut off distance was applied for 100 and 300 K. The area under the curve of the self-part of the Van Hove function is proportional to the number of particles in cages, and the area is normalized to 1 at the initial time. As shown, Ncage(t) at low temperatures do not decay even after long time durations, which is consistent with the behavior in MSD. From the definition of the cage decay function (eq 2), the nondecay part is a measure of the caging of ions. It is noteworthy that the observed lengthening of the NCL frequency region at low temperature regimes is quite similar to those observed in nonionic systems.3,4,32,34 It suggests that treatment and serious consideration of this region found by molecular dynamics (MD) simulations in glass forming materials is essential for a fundamental understanding of the glass transition. Temperature Dependence of Characteristic Times. To learn how each region is lengthened with decreasing temper-

N

Gs( r ⃗ , t ) = (1/N ) ∑ ⟨δ( ri (⃗ t ) − ri (0) − r ⃗ )⟩ ⃗

∫0

(1)

In Figure 2(a,b), the self-part of the Van Hove function at 300 K is shown, where r ≡ |r|⃗ . At temperatures ∼500 K and above, development of the second and further shells by jump motion is clear5,6 within a time scale of observation (∼several ns), while at 300 K the decay of the function is hardly observed in the linear scale plot (Figure 2a), although in the logarithmic plot (Figure 2b), contribution at around r = 3 Å is found. In Figure 2a,b, the second curve from the right shows the functions for the time ≈ tx1. A development of the next peak begins after this time. Therefore, Li ions are almost trapped within original sites during the observation time, and this time region becomes considerably longer at low temperatures than those at higher temperatures. This situation is also clarified by the cage decay function in the following section. C

DOI: 10.1021/acs.jpcc.7b03132 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

of fast and slow ions,5,6,17 as proven by the changes in the fractal dimension of the random walks.4−6 The value at 100 K is out of the range of this figure. If the line for 1/tx1 is extrapolated, the estimated value of tx1 at 100 K is found to be ∼0.1 s. Thus, within the time scale of conventional MD (typically ≈ ns), the flat part of the MSD found is for the NCL in the caged ions regime. These findings are consistent with our explanations of the time development of MSD and the Van Hove functions. In the following sections, further details of individual ionic motions as well as the characteristics of the cages are examined to understand the dynamics of caged ions. Squared Displacement of Each Ionic Motion. Squared displacement of arbitrarily chosen 6 ions at 300 K and those at 500 K are shown in Figure 5. Each ion has a different color in Figure 3. Cage decay measure by the area under the self-part of the Van Hove function normalized to 1 for 900, 800, 700, 600, 500, 300, and 100 K from the bottom to top. The cut off length was chosen to be 1.7 Å (The first minimum of the function). The value decays by the net jump to the neighbor or further shells.

ature, an inverse of each characteristic time is plotted against reciprocal temperature in Figure 4. The lines representing the

Figure 4. Inverse of the characteristic time scales plotted against inverse of temperatures. Lines are for fitted ones by assuming exponential. Absolute values of the slope increase in the order of 1/tx1, 1/tx2, and 1/tdif and therefore the slow down is enhanced at long time scale. The decrease of the value of the 1/tx1 corresponds to the lengthening of the NCL region. The value estimated from the MSD at 300 K is located on the line. Slope for 1/tdif is comparable to the slope for the temperature dependence of the diffusion coefficients (dotted and dashed line, vertical position is arbitrary shifted).

Figure 5. Squared displacements of arbitrary chosen 6 Li ions during the time scale of the NCL region. (a) 300 K and (b) 300 K (data interval is 10 times larger than the upper one. 125 points are used for each ion.) (c) 500 K. (The same number of data points are used as in the medium one.) Large displacements are found intermittently. The motion is not a harmonic vibration. Some ions show large squared displacements, while other ions are strongly localized.

Arrhenius temperature dependences of 1/tx1, 1/tx2, and 1/tdif, as well as the temperature dependence of the diffusion coefficients (a dotted and dashed line). The latter has a different dimension [m2 s−1] and has been shifted vertically for the sake of comparison of its slope with the other quantities. The absolute values of the slopes increase in the order of 1/tx1, 1/tx2, and 1/ tdif, and therefore the slowing down is enhanced at longer time scales. The activation energy in kJ/mol for 1/tx1 is 23.6, whereas that for 1/tx2 is 36.6. The activation energy for 1/tdif is 48.3, and this slope is comparable to (or slightly smaller than) that for diffusion coefficient, 52.7 (obtained for the range from 600 to 2000 K), as shown in Figure 4. The larger activation energy for the diffusion coefficients and tdif than that for individual jumps (observed between tx1 and tx2) is explained by the strong back correlated motion as an average

each figure. The impression gathered from the figure may be different depending on the time intervals of data points used. Figure 5(a,b) shows motions at 300 K during the time scale of NCL for the same ions but with different time intervals of plots. That is in the top figure, 1250 data points are connected by lines, while in the middle Figure 125 data points for each ion are connected. In both parts, the active and quiet states intermittently appeared for the two ions shown in red and blue. Other ions are more localized. Strong heterogeneity of the motion is also observed at 500 K during the time scale of NCL (10 ps) in Figure 5(c). Here strong heterogeneity means coexistence of ions with large displacements and immobile ions. The number of data points is 125 for each ion. Again, strong intermittency35,36 is found. As found, the motion within a cage is not a simple harmonic motion. The intermittent nature of the motion can be understood from the characteristics of jump motions related D

DOI: 10.1021/acs.jpcc.7b03132 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C to a deterministic chaos.4−6,37,38 The deterministic nature of the motion is previously clarified by the phase-space plot of each ionic motion in lithium metasilicate37 and also in ionic liquid.38 The anharmonicity of the caged ion motion was suggested to be the cause of the NCL.4 Shapes of the Density Profile of Localized Ions. The positions of atoms and ions at 300 K during 500 ps (within the NCL region) in an arbitrarily chosen slice with a width of 6 Å in the system are shown in Figure 6(a). Positions of Si and O

some distortions, especially for the sites of Li ions. This situation is shown clearer in Figure 7(a,b) for two examples of

Figure 6. (a) Positions of atoms at 300 K for the slice with a width of 6 Å in Li2SiO3. Positions are accumulated for 500 ps (NCL region). The point is plotted at every 0.4 ps. (b) Positions of atoms at 500 K for the slice with a width of 6 Å in Li2SiO3. Positions of atoms are accumulated for 10 ps (NCL region). The point is plotted at every 0.08 ps. The same number of data points (125) as in (a) for each site is included.

Figure 7. (a,b) Two examples of the set of the data points during NCL region of Li ions with oxygen atoms within neighboring distance 3.1 Å at initial time at 300 K for the slice with a width of 6 Å. (c) An example of the three-dimensional plot of the positions of Li ions during the NCL region, which exhibited the nonspherical shape. Colors in (a,b) are the same as that in Figure 6.

atoms are essentially localized except for some fluctuation to be shown later. The set of points (125 points in each site) for Li in the site is separated from the other ones but tends to show a pattern related to the jump paths observed at higher temperatures or at longer time scale. A similar plot of the data found at 500 K during 10 ps (NCL region) is shown in Figure 6(b). Since the observation time chosen in the two cases scale with the corresponding characteristic times, tx1 and the same number of data points are used, situation in both cases (at 300 and 500 K) are comparable in spite of quite different time scales of the observation times. The shape of the set of points in each potential well (site) in these figures is not necessarily spherical but is ellipsoid-like with

positions visited by Li ions at 300 K in the NCL regime projected onto the xy-planes, where surrounding oxygen atoms are also shown. Obviously, the shape determined is related to the coordination polyhedron made by oxygen atoms and its fluctuation. A three-dimensional plot of positions of Li ions within a cage for another case at 300 K is shown in Figure 7(c), where the ellipsoid-like shape of the cage is clear. In Figure 8(a,b), trajectories of Li ions in cages of other examples at 300 K are shown. The ion near the edge of the site tends to be sharply drawn back within the site. Thicker curves E

DOI: 10.1021/acs.jpcc.7b03132 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C M

φm(r) =

M

∑ φm(|rk − R|) = ∑ k=1

k=1

qk |rk − R|m

(3)

The vector R = r − rc is a relative position of the point P measured from the center of the subregion and qk is a point charge if it is concerned with Coulombic force. Here the interaction potential has the form of inversed power law ∝(1/ rm), including the case of Coulombic potential. When we represent the term within the () in the following expression by εk, 1/2 ⎡ ⎛ 2R·r r 2 ⎞⎤ |rk − R| = R ⎢1 − ⎜ 2 k − k2 ⎟⎥ , ⎢⎣ R ⎠⎥⎦ ⎝ R

R = |R|

Equation 3 can be expanded as follows when |εk| < 1. φm(r) =

∑ k

+ Figure 8. (a,b) Trajectories of localized ions at 300 K for arbitrary chosen two Li ions. Three dimensional positions during 500 ps runs (NCL region) are connected with the interval of 100 steps (time step = 4 fs). For each ion, 1250 positions are included. Thicker curves are results of local polynomial fittings in 3D. Numerical values for axes are in Å.

qk ⎧ ⎞ m 1 m ⎛⎜ m ⎨1 + εk + + 1⎟εk 2 ⎠ Rm ⎩ 2 2! 2 ⎝ 2

⎫ ⎞⎛ m ⎞ 1 m ⎛⎜ m + 1⎟⎜ + 2⎟εk 3 + ···⎬ ⎠⎝ 2 ⎠ ⎭ 3! 2 ⎝ 2

(3′)

At the arbitrarily chosen position r, the potential function between the cell and the particles at the position can be approximately represented as follows: 2R·rk

rk R

( ) = O( ):

Using the relation O φm(r) =

are results after denoising using a local (for each 100 data points) polynomial fitting. A nearly one-dimensional motion, along the direction of the longer axis of the nonspherical shapes, was observed in both cases and other cases. These results mean that the intermittent spike like motions found in Figure 5 are due to such motions along the long axis of the ellipsoid. Therefore, in a longer time scale than vibrations, the ion exhibited an-harmonic motion, and followed by further jumps having the deterministic nature. In both trajectories of Li ions and squared displacements, one can clearly see trials before jump motions due to deterministic motion related to chaos. This description of jumps is not necessarily the same as that of jump motions characterized as the statistical process39,40 without considering either the deformation of lattice or the dynamics of cages. Although a clear double well structure, which is often assumed, was not found, the shape of density of ions means a contribution of dipole (for Li−O or center of mass positions of these species) and quadrupole of the ionic interactions as well as charge−charge interactions. Therefore, the origin and properties of the NCL are closely related to the anharmonicity represented by multipole expansions of the potential field as discussed below. Multipole Expansion of the Potential Field. The treatment of the electric field by a multipole expansion41−43 is general, and the method is effectively used for acceleration of MD simulations.44,45 When the “multipole method” is applied to MD simulation, it is typically used with the “tree method” or the “hierarchy tree method”.46 Here we considered a case of interaction from M particles in a subregion (cell) to the outer point P for the expansion using multipoles as in the literatures.40−48 The center of the cell is represented as rc and the relative positions of particles k (k = 1, 2,...M) to the center are taken as rk. The potential formed at the point P is represented as follows:

R2

Q αβR αRβ Oαβγ R αRβR γ μα R α Z + + ··· m + m+4 m+2 R R R Rm + 6 (4)

In the right-hand side of the equation, Z is a sum of charge within a cell, Z = Σkqk, μα is a dipole moment, μα = Σkqkrkα, and m Qαβ is a quadrupole, Q αβ = ∑k qk 2 {(3rkα − δαβrk2)}. In this expression, μαRα is for a sum of terms with α = x,y,z; while rkα is a α = x,y,z component of vector rk. Thus, one can treat forces from some groups of atoms/ions by multipole expansion. The multipole expansion can also be represented by using spherical harmonics.49 As pointed out by Datta,50−52 the method of expansion is applicable for both interior and exterior (electronic and/or magnetic) fields (see refs 49 and 50 for the comparison of the interior and exterior cases). In the present ionic system, each Li ion is trapped in the cage formed by oxygen atoms and further shells and is rocked by fluctuation of the coordination polyhedron. If we considered the field formed by the cage, then the potential of the cage consists of a sum of multipole terms, and hence is anharmonic. The contents of the potential in NCL regime can be represented by the multipoles. When considering the motion of ions in the cage, anharmonicity of the potential was considered as an origin of nearly constant loss (NCL).3−6 Although the dipole contributes to the dynamics as pointed out in ref 1, the point charge term is not negligible when the total sum of q is not 0. A larger number of terms will be necessary to represent the situation more precisely. Generally, the contribution of dipole (∼1/r2) decays faster than point charge (∼1/r), and that of quadrupole (∼1/r3) decays faster than dipole. Nonspherical shapes of the positions visited by Li ions reflects the contribution of dipoles and quadrupoles, which were formed by the distribution of atoms/ions in the shells. Fluctuation of the Local Environment of the Ion Sites. The local field considered for the caging ion regime is not static but rather dynamic. F

DOI: 10.1021/acs.jpcc.7b03132 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C As previously shown in several works for the present system53 as well as an ionic liquid,54,55 the large displacement of Li ions to the next site is accompanied by the changes of coordination number of surrounding oxygen atoms. A similar result is also observed in a colloidal system56 in spite of the fact that the system is quite different from ionic conductor examined here. This fact implies the generality of these findings. In Figure 9, fluctuations of the coordination number for some arbitrary chosen Li ions at 500 K are shown. The motion

Figure 10. (a) Fluctuation of the coordination number around an arbitrary chosen Li ion at 100 K during 200 ps. (b) similar plot for the different Li ion. Numbers tended to fluctuate around a mean value in both cases. The fluctuation pattern shows an intermittency suggesting a deterministic nature of the fluctuation.

compensated each other. Again, strong intermittency suggests the existence of deterministic chaos governing the caged ion dynamics. It is noteworthy that the behavior of the cage is also strongly related to the more macroscopic mechanical properties of the systems as demonstrated in nanoporous Li disilicate glass, and also that looser caging of Li ions causes larger mobility of ions.13 In many metallic glasses, mechanical properties, such as strain rate and ductility, have been shown to have correlation to the relaxation time of the secondary β-relaxation.14 While the Poisson ratio (the ratio between transverse strain and longitudinal strain) and/or K/G (where, K is elastic bulk modulus and G is the shear modulus) is related to the Debye− Waller factor, which is also known to be related to the nearly constant loss (NCL).56 Ngai et al.14 have discussed the interrelation among these properties based on the Coupling model. Recent MD work for the hard sphere systems treated the exchanges of caging and escaping57 to characterize the dense fluid. Although our attention in the present work is mainly focused on the localization within a cage, it should not be forgotten that the ion dynamics are complex, heterogeneous, and have hierarchy in the dynamics of fast and slow ions,17,58,59 where the cooperativity of ions play important roles.3,4 Ionic motions at short and long times (of short length scale and long

Figure 9. (a) Fluctuation of the coordination number around an arbitrary chosen Li ion at 500 K during 6 ns. The coordination number changes in the beginning and after that it is fluctuated mainly between 5 and 6. (b) Similar plot for the different Li ion during 200 ps. A larger fluctuation was found in this case. During the runs, the position of the Li ion moves within the cage.

is dynamically heterogeneous. Both long time localizations and large fluctuations are found. In the example of Figure 9(a), the mean value of the fluctuation is shifted at an early stage of the observation but a stable state with fluctuation occuring mainly between 5 and 6 continues to hold up after long timespans (∼6 ns). In the example of the mobile central ion over an even shorter time scale of 200 ps, shown in Figure 9(b), the coordination number fluctuates between 3−8. In Figure 10(a,b), similar plots for ionic motions at 100 K are shown. As can be seen by inspection, the cages (coordination polyhedra) fluctuate even at 100 K. In contrast, the self-part of the Van Hove function seems to be rather stable even at higher temperatures (see Figure 2). This is because the Van Hove function is for a net motion among shells, where the forwardcorrelated motion and back correlated motions have G

DOI: 10.1021/acs.jpcc.7b03132 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

(2) Angell, C. A. Mobile Ions in Amorphous Solids. Annu. Rev. Phys. Chem. 1992, 43, 693−717. (3) Ngai, K L. Relaxation and Diffusion in Complex Systems; Springer: New York, 2011. (4) Habasaki, J.; León, C.; Ngai, K. L. Dynamics of Glassy, Crystalline and Liquid Ionic Conductors: Experiments, Theories, Simulations; Springer International: Switzerland, 2017. (5) Habasaki, J.; Ngai, K. L.; Hiwatari, Y. Molecular Dynamics Study of Cage Decay, Near Constant Loss, and Crossover to Cooperative Ion Hopping in Lithium Metasilicate. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2002, 66, 021205. (6) Habasaki, J.; Ngai, K. L.; Hiwatari, Y. Dynamics of Caged Ions in Glassy Ionic Conductors. J. Chem. Phys. 2004, 120, 8195−8200. (7) Ngai, K. L. Properties of the Constant Loss in Ionically Conducting Glasses, Melts, and Crystals. J. Chem. Phys. 1999, 110, 10576−10584. (8) Rivera, A.; León, C.; Varsamis, C. P. E.; Chryssikos, G. D.; Ngai, K. L.; Roland, C. M.; Buckley, L. J. Cation Mass Dependence of the Nearly Constant Dielectric Loss in Alkali Triborate Glasses. Phys. Rev. Lett. 2002, 88, 125902. (9) Lee, W. K.; Liu, J. F.; Nowick, A. S. Limiting Behavior of AC Conductivity in Ionically Conducting Crystals and Glasses: A New Universality. Phys. Rev. Lett. 1991, 67, 1559. (10) León, C.; Rivera, A.; Varez, J.; Sanz, J.; Santamaria, J.; Ngai, K. L. Origin of Constant Loss in Ionic Conductors. Phys. Rev. Lett. 2001, 86, 1279. (11) Habasaki, J.; Ngai, K. L. Molecular Dynamics Simulation of Ion Dynamics in Glassy Ionic Conductors: Evidence of the Primitive Ion Hopping Process. J. Non-Cryst. Solids 2006, 352, 5170−5177. (12) Laughman, D. M.; Banhatti, R. D.; Funke, K. Nearly constant loss effects in borate glasses. Phys. Chem. Chem. Phys. 2009, 11, 3158− 3167. (13) Habasaki, J. Molecular Dynamics Study of Nano-Porous MaterialsEnhancement of Mobility of Li ions in Lithium Disilicate. J. Chem. Phys. 2016, 145, 204503. (14) Ngai, K. L.; Wang, L.-M.; Liu, R.; Wang, W. H. Microscopic Dynamics Perspective on the Relationship between Poisson’s Ratio and Ductility of Metallic Glasses. J. Chem. Phys. 2014, 140, 044511. (15) Ngai, K. L.; Wang, L.-M.; Liu, R.; Wang, W. H. Erratum:“Microscopic Dynamics Perspective on the Relationship between Poisson’s Ratio and Ductility of Metallic Glasses. J. Chem. Phys. 2014, 140, 044511. (16) Yu, H. B.; Shen, X.; Wang, Z.; Gu, L.; Wang, W. H.; Bai, H. Y. Tensile Plasticity in Metallic Glasses with Pronounced β Relaxations. Phys. Rev. Lett. 2012, 108, 015504. (17) Habasaki, J.; Okada, I.; Hiwatari, Y. Origins of the Two-Step Relaxation and the Boson Peak in an Alkali Silicate Glass Studied by Molecular-Dynamics Simulation. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 1995, 52, 2681−2687. (18) Ngai, K. L.; Habasaki, J.; Hiwatari, Y.; León, C. A Combined Molecular Dynamics Simulation, Experimental and Coupling Model Study of the Ion Dynamics in Glassy Ionic Conductors. J. Phys.: Condens. Matter 2003, 15, S1607−1632. (19) Heuer, A.; Kunow, M.; Vogel, M.; Banhatti, R. D. Characterization of the Complex Ion Dynamics in Lithium Silicate Glasses via Computer Simulations. Phys. Chem. Chem. Phys. 2002, 4, 3185−3192. (20) Habasaki, J.; Ngai, K. L. The Mixed Alkali Effect in Ionically Conducting Glasses Revisited: A Study by Molecular Dynamics Simulation. Phys. Chem. Chem. Phys. 2007, 9, 4673−4689. (21) Habasaki, J.; Okada, I.; Hiwatari, Y. MD Study of the Mixed Alkali Effect in a Lithium-Potassium Metasilicate Glass. J. Non-Cryst. Solids 1995, 183, 12−21. (22) Habasaki, J.; Okada, I.; Hiwatari, Y. MD Study of the Mixed Alkali Effect in Terms of the Potential Surface in the LithiumPotassium Metasilicate Glass. J. Non-Cryst. Solids 1996, 208, 181−190. (23) Habasaki, J.; Ngai, K. L. Molecular Dynamics Study of Network Statistics in Lithium Disilicate: Qn Distribution and the PressureVolume Diagram. J. Chem. Phys. 2013, 139, 064503.

scale), which have different exponents to characterize them, tend to mix in multifractal manners. Such a multifractal nature characterizes well the hierarchical dynamics with heterogeneity.59,4 It is noteworthy that multifractal analysis is one of the possible ways to give a quantitative measure of the heterogeneity (see ref 4 and references therein).

■

CONCLUSIONS Explanation of the origin and properties of the caged ion dynamics exemplified as the NCL, and the manner of lengthening of different time regimes of ion dynamics reflected by the mean square displacement (MSD) with decreasing temperature in lithium metasilicate were proposed in previous works,5,6,11 and related to experiments.4,59−62 The new molecular dynamics simulation we performed at lower temperatures than ever before confirm the explanations and interpretations of the nature of the NCL given in the previous works. Observed lengthening of the NCL region is quite similar to that observed in nonionic systems near the glass transition temperature. However, it does not indicate the complete arrest of ions because the Li ion has the size of the diffusivity that is experimentally measurable, even near room temperature or below. For the case of nonionic system, the situation with this nondecaying region may be regarded as the computational glass transition. This behavior should be distinguished from the thermodynamic glass transition because this time scale (typically ≈ ns) is shorter by many orders than the characteristic time scale of the real glass transition experimentally observed. Moreover, motions of Li ions in the NCL region were examined in details. Nonspherical shapes of the density profile of Li ions in the NCL region were visualized, and the relation to a multipole expansion was argued and presented. Several characteristics such as heterogeneous and anharmonic motion with intermittency connected to deterministic chaos, low dimensional anisotropic motion along the long axis of the cage with the ellipsoid shape were visualized and discussed to elucidate the characteristics of the cages in the present work. The central Li ion in the coordination polyhedron (cage) has a tendency to move along with the fluctuation of the coordination number along the long axis with time in the NCL region. For further studies of caging dynamics, more quantitative investigation for contributions of multiple terms in the NCL regime would be useful.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected] (J.H.). ORCID

Junko Habasaki: 0000-0002-2887-2340 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS A part of this work was supported by JSPS KAKENHI Grant Number JP17K05570.

■

REFERENCES

(1) Ngai, K. L. A review of critical experimental facts in electrical relaxation and ionic diffusion in ionically conducting glasses and melts. J. Non-Cryst. Solids 1996, 203, 232−245. H

DOI: 10.1021/acs.jpcc.7b03132 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

(49) Thompson, W. J. Angular Momentum; John Wiley & Sons Inc.: New York, 2008. (50) Datta, S. Multipole Expansion of the Interaction Hamiltonian between a Charged Particle and a Non-Uniform Magnetic Field. Eur. J. Phys. 1984, 5, 243. (51) Datta, S. Multipole Expansion of the Interaction Hamiltonian between a Charged Particle and a Non-Uniform Magnetic Field: II. Exterior Field. Eur. J. Phys. 1985, 6, 257. (52) Dieterich, W.; Maass, P. Constant Dielectric Loss in Disordered Ionic Conductors: Theoretical Aspects. Solid State Ionics 2009, 180, 446. (53) Habasaki, J. Molecular-Dynamics Study of Glass Formation in the Li2SiO3 System. Mol. Phys. 1990, 70, 513−528. (54) Habasaki, J.; Ngai, K. L. Rigidity and Soft Percolation in the Glass Transition of an Atomistic Model of Ionic Liquid, 1-Ethyl-3Methyl Imidazolium Nitrate, from Molecular Dynamics Simulations Existence of Infinite Overlapping Networks in a Fragile Ionic Liquid. J. Chem. Phys. 2015, 142, 164501. (55) Zhang, Y.; Maginn, E. J. Direct Correlation between Ionic Liquid Transport Properties and Ion Pair Lifetimes: A Molecular Dynamics Study. J. Phys. Chem. Lett. 2015, 6, 700. (56) Weeks, E. R.; Weitz, D. A. Properties of Cage Rearrangements Observed near the Colloidal Glass Transition. Phys. Rev. Lett. 2002, 89 (2002), 095704. (57) Ngai, K. L.; Habasaki, J. An Alternative Explanation of the Change in T-Dependence of the Effective Debye-Waller Factor at Tc or TB. J. Chem. Phys. 2014, 141, 114502. (58) van Megen, W.; Schöpe, H. J. The Cage Effect in Systems of Hard Spheres. J. Chem. Phys. 2017, 146, 104503. (59) Habasaki, J.; Okada, I.; Hiwatari, Y. Fracton Excitation and Lévy Flight Dynamics in Alkali Silicate Glasses. Phys. Rev. B: Condens. Matter Mater. Phys. 1997, 55, 6309. (60) Kanert, O.; Steinert, J.; Jain, H.; Ngai, K. L. Nuclear Spin Relaxation and Atomic Motion in Inorganic Glasses. J. Non-Cryst. Solids 1991, 131−133, 1001−1010. (61) Rivera, A.; León, C.; Sanz, J.; Santamaria, J.; Moynihan, C. T.; Ngai, K. L. Crossover from Ionic Hopping to Nearly Constant Loss in the Fast Ionic Conductor Li0.18La0.61TiO3. Phys. Rev. B: Condens. Matter Mater. Phys. 2002, 65, 224302. (62) León, C.; Rivera, A.; Ngai, K. L. Correlation between Ion Hopping Conductivity and Near Constant Loss in Ionic Conductors. Phys. Rev. B: Condens. Matter Mater. Phys. 2004, 69, 134303.

(24) Habasaki, J.; Okada, I. Molecular Dynamics Simulation of Alkali Silicates Based on the Quantum Mechanical Potential Surfaces. Mol. Simul. 1992, 9, 319−326. (25) Gilbert, T. L. Soft-Sphere Model for Closed-Shell Atoms and Ions. J. Chem. Phys. 1968, 49, 2640−2641. (26) Ida, Y. Interionic Repulsive Force and Compressibility of Ions. Phys. Earth Planet. Inter. 1976, 13, 97−104. (27) Gonçalves, G. V.; Rino, J. P. Diffusion Mechanisms in Lithium Disilicate Melt by Molecular Dynamics Simulation. J. Non-Cryst. Solids 2014, 402, 91−95. (28) Nascimento, M. L. F.; Fokin, V. M.; Zanotto, E. D.; Abyzov, A. S. Abyzov. Dynamic Processes in a Silicate Liquid from above Melting to below the Glass Transition. J. Chem. Phys. 2011, 135, 194703. (29) Habasaki, J.; Ueda, A. Molecular Dynamics Study of One Component Soft-Core system - Analytical Expression of NonEquilibrium Relaxation in Constant Pressure Conditions. J. NonCryst. Solids 2016, 447, 212. (30) Habasaki, J.; Ngai, K. L. Molecular Dynamics Simulation of Ion Dynamics in Glassy Ionic Conductor: Evidence of the Primitive Ion Hopping Process. J. Non-Cryst. Solids 2006, 352, 5170−5177. (31) Ngai, K. L.; Habasaki, J.; Prevosto, D.; Capaccioli, S.; Paluch, M. Thermodynamic Scaling of α-Relaxation Time and Viscosity Stems from the Johari-Goldstein β-Relaxation or the Primitive Relaxation of the Coupling Model. J. Chem. Phys. 2012, 137, 034511. (32) Ngai, K. L.; Habasaki, J.; Prevosto, D.; Capaccioli, S.; Paluch, M. Erratum: “Thermodynamic Scaling of α-Relaxation Time and Viscosity Stems from the Johari-Goldstein β-Relaxation or the Primitive Relaxation of the Coupling Model” [2012, J. Chem. Phys. 2012, 137, 034511]. J. Chem. Phys. 2014, 140, 019901. (33) Van Hove, L. Correlations in Space and Time and Born Approximation Scattering in Systems of Interacting Particles. Phys. Rev. 1954, 95, 249−262. (34) Kubo, E.; Minoguchi, A.; Sotokawa, H.; Nozaki, R. Nearly Constant Dielectric Loss of Polyhydric Alcohols. J. Non-Cryst. Solids 2006, 352, 4724−4728. (35) Townnsend, A. A. Local Isotropy in the Turbulent Wake of a Cylinder. Aust. J. Chem. 1948, 1, 161. (36) Wang, X.-J. Dynamical Sporadicity and Anomalous Diffusion in the Lévy Motion. Phys. Rev. A: At., Mol., Opt. Phys. 1992, 45, 8407. (37) Habasaki, J.; Ngai, K. L.; Hiwatari, Y. Time Series Analysis of Ion Dynamics in Glassy Ionic Conductors Obtained by a Molecular Dynamics Simulation. J. Chem. Phys. 2005, 122, 054507. (38) Habasaki, J.; Ngai, K. L. Heterogeneous Dynamics of Ionic Liquids from Molecular Dynamics Simulations. J. Chem. Phys. 2008, 129, 194501. (39) Eyring, H. The Activated Complex in Chemical Reactions. J. Chem. Phys. 1935, 3, 107−115. (40) Vineyard, G. H. Frequency Factors and Isotope Effects in Solid State Rate Processes. J. Phys. Chem. Solids 1957, 3, 121. (41) Thorne, K. S. Multipole expansions of gravitational radiation. Rev. Mod. Phys. 1980, 52, 299. (42) Newman, E. T.; Penrose, R. 10 Exact Gravitationally-Conserved Quantities. Phys. Rev. Lett. 1965, 15, 231. (43) Kosov, D. S.; Popelier, P. L. A. Convergence of the Multipole Expansion for Electrostatic Potentials of Finite Topological Atoms. J. Chem. Phys. 2000, 113, 3969. (44) Barnes, J. E. L.; Hut, P. A Hierarchical O(N log N) ForceCalculation Algorithm. Nature 1986, 324, 446−449. (45) Greengard, L.; Rokhlin, V. A Fast Algorithm for Particle Simulations. J. Comput. Phys. 1987, 73, 325. (46) Giese, T. J.; York, D. M. Extension of Adaptive Tree Code and Fast Multipole Methods to High Angular Momentum Particle Charge Densities. J. Comput. Chem. 2008, 29, 1895. (47) Carrier, J.; Greengard, L.; Rokhlin, V. A Fast Adaptive Multipole Algorithm for Particle Simulations. SIAM, J. Sci. Stat. Comput. 1988, 9, 669. (48) Ueda, A., Molecular SimulationsFrom Classical to Quantum Methods-, Shoka-bo 2003 (in Japanese). I

DOI: 10.1021/acs.jpcc.7b03132 J. Phys. Chem. C XXXX, XXX, XXX−XXX