Relaxation process of hopping ionic conduction in lattice gas models

SOLID STATE IONICS

ELSEVIER

Solid State Ionics 79 (1995) 3-8

Relaxation process of hopping ionic conduction in lattice gas models Hiroshi Sato, Anuradha Datta, Takuma Ishikawa

*

School of Materials Engineering, Purdue University, W. Lafayette, IN 47907-1289, USA

Abstract An analytical, microscopic treatment has been developed for the relaxation process of hopping ionic conduction in interacting lattice gas systems to gain understanding of the origin of the universal dynamic response in solid electrolytes. Here, it is demonstrated that the relaxation motion of mobile ions is intrinsically non-Debye and follows the KohlrauschWilliams-Watts relaxation behavior even without interactions among mobile ions. This is an apparent contrast to the current view that the anomaly is “created” by interaction and disorder. Keyworak

Hopping ionic conduction; Kohlrausch-Williams-Watts

1. Introduction

A major interest of this treatment is to understand the “universal dynamic response” observed in the ionic conductivity in solid electrolytes by an analytical treatment from first principles. The phenomenon “ universal dynamic response’ ’ is succinctly represented by the following two empirical relations [1,2]; the first indicates that the frequency dependent ionic

conductivity power law

a(w)

U( 0) u UP

is represented

by Jonscher’s

(O
(1) with p limited between 0 and 1, and second, the relaxation process in the induced current when a constant external field is switched on at t = 0 in such systems is non-Debye (non-exponential) with the

* Faculty of Engineering, Tokyo Institute of Polytechnics. Atsugi, Kanagawa 243-02, Japan. 0167-2738/95/$09.50

relaxation behavior; Lattice gas models; Bethe method

behavior to reach the stationary state represented by the Kohlrausch-Williams-Watts (KWW) stretched exponential function B(t) Bl a(t)

(1-p=p,

-(t/r)“], o
)*

(2)

The relaxation process in the hopping conduction in a lattice gas is generally considered to be of Debye type [3], and prior theoretical efforts such as the coupling model, Ngai et al. [4], the cage effect jump relaxation model, Funke et al. 121, and the diffusion controlled relaxation model, Elliott et al. [5], have ascribed the origin of the nonlinear relaxation process to a coupling with the relaxation process of the surroundings through interactions with neighbors. It has been argued that the effect should arise from interactions and disorder [6]. However, questions still remain as to what specifically creates the KWW behavior. We have developed a microscopic, analytical treatment of the relaxation process

0 1995 Elsevier Science B.V. All rights reserved

SSDI 0167-2738(95)00020-8

-exp[

H. Sate et al. /Solid State Ionics 79 (1995) 3-8

4

of ionic conduction in a lattice gas model which can explain characteristics of universal dynamic response, and this shows that the relaxation motion of a single particle in the assembly of particles of the same species with negligibly small amount of vacant sites (the simplest lattice gas model) is intrinsically non-Debye, and follows the KWW relaxation behavior even without interaction in contrast with the current view indicated above. We try to show the origin of why such a discrepancy should arise here.

2. Derivation of the relaxation curve The task involved here is to derive the relaxation curve for the induced current upon the application of a constant external field at t = 0 in the lattice gas model. Here, we limit ourselves to the derivation of flow of a single particle in an assembly rather than that of the assembly, because the treatment of the former is far simpler. The relation between the two is now well known [7,8], but, as far as the mechanism of the relaxation process is concerned, the former is satisfactory as well. Typically, the relaxation curve can be divided into three time regions; adiabatic, intermediate and isothermal [9] as shown in Fig. 1. If the relaxation curve is normalized, the curve represents the timedependent correlation factor as introduced by Funke and the relaxed value in the isothermal region for flow of a single particle is commonly called the tracer correlation factor. The difficulty in deriving

I

0

0

'

correlation factor t

Fig. 1. Normalized relaxation curve. The quantity f*(t) sents the time-dependent correlation factor.

repre-

the relaxation behavior theoretically is to deal with all the three time regions with equal rigor. By using a combination of the Path Probability Method (PPM) of irreversible statistical mechanics [lo] (which is a systematic method to deal with the microscopic master equation approach) and the Bethe method (applied in time space) [ll], we show how the relaxation behavior can be consistently and analytically treated in the pair approximation from t = 0 to t --f ~0. In a conventional linear response theory, the flow is evaluated based on the “reference state” (usually taken ai the equilibrium state in the isothermal range). If the “reference state” does not change with time, the relaxation of flow is essentially exponential. The appearance of the KWW behavior must somehow be related to the change of the reference state in the intermediate region. This change cannot ordinarily be calculated analytically, because the evaluation of the nonequilibrium state is involved. Therefore, earlier workers calculated this as a realistic change of the state from the adiabatic to the isothermal time range as a result of coupling with the relaxation process of surroundings through interactions with mobile ions, and the induced current is evaluated as if the whole range were the isothermal range. This process is especially clearly shown in the treatment of Funke [2]. This kind of treatment thus generally requires models for the coupling as well as interactions among mobile ions. However, as it is shown, KWW behavior is created by the development of fluctuation from t = 0 to the equilibrium distribution as t -+ ~0,and the KWW behavior can be considered to be intrinsic for the motion of a particle in a lattice gas whether there are interactions among mobile ions or not. Here, we show how the relaxation curve can be directly calculated without going through the calculation of the reference state. We start from the derivation of the flow of a small number of tracer particles in an assembly of particles of a multicomponent system in a lattice gas with nearest neighbor pairwise interactions among particles to derive the relaxation process by the PPM. The hopping of particles is allowed only to nearest neighboring vacant sites. When a constant driving force for the ith species (which is identified to be the tagged species) is applied at t = 0, the normalized flow of particles of the ith species in the disordered distribution (without the long range order), ?Pj(t), to

5

H. Sato et al. /Solid State Ionics 79 (1995) 3-8

the first order of the driving force in the pair approximation, is expressed as Y&(t) = -hi+

xQj$ji(t).

(3)

Here, Qj specifies the equilibrium distribution of particles represented by the concentration of the jth species around the lattice site on the average, qji indicates the gradient in the distribution of particles created by the driving force for the ith species, ki. Thus, Qj essentially represents the equilibrium distribution taken as the reference state for the linear expansion. In addition, in the PPM, Qj represents the equilibrium distribution of particles calculated by the Cluster Variation Method of statistical thermodynamics (CVM, a static version of the PPM) and is independent of time. The equation has been solved in the limit t + 03 (stationary state) analytically by Kikuchi for the stationary values of Ti and +ji in terms of a time independent function Zi derived as [I21 qi=

-{[(2&-1)zi]/[2+(2ij-3)zi]}~~i, (4)

where Zi is expressed as Zi = c ( QjGj)/( gi + Gj).

(5)

Here rGi, etc. are defined as the jump frequency of ith species wi with the bond breaking factor due to the nearest neighbor pair interactions among particles, and 21% represents the coordination number of the lattice. Then, the tracer correlation factor fi is defined as the coefficient of ki in Eq. (4), f,= [(2&l)Z,]/[(l-Zi)

+(2Gl)Z,].

(6)

The physical meaning of the calculation of flow by the pair approximation of the PPM is clear from Fig. 2 which represents Eqs. (3), (4) and (5). Here, the situation after the tagged atom has jumped out (at t = O), is represented and its return probability into the vacancy is calculated in competition with the other particles around the vacancy whose species are specified by Qj. The quantity (1 - Zi) represents the jump back of the tagged particle while (2ij - l)Zi represents the escape of the vacancy. Therefore, the normalized jump back probability in the stationary state defined as Ti is Ti=

(1 _Zi)/[(l

-zi)

and in terms of

TV, fi

f, = (1 - Ti)/(l

+ (2ij-

P,=l-z

Fig. 2. Definition of Z in terms of I’, and P,.

(7)

is given as

+ Ti).

(8)

Fig. 2 thus confirms that the calculation of flow by the PPM is equivalent to the time correlation of the motion of the tagged particle from the time the particle has jumped out leaving the vacancy behind. Since Qj in Eq. (3) is time independent, Eq. (3) is equivalent to the equation for the stationary state at any time instant t. Therefore, as in Eqs. (7) and (8), we can define time dependent quantities Ti*(t), fi* (t) and Zi*(t) which tend to Ti, fi and Zi respectively as t -+ w as Ti*(t)

=

[l -z;(t)]/[(l +(2G,-

-z;(t))

l)z;(t)]

(9)

and fi*(t)

= [l-

Ti*(t)]/[l

+ Ti*(t)].

(IO)

The quantity f *(t) thus defined in Eq. (10) is called the time-dependent correlation factor and represents the normalized relaxation curve which we would like to derive. Among the time dependent quantities T * (t), f * (t) and Z *(t), the relation between T * (t) and T can be directly calculated based on Fig. 2 in an elementary fashion. The result yields T*(t)=Ti[l-exp(-t?it/~i)].

2;yz,

l)Zi]

(II)

Eq. (11) indicates that the relaxation process T * (t) is exponential if 7i is time independent as in the present case while, if 7i increases with time, it becomes KWW type. Since the value of 7i is determined by Qj, 7i can be indirectly time dependent if the quantity corresponding to Qj changes with time.

6

H. Sato et al. /Solid State lonics 79 (1995) 3-8

The analysis of the time dependence of ri is thus required to determine the cause of the KWW behavior. In the PPM, since Qj is constant from t = 0 to t = ~0,the same 7i represents both the adiabatic state and the isothermal state. As is clear from Fig. 2, on the other hand, the pair approximation of the PPM represents a correct adiabatic (t = 0) value of ri, because, here, only the knowledge of the direct neighbor of the vacancy is involved. In addition, the fact that Qj in Eq. (3) is an averaged value and is the same from t = 0 to t = 00in the PPM corresponds to the mean field concept in r-space. Therefore, the task is to improve the adiabatic value of 7i at t = 0 obtained by the PPM to the correct isothermal value corresponding to the pair approximation as t + ~0 through the intermediate region in the form of improving the isothermal value of the mean field approximation. The original Bethe method [lo] is known to correspond to the pair approximation in the Cluster Variation Method of statistical thermodynamics on which the PPM is based. In addition, the consistency relation [lo], which connects the thermodynamical properties of a small cluster to that of the outside gives the propagation of order with distance from the local order to the long-range order in the treatment of the order-disorder transformations [13]. Therefore, if this concept is applied to the calculation of 7i in t space as in this case, it is expected that the isothermal value of ri at t = CCcan be obtained from the adiabatic value of ri at t = 0 in the pair approximation in the form of the propagation of order. The introduction of the Bethe method to improve the degree of approximation of 7i through Z in Eq. (4) has been called the conversion process because this process corresponds to the conversion of the ensemble-averaging process taken in the PPM into the time-averaging process [12,14]. This is accomplished by following the fluctuation in Q, with time. We show in the following specifically the motion of a single particle in an assembly of particles of the same species without interactions and with a negligible amount of vacant sites. This illustrates the process of how the KWW behavior is introduced. The treatment for general cases such as for a dilute (a large number of vacancies) case or for binary systems can be carried out in a similar fashion [14]. A

further improvement of this treatment was made recently [15]. We start from a group of particles of nearest neighbor sites as in Fig. 2 (but of the same species) which represents the cluster in the pair approximation of the PPM at t = 0. Then, Z in Eq. (5) for the assembly of particles of the same species is Z=w/(w+w),

(12)

where w represents the jump frequency of the particle. From this, the adiabatic value of r is obtained as r = w/[ w + (2 0 - l)w]. In the Bethe treatment, since particles other than the tagged particle in the group are thought to be averaged particles having the same jump efficiency T as the tagged particle, we obtain for Z Z=[w(l-7)]/[w+w(l-T)] = (1-

7)/(2-

T).

(13)

From Eq. (13), using Eq. (71, the consistency relation in T is obtained in the form of an iteration formula as 7= 1/[2&-

(20-

l)r].

(14)

Eq. (14) can be written in a differential form as dr dN

2[2i,-(2&1)T][1-(2c5-1)T][1-s] (29-1)+[2;-(20-1)T]2

’

(15)

where N represents the number of iterations. As the number of iterations, N, increases, T approaches the isothermal value r + l/(2 L - 1) (N + a) from the adiabatic value r= l/(2&) (N = 1) in Eq. (14). Therefore, if N is related to time t, Eq. (15) represents the propagation of order, dr/dt, in t space through dr/d N. The number of iterations N represents the distance of propagation of information, and, hence, represents the number of jumps n which effectively covers N unit distance by the vacancy which makes a random walk motion. If so, N and the number of jumps n are related as N = 6 and t = n/w. Because N = 1 corresponds to t = 0, N - 1 = & [16]. This is the relation between N and t. The propagation of order in Eq. (14) means that information at the center of the group in Fig. 1 is accumulated over a distance with time by hopping of particles [17] represented by the motion of the vacancy. Here, specifically, the increase of r indicated in Fig. 3

H. Sato et al. /Solid State Ionics 79 (1995) 3-8

represents the return of the escaped vacancy with time. The presence of interactions among conduction ions, especially those which tend to create an ordering of conduction ions enhances the KWW behavior. In Fig. 3, r-t relations are plotted for different values of 2 ij with 26, = 2 for the one dimensional lattice, 3 for the two-dimensional honeycomb lattice, 4 for the two-dimensional square net and 8 for the body-centered cubic lattice. The increase of 7 with t means that 7 *(t) and, hence f*(t) also have the KWW behavior (Eq. (11)). In other words, the non Debye behavior is inherent although the effect is almost negligible except for the case, 2ij = 2. Based on Eq. (1.5), the increase of 7 as t + m for cases with 26, > 2 is exponential, but for the one-dimensional case (2 ij = 2), where the percolation limit is reached as t + 03, the decay becomes proportional to t-‘/’ as the percolation limit is approached. Therefore, in this case, there is a changeover in slope in the 7-t relation as is clear in Fig. 4, which shows d(log f * (t))/d(log t>. The extensive linearity in the relation log f * (t) - log t such as that observed for the case 2 ij = 2 is also observed in systems where f*(t) decreases with time to a small enough value such as in mixed alkali systems. This is due to the existence of a scaling law for the decrease of fi* (t) with time near the percolation limit. The relation log f * (t) - log t can be converted to log f * ( o) - log o (o represents the frequency) by simply converting t into 27r/o. The Jonscher power law exponent, p, can be obtained from either

1.00

,

0.20

-

0.00

’

I

f

0.0

2ti=8

, 10.0

I

20.0

30.0

40.0

I 50.0

wt Fig. 3. 7-t relations numbers, 2 6.

for structures

with different

coordination

3 *-ru 2

0.1 0.01

0.1

1

10

100

loo0

log wt Fig. 4. Comparison of Monte Carlo simulation analytical theory for 2 6 = 2 and for 2 cj = 3.

results with the

d(log f *(t))/d(log t) = -p or d(log f * Cm))/ d(log w) = p. The changeover of p with log t as a function of f* such as that observed for the case 26, = 2 is commonly observed also in the calculated results of general cases with interaction. The importance of this result in understanding the “universal dynamic response” will be discussed elsewhere. The relation, log f *(t) - log t, is shown for 2 G, = 2 and for 26, = 3 in Fig. 4 along with the results of Monte Carlo simulation of the time dependent correlation factor [18]. The simulation calculates f*(t) = -/d is excellent. In the case 2 G, = 3, where the tracer correlation factor obtained by the pair approximation coincides with the exact value of l/3, an almost exact agreement of the results is observed (Fig. 4). This shows that the pair approximation is quite acceptable for the calculation of

8

H. Sate et al. /Solid

f *(t). In the one-dimensional case, a change in slope in the log f*(t) - log t curve is observed at the value f*(t) N 0.1 (Fig. 4) in agreement with the theory. The values of p obtained from log f*(t) log t for 26, = 2, 3, 4, and 8 are 0.5, 0.32, 0.22 and 0.12 respectively although these values do not have any practical importance. The value 0.5 for the onedimensional lattice is, however, proven to be the exact value for this case. The agreement of the results of the simulation and the analytical results for the assembly of particles of the same species let us believe that the KWW relaxation behavior is an intrinsic characteristic of the motion of particles in lattice gas models and is a result of the development of fluctuation with time from t = 0. Therefore, no specific model for the appearance of the KWW behavior is required.

Acknowledgements This work is a result of cooperation with K. Funke of the University of Miinster, Germany. The authors appreciate his cooperation and comments. The authors also owe valuable discussions to R. Kikuchi of UCLA and G. L. Lied1 of Purdue University. This work was supported by the U.S. Depart-

State Ionics 79 (1995) 3-8

ment of Energy 84ER45133.

under

grant

number

DE-EGOZ

References [11 A.K. Jonscher, Nature 267 (1977) 673. (21 K. Funke, Prog. Solid State Chem. 22 (1993) 195. [31 J.C. Dyre, J. Appl. Phys. 64 (1988) 2456. [41 K.L. Ngai, A.K. Rajagopal and S. Teitler, J. Chem. Phys. 88 (1988) 5086. El S.R. Elliott and A.P. Owens, Phys. Rev. B 44 (1991) 47. 161J.P. Bouchard and A. Georges, Phys. Rep. 195 (1990) 127. [71 G.E. Murch and J.C. Dyre, CRC Crit. Rev. 1.5 (1989) 348. (81H. Sato and T. Ishikawa, to be published. 191 C. Zener, Elasticity and Anelasticity of Metals (University of Chicago Press, Chicago, 1948). [lOI R. Kikuchi, Prog. Theor. Phys. (Kyoto), Supplement No. 35 (1966) 1. WI H.A. Bethe, Proc. Roy. Sot. 150 (1935) 552. 1121 H. Sato and R. Kikuchi, Phys. Rev. B 28 (1983) 648. 1131 F. Zemike, Physica 1 (1940) 565. [141 H. Sato, in: Nontraditional Methods in Diffusion, eds. G.E. Murch, H.K. Bimbaum and J.R. Cost (The Metallurgical Society of AIME, 1984) pp. 203-235. [151 H. Sato, A. Datta and T. Ishikawa, to be published. (161H. Sato and A. Datta, Solid State Ionics 72 (1994) 19. [171 H. Sato, T. Ishikawa and K. Funke, Solid State lonics 53-56 (1992) 907. [I81T. Ishikawa and H. Sato, Lo be published.