A study of diffusion in a disordered chain: Computer simulation

JOURNAL OF SCIENCE OF HNUE Mathematical and Physical Sci., 2012, Vol. 57, No. 7, pp. 142-150 This paper is available online at A STUDY OF DIFFUSION IN A DISORDEREDCHAIN: COMPUTER SIMULATION Trinh Van Mung1, Pham Khac Hung2 and Nguyen Thi Thanh Ha2 1Vinh Phuc College, Vinh Phuc Province 2Department of Computational Physics, Institute of Engineering Physics, Hanoi University of Science & Technology Abstract.We have examined two simulation techniques for diffusion on disordered systems: the non-force and force method. The latter enables one to study diffusion on the system with periodic boundary conditions and size which is 20 times shorter than the non-force method. Furthermore, the force method is an appropriate way to study impurity diffusion in a more real model of amorphous solid constructed by molecular dynamic or statistic relaxation techniques. Using the non-force method we simulate diffusion on a disordered chain consisting of 2000 sites for a wide range of temperature. Three types of disordered chains, site, saddle and site-saddle system are considered. The influence of site and saddle disorders and Arrhenius behavior is also investigated and discussed in this work. Keywords: Diffusion, random walk, disordered system, simulation, non-force method. 1. Introduction Despite great effort over a long period of time, many aspects of diffusion in amorphous alloys (AM) are still poorly understood [1-3]. For example, it is still not clear why the existence of a wide, continuous spectrum of site and saddle point energy for diffusion process could also give rise to the Arrhenius behavior experimentally observed in certain AMs [4-6]. In comparison with crystal counterparts, the pre-exponential factor as well as activation energy exhibits certain specific features. For instance, the factor Do for diffusion of some impurities in AMs changes in the very wide range of 10±7cm2/s [2-3]. Several models have been suggested to clarify those observations. One, the Fisher random walk model, states that, for not too strong disorder, random walk in dimensions Received October 15, 2011. Accepted September 30, 2012. Physics Subject Classification: 62 44 01 01. Contact Trinh Van Mung, e-mail address: mungtv76@gmail.com 142 A study of diffusion in a disordered chain: computer simulation higher than two gives rise to classical diffusion behavior at large time [5]. A simulation of diffusion in regular disordered lattice in which the sites for target particle remain ordered but the potential barriers vary from site to site could provide useful information [7-8]. As shown in ref. [7], the mean square displacement 〈x2n〉 of the target particle after n jumps is proportional to fna2. Here a is the distance between the two nearest sites and f is correlation factor which is a function of temperature and may be less than 10−3 for strong disordered lattice. The Monte-Carlo (MC) simulation in [8] found a compensatory effect between disorders in the site and saddle point energies which originates in the Arrhenius behavior for diffusion of strong disordered lattice. However, it is very useful to conduct a simulation of the diffusion process in a more real model of AMs which is constructed by techniques like molecular dynamic (MD) or statistic relaxation method [9-10] because these models reproduce well both the dynamic and structural properties of AMs. The network of sites that are a link or a potential barrier between them can be directly calculated using MD models. Nevertheless, previous MC technique reported in refs. [8] can not be applied for that network due to the very long computation required and the imprecise result obtained. In the case of low temperature, for example, it is very impractical because one must run a very large number of particle hops. The aim of this article is twofold. First, it is to find a simulation technique which could be use to study systems with periodic boundary conditions and are of moderate size. Two simulation techniques based on target particle probability density and the Einstein diffusion equation have been investigated. We demonstrate the effectiveness of suggested simulation techniques on diffusion in a definite disordered chain with and without periodic boundary conditions, and compare this to the Monte-Carlo method. Second, the influence of two types of energetic disorders on diffusion in disordered chains is also considered and examined here. 2. Content 2.1. Calculation method As mentioned above, it is useful to study the diffusion of impurities on MD model sites. A simulation for this system can be performed by means of the MC method. Accordingly, the target particle that is initially located somewhere in the MD model moves from one site to another and visits a large number of sites along a random path. The position of the target particle is recorded during every run and then an average square displacement and diffusion time is calculated after n jumps. Such a simulation is not accurate and also not practical at low temperatures when the particle is captured by "deep trap" and moves many times between two neighboring sites. Now we examine two simulation techniques which consume much less computing time. Consider a chain of sites labeled by integer numbers 1, 2, 3, ...N , where each ith site is assigned a site energy εi and two neighboring sites i, i + 1 the transition energy εi,i+1 (saddle point 143 Trinh Van Mung, Pham Khac Hung and Nguyen Thi Thu Ha energy). The target particle locates initially at middle siteN/2 and then will be allowed to randomly move to left or right sites by a thermally activated Poissonian jump. The relative probability to jump from ith site into left and right neighboring site is given by pi = 1− qi = exp(−εi,i−1β) exp(−εi,i+1β) + exp(−εi,i−1β) . (2.1) Since the time period between hops from site i to site i+1 is τi,i+1 = τ0 exp((εi,i+1− εi)β) τi,i+1 = τ0 exp((εi,i+1 − εi)β), the average time that the particle stays at ith site can be given as τi = piτi,i−1 + qiτi,i+1 → τi = 2τ0 exp(−εiβ) exp(−εi,i+1β) + exp(−εi,i−1β) (2.2) here τ0 is the period of atomic vibration at a lattice site. Let us define c(i, n) as the probability of just arriving on the ith site as a result of n jumps. The probability at the next step n + 1 is given by c(i, n+ 1) = c(i+ 1, n)pi+1 + c(i− 1, n)qi−1. (2.3) As a initial configuration related to step n = 0, the probability c(i, 0) is equal to 1 for site N/2 and c(i, 0) = 0 for other sites. Eq. 2.3 allows the calculation c(i, n) at each step n and as such simulates the diffusion process. We also calculate average square displacement 〈x2n〉 and the time period to realize n steps tn obtained as follows = a 2 N∑ i ci(i, n)(i−N/2)2 (2.4) and tn = n∑ k N∑ i ci(i, k)τi (2.5) where a is the distance between two neighboring sites. It is now easy to calculate diffusion coefficient D that is determined as a slope of dependence of 〈x2n〉 versus tn. The algorithm above which described and denotes the non-force method can be used to simulate very long chain of sites, but it is not appropriate for systems of moderate size or periodic boundary conditions. To take this into account we employ a force method. Taking a chain withN sites and periodic boundary conditions, 144 A study of diffusion in a disordered chain: computer simulation the transition energy of site 1 and N will be equal to each other, i.e. ε0,1 = εN,N+1. Due to periodic boundary conditions, the target particle arriving at site N can jump to site N − 1 and 1. Analogously, after locating at the 1th site the target particle will jump into site 2 or N at next step. For a large number of steps n, the probability c(i, n) approaches equilibrium value, i.e. the process is stationary, given as ceq(i, n) = exp(−εi−1,iβ) + exp(−εi,i+1β) 2 N∑ 1 exp(−εi,i+1β) . (2.6) If we apply external field force g, then the site and transition energies of the considered system change according to εfi = εi + gia ε f i,i+1 = εi,i+1 + g(i+ 0.5)a (2.7) here the index f indicates the corresponding parameter of applied force g. Obviously the probability c(i, n) varies with n in accordance to Eq. 2.1 with parameters pfi , q f i determined by Eq. 2.7. Force g gives rise to drift flow in the system for which velocity follows Einstein diffusion relation v = Dg.β (2.8) To calculate the drift velocity we determine the averaged mean displacement 〈xn〉 given as 〈xn〉 = a N∑ i c(i, n)pi − c(i, n)qi (2.9) In the absence of force g the quantity 〈xn〉 is almost equal to zero. The drift velocity can be obtained by determining the slope of the dependence of 〈xn〉 versus tn. As such, the force method consists of following steps: (i) Set up set of site and transition energies ǫi, ǫi,i+1 on the chain; i = 1, 2, ...N . Then the initial probabilities c(i, n) are determined by Eq. 2.6. (ii) Set up set of site and transition energies ǫfi , ǫ f i,i+1 using Eq. 2.7. (iii) Calculate the probability c(i, n) at each step n according to Eq. 2.3 and using the parameters pfi , q f i . (iv) Calculate the mean displacement 〈xn〉 and compute the drift velocity, and finally the diffusion coefficient is determined by the Einstein relation 2.8. This method can be easily generalized for higher-dimension systems using a similar computing procedure. Hence, for a network of sites that consists of several thousand sites and is calculated using MDmodels, the force method is a more effective way to calculate the diffusion coefficient of such a system than the non-force or MC method. 145 Trinh Van Mung, Pham Khac Hung and Nguyen Thi Thu Ha 2.2. Results and discussion In order to check the validity of the force method, two chains consisting of 4.105 and 2.103 sites are constructed. A simulation is carried out for the long chain using the non-force method. The force method is applied for the short chain with periodic boundary conditions. Both sites and saddle points of the chain will be assigned to energies chosen at random or constant lowest value ǫ1x. We consider two types of energy distribution: two-level distribution where site/transition energy for every site amounts to either values ǫ1x or ǫ2x (ǫ1x < ǫ2x, x is denoted to s or t) and uniform distribution where site/transition energies are uniformly distributed in the interval [ε1x, ε2x]. Figure 1. The mean displacement (left) and square displacement (right) Figure 1 shows averaged square displacement 〈x2n〉 and mean displacement 〈xn〉 as a function of diffusion time tn. Here α is a fraction of low-level energies ǫ1x. Very straight lines are observed, the slope of which determines the diffusion coefficient and drift velocity. In what follow we shall use the quantities t1 = τ0 exp((ε1t − ε1s)β) and D1 = a 2/2t1 = D01 exp((ε1t − ε1s)β) which is the diffusion coefficient for an ordered chain with constant site and transition energies ǫ1s and ǫ1t, respectively. The result of the simulation for both long and short chains is shown in Figure 2. It can be seen that while both simulation methods provide an identical result the size of the chains employed here is significantly different. The discrepancy obtained here is less than 0.1%. Therefore, the non-force method allows one to attaining a precise result based on a system of relatively small size and it can be employed for simulating the diffusion on a network of sites in a real AM model which consists of several thousand atoms. An important point with respect to matter diffusion is that AMs are characterized by two kinds of disorder which have deeply different properties: on the one hand, there are site disorders corresponding to the so-called random-trapping model, on the other hand, there are saddle point disorders, corresponding to the random-hoping model [1]. 146 A study of diffusion in a disordered chain: computer simulation Figure 2. The diffusion coefficient of a two-level distribution system at temperature (ε2t − ε1t)β = 3 (left) in a uniform distribution system (right) Figure 3. Temperature dependence of the diffusion coefficient for saddle (left) and site systems (right) In the present work these disorders are examined separately on three corresponding chains: (1) the site system, where site energies are constant and equal to ǫ1s, but transition energies amount to random value, (2) the saddle system, where transition energies have constant the value of ǫ1t and site energies is chosen randomly, (3) the saddle-site system, where both site and transition energies are chosen randomly but their type of energy distribution is the same, i.e. it is two-level or uniform distribution with the same ε2x−ε1x. Figure 3 presents the temperature dependence of the diffusion coefficient for site and 147 Trinh Van Mung, Pham Khac Hung and Nguyen Thi Thu Ha saddle systems. The notable feature of these plots is the linearity for a limited interval of temperature in each case which could be fitted to following formula ln(D/D1) = ax(ε2x − ε1x)β + bx (2.10) We examine two temperature intervals (ε2t − ε1x)β = (1.5 - 3.5) and (7 - 9). The fitting parameters are listed in Table 1 show that parameter at at low temperatures is higher than at high temperatures which implies a deviation from Arrhenius behavior along the entire interval (1.5 - 9). Table 1. Parameters at and bt of saddle point and site systems Two-level distribution Uniform distributionα = 0.2 α = 0.4 α = 0.6 α = 0.8 Interval of (ε2t − ε1t)β 1.5-3.5 1.5-3.5 1.5-3.5 1.5-3.5 1.5-3.5 at -0.9773 -0.9431 -0.8798 -0.7281 -0.6867 bt 0.1395 0.2936 0.4736 0.6343 0.2044 Interval of (ε2t − ε1t)β 7.0-9.0 7.0-9.0 7.0-9.0 7.0-9.0 7.0-9.0 at -0.9995 -0.9979 -0.9953 -0.9906 -0.8743 bt 0.2184 0.4885 0.8921 1.6069 1.073 Interval of (ε2s − ε1s)β 7.0-9.0 7.0-9.0 7.0-9.0 7.0-9.0 7.0-9.0 as 0 0 0 0 0.1276 bs 1.6405 0.9243 0.5167 0.2193 1.082 Interval of (ε2s − ε1s)β 7.0-9.0 7.0-9.0 7.0-9.0 7.0-9.0 7.0-9.0 ats -0.9982 -0.9982 -0.9967 -0.9903 -0.7471 bts 1.797 1.382 1.4213 1.7954 2.1428 In the case of α = 0.8, parameter at changes from -0.7281 to -0.9906, i.e. about 1.36 times, indicating strong deviation from Arrhenius behavior. From Eq. 2.10 the diffusion coefficient can be reduced to D =D1 exp(ax(ε2x − ε1x)β + bx) =D01 exp(bx) exp(ax(ε2x − ε1x)β − (ε1t − ε1s)β) (2.11) Comparing this to an ordered chain diffusion, where site and transition energies are equal to the ǫ1s and ǫ1t, the effective activation energy changes by −ax(ε2t − ε1x) and the pre-exponential factor by a factor of exp(bx). Because the value of ax for the saddle system is negative, its effective activation energy increases. Meanwhile, the parameter bx is always positive so the pre-exponential factor increases 1.2 to 5.6 times in comparison with the ordered chain. In the case of two-level distribution and at low temperature, the effective activation energy approaches the value ε2t − ε1s or ε1t − ε1s for the saddle and site system, respectively, i.e., the highest energy level of the saddle point and the lowest energy level of the site have a dominant contribution to the activation energy in the low 148 A study of diffusion in a disordered chain: computer simulation temperature region (see Table 1). Figure 4 presents the results for the saddle-site system. It is noted that in the interval of (7 - 9), linearity is observed in all plots in Figure 4. Parameters at and bt are almost equal to the sum of corresponding site and saddle system parameters (see Table). Therefore, the pre-exponential factor D0 of a disordered chain is bigger when both types of energetic disorders exist. The effect of compensation proposed in ref.[1] is not observed in our simulation. Figure 4. Temperature dependence of the diffusion coefficient in a saddle-site system 3. Conclusion A new simulation technique is proposed which allows the study of diffusion in a disordered system of moderate size and with periodic boundary conditions. A simulation carried out which makes use of this method on a chain which is 20 times smaller provides an accuracy equal to that using the usual method. The influence of both types of disorders has been investigated in relation to disordered chain. Some features of diffusion in those systems have been revealed as follows: 1-Despite the existence of strong energetic disorder in the system, including site and saddle types, Arrhenius behavior is always observed within a limited temperature range, particularly at low temperatures. Both types of disorder increases the pre-exponential coefficient D0 by a factor of 1.2 to 9.1 depending on the actual distribution of site and transition energy. On the other hand, along the entire temperature interval investigated we observe strong deviation from Arrhenius law except when there was two-level distribution with α = 0.2. 2-Concerning the effective activation energy, it approaches the value of (ε2t − ε1s) or (ε1t − ε1s) for the saddle and site system, respectively, as temperature decreases. This means that the highest energy level of a saddle point and the lowest energy level of a site have a dominant contribution in determining the activation energy at low temperatures. 149 Trinh Van Mung, Pham Khac Hung and Nguyen Thi Thu Ha 3-When there exist both site and saddle disorders in the chain, parameters ats and bts derived from Arrhenius plot are almost equal to the sum of the parameters of the site and saddle systems. A compensation effect for pre-exponential factors was not observed in our simulation. REFERENCES [1] Y.Limoge, Science and Engineering, A226-228 228. [2] W. Frank, A. Horner, P. Scharwaechter and H. Kronmiiller, 1994. Materials Science and Engineering, A 179/A 180 36. [3] V.Naundorf, M.P.Macht, A.S.Bakai, N.Lazarev, 1998. J. of Non-Crystalline Solids. 224 122. [4] N. Eliaz, D. Fuks and D. Eliezer, 1999. Acta Materialia. 47, 10 2981. [5] M.E. Fisher, 1984. J. Stat. Phys. 34 667. [6] J.W. van de Leur, A.Yu. Orlov, 2009. Physics Letters A, 373, 31 2675. [7] Belashchenco D.K.,P.K.Hung, 1986. Izv.Akad.Nauk. SSSR, Metally, 2 57. [8] Limoge Y., Bocquet J. L., 1990. Phys.Rev.Lett. 65 60. [9] Vo Van Hoang, 2004. J. of Phys. B 348 347. [10] P.K. Hung, L.T. Vinh, P.H. Kien, 2010. J. of Non-Crystal. Sol., 356 1000. 150